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

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

for lou and callable cbcmain

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