source: stable/2.4/Cbc/src/CbcNode.cpp @ 1271

Last change on this file since 1271 was 1271, checked in by forrest, 10 years ago

Creating new stable branch 2.4 from trunk (rev 1270)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 171.3 KB
Line 
1/* $Id: CbcNode.cpp 1271 2009-11-05 15:57:25Z 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      } else if (iPass==1) {
3111        // looks like a solution - get paranoid
3112        bool roundAgain=false;
3113        // get basis
3114        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
3115        if (!ws)
3116          break;
3117        double tolerance;
3118        solver->getDblParam(OsiPrimalTolerance,tolerance);
3119        for (i=0;i<numberColumns;i++) {
3120          double value = saveSolution[i];
3121          if (value<lower[i]-tolerance) {
3122            saveSolution[i]=lower[i];
3123            roundAgain=true;
3124            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
3125          } else if (value>upper[i]+tolerance) {
3126            saveSolution[i]=upper[i];
3127            roundAgain=true;
3128            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
3129          } 
3130        }
3131        if (roundAgain) {
3132          // restore basis
3133          solver->setWarmStart(ws);
3134          solver->setColSolution(saveSolution);
3135          delete ws;
3136          bool takeHint;
3137          OsiHintStrength strength;
3138          solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
3139          solver->setHintParam(OsiDoDualInResolve,false,OsiHintDo) ;
3140          model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3141          //double newObjValue = solver->getObjSense()*solver->getObjValue();
3142          //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3143          solver->setHintParam(OsiDoDualInResolve,takeHint,strength) ;
3144          if (!solver->isProvenOptimal()) {
3145            // infeasible
3146            anyAction=-2;
3147            break;
3148          }
3149        } else {
3150          delete ws;
3151          break;
3152        }
3153      }
3154    }
3155    if (anyAction==-2) {
3156      break;
3157    }
3158    // skip if solution
3159    if (!numberUnsatisfied_)
3160      break;
3161    int skipAll = (numberNotTrusted==0||numberToDo==1) ? 1 : 0;
3162    bool doneHotStart=false;
3163    //DEPRECATED_STRATEGYint searchStrategy = saveSearchStrategy>=0 ? (saveSearchStrategy%10) : -1;
3164    int searchStrategy = model->searchStrategy();
3165    // But adjust depending on ratio of iterations
3166    if (searchStrategy>0) {
3167      if (numberBeforeTrust>=5&&numberBeforeTrust<=10) {
3168        if (searchStrategy!=2) {
3169          assert (searchStrategy==1);
3170          if (depth_>5) {
3171            int numberIterations = model->getIterationCount();
3172            int numberStrongIterations = model->numberStrongIterations();
3173            if (numberStrongIterations>numberIterations+10000) {
3174              searchStrategy=2;
3175              skipAll=1;
3176            } else if (numberStrongIterations*4+1000<numberIterations) {
3177              searchStrategy=3;
3178              skipAll=0;
3179            }
3180          } else {
3181            searchStrategy=3;
3182            skipAll=0;
3183          }
3184        }
3185      }
3186    }
3187    // worth trying if too many times
3188    // Save basis
3189    CoinWarmStart * ws = NULL;
3190    // save limit
3191    int saveLimit=0;
3192    solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
3193    if (!numberPassesLeft)
3194      skipAll=1;
3195    if (!skipAll) {
3196      ws = solver->getWarmStart();
3197      int limit=100;
3198      if (!saveStateOfSearch&&saveLimit<limit&&saveLimit==100)
3199        solver->setIntParam(OsiMaxNumIterationHotStart,limit); 
3200    }
3201    // Say which one will be best
3202    int whichChoice=0;
3203    int bestChoice;
3204    if (iBestGot>=0)
3205      bestChoice=iBestGot;
3206    else
3207      bestChoice=iBestNot;
3208    assert (bestChoice>=0);
3209    // If we have hit max time don't do strong branching
3210    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
3211                        model->getDblParam(CbcModel::CbcMaximumSeconds));
3212    // also give up if we are looping round too much
3213    if (hitMaxTime||numberPassesLeft<=0||useShadow==2) {
3214      int iObject = whichObject[bestChoice];
3215      OsiObject * object = model->modifiableObject(iObject);
3216      int preferredWay;
3217      object->infeasibility(&usefulInfo,preferredWay);
3218      CbcObject * obj =
3219        dynamic_cast <CbcObject *>(object) ;
3220      assert (obj);
3221      branch_=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3222      {
3223        CbcBranchingObject * branchObj =
3224          dynamic_cast <CbcBranchingObject *>(branch_) ;
3225        assert (branchObj);
3226        branchObj->way(preferredWay);
3227      }
3228      delete ws;
3229      ws=NULL;
3230      break;
3231    } else {
3232      // say do fast
3233      int easy=1;
3234      if (!skipAll)
3235        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
3236      int iDo;
3237#ifdef RANGING
3238      bool useRanging = model->allDynamic()&&!skipAll;
3239      if (useRanging) {
3240        double currentObjective = solver->getObjValue()*solver->getObjSense();
3241        double gap = cutoff-currentObjective;
3242        // relax a bit
3243        gap *= 1.0000001;
3244        gap = CoinMax(1.0e-5,gap);
3245        // off penalties if too much
3246        double needed = neededPenalties;
3247        needed *= numberRows;
3248        if (numberNodes) {
3249          if (needed>1.0e6) {
3250            neededPenalties=0;
3251          } else if (gap<1.0e5) {
3252            // maybe allow some not needed
3253            int extra = static_cast<int> ((1.0e6-needed)/numberRows);
3254            int nStored = numberObjects-optionalPenalties;
3255            extra = CoinMin(extra,nStored);
3256            for (int i=0;i<extra;i++) {
3257              objectMark[neededPenalties]=objectMark[optionalPenalties+i];
3258                which[neededPenalties++]=which[optionalPenalties+i];;
3259            }
3260          }
3261        }
3262        if (osiclp&&neededPenalties) {
3263          assert (!doneHotStart);
3264          xPen += neededPenalties;
3265          which--;
3266          which[0]=neededPenalties;
3267          osiclp->passInRanges(which);
3268          // Mark hot start and get ranges
3269          if (kPass) {
3270            // until can work out why solution can go funny
3271            int save = osiclp->specialOptions();
3272            osiclp->setSpecialOptions(save|256);
3273            solver->markHotStart();
3274            osiclp->setSpecialOptions(save);
3275          } else {
3276            solver->markHotStart();
3277          }
3278          doneHotStart=true;
3279          xMark++;
3280          kPass++;
3281          osiclp->passInRanges(NULL);
3282          const double * downCost=osiclp->upRange();
3283          const double * upCost=osiclp->downRange();
3284          bool problemFeasible=true;
3285          int numberFixed=0;
3286          for (int i=0;i<neededPenalties;i++) {
3287            int j = objectMark[i];
3288            int iObject = whichObject[j];
3289            OsiObject * object = model->modifiableObject(iObject);
3290            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3291              dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3292            // Use this object's numberBeforeTrust
3293            int numberBeforeTrust = dynamicObject->numberBeforeTrust();
3294            int iSequence=dynamicObject->columnNumber();
3295            double value = saveSolution[iSequence];
3296            value -= floor(value);
3297            double upPenalty = CoinMin(upCost[i],1.0e110)*(1.0-value);
3298            double downPenalty = CoinMin(downCost[i],1.0e110)*value;
3299            int numberThisDown = dynamicObject->numberTimesDown();
3300            int numberThisUp = dynamicObject->numberTimesUp();
3301            if (!numberBeforeTrust) {
3302              // override
3303              downEstimate[iObject]=downPenalty;
3304              upEstimate[iObject]=upPenalty;
3305              double min1 = CoinMin(downEstimate[iObject],
3306                                    upEstimate[iObject]);
3307              double max1 = CoinMax(downEstimate[iObject],
3308                                    upEstimate[iObject]);
3309              min1 = 0.8*min1+0.2*max1;
3310              sort[j] = - min1;
3311            } else if (numberThisDown<numberBeforeTrust||
3312                       numberThisUp<numberBeforeTrust) {
3313              double invTrust = 1.0/static_cast<double> (numberBeforeTrust);
3314              if (numberThisDown<numberBeforeTrust) {
3315                double fraction = numberThisDown*invTrust;
3316                downEstimate[iObject] = fraction*downEstimate[iObject]+(1.0-fraction)*downPenalty;
3317              }
3318              if (numberThisUp<numberBeforeTrust) {
3319                double fraction = numberThisUp*invTrust;
3320                upEstimate[iObject] = fraction*upEstimate[iObject]+(1.0-fraction)*upPenalty;
3321              }
3322              double min1 = CoinMin(downEstimate[iObject],
3323                                    upEstimate[iObject]);
3324              double max1 = CoinMax(downEstimate[iObject],
3325                                    upEstimate[iObject]);
3326              min1 = 0.8*min1+0.2*max1;
3327              min1 *= 10.0;
3328              if (!(numberThisDown+numberThisUp))
3329                min1 *= 100.0;
3330              sort[j] = - min1;
3331            }
3332            if (CoinMax(downPenalty,upPenalty)>gap) {
3333              printf("gap %g object %d has down range %g, up %g\n",
3334                     gap,i,downPenalty,upPenalty);
3335              //sort[j] -= 1.0e50; // make more likely to be chosen
3336              int number;
3337              if (downPenalty>gap) {
3338                number = dynamicObject->numberTimesDown();
3339                if (upPenalty>gap) 
3340                  problemFeasible=false;
3341                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver,&usefulInfo,1);
3342                //branch->fix(solver,saveLower,saveUpper,1);
3343                delete branch;
3344              } else {
3345                number = dynamicObject->numberTimesUp();
3346                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver,&usefulInfo,1);
3347                //branch->fix(solver,saveLower,saveUpper,-1);
3348                delete branch;
3349              }
3350              if(number>=numberBeforeTrust)
3351                dynamicObject->setNumberBeforeTrust(number+1);
3352              numberFixed++;
3353            }
3354            if (!numberNodes)
3355              printf("%d pen down ps %g -> %g up ps %g -> %g\n",
3356                     iObject,downPenalty,downPenalty,upPenalty,upPenalty);
3357          }
3358          if(numberFixed&&problemFeasible) {
3359            assert(doneHotStart);
3360            solver->unmarkHotStart();
3361            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3362            //double newObjValue = solver->getObjSense()*solver->getObjValue();
3363            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3364            solver->markHotStart();
3365            problemFeasible = solver->isProvenOptimal();
3366          }
3367          if(!problemFeasible) {
3368            fprintf(stdout,"both ways infeas on ranging - code needed\n");
3369            anyAction=-2;
3370            if (!choiceObject) {
3371              delete choice.possibleBranch;
3372              choice.possibleBranch=NULL;
3373            }
3374            //printf("Both infeasible for choice %d sequence %d\n",i,
3375            // model->object(choice.objectNumber)->columnNumber());
3376            // Delete the snapshot
3377            solver->unmarkHotStart();
3378            // back to normal
3379            solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3380            // restore basis
3381            solver->setWarmStart(ws);
3382            doneHotStart=false;
3383            delete ws;
3384            ws=NULL;
3385            break;
3386          }
3387        }
3388      }
3389#endif          /* RANGING */
3390      {
3391        int numberIterations = model->getIterationCount();
3392        //numberIterations += (model->numberExtraIterations()>>2);
3393        const int * strongInfo = model->strongInfo();
3394        //int numberDone = strongInfo[0]-strongInfo[3];
3395        int numberFixed = strongInfo[1]-strongInfo[4];
3396        int numberInfeasible = strongInfo[2]-strongInfo[5];
3397        assert (!strongInfo[3]);
3398        assert (!strongInfo[4]);
3399        assert (!strongInfo[5]);
3400        int numberStrongIterations = model->numberStrongIterations();
3401        int numberRows = solver->getNumRows();
3402        if (numberStrongIterations>numberIterations+CoinMin(100,10*numberRows)&&depth_>=4&&numberNodes>100) {
3403          if (20*numberInfeasible+4*numberFixed<numberNodes) {
3404            // Say never do
3405            skipAll=-1;
3406          }
3407        } 
3408      }
3409      // make sure best will be first
3410      if (iBestGot>=0)
3411        sort[iBestGot]=-COIN_DBL_MAX;
3412      // Actions 0 - exit for repeat, 1 resolve and try old choice,2 exit for continue
3413      if (anyAction)
3414        numberToDo=0; // skip as we will be trying again
3415      // Sort
3416      CoinSort_2(sort,sort+numberToDo,whichObject);
3417      // Change in objective opposite infeasible
3418      double worstFeasible=0.0;
3419      // Just first if strong off
3420      if (!numberStrong)
3421        numberToDo=CoinMin(numberToDo,1);
3422      if (searchStrategy==2)
3423        numberToDo=CoinMin(numberToDo,20);
3424      iDo=0;
3425      int saveLimit2;
3426      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit2);
3427      bool doQuickly = false; // numberToDo>2*numberStrong;
3428      if (searchStrategy==2)
3429        doQuickly=true;
3430      int numberTest=numberNotTrusted>0 ? numberStrong : (numberStrong+1)/2;
3431      if (searchStrategy==3) {
3432        // Previously decided we need strong
3433        doQuickly=false;
3434        numberTest = numberStrong;
3435      }
3436      // Try nearly always off
3437      if (skipAll>=0) {
3438        if (searchStrategy<2) {
3439          //if ((numberNodes%20)!=0) {
3440            if ((model->specialOptions()&8)==0) {
3441              numberTest=0;
3442              doQuickly=true;
3443            }
3444            //} else {
3445            //doQuickly=false;
3446            //numberTest=2*numberStrong;
3447            //skipAll=0;
3448            //}
3449        }
3450      } else {
3451        // Just take first
3452        doQuickly=true;
3453        numberTest=1;
3454      }
3455      int testDepth = (skipAll>=0) ? 8 : 4;
3456      if (depth_<testDepth&&numberStrong) {
3457        if (searchStrategy!=2) {
3458          doQuickly=false;
3459          int numberRows = solver->getNumRows();
3460          // whether to do this or not is important - think
3461          if (numberRows<300||numberRows+numberColumns<2500) {
3462            if (depth_<7)
3463              numberStrong = CoinMin(3*numberStrong,numberToDo);
3464            if (!depth_) 
3465              numberStrong=CoinMin(6*numberStrong,numberToDo);
3466          }
3467          numberTest=numberStrong;
3468          skipAll=0;
3469        }
3470      }
3471      // Do at least 5 strong
3472      if (numberColumns<1000&&(depth_<15||numberNodes<1000000))
3473        numberTest = CoinMax(numberTest,5);
3474      if ((model->specialOptions()&8)==0) {
3475        if (skipAll) {
3476          numberTest=0;
3477          doQuickly=true;
3478        }
3479      } else {
3480        // do 5 as strong is fixing
3481        numberTest = CoinMax(numberTest,5);
3482      }
3483      // see if switched off
3484      if (skipAll<0) {
3485        numberTest=0;
3486        doQuickly=true;
3487      }
3488      int realMaxHotIterations=999999;
3489      if (skipAll<0)
3490        numberToDo=1;
3491#ifdef DO_ALL_AT_ROOT
3492      if (!numberNodes) {
3493        printf("DOX %d test %d unsat %d - nobj %d\n",
3494               numberToDo,numberTest,numberUnsatisfied_,
3495               numberObjects);
3496        numberTest=numberToDo;
3497      }
3498#endif
3499      for ( iDo=0;iDo<numberToDo;iDo++) {
3500        int iObject = whichObject[iDo];
3501        OsiObject * object = model->modifiableObject(iObject);
3502        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3503          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3504        int iColumn = dynamicObject ? dynamicObject->columnNumber() : numberColumns+iObject;
3505        int preferredWay;
3506        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
3507        // may have become feasible
3508        if (!infeasibility)
3509          continue;
3510#ifndef NDEBUG
3511        if (iColumn<numberColumns) {
3512            const double * solution = model->testSolution();
3513            assert (saveSolution[iColumn]==solution[iColumn]);
3514          }
3515#endif
3516        CbcSimpleInteger * obj =
3517          dynamic_cast <CbcSimpleInteger *>(object) ;
3518        if (obj) {
3519          if (choiceObject) {
3520            obj->fillCreateBranch(choiceObject,&usefulInfo,preferredWay);
3521            choiceObject->setObject(dynamicObject);
3522          } else {
3523            choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3524          }
3525        } else {
3526          CbcObject * obj =
3527            dynamic_cast <CbcObject *>(object) ;
3528          assert (obj);
3529          choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3530        }
3531        // Save which object it was
3532        choice.objectNumber=iObject;
3533        choice.numIntInfeasUp = numberUnsatisfied_;
3534        choice.numIntInfeasDown = numberUnsatisfied_;
3535        choice.upMovement = upEstimate[iObject];
3536        choice.downMovement = downEstimate[iObject];
3537        assert (choice.upMovement>=0.0);
3538        assert (choice.downMovement>=0.0);
3539        choice.fix=0; // say not fixed
3540        // see if can skip strong branching
3541        int canSkip = choice.possibleBranch->fillStrongInfo(choice);
3542        if ((numberTest<=0||skipAll)) {
3543          if (iDo>20) {
3544#ifdef DO_ALL_AT_ROOT
3545            if (!numberNodes)
3546              printf("DOY test %d - iDo %d\n",
3547                     numberTest,iDo);
3548#endif
3549            if (!choiceObject) {
3550              delete choice.possibleBranch;
3551              choice.possibleBranch=NULL;
3552            }
3553            break; // give up anyway
3554          }
3555        }
3556        if (model->messageHandler()->logLevel()>3&&numberBeforeTrust&&dynamicObject) 
3557          dynamicObject->print(1,choice.possibleBranch->value());
3558        if (skipAll<0)
3559          canSkip=true;
3560        if (!canSkip) {
3561          if (!doneHotStart) {
3562            // Mark hot start
3563            doneHotStart=true;
3564            solver->markHotStart();
3565            xMark++;
3566          }
3567          numberTest--;
3568          // just do a few
3569          if (searchStrategy==2)
3570            solver->setIntParam(OsiMaxNumIterationHotStart,10); 
3571          double objectiveChange ;
3572          double newObjectiveValue=1.0e100;
3573          int j;
3574          // status is 0 finished, 1 infeasible and other
3575          int iStatus;
3576          /*
3577            Try the down direction first. (Specify the initial branching alternative as
3578            down with a call to way(-1). Each subsequent call to branch() performs the
3579            specified branch and advances the branch object state to the next branch
3580            alternative.)
3581          */
3582          choice.possibleBranch->way(-1) ;
3583          choice.possibleBranch->branch() ;
3584          solver->solveFromHotStart() ;
3585          bool needHotStartUpdate=false;
3586          numberStrongDone++;
3587          numberStrongIterations += solver->getIterationCount();
3588          /*
3589            We now have an estimate of objective degradation that we can use for strong
3590            branching. If we're over the cutoff, the variable is monotone up.
3591            If we actually made it to optimality, check for a solution, and if we have
3592            a good one, call setBestSolution to process it. Note that this may reduce the
3593            cutoff, so we check again to see if we can declare this variable monotone.
3594          */
3595          if (solver->isProvenOptimal())
3596            iStatus=0; // optimal
3597          else if (solver->isIterationLimitReached()
3598                   &&!solver->isDualObjectiveLimitReached())
3599            iStatus=2; // unknown
3600          else
3601            iStatus=1; // infeasible
3602          if (iStatus!=2&&solver->getIterationCount()>
3603              realMaxHotIterations)
3604            numberUnfinished++;
3605          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3606          choice.numItersDown = solver->getIterationCount();
3607          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3608          // Update branching information if wanted
3609          CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3610          if (cbcobj) {
3611            CbcObject * object = cbcobj->object();
3612            assert (object) ;
3613            CbcObjectUpdateData update = object->createUpdateInformation(solver,this,cbcobj);
3614            update.objectNumber_ = choice.objectNumber;
3615            model->addUpdateInformation(update);
3616          } else {
3617            decision->updateInformation( solver,this);
3618          }
3619          if (!iStatus) {
3620            choice.finishedDown = true ;
3621            if (newObjectiveValue>=cutoff) {
3622              objectiveChange = 1.0e100; // say infeasible
3623              numberStrongInfeasible++;
3624            } else {
3625              //#define TIGHTEN_BOUNDS
3626#ifdef TIGHTEN_BOUNDS
3627              // Can we tighten bounds?
3628              if (iColumn<numberColumns&&cutoff<1.0e20
3629                  &&objectiveChange>1.0e-5) {
3630                double value = saveSolution[iColumn];
3631                double down = value-floor(value);
3632                double changePer = objectiveChange/(down+1.0e-7);
3633                double distance = (cutoff-objectiveValue_)/changePer;
3634                distance += 1.0e-3;
3635                if (distance<5.0) {
3636                  double newLower = ceil(value-distance);
3637                  if (newLower>saveLower[iColumn]) {
3638                    //printf("Could increase lower bound on %d from %g to %g\n",
3639                    //   iColumn,saveLower[iColumn],newLower);
3640                    saveLower[iColumn]=newLower;
3641                    solver->setColLower(iColumn,newLower);
3642                  }
3643                }
3644                }
3645#endif
3646              // See if integer solution
3647              if (model->feasibleSolution(choice.numIntInfeasDown,
3648                                          choice.numObjInfeasDown)
3649                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
3650                if (auxiliaryInfo->solutionAddsCuts()) {
3651                  needHotStartUpdate=true;
3652                  solver->unmarkHotStart();
3653                }
3654                model->setBestSolution(CBC_STRONGSOL,
3655                                       newObjectiveValue,
3656                                       solver->getColSolution()) ;
3657                if (needHotStartUpdate) {
3658                  model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3659                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3660                  //objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
3661                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3662                  model->feasibleSolution(choice.numIntInfeasDown,
3663                                          choice.numObjInfeasDown);
3664                }
3665                model->setLastHeuristic(NULL);
3666                model->incrementUsed(solver->getColSolution());
3667                cutoff =model->getCutoff();
3668                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3669                  objectiveChange = 1.0e100 ;
3670                  numberStrongInfeasible++;
3671                }
3672              }
3673            }
3674          } else if (iStatus==1) {
3675            objectiveChange = 1.0e100 ;
3676            numberStrongInfeasible++;
3677          } else {
3678            // Can't say much as we did not finish
3679            choice.finishedDown = false ;
3680            numberUnfinished++;
3681          }
3682          choice.downMovement = objectiveChange ;
3683         
3684          // restore bounds
3685          for ( j=0;j<numberColumns;j++) {
3686            if (saveLower[j] != lower[j])
3687              solver->setColLower(j,saveLower[j]);
3688            if (saveUpper[j] != upper[j])
3689              solver->setColUpper(j,saveUpper[j]);
3690          }
3691          if(needHotStartUpdate) {
3692            needHotStartUpdate = false;
3693            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3694            //double newObjValue = solver->getObjSense()*solver->getObjValue();
3695            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3696            //we may again have an integer feasible solution
3697            int numberIntegerInfeasibilities;
3698            int numberObjectInfeasibilities;
3699            if (model->feasibleSolution(
3700                                        numberIntegerInfeasibilities,
3701                                        numberObjectInfeasibilities)) {
3702#ifdef BONMIN
3703              //In this case node has become integer feasible, let us exit the loop
3704              std::cout<<"Node has become integer feasible"<<std::endl;
3705              numberUnsatisfied_ = 0;
3706              break;
3707#endif
3708              double objValue = solver->getObjValue();
3709              model->setBestSolution(CBC_STRONGSOL,
3710                                     objValue,
3711                                     solver->getColSolution()) ;
3712              model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3713              //double newObjValue = solver->getObjSense()*solver->getObjValue();
3714              //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3715              cutoff =model->getCutoff();
3716            }
3717            solver->markHotStart();
3718          }
3719#ifdef DO_ALL_AT_ROOT
3720          if (!numberNodes)
3721            printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3722                   choice.objectNumber,iStatus,newObjectiveValue,choice.numItersDown,
3723                   choice.downMovement,choice.finishedDown,choice.numIntInfeasDown,
3724                   choice.numObjInfeasDown);
3725#endif
3726         
3727          // repeat the whole exercise, forcing the variable up
3728          choice.possibleBranch->branch();
3729          solver->solveFromHotStart() ;
3730          numberStrongDone++;
3731          numberStrongIterations += solver->getIterationCount();
3732          /*
3733            We now have an estimate of objective degradation that we can use for strong
3734            branching. If we're over the cutoff, the variable is monotone up.
3735            If we actually made it to optimality, check for a solution, and if we have
3736            a good one, call setBestSolution to process it. Note that this may reduce the
3737            cutoff, so we check again to see if we can declare this variable monotone.
3738          */
3739          if (solver->isProvenOptimal())
3740            iStatus=0; // optimal
3741          else if (solver->isIterationLimitReached()
3742                   &&!solver->isDualObjectiveLimitReached())
3743            iStatus=2; // unknown
3744          else
3745            iStatus=1; // infeasible
3746          if (iStatus!=2&&solver->getIterationCount()>
3747              realMaxHotIterations)
3748            numberUnfinished++;
3749          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3750          choice.numItersUp = solver->getIterationCount();
3751          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3752          // Update branching information if wanted
3753          cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3754          if (cbcobj) {
3755            CbcObject * object = cbcobj->object();
3756            assert (object) ;
3757            CbcObjectUpdateData update = object->createUpdateInformation(solver,this,cbcobj);
3758            update.objectNumber_ = choice.objectNumber;
3759            model->addUpdateInformation(update);
3760          } else {
3761            decision->updateInformation( solver,this);
3762          }
3763          if (!iStatus) {
3764            choice.finishedUp = true ;
3765            if (newObjectiveValue>=cutoff) {
3766              objectiveChange = 1.0e100; // say infeasible
3767              numberStrongInfeasible++;
3768            } else {
3769#ifdef TIGHTEN_BOUNDS
3770              // Can we tighten bounds?
3771              if (iColumn<numberColumns&&cutoff<1.0e20
3772                  &&objectiveChange>1.0e-5) {
3773                double value = saveSolution[iColumn];
3774                double up = ceil(value)-value;
3775                double changePer = objectiveChange/(up+1.0e-7);
3776                double distance = (cutoff-objectiveValue_)/changePer;
3777                distance += 1.0e-3;
3778                if (distance<5.0) {
3779                  double newUpper = floor(value+distance);
3780                  if (newUpper<saveUpper[iColumn]) {
3781                    //printf("Could decrease upper bound on %d from %g to %g\n",
3782                    //   iColumn,saveUpper[iColumn],newUpper);
3783                    saveUpper[iColumn]=newUpper;
3784                    solver->setColUpper(iColumn,newUpper);
3785                  }
3786                }
3787                }
3788#endif
3789              // See if integer solution
3790              if (model->feasibleSolution(choice.numIntInfeasUp,
3791                                          choice.numObjInfeasUp)
3792                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
3793#ifdef BONMIN
3794                std::cout<<"Node has become integer feasible"<<std::endl;
3795                numberUnsatisfied_ = 0;
3796                break;
3797#endif
3798                if (auxiliaryInfo->solutionAddsCuts()) {
3799                  needHotStartUpdate=true;
3800                  solver->unmarkHotStart();
3801                }
3802                model->setBestSolution(CBC_STRONGSOL,
3803                                       newObjectiveValue,
3804                                       solver->getColSolution()) ;
3805                if (needHotStartUpdate) {
3806                  model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3807                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3808                  //objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
3809                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3810                  model->feasibleSolution(choice.numIntInfeasDown,
3811                                          choice.numObjInfeasDown);
3812                }
3813                model->setLastHeuristic(NULL);
3814                model->incrementUsed(solver->getColSolution());
3815                cutoff =model->getCutoff();
3816                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3817                  objectiveChange = 1.0e100 ;
3818                  numberStrongInfeasible++;
3819                }
3820              }
3821            }
3822          } else if (iStatus==1) {
3823            objectiveChange = 1.0e100 ;
3824            numberStrongInfeasible++;
3825          } else {
3826            // Can't say much as we did not finish
3827            choice.finishedUp = false ;
3828            numberUnfinished++;
3829          }
3830          choice.upMovement = objectiveChange ;
3831         
3832          // restore bounds
3833          for ( j=0;j<numberColumns;j++) {
3834            if (saveLower[j] != lower[j])
3835              solver->setColLower(j,saveLower[j]);
3836            if (saveUpper[j] != upper[j])
3837              solver->setColUpper(j,saveUpper[j]);
3838          }
3839          if(needHotStartUpdate) {
3840            needHotStartUpdate = false;
3841            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3842            //double newObjValue = solver->getObjSense()*solver->getObjValue();
3843            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3844            //we may again have an integer feasible solution
3845            int numberIntegerInfeasibilities;
3846            int numberObjectInfeasibilities;
3847            if (model->feasibleSolution(
3848                                        numberIntegerInfeasibilities,
3849                                        numberObjectInfeasibilities)) {
3850              double objValue = solver->getObjValue();
3851              model->setBestSolution(CBC_STRONGSOL,
3852                                     objValue,
3853                                     solver->getColSolution()) ;
3854              model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3855              //double newObjValue = solver->getObjSense()*solver->getObjValue();
3856              //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3857              cutoff =model->getCutoff();
3858            }
3859            solver->markHotStart();
3860          }
3861         
3862#ifdef DO_ALL_AT_ROOT
3863          if (!numberNodes)
3864            printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3865                   choice.objectNumber,iStatus,newObjectiveValue,choice.numItersUp,
3866                   choice.upMovement,choice.finishedUp,choice.numIntInfeasUp,
3867                   choice.numObjInfeasUp);
3868#endif
3869        }
3870       
3871        solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit2); 
3872        /*
3873          End of evaluation for this candidate variable. Possibilities are:
3874          * Both sides below cutoff; this variable is a candidate for branching.
3875          * Both sides infeasible or above the objective cutoff: no further action
3876          here. Break from the evaluation loop and assume the node will be purged
3877          by the caller.
3878          * One side below cutoff: Install the branch (i.e., fix the variable). Break
3879          from the evaluation loop and assume the node will be reoptimised by the
3880          caller.
3881        */
3882        // reset
3883        choice.possibleBranch->resetNumberBranchesLeft();
3884        if (choice.upMovement<1.0e100) {
3885          if(choice.downMovement<1.0e100) {
3886            // In case solution coming in was odd
3887            choice.upMovement = CoinMax(0.0,choice.upMovement);
3888            choice.downMovement = CoinMax(0.0,choice.downMovement);
3889#if ZERO_ONE==2
3890            // branch on 0-1 first (temp)
3891            if (fabs(choice.possibleBranch->value())<1.0) {
3892              choice.upMovement *= ZERO_FAKE;
3893              choice.downMovement *= ZERO_FAKE;
3894            }
3895#endif
3896            // feasible - see which best
3897            if (!canSkip) {
3898              if (model->messageHandler()->logLevel()>3) 
3899                printf("sort %g downest %g upest %g ",sort[iDo],downEstimate[iObject],
3900                       upEstimate[iObject]);
3901              model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
3902                << iObject << iColumn
3903                <<choice.downMovement<<choice.numIntInfeasDown
3904                <<choice.upMovement<<choice.numIntInfeasUp
3905                <<choice.possibleBranch->value()
3906                <<CoinMessageEol;
3907            }
3908            int betterWay;
3909            {
3910              CbcBranchingObject * branchObj =
3911                dynamic_cast <CbcBranchingObject *>(branch_) ;
3912              if (branch_)
3913                assert (branchObj);
3914              betterWay = decision->betterBranch(choice.possibleBranch,
3915                                                 branchObj,
3916                                                 choice.upMovement,
3917                                                 choice.numIntInfeasUp ,
3918                                                 choice.downMovement,
3919                                                 choice.numIntInfeasDown );
3920            }
3921            if (betterWay) {
3922              // C) create branching object
3923              if (choiceObject) {
3924                delete branch_;
3925                branch_ = choice.possibleBranch->clone();
3926              } else {
3927                delete branch_;
3928                branch_ = choice.possibleBranch;
3929                choice.possibleBranch=NULL;
3930              }
3931              {
3932                CbcBranchingObject * branchObj =
3933                  dynamic_cast <CbcBranchingObject *>(branch_) ;
3934                assert (branchObj);
3935                //branchObj->way(preferredWay);
3936                branchObj->way(betterWay);
3937              }
3938              bestChoice = choice.objectNumber;
3939              whichChoice = iDo;
3940              if (numberStrong<=1) {
3941                delete ws;
3942                ws=NULL;
3943                break;
3944              }
3945            } else {
3946              if (!choiceObject) {
3947                delete choice.possibleBranch;
3948                choice.possibleBranch=NULL;
3949              }
3950              if (iDo>=2*numberStrong) {
3951                delete ws;
3952                ws=NULL;
3953                break;
3954              }
3955              if (!dynamicObject||dynamicObject->numberTimesUp()>1) {
3956                if (iDo-whichChoice>=numberStrong) {
3957                  if (!choiceObject) {
3958                    delete choice.possibleBranch;
3959                    choice.possibleBranch=NULL;
3960                  }
3961                  break; // give up
3962                }
3963              } else {
3964                if (iDo-whichChoice>=2*numberStrong) {
3965                  delete ws;
3966                  ws=NULL;
3967                  if (!choiceObject) {
3968                    delete choice.possibleBranch;
3969                    choice.possibleBranch=NULL;
3970                  }
3971                  break; // give up
3972                }
3973              }
3974            }
3975          } else {
3976            // up feasible, down infeasible
3977            anyAction=-1;
3978            worstFeasible = CoinMax(worstFeasible,choice.upMovement);
3979            model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
3980              << iObject << iColumn
3981              <<choice.downMovement<<choice.numIntInfeasDown
3982              <<choice.upMovement<<choice.numIntInfeasUp
3983              <<choice.possibleBranch->value()
3984              <<CoinMessageEol;
3985            //printf("Down infeasible for choice %d sequence %d\n",i,
3986            // model->object(choice.objectNumber)->columnNumber());
3987            choice.fix=1;
3988            numberToFix++;
3989            choice.possibleBranch->fix(solver,saveLower,saveUpper,1);
3990            if (!choiceObject) {
3991              choice.possibleBranch=NULL;
3992            } else {
3993              //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3994              choice.possibleBranch=choiceObject;
3995            }
3996            assert(doneHotStart);
3997            solver->unmarkHotStart();
3998            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3999            //double newObjValue = solver->getObjSense()*solver->getObjValue();
4000            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
4001            solver->markHotStart();
4002            // may be infeasible (if other way stopped on iterations)
4003            if (!solver->isProvenOptimal()) {
4004              // neither side feasible
4005              anyAction=-2;
4006              if (!choiceObject) {
4007                delete choice.possibleBranch;
4008                choice.possibleBranch=NULL;
4009              }
4010              //printf("Both infeasible for choice %d sequence %d\n",i,
4011              // model->object(choice.objectNumber)->columnNumber());
4012              delete ws;
4013              ws=NULL;
4014              break;
4015            }
4016          }
4017        } else {
4018          if(choice.downMovement<1.0e100) {
4019            // down feasible, up infeasible
4020            anyAction=-1;
4021            worstFeasible = CoinMax(worstFeasible,choice.downMovement);
4022            model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
4023              << iObject << iColumn
4024              <<choice.downMovement<<choice.numIntInfeasDown
4025              <<choice.upMovement<<choice.numIntInfeasUp
4026              <<choice.possibleBranch->value()
4027              <<CoinMessageEol;
4028            choice.fix=-1;
4029            numberToFix++;
4030            choice.possibleBranch->fix(solver,saveLower,saveUpper,-1);
4031            if (!choiceObject) {
4032              choice.possibleBranch=NULL;
4033            } else {
4034              //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
4035              choice.possibleBranch=choiceObject;
4036            }
4037            assert(doneHotStart);
4038            solver->unmarkHotStart();
4039            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
4040            //double newObjValue = solver->getObjSense()*solver->getObjValue();
4041            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
4042            solver->markHotStart();
4043            // may be infeasible (if other way stopped on iterations)
4044            if (!solver->isProvenOptimal()) {
4045              // neither side feasible
4046              anyAction=-2;
4047              if (!choiceObject) {
4048                delete choice.possibleBranch;
4049                choice.possibleBranch=NULL;
4050              }
4051              delete ws;
4052              ws=NULL;
4053              break;
4054            }
4055          } else {
4056            // neither side feasible
4057            anyAction=-2;
4058            if (!choiceObject) {
4059              delete choice.possibleBranch;
4060              choice.possibleBranch=NULL;
4061            }
4062            delete ws;
4063            ws=NULL;
4064            break;
4065          }
4066        }
4067        // Check max time
4068        hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
4069                       model->getDblParam(CbcModel::CbcMaximumSeconds));
4070        if (hitMaxTime) {
4071          // make sure rest are fast
4072          doQuickly=true;
4073          for ( int jDo=iDo+1;jDo<numberToDo;jDo++) {
4074            int iObject = whichObject[iDo];
4075            OsiObject * object = model->modifiableObject(iObject);
4076            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4077              dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4078            if (dynamicObject) 
4079              dynamicObject->setNumberBeforeTrust(0);
4080          }
4081          numberTest=0;
4082        }
4083        if (!choiceObject) {
4084          delete choice.possibleBranch;
4085        }
4086      }
4087      if (model->messageHandler()->logLevel()>3) { 
4088        if (anyAction==-2) {
4089          printf("infeasible\n");
4090        } else if(anyAction==-1) {
4091          printf("%d fixed AND choosing %d iDo %d iChosenWhen %d numberToDo %d\n",numberToFix,bestChoice,
4092                   iDo,whichChoice,numberToDo);
4093        } else {
4094          int iObject = whichObject[whichChoice];
4095          OsiObject * object = model->modifiableObject(iObject);
4096          CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4097            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4098          if (dynamicObject) {
4099            int iColumn = dynamicObject->columnNumber();
4100            printf("choosing %d (column %d) iChosenWhen %d numberToDo %d\n",bestChoice,
4101                   iColumn,whichChoice,numberToDo);
4102          }
4103        }
4104      }
4105      if (doneHotStart) {
4106        // Delete the snapshot
4107        solver->unmarkHotStart();
4108        // back to normal
4109        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
4110        // restore basis
4111        solver->setWarmStart(ws);
4112      }
4113      solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
4114      // Unless infeasible we will carry on
4115      // But we could fix anyway
4116      if (numberToFix&&!hitMaxTime) {
4117        if (anyAction!=-2) {
4118          // apply and take off
4119          bool feasible=true;
4120          // can do quick optimality check
4121          int easy=2;
4122          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
4123          model->resolve(NULL,11,saveSolution,saveLower,saveUpper) ;
4124          //double newObjValue = solver->getObjSense()*solver->getObjValue();
4125          //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
4126          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
4127          feasible = solver->isProvenOptimal();
4128          if (feasible) {
4129            anyAction=0;
4130            // See if candidate still possible
4131            if (branch_) {
4132              const OsiObject * object = model->object(bestChoice);
4133              int preferredWay;
4134              double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
4135              if (!infeasibility) {
4136                // take out
4137                delete branch_;
4138                branch_=NULL;
4139              } else {
4140                CbcBranchingObject * branchObj =
4141                  dynamic_cast <CbcBranchingObject *>(branch_) ;
4142                assert (branchObj);
4143                branchObj->way(preferredWay);
4144              }
4145            }
4146          } else {
4147            anyAction=-2;
4148            finished=true;
4149          }
4150        }
4151      }
4152      // If  fixed then round again
4153      if (!branch_&&anyAction!=-2) {
4154        finished=false;
4155      }
4156      delete ws;
4157    }
4158  }
4159  // update number of strong iterations etc
4160  model->incrementStrongInfo(numberStrongDone,numberStrongIterations,
4161                             anyAction==-2 ? 0:numberToFix,anyAction==-2);
4162  if (model->searchStrategy()==-1) {
4163#ifndef COIN_DEVELOP
4164    if (solver->messageHandler()->logLevel()>1)
4165#endif
4166      printf("%d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
4167             numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
4168             numberNotTrusted);
4169    // decide what to do
4170    int strategy=1;
4171    if (((numberUnfinished*4>numberStrongDone&&
4172          numberStrongInfeasible*40<numberStrongDone)||
4173         numberStrongInfeasible<0)&&model->numberStrong()<10&&model->numberBeforeTrust()<=20&&model->numberObjects()>CoinMax(1000,solver->getNumRows())) {
4174      strategy=2;
4175#ifdef COIN_DEVELOP
4176      //if (model->logLevel()>1)
4177      printf("going to strategy 2\n");
4178#endif
4179      // Weaken
4180      model->setNumberStrong(2);
4181      model->setNumberBeforeTrust(1);
4182      model->synchronizeNumberBeforeTrust();
4183    }
4184    if (numberNodes)
4185      strategy=1;  // should only happen after hot start
4186    model->setSearchStrategy(strategy);
4187  } else if (numberStrongDone) {
4188    //printf("%d strongB, %d iters, %d inf, %d not finished, %d not trusted\n",
4189    //   numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
4190    //   numberNotTrusted);
4191  }
4192  if (model->searchStrategy()==1&&numberNodes>500&&numberNodes<-510) {
4193#ifndef COIN_DEVELOP
4194    if (solver->messageHandler()->logLevel()>1)
4195#endif
4196      printf("after %d nodes - %d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
4197             numberNodes,numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
4198             numberNotTrusted);
4199    // decide what to do
4200    if (numberUnfinished*10>numberStrongDone+1||
4201        !numberStrongInfeasible) {
4202      printf("going to strategy 2\n");
4203      // Weaken
4204      model->setNumberStrong(2);
4205      model->setNumberBeforeTrust(1);
4206      model->synchronizeNumberBeforeTrust();
4207      model->setSearchStrategy(2);
4208    }
4209  }
4210  // Set guessed solution value
4211  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
4212 
4213  /*
4214    Cleanup, then we're finished
4215  */
4216  if (!model->branchingMethod())
4217    delete decision;
4218 
4219  delete choiceObject;
4220  delete [] sort;
4221  delete [] whichObject;
4222#ifdef RANGING
4223  delete [] objectMark;
4224#endif
4225  delete [] saveLower;
4226  delete [] saveUpper;
4227  delete [] upEstimate;
4228  delete [] downEstimate;
4229# ifdef COIN_HAS_CLP
4230  if (osiclp) {
4231    osiclp->setSpecialOptions(saveClpOptions);
4232  }
4233# endif
4234  // restore solution
4235  solver->setColSolution(saveSolution);
4236  model->reserveCurrentSolution(saveSolution);
4237  delete [] saveSolution;
4238  model->setStateOfSearch(saveStateOfSearch);
4239  model->setLogLevel(saveLogLevel);
4240  // delete extra regions
4241  if (usefulInfo.usefulRegion_) {
4242    delete [] usefulInfo.usefulRegion_;
4243    delete [] usefulInfo.indexRegion_;
4244    delete [] usefulInfo.pi_;
4245    usefulInfo.usefulRegion_ = NULL;
4246    usefulInfo.indexRegion_ = NULL;
4247    usefulInfo.pi_ = NULL;
4248  }
4249  useShadow = model->moreSpecialOptions()&7;
4250  if ((useShadow==5&&model->getSolutionCount())||useShadow==6) {
4251    // zap pseudo shadow prices
4252    model->pseudoShadow(-1);
4253    // and switch off
4254    model->setMoreSpecialOptions(model->moreSpecialOptions()&(~1023));
4255  }
4256  return anyAction;
4257}
4258int CbcNode::analyze (CbcModel *model, double * results)
4259{
4260  int i;
4261  int numberIterationsAllowed = model->numberAnalyzeIterations();
4262  OsiSolverInterface * solver = model->solver();
4263  objectiveValue_ = solver->getObjSense()*solver->getObjValue();
4264  double cutoff =model->getCutoff();
4265  const double * lower = solver->getColLower();
4266  const double * upper = solver->getColUpper();
4267  const double * dj = solver->getReducedCost();
4268  int numberObjects = model->numberObjects();
4269  int numberColumns = model->getNumCols();
4270  // Initialize arrays
4271  int numberIntegers = model->numberIntegers();
4272  int * back = new int[numberColumns];
4273  const int * integerVariable = model->integerVariable();
4274  for (i=0;i<numberColumns;i++) 
4275    back[i]=-1;
4276  // What results is
4277  double * newLower = results;
4278  double * objLower = newLower+numberIntegers;
4279  double * newUpper = objLower+numberIntegers;
4280  double * objUpper = newUpper+numberIntegers;
4281  for (i=0;i<numberIntegers;i++) {
4282    int iColumn = integerVariable[i];
4283    back[iColumn]=i;
4284    newLower[i]=0.0;
4285    objLower[i]=-COIN_DBL_MAX;
4286    newUpper[i]=0.0;
4287    objUpper[i]=-COIN_DBL_MAX;
4288  }
4289  double * saveUpper = new double[numberColumns];
4290  double * saveLower = new double[numberColumns];
4291  int anyAction=0;
4292  // Save solution in case heuristics need good solution later
4293 
4294  double * saveSolution = new double[numberColumns];
4295  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
4296  model->reserveCurrentSolution(saveSolution);
4297  for (i=0;i<numberColumns;i++) {
4298    saveLower[i] = lower[i];
4299    saveUpper[i] = upper[i];
4300  }
4301  // Get arrays to sort
4302  double * sort = new double[numberObjects];
4303  int * whichObject = new int[numberObjects];
4304  int numberToFix=0;
4305  int numberToDo=0;
4306  double integerTolerance = 
4307    model->getDblParam(CbcModel::CbcIntegerTolerance);
4308  // point to useful information
4309  OsiBranchingInformation usefulInfo = model->usefulInformation();
4310  // and modify
4311  usefulInfo.depth_=depth_;
4312 
4313  // compute current state
4314  int numberObjectInfeasibilities; // just odd ones
4315  int numberIntegerInfeasibilities;
4316  model->feasibleSolution(
4317                          numberIntegerInfeasibilities,
4318                          numberObjectInfeasibilities);
4319# ifdef COIN_HAS_CLP
4320  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
4321  int saveClpOptions=0;
4322  bool fastIterations = (model->specialOptions()&8)!=0;
4323  if (osiclp&&fastIterations) {
4324    // for faster hot start
4325    saveClpOptions = osiclp->specialOptions();
4326    osiclp->setSpecialOptions(saveClpOptions|8192);
4327  }
4328# else
4329  bool fastIterations = false ;
4330# endif
4331  /*
4332    Scan for branching objects that indicate infeasibility. Choose candidates
4333    using priority as the first criteria, then integer infeasibility.
4334   
4335    The algorithm is to fill the array with a set of good candidates (by
4336    infeasibility) with priority bestPriority.  Finding a candidate with
4337    priority better (less) than bestPriority flushes the choice array. (This
4338    serves as initialization when the first candidate is found.)
4339   
4340  */
4341  numberToDo=0;
4342  for (i=0;i<numberObjects;i++) {
4343    OsiObject * object = model->modifiableObject(i);
4344    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4345      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4346    if(!dynamicObject)
4347      continue;
4348    int preferredWay;
4349    double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
4350    int iColumn = dynamicObject->columnNumber();
4351    if (saveUpper[iColumn]==saveLower[iColumn])
4352      continue;
4353    if (infeasibility)
4354      sort[numberToDo]=-1.0e10-infeasibility;
4355    else
4356      sort[numberToDo]=-fabs(dj[iColumn]);
4357    whichObject[numberToDo++]=i;
4358  }
4359  // Save basis
4360  CoinWarmStart * ws = solver->getWarmStart();
4361  int saveLimit;
4362  solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
4363  int targetIterations = CoinMax(500,numberIterationsAllowed/numberObjects);
4364  if (saveLimit<targetIterations)
4365    solver->setIntParam(OsiMaxNumIterationHotStart,targetIterations); 
4366  // Mark hot start
4367  solver->markHotStart();
4368  // Sort
4369  CoinSort_2(sort,sort+numberToDo,whichObject);
4370  double * currentSolution = model->currentSolution();
4371  double objMin = 1.0e50;
4372  double objMax = -1.0e50;
4373  bool needResolve=false;
4374  int iDo;
4375  for (iDo=0;iDo<numberToDo;iDo++) {
4376    CbcStrongInfo choice;
4377    int iObject = whichObject[iDo];
4378    OsiObject * object = model->modifiableObject(iObject);
4379    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4380      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4381    if (!dynamicObject)
4382      continue;
4383    int iColumn = dynamicObject->columnNumber();
4384    int preferredWay;
4385    object->infeasibility(&usefulInfo,preferredWay);
4386    double value = currentSolution[iColumn];
4387    double nearest = floor(value+0.5);
4388    double lowerValue = floor(value);
4389    bool satisfied=false;
4390    if (fabs(value-nearest)<=integerTolerance||value<saveLower[iColumn]||value>saveUpper[iColumn]) {
4391      satisfied=true;
4392      double newValue;
4393      if (nearest<saveUpper[iColumn]) {
4394        newValue = nearest + 1.0001*integerTolerance;
4395        lowerValue = nearest;
4396      } else {
4397        newValue = nearest - 1.0001*integerTolerance;
4398        lowerValue = nearest-1;
4399      }
4400      currentSolution[iColumn]=newValue;
4401    }
4402    double upperValue = lowerValue+1.0;
4403    //CbcSimpleInteger * obj =
4404    //dynamic_cast <CbcSimpleInteger *>(object) ;
4405    //if (obj) {
4406    //choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
4407    //} else {
4408      CbcObject * obj =
4409        dynamic_cast <CbcObject *>(object) ;
4410      assert (obj);
4411      choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
4412      //}
4413    currentSolution[iColumn]=value;
4414    // Save which object it was
4415    choice.objectNumber=iObject;
4416    choice.numIntInfeasUp = numberUnsatisfied_;
4417    choice.numIntInfeasDown = numberUnsatisfied_;
4418    choice.downMovement = 0.0;
4419    choice.upMovement = 0.0;
4420    choice.numItersDown = 0;
4421    choice.numItersUp = 0;
4422    choice.fix=0; // say not fixed
4423    double objectiveChange ;
4424    double newObjectiveValue=1.0e100;
4425    int j;
4426    // status is 0 finished, 1 infeasible and other
4427    int iStatus;
4428    /*
4429      Try the down direction first. (Specify the initial branching alternative as
4430      down with a call to way(-1). Each subsequent call to branch() performs the
4431      specified branch and advances the branch object state to the next branch
4432      alternative.)
4433    */
4434    choice.possibleBranch->way(-1) ;
4435    choice.possibleBranch->branch() ;
4436    if (fabs(value-lowerValue)>integerTolerance) {
4437      solver->solveFromHotStart() ;
4438      /*
4439        We now have an estimate of objective degradation that we can use for strong
4440        branching. If we're over the cutoff, the variable is monotone up.
4441        If we actually made it to optimality, check for a solution, and if we have
4442        a good one, call setBestSolution to process it. Note that this may reduce the
4443        cutoff, so we check again to see if we can declare this variable monotone.
4444      */
4445      if (solver->isProvenOptimal())
4446        iStatus=0; // optimal
4447      else if (solver->isIterationLimitReached()
4448               &&!solver->isDualObjectiveLimitReached())
4449        iStatus=2; // unknown
4450      else
4451        iStatus=1; // infeasible
4452      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
4453      choice.numItersDown = solver->getIterationCount();
4454      numberIterationsAllowed -= choice.numItersDown;
4455      objectiveChange = newObjectiveValue  - objectiveValue_;
4456      if (!iStatus) {
4457        choice.finishedDown = true ;
4458        if (newObjectiveValue>=cutoff) {
4459          objectiveChange = 1.0e100; // say infeasible
4460        } else {
4461          // See if integer solution
4462          if (model->feasibleSolution(choice.numIntInfeasDown,
4463                                      choice.numObjInfeasDown)
4464              &&model->problemFeasibility()->feasible(model,-1)>=0) {
4465            model->setBestSolution(CBC_STRONGSOL,
4466                                   newObjectiveValue,
4467                                   solver->getColSolution()) ;
4468            model->setLastHeuristic(NULL);
4469            model->incrementUsed(solver->getColSolution());
4470            cutoff =model->getCutoff();
4471            if (newObjectiveValue >= cutoff)    //  *new* cutoff
4472              objectiveChange = 1.0e100 ;
4473          }
4474        }
4475      } else if (iStatus==1) {
4476        objectiveChange = 1.0e100 ;
4477      } else {
4478        // Can't say much as we did not finish
4479        choice.finishedDown = false ;
4480      }
4481      choice.downMovement = objectiveChange ;
4482    }
4483    // restore bounds
4484    for ( j=0;j<numberColumns;j++) {
4485      if (saveLower[j] != lower[j])
4486        solver->setColLower(j,saveLower[j]);
4487      if (saveUpper[j] != upper[j])
4488        solver->setColUpper(j,saveUpper[j]);
4489    }
4490    // repeat the whole exercise, forcing the variable up
4491    choice.possibleBranch->branch();
4492    if (fabs(value-upperValue)>integerTolerance) {
4493      solver->solveFromHotStart() ;
4494      /*
4495        We now have an estimate of objective degradation that we can use for strong
4496        branching. If we're over the cutoff, the variable is monotone up.
4497        If we actually made it to optimality, check for a solution, and if we have
4498        a good one, call setBestSolution to process it. Note that this may reduce the
4499        cutoff, so we check again to see if we can declare this variable monotone.
4500      */
4501      if (solver->isProvenOptimal())
4502        iStatus=0; // optimal
4503      else if (solver->isIterationLimitReached()
4504               &&!solver->isDualObjectiveLimitReached())
4505        iStatus=2; // unknown
4506      else
4507        iStatus=1; // infeasible
4508      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
4509      choice.numItersUp = solver->getIterationCount();
4510      numberIterationsAllowed -= choice.numItersUp;
4511      objectiveChange = newObjectiveValue  - objectiveValue_;
4512      if (!iStatus) {
4513        choice.finishedUp = true ;
4514        if (newObjectiveValue>=cutoff) {
4515          objectiveChange = 1.0e100; // say infeasible
4516        } else {
4517          // See if integer solution
4518          if (model->feasibleSolution(choice.numIntInfeasUp,
4519                                      choice.numObjInfeasUp)
4520              &&model->problemFeasibility()->feasible(model,-1)>=0) {
4521            model->setBestSolution(CBC_STRONGSOL,
4522                                   newObjectiveValue,
4523                                   solver->getColSolution()) ;
4524            model->setLastHeuristic(NULL);
4525            model->incrementUsed(solver->getColSolution());
4526            cutoff =model->getCutoff();
4527            if (newObjectiveValue >= cutoff)    //  *new* cutoff
4528              objectiveChange = 1.0e100 ;
4529          }
4530        }
4531      } else if (iStatus==1) {
4532        objectiveChange = 1.0e100 ;
4533      } else {
4534        // Can't say much as we did not finish
4535        choice.finishedUp = false ;
4536      }
4537      choice.upMovement = objectiveChange ;
4538     
4539      // restore bounds
4540      for ( j=0;j<numberColumns;j++) {
4541        if (saveLower[j] != lower[j])
4542          solver->setColLower(j,saveLower[j]);
4543        if (saveUpper[j] != upper[j])
4544          solver->setColUpper(j,saveUpper[j]);
4545      }
4546    }
4547    // If objective goes above certain amount we can set bound
4548    int jInt = back[iColumn];
4549    newLower[jInt]=upperValue;
4550    if (choice.finishedDown)
4551      objLower[jInt]=choice.downMovement+objectiveValue_;
4552    else
4553      objLower[jInt]=objectiveValue_;
4554    newUpper[jInt]=lowerValue;
4555    if (choice.finishedUp)
4556      objUpper[jInt]=choice.upMovement+objectiveValue_;
4557    else
4558      objUpper[jInt]=objectiveValue_;
4559    objMin = CoinMin(CoinMin(objLower[jInt],objUpper[jInt]),objMin);
4560    /*
4561      End of evaluation for this candidate variable. Possibilities are:
4562      * Both sides below cutoff; this variable is a candidate for branching.
4563      * Both sides infeasible or above the objective cutoff: no further action
4564      here. Break from the evaluation loop and assume the node will be purged
4565      by the caller.
4566      * One side below cutoff: Install the branch (i.e., fix the variable). Break
4567      from the evaluation loop and assume the node will be reoptimised by the
4568      caller.
4569    */
4570    if (choice.upMovement<1.0e100) {
4571      if(choice.downMovement<1.0e100) {
4572        objMax = CoinMax(CoinMax(objLower[jInt],objUpper[jInt]),objMax);
4573        // In case solution coming in was odd
4574        choice.upMovement = CoinMax(0.0,choice.upMovement);
4575        choice.downMovement = CoinMax(0.0,choice.downMovement);
4576        // feasible -
4577        model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
4578          << iObject << iColumn
4579          <<choice.downMovement<<choice.numIntInfeasDown
4580          <<choice.upMovement<<choice.numIntInfeasUp
4581          <<value
4582          <<CoinMessageEol;
4583      } else {
4584        // up feasible, down infeasible
4585        anyAction=-1;
4586        if (!satisfied)
4587          needResolve=true;
4588        choice.fix=1;
4589        numberToFix++;
4590        saveLower[iColumn]=upperValue;
4591        solver->setColLower(iColumn,upperValue);
4592      }
4593    } else {
4594      if(choice.downMovement<1.0e100) {
4595        // down feasible, up infeasible
4596        anyAction=-1;
4597        if (!satisfied)
4598          needResolve=true;
4599        choice.fix=-1;
4600        numberToFix++;
4601        saveUpper[iColumn]=lowerValue;
4602        solver->setColUpper(iColumn,lowerValue);
4603      } else {
4604        // neither side feasible
4605        anyAction=-2;
4606        printf("Both infeasible for choice %d sequence %d\n",i,
4607               model->object(choice.objectNumber)->columnNumber());
4608        delete ws;
4609        ws=NULL;
4610        //solver->writeMps("bad");
4611        numberToFix=-1;
4612        delete choice.possibleBranch;
4613        choice.possibleBranch=NULL;
4614        break;
4615      }
4616    }
4617    delete choice.possibleBranch;
4618    if (numberIterationsAllowed<=0)
4619      break;
4620    //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
4621    //     choice.downMovement,choice.upMovement,value);
4622  }
4623  printf("Best possible solution %g, can fix more if solution of %g found - looked at %d variables in %d iterations\n",
4624         objMin,objMax,iDo,model->numberAnalyzeIterations()-numberIterationsAllowed);
4625  model->setNumberAnalyzeIterations(numberIterationsAllowed);
4626  // Delete the snapshot
4627  solver->unmarkHotStart();
4628  // back to normal
4629  solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
4630  solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
4631  // restore basis
4632  solver->setWarmStart(ws);
4633  delete ws;
4634 
4635  delete [] sort;
4636  delete [] whichObject;
4637  delete [] saveLower;
4638  delete [] saveUpper;
4639  delete [] back;
4640  // restore solution
4641  solver->setColSolution(saveSolution);
4642# ifdef COIN_HAS_CLP
4643  if (osiclp) 
4644    osiclp->setSpecialOptions(saveClpOptions);
4645# endif
4646  model->reserveCurrentSolution(saveSolution);
4647  delete [] saveSolution;
4648  if (needResolve)
4649    solver->resolve();
4650  return numberToFix;
4651}
4652
4653
4654CbcNode::CbcNode(const CbcNode & rhs)
4655  : CoinTreeNode(rhs) 
4656{ 
4657#ifdef CHECK_NODE
4658  printf("CbcNode %x Constructor from rhs %x\n",this,&rhs);
4659#endif
4660  if (rhs.nodeInfo_)
4661    nodeInfo_ = rhs.nodeInfo_->clone();
4662  else
4663    nodeInfo_=NULL;
4664  objectiveValue_=rhs.objectiveValue_;
4665  guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4666  sumInfeasibilities_ = rhs.sumInfeasibilities_;
4667  if (rhs.branch_)
4668    branch_=rhs.branch_->clone();
4669  else
4670    branch_=NULL;
4671  depth_ = rhs.depth_;
4672  numberUnsatisfied_ = rhs.numberUnsatisfied_;
4673  nodeNumber_ = rhs.nodeNumber_;
4674  state_ = rhs.state_;
4675  if (nodeInfo_)
4676    assert ((state_&2)!=0);
4677  else
4678    assert ((state_&2)==0);
4679}
4680
4681CbcNode &
4682CbcNode::operator=(const CbcNode & rhs)
4683{
4684  if (this != &rhs) {
4685    delete nodeInfo_;
4686    if (rhs.nodeInfo_)
4687      nodeInfo_ = rhs.nodeInfo_->clone();
4688    else
4689      nodeInfo_ = NULL;
4690    objectiveValue_=rhs.objectiveValue_;
4691    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4692    sumInfeasibilities_ = rhs.sumInfeasibilities_;
4693    if (rhs.branch_)
4694      branch_=rhs.branch_->clone();
4695    else
4696      branch_=NULL,
4697        depth_ = rhs.depth_;
4698    numberUnsatisfied_ = rhs.numberUnsatisfied_;
4699    nodeNumber_ = rhs.nodeNumber_;
4700    state_ = rhs.state_;
4701    if (nodeInfo_)
4702      assert ((state_&2)!=0);
4703    else
4704      assert ((state_&2)==0);
4705  }
4706  return *this;
4707}
4708CbcNode::~CbcNode ()
4709{
4710#ifdef CHECK_NODE
4711  if (nodeInfo_) {
4712    printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
4713           this,nodeInfo_,nodeInfo_->numberPointingToThis());
4714    //assert(nodeInfo_->numberPointingToThis()>=0);
4715  } else {
4716    printf("CbcNode %x Destructor nodeInfo %x (?)\n",
4717           this,nodeInfo_);
4718  }
4719#endif
4720  if (nodeInfo_) {
4721    // was if (nodeInfo_&&(state_&2)!=0) {
4722    nodeInfo_->nullOwner();
4723    int numberToDelete=nodeInfo_->numberBranchesLeft();
4724    //    CbcNodeInfo * parent = nodeInfo_->parent();
4725    //assert (nodeInfo_->numberPointingToThis()>0);
4726    if (nodeInfo_->decrement(numberToDelete)==0||(state_&2)==0) {
4727      if ((state_&2)==0) 
4728        nodeInfo_->nullParent();
4729      delete nodeInfo_;
4730    } else {
4731      //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,nodeInfo_->parent());
4732      // anyway decrement parent
4733      //if (parent)
4734      ///parent->decrement(1);
4735    }
4736  }
4737  delete branch_;
4738}
4739// Decrement  active cut counts
4740void 
4741CbcNode::decrementCuts(int change)
4742{
4743  if (nodeInfo_)
4744    assert ((state_&2)!=0);
4745  else
4746    assert ((state_&2)==0);
4747  if(nodeInfo_) {
4748    nodeInfo_->decrementCuts(change);
4749  }
4750}
4751void 
4752CbcNode::decrementParentCuts(CbcModel * model, int change)
4753{
4754  if (nodeInfo_)
4755    assert ((state_&2)!=0);
4756  else
4757    assert ((state_&2)==0);
4758  if(nodeInfo_) {
4759    nodeInfo_->decrementParentCuts(model, change);
4760  }
4761}
4762
4763/*
4764  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
4765  in the attached nodeInfo_.
4766*/
4767void
4768CbcNode::initializeInfo()
4769{
4770  assert(nodeInfo_ && branch_) ;
4771  nodeInfo_->initializeInfo(branch_->numberBranches());
4772  assert ((state_&2)!=0);
4773  assert (nodeInfo_->numberBranchesLeft()==
4774          branch_->numberBranchesLeft());
4775}
4776// Nulls out node info
4777void 
4778CbcNode::nullNodeInfo()
4779{
4780  nodeInfo_=NULL;
4781  // say not active
4782  state_ &= ~2;
4783}
4784
4785int
4786CbcNode::branch(OsiSolverInterface * solver)
4787{
4788  double changeInGuessed;
4789  assert (nodeInfo_->numberBranchesLeft()==
4790          branch_->numberBranchesLeft());
4791  if (!solver)
4792    changeInGuessed=branch_->branch();
4793  else
4794    changeInGuessed=branch_->branch(solver);
4795  guessedObjectiveValue_+= changeInGuessed;
4796  //#define PRINTIT
4797#ifdef PRINTIT
4798  int numberLeft = nodeInfo_->numberBranchesLeft();
4799  CbcNodeInfo * parent = nodeInfo_->parent();
4800  int parentNodeNumber = -1;
4801  CbcBranchingObject * object1 = 
4802    dynamic_cast<CbcBranchingObject *>(branch_) ;
4803  //OsiObject * object = object1->
4804  //int sequence = object->columnNumber);
4805  int id=-1;
4806  double value=0.0;
4807  if (object1) {
4808    id = object1->variable();
4809    value = object1->value();
4810  }
4811  printf("id %d value %g objvalue %g\n",id,value,objectiveValue_);
4812  if (parent)
4813    parentNodeNumber = parent->nodeNumber();
4814  printf("Node number %d, %s, way %d, depth %d, parent node number %d\n",
4815         nodeInfo_->nodeNumber(),(numberLeft==2) ? "leftBranch" : "rightBranch",
4816         way(),depth_,parentNodeNumber);
4817  assert (parentNodeNumber!=nodeInfo_->nodeNumber());
4818#endif
4819  return nodeInfo_->branchedOn();
4820}
4821/* Active arm of the attached OsiBranchingObject.
4822   
4823   In the simplest instance, coded -1 for the down arm of the branch, +1 for
4824   the up arm. But see OsiBranchingObject::way()
4825   Use nodeInfo--.numberBranchesLeft_ to see how active
4826*/
4827int 
4828CbcNode::way() const
4829{
4830  if (branch_) {
4831    CbcBranchingObject * obj =
4832      dynamic_cast <CbcBranchingObject *>(branch_) ;
4833    if (obj) {
4834      return obj->way();
4835    } else {
4836      OsiTwoWayBranchingObject * obj2 =
4837        dynamic_cast <OsiTwoWayBranchingObject *>(branch_) ;
4838      assert (obj2);
4839      return obj2->way();
4840    }
4841  } else {
4842    return 0;
4843  }
4844}
4845/* Create a branching object for the node
4846   
4847   The routine scans the object list of the model and selects a set of
4848   unsatisfied objects as candidates for branching. The candidates are
4849   evaluated, and an appropriate branch object is installed.
4850   
4851   The numberPassesLeft is decremented to stop fixing one variable each time
4852   and going on and on (e.g. for stock cutting, air crew scheduling)
4853   
4854   If evaluation determines that an object is monotone or infeasible,
4855   the routine returns immediately. In the case of a monotone object,
4856   the branch object has already been called to modify the model.
4857   
4858   Return value:
4859   <ul>
4860   <li>  0: A branching object has been installed
4861   <li> -1: A monotone object was discovered
4862   <li> -2: An infeasible object was discovered
4863   </ul>
4864   Branch state:
4865   <ul>
4866   <li> -1: start
4867   <li> -1: A monotone object was discovered
4868   <li> -2: An infeasible object was discovered
4869   </ul>
4870*/
4871int 
4872CbcNode::chooseOsiBranch (CbcModel * model,
4873                          CbcNode * lastNode,
4874                          OsiBranchingInformation * usefulInfo,
4875                          int branchState)
4876{
4877  int returnStatus=0;
4878  if (lastNode)
4879    depth_ = lastNode->depth_+1;
4880  else
4881    depth_ = 0;
4882  OsiSolverInterface * solver = model->solver();
4883  objectiveValue_ = solver->getObjValue()*solver->getObjSense();
4884  usefulInfo->objectiveValue_ = objectiveValue_;
4885  usefulInfo->depth_ = depth_;
4886  const double * saveInfoSol = usefulInfo->solution_;
4887  double * saveSolution = new double[solver->getNumCols()];
4888  memcpy(saveSolution,solver->getColSolution(),solver->getNumCols()*sizeof(double));
4889  usefulInfo->solution_ = saveSolution;
4890  OsiChooseVariable * choose = model->branchingMethod()->chooseMethod();
4891  int numberUnsatisfied=-1;
4892  if (branchState<0) {
4893    // initialize
4894    // initialize sum of "infeasibilities"
4895    sumInfeasibilities_ = 0.0;
4896    numberUnsatisfied = choose->setupList(usefulInfo,true);
4897    numberUnsatisfied_ = numberUnsatisfied;
4898    branchState=0;
4899    if (numberUnsatisfied_<0) {
4900      // infeasible
4901      delete [] saveSolution;
4902      return -2;
4903    }
4904  }
4905  // unset best
4906  int best=-1;
4907  choose->setBestObjectIndex(-1);
4908  if (numberUnsatisfied) {
4909    if (branchState>0||!choose->numberOnList()) {
4910      // we need to return at once - don't do strong branching or anything
4911      if (choose->numberOnList()||!choose->numberStrong()) {
4912        best = choose->candidates()[0];
4913        choose->setBestObjectIndex(best);
4914      } else {
4915        // nothing on list - need to try again - keep any solution
4916        numberUnsatisfied = choose->setupList(usefulInfo, false);
4917        numberUnsatisfied_ = numberUnsatisfied;
4918        if (numberUnsatisfied) {
4919          best = choose->candidates()[0];
4920          choose->setBestObjectIndex(best);
4921        }
4922      }
4923    } else {
4924      // carry on with strong branching or whatever
4925      int returnCode = choose->chooseVariable(solver, usefulInfo,true);
4926      // update number of strong iterations etc
4927      model->incrementStrongInfo(choose->numberStrongDone(),choose->numberStrongIterations(),
4928                                 returnCode==-1 ? 0:choose->numberStrongFixed(),returnCode==-1);
4929      if (returnCode>1) {
4930        // has fixed some
4931        returnStatus=-1;
4932      } else if (returnCode==-1) {
4933        // infeasible
4934        returnStatus=-2;
4935      } else if (returnCode==0) {
4936        // normal
4937        returnStatus=0;
4938        numberUnsatisfied=1;
4939      } else {
4940        // ones on list satisfied - double check
4941        numberUnsatisfied = choose->setupList(usefulInfo, false);
4942        numberUnsatisfied_ = numberUnsatisfied;
4943        if (numberUnsatisfied) {
4944          best = choose->candidates()[0];
4945          choose->setBestObjectIndex(best);
4946        }
4947      }
4948    }
4949  } 
4950  delete branch_;
4951  branch_ = NULL;
4952  guessedObjectiveValue_ = COIN_DBL_MAX;//objectiveValue_; // for now
4953  if (!returnStatus) {
4954    if (numberUnsatisfied) {
4955      // create branching object
4956      const OsiObject * obj = model->solver()->object(choose->bestObjectIndex());
4957      //const OsiSolverInterface * solver = usefulInfo->solver_;
4958      branch_ = obj->createBranch(model->solver(),usefulInfo,obj->whichWay());
4959    }
4960  }
4961  usefulInfo->solution_=saveInfoSol;
4962  delete [] saveSolution;
4963  // may have got solution
4964  if (choose->goodSolution()
4965      &&model->problemFeasibility()->feasible(model,-1)>=0) {
4966    // yes
4967    double objValue = choose->goodObjectiveValue();
4968    model->setBestSolution(CBC_STRONGSOL,
4969                           objValue,
4970                           choose->goodSolution()) ;
4971    model->setLastHeuristic(NULL);
4972    model->incrementUsed(choose->goodSolution());
4973    choose->clearGoodSolution();
4974  }
4975  return returnStatus;
4976}
4977int 
4978CbcNode::chooseClpBranch (CbcModel * model,
4979                          CbcNode * lastNode)
4980{ 
4981  assert(lastNode);
4982  depth_ = lastNode->depth_+1;
4983  delete branch_;
4984  branch_=NULL;
4985  OsiSolverInterface * solver = model->solver();
4986  //double saveObjectiveValue = solver->getObjValue();
4987  //double objectiveValue = CoinMax(solver->getObjSense()*saveObjectiveValue,objectiveValue_);
4988  const double * lower = solver->getColLower();
4989  const double * upper = solver->getColUpper();
4990  // point to useful information
4991  OsiBranchingInformation usefulInfo = model->usefulInformation();
4992  // and modify
4993  usefulInfo.depth_=depth_;
4994  int i;
4995  //bool beforeSolution = model->getSolutionCount()==0;
4996  int numberObjects = model->numberObjects();
4997  int numberColumns = model->getNumCols();
4998  double * saveUpper = new double[numberColumns];
4999  double * saveLower = new double[numberColumns];
5000 
5001  // Save solution in case heuristics need good solution later
5002 
5003  double * saveSolution = new double[numberColumns];
5004  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
5005  model->reserveCurrentSolution(saveSolution);
5006  for (i=0;i<numberColumns;i++) {
5007    saveLower[i] = lower[i];
5008    saveUpper[i] = upper[i];
5009  }
5010  // Save basis
5011  CoinWarmStart * ws = solver->getWarmStart();
5012  numberUnsatisfied_ = 0;
5013  // initialize sum of "infeasibilities"
5014  sumInfeasibilities_ = 0.0;
5015  // Note looks as if off end (hidden one)
5016  OsiObject * object = model->modifiableObject(numberObjects);
5017  CbcGeneralDepth * thisOne = dynamic_cast <CbcGeneralDepth *> (object);
5018  assert (thisOne);
5019  OsiClpSolverInterface * clpSolver
5020    = dynamic_cast<OsiClpSolverInterface *> (solver);
5021  assert (clpSolver);
5022  ClpSimplex * simplex = clpSolver->getModelPtr();
5023  int preferredWay;
5024  double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
5025  if (thisOne->whichSolution()>=0) {
5026    ClpNode * nodeInfo = thisOne->nodeInfo(thisOne->whichSolution());
5027    nodeInfo->applyNode(simplex,2);
5028    int saveLogLevel = simplex->logLevel();
5029    simplex->setLogLevel(0);
5030    simplex->dual();
5031    simplex->setLogLevel(saveLogLevel);
5032    double cutoff = model->getCutoff();
5033    bool goodSolution=true;
5034    if (simplex->status()) {
5035      //simplex->writeMps("bad7.mps",2);
5036      if (nodeInfo->objectiveValue()>cutoff-1.0e-2)
5037        goodSolution=false;
5038      else
5039        assert (!simplex->status());
5040    }
5041    if (goodSolution) {
5042      double newObjectiveValue = solver->getObjSense()*solver->getObjValue();
5043      // See if integer solution
5044      int numInf;
5045      int numInf2;
5046      bool gotSol = model->feasibleSolution(numInf,numInf2);
5047      if (!gotSol) {
5048        printf("numinf %d\n",numInf);
5049        double * sol = simplex->primalColumnSolution();
5050        for (int i=0;i<numberColumns;i++) {
5051          if (simplex->isInteger(i)) {
5052            double value = floor(sol[i]+0.5);
5053            if (fabs(value-sol[i])>1.0e-7) {
5054              printf("%d value %g\n",i,sol[i]);
5055              if (fabs(value-sol[i])<1.0e-3) {
5056                sol[i]=value;
5057              }
5058            }
5059          }
5060        }
5061        simplex->writeMps("bad8.mps",2);
5062        bool gotSol = model->feasibleSolution(numInf,numInf2);
5063        if (!gotSol) 
5064          assert (gotSol);
5065      }
5066      model->setBestSolution(CBC_STRONGSOL,
5067                             newObjectiveValue,
5068                             solver->getColSolution()) ;
5069      model->setLastHeuristic(NULL);
5070      model->incrementUsed(solver->getColSolution());
5071    }
5072  }
5073  // restore bounds
5074  { for (int j=0;j<numberColumns;j++) {
5075      if (saveLower[j] != lower[j])
5076        solver->setColLower(j,saveLower[j]);
5077      if (saveUpper[j] != upper[j])
5078        solver->setColUpper(j,saveUpper[j]);
5079    }
5080  }
5081  // restore basis
5082  solver->setWarmStart(ws);
5083  delete ws;
5084  int anyAction;
5085  //#define CHECK_PATH
5086#ifdef CHECK_PATH
5087  extern int gotGoodNode_Z;
5088  if (gotGoodNode_Z>=0)
5089    printf("good node %d %g\n",gotGoodNode_Z,infeasibility);
5090#endif
5091  if (infeasibility>0.0) {
5092    if (infeasibility==COIN_DBL_MAX) {
5093      anyAction=-2; // infeasible
5094    } else {
5095      branch_=thisOne->createCbcBranch(solver,&usefulInfo,preferredWay);
5096      // Set to firat one (and change when re-pushing)
5097      CbcGeneralBranchingObject * branch = 
5098        dynamic_cast <CbcGeneralBranchingObject *> (branch_);
5099      branch->state(objectiveValue_,sumInfeasibilities_,
5100                    numberUnsatisfied_,0);
5101      branch->setNode(this);
5102      anyAction=0;
5103    }
5104  } else {
5105    anyAction=-1;
5106  }
5107#ifdef CHECK_PATH
5108  gotGoodNode_Z=-1;
5109#endif
5110  // Set guessed solution value
5111  guessedObjectiveValue_ = objectiveValue_+1.0e-5;
5112  delete [] saveLower;
5113  delete [] saveUpper;
5114 
5115  // restore solution
5116  solver->setColSolution(saveSolution);
5117  delete [] saveSolution;
5118  return anyAction;
5119}
5120/* Double checks in case node can change its mind!
5121   Returns objective value
5122   Can change objective etc */
5123double
5124CbcNode::checkIsCutoff(double cutoff)
5125{
5126  branch_->checkIsCutoff(cutoff);
5127  return objectiveValue_;
5128}
Note: See TracBrowser for help on using the repository browser.