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

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

for osi methods

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