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

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

for osibranching

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