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

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

for nonlinear and start moving to OsiTree?
afor n

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