source: trunk/Cbc/src/CbcNode.cpp @ 833

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

changes to heuristics and allow collection of more statistics

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