source: branches/devel/Cbc/src/CbcNode.cpp @ 445

Last change on this file since 445 was 445, checked in by forrest, 13 years ago

for OsiSOS and add more stats to cut generators

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