source: trunk/CbcNode.cpp @ 268

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