source: trunk/CbcNode.cpp @ 267

Last change on this file since 267 was 267, checked in by forrest, 14 years ago

for nonlinear

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