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

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

probably breaking everything

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