source: trunk/Cbc/src/CbcNode.cpp @ 1240

Last change on this file since 1240 was 1240, checked in by forrest, 10 years ago

adding stuff for heuristics

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 171.3 KB
Line 
1/* $Id: CbcNode.cpp 1240 2009-10-02 18:41:44Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8
9#include "CbcConfig.h"
10
11#include <string>
12//#define CBC_DEBUG 1
13//#define CHECK_CUT_COUNTS
14//#define CHECK_NODE
15//#define CBC_CHECK_BASIS
16#include <cassert>
17#include <cfloat>
18#define CUTS
19#include "OsiSolverInterface.hpp"
20#include "OsiChooseVariable.hpp"
21#include "OsiAuxInfo.hpp"
22#include "OsiSolverBranch.hpp"
23#include "CoinWarmStartBasis.hpp"
24#include "CoinTime.hpp"
25#include "CbcModel.hpp"
26#include "CbcNode.hpp"
27#include "CbcStatistics.hpp"
28#include "CbcStrategy.hpp"
29#include "CbcBranchActual.hpp"
30#include "CbcBranchDynamic.hpp"
31#include "OsiRowCut.hpp"
32#include "OsiRowCutDebugger.hpp"
33#include "OsiCuts.hpp"
34#include "CbcCountRowCut.hpp"
35#include "CbcFeasibilityBase.hpp"
36#include "CbcMessage.hpp"
37#ifdef COIN_HAS_CLP
38#include "OsiClpSolverInterface.hpp"
39#include "ClpSimplexOther.hpp"
40#endif
41using namespace std;
42#include "CglCutGenerator.hpp"
43// Default Constructor
44CbcNodeInfo::CbcNodeInfo ()
45  :
46  numberPointingToThis_(0),
47  parent_(NULL),
48  parentBranch_(NULL),
49  owner_(NULL),
50  numberCuts_(0),
51  nodeNumber_(0),
52  cuts_(NULL),
53  numberRows_(0),
54  numberBranchesLeft_(0),
55  active_(7)
56{
57#ifdef CHECK_NODE
58  printf("CbcNodeInfo %x Constructor\n",this);
59#endif
60}
61
62void
63CbcNodeInfo::setParentBasedData()
64{
65  if (parent_) {
66    numberRows_ = parent_->numberRows_+parent_->numberCuts_;
67    //parent_->increment();
68    if (parent_->owner()) {
69      const OsiBranchingObject* br = parent_->owner()->branchingObject();
70      assert(br);
71      parentBranch_ = br->clone();
72    }
73  }
74}
75
76void
77CbcNodeInfo::unsetParentBasedData()
78{
79  if (parent_) {
80    numberRows_ = 0;
81    if (parent_->owner()) {
82      delete parentBranch_;
83      parentBranch_ = NULL;
84    }
85  }
86}
87
88#if 0
89// Constructor given parent
90CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent)
91  :
92  numberPointingToThis_(2),
93  parent_(parent),
94  parentBranch_(NULL),
95  owner_(NULL),
96  numberCuts_(0),
97  nodeNumber_(0),
98  cuts_(NULL),
99  numberRows_(0),
100  numberBranchesLeft_(2),
101  active_(7)
102{
103#ifdef CHECK_NODE
104  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
105#endif
106  //setParentBasedData();
107}
108#endif
109
110// Copy Constructor
111CbcNodeInfo::CbcNodeInfo (const CbcNodeInfo & rhs)
112  :
113  numberPointingToThis_(rhs.numberPointingToThis_),
114  parent_(rhs.parent_),
115  parentBranch_(NULL),
116  owner_(rhs.owner_),
117  numberCuts_(rhs.numberCuts_),
118  nodeNumber_(rhs.nodeNumber_),
119  cuts_(NULL),
120  numberRows_(rhs.numberRows_),
121  numberBranchesLeft_(rhs.numberBranchesLeft_),
122  active_(rhs.active_)
123{
124#ifdef CHECK_NODE
125  printf("CbcNodeInfo %x Copy constructor\n",this);
126#endif
127  if (numberCuts_) {
128    cuts_ = new CbcCountRowCut * [numberCuts_];
129    int n=0;
130    for (int i=0;i<numberCuts_;i++) {
131      CbcCountRowCut * thisCut = rhs.cuts_[i];
132      if (thisCut) {
133        // I think this is correct - new one should take priority
134        thisCut->setInfo(this,n);
135        thisCut->increment(numberBranchesLeft_); 
136        cuts_[n++] = thisCut;
137      }
138    }
139    numberCuts_=n;
140  }
141  if (rhs.parentBranch_) {
142    parentBranch_ = rhs.parentBranch_->clone();
143  }
144}
145// Constructor given parent and owner
146CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent, CbcNode * owner)
147  :
148  numberPointingToThis_(2),
149  parent_(parent),
150  parentBranch_(NULL),
151  owner_(owner),
152  numberCuts_(0),
153  nodeNumber_(0),
154  cuts_(NULL),
155  numberRows_(0),
156  numberBranchesLeft_(2),
157  active_(7)
158{
159#ifdef CHECK_NODE
160  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
161#endif
162  //setParentBasedData();
163}
164
165/**
166   Take care to detach from the owning CbcNode and decrement the reference
167   count in the parent.  If this is the last nodeInfo object pointing to the
168   parent, make a recursive call to delete the parent.
169*/
170CbcNodeInfo::~CbcNodeInfo()
171{
172#ifdef CHECK_NODE
173  printf("CbcNodeInfo %x Destructor parent %x\n",this,parent_);
174#endif
175 
176  assert(!numberPointingToThis_);
177  // But there may be some left (max nodes?)
178  for (int i=0;i<numberCuts_;i++) {
179    if (cuts_[i]) {
180#ifndef GLOBAL_CUTS_JUST_POINTERS
181      delete cuts_[i];
182#else
183      if (cuts_[i]->globallyValidAsInteger()!=2)
184        delete cuts_[i];
185#endif
186    }
187  }
188  delete [] cuts_;
189  if (owner_) 
190    owner_->nullNodeInfo();
191  if (parent_) {
192    int numberLinks = parent_->decrement();
193    if (!numberLinks) delete parent_;
194  }
195  delete parentBranch_;
196}
197
198
199//#define ALLCUTS
200void
201CbcNodeInfo::decrementCuts(int change)
202{
203  int i;
204  // get rid of all remaining if negative
205  int changeThis;
206  if (change<0)
207    changeThis = numberBranchesLeft_;
208  else
209    changeThis = change;
210  // decrement cut counts
211  for (i=0;i<numberCuts_;i++) {
212    if (cuts_[i]) {
213      int number = cuts_[i]->decrement(changeThis);
214      if (!number) {
215        //printf("info %x del cut %d %x\n",this,i,cuts_[i]);
216#ifndef GLOBAL_CUTS_JUST_POINTERS
217        delete cuts_[i];
218#else
219        if (cuts_[i]->globallyValidAsInteger()!=2)
220          delete cuts_[i];
221#endif
222        cuts_[i]=NULL;
223      }
224    }
225  }
226}
227void
228CbcNodeInfo::incrementCuts(int change)
229{
230  int i;
231  assert (change>0);
232  // increment cut counts
233  for (i=0;i<numberCuts_;i++) {
234    if (cuts_[i]) 
235      cuts_[i]->increment(change);
236  }
237}
238void
239CbcNodeInfo::decrementParentCuts(CbcModel * model,int change)
240{
241  if (parent_) {
242    // get rid of all remaining if negative
243    int changeThis;
244    if (change<0)
245      changeThis = numberBranchesLeft_;
246    else
247      changeThis = change;
248    int i;
249    // Get over-estimate of space needed for basis
250    CoinWarmStartBasis & dummy = model->workingBasis();
251    dummy.setSize(0,numberRows_+numberCuts_);
252    buildRowBasis(dummy);
253    /* everything is zero (i.e. free) so we can use to see
254       if latest basis */
255    CbcNodeInfo * thisInfo = parent_;
256    while (thisInfo) 
257      thisInfo = thisInfo->buildRowBasis(dummy);
258    // decrement cut counts
259    thisInfo = parent_;
260    int numberRows=numberRows_;
261    while (thisInfo) {
262      for (i=thisInfo->numberCuts_-1;i>=0;i--) {
263        CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
264#ifdef ALLCUTS
265        status = CoinWarmStartBasis::isFree;
266#endif
267        if (thisInfo->cuts_[i]) {
268          int number=1;
269          if (status!=CoinWarmStartBasis::basic) {
270            // tight - drop 1 or 2
271            if (change<0)
272              number = thisInfo->cuts_[i]->decrement(changeThis);
273            else
274              number = thisInfo->cuts_[i]->decrement(change);
275          }
276          if (!number) {
277#ifndef GLOBAL_CUTS_JUST_POINTERS
278            delete thisInfo->cuts_[i];
279#else
280            if (thisInfo->cuts_[i]->globallyValidAsInteger()!=2)
281              delete thisInfo->cuts_[i];
282#endif
283            thisInfo->cuts_[i]=NULL;
284          }
285        }
286      }
287      thisInfo = thisInfo->parent_;
288    }
289  }
290}
291#if 0
292void
293CbcNodeInfo::incrementParentCuts(CbcModel * model, int change)
294{
295  if (parent_) {
296    int i;
297    // Get over-estimate of space needed for basis
298    CoinWarmStartBasis & dummy = model->workingBasis();
299    dummy.setSize(0,numberRows_+numberCuts_);
300    /* everything is zero (i.e. free) so we can use to see
301       if latest basis */
302    buildRowBasis(dummy);
303    CbcNodeInfo * thisInfo = parent_;
304    while (thisInfo)
305      thisInfo = thisInfo->buildRowBasis(dummy);
306    // increment cut counts
307    thisInfo = parent_;
308    int numberRows=numberRows_;
309    while (thisInfo) {
310      for (i=thisInfo->numberCuts_-1;i>=0;i--) {
311        CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
312#ifdef ALLCUTS
313        status = CoinWarmStartBasis::isFree;
314#endif
315        if (thisInfo->cuts_[i]&&status!=CoinWarmStartBasis::basic) {
316          thisInfo->cuts_[i]->increment(change);
317        }
318      }
319      thisInfo = thisInfo->parent_;
320    }
321  }
322}
323#endif
324/*
325  Append cuts to the cuts_ array in a nodeInfo. The initial reference count
326  is set to numberToBranchOn, which will normally be the number of arms
327  defined for the CbcBranchingObject attached to the CbcNode that owns this
328  CbcNodeInfo.
329*/
330void
331CbcNodeInfo::addCuts (OsiCuts & cuts, int numberToBranchOn,
332                      /*int * whichGenerator,*/int numberPointingToThis)
333{
334  int numberCuts = cuts.sizeRowCuts();
335  if (numberCuts) {
336    int i;
337    if (!numberCuts_) {
338      cuts_ = new CbcCountRowCut * [numberCuts];
339    } else {
340      CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
341      memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *));
342      delete [] cuts_;
343      cuts_ = temp;
344    }
345    for (i=0;i<numberCuts;i++) {
346      CbcCountRowCut * thisCut = new CbcCountRowCut(*cuts.rowCutPtr(i),
347                                                    this,numberCuts_,
348                                                    -1,numberPointingToThis);
349      thisCut->increment(numberToBranchOn); 
350      cuts_[numberCuts_++] = thisCut;
351#ifdef CBC_DEBUG
352#if CBC_DEBUG>1
353      int n=thisCut->row().getNumElements();
354      printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
355             thisCut->ub());
356      int j;
357      const int * index = thisCut->row().getIndices();
358      const double * element = thisCut->row().getElements();
359      for (j=0;j<n;j++) {
360        printf(" (%d,%g)",index[j],element[j]);
361        assert(fabs(element[j])>1.00e-12);
362      }
363      printf("\n");
364#else
365      int n=thisCut->row().getNumElements();
366      int j;
367      const double * element = thisCut->row().getElements();
368      for (j=0;j<n;j++) {
369        assert(fabs(element[j])>1.00e-12);
370      }
371#endif
372#endif
373    }
374  }
375}
376
377void
378CbcNodeInfo::addCuts(int numberCuts, CbcCountRowCut ** cut, 
379                     int numberToBranchOn)
380{
381  if (numberCuts) {
382    int i;
383    if (!numberCuts_) {
384      cuts_ = new CbcCountRowCut * [numberCuts];
385    } else {
386      CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
387      memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *));
388      delete [] cuts_;
389      cuts_ = temp;
390    }
391    for (i=0;i<numberCuts;i++) {
392      CbcCountRowCut * thisCut = cut[i];
393      thisCut->setInfo(this,numberCuts_);
394      //printf("info %x cut %d %x\n",this,i,thisCut);
395      thisCut->increment(numberToBranchOn); 
396      cuts_[numberCuts_++] = thisCut;
397#ifdef CBC_DEBUG
398      int n=thisCut->row().getNumElements();
399#if CBC_DEBUG>1
400      printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
401             thisCut->ub());
402#endif
403      int j;
404#if CBC_DEBUG>1
405      const int * index = thisCut->row().getIndices();
406#endif
407      const double * element = thisCut->row().getElements();
408      for (j=0;j<n;j++) {
409#if CBC_DEBUG>1
410        printf(" (%d,%g)",index[j],element[j]);
411#endif
412        assert(fabs(element[j])>1.00e-12);
413      }
414      printf("\n");
415#endif
416    }
417  }
418}
419
420// delete cuts
421void
422CbcNodeInfo::deleteCuts(int numberToDelete, CbcCountRowCut ** cuts)
423{
424  int i;
425  int j;
426  int last=-1;
427  for (i=0;i<numberToDelete;i++) {
428    CbcCountRowCut * next = cuts[i];
429    for (j=last+1;j<numberCuts_;j++) {
430      if (next==cuts_[j])
431        break;
432    }
433    if (j==numberCuts_) {
434      // start from beginning
435      for (j=0;j<last;j++) {
436        if (next==cuts_[j])
437          break;
438      }
439      assert(j<last);
440    }
441    last=j;
442    int number = cuts_[j]->decrement();
443    if (!number) {
444#ifndef GLOBAL_CUTS_JUST_POINTERS
445      delete cuts_[j];
446#else
447      if (cuts_[j]->globallyValidAsInteger()!=2)
448        delete cuts_[j];
449#endif
450    }
451    cuts_[j]=NULL;
452  }
453  j=0;
454  for (i=0;i<numberCuts_;i++) {
455    if (cuts_[i])
456      cuts_[j++]=cuts_[i];
457  }
458  numberCuts_ = j;
459}
460
461// delete cuts
462void
463CbcNodeInfo::deleteCuts(int numberToDelete, int * which)
464{
465  int i;
466  for (i=0;i<numberToDelete;i++) {
467    int iCut=which[i];
468    int number = cuts_[iCut]->decrement();
469    if (!number) {
470#ifndef GLOBAL_CUTS_JUST_POINTERS
471      delete cuts_[iCut];
472#else
473      if (cuts_[iCut]->globallyValidAsInteger()!=2)
474        delete cuts_[iCut];
475#endif
476    }
477    cuts_[iCut]=NULL;
478  }
479  int j=0;
480  for (i=0;i<numberCuts_;i++) {
481    if (cuts_[i])
482      cuts_[j++]=cuts_[i];
483  }
484  numberCuts_ = j;
485}
486
487// Really delete a cut
488void 
489CbcNodeInfo::deleteCut(int whichOne)
490{
491  assert(whichOne<numberCuts_);
492  cuts_[whichOne]=NULL;
493}
494/* Deactivate node information.
495   1 - bounds
496   2 - cuts
497   4 - basis!
498*/
499void 
500CbcNodeInfo::deactivate(int mode)
501{
502  active_ &= (~mode);
503}
504
505CbcFullNodeInfo::CbcFullNodeInfo() :
506  CbcNodeInfo(),
507  basis_(),
508  numberIntegers_(0),
509  lower_(NULL),
510  upper_(NULL)
511{
512}
513CbcFullNodeInfo::CbcFullNodeInfo(CbcModel * model,
514                                 int numberRowsAtContinuous) :
515  CbcNodeInfo(NULL, model->currentNode())
516{
517  OsiSolverInterface * solver = model->solver();
518  numberRows_ = numberRowsAtContinuous;
519  numberIntegers_ = model->numberIntegers();
520  int numberColumns = model->getNumCols();
521  lower_ = new double [numberColumns];
522  upper_ = new double [numberColumns];
523  const double * lower = solver->getColLower();
524  const double * upper = solver->getColUpper();
525  int i;
526 
527  for (i=0;i<numberColumns;i++) {
528    lower_[i]=lower[i];
529    upper_[i]=upper[i];
530  }
531 
532  basis_ =  dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
533}
534
535CbcFullNodeInfo::CbcFullNodeInfo(const CbcFullNodeInfo & rhs) :
536  CbcNodeInfo(rhs)
537{ 
538  basis_= dynamic_cast<CoinWarmStartBasis *>(rhs.basis_->clone()) ;
539  numberIntegers_=rhs.numberIntegers_;
540  lower_=NULL;
541  upper_=NULL;
542  if (rhs.lower_!=NULL) {
543    int numberColumns = basis_->getNumStructural();
544    lower_ = new double [numberColumns];
545    upper_ = new double [numberColumns];
546    assert (upper_!=NULL);
547    memcpy(lower_,rhs.lower_,numberColumns*sizeof(double));
548    memcpy(upper_,rhs.upper_,numberColumns*sizeof(double));
549  }
550}
551
552CbcNodeInfo * 
553CbcFullNodeInfo::clone() const
554{ 
555  return (new CbcFullNodeInfo(*this));
556}
557
558CbcFullNodeInfo::~CbcFullNodeInfo ()
559{
560  delete basis_ ;
561  delete [] lower_;
562  delete [] upper_;
563}
564
565/*
566  The basis supplied as a parameter is deleted and replaced with a new basis
567  appropriate for the node, and lower and upper bounds on variables are
568  reset according to the stored bounds arrays. Any cuts associated with this
569  node are added to the list in addCuts, but not actually added to the
570  constraint system in the model.
571 
572  Why pass in a basis at all? The short answer is ``We need the parameter to
573  pass out a basis, so might as well use it to pass in the size.''
574 
575  A longer answer is that in practice we take a memory allocation hit up in
576  addCuts1 (the only place applyToModel is called) when we setSize() the
577  basis that's passed in. It's immediately tossed here in favour of a clone
578  of the basis attached to this nodeInfo. This can probably be fixed, given
579  a bit of thought.
580*/
581
582void CbcFullNodeInfo::applyToModel (CbcModel *model,
583                                    CoinWarmStartBasis *&basis,
584                                    CbcCountRowCut **addCuts,
585                                    int &currentNumberCuts) const 
586 
587{ OsiSolverInterface *solver = model->solver() ;
588 
589  // branch - do bounds
590  assert (active_==7||active_==15);
591  int i;
592  solver->setColLower(lower_);
593  solver->setColUpper(upper_);
594  int numberColumns = model->getNumCols();
595  // move basis - but make sure size stays
596  // for bon-min - should not be needed int numberRows = model->getNumRows();
597  int numberRows=basis->getNumArtificial();
598  delete basis ;
599  if (basis_) {
600    basis = dynamic_cast<CoinWarmStartBasis *>(basis_->clone()) ;
601    basis->resize(numberRows,numberColumns);
602#ifdef CBC_CHECK_BASIS
603    std::cout << "Basis (after applying root " <<this<<") "<< std::endl ;
604    basis->print() ;
605#endif   
606  } else {
607    // We have a solver without a basis
608    basis=NULL;
609  }
610  for (i=0;i<numberCuts_;i++) 
611    addCuts[currentNumberCuts+i]= cuts_[i];
612  currentNumberCuts += numberCuts_;
613  assert(!parent_);
614  return ;
615}
616// Just apply bounds to one variable (1=>infeasible)
617int 
618CbcFullNodeInfo::applyBounds(int iColumn, double & lower, double & upper,int force) 
619{
620  if ((force&&1)==0) {
621    if (lower>lower_[iColumn])
622      printf("%d odd lower going from %g to %g\n",iColumn,lower,lower_[iColumn]);
623    lower = lower_[iColumn];
624  } else {
625    lower_[iColumn]=lower;
626  }
627  if ((force&&2)==0) {
628    if (upper<upper_[iColumn])
629      printf("%d odd upper going from %g to %g\n",iColumn,upper,upper_[iColumn]);
630    upper = upper_[iColumn];
631  } else {
632    upper_[iColumn]=upper;
633  }
634  return (upper_[iColumn]>=lower_[iColumn]) ? 0 : 1;
635}
636
637/* Builds up row basis backwards (until original model).
638   Returns NULL or previous one to apply .
639   Depends on Free being 0 and impossible for cuts
640*/
641CbcNodeInfo * 
642CbcFullNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
643{
644  const unsigned int * saved = 
645    reinterpret_cast<const unsigned int *> (basis_->getArtificialStatus());
646  unsigned int * now = 
647    reinterpret_cast<unsigned int *> (basis.getArtificialStatus());
648  int number=basis_->getNumArtificial()>>4;;
649  int i;
650  for (i=0;i<number;i++) { 
651    if (!now[i])
652      now[i] = saved[i];
653  }
654  return NULL;
655}
656
657
658// Default constructor
659CbcPartialNodeInfo::CbcPartialNodeInfo()
660 
661  : CbcNodeInfo(),
662    basisDiff_(NULL),
663    variables_(NULL),
664    newBounds_(NULL),
665    numberChangedBounds_(0)
666   
667{ /* this space intentionally left blank */ }
668
669// Constructor from current state
670CbcPartialNodeInfo::CbcPartialNodeInfo (CbcNodeInfo *parent, CbcNode *owner,
671                                        int numberChangedBounds,
672                                        const int *variables,
673                                        const double *boundChanges,
674                                        const CoinWarmStartDiff *basisDiff)
675  : CbcNodeInfo(parent,owner)
676{
677  basisDiff_ = basisDiff->clone() ;
678#ifdef CBC_CHECK_BASIS
679  std::cout << "Constructor (" <<this<<") "<< std::endl ;
680#endif   
681 
682  numberChangedBounds_ = numberChangedBounds;
683  int size = numberChangedBounds_*(sizeof(double)+sizeof(int));
684  char * temp = new char [size];
685  newBounds_ = reinterpret_cast<double *> (temp);
686  variables_ = reinterpret_cast<int *> (newBounds_+numberChangedBounds_);
687 
688  int i ;
689  for (i=0;i<numberChangedBounds_;i++) {
690    variables_[i]=variables[i];
691    newBounds_[i]=boundChanges[i];
692  }
693}
694
695CbcPartialNodeInfo::CbcPartialNodeInfo (const CbcPartialNodeInfo & rhs)
696 
697  : CbcNodeInfo(rhs)
698   
699{ basisDiff_ = rhs.basisDiff_->clone() ;
700 
701#ifdef CBC_CHECK_BASIS
702  std::cout << "Copy constructor (" <<this<<") from "<<this<< std::endl ;
703#endif   
704  numberChangedBounds_ = rhs.numberChangedBounds_;
705  int size = numberChangedBounds_*(sizeof(double)+sizeof(int));
706  char * temp = new char [size];
707  newBounds_ = reinterpret_cast<double *> (temp);
708  variables_ = reinterpret_cast<int *> (newBounds_+numberChangedBounds_);
709 
710  int i ;
711  for (i=0;i<numberChangedBounds_;i++) {
712    variables_[i]=rhs.variables_[i];
713    newBounds_[i]=rhs.newBounds_[i];
714  }
715}
716
717CbcNodeInfo * 
718CbcPartialNodeInfo::clone() const
719{ 
720  return (new CbcPartialNodeInfo(*this));
721}
722
723
724CbcPartialNodeInfo::~CbcPartialNodeInfo ()
725{
726  delete basisDiff_ ;
727  delete [] newBounds_;
728}
729
730
731/**
732   The basis supplied as a parameter is incrementally modified, and lower and
733   upper bounds on variables in the model are incrementally modified. Any
734   cuts associated with this node are added to the list in addCuts.
735*/
736
737void CbcPartialNodeInfo::applyToModel (CbcModel *model,
738                                       CoinWarmStartBasis *&basis,
739                                       CbcCountRowCut **addCuts,
740                                       int &currentNumberCuts) const 
741 
742{ OsiSolverInterface *solver = model->solver();
743  if ((active_&4)!=0) {
744    basis->applyDiff(basisDiff_) ;
745#ifdef CBC_CHECK_BASIS
746    std::cout << "Basis (after applying " <<this<<") "<< std::endl ;
747    basis->print() ;
748#endif   
749  }
750 
751  // branch - do bounds
752  int i;
753  if ((active_&1)!=0) {
754    for (i=0;i<numberChangedBounds_;i++) {
755      int variable = variables_[i];
756      int k = variable&0x3fffffff;
757      if ((variable&0x80000000)==0) {
758        // lower bound changing
759        //#define CBC_PRINT2
760#ifdef CBC_PRINT2
761        if(solver->getColLower()[k]!=newBounds_[i])
762          printf("lower change for column %d - from %g to %g\n",k,solver->getColLower()[k],newBounds_[i]);
763#endif
764#ifndef NDEBUG
765        if ((variable&0x40000000)==0&&false) {
766          double oldValue = solver->getColLower()[k];
767          assert (newBounds_[i]>oldValue-1.0e-8);
768          if (newBounds_[i]<oldValue+1.0e-8)
769            printf("bad null lower change for column %d - bound %g\n",k,oldValue);
770        }
771#endif
772        solver->setColLower(k,newBounds_[i]);
773      } else {
774        // upper bound changing
775#ifdef CBC_PRINT2
776        if(solver->getColUpper()[k]!=newBounds_[i])
777          printf("upper change for column %d - from %g to %g\n",k,solver->getColUpper()[k],newBounds_[i]);
778#endif
779#ifndef NDEBUG
780        if ((variable&0x40000000)==0&&false) {
781          double oldValue = solver->getColUpper()[k];
782          assert (newBounds_[i]<oldValue+1.0e-8);
783          if (newBounds_[i]>oldValue-1.0e-8)
784            printf("bad null upper change for column %d - bound %g\n",k,oldValue);
785        }
786#endif
787        solver->setColUpper(k,newBounds_[i]);
788      }
789    }
790  }
791  if ((active_&2)!=0) {
792    for (i=0;i<numberCuts_;i++) {
793      addCuts[currentNumberCuts+i]= cuts_[i];
794      if (cuts_[i]&&model->messageHandler()->logLevel()>4) {
795        cuts_[i]->print();
796      }
797    }
798   
799    currentNumberCuts += numberCuts_;
800  }
801  return ;
802}
803// Just apply bounds to one variable (1=>infeasible)
804int
805CbcPartialNodeInfo::applyBounds(int iColumn, double & lower, double & upper,int force) 
806{
807  // branch - do bounds
808  int i;
809  int found=0;
810  double newLower = -COIN_DBL_MAX;
811  double newUpper = COIN_DBL_MAX;
812  for (i=0;i<numberChangedBounds_;i++) {
813    int variable = variables_[i];
814    int k = variable&0x3fffffff;
815    if (k==iColumn) {
816      if ((variable&0x80000000)==0) {
817        // lower bound changing
818        found |= 1;
819        newLower = CoinMax(newLower,newBounds_[i]);
820        if ((force&1)==0) {
821          if (lower>newBounds_[i])
822            printf("%d odd lower going from %g to %g\n",iColumn,lower,newBounds_[i]);
823          lower = newBounds_[i];
824        } else {
825          newBounds_[i]=lower;
826          variables_[i] |= 0x40000000; // say can go odd way
827        }
828      } else {
829        // upper bound changing
830        found |= 2;
831        newUpper = CoinMin(newUpper,newBounds_[i]);
832        if ((force&2)==0) {
833          if (upper<newBounds_[i])
834            printf("%d odd upper going from %g to %g\n",iColumn,upper,newBounds_[i]);
835          upper = newBounds_[i];
836        } else {
837          newBounds_[i]=upper;
838          variables_[i] |= 0x40000000; // say can go odd way
839        }
840      }
841    }
842  }
843  newLower = CoinMax(newLower,lower);
844  newUpper = CoinMin(newUpper,upper);
845  int nAdd=0;
846  if ((force&2)!=0&&(found&2)==0) {
847    // need to add new upper
848    nAdd++;
849  }
850  if ((force&1)!=0&&(found&1)==0) {
851    // need to add new lower
852    nAdd++;
853  }
854  if (nAdd) { 
855    int size = (numberChangedBounds_+nAdd)*(sizeof(double)+sizeof(int));
856    char * temp = new char [size];
857    double * newBounds = reinterpret_cast<double *> (temp);
858    int * variables = reinterpret_cast<int *> (newBounds+numberChangedBounds_+nAdd);
859   
860    int i ;
861    for (i=0;i<numberChangedBounds_;i++) {
862      variables[i]=variables_[i];
863      newBounds[i]=newBounds_[i];
864    }
865    delete [] newBounds_;
866    newBounds_ = newBounds;
867    variables_ = variables;
868    if ((force&2)!=0&&(found&2)==0) {
869      // need to add new upper
870      int variable = iColumn | 0x80000000;
871      variables_[numberChangedBounds_]=variable;
872      newBounds_[numberChangedBounds_++]=newUpper;
873    }
874    if ((force&1)!=0&&(found&1)==0) {
875      // need to add new lower
876      int variable = iColumn;
877      variables_[numberChangedBounds_]=variable;
878      newBounds_[numberChangedBounds_++]=newLower;
879    }
880  }
881 
882  return (newUpper>=newLower) ? 0 : 1;
883}
884
885/* Builds up row basis backwards (until original model).
886   Returns NULL or previous one to apply .
887   Depends on Free being 0 and impossible for cuts
888*/
889
890CbcNodeInfo * 
891CbcPartialNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
892 
893{ basis.applyDiff(basisDiff_) ;
894 
895  return parent_ ; }
896CbcNode::CbcNode() :
897  nodeInfo_(NULL),
898  objectiveValue_(1.0e100),
899  guessedObjectiveValue_(1.0e100),
900  sumInfeasibilities_(0.0),
901  branch_(NULL),
902  depth_(-1),
903  numberUnsatisfied_(0),
904  nodeNumber_(-1),
905  state_(0)
906{
907#ifdef CHECK_NODE
908  printf("CbcNode %x Constructor\n",this);
909#endif
910}
911// Print
912void 
913CbcNode::print() const
914{
915  printf("number %d obj %g depth %d sumun %g nunsat %d state %d\n",
916         nodeNumber_,objectiveValue_,depth_,sumInfeasibilities_,numberUnsatisfied_,state_);
917}
918CbcNode::CbcNode(CbcModel * model,
919                 CbcNode * lastNode) :
920  nodeInfo_(NULL),
921  objectiveValue_(1.0e100),
922  guessedObjectiveValue_(1.0e100),
923  sumInfeasibilities_(0.0),
924  branch_(NULL),
925  depth_(-1),
926  numberUnsatisfied_(0),
927  nodeNumber_(-1),
928  state_(0)
929{
930#ifdef CHECK_NODE
931  printf("CbcNode %x Constructor from model\n",this);
932#endif
933  model->setObjectiveValue(this,lastNode);
934 
935  if (lastNode) {
936    if (lastNode->nodeInfo_) {
937      lastNode->nodeInfo_->increment();
938    }
939  }
940  nodeNumber_= model->getNodeCount();
941}
942
943#define CBC_NEW_CREATEINFO
944#ifdef CBC_NEW_CREATEINFO
945
946/*
947  New createInfo, with basis manipulation hidden inside mergeBasis. Allows
948  solvers to override and carry over all information from one basis to
949  another.
950*/
951
952void
953CbcNode::createInfo (CbcModel *model,
954                     CbcNode *lastNode,
955                     const CoinWarmStartBasis *lastws,
956                     const double *lastLower, const double *lastUpper,
957                     int numberOldActiveCuts, int numberNewCuts)
958 
959{ OsiSolverInterface *solver = model->solver();
960  CbcStrategy *strategy = model->strategy();
961  /*
962    The root --- no parent. Create full basis and bounds information.
963  */
964  if (!lastNode)
965    { 
966      if (!strategy)
967        nodeInfo_=new CbcFullNodeInfo(model,solver->getNumRows());
968      else
969        nodeInfo_ = strategy->fullNodeInfo(model,solver->getNumRows());
970    } else {
971    /*
972      Not the root. Create an edit from the parent's basis & bound information.
973      This is not quite as straightforward as it seems. We need to reintroduce
974      cuts we may have dropped out of the basis, in the correct position, because
975      this whole process is strictly positional. Start by grabbing the current
976      basis.
977    */
978    bool mustDeleteBasis;
979    const CoinWarmStartBasis *ws =
980      dynamic_cast<const CoinWarmStartBasis*>(solver->getPointerToWarmStart(mustDeleteBasis));
981    assert(ws!=NULL); // make sure not volume
982    //int numberArtificials = lastws->getNumArtificial();
983    int numberColumns = solver->getNumCols();
984    int numberRowsAtContinuous = model->numberRowsAtContinuous();
985    int currentNumberCuts = model->currentNumberCuts();
986#   ifdef CBC_CHECK_BASIS
987    std::cout
988      << "Before expansion: orig " << numberRowsAtContinuous
989      << ", old " << numberOldActiveCuts
990      << ", new " << numberNewCuts
991      << ", current " << currentNumberCuts << "." << std::endl ;
992    ws->print();
993#   endif
994    /*
995      Clone the basis and resize it to hold the structural constraints, plus
996      all the cuts: old cuts, both active and inactive (currentNumberCuts),
997      and new cuts (numberNewCuts). This will become the expanded basis.
998    */
999    CoinWarmStartBasis *expanded = 
1000      dynamic_cast<CoinWarmStartBasis *>(ws->clone()) ;
1001    int iCompact = numberRowsAtContinuous+numberOldActiveCuts+numberNewCuts ;
1002    // int nPartial = numberRowsAtContinuous+currentNumberCuts;
1003    int iFull = numberRowsAtContinuous+currentNumberCuts+numberNewCuts;
1004    // int maxBasisLength = ((iFull+15)>>4)+((numberColumns+15)>>4);
1005    // printf("l %d full %d\n",maxBasisLength,iFull);
1006    expanded->resize(iFull,numberColumns);
1007#   ifdef CBC_CHECK_BASIS
1008    std::cout
1009      << "\tFull basis " << iFull << " rows, "
1010      << numberColumns << " columns; compact "
1011      << iCompact << " rows." << std::endl ;
1012#   endif
1013    /*
1014      Now flesh out the expanded basis. The clone already has the
1015      correct status information for the variables and for the structural
1016      (numberRowsAtContinuous) constraints. Any indices beyond nPartial must be
1017      cuts created while processing this node --- they can be copied en bloc
1018      into the correct position in the expanded basis. The space reserved for
1019      xferRows is a gross overestimate.
1020    */
1021    CoinWarmStartBasis::XferVec xferRows ;
1022    xferRows.reserve(iFull-numberRowsAtContinuous+1) ;
1023    if (numberNewCuts) {
1024      xferRows.push_back(
1025                         CoinWarmStartBasis::XferEntry(iCompact-numberNewCuts,
1026                                                       iFull-numberNewCuts,numberNewCuts)) ;
1027    }
1028    /*
1029      From nPartial down, record the entries we want to copy from the current
1030      basis (the entries for the active cuts; non-zero in the list returned
1031      by addedCuts). Fill the expanded basis with entries showing a status of
1032      basic for the deactivated (loose) cuts.
1033    */
1034    CbcCountRowCut **cut = model->addedCuts();
1035    iFull -= (numberNewCuts+1) ;
1036    iCompact -= (numberNewCuts+1) ;
1037    int runLen = 0 ;
1038    CoinWarmStartBasis::XferEntry entry(-1,-1,-1) ;
1039    while (iFull >= numberRowsAtContinuous) { 
1040      for ( ; iFull >= numberRowsAtContinuous &&
1041              cut[iFull-numberRowsAtContinuous] ; iFull--)
1042        runLen++ ;
1043      if (runLen) {
1044        iCompact -= runLen ;
1045        entry.first = iCompact+1 ;
1046        entry.second = iFull+1 ;
1047        entry.third = runLen ;
1048        runLen = 0 ;
1049        xferRows.push_back(entry) ;
1050      }
1051      for ( ; iFull >= numberRowsAtContinuous &&
1052              !cut[iFull-numberRowsAtContinuous] ; iFull--)
1053        expanded->setArtifStatus(iFull,CoinWarmStartBasis::basic);
1054    }
1055    /*
1056      Finally, call mergeBasis to copy over entries from the current basis to
1057      the expanded basis. Since we cloned the expanded basis from the active basis
1058      and haven't changed the number of variables, only row status entries need
1059      to be copied.
1060    */
1061    expanded->mergeBasis(ws,&xferRows,0) ;
1062   
1063#ifdef CBC_CHECK_BASIS
1064    std::cout << "Expanded basis:" << std::endl ;
1065    expanded->print() ;
1066    std::cout << "Diffing against:" << std::endl ;
1067    lastws->print() ;
1068#endif   
1069    assert (expanded->getNumArtificial()>=lastws->getNumArtificial());
1070#ifdef CLP_INVESTIGATE
1071    if (!expanded->fullBasis()) {
1072      int iFull = numberRowsAtContinuous+currentNumberCuts+numberNewCuts;
1073      printf("cont %d old %d new %d current %d full inc %d full %d\n",
1074             numberRowsAtContinuous,numberOldActiveCuts,numberNewCuts,
1075             currentNumberCuts,iFull,iFull-numberNewCuts);
1076    }
1077#endif
1078   
1079    /*
1080      Now that we have two bases in proper positional correspondence, creating
1081      the actual diff is dead easy.
1082     
1083      Note that we're going to compare the expanded basis here to the stripped
1084      basis (lastws) produced by addCuts. It doesn't affect the correctness (the
1085      diff process has no knowledge of the meaning of an entry) but it does
1086      mean that we'll always generate a whack of diff entries because the expanded
1087      basis is considerably larger than the stripped basis.
1088    */
1089    CoinWarmStartDiff *basisDiff = expanded->generateDiff(lastws) ;
1090   
1091    /*
1092      Diff the bound vectors. It's assumed the number of structural variables
1093      is not changing. For branching objects that change bounds on integer
1094      variables, we should see at least one bound change as a consequence
1095      of applying the branch that generated this subproblem from its parent.
1096      This need not hold for other types of branching objects (hyperplane
1097      branches, for example).
1098    */
1099    const double * lower = solver->getColLower();
1100    const double * upper = solver->getColUpper();
1101   
1102    double *boundChanges = new double [2*numberColumns] ;
1103    int *variables = new int [2*numberColumns] ;
1104    int numberChangedBounds=0;
1105   
1106    int i;
1107    for (i=0;i<numberColumns;i++) {
1108      if (lower[i]!=lastLower[i]) {
1109        variables[numberChangedBounds]=i;
1110        boundChanges[numberChangedBounds++]=lower[i];
1111      }
1112      if (upper[i]!=lastUpper[i]) {
1113        variables[numberChangedBounds]=i|0x80000000;
1114        boundChanges[numberChangedBounds++]=upper[i];
1115      }
1116#ifdef CBC_DEBUG
1117      if (lower[i] != lastLower[i]) {
1118        std::cout
1119          << "lower on " << i << " changed from "
1120          << lastLower[i] << " to " << lower[i] << std::endl ;
1121      }
1122      if (upper[i] != lastUpper[i]) {
1123        std::cout
1124          << "upper on " << i << " changed from "
1125          << lastUpper[i] << " to " << upper[i] << std::endl ;
1126      }
1127#endif
1128    }
1129#ifdef CBC_DEBUG
1130    std::cout << numberChangedBounds << " changed bounds." << std::endl ;
1131#endif
1132    //if (lastNode->branchingObject()->boundBranch())
1133    //assert (numberChangedBounds);
1134    /*
1135      Hand the lot over to the CbcPartialNodeInfo constructor, then clean up and
1136      return.
1137    */
1138    if (!strategy)
1139      nodeInfo_ =
1140        new CbcPartialNodeInfo(lastNode->nodeInfo_,this,numberChangedBounds,
1141                               variables,boundChanges,basisDiff) ;
1142    else
1143      nodeInfo_ =
1144        strategy->partialNodeInfo(model,lastNode->nodeInfo_,this,
1145                                  numberChangedBounds,variables,boundChanges,
1146                                  basisDiff) ;
1147    delete basisDiff ;
1148    delete [] boundChanges;
1149    delete [] variables;
1150    delete expanded ;
1151    if  (mustDeleteBasis)
1152      delete ws;
1153  }
1154  // Set node number
1155  nodeInfo_->setNodeNumber(model->getNodeCount2());
1156  state_ |= 2; // say active
1157}
1158
1159#else   // CBC_NEW_CREATEINFO
1160
1161/*
1162  Original createInfo, with bare manipulation of basis vectors. Fails if solver
1163  maintains additional information in basis.
1164*/
1165
1166void
1167CbcNode::createInfo (CbcModel *model,
1168                     CbcNode *lastNode,
1169                     const CoinWarmStartBasis *lastws,
1170                     const double *lastLower, const double *lastUpper,
1171                     int numberOldActiveCuts,int numberNewCuts)
1172{ OsiSolverInterface * solver = model->solver();
1173  CbcStrategy * strategy = model->strategy();
1174  /*
1175    The root --- no parent. Create full basis and bounds information.
1176  */
1177  if (!lastNode)
1178    { 
1179      if (!strategy)
1180        nodeInfo_=new CbcFullNodeInfo(model,solver->getNumRows());
1181      else
1182        nodeInfo_ = strategy->fullNodeInfo(model,solver->getNumRows());
1183    }
1184  /*
1185    Not the root. Create an edit from the parent's basis & bound information.
1186    This is not quite as straightforward as it seems. We need to reintroduce
1187    cuts we may have dropped out of the basis, in the correct position, because
1188    this whole process is strictly positional. Start by grabbing the current
1189    basis.
1190  */
1191  else
1192    { 
1193      bool mustDeleteBasis;
1194      const CoinWarmStartBasis* ws =
1195        dynamic_cast<const CoinWarmStartBasis*>(solver->getPointerToWarmStart(mustDeleteBasis));
1196      assert(ws!=NULL); // make sure not volume
1197      //int numberArtificials = lastws->getNumArtificial();
1198      int numberColumns = solver->getNumCols();
1199     
1200      const double * lower = solver->getColLower();
1201      const double * upper = solver->getColUpper();
1202     
1203      int i;
1204      /*
1205        Create a clone and resize it to hold all the structural constraints, plus
1206        all the cuts: old cuts, both active and inactive (currentNumberCuts), and
1207        new cuts (numberNewCuts).
1208       
1209        TODO: You'd think that the set of constraints (logicals) in the expanded
1210        basis should match the set represented in lastws. At least, that's
1211        what I thought. But at the point I first looked hard at this bit of
1212        code, it turned out that lastws was the stripped basis produced at
1213        the end of addCuts(), rather than the raw basis handed back by
1214        addCuts1(). The expanded basis here is equivalent to the raw basis of
1215        addCuts1(). I said ``whoa, that's not good, I must have introduced a
1216        bug'' and went back to John's code to see where I'd gone wrong.
1217        And discovered the same `error' in his code.
1218       
1219        After a bit of thought, my conclusion is that correctness is not
1220        affected by whether lastws is the stripped or raw basis. The diffs
1221        have no semantics --- just a set of changes that need to be made
1222        to convert lastws into expanded. I think the only effect is that we
1223        store a lot more diffs (everything in expanded that's not covered by
1224        the stripped basis). But I need to give this more thought. There
1225        may well be some subtle error cases.
1226       
1227        In the mean time, I've twiddled addCuts() to set lastws to the raw
1228        basis. Makes me (Lou) less nervous to compare apples to apples.
1229      */
1230      CoinWarmStartBasis *expanded = 
1231        dynamic_cast<CoinWarmStartBasis *>(ws->clone()) ;
1232      int numberRowsAtContinuous = model->numberRowsAtContinuous();
1233      int iFull = numberRowsAtContinuous+model->currentNumberCuts()+
1234        numberNewCuts;
1235      //int numberArtificialsNow = iFull;
1236      //int maxBasisLength = ((iFull+15)>>4)+((numberColumns+15)>>4);
1237      //printf("l %d full %d\n",maxBasisLength,iFull);
1238      if (expanded) 
1239        expanded->resize(iFull,numberColumns);
1240#ifdef CBC_CHECK_BASIS
1241      printf("Before expansion: orig %d, old %d, new %d, current %d\n",
1242             numberRowsAtContinuous,numberOldActiveCuts,numberNewCuts,
1243             model->currentNumberCuts()) ;
1244      ws->print();
1245#endif
1246      /*
1247        Now fill in the expanded basis. Any indices beyond nPartial must
1248        be cuts created while processing this node --- they can be copied directly
1249        into the expanded basis. From nPartial down, pull the status of active cuts
1250        from ws, interleaving with a B entry for the deactivated (loose) cuts.
1251      */
1252      int numberDropped = model->currentNumberCuts()-numberOldActiveCuts;
1253      int iCompact=iFull-numberDropped;
1254      CbcCountRowCut ** cut = model->addedCuts();
1255      int nPartial = model->currentNumberCuts()+numberRowsAtContinuous;
1256      iFull--;
1257      for (;iFull>=nPartial;iFull--) {
1258        CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
1259        //assert (status != CoinWarmStartBasis::basic); // may be permanent cut
1260        expanded->setArtifStatus(iFull,status);
1261      }
1262      for (;iFull>=numberRowsAtContinuous;iFull--) {
1263        if (cut[iFull-numberRowsAtContinuous]) {
1264          CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
1265          // If no cut generator being used then we may have basic variables
1266          //if (model->getMaximumCutPasses()&&
1267          //  status == CoinWarmStartBasis::basic)
1268          //printf("cut basic\n");
1269          expanded->setArtifStatus(iFull,status);
1270        } else {
1271          expanded->setArtifStatus(iFull,CoinWarmStartBasis::basic);
1272        }
1273      }
1274#ifdef CBC_CHECK_BASIS
1275      printf("Expanded basis\n");
1276      expanded->print() ;
1277      printf("Diffing against\n") ;
1278      lastws->print() ;
1279#endif   
1280      /*
1281        Now that we have two bases in proper positional correspondence, creating
1282        the actual diff is dead easy.
1283      */
1284     
1285      CoinWarmStartDiff *basisDiff = expanded->generateDiff(lastws) ;
1286      /*
1287        Diff the bound vectors. It's assumed the number of structural variables is
1288        not changing. Assuming that branching objects all involve integer variables,
1289        we should see at least one bound change as a consequence of processing this
1290        subproblem. Different types of branching objects could break this assertion.
1291        Not true at all - we have not applied current branch - JJF.
1292      */
1293      double *boundChanges = new double [2*numberColumns] ;
1294      int *variables = new int [2*numberColumns] ;
1295      int numberChangedBounds=0;
1296      for (i=0;i<numberColumns;i++) {
1297        if (lower[i]!=lastLower[i]) {
1298          variables[numberChangedBounds]=i;
1299          boundChanges[numberChangedBounds++]=lower[i];
1300        }
1301        if (upper[i]!=lastUpper[i]) {
1302          variables[numberChangedBounds]=i|0x80000000;
1303          boundChanges[numberChangedBounds++]=upper[i];
1304        }
1305#ifdef CBC_DEBUG
1306        if (lower[i]!=lastLower[i])
1307          printf("lower on %d changed from %g to %g\n",
1308                 i,lastLower[i],lower[i]);
1309        if (upper[i]!=lastUpper[i])
1310          printf("upper on %d changed from %g to %g\n",
1311                 i,lastUpper[i],upper[i]);
1312#endif
1313      }
1314#ifdef CBC_DEBUG
1315      printf("%d changed bounds\n",numberChangedBounds) ;
1316#endif
1317      //if (lastNode->branchingObject()->boundBranch())
1318      //assert (numberChangedBounds);
1319      /*
1320        Hand the lot over to the CbcPartialNodeInfo constructor, then clean up and
1321        return.
1322      */
1323      if (!strategy)
1324        nodeInfo_ =
1325          new CbcPartialNodeInfo(lastNode->nodeInfo_,this,numberChangedBounds,
1326                                 variables,boundChanges,basisDiff) ;
1327      else
1328        nodeInfo_ = strategy->partialNodeInfo(model, lastNode->nodeInfo_,this,numberChangedBounds,
1329                                              variables,boundChanges,basisDiff) ;
1330      delete basisDiff ;
1331      delete [] boundChanges;
1332      delete [] variables;
1333      delete expanded ;
1334      if  (mustDeleteBasis)
1335        delete ws;
1336    }
1337  // Set node number
1338  nodeInfo_->setNodeNumber(model->getNodeCount2());
1339  state_ |= 2; // say active
1340}
1341
1342#endif  // CBC_NEW_CREATEINFO
1343/*
1344  The routine scans through the object list of the model looking for objects
1345  that indicate infeasibility. It tests each object using strong branching
1346  and selects the one with the least objective degradation.  A corresponding
1347  branching object is left attached to lastNode.
1348 
1349  If strong branching is disabled, a candidate object is chosen essentially
1350  at random (whatever object ends up in pos'n 0 of the candidate array).
1351 
1352  If a branching candidate is found to be monotone, bounds are set to fix the
1353  variable and the routine immediately returns (the caller is expected to
1354  reoptimize).
1355 
1356  If a branching candidate is found to result in infeasibility in both
1357  directions, the routine immediately returns an indication of infeasibility.
1358 
1359  Returns:  0   both branch directions are feasible
1360  -1    branching variable is monotone
1361  -2    infeasible
1362 
1363  Original comments:
1364  Here could go cuts etc etc
1365  For now just fix on objective from strong branching.
1366*/
1367
1368int CbcNode::chooseBranch (CbcModel *model, CbcNode *lastNode,int numberPassesLeft)
1369 
1370{ if (lastNode)
1371    depth_ = lastNode->depth_+1;
1372  else
1373    depth_ = 0;
1374  delete branch_;
1375  branch_=NULL;
1376  OsiSolverInterface * solver = model->solver();
1377# ifdef COIN_HAS_CLP
1378  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
1379  int saveClpOptions=0;
1380  if (osiclp) {
1381    // for faster hot start
1382    saveClpOptions = osiclp->specialOptions();
1383    osiclp->setSpecialOptions(saveClpOptions|8192);
1384  }
1385# else
1386  OsiSolverInterface *osiclp = NULL ;
1387# endif
1388  double saveObjectiveValue = solver->getObjValue();
1389  double objectiveValue = CoinMax(solver->getObjSense()*saveObjectiveValue,objectiveValue_);
1390  const double * lower = solver->getColLower();
1391  const double * upper = solver->getColUpper();
1392  // See what user thinks
1393  int anyAction=model->problemFeasibility()->feasible(model,0);
1394  if (anyAction) {
1395    // will return -2 if infeasible , 0 if treat as integer
1396    return anyAction-1;
1397  }
1398  double integerTolerance = 
1399    model->getDblParam(CbcModel::CbcIntegerTolerance);
1400  // point to useful information
1401  OsiBranchingInformation usefulInfo = model->usefulInformation();
1402  // and modify
1403  usefulInfo.depth_=depth_;
1404  int i;
1405  bool beforeSolution = model->getSolutionCount()==0;
1406  int numberStrong=model->numberStrong();
1407  // switch off strong if hotstart
1408  const double * hotstartSolution = model->hotstartSolution();
1409  const int * hotstartPriorities = model->hotstartPriorities();
1410  int numberObjects = model->numberObjects();
1411  int numberColumns = model->getNumCols();
1412  double * saveUpper = new double[numberColumns];
1413  double * saveLower = new double[numberColumns];
1414  for (i=0;i<numberColumns;i++) {
1415    saveLower[i] = lower[i];
1416    saveUpper[i] = upper[i];
1417  }
1418 
1419  // Save solution in case heuristics need good solution later
1420 
1421  double * saveSolution = new double[numberColumns];
1422  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1423  model->reserveCurrentSolution(saveSolution);
1424  if (hotstartSolution) {
1425    numberStrong=0;
1426    if((model->moreSpecialOptions()&1024)!=0) {
1427      int nBad=0;
1428      int nUnsat=0;
1429      int nDiff=0;
1430      for (int i=0;i<numberObjects;i++) {
1431        OsiObject * object = model->modifiableObject(i);
1432        const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
1433        if (thisOne) {
1434          int iColumn = thisOne->columnNumber();
1435          double targetValue = hotstartSolution[iColumn];
1436          double value = saveSolution[iColumn];
1437          if (fabs(value-floor(value+0.5))>1.0e-6) {
1438            nUnsat++;
1439#ifdef CLP_INVESTIGATE
1440            printf("H %d is %g target %g\n",iColumn,value,targetValue);
1441#endif
1442          } else if (fabs(targetValue-value)>1.0e-6) {
1443            nDiff++;
1444          }
1445          if (targetValue<saveLower[iColumn]||
1446              targetValue>saveUpper[iColumn]) {
1447#ifdef CLP_INVESTIGATE
1448            printf("%d has target %g and current bounds %g and %g\n",
1449                   iColumn,targetValue,saveLower[iColumn],saveUpper[iColumn]);
1450#endif
1451            nBad++;
1452          }
1453        }
1454      }
1455#ifdef CLP_INVESTIGATE
1456      printf("Hot %d unsatisfied, %d outside limits, %d different\n",
1457             nUnsat,nBad,nDiff); 
1458#endif
1459      if (nBad) {
1460        // switch off as not possible
1461          hotstartSolution=NULL;
1462          model->setHotstartSolution(NULL,NULL);
1463          usefulInfo.hotstartSolution_=NULL;
1464      }
1465    }
1466  }
1467  int numberStrongDone=0;
1468  int numberUnfinished=0;
1469  int numberStrongInfeasible=0;
1470  int numberStrongIterations=0;
1471  int saveNumberStrong=numberStrong;
1472  bool checkFeasibility = numberObjects>model->numberIntegers();
1473  int maximumStrong = CoinMax(CoinMin(numberStrong,numberObjects),1);
1474  /*
1475    Get a branching decision object. Use the default decision criteria unless
1476    the user has loaded a decision method into the model.
1477  */
1478  CbcBranchDecision *decision = model->branchingMethod();
1479  CbcDynamicPseudoCostBranchingObject * dynamicBranchingObject =
1480    dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(decision);
1481  if (!decision||dynamicBranchingObject)
1482    decision = new CbcBranchDefaultDecision();
1483  decision->initialize(model);
1484  CbcStrongInfo * choice = new CbcStrongInfo[maximumStrong];
1485  // May go round twice if strong branching fixes all local candidates
1486  bool finished=false;
1487  double estimatedDegradation=0.0; 
1488  while(!finished) {
1489    finished=true;
1490    // Some objects may compute an estimate of best solution from here
1491    estimatedDegradation=0.0; 
1492    //int numberIntegerInfeasibilities=0; // without odd ones
1493    numberStrongDone=0;
1494    numberUnfinished=0;
1495    numberStrongInfeasible=0;
1496    numberStrongIterations=0;
1497   
1498    // We may go round this loop twice (only if we think we have solution)
1499    for (int iPass=0;iPass<2;iPass++) {
1500     
1501      // compute current state
1502      //int numberObjectInfeasibilities; // just odd ones
1503      //model->feasibleSolution(
1504      //                      numberIntegerInfeasibilities,
1505      //                      numberObjectInfeasibilities);
1506      // Some objects may compute an estimate of best solution from here
1507      estimatedDegradation=0.0; 
1508      numberUnsatisfied_ = 0;
1509      // initialize sum of "infeasibilities"
1510      sumInfeasibilities_ = 0.0;
1511      int bestPriority=COIN_INT_MAX;
1512      /*
1513        Scan for branching objects that indicate infeasibility. Choose the best
1514        maximumStrong candidates, using priority as the first criteria, then
1515        integer infeasibility.
1516       
1517        The algorithm is to fill the choice array with a set of good candidates (by
1518        infeasibility) with priority bestPriority.  Finding a candidate with
1519        priority better (less) than bestPriority flushes the choice array. (This
1520        serves as initialization when the first candidate is found.)
1521       
1522        A new candidate is added to choices only if its infeasibility exceeds the
1523        current max infeasibility (mostAway). When a candidate is added, it
1524        replaces the candidate with the smallest infeasibility (tracked by
1525        iSmallest).
1526      */
1527      int iSmallest = 0;
1528      double mostAway = 1.0e-100;
1529      for (i = 0 ; i < maximumStrong ; i++)
1530        choice[i].possibleBranch = NULL ;
1531      numberStrong=0;
1532      bool canDoOneHot=false;
1533      for (i=0;i<numberObjects;i++) {
1534        OsiObject * object = model->modifiableObject(i);
1535        int preferredWay;
1536        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
1537        int priorityLevel = object->priority();
1538        if (hotstartSolution) {
1539          // we are doing hot start
1540          const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
1541          if (thisOne) {
1542            int iColumn = thisOne->columnNumber();
1543            bool canDoThisHot=true;
1544            double targetValue = hotstartSolution[iColumn];
1545            if (saveUpper[iColumn]>saveLower[iColumn]) {
1546              double value = saveSolution[iColumn];
1547              if (hotstartPriorities)
1548                priorityLevel=hotstartPriorities[iColumn]; 
1549              //double originalLower = thisOne->originalLower();
1550              //double originalUpper = thisOne->originalUpper();
1551              // switch off if not possible
1552              if (targetValue>=saveLower[iColumn]&&targetValue<=saveUpper[iColumn]) {
1553                /* priority outranks rest always if negative
1554                   otherwise can be downgraded if at correct level.
1555                   Infeasibility may be increased to choose 1.0 values first.
1556                   choose one near wanted value
1557                */
1558                if (fabs(value-targetValue)>integerTolerance) {
1559                  //if (infeasibility>0.01)
1560                  //infeasibility = fabs(1.0e6-fabs(value-targetValue));
1561                  //else
1562                  infeasibility = fabs(value-targetValue);
1563                  //if (targetValue==1.0)
1564                  //infeasibility += 1.0;
1565                  if (value>targetValue) {
1566                    preferredWay=-1;
1567                  } else {
1568                    preferredWay=1;
1569                  }
1570                  priorityLevel = CoinAbs(priorityLevel);
1571                } else if (priorityLevel<0) {
1572                  priorityLevel = CoinAbs(priorityLevel);
1573                  if (targetValue==saveLower[iColumn]) {
1574                    infeasibility = integerTolerance+1.0e-12;
1575                    preferredWay=-1;
1576                  } else if (targetValue==saveUpper[iColumn]) {
1577                    infeasibility = integerTolerance+1.0e-12;
1578                    preferredWay=1;
1579                  } else {
1580                    // can't
1581                    priorityLevel += 10000000;
1582                    canDoThisHot=false;
1583                  }
1584                } else {
1585                  priorityLevel += 10000000;
1586                  canDoThisHot=false;
1587                }
1588              } else {
1589                // switch off if not possible
1590                canDoThisHot=false;
1591              }
1592              if (canDoThisHot)
1593                canDoOneHot=true;
1594            } else if (targetValue<saveLower[iColumn]||targetValue>saveUpper[iColumn]) {
1595            }
1596          } else {
1597            priorityLevel += 10000000;
1598          }
1599        }
1600        if (infeasibility) {
1601          // Increase estimated degradation to solution
1602          estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate());
1603          numberUnsatisfied_++;
1604          sumInfeasibilities_ += infeasibility;
1605          // Better priority? Flush choices.
1606          if (priorityLevel<bestPriority) {
1607            int j;
1608            iSmallest=0;
1609            for (j=0;j<maximumStrong;j++) {
1610              choice[j].upMovement=0.0;
1611              delete choice[j].possibleBranch;
1612              choice[j].possibleBranch=NULL;
1613            }
1614            bestPriority = priorityLevel;
1615            mostAway=1.0e-100;
1616            numberStrong=0;
1617          } else if (priorityLevel>bestPriority) {
1618            continue;
1619          }
1620          // Check for suitability based on infeasibility.
1621          if (infeasibility>mostAway) {
1622            //add to list
1623            choice[iSmallest].upMovement=infeasibility;
1624            delete choice[iSmallest].possibleBranch;
1625            CbcObject * obj =
1626              dynamic_cast <CbcObject *>(object) ;
1627            assert (obj);
1628            choice[iSmallest].possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
1629            numberStrong = CoinMax(numberStrong,iSmallest+1);
1630            // Save which object it was
1631            choice[iSmallest].objectNumber=i;
1632            int j;
1633            iSmallest=-1;
1634            mostAway = 1.0e50;
1635            for (j=0;j<maximumStrong;j++) {
1636              if (choice[j].upMovement<mostAway) {
1637                mostAway=choice[j].upMovement;
1638                iSmallest=j;
1639              }
1640            }
1641          }
1642        }
1643      }
1644      if (!canDoOneHot&&hotstartSolution) {
1645        // switch off as not possible
1646        hotstartSolution=NULL;
1647        model->setHotstartSolution(NULL,NULL);
1648        usefulInfo.hotstartSolution_=NULL;
1649      }
1650      if (numberUnsatisfied_) {
1651        // some infeasibilities - go to next steps
1652#ifdef CLP_INVESTIGATE
1653        if (hotstartSolution) {
1654          int k=choice[0].objectNumber;
1655          OsiObject * object = model->modifiableObject(k);
1656          const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
1657          assert (thisOne);
1658          int iColumn = thisOne->columnNumber();
1659          double targetValue = hotstartSolution[iColumn];
1660          double value = saveSolution[iColumn];
1661          printf("Branch on %d has target %g (value %g) and current bounds %g and %g\n",
1662                 iColumn,targetValue,value,saveLower[iColumn],saveUpper[iColumn]);
1663        }
1664#endif
1665        break;
1666      } else if (!iPass) {
1667        // looks like a solution - get paranoid
1668        bool roundAgain=false;
1669        // get basis
1670        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
1671        if (!ws)
1672          break;
1673        for (i=0;i<numberColumns;i++) {
1674          double value = saveSolution[i];
1675          if (value<lower[i]) {
1676            saveSolution[i]=lower[i];
1677            roundAgain=true;
1678            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
1679          } else if (value>upper[i]) {
1680            saveSolution[i]=upper[i];
1681            roundAgain=true;
1682            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
1683          } 
1684        }
1685        if (roundAgain&&saveNumberStrong) {
1686          // restore basis
1687          solver->setWarmStart(ws);
1688          delete ws;
1689          solver->resolve();
1690          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1691          model->reserveCurrentSolution(saveSolution);
1692          if (!solver->isProvenOptimal()) {
1693            // infeasible
1694            anyAction=-2;
1695            break;
1696          }
1697        } else {
1698          delete ws;
1699          break;
1700        }
1701      }
1702    }
1703    /* Some solvers can do the strong branching calculations faster if
1704       they do them all at once.  At present only Clp does for ordinary
1705       integers but I think this coding would be easy to modify
1706    */
1707    bool allNormal=true; // to say if we can do fast strong branching
1708    // Say which one will be best
1709    int bestChoice=0;
1710    double worstInfeasibility=0.0;
1711    for (i=0;i<numberStrong;i++) {
1712      choice[i].numIntInfeasUp = numberUnsatisfied_;
1713      choice[i].numIntInfeasDown = numberUnsatisfied_;
1714      choice[i].fix=0; // say not fixed
1715      if (!dynamic_cast <const CbcSimpleInteger *> (model->object(choice[i].objectNumber)))
1716        allNormal=false; // Something odd so lets skip clever fast branching
1717      if ( !model->object(choice[i].objectNumber)->boundBranch())
1718        numberStrong=0; // switch off
1719      if ( choice[i].possibleBranch->numberBranches()>2)
1720        numberStrong=0; // switch off
1721      // Do best choice in case switched off
1722      if (choice[i].upMovement>worstInfeasibility) {
1723        worstInfeasibility=choice[i].upMovement;
1724        bestChoice=i;
1725      }
1726    }
1727    // If we have hit max time don't do strong branching
1728    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
1729                        model->getDblParam(CbcModel::CbcMaximumSeconds));
1730    // also give up if we are looping round too much
1731    if (hitMaxTime||numberPassesLeft<=0)
1732      numberStrong=0;
1733    /*
1734      Is strong branching enabled? If so, set up and do it. Otherwise, we'll
1735      fall through to simple branching.
1736     
1737      Setup for strong branching involves saving the current basis (for restoration
1738      afterwards) and setting up for hot starts.
1739    */
1740    if (numberStrong&&saveNumberStrong) {
1741     
1742      bool solveAll=false; // set true to say look at all even if some fixed (experiment)
1743      solveAll=true;
1744      // worth trying if too many times
1745      // Save basis
1746      CoinWarmStart * ws = solver->getWarmStart();
1747      // save limit
1748      int saveLimit;
1749      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
1750      if (beforeSolution&&saveLimit<100)
1751        solver->setIntParam(OsiMaxNumIterationHotStart,100); // go to end
1752#     ifdef COIN_HAS_CLP     
1753      /* If we are doing all strong branching in one go then we create new arrays
1754         to store information.  If clp NULL then doing old way.
1755         Going down -
1756         outputSolution[2*i] is final solution.
1757         outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown
1758         outputStuff[2*i+numberStrong] is number iterations
1759         On entry newUpper[i] is new upper bound, on exit obj change
1760         Going up -
1761         outputSolution[2*i+1] is final solution.
1762         outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
1763         outputStuff[2*i+1+numberStrong] is number iterations
1764         On entry newLower[i] is new lower bound, on exit obj change
1765      */
1766      ClpSimplex * clp=NULL;
1767      double * newLower = NULL;
1768      double * newUpper = NULL;
1769      double ** outputSolution=NULL;
1770      int * outputStuff=NULL;
1771      // Go back to normal way if user wants it
1772      if (osiclp&&(osiclp->specialOptions()&16)!=0&&osiclp->specialOptions()>0)
1773        allNormal=false;
1774      if (osiclp&&!allNormal) {
1775        // say do fast
1776        int easy=1;
1777        osiclp->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
1778      }
1779      if (osiclp&& allNormal) {
1780        clp = osiclp->getModelPtr();
1781        // Clp - do a different way
1782        newLower = new double[numberStrong];
1783        newUpper = new double[numberStrong];
1784        outputSolution = new double * [2*numberStrong];
1785        outputStuff = new int [4*numberStrong];
1786        int * which = new int[numberStrong];
1787        int startFinishOptions;
1788        int specialOptions = osiclp->specialOptions();
1789        int clpOptions = clp->specialOptions();
1790        int returnCode=0;
1791#define CRUNCH
1792#ifdef CRUNCH
1793        // Crunch down problem
1794        int numberRows = clp->numberRows();
1795        // Use dual region
1796        double * rhs = clp->dualRowSolution();
1797        int * whichRow = new int[3*numberRows];
1798        int * whichColumn = new int[2*numberColumns];
1799        int nBound;
1800        ClpSimplex * small = static_cast<ClpSimplexOther *> (clp)->crunch(rhs,whichRow,whichColumn,nBound,true);
1801        if (!small) {
1802          anyAction=-2;
1803          //printf("XXXX Inf by inspection\n");
1804          delete [] whichColumn;
1805          whichColumn=NULL;
1806          delete [] whichRow;
1807          whichRow=NULL;
1808          break;
1809        } else {
1810          clp = small;
1811        }
1812#else
1813        int saveLogLevel = clp->logLevel();
1814        int saveMaxIts = clp->maximumIterations();
1815#endif
1816        clp->setLogLevel(0);
1817        if((specialOptions&1)==0) {
1818          startFinishOptions=0;
1819          clp->setSpecialOptions(clpOptions|(64|1024));
1820        } else {
1821          startFinishOptions=1+2+4;
1822          //startFinishOptions=1+4; // for moment re-factorize
1823          if((specialOptions&4)==0) 
1824            clp->setSpecialOptions(clpOptions|(64|128|512|1024|4096));
1825          else
1826            clp->setSpecialOptions(clpOptions|(64|128|512|1024|2048|4096));
1827        }
1828        // User may want to clean up before strong branching
1829        if ((clp->specialOptions()&32)!=0) {
1830          clp->primal(1);
1831          if (clp->numberIterations())
1832            model->messageHandler()->message(CBC_ITERATE_STRONG,*model->messagesPointer())
1833              << clp->numberIterations()
1834              <<CoinMessageEol;
1835        }
1836        clp->setMaximumIterations(saveLimit);
1837#ifdef CRUNCH
1838        int * backColumn = whichColumn+numberColumns;
1839#endif
1840        for (i=0;i<numberStrong;i++) {
1841          int iObject = choice[i].objectNumber;
1842          const OsiObject * object = model->object(iObject);
1843          const CbcSimpleInteger * simple = static_cast <const CbcSimpleInteger *> (object);
1844          int iSequence = simple->columnNumber();
1845          newLower[i]= ceil(saveSolution[iSequence]);
1846          newUpper[i]= floor(saveSolution[iSequence]);
1847#ifdef CRUNCH
1848          iSequence = backColumn[iSequence];
1849          assert (iSequence>=0);
1850#endif
1851          which[i]=iSequence;
1852          outputSolution[2*i]= new double [numberColumns];
1853          outputSolution[2*i+1]= new double [numberColumns];
1854        }
1855        //clp->writeMps("bad");
1856        returnCode=clp->strongBranching(numberStrong,which,
1857                                        newLower, newUpper,outputSolution,
1858                                        outputStuff,outputStuff+2*numberStrong,!solveAll,false,
1859                                        startFinishOptions);
1860#ifndef CRUNCH
1861        clp->setSpecialOptions(clpOptions); // restore
1862        clp->setMaximumIterations(saveMaxIts);
1863        clp->setLogLevel(saveLogLevel);
1864#endif
1865        if (returnCode==-2) {
1866          // bad factorization!!!
1867          // Doing normal way
1868          // Mark hot start
1869          solver->markHotStart();
1870          clp = NULL;
1871        } else {
1872#ifdef CRUNCH
1873          // extract solution
1874          //bool checkSol=true;
1875          for (i=0;i<numberStrong;i++) {
1876            int iObject = choice[i].objectNumber;
1877            const OsiObject * object = model->object(iObject);
1878            const CbcSimpleInteger * simple = static_cast <const CbcSimpleInteger *> (object);
1879            int iSequence = simple->columnNumber();
1880            which[i]=iSequence;
1881            double * sol = outputSolution[2*i];
1882            double * sol2 = outputSolution[2*i+1];
1883            //bool x=true;
1884            //bool x2=true;
1885            for (int iColumn=numberColumns-1;iColumn>=0;iColumn--) {
1886              int jColumn = backColumn[iColumn];
1887              if (jColumn>=0) {
1888                sol[iColumn]=sol[jColumn];
1889                sol2[iColumn]=sol2[jColumn];
1890              } else {
1891                sol[iColumn]=saveSolution[iColumn];
1892                sol2[iColumn]=saveSolution[iColumn];
1893              }
1894            }
1895          }
1896#endif
1897        }
1898#ifdef CRUNCH
1899        delete [] whichColumn;
1900        delete [] whichRow;
1901        delete small;
1902#endif
1903        delete [] which;
1904      } else {
1905        // Doing normal way
1906        // Mark hot start
1907        solver->markHotStart();
1908      }
1909#     else      /* COIN_HAS_CLP */
1910     
1911      OsiSolverInterface *clp = NULL ;
1912      double **outputSolution = NULL ;
1913      int *outputStuff = NULL ;
1914      double * newLower = NULL ;
1915      double * newUpper = NULL ;
1916     
1917      solver->markHotStart();
1918     
1919#     endif     /* COIN_HAS_CLP */
1920      /*
1921        Open a loop to do the strong branching LPs. For each candidate variable,
1922        solve an LP with the variable forced down, then up. If a direction turns
1923        out to be infeasible or monotonic (i.e., over the dual objective cutoff),
1924        force the objective change to be big (1.0e100). If we determine the problem
1925        is infeasible, or find a monotone variable, escape the loop.
1926       
1927        TODO: The `restore bounds' part might be better encapsulated as an
1928        unbranch() method. Branching objects more exotic than simple integers
1929        or cliques might not restrict themselves to variable bounds.
1930       
1931        TODO: Virtuous solvers invalidate the current solution (or give bogus
1932        results :-) when the bounds are changed out from under them. So we
1933        need to do all the work associated with finding a new solution before
1934        restoring the bounds.
1935      */
1936      for (i = 0 ; i < numberStrong ; i++)
1937        { double objectiveChange ;
1938          double newObjectiveValue=1.0e100;
1939          // status is 0 finished, 1 infeasible and other
1940          int iStatus;
1941          /*
1942            Try the down direction first. (Specify the initial branching alternative as
1943            down with a call to way(-1). Each subsequent call to branch() performs the
1944            specified branch and advances the branch object state to the next branch
1945            alternative.)
1946          */
1947          if (!clp) {
1948            choice[i].possibleBranch->way(-1) ;
1949            choice[i].possibleBranch->branch() ;
1950            bool feasible=true;
1951            if (checkFeasibility) {
1952              // check branching did not make infeasible
1953              int iColumn;
1954              int numberColumns = solver->getNumCols();
1955              const double * columnLower = solver->getColLower();
1956              const double * columnUpper = solver->getColUpper();
1957              for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1958                if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1959                  feasible=false;
1960              }
1961            }
1962            if (feasible) {
1963              solver->solveFromHotStart() ;
1964              numberStrongDone++;
1965              numberStrongIterations += solver->getIterationCount();
1966              /*
1967                We now have an estimate of objective degradation that we can use for strong
1968                branching. If we're over the cutoff, the variable is monotone up.
1969                If we actually made it to optimality, check for a solution, and if we have
1970                a good one, call setBestSolution to process it. Note that this may reduce the
1971                cutoff, so we check again to see if we can declare this variable monotone.
1972              */
1973              if (solver->isProvenOptimal())
1974                iStatus=0; // optimal
1975              else if (solver->isIterationLimitReached()
1976                       &&!solver->isDualObjectiveLimitReached())
1977                iStatus=2; // unknown
1978              else
1979                iStatus=1; // infeasible
1980              newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1981              choice[i].numItersDown = solver->getIterationCount();
1982            } else {
1983              iStatus=1; // infeasible
1984              newObjectiveValue = 1.0e100;
1985              choice[i].numItersDown = 0;
1986            }
1987          } else {
1988            iStatus = outputStuff[2*i];
1989            choice[i].numItersDown = outputStuff[2*numberStrong+2*i];
1990            numberStrongDone++;
1991            numberStrongIterations += choice[i].numItersDown;
1992            newObjectiveValue = objectiveValue+newUpper[i];
1993            solver->setColSolution(outputSolution[2*i]);
1994          }
1995          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
1996          if (!iStatus) {
1997            choice[i].finishedDown = true ;
1998            if (newObjectiveValue>=model->getCutoff()) {
1999              objectiveChange = 1.0e100; // say infeasible
2000              numberStrongInfeasible++;
2001            } else {
2002              // See if integer solution
2003              if (model->feasibleSolution(choice[i].numIntInfeasDown,
2004                                          choice[i].numObjInfeasDown)
2005                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
2006                model->setBestSolution(CBC_STRONGSOL,
2007                                       newObjectiveValue,
2008                                       solver->getColSolution()) ;
2009                // only needed for odd solvers
2010                newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2011                objectiveChange = CoinMax(newObjectiveValue-objectiveValue_,0.0) ;
2012                model->setLastHeuristic(NULL);
2013                model->incrementUsed(solver->getColSolution());
2014                if (newObjectiveValue >= model->getCutoff()) {  //  *new* cutoff
2015                  objectiveChange = 1.0e100 ;
2016                  numberStrongInfeasible++;
2017                }
2018              }
2019            }
2020          } else if (iStatus==1) {
2021            objectiveChange = 1.0e100 ;
2022            numberStrongInfeasible++;
2023          } else {
2024            // Can't say much as we did not finish
2025            choice[i].finishedDown = false ;
2026            numberUnfinished++;
2027          }
2028          choice[i].downMovement = objectiveChange ;
2029         
2030          // restore bounds
2031          if (!clp)
2032            { for (int j=0;j<numberColumns;j++) {
2033                if (saveLower[j] != lower[j])
2034                  solver->setColLower(j,saveLower[j]);
2035                if (saveUpper[j] != upper[j])
2036                  solver->setColUpper(j,saveUpper[j]);
2037              }
2038            }
2039          //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2040          //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersDown,
2041          //     choice[i].downMovement,choice[i].finishedDown,choice[i].numIntInfeasDown,
2042          //     choice[i].numObjInfeasDown);
2043         
2044          // repeat the whole exercise, forcing the variable up
2045          if (!clp) {
2046            bool feasible=true;
2047            // If odd branching then maybe just one possibility
2048            if(choice[i].possibleBranch->numberBranchesLeft()>0) {
2049              choice[i].possibleBranch->branch();
2050              if (checkFeasibility) {
2051                // check branching did not make infeasible
2052                int iColumn;
2053                int numberColumns = solver->getNumCols();
2054                const double * columnLower = solver->getColLower();
2055                const double * columnUpper = solver->getColUpper();
2056                for (iColumn= 0;iColumn<numberColumns;iColumn++) {
2057                  if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
2058                    feasible=false;
2059                }
2060              }
2061            } else {
2062              // second branch infeasible
2063              feasible=false;
2064            }
2065            if (feasible) {
2066              solver->solveFromHotStart() ;
2067              numberStrongDone++;
2068              numberStrongIterations += solver->getIterationCount();
2069              /*
2070                We now have an estimate of objective degradation that we can use for strong
2071                branching. If we're over the cutoff, the variable is monotone up.
2072                If we actually made it to optimality, check for a solution, and if we have
2073                a good one, call setBestSolution to process it. Note that this may reduce the
2074                cutoff, so we check again to see if we can declare this variable monotone.
2075              */
2076              if (solver->isProvenOptimal())
2077                iStatus=0; // optimal
2078              else if (solver->isIterationLimitReached()
2079                       &&!solver->isDualObjectiveLimitReached())
2080                iStatus=2; // unknown
2081              else
2082                iStatus=1; // infeasible
2083              newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2084              choice[i].numItersUp = solver->getIterationCount();
2085            } else {
2086              iStatus=1; // infeasible
2087              newObjectiveValue = 1.0e100;
2088              choice[i].numItersDown = 0;
2089            }
2090          } else {
2091            iStatus = outputStuff[2*i+1];
2092            choice[i].numItersUp = outputStuff[2*numberStrong+2*i+1];
2093            numberStrongDone++;
2094            numberStrongIterations += choice[i].numItersUp;
2095            newObjectiveValue = objectiveValue+newLower[i];
2096            solver->setColSolution(outputSolution[2*i+1]);
2097          }
2098          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
2099          if (!iStatus) {
2100            choice[i].finishedUp = true ;
2101            if (newObjectiveValue>=model->getCutoff()) {
2102              objectiveChange = 1.0e100; // say infeasible
2103              numberStrongInfeasible++;
2104            } else {
2105              // See if integer solution
2106              if (model->feasibleSolution(choice[i].numIntInfeasUp,
2107                                          choice[i].numObjInfeasUp)
2108                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
2109                model->setBestSolution(CBC_STRONGSOL,
2110                                       newObjectiveValue,
2111                                       solver->getColSolution()) ;
2112                // only needed for odd solvers
2113                newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2114                objectiveChange = CoinMax(newObjectiveValue-objectiveValue_,0.0) ;
2115                model->setLastHeuristic(NULL);
2116                model->incrementUsed(solver->getColSolution());
2117                if (newObjectiveValue >= model->getCutoff()) {  //  *new* cutoff
2118                  objectiveChange = 1.0e100 ;
2119                  numberStrongInfeasible++;
2120                }
2121              }
2122            }
2123          } else if (iStatus==1) {
2124            objectiveChange = 1.0e100 ;
2125            numberStrongInfeasible++;
2126          } else {
2127            // Can't say much as we did not finish
2128            choice[i].finishedUp = false ;
2129            numberUnfinished++;
2130          }
2131          choice[i].upMovement = objectiveChange ;
2132         
2133          // restore bounds
2134          if (!clp)
2135            { for (int j=0;j<numberColumns;j++) {
2136                if (saveLower[j] != lower[j])
2137                  solver->setColLower(j,saveLower[j]);
2138                if (saveUpper[j] != upper[j])
2139                  solver->setColUpper(j,saveUpper[j]);
2140              }
2141            }
2142         
2143          //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2144          //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersUp,
2145          //     choice[i].upMovement,choice[i].finishedUp,choice[i].numIntInfeasUp,
2146          //     choice[i].numObjInfeasUp);
2147         
2148          /*
2149            End of evaluation for this candidate variable. Possibilities are:
2150            * Both sides below cutoff; this variable is a candidate for branching.
2151            * Both sides infeasible or above the objective cutoff: no further action
2152            here. Break from the evaluation loop and assume the node will be purged
2153            by the caller.
2154            * One side below cutoff: Install the branch (i.e., fix the variable). Break
2155            from the evaluation loop and assume the node will be reoptimised by the
2156            caller.
2157          */
2158          // reset
2159          choice[i].possibleBranch->resetNumberBranchesLeft();
2160          if (choice[i].upMovement<1.0e100) {
2161            if(choice[i].downMovement<1.0e100) {
2162              // feasible - no action
2163            } else {
2164              // up feasible, down infeasible
2165              anyAction=-1;
2166              //printf("Down infeasible for choice %d sequence %d\n",i,
2167              // model->object(choice[i].objectNumber)->columnNumber());
2168              if (!solveAll) {
2169                choice[i].possibleBranch->way(1);
2170                choice[i].possibleBranch->branch();
2171                break;
2172              } else {
2173                choice[i].fix=1;
2174              }
2175            }
2176          } else {
2177            if(choice[i].downMovement<1.0e100) {
2178              // down feasible, up infeasible
2179              anyAction=-1;
2180              //printf("Up infeasible for choice %d sequence %d\n",i,
2181              // model->object(choice[i].objectNumber)->columnNumber());
2182              if (!solveAll) {
2183                choice[i].possibleBranch->way(-1);
2184                choice[i].possibleBranch->branch();
2185                break;
2186              } else {
2187                choice[i].fix=-1;
2188              }
2189            } else {
2190              // neither side feasible
2191              anyAction=-2;
2192              //printf("Both infeasible for choice %d sequence %d\n",i,
2193              // model->object(choice[i].objectNumber)->columnNumber());
2194              break;
2195            }
2196          }
2197          bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
2198                              model->getDblParam(CbcModel::CbcMaximumSeconds));
2199          if (hitMaxTime) {
2200            numberStrong=i+1;
2201            break;
2202          }
2203        }
2204      if (!clp) {
2205        // Delete the snapshot
2206        solver->unmarkHotStart();
2207      } else {
2208        delete [] newLower;
2209        delete [] newUpper;
2210        delete [] outputStuff;
2211        int i;
2212        for (i=0;i<2*numberStrong;i++)
2213          delete [] outputSolution[i];
2214        delete [] outputSolution;
2215      }
2216      solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
2217      // restore basis
2218      solver->setWarmStart(ws);
2219      // Unless infeasible we will carry on
2220      // But we could fix anyway
2221      if (anyAction==-1&&solveAll) {
2222        // apply and take off
2223        for (i = 0 ; i < numberStrong ; i++) {
2224          if (choice[i].fix) {
2225            choice[i].possibleBranch->way(choice[i].fix) ;
2226            choice[i].possibleBranch->branch() ;
2227          }
2228        }
2229        bool feasible=true;
2230        if (checkFeasibility) {
2231          // check branching did not make infeasible
2232          int iColumn;
2233          int numberColumns = solver->getNumCols();
2234          const double * columnLower = solver->getColLower();
2235          const double * columnUpper = solver->getColUpper();
2236          for (iColumn= 0;iColumn<numberColumns;iColumn++) {
2237            if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
2238              feasible=false;
2239          }
2240        }
2241        if (feasible) {
2242          // can do quick optimality check
2243          int easy=2;
2244          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
2245          solver->resolve() ;
2246          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
2247          feasible = solver->isProvenOptimal();
2248        }
2249        if (feasible) {
2250          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
2251          model->reserveCurrentSolution(saveSolution);
2252          memcpy(saveLower,solver->getColLower(),numberColumns*sizeof(double));
2253          memcpy(saveUpper,solver->getColUpper(),numberColumns*sizeof(double));
2254          // Clean up all candidates whih are fixed
2255          int numberLeft=0;
2256          for (i = 0 ; i < numberStrong ; i++) {
2257            CbcStrongInfo thisChoice = choice[i];
2258            choice[i].possibleBranch=NULL;
2259            const OsiObject * object = model->object(thisChoice.objectNumber);
2260            int preferredWay;
2261            double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
2262            if (!infeasibility) {
2263              // take out
2264              delete thisChoice.possibleBranch;
2265            } else {
2266              choice[numberLeft++]=thisChoice;
2267            }
2268          }
2269          numberStrong=numberLeft;
2270          for (;i<maximumStrong;i++) {
2271            delete choice[i].possibleBranch;
2272            choice[i].possibleBranch=NULL;
2273          }
2274          // If all fixed then round again
2275          if (!numberLeft) {
2276            finished=false;
2277            numberStrong=0;
2278            saveNumberStrong=0;
2279            maximumStrong=1;
2280          } else {
2281            anyAction=0;
2282          }
2283          // If these two uncommented then different action
2284          anyAction=-1;
2285          finished=true;
2286          //printf("some fixed but continuing %d left\n",numberLeft);
2287        } else {
2288          anyAction=-2; // say infeasible
2289        }
2290      }
2291      delete ws;
2292      int numberNodes = model->getNodeCount();
2293      // update number of strong iterations etc
2294      model->incrementStrongInfo(numberStrongDone,numberStrongIterations,
2295                                 anyAction==-2 ? 0:numberStrongInfeasible,anyAction==-2);
2296     
2297      /*
2298        anyAction >= 0 indicates that strong branching didn't produce any monotone
2299        variables. Sift through the candidates for the best one.
2300       
2301        QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is
2302        local to this code block. Perhaps should be numberNodes_ from model?
2303        Unclear what this calculation is doing.
2304      */
2305      if (anyAction>=0) {
2306       
2307        // get average cost per iteration and assume stopped ones
2308        // would stop after 50% more iterations at average cost??? !!! ???
2309        double averageCostPerIteration=0.0;
2310        double totalNumberIterations=1.0;
2311        int smallestNumberInfeasibilities=COIN_INT_MAX;
2312        for (i=0;i<numberStrong;i++) {
2313          totalNumberIterations += choice[i].numItersDown +
2314            choice[i].numItersUp ;
2315          averageCostPerIteration += choice[i].downMovement +
2316            choice[i].upMovement;
2317          smallestNumberInfeasibilities= 
2318            CoinMin(CoinMin(choice[i].numIntInfeasDown ,
2319                            choice[i].numIntInfeasUp ),
2320                    smallestNumberInfeasibilities);
2321        }
2322        //if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities)
2323        //numberNodes=1000000; // switch off search for better solution
2324        numberNodes=1000000; // switch off anyway
2325        averageCostPerIteration /= totalNumberIterations;
2326        // all feasible - choose best bet
2327       
2328        // New method does all at once so it can be more sophisticated
2329        // in deciding how to balance actions.
2330        // But it does need arrays
2331        double * changeUp = new double [numberStrong];
2332        int * numberInfeasibilitiesUp = new int [numberStrong];
2333        double * changeDown = new double [numberStrong];
2334        int * numberInfeasibilitiesDown = new int [numberStrong];
2335        CbcBranchingObject ** objects = new CbcBranchingObject * [ numberStrong];
2336        for (i = 0 ; i < numberStrong ; i++) {
2337          int iColumn = choice[i].possibleBranch->variable() ;
2338          model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
2339            << i << iColumn
2340            <<choice[i].downMovement<<choice[i].numIntInfeasDown
2341            <<choice[i].upMovement<<choice[i].numIntInfeasUp
2342            <<choice[i].possibleBranch->value()
2343            <<CoinMessageEol;
2344          changeUp[i]=choice[i].upMovement;
2345          numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
2346          changeDown[i]=choice[i].downMovement;
2347          numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
2348          objects[i] = choice[i].possibleBranch;
2349        }
2350        int whichObject = decision->bestBranch(objects,numberStrong,numberUnsatisfied_,
2351                                               changeUp,numberInfeasibilitiesUp,
2352                                               changeDown,numberInfeasibilitiesDown,
2353                                               objectiveValue_);
2354        // move branching object and make sure it will not be deleted
2355        if (whichObject>=0) {
2356          branch_ = objects[whichObject];
2357          if (model->messageHandler()->logLevel()>3) 
2358            printf("Choosing column %d\n",choice[whichObject].possibleBranch->variable()) ;
2359          choice[whichObject].possibleBranch=NULL;
2360        }
2361        delete [] changeUp;
2362        delete [] numberInfeasibilitiesUp;
2363        delete [] changeDown;
2364        delete [] numberInfeasibilitiesDown;
2365        delete [] objects;
2366      }
2367#     ifdef COIN_HAS_CLP
2368      if (osiclp&&!allNormal) {
2369        // back to normal
2370        osiclp->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
2371      }
2372#     endif 
2373    }
2374    /*
2375      Simple branching. Probably just one, but we may have got here
2376      because of an odd branch e.g. a cut
2377    */
2378    else {
2379      // not strong
2380      // C) create branching object
2381      branch_ = choice[bestChoice].possibleBranch;
2382      choice[bestChoice].possibleBranch=NULL;
2383    }
2384  }
2385  // Set guessed solution value
2386  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
2387  /*
2388    Cleanup, then we're outta here.
2389  */
2390  if (!model->branchingMethod()||dynamicBranchingObject)
2391    delete decision;
2392 
2393  for (i=0;i<maximumStrong;i++)
2394    delete choice[i].possibleBranch;
2395  delete [] choice;
2396  delete [] saveLower;
2397  delete [] saveUpper;
2398 
2399  // restore solution
2400  solver->setColSolution(saveSolution);
2401  delete [] saveSolution;
2402# ifdef COIN_HAS_CLP
2403  if (osiclp) 
2404    osiclp->setSpecialOptions(saveClpOptions);
2405# endif
2406  return anyAction;
2407}
2408
2409/*
2410  Version for dynamic pseudo costs.
2411 
2412  **** For now just return if anything odd
2413  later allow even if odd
2414 
2415  The routine scans through the object list of the model looking for objects
2416  that indicate infeasibility. It tests each object using strong branching
2417  and selects the one with the least objective degradation.  A corresponding
2418  branching object is left attached to lastNode.
2419  This version gives preference in evaluation to variables which
2420  have not been evaluated many times.  It also uses numberStrong
2421  to say give up if last few tries have not changed incumbent.
2422  See Achterberg, Koch and Martin.
2423 
2424  If strong branching is disabled, a candidate object is chosen essentially
2425  at random (whatever object ends up in pos'n 0 of the candidate array).
2426 
2427  If a branching candidate is found to be monotone, bounds are set to fix the
2428  variable and the routine immediately returns (the caller is expected to
2429  reoptimize).
2430 
2431  If a branching candidate is found to result in infeasibility in both
2432  directions, the routine immediately returns an indication of infeasibility.
2433 
2434  Returns:  0   both branch directions are feasible
2435  -1    branching variable is monotone
2436  -2    infeasible
2437  -3   Use another method
2438 
2439  For now just fix on objective from strong branching.
2440*/
2441
2442int CbcNode::chooseDynamicBranch (CbcModel *model, CbcNode *lastNode,
2443                                  OsiSolverBranch * & /*branches*/,
2444                                  int numberPassesLeft)
2445 
2446{ if (lastNode)
2447    depth_ = lastNode->depth_+1;
2448  else
2449    depth_ = 0;
2450  // Go to other choose if hot start
2451  if (model->hotstartSolution()&&
2452      (((model->moreSpecialOptions()&1024)==0)||false))
2453    return -3;
2454  delete branch_;
2455  branch_=NULL;
2456  OsiSolverInterface * solver = model->solver();
2457  // get information on solver type
2458  const OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
2459  const OsiBabSolver * auxiliaryInfo = dynamic_cast<const OsiBabSolver *> (auxInfo);
2460  if (!auxiliaryInfo) {
2461    // use one from CbcModel
2462    auxiliaryInfo = model->solverCharacteristics();
2463  }
2464  int numberObjects = model->numberObjects();
2465  // If very odd set of objects then use older chooseBranch
2466  bool useOldWay = false;
2467  // point to useful information
2468  OsiBranchingInformation usefulInfo = model->usefulInformation();
2469  if (numberObjects>model->numberIntegers()) {
2470    for (int i=model->numberIntegers();i<numberObjects;i++) {
2471      OsiObject * object = model->modifiableObject(i);
2472      CbcObject * obj = dynamic_cast <CbcObject *>(object) ;
2473      if (!obj || !obj->optionalObject()) {
2474        int preferredWay;
2475        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
2476        if (infeasibility) {
2477          useOldWay=true;
2478          break;
2479        }
2480      }
2481    }
2482  }
2483  if ((model->specialOptions()&128)!=0)
2484    useOldWay=false; // allow
2485  // For now return if not simple
2486  if (useOldWay)
2487    return -3;
2488  // Modify useful info
2489  usefulInfo.depth_=depth_;
2490  if ((model->specialOptions()&128)!=0) {
2491    // SOS - shadow prices
2492    int numberRows = solver->getNumRows();
2493    const double * pi = usefulInfo.pi_;
2494    double sumPi=0.0;
2495    for (int i=0;i<numberRows;i++) 
2496      sumPi += fabs(pi[i]);
2497    sumPi /= static_cast<double> (numberRows);
2498    // and scale back
2499    sumPi *= 0.01;
2500    usefulInfo.defaultDual_ = sumPi; // switch on
2501    int numberColumns = solver->getNumCols();
2502    int size = CoinMax(numberColumns,2*numberRows);
2503    usefulInfo.usefulRegion_ = new double [size];
2504    CoinZeroN(usefulInfo.usefulRegion_,size);
2505    usefulInfo.indexRegion_ = new int [size];
2506    // pi may change
2507    usefulInfo.pi_=CoinCopyOfArray(usefulInfo.pi_,numberRows);
2508  }
2509  assert (auxiliaryInfo);
2510  double cutoff =model->getCutoff();
2511  const double * lower = solver->getColLower();
2512  const double * upper = solver->getColUpper();
2513  // See if user thinks infeasible
2514  int anyAction=model->problemFeasibility()->feasible(model,0);
2515  if (anyAction) {
2516    // will return -2 if infeasible , 0 if treat as integer
2517    return anyAction-1;
2518  }
2519  int i;
2520  int saveStateOfSearch = model->stateOfSearch()%10;
2521  int numberStrong=model->numberStrong();
2522  /* Ranging is switched off.
2523     The idea is that you can find out the effect of one iteration
2524     on each unsatisfied variable cheaply.  Then use this
2525     if you have not got much else to go on.
2526  */
2527  //#define RANGING
2528#ifdef RANGING
2529  // must have clp
2530#ifndef COIN_HAS_CLP
2531#  warning("Ranging switched off as not Clp");
2532#undef RANGING
2533#endif
2534  // Pass number
2535  int kPass=0;
2536  int numberRows = solver->getNumRows();
2537#endif
2538  int numberColumns = model->getNumCols();
2539  double * saveUpper = new double[numberColumns];
2540  double * saveLower = new double[numberColumns];
2541  for (i=0;i<numberColumns;i++) {
2542    saveLower[i] = lower[i];
2543    saveUpper[i] = upper[i];
2544  }
2545 
2546  // Save solution in case heuristics need good solution later
2547 
2548  double * saveSolution = new double[numberColumns];
2549  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
2550  model->reserveCurrentSolution(saveSolution);
2551  const double * hotstartSolution = model->hotstartSolution();
2552  const int * hotstartPriorities = model->hotstartPriorities();
2553  double integerTolerance = 
2554    model->getDblParam(CbcModel::CbcIntegerTolerance);
2555  if (hotstartSolution) {
2556    if((model->moreSpecialOptions()&1024)!=0) {
2557      int nBad=0;
2558      int nUnsat=0;
2559      int nDiff=0;
2560      for (int i=0;i<numberObjects;i++) {
2561        OsiObject * object = model->modifiableObject(i);
2562        const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
2563        if (thisOne) {
2564          int iColumn = thisOne->columnNumber();
2565          double targetValue = hotstartSolution[iColumn];
2566          double value = saveSolution[iColumn];
2567          if (fabs(value-floor(value+0.5))>1.0e-6) {
2568            nUnsat++;
2569#ifdef CLP_INVESTIGATE
2570            printf("H %d is %g target %g\n",iColumn,value,targetValue);
2571#endif
2572          } else if (fabs(targetValue-value)>1.0e-6) {
2573            nDiff++;
2574          }
2575          if (targetValue<saveLower[iColumn]||
2576              targetValue>saveUpper[iColumn]) {
2577#ifdef CLP_INVESTIGATE
2578            printf("%d has target %g and current bounds %g and %g\n",
2579                   iColumn,targetValue,saveLower[iColumn],saveUpper[iColumn]);
2580#endif
2581            nBad++;
2582          }
2583        }
2584      }
2585#ifdef CLP_INVESTIGATE
2586      printf("Hot %d unsatisfied, %d outside limits, %d different\n",
2587             nUnsat,nBad,nDiff); 
2588#endif
2589      if (nBad) {
2590        // switch off as not possible
2591          hotstartSolution=NULL;
2592          model->setHotstartSolution(NULL,NULL);
2593          usefulInfo.hotstartSolution_=NULL;
2594      }
2595    }
2596  }
2597  /*
2598    Get a branching decision object. Use the default dynamic decision criteria unless
2599    the user has loaded a decision method into the model.
2600  */
2601  CbcBranchDecision *decision = model->branchingMethod();
2602  if (!decision)
2603    decision = new CbcBranchDynamicDecision();
2604  int xMark=0;
2605  // Get arrays to sort
2606  double * sort = new double[numberObjects];
2607  int * whichObject = new int[numberObjects];
2608#ifdef RANGING
2609  int xPen=0;
2610  int * objectMark = new int[2*numberObjects+1];
2611#endif
2612  // Arrays with movements
2613  double * upEstimate = new double[numberObjects];
2614  double * downEstimate = new double[numberObjects];
2615  double estimatedDegradation=0.0; 
2616  int numberNodes=model->getNodeCount();
2617  int saveLogLevel = model->logLevel();
2618#if 0
2619  if ((numberNodes%500)==0) {
2620    model->setLogLevel(6);
2621    // Get average up and down costs
2622    double averageUp=0.0;
2623    double averageDown=0.0;
2624    int numberUp=0;
2625    int numberDown=0;
2626    int i;
2627    for ( i=0;i<numberObjects;i++) {
2628      OsiObject * object = model->modifiableObject(i);
2629      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2630        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2631      assert(dynamicObject);
2632      int  numberUp2=0;
2633      int numberDown2=0;
2634      double up=0.0;
2635      double down=0.0;
2636      if (dynamicObject->numberTimesUp()) {
2637        numberUp++;
2638        averageUp += dynamicObject->upDynamicPseudoCost();
2639        numberUp2 += dynamicObject->numberTimesUp();
2640        up = dynamicObject->upDynamicPseudoCost();
2641      }
2642      if (dynamicObject->numberTimesDown()) {
2643        numberDown++;
2644        averageDown += dynamicObject->downDynamicPseudoCost();
2645        numberDown2 += dynamicObject->numberTimesDown();
2646        down = dynamicObject->downDynamicPseudoCost();
2647      }
2648      if (numberUp2||numberDown2)
2649        printf("col %d - up %d times cost %g, - down %d times cost %g\n",
2650               dynamicObject->columnNumber(),numberUp2,up,numberDown2,down);
2651    }
2652    if (numberUp)
2653      averageUp /= static_cast<double> (numberUp);
2654    else
2655      averageUp=1.0;
2656    if (numberDown)
2657      averageDown /= static_cast<double> (numberDown);
2658    else
2659      averageDown=1.0;
2660    printf("total - up %d vars average %g, - down %d vars average %g\n",
2661           numberUp,averageUp,numberDown,averageDown);
2662  }
2663#endif
2664  int numberBeforeTrust = model->numberBeforeTrust();
2665  // May go round twice if strong branching fixes all local candidates
2666  bool finished=false;
2667  int numberToFix=0;
2668# ifdef COIN_HAS_CLP
2669  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
2670  int saveClpOptions=0;
2671  if (osiclp) {
2672    // for faster hot start
2673    saveClpOptions = osiclp->specialOptions();
2674    osiclp->setSpecialOptions(saveClpOptions|8192);
2675  }
2676# else
2677  OsiSolverInterface *osiclp = NULL ;
2678# endif
2679  const CglTreeProbingInfo * probingInfo = NULL; //model->probingInfo();
2680  // Old code left in with DEPRECATED_STRATEGY
2681  assert (model->searchStrategy()==-1||
2682          model->searchStrategy()==1||
2683          model->searchStrategy()==2);
2684#ifdef DEPRECATED_STRATEGY
2685  int saveSearchStrategy2 = model->searchStrategy();
2686#endif
2687  // Get average up and down costs
2688  {
2689    double averageUp=0.0;
2690    double averageDown=0.0;
2691    int numberUp=0;
2692    int numberDown=0;
2693    int i;
2694    for ( i=0;i<numberObjects;i++) {
2695      OsiObject * object = model->modifiableObject(i);
2696      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2697        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2698      if (dynamicObject) {
2699        if (dynamicObject->numberTimesUp()) {
2700          numberUp++;
2701          averageUp += dynamicObject->upDynamicPseudoCost();
2702        }
2703        if (dynamicObject->numberTimesDown()) {
2704          numberDown++;
2705          averageDown += dynamicObject->downDynamicPseudoCost();
2706        }
2707      }
2708    }
2709    if (numberUp) 
2710      averageUp /= static_cast<double> (numberUp);
2711    else
2712      averageUp=1.0;
2713    if (numberDown) 
2714      averageDown /= static_cast<double> (numberDown);
2715    else
2716      averageDown=1.0;
2717    for ( i=0;i<numberObjects;i++) {
2718      OsiObject * object = model->modifiableObject(i);
2719      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2720        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2721      if (dynamicObject) {
2722        if (!dynamicObject->numberTimesUp()) 
2723          dynamicObject->setUpDynamicPseudoCost(averageUp);
2724        if (!dynamicObject->numberTimesDown()) 
2725          dynamicObject->setDownDynamicPseudoCost(averageDown);
2726      }
2727    }
2728  }
2729  /*
2730    1 strong
2731    2 no strong
2732    3 strong just before solution
2733    4 no strong just before solution
2734    5 strong first time or before solution
2735    6 strong first time
2736  */
2737  int useShadow = model->moreSpecialOptions()&7;
2738  if (useShadow>2) {
2739    if (model->getSolutionCount()) {
2740      if (numberNodes||useShadow<5) {
2741        useShadow=0;
2742        // zap pseudo shadow prices
2743        model->pseudoShadow(-1);
2744        // and switch off
2745        model->setMoreSpecialOptions(model->moreSpecialOptions()&(~1023));
2746      } else {
2747        useShadow=1;
2748      }
2749    } else if (useShadow<5) {
2750      useShadow -= 2;
2751    } else {
2752      useShadow=1;
2753    }
2754  } 
2755  if (useShadow) {
2756    // pseudo shadow prices
2757    model->pseudoShadow((model->moreSpecialOptions()>>3)&63);
2758  }
2759#ifdef DEPRECATED_STRATEGY
2760  { // in for tabbing
2761  } else if (saveSearchStrategy2<1999) {
2762    // pseudo shadow prices
2763    model->pseudoShadow(NULL,NULL);
2764  } else if (saveSearchStrategy2<2999) {
2765    // leave old ones
2766  } else if (saveSearchStrategy2<3999) {
2767    // pseudo shadow prices at root
2768    if (!numberNodes)
2769      model->pseudoShadow(NULL,NULL);
2770  } else {
2771    abort();
2772  }
2773  if (saveSearchStrategy2>=0)
2774    saveSearchStrategy2 = saveSearchStrategy2 % 1000;
2775  if (saveSearchStrategy2==999)
2776    saveSearchStrategy2=-1;
2777  int saveSearchStrategy = saveSearchStrategy2<99 ? saveSearchStrategy2 : saveSearchStrategy2-100;
2778#endif //DEPRECATED_STRATEGY
2779  int numberNotTrusted=0;
2780  int numberStrongDone=0;
2781  int numberUnfinished=0;
2782  int numberStrongInfeasible=0;
2783  int numberStrongIterations=0;
2784  // so we can save lots of stuff
2785  CbcStrongInfo choice;
2786  CbcDynamicPseudoCostBranchingObject * choiceObject = NULL;
2787  if (model->allDynamic()) {
2788    CbcSimpleIntegerDynamicPseudoCost * object = NULL;
2789    choiceObject=new CbcDynamicPseudoCostBranchingObject(model,0,-1,0.5,object);
2790  }
2791  choice.possibleBranch=choiceObject;
2792  numberPassesLeft = CoinMax(numberPassesLeft,2);
2793  while(!finished) {
2794    numberPassesLeft--;
2795    finished=true;
2796    decision->initialize(model);
2797    // Some objects may compute an estimate of best solution from here
2798    estimatedDegradation=0.0; 
2799    numberToFix=0;
2800    int numberToDo=0;
2801    int iBestNot=-1;
2802    int iBestGot=-1;
2803    double best=0.0;
2804    numberNotTrusted=0;
2805    numberStrongDone=0;
2806    numberUnfinished=0;
2807    numberStrongInfeasible=0;
2808    numberStrongIterations=0;
2809#ifdef RANGING
2810    int * which = objectMark+numberObjects+1;
2811    int neededPenalties;
2812    int optionalPenalties;
2813#endif
2814    // We may go round this loop three times (only if we think we have solution)
2815    for (int iPass=0;iPass<3;iPass++) {
2816     
2817      // Some objects may compute an estimate of best solution from here
2818      estimatedDegradation=0.0; 
2819      numberUnsatisfied_ = 0;
2820      // initialize sum of "infeasibilities"
2821      sumInfeasibilities_ = 0.0;
2822      int bestPriority=COIN_INT_MAX;
2823      int number01 = 0;
2824      const cliqueEntry * entry = NULL;
2825      const int * toZero = NULL;
2826      const int * toOne = NULL;
2827      const int * backward = NULL;
2828      int numberUnsatisProbed=0;
2829      int numberUnsatisNotProbed=0; // 0-1
2830      if (probingInfo) {
2831        number01 = probingInfo->numberIntegers();
2832        entry = probingInfo->fixEntries();
2833        toZero = probingInfo->toZero();
2834        toOne = probingInfo->toOne();
2835        backward = probingInfo->backward();
2836        if (!toZero[number01]||number01<numberObjects||true) {
2837          // no info
2838          probingInfo=NULL;
2839        }
2840      }
2841      /*
2842        Scan for branching objects that indicate infeasibility. Choose candidates
2843        using priority as the first criteria, then integer infeasibility.
2844       
2845        The algorithm is to fill the array with a set of good candidates (by
2846        infeasibility) with priority bestPriority.  Finding a candidate with
2847        priority better (less) than bestPriority flushes the choice array. (This
2848        serves as initialization when the first candidate is found.)
2849       
2850      */
2851      numberToDo=0;
2852#ifdef RANGING
2853      neededPenalties=0;
2854      optionalPenalties=numberObjects;
2855#endif
2856      iBestNot=-1;
2857      double bestNot=0.0;
2858      iBestGot=-1;
2859      best=0.0;
2860      /* Problem type as set by user or found by analysis.  This will be extended
2861         0 - not known
2862         1 - Set partitioning <=
2863         2 - Set partitioning ==
2864         3 - Set covering
2865         4 - all +- 1 or all +1 and odd
2866      */
2867      int problemType = model->problemType();
2868      bool canDoOneHot=false;
2869      for (i=0;i<numberObjects;i++) {
2870        OsiObject * object = model->modifiableObject(i);
2871        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2872          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2873        int preferredWay;
2874        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
2875        int priorityLevel = object->priority();
2876        if (hotstartSolution) {
2877          // we are doing hot start
2878          const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
2879          if (thisOne) {
2880            int iColumn = thisOne->columnNumber();
2881            bool canDoThisHot=true;
2882            double targetValue = hotstartSolution[iColumn];
2883            if (saveUpper[iColumn]>saveLower[iColumn]) {
2884              double value = saveSolution[iColumn];
2885              if (hotstartPriorities)
2886                priorityLevel=hotstartPriorities[iColumn]; 
2887              //double originalLower = thisOne->originalLower();
2888              //double originalUpper = thisOne->originalUpper();
2889              // switch off if not possible
2890              if (targetValue>=saveLower[iColumn]&&targetValue<=saveUpper[iColumn]) {
2891                /* priority outranks rest always if negative
2892                   otherwise can be downgraded if at correct level.
2893                   Infeasibility may be increased to choose 1.0 values first.
2894                   choose one near wanted value
2895                */
2896                if (fabs(value-targetValue)>integerTolerance) {
2897                  //if (infeasibility>0.01)
2898                  //infeasibility = fabs(1.0e6-fabs(value-targetValue));
2899                  //else
2900                  infeasibility = fabs(value-targetValue);
2901                  //if (targetValue==1.0)
2902                  //infeasibility += 1.0;
2903                  if (value>targetValue) {
2904                    preferredWay=-1;
2905                  } else {
2906                    preferredWay=1;
2907                  }
2908                  priorityLevel = CoinAbs(priorityLevel);
2909                } else if (priorityLevel<0) {
2910                  priorityLevel = CoinAbs(priorityLevel);
2911                  if (targetValue==saveLower[iColumn]) {
2912                    infeasibility = integerTolerance+1.0e-12;
2913                    preferredWay=-1;
2914                  } else if (targetValue==saveUpper[iColumn]) {
2915                    infeasibility = integerTolerance+1.0e-12;
2916                    preferredWay=1;
2917                  } else {
2918                    // can't
2919                    priorityLevel += 10000000;
2920                    canDoThisHot=false;
2921                  }
2922                } else {
2923                  priorityLevel += 10000000;
2924                  canDoThisHot=false;
2925                }
2926              } else {
2927                // switch off if not possible
2928                canDoThisHot=false;
2929              }
2930              if (canDoThisHot)
2931                canDoOneHot=true;
2932            } else if (targetValue<saveLower[iColumn]||targetValue>saveUpper[iColumn]) {
2933            }
2934          } else {
2935            priorityLevel += 10000000;
2936          }
2937        }
2938#define ZERO_ONE 0
2939#define ZERO_FAKE 1.0e20;
2940#if ZERO_ONE==1
2941        // branch on 0-1 first (temp)
2942        if (fabs(saveSolution[dynamicObject->columnNumber()])<1.0)
2943          priorityLevel--;
2944#endif
2945#if ZERO_ONE==2
2946        if (fabs(saveSolution[dynamicObject->columnNumber()])<1.0)
2947          infeasibility *= ZERO_FAKE;
2948#endif
2949        if (infeasibility) {
2950          int iColumn = numberColumns+i;
2951          bool gotDown=false;
2952          int numberThisDown = 0;
2953          bool gotUp=false;
2954          int numberThisUp = 0;
2955          double downGuess=object->downEstimate();
2956          double upGuess=object->upEstimate();
2957          if (dynamicObject) {
2958            // Use this object's numberBeforeTrust
2959            int numberBeforeTrust = dynamicObject->numberBeforeTrust();
2960            iColumn = dynamicObject->columnNumber();
2961            gotDown=false;
2962            numberThisDown = dynamicObject->numberTimesDown();
2963            if (numberThisDown>=numberBeforeTrust)
2964              gotDown=true;
2965            gotUp=false;
2966            numberThisUp = dynamicObject->numberTimesUp();
2967            if (numberThisUp>=numberBeforeTrust)
2968              gotUp=true;
2969            if (!depth_&&false) {
2970              // try closest to 0.5
2971              double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2972              infeasibility = fabs(0.5-part);
2973            }
2974            if (problemType>0&&problemType<4&&false) {
2975              // try closest to 0.5
2976              double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2977              infeasibility = 0.5-fabs(0.5-part);
2978            }
2979            if (probingInfo) {
2980              int iSeq = backward[iColumn];
2981              assert (iSeq>=0);
2982              infeasibility = 1.0 + (toZero[iSeq+1]-toZero[iSeq])+ 
2983                5.0*CoinMin(toOne[iSeq]-toZero[iSeq],toZero[iSeq+1]-toOne[iSeq]);
2984              if (toZero[iSeq+1]>toZero[iSeq]) {
2985                numberUnsatisProbed++;
2986              } else {
2987                numberUnsatisNotProbed++;
2988              }
2989            }
2990          } else {
2991            // see if SOS
2992            CbcSOS * sosObject =
2993              dynamic_cast <CbcSOS *>(object) ;
2994            if (sosObject) {
2995              gotDown=false;
2996              numberThisDown = sosObject->numberTimesDown();
2997              if (numberThisDown>=numberBeforeTrust)
2998                gotDown=true;
2999              gotUp=false;
3000              numberThisUp = sosObject->numberTimesUp();
3001              if (numberThisUp>=numberBeforeTrust)
3002                gotUp=true;
3003            } else {
3004              gotDown=true;
3005              numberThisDown=999999;
3006              downGuess=1.0e20;
3007              gotUp=true;
3008              numberThisUp=999999;
3009              upGuess=1.0e20;
3010              numberPassesLeft=0;
3011            }
3012          }
3013          // Increase estimated degradation to solution
3014          estimatedDegradation += CoinMin(downGuess,upGuess);
3015          downEstimate[i]=downGuess;
3016          upEstimate[i]=upGuess;
3017          numberUnsatisfied_++;
3018          sumInfeasibilities_ += infeasibility;
3019          // Better priority? Flush choices.
3020          if (priorityLevel<bestPriority) {
3021            numberToDo=0;
3022            bestPriority = priorityLevel;
3023            iBestGot=-1;
3024            best=0.0;
3025            numberNotTrusted=0;
3026#ifdef RANGING
3027            neededPenalties=0;
3028            optionalPenalties=numberObjects;
3029#endif
3030          } else if (priorityLevel>bestPriority) {
3031            continue;
3032          }
3033          if (!gotUp||!gotDown)
3034            numberNotTrusted++;
3035          // Check for suitability based on infeasibility.
3036          if ((gotDown&&gotUp)&&numberStrong>0) {
3037            sort[numberToDo]=-infeasibility;
3038            if (infeasibility>best) {
3039              best=infeasibility;
3040              iBestGot=numberToDo;
3041            }
3042#ifdef RANGING
3043            if (dynamicObject) {
3044              objectMark[--optionalPenalties]=numberToDo;
3045              which[optionalPenalties]=iColumn;
3046            }
3047#endif
3048          } else {
3049#ifdef RANGING
3050            if (dynamicObject) {
3051              objectMark[neededPenalties]=numberToDo;
3052              which[neededPenalties++]=iColumn;
3053            }
3054#endif
3055            sort[numberToDo]=-10.0*infeasibility;
3056            if (!(numberThisUp+numberThisDown))
3057              sort[numberToDo] *= 100.0; // make even more likely
3058            if (iColumn<numberColumns) {
3059              double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
3060              if (1.0-fabs(part-0.5)>bestNot) {
3061                iBestNot=numberToDo;
3062                bestNot = 1.0-fabs(part-0.5);
3063              }
3064            } else {
3065              // SOS
3066              if (-sort[numberToDo]>bestNot) {
3067                iBestNot=numberToDo;
3068                bestNot = -sort[numberToDo];
3069              }
3070            }
3071          }
3072          if (model->messageHandler()->logLevel()>3) { 
3073            printf("%d (%d) down %d %g up %d %g - infeas %g - sort %g solution %g\n",
3074                   i,iColumn,numberThisDown,object->downEstimate(),numberThisUp,object->upEstimate(),
3075                   infeasibility,sort[numberToDo],saveSolution[iColumn]);
3076          }
3077          whichObject[numberToDo++]=i;
3078        } else {
3079          // for debug
3080          downEstimate[i]=-1.0;
3081          upEstimate[i]=-1.0;
3082        }
3083      }
3084      if (numberUnsatisfied_) {
3085        if (probingInfo&&false)
3086          printf("nunsat %d, %d probed, %d other 0-1\n",numberUnsatisfied_,
3087                 numberUnsatisProbed,numberUnsatisNotProbed);
3088        // some infeasibilities - go to next steps
3089        if (!canDoOneHot&&hotstartSolution) {
3090          // switch off as not possible
3091          hotstartSolution=NULL;
3092          model->setHotstartSolution(NULL,NULL);
3093          usefulInfo.hotstartSolution_=NULL;
3094        }
3095        break;
3096      } else if (!iPass) {
3097        // may just need resolve
3098        model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3099        //double newObjValue = solver->getObjSense()*solver->getObjValue();
3100        //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3101        if (!solver->isProvenOptimal()) {
3102          // infeasible
3103          anyAction=-2;
3104          break;
3105        }
3106      } else if (iPass==1) {
3107        // looks like a solution - get paranoid
3108        bool roundAgain=false;
3109        // get basis
3110        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
3111        if (!ws)
3112          break;
3113        double tolerance;
3114        solver->getDblParam(OsiPrimalTolerance,tolerance);
3115        for (i=0;i<numberColumns;i++) {
3116          double value = saveSolution[i];
3117          if (value<lower[i]-tolerance) {
3118            saveSolution[i]=lower[i];
3119            roundAgain=true;
3120            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
3121          } else if (value>upper[i]+tolerance) {
3122            saveSolution[i]=upper[i];
3123            roundAgain=true;
3124            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
3125          } 
3126        }
3127        if (roundAgain) {
3128          // restore basis
3129          solver->setWarmStart(ws);
3130          solver->setColSolution(saveSolution);
3131          delete ws;
3132          bool takeHint;
3133          OsiHintStrength strength;
3134          solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
3135          solver->setHintParam(OsiDoDualInResolve,false,OsiHintDo) ;
3136          model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3137          //double newObjValue = solver->getObjSense()*solver->getObjValue();
3138          //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3139          solver->setHintParam(OsiDoDualInResolve,takeHint,strength) ;
3140          if (!solver->isProvenOptimal()) {
3141            // infeasible
3142            anyAction=-2;
3143            break;
3144          }
3145        } else {
3146          delete ws;
3147          break;
3148        }
3149      }
3150    }
3151    if (anyAction==-2) {
3152      break;
3153    }
3154    // skip if solution
3155    if (!numberUnsatisfied_)
3156      break;
3157    int skipAll = (numberNotTrusted==0||numberToDo==1) ? 1 : 0;
3158    bool doneHotStart=false;
3159    //DEPRECATED_STRATEGYint searchStrategy = saveSearchStrategy>=0 ? (saveSearchStrategy%10) : -1;
3160    int searchStrategy = model->searchStrategy();
3161    // But adjust depending on ratio of iterations
3162    if (searchStrategy>0) {
3163      if (numberBeforeTrust>=5&&numberBeforeTrust<=10) {
3164        if (searchStrategy!=2) {
3165          assert (searchStrategy==1);
3166          if (depth_>5) {
3167            int numberIterations = model->getIterationCount();
3168            int numberStrongIterations = model->numberStrongIterations();
3169            if (numberStrongIterations>numberIterations+10000) {
3170              searchStrategy=2;
3171              skipAll=1;
3172            } else if (numberStrongIterations*4+1000<numberIterations) {
3173              searchStrategy=3;
3174              skipAll=0;
3175            }
3176          } else {
3177            searchStrategy=3;
3178            skipAll=0;
3179          }
3180        }
3181      }
3182    }
3183    // worth trying if too many times
3184    // Save basis
3185    CoinWarmStart * ws = NULL;
3186    // save limit
3187    int saveLimit=0;
3188    solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
3189    if (!numberPassesLeft)
3190      skipAll=1;
3191    if (!skipAll) {
3192      ws = solver->getWarmStart();
3193      int limit=100;
3194      if (!saveStateOfSearch&&saveLimit<limit&&saveLimit==100)
3195        solver->setIntParam(OsiMaxNumIterationHotStart,limit); 
3196    }
3197    // Say which one will be best
3198    int whichChoice=0;
3199    int bestChoice;
3200    if (iBestGot>=0)
3201      bestChoice=iBestGot;
3202    else
3203      bestChoice=iBestNot;
3204    assert (bestChoice>=0);
3205    // If we have hit max time don't do strong branching
3206    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
3207                        model->getDblParam(CbcModel::CbcMaximumSeconds));
3208    // also give up if we are looping round too much
3209    if (hitMaxTime||numberPassesLeft<=0||useShadow==2) {
3210      int iObject = whichObject[bestChoice];
3211      OsiObject * object = model->modifiableObject(iObject);
3212      int preferredWay;
3213      object->infeasibility(&usefulInfo,preferredWay);
3214      CbcObject * obj =
3215        dynamic_cast <CbcObject *>(object) ;
3216      assert (obj);
3217      branch_=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3218      {
3219        CbcBranchingObject * branchObj =
3220          dynamic_cast <CbcBranchingObject *>(branch_) ;
3221        assert (branchObj);
3222        branchObj->way(preferredWay);
3223      }
3224      delete ws;
3225      ws=NULL;
3226      break;
3227    } else {
3228      // say do fast
3229      int easy=1;
3230      if (!skipAll)
3231        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
3232      int iDo;
3233#ifdef RANGING
3234      bool useRanging = model->allDynamic()&&!skipAll;
3235      if (useRanging) {
3236        double currentObjective = solver->getObjValue()*solver->getObjSense();
3237        double gap = cutoff-currentObjective;
3238        // relax a bit
3239        gap *= 1.0000001;
3240        gap = CoinMax(1.0e-5,gap);
3241        // off penalties if too much
3242        double needed = neededPenalties;
3243        needed *= numberRows;
3244        if (numberNodes) {
3245          if (needed>1.0e6) {
3246            neededPenalties=0;
3247          } else if (gap<1.0e5) {
3248            // maybe allow some not needed
3249            int extra = static_cast<int> ((1.0e6-needed)/numberRows);
3250            int nStored = numberObjects-optionalPenalties;
3251            extra = CoinMin(extra,nStored);
3252            for (int i=0;i<extra;i++) {
3253              objectMark[neededPenalties]=objectMark[optionalPenalties+i];
3254                which[neededPenalties++]=which[optionalPenalties+i];;
3255            }
3256          }
3257        }
3258        if (osiclp&&neededPenalties) {
3259          assert (!doneHotStart);
3260          xPen += neededPenalties;
3261          which--;
3262          which[0]=neededPenalties;
3263          osiclp->passInRanges(which);
3264          // Mark hot start and get ranges
3265          if (kPass) {
3266            // until can work out why solution can go funny
3267            int save = osiclp->specialOptions();
3268            osiclp->setSpecialOptions(save|256);
3269            solver->markHotStart();
3270            osiclp->setSpecialOptions(save);
3271          } else {
3272            solver->markHotStart();
3273          }
3274          doneHotStart=true;
3275          xMark++;
3276          kPass++;
3277          osiclp->passInRanges(NULL);
3278          const double * downCost=osiclp->upRange();
3279          const double * upCost=osiclp->downRange();
3280          bool problemFeasible=true;
3281          int numberFixed=0;
3282          for (int i=0;i<neededPenalties;i++) {
3283            int j = objectMark[i];
3284            int iObject = whichObject[j];
3285            OsiObject * object = model->modifiableObject(iObject);
3286            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3287              dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3288            // Use this object's numberBeforeTrust
3289            int numberBeforeTrust = dynamicObject->numberBeforeTrust();
3290            int iSequence=dynamicObject->columnNumber();
3291            double value = saveSolution[iSequence];
3292            value -= floor(value);
3293            double upPenalty = CoinMin(upCost[i],1.0e110)*(1.0-value);
3294            double downPenalty = CoinMin(downCost[i],1.0e110)*value;
3295            int numberThisDown = dynamicObject->numberTimesDown();
3296            int numberThisUp = dynamicObject->numberTimesUp();
3297            if (!numberBeforeTrust) {
3298              // override
3299              downEstimate[iObject]=downPenalty;
3300              upEstimate[iObject]=upPenalty;
3301              double min1 = CoinMin(downEstimate[iObject],
3302                                    upEstimate[iObject]);
3303              double max1 = CoinMax(downEstimate[iObject],
3304                                    upEstimate[iObject]);
3305              min1 = 0.8*min1+0.2*max1;
3306              sort[j] = - min1;
3307            } else if (numberThisDown<numberBeforeTrust||
3308                       numberThisUp<numberBeforeTrust) {
3309              double invTrust = 1.0/static_cast<double> (numberBeforeTrust);
3310              if (numberThisDown<numberBeforeTrust) {
3311                double fraction = numberThisDown*invTrust;
3312                downEstimate[iObject] = fraction*downEstimate[iObject]+(1.0-fraction)*downPenalty;
3313              }
3314              if (numberThisUp<numberBeforeTrust) {
3315                double fraction = numberThisUp*invTrust;
3316                upEstimate[iObject] = fraction*upEstimate[iObject]+(1.0-fraction)*upPenalty;
3317              }
3318              double min1 = CoinMin(downEstimate[iObject],
3319                                    upEstimate[iObject]);
3320              double max1 = CoinMax(downEstimate[iObject],
3321                                    upEstimate[iObject]);
3322              min1 = 0.8*min1+0.2*max1;
3323              min1 *= 10.0;
3324              if (!(numberThisDown+numberThisUp))
3325                min1 *= 100.0;
3326              sort[j] = - min1;
3327            }
3328            if (CoinMax(downPenalty,upPenalty)>gap) {
3329              printf("gap %g object %d has down range %g, up %g\n",
3330                     gap,i,downPenalty,upPenalty);
3331              //sort[j] -= 1.0e50; // make more likely to be chosen
3332              int number;
3333              if (downPenalty>gap) {
3334                number = dynamicObject->numberTimesDown();
3335                if (upPenalty>gap) 
3336                  problemFeasible=false;
3337                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver,&usefulInfo,1);
3338                //branch->fix(solver,saveLower,saveUpper,1);
3339                delete branch;
3340              } else {
3341                number = dynamicObject->numberTimesUp();
3342                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver,&usefulInfo,1);
3343                //branch->fix(solver,saveLower,saveUpper,-1);
3344                delete branch;
3345              }
3346              if(number>=numberBeforeTrust)
3347                dynamicObject->setNumberBeforeTrust(number+1);
3348              numberFixed++;
3349            }
3350            if (!numberNodes)
3351              printf("%d pen down ps %g -> %g up ps %g -> %g\n",
3352                     iObject,downPenalty,downPenalty,upPenalty,upPenalty);
3353          }
3354          if(numberFixed&&problemFeasible) {
3355            assert(doneHotStart);
3356            solver->unmarkHotStart();
3357            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3358            //double newObjValue = solver->getObjSense()*solver->getObjValue();
3359            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3360            solver->markHotStart();
3361            problemFeasible = solver->isProvenOptimal();
3362          }
3363          if(!problemFeasible) {
3364            fprintf(stdout,"both ways infeas on ranging - code needed\n");
3365            anyAction=-2;
3366            if (!choiceObject) {
3367              delete choice.possibleBranch;
3368              choice.possibleBranch=NULL;
3369            }
3370            //printf("Both infeasible for choice %d sequence %d\n",i,
3371            // model->object(choice.objectNumber)->columnNumber());
3372            // Delete the snapshot
3373            solver->unmarkHotStart();
3374            // back to normal
3375            solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3376            // restore basis
3377            solver->setWarmStart(ws);
3378            doneHotStart=false;
3379            delete ws;
3380            ws=NULL;
3381            break;
3382          }
3383        }
3384      }
3385#endif          /* RANGING */
3386      {
3387        int numberIterations = model->getIterationCount();
3388        //numberIterations += (model->numberExtraIterations()>>2);
3389        const int * strongInfo = model->strongInfo();
3390        //int numberDone = strongInfo[0]-strongInfo[3];
3391        int numberFixed = strongInfo[1]-strongInfo[4];
3392        int numberInfeasible = strongInfo[2]-strongInfo[5];
3393        assert (!strongInfo[3]);
3394        assert (!strongInfo[4]);
3395        assert (!strongInfo[5]);
3396        int numberStrongIterations = model->numberStrongIterations();
3397        int numberRows = solver->getNumRows();
3398        if (numberStrongIterations>numberIterations+CoinMin(100,10*numberRows)&&depth_>=4&&numberNodes>100) {
3399          if (20*numberInfeasible+4*numberFixed<numberNodes) {
3400            // Say never do
3401            skipAll=-1;
3402          }
3403        } 
3404      }
3405      // make sure best will be first
3406      if (iBestGot>=0)
3407        sort[iBestGot]=-COIN_DBL_MAX;
3408      // Actions 0 - exit for repeat, 1 resolve and try old choice,2 exit for continue
3409      if (anyAction)
3410        numberToDo=0; // skip as we will be trying again
3411      // Sort
3412      CoinSort_2(sort,sort+numberToDo,whichObject);
3413      // Change in objective opposite infeasible
3414      double worstFeasible=0.0;
3415      // Just first if strong off
3416      if (!numberStrong)
3417        numberToDo=CoinMin(numberToDo,1);
3418      if (searchStrategy==2)
3419        numberToDo=CoinMin(numberToDo,20);
3420      iDo=0;
3421      int saveLimit2;
3422      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit2);
3423      bool doQuickly = false; // numberToDo>2*numberStrong;
3424      if (searchStrategy==2)
3425        doQuickly=true;
3426      int numberTest=numberNotTrusted>0 ? numberStrong : (numberStrong+1)/2;
3427      if (searchStrategy==3) {
3428        // Previously decided we need strong
3429        doQuickly=false;
3430        numberTest = numberStrong;
3431      }
3432      // Try nearly always off
3433      if (skipAll>=0) {
3434        if (searchStrategy<2) {
3435          //if ((numberNodes%20)!=0) {
3436            if ((model->specialOptions()&8)==0) {
3437              numberTest=0;
3438              doQuickly=true;
3439            }
3440            //} else {
3441            //doQuickly=false;
3442            //numberTest=2*numberStrong;
3443            //skipAll=0;
3444            //}
3445        }
3446      } else {
3447        // Just take first
3448        doQuickly=true;
3449        numberTest=1;
3450      }
3451      int testDepth = (skipAll>=0) ? 8 : 4;
3452      if (depth_<testDepth&&numberStrong) {
3453        if (searchStrategy!=2) {
3454          doQuickly=false;
3455          int numberRows = solver->getNumRows();
3456          // whether to do this or not is important - think
3457          if (numberRows<300||numberRows+numberColumns<2500) {
3458            if (depth_<7)
3459              numberStrong = CoinMin(3*numberStrong,numberToDo);
3460            if (!depth_) 
3461              numberStrong=CoinMin(6*numberStrong,numberToDo);
3462          }
3463          numberTest=numberStrong;
3464          skipAll=0;
3465        }
3466      }
3467      // Do at least 5 strong
3468      if (numberColumns<1000&&(depth_<15||numberNodes<1000000))
3469        numberTest = CoinMax(numberTest,5);
3470      if ((model->specialOptions()&8)==0) {
3471        if (skipAll) {
3472          numberTest=0;
3473          doQuickly=true;
3474        }
3475      } else {
3476        // do 5 as strong is fixing
3477        numberTest = CoinMax(numberTest,5);
3478      }
3479      // see if switched off
3480      if (skipAll<0) {
3481        numberTest=0;
3482        doQuickly=true;
3483      }
3484      int realMaxHotIterations=999999;
3485      if (skipAll<0)
3486        numberToDo=1;
3487#ifdef DO_ALL_AT_ROOT
3488      if (!numberNodes) {
3489        printf("DOX %d test %d unsat %d - nobj %d\n",
3490               numberToDo,numberTest,numberUnsatisfied_,
3491               numberObjects);
3492        numberTest=numberToDo;
3493      }
3494#endif
3495      for ( iDo=0;iDo<numberToDo;iDo++) {
3496        int iObject = whichObject[iDo];
3497        OsiObject * object = model->modifiableObject(iObject);
3498        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3499          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3500        int iColumn = dynamicObject ? dynamicObject->columnNumber() : numberColumns+iObject;
3501        int preferredWay;
3502        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
3503        // may have become feasible
3504        if (!infeasibility)
3505          continue;
3506#ifndef NDEBUG
3507        if (iColumn<numberColumns) {
3508            const double * solution = model->testSolution();
3509            assert (saveSolution[iColumn]==solution[iColumn]);
3510          }
3511#endif
3512        CbcSimpleInteger * obj =
3513          dynamic_cast <CbcSimpleInteger *>(object) ;
3514        if (obj) {
3515          if (choiceObject) {
3516            obj->fillCreateBranch(choiceObject,&usefulInfo,preferredWay);
3517            choiceObject->setObject(dynamicObject);
3518          } else {
3519            choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3520          }
3521        } else {
3522          CbcObject * obj =
3523            dynamic_cast <CbcObject *>(object) ;
3524          assert (obj);
3525          choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3526        }
3527        // Save which object it was
3528        choice.objectNumber=iObject;
3529        choice.numIntInfeasUp = numberUnsatisfied_;
3530        choice.numIntInfeasDown = numberUnsatisfied_;
3531        choice.upMovement = upEstimate[iObject];
3532        choice.downMovement = downEstimate[iObject];
3533        assert (choice.upMovement>=0.0);
3534        assert (choice.downMovement>=0.0);
3535        choice.fix=0; // say not fixed
3536        // see if can skip strong branching
3537        int canSkip = choice.possibleBranch->fillStrongInfo(choice);
3538        if ((numberTest<=0||skipAll)) {
3539          if (iDo>20) {
3540#ifdef DO_ALL_AT_ROOT
3541            if (!numberNodes)
3542              printf("DOY test %d - iDo %d\n",
3543                     numberTest,iDo);
3544#endif
3545            if (!choiceObject) {
3546              delete choice.possibleBranch;
3547              choice.possibleBranch=NULL;
3548            }
3549            break; // give up anyway
3550          }
3551        }
3552        if (model->messageHandler()->logLevel()>3&&numberBeforeTrust&&dynamicObject) 
3553          dynamicObject->print(1,choice.possibleBranch->value());
3554        if (skipAll<0)
3555          canSkip=true;
3556        if (!canSkip) {
3557          if (!doneHotStart) {
3558            // Mark hot start
3559            doneHotStart=true;
3560            solver->markHotStart();
3561            xMark++;
3562          }
3563          numberTest--;
3564          // just do a few
3565          if (searchStrategy==2)
3566            solver->setIntParam(OsiMaxNumIterationHotStart,10); 
3567          double objectiveChange ;
3568          double newObjectiveValue=1.0e100;
3569          int j;
3570          // status is 0 finished, 1 infeasible and other
3571          int iStatus;
3572          /*
3573            Try the down direction first. (Specify the initial branching alternative as
3574            down with a call to way(-1). Each subsequent call to branch() performs the
3575            specified branch and advances the branch object state to the next branch
3576            alternative.)
3577          */
3578          choice.possibleBranch->way(-1) ;
3579          choice.possibleBranch->branch() ;
3580          solver->solveFromHotStart() ;
3581          bool needHotStartUpdate=false;
3582          numberStrongDone++;
3583          numberStrongIterations += solver->getIterationCount();
3584          /*
3585            We now have an estimate of objective degradation that we can use for strong
3586            branching. If we're over the cutoff, the variable is monotone up.
3587            If we actually made it to optimality, check for a solution, and if we have
3588            a good one, call setBestSolution to process it. Note that this may reduce the
3589            cutoff, so we check again to see if we can declare this variable monotone.
3590          */
3591          if (solver->isProvenOptimal())
3592            iStatus=0; // optimal
3593          else if (solver->isIterationLimitReached()
3594                   &&!solver->isDualObjectiveLimitReached())
3595            iStatus=2; // unknown
3596          else
3597            iStatus=1; // infeasible
3598          if (iStatus!=2&&solver->getIterationCount()>
3599              realMaxHotIterations)
3600            numberUnfinished++;
3601          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3602          choice.numItersDown = solver->getIterationCount();
3603          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3604          // Update branching information if wanted
3605          CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3606          if (cbcobj) {
3607            CbcObject * object = cbcobj->object();
3608            assert (object) ;
3609            CbcObjectUpdateData update = object->createUpdateInformation(solver,this,cbcobj);
3610            update.objectNumber_ = choice.objectNumber;
3611            model->addUpdateInformation(update);
3612          } else {
3613            decision->updateInformation( solver,this);
3614          }
3615          if (!iStatus) {
3616            choice.finishedDown = true ;
3617            if (newObjectiveValue>=cutoff) {
3618              objectiveChange = 1.0e100; // say infeasible
3619              numberStrongInfeasible++;
3620            } else {
3621              //#define TIGHTEN_BOUNDS
3622#ifdef TIGHTEN_BOUNDS
3623              // Can we tighten bounds?
3624              if (iColumn<numberColumns&&cutoff<1.0e20
3625                  &&objectiveChange>1.0e-5) {
3626                double value = saveSolution[iColumn];
3627                double down = value-floor(value);
3628                double changePer = objectiveChange/(down+1.0e-7);
3629                double distance = (cutoff-objectiveValue_)/changePer;
3630                distance += 1.0e-3;
3631                if (distance<5.0) {
3632                  double newLower = ceil(value-distance);
3633                  if (newLower>saveLower[iColumn]) {
3634                    //printf("Could increase lower bound on %d from %g to %g\n",
3635                    //   iColumn,saveLower[iColumn],newLower);
3636                    saveLower[iColumn]=newLower;
3637                    solver->setColLower(iColumn,newLower);
3638                  }
3639                }
3640                }
3641#endif
3642              // See if integer solution
3643              if (model->feasibleSolution(choice.numIntInfeasDown,
3644                                          choice.numObjInfeasDown)
3645                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
3646                if (auxiliaryInfo->solutionAddsCuts()) {
3647                  needHotStartUpdate=true;
3648                  solver->unmarkHotStart();
3649                }
3650                model->setBestSolution(CBC_STRONGSOL,
3651                                       newObjectiveValue,
3652                                       solver->getColSolution()) ;
3653                if (needHotStartUpdate) {
3654                  model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3655                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3656                  //objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
3657                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3658                  model->feasibleSolution(choice.numIntInfeasDown,
3659                                          choice.numObjInfeasDown);
3660                }
3661                model->setLastHeuristic(NULL);
3662                model->incrementUsed(solver->getColSolution());
3663                cutoff =model->getCutoff();
3664                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3665                  objectiveChange = 1.0e100 ;
3666                  numberStrongInfeasible++;
3667                }
3668              }
3669            }
3670          } else if (iStatus==1) {
3671            objectiveChange = 1.0e100 ;
3672            numberStrongInfeasible++;
3673          } else {
3674            // Can't say much as we did not finish
3675            choice.finishedDown = false ;
3676            numberUnfinished++;
3677          }
3678          choice.downMovement = objectiveChange ;
3679         
3680          // restore bounds
3681          for ( j=0;j<numberColumns;j++) {
3682            if (saveLower[j] != lower[j])
3683              solver->setColLower(j,saveLower[j]);
3684            if (saveUpper[j] != upper[j])
3685              solver->setColUpper(j,saveUpper[j]);
3686          }
3687          if(needHotStartUpdate) {
3688            needHotStartUpdate = false;
3689            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3690            //double newObjValue = solver->getObjSense()*solver->getObjValue();
3691            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3692            //we may again have an integer feasible solution
3693            int numberIntegerInfeasibilities;
3694            int numberObjectInfeasibilities;
3695            if (model->feasibleSolution(
3696                                        numberIntegerInfeasibilities,
3697                                        numberObjectInfeasibilities)) {
3698#ifdef BONMIN
3699              //In this case node has become integer feasible, let us exit the loop
3700              std::cout<<"Node has become integer feasible"<<std::endl;
3701              numberUnsatisfied_ = 0;
3702              break;
3703#endif
3704              double objValue = solver->getObjValue();
3705              model->setBestSolution(CBC_STRONGSOL,
3706                                     objValue,
3707                                     solver->getColSolution()) ;
3708              model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3709              //double newObjValue = solver->getObjSense()*solver->getObjValue();
3710              //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3711              cutoff =model->getCutoff();
3712            }
3713            solver->markHotStart();
3714          }
3715#ifdef DO_ALL_AT_ROOT
3716          if (!numberNodes)
3717            printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3718                   choice.objectNumber,iStatus,newObjectiveValue,choice.numItersDown,
3719                   choice.downMovement,choice.finishedDown,choice.numIntInfeasDown,
3720                   choice.numObjInfeasDown);
3721#endif
3722         
3723          // repeat the whole exercise, forcing the variable up
3724          choice.possibleBranch->branch();
3725          solver->solveFromHotStart() ;
3726          numberStrongDone++;
3727          numberStrongIterations += solver->getIterationCount();
3728          /*
3729            We now have an estimate of objective degradation that we can use for strong
3730            branching. If we're over the cutoff, the variable is monotone up.
3731            If we actually made it to optimality, check for a solution, and if we have
3732            a good one, call setBestSolution to process it. Note that this may reduce the
3733            cutoff, so we check again to see if we can declare this variable monotone.
3734          */
3735          if (solver->isProvenOptimal())
3736            iStatus=0; // optimal
3737          else if (solver->isIterationLimitReached()
3738                   &&!solver->isDualObjectiveLimitReached())
3739            iStatus=2; // unknown
3740          else
3741            iStatus=1; // infeasible
3742          if (iStatus!=2&&solver->getIterationCount()>
3743              realMaxHotIterations)
3744            numberUnfinished++;
3745          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3746          choice.numItersUp = solver->getIterationCount();
3747          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3748          // Update branching information if wanted
3749          cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3750          if (cbcobj) {
3751            CbcObject * object = cbcobj->object();
3752            assert (object) ;
3753            CbcObjectUpdateData update = object->createUpdateInformation(solver,this,cbcobj);
3754            update.objectNumber_ = choice.objectNumber;
3755            model->addUpdateInformation(update);
3756          } else {
3757            decision->updateInformation( solver,this);
3758          }
3759          if (!iStatus) {
3760            choice.finishedUp = true ;
3761            if (newObjectiveValue>=cutoff) {
3762              objectiveChange = 1.0e100; // say infeasible
3763              numberStrongInfeasible++;
3764            } else {
3765#ifdef TIGHTEN_BOUNDS
3766              // Can we tighten bounds?
3767              if (iColumn<numberColumns&&cutoff<1.0e20
3768                  &&objectiveChange>1.0e-5) {
3769                double value = saveSolution[iColumn];
3770                double up = ceil(value)-value;
3771                double changePer = objectiveChange/(up+1.0e-7);
3772                double distance = (cutoff-objectiveValue_)/changePer;
3773                distance += 1.0e-3;
3774                if (distance<5.0) {
3775                  double newUpper = floor(value+distance);
3776                  if (newUpper<saveUpper[iColumn]) {
3777                    //printf("Could decrease upper bound on %d from %g to %g\n",
3778                    //   iColumn,saveUpper[iColumn],newUpper);
3779                    saveUpper[iColumn]=newUpper;
3780                    solver->setColUpper(iColumn,newUpper);
3781                  }
3782                }
3783                }
3784#endif
3785              // See if integer solution
3786              if (model->feasibleSolution(choice.numIntInfeasUp,
3787                                          choice.numObjInfeasUp)
3788                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
3789#ifdef BONMIN
3790                std::cout<<"Node has become integer feasible"<<std::endl;
3791                numberUnsatisfied_ = 0;
3792                break;
3793#endif
3794                if (auxiliaryInfo->solutionAddsCuts()) {
3795                  needHotStartUpdate=true;
3796                  solver->unmarkHotStart();
3797                }
3798                model->setBestSolution(CBC_STRONGSOL,
3799                                       newObjectiveValue,
3800                                       solver->getColSolution()) ;
3801                if (needHotStartUpdate) {
3802                  model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3803                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3804                  //objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
3805                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3806                  model->feasibleSolution(choice.numIntInfeasDown,
3807                                          choice.numObjInfeasDown);
3808                }
3809                model->setLastHeuristic(NULL);
3810                model->incrementUsed(solver->getColSolution());
3811                cutoff =model->getCutoff();
3812                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3813                  objectiveChange = 1.0e100 ;
3814                  numberStrongInfeasible++;
3815                }
3816              }
3817            }
3818          } else if (iStatus==1) {
3819            objectiveChange = 1.0e100 ;
3820            numberStrongInfeasible++;
3821          } else {
3822            // Can't say much as we did not finish
3823            choice.finishedUp = false ;
3824            numberUnfinished++;
3825          }
3826          choice.upMovement = objectiveChange ;
3827         
3828          // restore bounds
3829          for ( j=0;j<numberColumns;j++) {
3830            if (saveLower[j] != lower[j])
3831              solver->setColLower(j,saveLower[j]);
3832            if (saveUpper[j] != upper[j])
3833              solver->setColUpper(j,saveUpper[j]);
3834          }
3835          if(needHotStartUpdate) {
3836            needHotStartUpdate = false;
3837            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3838            //double newObjValue = solver->getObjSense()*solver->getObjValue();
3839            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3840            //we may again have an integer feasible solution
3841            int numberIntegerInfeasibilities;
3842            int numberObjectInfeasibilities;
3843            if (model->feasibleSolution(
3844                                        numberIntegerInfeasibilities,
3845                                        numberObjectInfeasibilities)) {
3846              double objValue = solver->getObjValue();
3847              model->setBestSolution(CBC_STRONGSOL,
3848                                     objValue,
3849                                     solver->getColSolution()) ;
3850              model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3851              //double newObjValue = solver->getObjSense()*solver->getObjValue();
3852              //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3853              cutoff =model->getCutoff();
3854            }
3855            solver->markHotStart();
3856          }
3857         
3858#ifdef DO_ALL_AT_ROOT
3859          if (!numberNodes)
3860            printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3861                   choice.objectNumber,iStatus,newObjectiveValue,choice.numItersUp,
3862                   choice.upMovement,choice.finishedUp,choice.numIntInfeasUp,
3863                   choice.numObjInfeasUp);
3864#endif
3865        }
3866       
3867        solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit2); 
3868        /*
3869          End of evaluation for this candidate variable. Possibilities are:
3870          * Both sides below cutoff; this variable is a candidate for branching.
3871          * Both sides infeasible or above the objective cutoff: no further action
3872          here. Break from the evaluation loop and assume the node will be purged
3873          by the caller.
3874          * One side below cutoff: Install the branch (i.e., fix the variable). Break
3875          from the evaluation loop and assume the node will be reoptimised by the
3876          caller.
3877        */
3878        // reset
3879        choice.possibleBranch->resetNumberBranchesLeft();
3880        if (choice.upMovement<1.0e100) {
3881          if(choice.downMovement<1.0e100) {
3882            // In case solution coming in was odd
3883            choice.upMovement = CoinMax(0.0,choice.upMovement);
3884            choice.downMovement = CoinMax(0.0,choice.downMovement);
3885#if ZERO_ONE==2
3886            // branch on 0-1 first (temp)
3887            if (fabs(choice.possibleBranch->value())<1.0) {
3888              choice.upMovement *= ZERO_FAKE;
3889              choice.downMovement *= ZERO_FAKE;
3890            }
3891#endif
3892            // feasible - see which best
3893            if (!canSkip) {
3894              if (model->messageHandler()->logLevel()>3) 
3895                printf("sort %g downest %g upest %g ",sort[iDo],downEstimate[iObject],
3896                       upEstimate[iObject]);
3897              model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
3898                << iObject << iColumn
3899                <<choice.downMovement<<choice.numIntInfeasDown
3900                <<choice.upMovement<<choice.numIntInfeasUp
3901                <<choice.possibleBranch->value()
3902                <<CoinMessageEol;
3903            }
3904            int betterWay;
3905            {
3906              CbcBranchingObject * branchObj =
3907                dynamic_cast <CbcBranchingObject *>(branch_) ;
3908              if (branch_)
3909                assert (branchObj);
3910              betterWay = decision->betterBranch(choice.possibleBranch,
3911                                                 branchObj,
3912                                                 choice.upMovement,
3913                                                 choice.numIntInfeasUp ,
3914                                                 choice.downMovement,
3915                                                 choice.numIntInfeasDown );
3916            }
3917            if (betterWay) {
3918              // C) create branching object
3919              if (choiceObject) {
3920                delete branch_;
3921                branch_ = choice.possibleBranch->clone();
3922              } else {
3923                delete branch_;
3924                branch_ = choice.possibleBranch;
3925                choice.possibleBranch=NULL;
3926              }
3927              {
3928                CbcBranchingObject * branchObj =
3929                  dynamic_cast <CbcBranchingObject *>(branch_) ;
3930                assert (branchObj);
3931                //branchObj->way(preferredWay);
3932                branchObj->way(betterWay);
3933              }
3934              bestChoice = choice.objectNumber;
3935              whichChoice = iDo;
3936              if (numberStrong<=1) {
3937                delete ws;
3938                ws=NULL;
3939                break;
3940              }
3941            } else {
3942              if (!choiceObject) {
3943                delete choice.possibleBranch;
3944                choice.possibleBranch=NULL;
3945              }
3946              if (iDo>=2*numberStrong) {
3947                delete ws;
3948                ws=NULL;
3949                break;
3950              }
3951              if (!dynamicObject||dynamicObject->numberTimesUp()>1) {
3952                if (iDo-whichChoice>=numberStrong) {
3953                  if (!choiceObject) {
3954                    delete choice.possibleBranch;
3955                    choice.possibleBranch=NULL;
3956                  }
3957                  break; // give up
3958                }
3959              } else {
3960                if (iDo-whichChoice>=2*numberStrong) {
3961                  delete ws;
3962                  ws=NULL;
3963                  if (!choiceObject) {
3964                    delete choice.possibleBranch;
3965                    choice.possibleBranch=NULL;
3966                  }
3967                  break; // give up
3968                }
3969              }
3970            }
3971          } else {
3972            // up feasible, down infeasible
3973            anyAction=-1;
3974            worstFeasible = CoinMax(worstFeasible,choice.upMovement);
3975            model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
3976              << iObject << iColumn
3977              <<choice.downMovement<<choice.numIntInfeasDown
3978              <<choice.upMovement<<choice.numIntInfeasUp
3979              <<choice.possibleBranch->value()
3980              <<CoinMessageEol;
3981            //printf("Down infeasible for choice %d sequence %d\n",i,
3982            // model->object(choice.objectNumber)->columnNumber());
3983            choice.fix=1;
3984            numberToFix++;
3985            choice.possibleBranch->fix(solver,saveLower,saveUpper,1);
3986            if (!choiceObject) {
3987              choice.possibleBranch=NULL;
3988            } else {
3989              //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3990              choice.possibleBranch=choiceObject;
3991            }
3992            assert(doneHotStart);
3993            solver->unmarkHotStart();
3994            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3995            //double newObjValue = solver->getObjSense()*solver->getObjValue();
3996            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3997            solver->markHotStart();
3998            // may be infeasible (if other way stopped on iterations)
3999            if (!solver->isProvenOptimal()) {
4000              // neither side feasible
4001              anyAction=-2;
4002              if (!choiceObject) {
4003                delete choice.possibleBranch;
4004                choice.possibleBranch=NULL;
4005              }
4006              //printf("Both infeasible for choice %d sequence %d\n",i,
4007              // model->object(choice.objectNumber)->columnNumber());
4008              delete ws;
4009              ws=NULL;
4010              break;
4011            }
4012          }
4013        } else {
4014          if(choice.downMovement<1.0e100) {
4015            // down feasible, up infeasible
4016            anyAction=-1;
4017            worstFeasible = CoinMax(worstFeasible,choice.downMovement);
4018            model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
4019              << iObject << iColumn
4020              <<choice.downMovement<<choice.numIntInfeasDown
4021              <<choice.upMovement<<choice.numIntInfeasUp
4022              <<choice.possibleBranch->value()
4023              <<CoinMessageEol;
4024            choice.fix=-1;
4025            numberToFix++;
4026            choice.possibleBranch->fix(solver,saveLower,saveUpper,-1);
4027            if (!choiceObject) {
4028              choice.possibleBranch=NULL;
4029            } else {
4030              //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
4031              choice.possibleBranch=choiceObject;
4032            }
4033            assert(doneHotStart);
4034            solver->unmarkHotStart();
4035            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
4036            //double newObjValue = solver->getObjSense()*solver->getObjValue();
4037            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
4038            solver->markHotStart();
4039            // may be infeasible (if other way stopped on iterations)
4040            if (!solver->isProvenOptimal()) {
4041              // neither side feasible
4042              anyAction=-2;
4043              if (!choiceObject) {
4044                delete choice.possibleBranch;
4045                choice.possibleBranch=NULL;
4046              }
4047              delete ws;
4048              ws=NULL;
4049              break;
4050            }
4051          } else {
4052            // neither side feasible
4053            anyAction=-2;
4054            if (!choiceObject) {
4055              delete choice.possibleBranch;
4056              choice.possibleBranch=NULL;
4057            }
4058            delete ws;
4059            ws=NULL;
4060            break;
4061          }
4062        }
4063        // Check max time
4064        hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
4065                       model->getDblParam(CbcModel::CbcMaximumSeconds));
4066        if (hitMaxTime) {
4067          // make sure rest are fast
4068          doQuickly=true;
4069          for ( int jDo=iDo+1;jDo<numberToDo;jDo++) {
4070            int iObject = whichObject[iDo];
4071            OsiObject * object = model->modifiableObject(iObject);
4072            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4073              dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4074            if (dynamicObject) 
4075              dynamicObject->setNumberBeforeTrust(0);
4076          }
4077          numberTest=0;
4078        }
4079        if (!choiceObject) {
4080          delete choice.possibleBranch;
4081        }
4082      }
4083      if (model->messageHandler()->logLevel()>3) { 
4084        if (anyAction==-2) {
4085          printf("infeasible\n");
4086        } else if(anyAction==-1) {
4087          printf("%d fixed AND choosing %d iDo %d iChosenWhen %d numberToDo %d\n",numberToFix,bestChoice,
4088                   iDo,whichChoice,numberToDo);
4089        } else {
4090          int iObject = whichObject[whichChoice];
4091          OsiObject * object = model->modifiableObject(iObject);
4092          CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4093            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4094          if (dynamicObject) {
4095            int iColumn = dynamicObject->columnNumber();
4096            printf("choosing %d (column %d) iChosenWhen %d numberToDo %d\n",bestChoice,
4097                   iColumn,whichChoice,numberToDo);
4098          }
4099        }
4100      }
4101      if (doneHotStart) {
4102        // Delete the snapshot
4103        solver->unmarkHotStart();
4104        // back to normal
4105        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
4106        // restore basis
4107        solver->setWarmStart(ws);
4108      }
4109      solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
4110      // Unless infeasible we will carry on
4111      // But we could fix anyway
4112      if (numberToFix&&!hitMaxTime) {
4113        if (anyAction!=-2) {
4114          // apply and take off
4115          bool feasible=true;
4116          // can do quick optimality check
4117          int easy=2;
4118          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
4119          model->resolve(NULL,11,saveSolution,saveLower,saveUpper) ;
4120          //double newObjValue = solver->getObjSense()*solver->getObjValue();
4121          //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
4122          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
4123          feasible = solver->isProvenOptimal();
4124          if (feasible) {
4125            anyAction=0;
4126            // See if candidate still possible
4127            if (branch_) {
4128              const OsiObject * object = model->object(bestChoice);
4129              int preferredWay;
4130              double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
4131              if (!infeasibility) {
4132                // take out
4133                delete branch_;
4134                branch_=NULL;
4135              } else {
4136                CbcBranchingObject * branchObj =
4137                  dynamic_cast <CbcBranchingObject *>(branch_) ;
4138                assert (branchObj);
4139                branchObj->way(preferredWay);
4140              }
4141            }
4142          } else {
4143            anyAction=-2;
4144            finished=true;
4145          }
4146        }
4147      }
4148      // If  fixed then round again
4149      if (!branch_&&anyAction!=-2) {
4150        finished=false;
4151      }
4152      delete ws;
4153    }
4154  }
4155  // update number of strong iterations etc
4156  model->incrementStrongInfo(numberStrongDone,numberStrongIterations,
4157                             anyAction==-2 ? 0:numberToFix,anyAction==-2);
4158  if (model->searchStrategy()==-1) {
4159#ifndef COIN_DEVELOP
4160    if (solver->messageHandler()->logLevel()>1)
4161#endif
4162      printf("%d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
4163             numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
4164             numberNotTrusted);
4165    // decide what to do
4166    int strategy=1;
4167    if (((numberUnfinished*4>numberStrongDone&&
4168          numberStrongInfeasible*40<numberStrongDone)||
4169         numberStrongInfeasible<0)&&model->numberStrong()<10&&model->numberBeforeTrust()<=20&&model->numberObjects()>CoinMax(1000,solver->getNumRows())) {
4170      strategy=2;
4171#ifdef COIN_DEVELOP
4172      //if (model->logLevel()>1)
4173      printf("going to strategy 2\n");
4174#endif
4175      // Weaken
4176      model->setNumberStrong(2);
4177      model->setNumberBeforeTrust(1);
4178      model->synchronizeNumberBeforeTrust();
4179    }
4180    if (numberNodes)
4181      strategy=1;  // should only happen after hot start
4182    model->setSearchStrategy(strategy);
4183  } else if (numberStrongDone) {
4184    //printf("%d strongB, %d iters, %d inf, %d not finished, %d not trusted\n",
4185    //   numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
4186    //   numberNotTrusted);
4187  }
4188  if (model->searchStrategy()==1&&numberNodes>500&&numberNodes<-510) {
4189#ifndef COIN_DEVELOP
4190    if (solver->messageHandler()->logLevel()>1)
4191#endif
4192      printf("after %d nodes - %d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
4193             numberNodes,numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
4194             numberNotTrusted);
4195    // decide what to do
4196    if (numberUnfinished*10>numberStrongDone+1||
4197        !numberStrongInfeasible) {
4198      printf("going to strategy 2\n");
4199      // Weaken
4200      model->setNumberStrong(2);
4201      model->setNumberBeforeTrust(1);
4202      model->synchronizeNumberBeforeTrust();
4203      model->setSearchStrategy(2);
4204    }
4205  }
4206  // Set guessed solution value
4207  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
4208 
4209  /*
4210    Cleanup, then we're finished
4211  */
4212  if (!model->branchingMethod())
4213    delete decision;
4214 
4215  delete choiceObject;
4216  delete [] sort;
4217  delete [] whichObject;
4218#ifdef RANGING
4219  delete [] objectMark;
4220#endif
4221  delete [] saveLower;
4222  delete [] saveUpper;
4223  delete [] upEstimate;
4224  delete [] downEstimate;
4225# ifdef COIN_HAS_CLP
4226  if (osiclp) {
4227    osiclp->setSpecialOptions(saveClpOptions);
4228  }
4229# endif
4230  // restore solution
4231  solver->setColSolution(saveSolution);
4232  model->reserveCurrentSolution(saveSolution);
4233  delete [] saveSolution;
4234  model->setStateOfSearch(saveStateOfSearch);
4235  model->setLogLevel(saveLogLevel);
4236  // delete extra regions
4237  if (usefulInfo.usefulRegion_) {
4238    delete [] usefulInfo.usefulRegion_;
4239    delete [] usefulInfo.indexRegion_;
4240    delete [] usefulInfo.pi_;
4241    usefulInfo.usefulRegion_ = NULL;
4242    usefulInfo.indexRegion_ = NULL;
4243    usefulInfo.pi_ = NULL;
4244  }
4245  useShadow = model->moreSpecialOptions()&7;
4246  if ((useShadow==5&&model->getSolutionCount())||useShadow==6) {
4247    // zap pseudo shadow prices
4248    model->pseudoShadow(-1);
4249    // and switch off
4250    model->setMoreSpecialOptions(model->moreSpecialOptions()&(~1023));
4251  }
4252  return anyAction;
4253}
4254int CbcNode::analyze (CbcModel *model, double * results)
4255{
4256  int i;
4257  int numberIterationsAllowed = model->numberAnalyzeIterations();
4258  OsiSolverInterface * solver = model->solver();
4259  objectiveValue_ = solver->getObjSense()*solver->getObjValue();
4260  double cutoff =model->getCutoff();
4261  const double * lower = solver->getColLower();
4262  const double * upper = solver->getColUpper();
4263  const double * dj = solver->getReducedCost();
4264  int numberObjects = model->numberObjects();
4265  int numberColumns = model->getNumCols();
4266  // Initialize arrays
4267  int numberIntegers = model->numberIntegers();
4268  int * back = new int[numberColumns];
4269  const int * integerVariable = model->integerVariable();
4270  for (i=0;i<numberColumns;i++) 
4271    back[i]=-1;
4272  // What results is
4273  double * newLower = results;
4274  double * objLower = newLower+numberIntegers;
4275  double * newUpper = objLower+numberIntegers;
4276  double * objUpper = newUpper+numberIntegers;
4277  for (i=0;i<numberIntegers;i++) {
4278    int iColumn = integerVariable[i];
4279    back[iColumn]=i;
4280    newLower[i]=0.0;
4281    objLower[i]=-COIN_DBL_MAX;
4282    newUpper[i]=0.0;
4283    objUpper[i]=-COIN_DBL_MAX;
4284  }
4285  double * saveUpper = new double[numberColumns];
4286  double * saveLower = new double[numberColumns];
4287  int anyAction=0;
4288  // Save solution in case heuristics need good solution later
4289 
4290  double * saveSolution = new double[numberColumns];
4291  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
4292  model->reserveCurrentSolution(saveSolution);
4293  for (i=0;i<numberColumns;i++) {
4294    saveLower[i] = lower[i];
4295    saveUpper[i] = upper[i];
4296  }
4297  // Get arrays to sort
4298  double * sort = new double[numberObjects];
4299  int * whichObject = new int[numberObjects];
4300  int numberToFix=0;
4301  int numberToDo=0;
4302  double integerTolerance = 
4303    model->getDblParam(CbcModel::CbcIntegerTolerance);
4304  // point to useful information
4305  OsiBranchingInformation usefulInfo = model->usefulInformation();
4306  // and modify
4307  usefulInfo.depth_=depth_;
4308 
4309  // compute current state
4310  int numberObjectInfeasibilities; // just odd ones
4311  int numberIntegerInfeasibilities;
4312  model->feasibleSolution(
4313                          numberIntegerInfeasibilities,
4314                          numberObjectInfeasibilities);
4315# ifdef COIN_HAS_CLP
4316  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
4317  int saveClpOptions=0;
4318  bool fastIterations = (model->specialOptions()&8)!=0;
4319  if (osiclp&&fastIterations) {
4320    // for faster hot start
4321    saveClpOptions = osiclp->specialOptions();
4322    osiclp->setSpecialOptions(saveClpOptions|8192);
4323  }
4324# else
4325  bool fastIterations = false ;
4326# endif
4327  /*
4328    Scan for branching objects that indicate infeasibility. Choose candidates
4329    using priority as the first criteria, then integer infeasibility.
4330   
4331    The algorithm is to fill the array with a set of good candidates (by
4332    infeasibility) with priority bestPriority.  Finding a candidate with
4333    priority better (less) than bestPriority flushes the choice array. (This
4334    serves as initialization when the first candidate is found.)
4335   
4336  */
4337  numberToDo=0;
4338  for (i=0;i<numberObjects;i++) {
4339    OsiObject * object = model->modifiableObject(i);
4340    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4341      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4342    if(!dynamicObject)
4343      continue;
4344    int preferredWay;
4345    double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
4346    int iColumn = dynamicObject->columnNumber();
4347    if (saveUpper[iColumn]==saveLower[iColumn])
4348      continue;
4349    if (infeasibility)
4350      sort[numberToDo]=-1.0e10-infeasibility;
4351    else
4352      sort[numberToDo]=-fabs(dj[iColumn]);
4353    whichObject[numberToDo++]=i;
4354  }
4355  // Save basis
4356  CoinWarmStart * ws = solver->getWarmStart();
4357  int saveLimit;
4358  solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
4359  int targetIterations = CoinMax(500,numberIterationsAllowed/numberObjects);
4360  if (saveLimit<targetIterations)
4361    solver->setIntParam(OsiMaxNumIterationHotStart,targetIterations); 
4362  // Mark hot start
4363  solver->markHotStart();
4364  // Sort
4365  CoinSort_2(sort,sort+numberToDo,whichObject);
4366  double * currentSolution = model->currentSolution();
4367  double objMin = 1.0e50;
4368  double objMax = -1.0e50;
4369  bool needResolve=false;
4370  int iDo;
4371  for (iDo=0;iDo<numberToDo;iDo++) {
4372    CbcStrongInfo choice;
4373    int iObject = whichObject[iDo];
4374    OsiObject * object = model->modifiableObject(iObject);
4375    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4376      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4377    if (!dynamicObject)
4378      continue;
4379    int iColumn = dynamicObject->columnNumber();
4380    int preferredWay;
4381    object->infeasibility(&usefulInfo,preferredWay);
4382    double value = currentSolution[iColumn];
4383    double nearest = floor(value+0.5);
4384    double lowerValue = floor(value);
4385    bool satisfied=false;
4386    if (fabs(value-nearest)<=integerTolerance||value<saveLower[iColumn]||value>saveUpper[iColumn]) {
4387      satisfied=true;
4388      double newValue;
4389      if (nearest<saveUpper[iColumn]) {
4390        newValue = nearest + 1.0001*integerTolerance;
4391        lowerValue = nearest;
4392      } else {
4393        newValue = nearest - 1.0001*integerTolerance;
4394        lowerValue = nearest-1;
4395      }
4396      currentSolution[iColumn]=newValue;
4397    }
4398    double upperValue = lowerValue+1.0;
4399    //CbcSimpleInteger * obj =
4400    //dynamic_cast <CbcSimpleInteger *>(object) ;
4401    //if (obj) {
4402    //choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
4403    //} else {
4404      CbcObject * obj =
4405        dynamic_cast <CbcObject *>(object) ;
4406      assert (obj);
4407      choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
4408      //}
4409    currentSolution[iColumn]=value;
4410    // Save which object it was
4411    choice.objectNumber=iObject;
4412    choice.numIntInfeasUp = numberUnsatisfied_;
4413    choice.numIntInfeasDown = numberUnsatisfied_;
4414    choice.downMovement = 0.0;
4415    choice.upMovement = 0.0;
4416    choice.numItersDown = 0;
4417    choice.numItersUp = 0;
4418    choice.fix=0; // say not fixed
4419    double objectiveChange ;
4420    double newObjectiveValue=1.0e100;
4421    int j;
4422    // status is 0 finished, 1 infeasible and other
4423    int iStatus;
4424    /*
4425      Try the down direction first. (Specify the initial branching alternative as
4426      down with a call to way(-1). Each subsequent call to branch() performs the
4427      specified branch and advances the branch object state to the next branch
4428      alternative.)
4429    */
4430    choice.possibleBranch->way(-1) ;
4431    choice.possibleBranch->branch() ;
4432    if (fabs(value-lowerValue)>integerTolerance) {
4433      solver->solveFromHotStart() ;
4434      /*
4435        We now have an estimate of objective degradation that we can use for strong
4436        branching. If we're over the cutoff, the variable is monotone up.
4437        If we actually made it to optimality, check for a solution, and if we have
4438        a good one, call setBestSolution to process it. Note that this may reduce the
4439        cutoff, so we check again to see if we can declare this variable monotone.
4440      */
4441      if (solver->isProvenOptimal())
4442        iStatus=0; // optimal
4443      else if (solver->isIterationLimitReached()
4444               &&!solver->isDualObjectiveLimitReached())
4445        iStatus=2; // unknown
4446      else
4447        iStatus=1; // infeasible
4448      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
4449      choice.numItersDown = solver->getIterationCount();
4450      numberIterationsAllowed -= choice.numItersDown;
4451      objectiveChange = newObjectiveValue  - objectiveValue_;
4452      if (!iStatus) {
4453        choice.finishedDown = true ;
4454        if (newObjectiveValue>=cutoff) {
4455          objectiveChange = 1.0e100; // say infeasible
4456        } else {
4457          // See if integer solution
4458          if (model->feasibleSolution(choice.numIntInfeasDown,
4459                                      choice.numObjInfeasDown)
4460              &&model->problemFeasibility()->feasible(model,-1)>=0) {
4461            model->setBestSolution(CBC_STRONGSOL,
4462                                   newObjectiveValue,
4463                                   solver->getColSolution()) ;
4464            model->setLastHeuristic(NULL);
4465            model->incrementUsed(solver->getColSolution());
4466            cutoff =model->getCutoff();
4467            if (newObjectiveValue >= cutoff)    //  *new* cutoff
4468              objectiveChange = 1.0e100 ;
4469          }
4470        }
4471      } else if (iStatus==1) {
4472        objectiveChange = 1.0e100 ;
4473      } else {
4474        // Can't say much as we did not finish
4475        choice.finishedDown = false ;
4476      }
4477      choice.downMovement = objectiveChange ;
4478    }
4479    // restore bounds
4480    for ( j=0;j<numberColumns;j++) {
4481      if (saveLower[j] != lower[j])
4482        solver->setColLower(j,saveLower[j]);
4483      if (saveUpper[j] != upper[j])
4484        solver->setColUpper(j,saveUpper[j]);
4485    }
4486    // repeat the whole exercise, forcing the variable up
4487    choice.possibleBranch->branch();
4488    if (fabs(value-upperValue)>integerTolerance) {
4489      solver->solveFromHotStart() ;
4490      /*
4491        We now have an estimate of objective degradation that we can use for strong
4492        branching. If we're over the cutoff, the variable is monotone up.
4493        If we actually made it to optimality, check for a solution, and if we have
4494        a good one, call setBestSolution to process it. Note that this may reduce the
4495        cutoff, so we check again to see if we can declare this variable monotone.
4496      */
4497      if (solver->isProvenOptimal())
4498        iStatus=0; // optimal
4499      else if (solver->isIterationLimitReached()
4500               &&!solver->isDualObjectiveLimitReached())
4501        iStatus=2; // unknown
4502      else
4503        iStatus=1; // infeasible
4504      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
4505      choice.numItersUp = solver->getIterationCount();
4506      numberIterationsAllowed -= choice.numItersUp;
4507      objectiveChange = newObjectiveValue  - objectiveValue_;
4508      if (!iStatus) {
4509        choice.finishedUp = true ;
4510        if (newObjectiveValue>=cutoff) {
4511          objectiveChange = 1.0e100; // say infeasible
4512        } else {
4513          // See if integer solution
4514          if (model->feasibleSolution(choice.numIntInfeasUp,
4515                                      choice.numObjInfeasUp)
4516              &&model->problemFeasibility()->feasible(model,-1)>=0) {
4517            model->setBestSolution(CBC_STRONGSOL,
4518                                   newObjectiveValue,
4519                                   solver->getColSolution()) ;
4520            model->setLastHeuristic(NULL);
4521            model->incrementUsed(solver->getColSolution());
4522            cutoff =model->getCutoff();
4523            if (newObjectiveValue >= cutoff)    //  *new* cutoff
4524              objectiveChange = 1.0e100 ;
4525          }
4526        }
4527      } else if (iStatus==1) {
4528        objectiveChange = 1.0e100 ;
4529      } else {
4530        // Can't say much as we did not finish
4531        choice.finishedUp = false ;
4532      }
4533      choice.upMovement = objectiveChange ;
4534     
4535      // restore bounds
4536      for ( j=0;j<numberColumns;j++) {
4537        if (saveLower[j] != lower[j])
4538          solver->setColLower(j,saveLower[j]);
4539        if (saveUpper[j] != upper[j])
4540          solver->setColUpper(j,saveUpper[j]);
4541      }
4542    }
4543    // If objective goes above certain amount we can set bound
4544    int jInt = back[iColumn];
4545    newLower[jInt]=upperValue;
4546    if (choice.finishedDown)
4547      objLower[jInt]=choice.downMovement+objectiveValue_;
4548    else
4549      objLower[jInt]=objectiveValue_;
4550    newUpper[jInt]=lowerValue;
4551    if (choice.finishedUp)
4552      objUpper[jInt]=choice.upMovement+objectiveValue_;
4553    else
4554      objUpper[jInt]=objectiveValue_;
4555    objMin = CoinMin(CoinMin(objLower[jInt],objUpper[jInt]),objMin);
4556    /*
4557      End of evaluation for this candidate variable. Possibilities are:
4558      * Both sides below cutoff; this variable is a candidate for branching.
4559      * Both sides infeasible or above the objective cutoff: no further action
4560      here. Break from the evaluation loop and assume the node will be purged
4561      by the caller.
4562      * One side below cutoff: Install the branch (i.e., fix the variable). Break
4563      from the evaluation loop and assume the node will be reoptimised by the
4564      caller.
4565    */
4566    if (choice.upMovement<1.0e100) {
4567      if(choice.downMovement<1.0e100) {
4568        objMax = CoinMax(CoinMax(objLower[jInt],objUpper[jInt]),objMax);
4569        // In case solution coming in was odd
4570        choice.upMovement = CoinMax(0.0,choice.upMovement);
4571        choice.downMovement = CoinMax(0.0,choice.downMovement);
4572        // feasible -
4573        model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
4574          << iObject << iColumn
4575          <<choice.downMovement<<choice.numIntInfeasDown
4576          <<choice.upMovement<<choice.numIntInfeasUp
4577          <<value
4578          <<CoinMessageEol;
4579      } else {
4580        // up feasible, down infeasible
4581        anyAction=-1;
4582        if (!satisfied)
4583          needResolve=true;
4584        choice.fix=1;
4585        numberToFix++;
4586        saveLower[iColumn]=upperValue;
4587        solver->setColLower(iColumn,upperValue);
4588      }
4589    } else {
4590      if(choice.downMovement<1.0e100) {
4591        // down feasible, up infeasible
4592        anyAction=-1;
4593        if (!satisfied)
4594          needResolve=true;
4595        choice.fix=-1;
4596        numberToFix++;
4597        saveUpper[iColumn]=lowerValue;
4598        solver->setColUpper(iColumn,lowerValue);
4599      } else {
4600        // neither side feasible
4601        anyAction=-2;
4602        printf("Both infeasible for choice %d sequence %d\n",i,
4603               model->object(choice.objectNumber)->columnNumber());
4604        delete ws;
4605        ws=NULL;
4606        //solver->writeMps("bad");
4607        numberToFix=-1;
4608        delete choice.possibleBranch;
4609        choice.possibleBranch=NULL;
4610        break;
4611      }
4612    }
4613    delete choice.possibleBranch;
4614    if (numberIterationsAllowed<=0)
4615      break;
4616    //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
4617    //     choice.downMovement,choice.upMovement,value);
4618  }
4619  printf("Best possible solution %g, can fix more if solution of %g found - looked at %d variables in %d iterations\n",
4620         objMin,objMax,iDo,model->numberAnalyzeIterations()-numberIterationsAllowed);
4621  model->setNumberAnalyzeIterations(numberIterationsAllowed);
4622  // Delete the snapshot
4623  solver->unmarkHotStart();
4624  // back to normal
4625  solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
4626  solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
4627  // restore basis
4628  solver->setWarmStart(ws);
4629  delete ws;
4630 
4631  delete [] sort;
4632  delete [] whichObject;
4633  delete [] saveLower;
4634  delete [] saveUpper;
4635  delete [] back;
4636  // restore solution
4637  solver->setColSolution(saveSolution);
4638# ifdef COIN_HAS_CLP
4639  if (osiclp) 
4640    osiclp->setSpecialOptions(saveClpOptions);
4641# endif
4642  model->reserveCurrentSolution(saveSolution);
4643  delete [] saveSolution;
4644  if (needResolve)
4645    solver->resolve();
4646  return numberToFix;
4647}
4648
4649
4650CbcNode::CbcNode(const CbcNode & rhs)
4651  : CoinTreeNode(rhs) 
4652{ 
4653#ifdef CHECK_NODE
4654  printf("CbcNode %x Constructor from rhs %x\n",this,&rhs);
4655#endif
4656  if (rhs.nodeInfo_)
4657    nodeInfo_ = rhs.nodeInfo_->clone();
4658  else
4659    nodeInfo_=NULL;
4660  objectiveValue_=rhs.objectiveValue_;
4661  guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4662  sumInfeasibilities_ = rhs.sumInfeasibilities_;
4663  if (rhs.branch_)
4664    branch_=rhs.branch_->clone();
4665  else
4666    branch_=NULL;
4667  depth_ = rhs.depth_;
4668  numberUnsatisfied_ = rhs.numberUnsatisfied_;
4669  nodeNumber_ = rhs.nodeNumber_;
4670  state_ = rhs.state_;
4671  if (nodeInfo_)
4672    assert ((state_&2)!=0);
4673  else
4674    assert ((state_&2)==0);
4675}
4676
4677CbcNode &
4678CbcNode::operator=(const CbcNode & rhs)
4679{
4680  if (this != &rhs) {
4681    delete nodeInfo_;
4682    if (rhs.nodeInfo_)
4683      nodeInfo_ = rhs.nodeInfo_->clone();
4684    else
4685      nodeInfo_ = NULL;
4686    objectiveValue_=rhs.objectiveValue_;
4687    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4688    sumInfeasibilities_ = rhs.sumInfeasibilities_;
4689    if (rhs.branch_)
4690      branch_=rhs.branch_->clone();
4691    else
4692      branch_=NULL,
4693        depth_ = rhs.depth_;
4694    numberUnsatisfied_ = rhs.numberUnsatisfied_;
4695    nodeNumber_ = rhs.nodeNumber_;
4696    state_ = rhs.state_;
4697    if (nodeInfo_)
4698      assert ((state_&2)!=0);
4699    else
4700      assert ((state_&2)==0);
4701  }
4702  return *this;
4703}
4704CbcNode::~CbcNode ()
4705{
4706#ifdef CHECK_NODE
4707  if (nodeInfo_) {
4708    printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
4709           this,nodeInfo_,nodeInfo_->numberPointingToThis());
4710    //assert(nodeInfo_->numberPointingToThis()>=0);
4711  } else {
4712    printf("CbcNode %x Destructor nodeInfo %x (?)\n",
4713           this,nodeInfo_);
4714  }
4715#endif
4716  if (nodeInfo_) {
4717    // was if (nodeInfo_&&(state_&2)!=0) {
4718    nodeInfo_->nullOwner();
4719    int numberToDelete=nodeInfo_->numberBranchesLeft();
4720    //    CbcNodeInfo * parent = nodeInfo_->parent();
4721    //assert (nodeInfo_->numberPointingToThis()>0);
4722    if (nodeInfo_->decrement(numberToDelete)==0||(state_&2)==0) {
4723      if ((state_&2)==0) 
4724        nodeInfo_->nullParent();
4725      delete nodeInfo_;
4726    } else {
4727      //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,nodeInfo_->parent());
4728      // anyway decrement parent
4729      //if (parent)
4730      ///parent->decrement(1);
4731    }
4732  }
4733  delete branch_;
4734}
4735// Decrement  active cut counts
4736void 
4737CbcNode::decrementCuts(int change)
4738{
4739  if (nodeInfo_)
4740    assert ((state_&2)!=0);
4741  else
4742    assert ((state_&2)==0);
4743  if(nodeInfo_) {
4744    nodeInfo_->decrementCuts(change);
4745  }
4746}
4747void 
4748CbcNode::decrementParentCuts(CbcModel * model, int change)
4749{
4750  if (nodeInfo_)
4751    assert ((state_&2)!=0);
4752  else
4753    assert ((state_&2)==0);
4754  if(nodeInfo_) {
4755    nodeInfo_->decrementParentCuts(model, change);
4756  }
4757}
4758
4759/*
4760  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
4761  in the attached nodeInfo_.
4762*/
4763void
4764CbcNode::initializeInfo()
4765{
4766  assert(nodeInfo_ && branch_) ;
4767  nodeInfo_->initializeInfo(branch_->numberBranches());
4768  assert ((state_&2)!=0);
4769  assert (nodeInfo_->numberBranchesLeft()==
4770          branch_->numberBranchesLeft());
4771}
4772// Nulls out node info
4773void 
4774CbcNode::nullNodeInfo()
4775{
4776  nodeInfo_=NULL;
4777  // say not active
4778  state_ &= ~2;
4779}
4780
4781int
4782CbcNode::branch(OsiSolverInterface * solver)
4783{
4784  double changeInGuessed;
4785  assert (nodeInfo_->numberBranchesLeft()==
4786          branch_->numberBranchesLeft());
4787  if (!solver)
4788    changeInGuessed=branch_->branch();
4789  else
4790    changeInGuessed=branch_->branch(solver);
4791  guessedObjectiveValue_+= changeInGuessed;
4792  //#define PRINTIT
4793#ifdef PRINTIT
4794  int numberLeft = nodeInfo_->numberBranchesLeft();
4795  CbcNodeInfo * parent = nodeInfo_->parent();
4796  int parentNodeNumber = -1;
4797  CbcBranchingObject * object1 = 
4798    dynamic_cast<CbcBranchingObject *>(branch_) ;
4799  //OsiObject * object = object1->
4800  //int sequence = object->columnNumber);
4801  int id=-1;
4802  double value=0.0;
4803  if (object1) {
4804    id = object1->variable();
4805    value = object1->value();
4806  }
4807  printf("id %d value %g objvalue %g\n",id,value,objectiveValue_);
4808  if (parent)
4809    parentNodeNumber = parent->nodeNumber();
4810  printf("Node number %d, %s, way %d, depth %d, parent node number %d\n",
4811         nodeInfo_->nodeNumber(),(numberLeft==2) ? "leftBranch" : "rightBranch",
4812         way(),depth_,parentNodeNumber);
4813  assert (parentNodeNumber!=nodeInfo_->nodeNumber());
4814#endif
4815  return nodeInfo_->branchedOn();
4816}
4817/* Active arm of the attached OsiBranchingObject.
4818   
4819   In the simplest instance, coded -1 for the down arm of the branch, +1 for
4820   the up arm. But see OsiBranchingObject::way()
4821   Use nodeInfo--.numberBranchesLeft_ to see how active
4822*/
4823int 
4824CbcNode::way() const
4825{
4826  if (branch_) {
4827    CbcBranchingObject * obj =
4828      dynamic_cast <CbcBranchingObject *>(branch_) ;
4829    if (obj) {
4830      return obj->way();
4831    } else {
4832      OsiTwoWayBranchingObject * obj2 =
4833        dynamic_cast <OsiTwoWayBranchingObject *>(branch_) ;
4834      assert (obj2);
4835      return obj2->way();
4836    }
4837  } else {
4838    return 0;
4839  }
4840}
4841/* Create a branching object for the node
4842   
4843   The routine scans the object list of the model and selects a set of
4844   unsatisfied objects as candidates for branching. The candidates are
4845   evaluated, and an appropriate branch object is installed.
4846   
4847   The numberPassesLeft is decremented to stop fixing one variable each time
4848   and going on and on (e.g. for stock cutting, air crew scheduling)
4849   
4850   If evaluation determines that an object is monotone or infeasible,
4851   the routine returns immediately. In the case of a monotone object,
4852   the branch object has already been called to modify the model.
4853   
4854   Return value:
4855   <ul>
4856   <li>  0: A branching object has been installed
4857   <li> -1: A monotone object was discovered
4858   <li> -2: An infeasible object was discovered
4859   </ul>
4860   Branch state:
4861   <ul>
4862   <li> -1: start
4863   <li> -1: A monotone object was discovered
4864   <li> -2: An infeasible object was discovered
4865   </ul>
4866*/
4867int 
4868CbcNode::chooseOsiBranch (CbcModel * model,
4869                          CbcNode * lastNode,
4870                          OsiBranchingInformation * usefulInfo,
4871                          int branchState)
4872{
4873  int returnStatus=0;
4874  if (lastNode)
4875    depth_ = lastNode->depth_+1;
4876  else
4877    depth_ = 0;
4878  OsiSolverInterface * solver = model->solver();
4879  objectiveValue_ = solver->getObjValue()*solver->getObjSense();
4880  usefulInfo->objectiveValue_ = objectiveValue_;
4881  usefulInfo->depth_ = depth_;
4882  const double * saveInfoSol = usefulInfo->solution_;
4883  double * saveSolution = new double[solver->getNumCols()];
4884  memcpy(saveSolution,solver->getColSolution(),solver->getNumCols()*sizeof(double));
4885  usefulInfo->solution_ = saveSolution;
4886  OsiChooseVariable * choose = model->branchingMethod()->chooseMethod();
4887  int numberUnsatisfied=-1;
4888  if (branchState<0) {
4889    // initialize
4890    // initialize sum of "infeasibilities"
4891    sumInfeasibilities_ = 0.0;
4892    numberUnsatisfied = choose->setupList(usefulInfo,true);
4893    numberUnsatisfied_ = numberUnsatisfied;
4894    branchState=0;
4895    if (numberUnsatisfied_<0) {
4896      // infeasible
4897      delete [] saveSolution;
4898      return -2;
4899    }
4900  }
4901  // unset best
4902  int best=-1;
4903  choose->setBestObjectIndex(-1);
4904  if (numberUnsatisfied) {
4905    if (branchState>0||!choose->numberOnList()) {
4906      // we need to return at once - don't do strong branching or anything
4907      if (choose->numberOnList()||!choose->numberStrong()) {
4908        best = choose->candidates()[0];
4909        choose->setBestObjectIndex(best);
4910      } else {
4911        // nothing on list - need to try again - keep any solution
4912        numberUnsatisfied = choose->setupList(usefulInfo, false);
4913        numberUnsatisfied_ = numberUnsatisfied;
4914        if (numberUnsatisfied) {
4915          best = choose->candidates()[0];
4916          choose->setBestObjectIndex(best);
4917        }
4918      }
4919    } else {
4920      // carry on with strong branching or whatever
4921      int returnCode = choose->chooseVariable(solver, usefulInfo,true);
4922      // update number of strong iterations etc
4923      model->incrementStrongInfo(choose->numberStrongDone(),choose->numberStrongIterations(),
4924                                 returnCode==-1 ? 0:choose->numberStrongFixed(),returnCode==-1);
4925      if (returnCode>1) {
4926        // has fixed some
4927        returnStatus=-1;
4928      } else if (returnCode==-1) {
4929        // infeasible
4930        returnStatus=-2;
4931      } else if (returnCode==0) {
4932        // normal
4933        returnStatus=0;
4934        numberUnsatisfied=1;
4935      } else {
4936        // ones on list satisfied - double check
4937        numberUnsatisfied = choose->setupList(usefulInfo, false);
4938        numberUnsatisfied_ = numberUnsatisfied;
4939        if (numberUnsatisfied) {
4940          best = choose->candidates()[0];
4941          choose->setBestObjectIndex(best);
4942        }
4943      }
4944    }
4945  } 
4946  delete branch_;
4947  branch_ = NULL;
4948  guessedObjectiveValue_ = COIN_DBL_MAX;//objectiveValue_; // for now
4949  if (!returnStatus) {
4950    if (numberUnsatisfied) {
4951      // create branching object
4952      const OsiObject * obj = model->solver()->object(choose->bestObjectIndex());
4953      //const OsiSolverInterface * solver = usefulInfo->solver_;
4954      branch_ = obj->createBranch(model->solver(),usefulInfo,obj->whichWay());
4955    }
4956  }
4957  usefulInfo->solution_=saveInfoSol;
4958  delete [] saveSolution;
4959  // may have got solution
4960  if (choose->goodSolution()
4961      &&model->problemFeasibility()->feasible(model,-1)>=0) {
4962    // yes
4963    double objValue = choose->goodObjectiveValue();
4964    model->setBestSolution(CBC_STRONGSOL,
4965                           objValue,
4966                           choose->goodSolution()) ;
4967    model->setLastHeuristic(NULL);
4968    model->incrementUsed(choose->goodSolution());
4969    choose->clearGoodSolution();
4970  }
4971  return returnStatus;
4972}
4973int 
4974CbcNode::chooseClpBranch (CbcModel * model,
4975                          CbcNode * lastNode)
4976{ 
4977  assert(lastNode);
4978  depth_ = lastNode->depth_+1;
4979  delete branch_;
4980  branch_=NULL;
4981  OsiSolverInterface * solver = model->solver();
4982  //double saveObjectiveValue = solver->getObjValue();
4983  //double objectiveValue = CoinMax(solver->getObjSense()*saveObjectiveValue,objectiveValue_);
4984  const double * lower = solver->getColLower();
4985  const double * upper = solver->getColUpper();
4986  // point to useful information
4987  OsiBranchingInformation usefulInfo = model->usefulInformation();
4988  // and modify
4989  usefulInfo.depth_=depth_;
4990  int i;
4991  //bool beforeSolution = model->getSolutionCount()==0;
4992  int numberObjects = model->numberObjects();
4993  int numberColumns = model->getNumCols();
4994  double * saveUpper = new double[numberColumns];
4995  double * saveLower = new double[numberColumns];
4996 
4997  // Save solution in case heuristics need good solution later
4998 
4999  double * saveSolution = new double[numberColumns];
5000  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
5001  model->reserveCurrentSolution(saveSolution);
5002  for (i=0;i<numberColumns;i++) {
5003    saveLower[i] = lower[i];
5004    saveUpper[i] = upper[i];
5005  }
5006  // Save basis
5007  CoinWarmStart * ws = solver->getWarmStart();
5008  numberUnsatisfied_ = 0;
5009  // initialize sum of "infeasibilities"
5010  sumInfeasibilities_ = 0.0;
5011  // Note looks as if off end (hidden one)
5012  OsiObject * object = model->modifiableObject(numberObjects);
5013  CbcGeneralDepth * thisOne = dynamic_cast <CbcGeneralDepth *> (object);
5014  assert (thisOne);
5015  OsiClpSolverInterface * clpSolver
5016    = dynamic_cast<OsiClpSolverInterface *> (solver);
5017  assert (clpSolver);
5018  ClpSimplex * simplex = clpSolver->getModelPtr();
5019  int preferredWay;
5020  double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
5021  if (thisOne->whichSolution()>=0) {
5022    ClpNode * nodeInfo = thisOne->nodeInfo(thisOne->whichSolution());
5023    nodeInfo->applyNode(simplex,2);
5024    int saveLogLevel = simplex->logLevel();
5025    simplex->setLogLevel(0);
5026    simplex->dual();
5027    simplex->setLogLevel(saveLogLevel);
5028    double cutoff = model->getCutoff();
5029    bool goodSolution=true;
5030    if (simplex->status()) {
5031      //simplex->writeMps("bad7.mps",2);
5032      if (nodeInfo->objectiveValue()>cutoff-1.0e-2)
5033        goodSolution=false;
5034      else
5035        assert (!simplex->status());
5036    }
5037    if (goodSolution) {
5038      double newObjectiveValue = solver->getObjSense()*solver->getObjValue();
5039      // See if integer solution
5040      int numInf;
5041      int numInf2;
5042      bool gotSol = model->feasibleSolution(numInf,numInf2);
5043      if (!gotSol) {
5044        printf("numinf %d\n",numInf);
5045        double * sol = simplex->primalColumnSolution();
5046        for (int i=0;i<numberColumns;i++) {
5047          if (simplex->isInteger(i)) {
5048            double value = floor(sol[i]+0.5);
5049            if (fabs(value-sol[i])>1.0e-7) {
5050              printf("%d value %g\n",i,sol[i]);
5051              if (fabs(value-sol[i])<1.0e-3) {
5052                sol[i]=value;
5053              }
5054            }
5055          }
5056        }
5057        simplex->writeMps("bad8.mps",2);
5058        bool gotSol = model->feasibleSolution(numInf,numInf2);
5059        if (!gotSol) 
5060          assert (gotSol);
5061      }
5062      model->setBestSolution(CBC_STRONGSOL,
5063                             newObjectiveValue,
5064                             solver->getColSolution()) ;
5065      model->setLastHeuristic(NULL);
5066      model->incrementUsed(solver->getColSolution());
5067    }
5068  }
5069  // restore bounds
5070  { for (int j=0;j<numberColumns;j++) {
5071      if (saveLower[j] != lower[j])
5072        solver->setColLower(j,saveLower[j]);
5073      if (saveUpper[j] != upper[j])
5074        solver->setColUpper(j,saveUpper[j]);
5075    }
5076  }
5077  // restore basis
5078  solver->setWarmStart(ws);
5079  delete ws;
5080  int anyAction;
5081  //#define CHECK_PATH
5082#ifdef CHECK_PATH
5083  extern int gotGoodNode_Z;
5084  if (gotGoodNode_Z>=0)
5085    printf("good node %d %g\n",gotGoodNode_Z,infeasibility);
5086#endif
5087  if (infeasibility>0.0) {
5088    if (infeasibility==COIN_DBL_MAX) {
5089      anyAction=-2; // infeasible
5090    } else {
5091      branch_=thisOne->createCbcBranch(solver,&usefulInfo,preferredWay);
5092      // Set to firat one (and change when re-pushing)
5093      CbcGeneralBranchingObject * branch = 
5094        dynamic_cast <CbcGeneralBranchingObject *> (branch_);
5095      branch->state(objectiveValue_,sumInfeasibilities_,
5096                    numberUnsatisfied_,0);
5097      branch->setNode(this);
5098      anyAction=0;
5099    }
5100  } else {
5101    anyAction=-1;
5102  }
5103#ifdef CHECK_PATH
5104  gotGoodNode_Z=-1;
5105#endif
5106  // Set guessed solution value
5107  guessedObjectiveValue_ = objectiveValue_+1.0e-5;
5108  delete [] saveLower;
5109  delete [] saveUpper;
5110 
5111  // restore solution
5112  solver->setColSolution(saveSolution);
5113  delete [] saveSolution;
5114  return anyAction;
5115}
5116/* Double checks in case node can change its mind!
5117   Returns objective value
5118   Can change objective etc */
5119double
5120CbcNode::checkIsCutoff(double cutoff)
5121{
5122  branch_->checkIsCutoff(cutoff);
5123  return objectiveValue_;
5124}
Note: See TracBrowser for help on using the repository browser.