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

Last change on this file since 778 was 778, checked in by andreasw, 12 years ago

setting HUGE dummy value for guessedObjectiveValue_ in osiChooseBranch; that way we can use guessed objectives in Bonmin

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