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

Last change on this file since 1419 was 1419, checked in by forrest, 11 years ago

take out a few redundant lines

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