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

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

fixes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 161.2 KB
Line 
1/* $Id: CbcNode.cpp 1212 2009-08-21 16:19:13Z 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  if (model->hotstartSolution())
1409    numberStrong=0;
1410  int numberStrongDone=0;
1411  int numberUnfinished=0;
1412  int numberStrongInfeasible=0;
1413  int numberStrongIterations=0;
1414  int saveNumberStrong=numberStrong;
1415  int numberObjects = model->numberObjects();
1416  bool checkFeasibility = numberObjects>model->numberIntegers();
1417  int maximumStrong = CoinMax(CoinMin(numberStrong,numberObjects),1);
1418  int numberColumns = model->getNumCols();
1419  double * saveUpper = new double[numberColumns];
1420  double * saveLower = new double[numberColumns];
1421 
1422  // Save solution in case heuristics need good solution later
1423 
1424  double * saveSolution = new double[numberColumns];
1425  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1426  model->reserveCurrentSolution(saveSolution);
1427  /*
1428    Get a branching decision object. Use the default decision criteria unless
1429    the user has loaded a decision method into the model.
1430  */
1431  CbcBranchDecision *decision = model->branchingMethod();
1432  CbcDynamicPseudoCostBranchingObject * dynamicBranchingObject =
1433    dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(decision);
1434  if (!decision||dynamicBranchingObject)
1435    decision = new CbcBranchDefaultDecision();
1436  decision->initialize(model);
1437  CbcStrongInfo * choice = new CbcStrongInfo[maximumStrong];
1438  for (i=0;i<numberColumns;i++) {
1439    saveLower[i] = lower[i];
1440    saveUpper[i] = upper[i];
1441  }
1442  // May go round twice if strong branching fixes all local candidates
1443  bool finished=false;
1444  double estimatedDegradation=0.0; 
1445  while(!finished) {
1446    finished=true;
1447    // Some objects may compute an estimate of best solution from here
1448    estimatedDegradation=0.0; 
1449    //int numberIntegerInfeasibilities=0; // without odd ones
1450    numberStrongDone=0;
1451    numberUnfinished=0;
1452    numberStrongInfeasible=0;
1453    numberStrongIterations=0;
1454   
1455    // We may go round this loop twice (only if we think we have solution)
1456    for (int iPass=0;iPass<2;iPass++) {
1457     
1458      // compute current state
1459      //int numberObjectInfeasibilities; // just odd ones
1460      //model->feasibleSolution(
1461      //                      numberIntegerInfeasibilities,
1462      //                      numberObjectInfeasibilities);
1463      const double * hotstartSolution = model->hotstartSolution();
1464      const int * hotstartPriorities = model->hotstartPriorities();
1465     
1466      // Some objects may compute an estimate of best solution from here
1467      estimatedDegradation=0.0; 
1468      numberUnsatisfied_ = 0;
1469      // initialize sum of "infeasibilities"
1470      sumInfeasibilities_ = 0.0;
1471      int bestPriority=COIN_INT_MAX;
1472      /*
1473        Scan for branching objects that indicate infeasibility. Choose the best
1474        maximumStrong candidates, using priority as the first criteria, then
1475        integer infeasibility.
1476       
1477        The algorithm is to fill the choice array with a set of good candidates (by
1478        infeasibility) with priority bestPriority.  Finding a candidate with
1479        priority better (less) than bestPriority flushes the choice array. (This
1480        serves as initialization when the first candidate is found.)
1481       
1482        A new candidate is added to choices only if its infeasibility exceeds the
1483        current max infeasibility (mostAway). When a candidate is added, it
1484        replaces the candidate with the smallest infeasibility (tracked by
1485        iSmallest).
1486      */
1487      int iSmallest = 0;
1488      double mostAway = 1.0e-100;
1489      for (i = 0 ; i < maximumStrong ; i++)
1490        choice[i].possibleBranch = NULL ;
1491      numberStrong=0;
1492      bool canDoOneHot=false;
1493      for (i=0;i<numberObjects;i++) {
1494        OsiObject * object = model->modifiableObject(i);
1495        int preferredWay;
1496        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
1497        int priorityLevel = object->priority();
1498        if (hotstartSolution) {
1499          // we are doing hot start
1500          const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
1501          if (thisOne) {
1502            int iColumn = thisOne->columnNumber();
1503            bool canDoThisHot=true;
1504            double targetValue = hotstartSolution[iColumn];
1505            if (saveUpper[iColumn]>saveLower[iColumn]) {
1506              double value = saveSolution[iColumn];
1507              if (hotstartPriorities)
1508                priorityLevel=hotstartPriorities[iColumn]; 
1509              //double originalLower = thisOne->originalLower();
1510              //double originalUpper = thisOne->originalUpper();
1511              // switch off if not possible
1512              if (targetValue>=saveLower[iColumn]&&targetValue<=saveUpper[iColumn]) {
1513                /* priority outranks rest always if negative
1514                   otherwise can be downgraded if at correct level.
1515                   Infeasibility may be increased to choose 1.0 values first.
1516                   choose one near wanted value
1517                */
1518                if (fabs(value-targetValue)>integerTolerance) {
1519                  infeasibility = fabs(value-targetValue);
1520                  //if (targetValue==1.0)
1521                  //infeasibility += 1.0;
1522                  if (value>targetValue) {
1523                    preferredWay=-1;
1524                  } else {
1525                    preferredWay=1;
1526                  }
1527                  priorityLevel = CoinAbs(priorityLevel);
1528                } else if (priorityLevel<0) {
1529                  priorityLevel = CoinAbs(priorityLevel);
1530                  if (targetValue==saveLower[iColumn]) {
1531                    infeasibility = integerTolerance+1.0e-12;
1532                    preferredWay=-1;
1533                  } else if (targetValue==saveUpper[iColumn]) {
1534                    infeasibility = integerTolerance+1.0e-12;
1535                    preferredWay=1;
1536                  } else {
1537                    // can't
1538                    priorityLevel += 10000000;
1539                    canDoThisHot=false;
1540                  }
1541                } else {
1542                  priorityLevel += 10000000;
1543                  canDoThisHot=false;
1544                }
1545              } else {
1546                // switch off if not possible
1547                canDoThisHot=false;
1548              }
1549              if (canDoThisHot)
1550                canDoOneHot=true;
1551            } else if (targetValue<saveLower[iColumn]||targetValue>saveUpper[iColumn]) {
1552            }
1553          } else {
1554            priorityLevel += 10000000;
1555          }
1556        }
1557        if (infeasibility) {
1558          // Increase estimated degradation to solution
1559          estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate());
1560          numberUnsatisfied_++;
1561          sumInfeasibilities_ += infeasibility;
1562          // Better priority? Flush choices.
1563          if (priorityLevel<bestPriority) {
1564            int j;
1565            iSmallest=0;
1566            for (j=0;j<maximumStrong;j++) {
1567              choice[j].upMovement=0.0;
1568              delete choice[j].possibleBranch;
1569              choice[j].possibleBranch=NULL;
1570            }
1571            bestPriority = priorityLevel;
1572            mostAway=1.0e-100;
1573            numberStrong=0;
1574          } else if (priorityLevel>bestPriority) {
1575            continue;
1576          }
1577          // Check for suitability based on infeasibility.
1578          if (infeasibility>mostAway) {
1579            //add to list
1580            choice[iSmallest].upMovement=infeasibility;
1581            delete choice[iSmallest].possibleBranch;
1582            CbcObject * obj =
1583              dynamic_cast <CbcObject *>(object) ;
1584            assert (obj);
1585            choice[iSmallest].possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
1586            numberStrong = CoinMax(numberStrong,iSmallest+1);
1587            // Save which object it was
1588            choice[iSmallest].objectNumber=i;
1589            int j;
1590            iSmallest=-1;
1591            mostAway = 1.0e50;
1592            for (j=0;j<maximumStrong;j++) {
1593              if (choice[j].upMovement<mostAway) {
1594                mostAway=choice[j].upMovement;
1595                iSmallest=j;
1596              }
1597            }
1598          }
1599        }
1600      }
1601      if (!canDoOneHot&&hotstartSolution) {
1602        // switch off as not possible
1603        hotstartSolution=NULL;
1604        model->setHotstartSolution(NULL,NULL);
1605      }
1606      if (numberUnsatisfied_) {
1607        // some infeasibilities - go to next steps
1608        break;
1609      } else if (!iPass) {
1610        // looks like a solution - get paranoid
1611        bool roundAgain=false;
1612        // get basis
1613        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
1614        if (!ws)
1615          break;
1616        for (i=0;i<numberColumns;i++) {
1617          double value = saveSolution[i];
1618          if (value<lower[i]) {
1619            saveSolution[i]=lower[i];
1620            roundAgain=true;
1621            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
1622          } else if (value>upper[i]) {
1623            saveSolution[i]=upper[i];
1624            roundAgain=true;
1625            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
1626          } 
1627        }
1628        if (roundAgain&&saveNumberStrong) {
1629          // restore basis
1630          solver->setWarmStart(ws);
1631          delete ws;
1632          solver->resolve();
1633          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1634          model->reserveCurrentSolution(saveSolution);
1635          if (!solver->isProvenOptimal()) {
1636            // infeasible
1637            anyAction=-2;
1638            break;
1639          }
1640        } else {
1641          delete ws;
1642          break;
1643        }
1644      }
1645    }
1646    /* Some solvers can do the strong branching calculations faster if
1647       they do them all at once.  At present only Clp does for ordinary
1648       integers but I think this coding would be easy to modify
1649    */
1650    bool allNormal=true; // to say if we can do fast strong branching
1651    // Say which one will be best
1652    int bestChoice=0;
1653    double worstInfeasibility=0.0;
1654    for (i=0;i<numberStrong;i++) {
1655      choice[i].numIntInfeasUp = numberUnsatisfied_;
1656      choice[i].numIntInfeasDown = numberUnsatisfied_;
1657      choice[i].fix=0; // say not fixed
1658      if (!dynamic_cast <const CbcSimpleInteger *> (model->object(choice[i].objectNumber)))
1659        allNormal=false; // Something odd so lets skip clever fast branching
1660      if ( !model->object(choice[i].objectNumber)->boundBranch())
1661        numberStrong=0; // switch off
1662      if ( choice[i].possibleBranch->numberBranches()>2)
1663        numberStrong=0; // switch off
1664      // Do best choice in case switched off
1665      if (choice[i].upMovement>worstInfeasibility) {
1666        worstInfeasibility=choice[i].upMovement;
1667        bestChoice=i;
1668      }
1669    }
1670    // If we have hit max time don't do strong branching
1671    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
1672                        model->getDblParam(CbcModel::CbcMaximumSeconds));
1673    // also give up if we are looping round too much
1674    if (hitMaxTime||numberPassesLeft<=0)
1675      numberStrong=0;
1676    /*
1677      Is strong branching enabled? If so, set up and do it. Otherwise, we'll
1678      fall through to simple branching.
1679     
1680      Setup for strong branching involves saving the current basis (for restoration
1681      afterwards) and setting up for hot starts.
1682    */
1683    if (numberStrong&&saveNumberStrong) {
1684     
1685      bool solveAll=false; // set true to say look at all even if some fixed (experiment)
1686      solveAll=true;
1687      // worth trying if too many times
1688      // Save basis
1689      CoinWarmStart * ws = solver->getWarmStart();
1690      // save limit
1691      int saveLimit;
1692      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
1693      if (beforeSolution&&saveLimit<100)
1694        solver->setIntParam(OsiMaxNumIterationHotStart,100); // go to end
1695#     ifdef COIN_HAS_CLP     
1696      /* If we are doing all strong branching in one go then we create new arrays
1697         to store information.  If clp NULL then doing old way.
1698         Going down -
1699         outputSolution[2*i] is final solution.
1700         outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown
1701         outputStuff[2*i+numberStrong] is number iterations
1702         On entry newUpper[i] is new upper bound, on exit obj change
1703         Going up -
1704         outputSolution[2*i+1] is final solution.
1705         outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
1706         outputStuff[2*i+1+numberStrong] is number iterations
1707         On entry newLower[i] is new lower bound, on exit obj change
1708      */
1709      ClpSimplex * clp=NULL;
1710      double * newLower = NULL;
1711      double * newUpper = NULL;
1712      double ** outputSolution=NULL;
1713      int * outputStuff=NULL;
1714      // Go back to normal way if user wants it
1715      if (osiclp&&(osiclp->specialOptions()&16)!=0&&osiclp->specialOptions()>0)
1716        allNormal=false;
1717      if (osiclp&&!allNormal) {
1718        // say do fast
1719        int easy=1;
1720        osiclp->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
1721      }
1722      if (osiclp&& allNormal) {
1723        clp = osiclp->getModelPtr();
1724        // Clp - do a different way
1725        newLower = new double[numberStrong];
1726        newUpper = new double[numberStrong];
1727        outputSolution = new double * [2*numberStrong];
1728        outputStuff = new int [4*numberStrong];
1729        int * which = new int[numberStrong];
1730        int startFinishOptions;
1731        int specialOptions = osiclp->specialOptions();
1732        int clpOptions = clp->specialOptions();
1733        int returnCode=0;
1734#define CRUNCH
1735#ifdef CRUNCH
1736        // Crunch down problem
1737        int numberRows = clp->numberRows();
1738        // Use dual region
1739        double * rhs = clp->dualRowSolution();
1740        int * whichRow = new int[3*numberRows];
1741        int * whichColumn = new int[2*numberColumns];
1742        int nBound;
1743        ClpSimplex * small = static_cast<ClpSimplexOther *> (clp)->crunch(rhs,whichRow,whichColumn,nBound,true);
1744        if (!small) {
1745          anyAction=-2;
1746          //printf("XXXX Inf by inspection\n");
1747          delete [] whichColumn;
1748          whichColumn=NULL;
1749          delete [] whichRow;
1750          whichRow=NULL;
1751          break;
1752        } else {
1753          clp = small;
1754        }
1755#else
1756        int saveLogLevel = clp->logLevel();
1757        int saveMaxIts = clp->maximumIterations();
1758#endif
1759        clp->setLogLevel(0);
1760        if((specialOptions&1)==0) {
1761          startFinishOptions=0;
1762          clp->setSpecialOptions(clpOptions|(64|1024));
1763        } else {
1764          startFinishOptions=1+2+4;
1765          //startFinishOptions=1+4; // for moment re-factorize
1766          if((specialOptions&4)==0) 
1767            clp->setSpecialOptions(clpOptions|(64|128|512|1024|4096));
1768          else
1769            clp->setSpecialOptions(clpOptions|(64|128|512|1024|2048|4096));
1770        }
1771        // User may want to clean up before strong branching
1772        if ((clp->specialOptions()&32)!=0) {
1773          clp->primal(1);
1774          if (clp->numberIterations())
1775            model->messageHandler()->message(CBC_ITERATE_STRONG,*model->messagesPointer())
1776              << clp->numberIterations()
1777              <<CoinMessageEol;
1778        }
1779        clp->setMaximumIterations(saveLimit);
1780#ifdef CRUNCH
1781        int * backColumn = whichColumn+numberColumns;
1782#endif
1783        for (i=0;i<numberStrong;i++) {
1784          int iObject = choice[i].objectNumber;
1785          const OsiObject * object = model->object(iObject);
1786          const CbcSimpleInteger * simple = static_cast <const CbcSimpleInteger *> (object);
1787          int iSequence = simple->columnNumber();
1788          newLower[i]= ceil(saveSolution[iSequence]);
1789          newUpper[i]= floor(saveSolution[iSequence]);
1790#ifdef CRUNCH
1791          iSequence = backColumn[iSequence];
1792          assert (iSequence>=0);
1793#endif
1794          which[i]=iSequence;
1795          outputSolution[2*i]= new double [numberColumns];
1796          outputSolution[2*i+1]= new double [numberColumns];
1797        }
1798        //clp->writeMps("bad");
1799        returnCode=clp->strongBranching(numberStrong,which,
1800                                        newLower, newUpper,outputSolution,
1801                                        outputStuff,outputStuff+2*numberStrong,!solveAll,false,
1802                                        startFinishOptions);
1803#ifndef CRUNCH
1804        clp->setSpecialOptions(clpOptions); // restore
1805        clp->setMaximumIterations(saveMaxIts);
1806        clp->setLogLevel(saveLogLevel);
1807#endif
1808        if (returnCode==-2) {
1809          // bad factorization!!!
1810          // Doing normal way
1811          // Mark hot start
1812          solver->markHotStart();
1813          clp = NULL;
1814        } else {
1815#ifdef CRUNCH
1816          // extract solution
1817          //bool checkSol=true;
1818          for (i=0;i<numberStrong;i++) {
1819            int iObject = choice[i].objectNumber;
1820            const OsiObject * object = model->object(iObject);
1821            const CbcSimpleInteger * simple = static_cast <const CbcSimpleInteger *> (object);
1822            int iSequence = simple->columnNumber();
1823            which[i]=iSequence;
1824            double * sol = outputSolution[2*i];
1825            double * sol2 = outputSolution[2*i+1];
1826            //bool x=true;
1827            //bool x2=true;
1828            for (int iColumn=numberColumns-1;iColumn>=0;iColumn--) {
1829              int jColumn = backColumn[iColumn];
1830              if (jColumn>=0) {
1831                sol[iColumn]=sol[jColumn];
1832                sol2[iColumn]=sol2[jColumn];
1833              } else {
1834                sol[iColumn]=saveSolution[iColumn];
1835                sol2[iColumn]=saveSolution[iColumn];
1836              }
1837            }
1838          }
1839#endif
1840        }
1841#ifdef CRUNCH
1842        delete [] whichColumn;
1843        delete [] whichRow;
1844        delete small;
1845#endif
1846        delete [] which;
1847      } else {
1848        // Doing normal way
1849        // Mark hot start
1850        solver->markHotStart();
1851      }
1852#     else      /* COIN_HAS_CLP */
1853     
1854      OsiSolverInterface *clp = NULL ;
1855      double **outputSolution = NULL ;
1856      int *outputStuff = NULL ;
1857      double * newLower = NULL ;
1858      double * newUpper = NULL ;
1859     
1860      solver->markHotStart();
1861     
1862#     endif     /* COIN_HAS_CLP */
1863      /*
1864        Open a loop to do the strong branching LPs. For each candidate variable,
1865        solve an LP with the variable forced down, then up. If a direction turns
1866        out to be infeasible or monotonic (i.e., over the dual objective cutoff),
1867        force the objective change to be big (1.0e100). If we determine the problem
1868        is infeasible, or find a monotone variable, escape the loop.
1869       
1870        TODO: The `restore bounds' part might be better encapsulated as an
1871        unbranch() method. Branching objects more exotic than simple integers
1872        or cliques might not restrict themselves to variable bounds.
1873       
1874        TODO: Virtuous solvers invalidate the current solution (or give bogus
1875        results :-) when the bounds are changed out from under them. So we
1876        need to do all the work associated with finding a new solution before
1877        restoring the bounds.
1878      */
1879      for (i = 0 ; i < numberStrong ; i++)
1880        { double objectiveChange ;
1881          double newObjectiveValue=1.0e100;
1882          // status is 0 finished, 1 infeasible and other
1883          int iStatus;
1884          /*
1885            Try the down direction first. (Specify the initial branching alternative as
1886            down with a call to way(-1). Each subsequent call to branch() performs the
1887            specified branch and advances the branch object state to the next branch
1888            alternative.)
1889          */
1890          if (!clp) {
1891            choice[i].possibleBranch->way(-1) ;
1892            choice[i].possibleBranch->branch() ;
1893            bool feasible=true;
1894            if (checkFeasibility) {
1895              // check branching did not make infeasible
1896              int iColumn;
1897              int numberColumns = solver->getNumCols();
1898              const double * columnLower = solver->getColLower();
1899              const double * columnUpper = solver->getColUpper();
1900              for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1901                if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1902                  feasible=false;
1903              }
1904            }
1905            if (feasible) {
1906              solver->solveFromHotStart() ;
1907              numberStrongDone++;
1908              numberStrongIterations += solver->getIterationCount();
1909              /*
1910                We now have an estimate of objective degradation that we can use for strong
1911                branching. If we're over the cutoff, the variable is monotone up.
1912                If we actually made it to optimality, check for a solution, and if we have
1913                a good one, call setBestSolution to process it. Note that this may reduce the
1914                cutoff, so we check again to see if we can declare this variable monotone.
1915              */
1916              if (solver->isProvenOptimal())
1917                iStatus=0; // optimal
1918              else if (solver->isIterationLimitReached()
1919                       &&!solver->isDualObjectiveLimitReached())
1920                iStatus=2; // unknown
1921              else
1922                iStatus=1; // infeasible
1923              newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1924              choice[i].numItersDown = solver->getIterationCount();
1925            } else {
1926              iStatus=1; // infeasible
1927              newObjectiveValue = 1.0e100;
1928              choice[i].numItersDown = 0;
1929            }
1930          } else {
1931            iStatus = outputStuff[2*i];
1932            choice[i].numItersDown = outputStuff[2*numberStrong+2*i];
1933            numberStrongDone++;
1934            numberStrongIterations += choice[i].numItersDown;
1935            newObjectiveValue = objectiveValue+newUpper[i];
1936            solver->setColSolution(outputSolution[2*i]);
1937          }
1938          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
1939          if (!iStatus) {
1940            choice[i].finishedDown = true ;
1941            if (newObjectiveValue>=model->getCutoff()) {
1942              objectiveChange = 1.0e100; // say infeasible
1943              numberStrongInfeasible++;
1944            } else {
1945              // See if integer solution
1946              if (model->feasibleSolution(choice[i].numIntInfeasDown,
1947                                          choice[i].numObjInfeasDown)
1948                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
1949                model->setBestSolution(CBC_STRONGSOL,
1950                                       newObjectiveValue,
1951                                       solver->getColSolution()) ;
1952                // only needed for odd solvers
1953                newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1954                objectiveChange = CoinMax(newObjectiveValue-objectiveValue_,0.0) ;
1955                model->setLastHeuristic(NULL);
1956                model->incrementUsed(solver->getColSolution());
1957                if (newObjectiveValue >= model->getCutoff()) {  //  *new* cutoff
1958                  objectiveChange = 1.0e100 ;
1959                  numberStrongInfeasible++;
1960                }
1961              }
1962            }
1963          } else if (iStatus==1) {
1964            objectiveChange = 1.0e100 ;
1965            numberStrongInfeasible++;
1966          } else {
1967            // Can't say much as we did not finish
1968            choice[i].finishedDown = false ;
1969            numberUnfinished++;
1970          }
1971          choice[i].downMovement = objectiveChange ;
1972         
1973          // restore bounds
1974          if (!clp)
1975            { for (int j=0;j<numberColumns;j++) {
1976                if (saveLower[j] != lower[j])
1977                  solver->setColLower(j,saveLower[j]);
1978                if (saveUpper[j] != upper[j])
1979                  solver->setColUpper(j,saveUpper[j]);
1980              }
1981            }
1982          //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1983          //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersDown,
1984          //     choice[i].downMovement,choice[i].finishedDown,choice[i].numIntInfeasDown,
1985          //     choice[i].numObjInfeasDown);
1986         
1987          // repeat the whole exercise, forcing the variable up
1988          if (!clp) {
1989            bool feasible=true;
1990            // If odd branching then maybe just one possibility
1991            if(choice[i].possibleBranch->numberBranchesLeft()>0) {
1992              choice[i].possibleBranch->branch();
1993              if (checkFeasibility) {
1994                // check branching did not make infeasible
1995                int iColumn;
1996                int numberColumns = solver->getNumCols();
1997                const double * columnLower = solver->getColLower();
1998                const double * columnUpper = solver->getColUpper();
1999                for (iColumn= 0;iColumn<numberColumns;iColumn++) {
2000                  if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
2001                    feasible=false;
2002                }
2003              }
2004            } else {
2005              // second branch infeasible
2006              feasible=false;
2007            }
2008            if (feasible) {
2009              solver->solveFromHotStart() ;
2010              numberStrongDone++;
2011              numberStrongIterations += solver->getIterationCount();
2012              /*
2013                We now have an estimate of objective degradation that we can use for strong
2014                branching. If we're over the cutoff, the variable is monotone up.
2015                If we actually made it to optimality, check for a solution, and if we have
2016                a good one, call setBestSolution to process it. Note that this may reduce the
2017                cutoff, so we check again to see if we can declare this variable monotone.
2018              */
2019              if (solver->isProvenOptimal())
2020                iStatus=0; // optimal
2021              else if (solver->isIterationLimitReached()
2022                       &&!solver->isDualObjectiveLimitReached())
2023                iStatus=2; // unknown
2024              else
2025                iStatus=1; // infeasible
2026              newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2027              choice[i].numItersUp = solver->getIterationCount();
2028            } else {
2029              iStatus=1; // infeasible
2030              newObjectiveValue = 1.0e100;
2031              choice[i].numItersDown = 0;
2032            }
2033          } else {
2034            iStatus = outputStuff[2*i+1];
2035            choice[i].numItersUp = outputStuff[2*numberStrong+2*i+1];
2036            numberStrongDone++;
2037            numberStrongIterations += choice[i].numItersUp;
2038            newObjectiveValue = objectiveValue+newLower[i];
2039            solver->setColSolution(outputSolution[2*i+1]);
2040          }
2041          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
2042          if (!iStatus) {
2043            choice[i].finishedUp = true ;
2044            if (newObjectiveValue>=model->getCutoff()) {
2045              objectiveChange = 1.0e100; // say infeasible
2046              numberStrongInfeasible++;
2047            } else {
2048              // See if integer solution
2049              if (model->feasibleSolution(choice[i].numIntInfeasUp,
2050                                          choice[i].numObjInfeasUp)
2051                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
2052                model->setBestSolution(CBC_STRONGSOL,
2053                                       newObjectiveValue,
2054                                       solver->getColSolution()) ;
2055                // only needed for odd solvers
2056                newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2057                objectiveChange = CoinMax(newObjectiveValue-objectiveValue_,0.0) ;
2058                model->setLastHeuristic(NULL);
2059                model->incrementUsed(solver->getColSolution());
2060                if (newObjectiveValue >= model->getCutoff()) {  //  *new* cutoff
2061                  objectiveChange = 1.0e100 ;
2062                  numberStrongInfeasible++;
2063                }
2064              }
2065            }
2066          } else if (iStatus==1) {
2067            objectiveChange = 1.0e100 ;
2068            numberStrongInfeasible++;
2069          } else {
2070            // Can't say much as we did not finish
2071            choice[i].finishedUp = false ;
2072            numberUnfinished++;
2073          }
2074          choice[i].upMovement = objectiveChange ;
2075         
2076          // restore bounds
2077          if (!clp)
2078            { for (int j=0;j<numberColumns;j++) {
2079                if (saveLower[j] != lower[j])
2080                  solver->setColLower(j,saveLower[j]);
2081                if (saveUpper[j] != upper[j])
2082                  solver->setColUpper(j,saveUpper[j]);
2083              }
2084            }
2085         
2086          //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2087          //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersUp,
2088          //     choice[i].upMovement,choice[i].finishedUp,choice[i].numIntInfeasUp,
2089          //     choice[i].numObjInfeasUp);
2090         
2091          /*
2092            End of evaluation for this candidate variable. Possibilities are:
2093            * Both sides below cutoff; this variable is a candidate for branching.
2094            * Both sides infeasible or above the objective cutoff: no further action
2095            here. Break from the evaluation loop and assume the node will be purged
2096            by the caller.
2097            * One side below cutoff: Install the branch (i.e., fix the variable). Break
2098            from the evaluation loop and assume the node will be reoptimised by the
2099            caller.
2100          */
2101          // reset
2102          choice[i].possibleBranch->resetNumberBranchesLeft();
2103          if (choice[i].upMovement<1.0e100) {
2104            if(choice[i].downMovement<1.0e100) {
2105              // feasible - no action
2106            } else {
2107              // up feasible, down infeasible
2108              anyAction=-1;
2109              //printf("Down infeasible for choice %d sequence %d\n",i,
2110              // model->object(choice[i].objectNumber)->columnNumber());
2111              if (!solveAll) {
2112                choice[i].possibleBranch->way(1);
2113                choice[i].possibleBranch->branch();
2114                break;
2115              } else {
2116                choice[i].fix=1;
2117              }
2118            }
2119          } else {
2120            if(choice[i].downMovement<1.0e100) {
2121              // down feasible, up infeasible
2122              anyAction=-1;
2123              //printf("Up infeasible for choice %d sequence %d\n",i,
2124              // model->object(choice[i].objectNumber)->columnNumber());
2125              if (!solveAll) {
2126                choice[i].possibleBranch->way(-1);
2127                choice[i].possibleBranch->branch();
2128                break;
2129              } else {
2130                choice[i].fix=-1;
2131              }
2132            } else {
2133              // neither side feasible
2134              anyAction=-2;
2135              //printf("Both infeasible for choice %d sequence %d\n",i,
2136              // model->object(choice[i].objectNumber)->columnNumber());
2137              break;
2138            }
2139          }
2140          bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
2141                              model->getDblParam(CbcModel::CbcMaximumSeconds));
2142          if (hitMaxTime) {
2143            numberStrong=i+1;
2144            break;
2145          }
2146        }
2147      if (!clp) {
2148        // Delete the snapshot
2149        solver->unmarkHotStart();
2150      } else {
2151        delete [] newLower;
2152        delete [] newUpper;
2153        delete [] outputStuff;
2154        int i;
2155        for (i=0;i<2*numberStrong;i++)
2156          delete [] outputSolution[i];
2157        delete [] outputSolution;
2158      }
2159      solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
2160      // restore basis
2161      solver->setWarmStart(ws);
2162      // Unless infeasible we will carry on
2163      // But we could fix anyway
2164      if (anyAction==-1&&solveAll) {
2165        // apply and take off
2166        for (i = 0 ; i < numberStrong ; i++) {
2167          if (choice[i].fix) {
2168            choice[i].possibleBranch->way(choice[i].fix) ;
2169            choice[i].possibleBranch->branch() ;
2170          }
2171        }
2172        bool feasible=true;
2173        if (checkFeasibility) {
2174          // check branching did not make infeasible
2175          int iColumn;
2176          int numberColumns = solver->getNumCols();
2177          const double * columnLower = solver->getColLower();
2178          const double * columnUpper = solver->getColUpper();
2179          for (iColumn= 0;iColumn<numberColumns;iColumn++) {
2180            if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
2181              feasible=false;
2182          }
2183        }
2184        if (feasible) {
2185          // can do quick optimality check
2186          int easy=2;
2187          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
2188          solver->resolve() ;
2189          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
2190          feasible = solver->isProvenOptimal();
2191        }
2192        if (feasible) {
2193          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
2194          model->reserveCurrentSolution(saveSolution);
2195          memcpy(saveLower,solver->getColLower(),numberColumns*sizeof(double));
2196          memcpy(saveUpper,solver->getColUpper(),numberColumns*sizeof(double));
2197          // Clean up all candidates whih are fixed
2198          int numberLeft=0;
2199          for (i = 0 ; i < numberStrong ; i++) {
2200            CbcStrongInfo thisChoice = choice[i];
2201            choice[i].possibleBranch=NULL;
2202            const OsiObject * object = model->object(thisChoice.objectNumber);
2203            int preferredWay;
2204            double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
2205            if (!infeasibility) {
2206              // take out
2207              delete thisChoice.possibleBranch;
2208            } else {
2209              choice[numberLeft++]=thisChoice;
2210            }
2211          }
2212          numberStrong=numberLeft;
2213          for (;i<maximumStrong;i++) {
2214            delete choice[i].possibleBranch;
2215            choice[i].possibleBranch=NULL;
2216          }
2217          // If all fixed then round again
2218          if (!numberLeft) {
2219            finished=false;
2220            numberStrong=0;
2221            saveNumberStrong=0;
2222            maximumStrong=1;
2223          } else {
2224            anyAction=0;
2225          }
2226          // If these two uncommented then different action
2227          anyAction=-1;
2228          finished=true;
2229          //printf("some fixed but continuing %d left\n",numberLeft);
2230        } else {
2231          anyAction=-2; // say infeasible
2232        }
2233      }
2234      delete ws;
2235      int numberNodes = model->getNodeCount();
2236      // update number of strong iterations etc
2237      model->incrementStrongInfo(numberStrongDone,numberStrongIterations,
2238                                 anyAction==-2 ? 0:numberStrongInfeasible,anyAction==-2);
2239     
2240      /*
2241        anyAction >= 0 indicates that strong branching didn't produce any monotone
2242        variables. Sift through the candidates for the best one.
2243       
2244        QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is
2245        local to this code block. Perhaps should be numberNodes_ from model?
2246        Unclear what this calculation is doing.
2247      */
2248      if (anyAction>=0) {
2249       
2250        // get average cost per iteration and assume stopped ones
2251        // would stop after 50% more iterations at average cost??? !!! ???
2252        double averageCostPerIteration=0.0;
2253        double totalNumberIterations=1.0;
2254        int smallestNumberInfeasibilities=COIN_INT_MAX;
2255        for (i=0;i<numberStrong;i++) {
2256          totalNumberIterations += choice[i].numItersDown +
2257            choice[i].numItersUp ;
2258          averageCostPerIteration += choice[i].downMovement +
2259            choice[i].upMovement;
2260          smallestNumberInfeasibilities= 
2261            CoinMin(CoinMin(choice[i].numIntInfeasDown ,
2262                            choice[i].numIntInfeasUp ),
2263                    smallestNumberInfeasibilities);
2264        }
2265        //if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities)
2266        //numberNodes=1000000; // switch off search for better solution
2267        numberNodes=1000000; // switch off anyway
2268        averageCostPerIteration /= totalNumberIterations;
2269        // all feasible - choose best bet
2270       
2271        // New method does all at once so it can be more sophisticated
2272        // in deciding how to balance actions.
2273        // But it does need arrays
2274        double * changeUp = new double [numberStrong];
2275        int * numberInfeasibilitiesUp = new int [numberStrong];
2276        double * changeDown = new double [numberStrong];
2277        int * numberInfeasibilitiesDown = new int [numberStrong];
2278        CbcBranchingObject ** objects = new CbcBranchingObject * [ numberStrong];
2279        for (i = 0 ; i < numberStrong ; i++) {
2280          int iColumn = choice[i].possibleBranch->variable() ;
2281          model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
2282            << i << iColumn
2283            <<choice[i].downMovement<<choice[i].numIntInfeasDown
2284            <<choice[i].upMovement<<choice[i].numIntInfeasUp
2285            <<choice[i].possibleBranch->value()
2286            <<CoinMessageEol;
2287          changeUp[i]=choice[i].upMovement;
2288          numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
2289          changeDown[i]=choice[i].downMovement;
2290          numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
2291          objects[i] = choice[i].possibleBranch;
2292        }
2293        int whichObject = decision->bestBranch(objects,numberStrong,numberUnsatisfied_,
2294                                               changeUp,numberInfeasibilitiesUp,
2295                                               changeDown,numberInfeasibilitiesDown,
2296                                               objectiveValue_);
2297        // move branching object and make sure it will not be deleted
2298        if (whichObject>=0) {
2299          branch_ = objects[whichObject];
2300          if (model->messageHandler()->logLevel()>3) 
2301            printf("Choosing column %d\n",choice[whichObject].possibleBranch->variable()) ;
2302          choice[whichObject].possibleBranch=NULL;
2303        }
2304        delete [] changeUp;
2305        delete [] numberInfeasibilitiesUp;
2306        delete [] changeDown;
2307        delete [] numberInfeasibilitiesDown;
2308        delete [] objects;
2309      }
2310#     ifdef COIN_HAS_CLP
2311      if (osiclp&&!allNormal) {
2312        // back to normal
2313        osiclp->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
2314      }
2315#     endif 
2316    }
2317    /*
2318      Simple branching. Probably just one, but we may have got here
2319      because of an odd branch e.g. a cut
2320    */
2321    else {
2322      // not strong
2323      // C) create branching object
2324      branch_ = choice[bestChoice].possibleBranch;
2325      choice[bestChoice].possibleBranch=NULL;
2326    }
2327  }
2328  // Set guessed solution value
2329  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
2330  /*
2331    Cleanup, then we're outta here.
2332  */
2333  if (!model->branchingMethod()||dynamicBranchingObject)
2334    delete decision;
2335 
2336  for (i=0;i<maximumStrong;i++)
2337    delete choice[i].possibleBranch;
2338  delete [] choice;
2339  delete [] saveLower;
2340  delete [] saveUpper;
2341 
2342  // restore solution
2343  solver->setColSolution(saveSolution);
2344  delete [] saveSolution;
2345# ifdef COIN_HAS_CLP
2346  if (osiclp) 
2347    osiclp->setSpecialOptions(saveClpOptions);
2348# endif
2349  return anyAction;
2350}
2351
2352/*
2353  Version for dynamic pseudo costs.
2354 
2355  **** For now just return if anything odd
2356  later allow even if odd
2357 
2358  The routine scans through the object list of the model looking for objects
2359  that indicate infeasibility. It tests each object using strong branching
2360  and selects the one with the least objective degradation.  A corresponding
2361  branching object is left attached to lastNode.
2362  This version gives preference in evaluation to variables which
2363  have not been evaluated many times.  It also uses numberStrong
2364  to say give up if last few tries have not changed incumbent.
2365  See Achterberg, Koch and Martin.
2366 
2367  If strong branching is disabled, a candidate object is chosen essentially
2368  at random (whatever object ends up in pos'n 0 of the candidate array).
2369 
2370  If a branching candidate is found to be monotone, bounds are set to fix the
2371  variable and the routine immediately returns (the caller is expected to
2372  reoptimize).
2373 
2374  If a branching candidate is found to result in infeasibility in both
2375  directions, the routine immediately returns an indication of infeasibility.
2376 
2377  Returns:  0   both branch directions are feasible
2378  -1    branching variable is monotone
2379  -2    infeasible
2380  -3   Use another method
2381 
2382  For now just fix on objective from strong branching.
2383*/
2384
2385int CbcNode::chooseDynamicBranch (CbcModel *model, CbcNode *lastNode,
2386                                  OsiSolverBranch * & /*branches*/,
2387                                  int numberPassesLeft)
2388 
2389{ if (lastNode)
2390    depth_ = lastNode->depth_+1;
2391  else
2392    depth_ = 0;
2393  // Go to other choose if hot start
2394  if (model->hotstartSolution()) 
2395    return -3;
2396  delete branch_;
2397  branch_=NULL;
2398  OsiSolverInterface * solver = model->solver();
2399  // get information on solver type
2400  const OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
2401  const OsiBabSolver * auxiliaryInfo = dynamic_cast<const OsiBabSolver *> (auxInfo);
2402  if (!auxiliaryInfo) {
2403    // use one from CbcModel
2404    auxiliaryInfo = model->solverCharacteristics();
2405  }
2406  int numberObjects = model->numberObjects();
2407  // If very odd set of objects then use older chooseBranch
2408  bool useOldWay = false;
2409  // point to useful information
2410  OsiBranchingInformation usefulInfo = model->usefulInformation();
2411  if (numberObjects>model->numberIntegers()) {
2412    for (int i=model->numberIntegers();i<numberObjects;i++) {
2413      OsiObject * object = model->modifiableObject(i);
2414      CbcObject * obj = dynamic_cast <CbcObject *>(object) ;
2415      if (!obj || !obj->optionalObject()) {
2416        int preferredWay;
2417        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
2418        if (infeasibility) {
2419          useOldWay=true;
2420          break;
2421        }
2422      }
2423    }
2424  }
2425  if ((model->specialOptions()&128)!=0)
2426    useOldWay=false; // allow
2427  // For now return if not simple
2428  if (useOldWay)
2429    return -3;
2430  // Modify useful info
2431  usefulInfo.depth_=depth_;
2432  if ((model->specialOptions()&128)!=0) {
2433    // SOS - shadow prices
2434    int numberRows = solver->getNumRows();
2435    const double * pi = usefulInfo.pi_;
2436    double sumPi=0.0;
2437    for (int i=0;i<numberRows;i++) 
2438      sumPi += fabs(pi[i]);
2439    sumPi /= static_cast<double> (numberRows);
2440    // and scale back
2441    sumPi *= 0.01;
2442    usefulInfo.defaultDual_ = sumPi; // switch on
2443    int numberColumns = solver->getNumCols();
2444    int size = CoinMax(numberColumns,2*numberRows);
2445    usefulInfo.usefulRegion_ = new double [size];
2446    CoinZeroN(usefulInfo.usefulRegion_,size);
2447    usefulInfo.indexRegion_ = new int [size];
2448    // pi may change
2449    usefulInfo.pi_=CoinCopyOfArray(usefulInfo.pi_,numberRows);
2450  }
2451  assert (auxiliaryInfo);
2452  double cutoff =model->getCutoff();
2453  const double * lower = solver->getColLower();
2454  const double * upper = solver->getColUpper();
2455  // See if user thinks infeasible
2456  int anyAction=model->problemFeasibility()->feasible(model,0);
2457  if (anyAction) {
2458    // will return -2 if infeasible , 0 if treat as integer
2459    return anyAction-1;
2460  }
2461  int i;
2462  int saveStateOfSearch = model->stateOfSearch()%10;
2463  int numberStrong=model->numberStrong();
2464  /* Ranging is switched off.
2465     The idea is that you can find out the effect of one iteration
2466     on each unsatisfied variable cheaply.  Then use this
2467     if you have not got much else to go on.
2468  */
2469  //#define RANGING
2470#ifdef RANGING
2471  // must have clp
2472#ifndef COIN_HAS_CLP
2473#  warning("Ranging switched off as not Clp");
2474#undef RANGING
2475#endif
2476  // Pass number
2477  int kPass=0;
2478  int numberRows = solver->getNumRows();
2479#endif
2480  int numberColumns = model->getNumCols();
2481  double * saveUpper = new double[numberColumns];
2482  double * saveLower = new double[numberColumns];
2483  for (i=0;i<numberColumns;i++) {
2484    saveLower[i] = lower[i];
2485    saveUpper[i] = upper[i];
2486  }
2487 
2488  // Save solution in case heuristics need good solution later
2489 
2490  double * saveSolution = new double[numberColumns];
2491  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
2492  model->reserveCurrentSolution(saveSolution);
2493  /*
2494    Get a branching decision object. Use the default dynamic decision criteria unless
2495    the user has loaded a decision method into the model.
2496  */
2497  CbcBranchDecision *decision = model->branchingMethod();
2498  if (!decision)
2499    decision = new CbcBranchDynamicDecision();
2500  int xMark=0;
2501  // Get arrays to sort
2502  double * sort = new double[numberObjects];
2503  int * whichObject = new int[numberObjects];
2504#ifdef RANGING
2505  int xPen=0;
2506  int * objectMark = new int[2*numberObjects+1];
2507#endif
2508  // Arrays with movements
2509  double * upEstimate = new double[numberObjects];
2510  double * downEstimate = new double[numberObjects];
2511  double estimatedDegradation=0.0; 
2512  int numberNodes=model->getNodeCount();
2513  int saveLogLevel = model->logLevel();
2514#if 0
2515  if ((numberNodes%500)==0) {
2516    model->setLogLevel(6);
2517    // Get average up and down costs
2518    double averageUp=0.0;
2519    double averageDown=0.0;
2520    int numberUp=0;
2521    int numberDown=0;
2522    int i;
2523    for ( i=0;i<numberObjects;i++) {
2524      OsiObject * object = model->modifiableObject(i);
2525      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2526        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2527      assert(dynamicObject);
2528      int  numberUp2=0;
2529      int numberDown2=0;
2530      double up=0.0;
2531      double down=0.0;
2532      if (dynamicObject->numberTimesUp()) {
2533        numberUp++;
2534        averageUp += dynamicObject->upDynamicPseudoCost();
2535        numberUp2 += dynamicObject->numberTimesUp();
2536        up = dynamicObject->upDynamicPseudoCost();
2537      }
2538      if (dynamicObject->numberTimesDown()) {
2539        numberDown++;
2540        averageDown += dynamicObject->downDynamicPseudoCost();
2541        numberDown2 += dynamicObject->numberTimesDown();
2542        down = dynamicObject->downDynamicPseudoCost();
2543      }
2544      if (numberUp2||numberDown2)
2545        printf("col %d - up %d times cost %g, - down %d times cost %g\n",
2546               dynamicObject->columnNumber(),numberUp2,up,numberDown2,down);
2547    }
2548    if (numberUp)
2549      averageUp /= static_cast<double> (numberUp);
2550    else
2551      averageUp=1.0;
2552    if (numberDown)
2553      averageDown /= static_cast<double> (numberDown);
2554    else
2555      averageDown=1.0;
2556    printf("total - up %d vars average %g, - down %d vars average %g\n",
2557           numberUp,averageUp,numberDown,averageDown);
2558  }
2559#endif
2560  int numberBeforeTrust = model->numberBeforeTrust();
2561  // May go round twice if strong branching fixes all local candidates
2562  bool finished=false;
2563  int numberToFix=0;
2564# ifdef COIN_HAS_CLP
2565  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
2566  int saveClpOptions=0;
2567  if (osiclp) {
2568    // for faster hot start
2569    saveClpOptions = osiclp->specialOptions();
2570    osiclp->setSpecialOptions(saveClpOptions|8192);
2571  }
2572# else
2573  OsiSolverInterface *osiclp = NULL ;
2574# endif
2575  const CglTreeProbingInfo * probingInfo = NULL; //model->probingInfo();
2576  // Old code left in with DEPRECATED_STRATEGY
2577  assert (model->searchStrategy()==-1||
2578          model->searchStrategy()==1||
2579          model->searchStrategy()==2);
2580#ifdef DEPRECATED_STRATEGY
2581  int saveSearchStrategy2 = model->searchStrategy();
2582#endif
2583  // Get average up and down costs
2584  double averageUp=0.0;
2585  double averageDown=0.0;
2586  {
2587    int numberUp=0;
2588    int numberDown=0;
2589    int i;
2590    for ( i=0;i<numberObjects;i++) {
2591      OsiObject * object = model->modifiableObject(i);
2592      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2593        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2594      if (dynamicObject) {
2595        if (dynamicObject->numberTimesUp()) {
2596          numberUp++;
2597          averageUp += dynamicObject->upDynamicPseudoCost();
2598        }
2599        if (dynamicObject->numberTimesDown()) {
2600          numberDown++;
2601          averageDown += dynamicObject->downDynamicPseudoCost();
2602        }
2603      }
2604    }
2605    if (numberUp) 
2606      averageUp /= static_cast<double> (numberUp);
2607    else
2608      averageUp=1.0;
2609    if (numberDown) 
2610      averageDown /= static_cast<double> (numberDown);
2611    else
2612      averageDown=1.0;
2613    for ( i=0;i<numberObjects;i++) {
2614      OsiObject * object = model->modifiableObject(i);
2615      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2616        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2617      if (dynamicObject) {
2618        if (!dynamicObject->numberTimesUp()) 
2619          dynamicObject->setUpDynamicPseudoCost(averageUp);
2620        if (!dynamicObject->numberTimesDown()) 
2621          dynamicObject->setDownDynamicPseudoCost(averageDown);
2622      }
2623    }
2624  }
2625#ifdef DEPRECATED_STRATEGY
2626  { // in for tabbing
2627  } else if (saveSearchStrategy2<1999) {
2628    // pseudo shadow prices
2629    model->pseudoShadow(NULL,NULL);
2630  } else if (saveSearchStrategy2<2999) {
2631    // leave old ones
2632  } else if (saveSearchStrategy2<3999) {
2633    // pseudo shadow prices at root
2634    if (!numberNodes)
2635      model->pseudoShadow(NULL,NULL);
2636  } else {
2637    abort();
2638  }
2639  if (saveSearchStrategy2>=0)
2640    saveSearchStrategy2 = saveSearchStrategy2 % 1000;
2641  if (saveSearchStrategy2==999)
2642    saveSearchStrategy2=-1;
2643  int saveSearchStrategy = saveSearchStrategy2<99 ? saveSearchStrategy2 : saveSearchStrategy2-100;
2644#endif //DEPRECATED_STRATEGY
2645  int numberNotTrusted=0;
2646  int numberStrongDone=0;
2647  int numberUnfinished=0;
2648  int numberStrongInfeasible=0;
2649  int numberStrongIterations=0;
2650  // so we can save lots of stuff
2651  CbcStrongInfo choice;
2652  CbcDynamicPseudoCostBranchingObject * choiceObject = NULL;
2653  if (model->allDynamic()) {
2654    CbcSimpleIntegerDynamicPseudoCost * object = NULL;
2655    choiceObject=new CbcDynamicPseudoCostBranchingObject(model,0,-1,0.5,object);
2656  }
2657  choice.possibleBranch=choiceObject;
2658  numberPassesLeft = CoinMax(numberPassesLeft,2);
2659  while(!finished) {
2660    numberPassesLeft--;
2661    finished=true;
2662    decision->initialize(model);
2663    // Some objects may compute an estimate of best solution from here
2664    estimatedDegradation=0.0; 
2665    numberToFix=0;
2666    int numberToDo=0;
2667    int iBestNot=-1;
2668    int iBestGot=-1;
2669    double best=0.0;
2670    numberNotTrusted=0;
2671    numberStrongDone=0;
2672    numberUnfinished=0;
2673    numberStrongInfeasible=0;
2674    numberStrongIterations=0;
2675#ifdef RANGING
2676    int * which = objectMark+numberObjects+1;
2677    int neededPenalties;
2678    int optionalPenalties;
2679#endif
2680    // We may go round this loop three times (only if we think we have solution)
2681    for (int iPass=0;iPass<3;iPass++) {
2682     
2683      // Some objects may compute an estimate of best solution from here
2684      estimatedDegradation=0.0; 
2685      numberUnsatisfied_ = 0;
2686      // initialize sum of "infeasibilities"
2687      sumInfeasibilities_ = 0.0;
2688      int bestPriority=COIN_INT_MAX;
2689      int number01 = 0;
2690      const cliqueEntry * entry = NULL;
2691      const int * toZero = NULL;
2692      const int * toOne = NULL;
2693      const int * backward = NULL;
2694      int numberUnsatisProbed=0;
2695      int numberUnsatisNotProbed=0; // 0-1
2696      if (probingInfo) {
2697        number01 = probingInfo->numberIntegers();
2698        entry = probingInfo->fixEntries();
2699        toZero = probingInfo->toZero();
2700        toOne = probingInfo->toOne();
2701        backward = probingInfo->backward();
2702        if (!toZero[number01]||number01<numberObjects||true) {
2703          // no info
2704          probingInfo=NULL;
2705        }
2706      }
2707      /*
2708        Scan for branching objects that indicate infeasibility. Choose candidates
2709        using priority as the first criteria, then integer infeasibility.
2710       
2711        The algorithm is to fill the array with a set of good candidates (by
2712        infeasibility) with priority bestPriority.  Finding a candidate with
2713        priority better (less) than bestPriority flushes the choice array. (This
2714        serves as initialization when the first candidate is found.)
2715       
2716      */
2717      numberToDo=0;
2718#ifdef RANGING
2719      neededPenalties=0;
2720      optionalPenalties=numberObjects;
2721#endif
2722      iBestNot=-1;
2723      double bestNot=0.0;
2724      iBestGot=-1;
2725      best=0.0;
2726      /* Problem type as set by user or found by analysis.  This will be extended
2727         0 - not known
2728         1 - Set partitioning <=
2729         2 - Set partitioning ==
2730         3 - Set covering
2731         4 - all +- 1 or all +1 and odd
2732      */
2733      int problemType = model->problemType();
2734      for (i=0;i<numberObjects;i++) {
2735        OsiObject * object = model->modifiableObject(i);
2736        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2737          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2738        int preferredWay;
2739        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
2740        int priorityLevel = object->priority();
2741#define ZERO_ONE 0
2742#define ZERO_FAKE 1.0e20;
2743#if ZERO_ONE==1
2744        // branch on 0-1 first (temp)
2745        if (fabs(saveSolution[dynamicObject->columnNumber()])<1.0)
2746          priorityLevel--;
2747#endif
2748#if ZERO_ONE==2
2749        if (fabs(saveSolution[dynamicObject->columnNumber()])<1.0)
2750          infeasibility *= ZERO_FAKE;
2751#endif
2752        if (infeasibility) {
2753          int iColumn = numberColumns+i;
2754          bool gotDown=false;
2755          int numberThisDown = 0;
2756          bool gotUp=false;
2757          int numberThisUp = 0;
2758          double downGuess=object->downEstimate();
2759          double upGuess=object->upEstimate();
2760          if (dynamicObject) {
2761            // Use this object's numberBeforeTrust
2762            int numberBeforeTrust = dynamicObject->numberBeforeTrust();
2763            iColumn = dynamicObject->columnNumber();
2764            gotDown=false;
2765            numberThisDown = dynamicObject->numberTimesDown();
2766            if (numberThisDown>=numberBeforeTrust)
2767              gotDown=true;
2768            gotUp=false;
2769            numberThisUp = dynamicObject->numberTimesUp();
2770            if (numberThisUp>=numberBeforeTrust)
2771              gotUp=true;
2772            if (!depth_&&false) {
2773              // try closest to 0.5
2774              double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2775              infeasibility = fabs(0.5-part);
2776            }
2777            if (problemType>0&&problemType<4&&false) {
2778              // try closest to 0.5
2779              double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2780              infeasibility = 0.5-fabs(0.5-part);
2781            }
2782            if (probingInfo) {
2783              int iSeq = backward[iColumn];
2784              assert (iSeq>=0);
2785              infeasibility = 1.0 + (toZero[iSeq+1]-toZero[iSeq])+ 
2786                5.0*CoinMin(toOne[iSeq]-toZero[iSeq],toZero[iSeq+1]-toOne[iSeq]);
2787              if (toZero[iSeq+1]>toZero[iSeq]) {
2788                numberUnsatisProbed++;
2789              } else {
2790                numberUnsatisNotProbed++;
2791              }
2792            }
2793          } else {
2794            // see if SOS
2795            CbcSOS * sosObject =
2796              dynamic_cast <CbcSOS *>(object) ;
2797            if (sosObject) {
2798              gotDown=false;
2799              numberThisDown = sosObject->numberTimesDown();
2800              if (numberThisDown>=numberBeforeTrust)
2801                gotDown=true;
2802              gotUp=false;
2803              numberThisUp = sosObject->numberTimesUp();
2804              if (numberThisUp>=numberBeforeTrust)
2805                gotUp=true;
2806            } else {
2807              gotDown=true;
2808              numberThisDown=999999;
2809              downGuess=1.0e20;
2810              gotUp=true;
2811              numberThisUp=999999;
2812              upGuess=1.0e20;
2813              numberPassesLeft=0;
2814            }
2815          }
2816          // Increase estimated degradation to solution
2817          estimatedDegradation += CoinMin(downGuess,upGuess);
2818          downEstimate[i]=downGuess;
2819          upEstimate[i]=upGuess;
2820          numberUnsatisfied_++;
2821          sumInfeasibilities_ += infeasibility;
2822          // Better priority? Flush choices.
2823          if (priorityLevel<bestPriority) {
2824            numberToDo=0;
2825            bestPriority = priorityLevel;
2826            iBestGot=-1;
2827            best=0.0;
2828            numberNotTrusted=0;
2829#ifdef RANGING
2830            neededPenalties=0;
2831            optionalPenalties=numberObjects;
2832#endif
2833          } else if (priorityLevel>bestPriority) {
2834            continue;
2835          }
2836          if (!gotUp||!gotDown)
2837            numberNotTrusted++;
2838          // Check for suitability based on infeasibility.
2839          if ((gotDown&&gotUp)&&numberStrong>0) {
2840            sort[numberToDo]=-infeasibility;
2841            if (infeasibility>best) {
2842              best=infeasibility;
2843              iBestGot=numberToDo;
2844            }
2845#ifdef RANGING
2846            if (dynamicObject) {
2847              objectMark[--optionalPenalties]=numberToDo;
2848              which[optionalPenalties]=iColumn;
2849            }
2850#endif
2851          } else {
2852#ifdef RANGING
2853            if (dynamicObject) {
2854              objectMark[neededPenalties]=numberToDo;
2855              which[neededPenalties++]=iColumn;
2856            }
2857#endif
2858            sort[numberToDo]=-10.0*infeasibility;
2859            if (!(numberThisUp+numberThisDown))
2860              sort[numberToDo] *= 100.0; // make even more likely
2861            if (iColumn<numberColumns) {
2862              double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2863              if (1.0-fabs(part-0.5)>bestNot) {
2864                iBestNot=numberToDo;
2865                bestNot = 1.0-fabs(part-0.5);
2866              }
2867            } else {
2868              // SOS
2869              if (-sort[numberToDo]>bestNot) {
2870                iBestNot=numberToDo;
2871                bestNot = -sort[numberToDo];
2872              }
2873            }
2874          }
2875          if (model->messageHandler()->logLevel()>3) { 
2876            printf("%d (%d) down %d %g up %d %g - infeas %g - sort %g solution %g\n",
2877                   i,iColumn,numberThisDown,object->downEstimate(),numberThisUp,object->upEstimate(),
2878                   infeasibility,sort[numberToDo],saveSolution[iColumn]);
2879          }
2880          whichObject[numberToDo++]=i;
2881        } else {
2882          // for debug
2883          downEstimate[i]=-1.0;
2884          upEstimate[i]=-1.0;
2885        }
2886      }
2887      if (numberUnsatisfied_) {
2888        if (probingInfo&&false)
2889          printf("nunsat %d, %d probed, %d other 0-1\n",numberUnsatisfied_,
2890                 numberUnsatisProbed,numberUnsatisNotProbed);
2891        // some infeasibilities - go to next steps
2892        break;
2893      } else if (!iPass) {
2894        // may just need resolve
2895        model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
2896        if (!solver->isProvenOptimal()) {
2897          // infeasible
2898          anyAction=-2;
2899          break;
2900        }
2901      } else if (iPass==1) {
2902        // looks like a solution - get paranoid
2903        bool roundAgain=false;
2904        // get basis
2905        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
2906        if (!ws)
2907          break;
2908        double tolerance;
2909        solver->getDblParam(OsiPrimalTolerance,tolerance);
2910        for (i=0;i<numberColumns;i++) {
2911          double value = saveSolution[i];
2912          if (value<lower[i]-tolerance) {
2913            saveSolution[i]=lower[i];
2914            roundAgain=true;
2915            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
2916          } else if (value>upper[i]+tolerance) {
2917            saveSolution[i]=upper[i];
2918            roundAgain=true;
2919            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
2920          } 
2921        }
2922        if (roundAgain) {
2923          // restore basis
2924          solver->setWarmStart(ws);
2925          solver->setColSolution(saveSolution);
2926          delete ws;
2927          bool takeHint;
2928          OsiHintStrength strength;
2929          solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
2930          solver->setHintParam(OsiDoDualInResolve,false,OsiHintDo) ;
2931          model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
2932          solver->setHintParam(OsiDoDualInResolve,takeHint,strength) ;
2933          if (!solver->isProvenOptimal()) {
2934            // infeasible
2935            anyAction=-2;
2936            break;
2937          }
2938        } else {
2939          delete ws;
2940          break;
2941        }
2942      }
2943    }
2944    if (anyAction==-2) {
2945      break;
2946    }
2947    // skip if solution
2948    if (!numberUnsatisfied_)
2949      break;
2950    int skipAll = (numberNotTrusted==0||numberToDo==1) ? 1 : 0;
2951    bool doneHotStart=false;
2952    //DEPRECATED_STRATEGYint searchStrategy = saveSearchStrategy>=0 ? (saveSearchStrategy%10) : -1;
2953    int searchStrategy = model->searchStrategy();
2954    // But adjust depending on ratio of iterations
2955    if (searchStrategy>0) {
2956      if (numberBeforeTrust>=5&&numberBeforeTrust<=10) {
2957        if (searchStrategy!=2) {
2958          assert (searchStrategy==1);
2959          if (depth_>5) {
2960            int numberIterations = model->getIterationCount();
2961            int numberStrongIterations = model->numberStrongIterations();
2962            if (numberStrongIterations>numberIterations+10000) {
2963              searchStrategy=2;
2964              skipAll=1;
2965            } else if (numberStrongIterations*4+1000<numberIterations) {
2966              searchStrategy=3;
2967              skipAll=0;
2968            }
2969          } else {
2970            searchStrategy=3;
2971            skipAll=0;
2972          }
2973        }
2974      }
2975    }
2976    // worth trying if too many times
2977    // Save basis
2978    CoinWarmStart * ws = NULL;
2979    // save limit
2980    int saveLimit=0;
2981    solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
2982    if (!numberPassesLeft)
2983      skipAll=1;
2984    if (!skipAll) {
2985      ws = solver->getWarmStart();
2986      int limit=100;
2987      if (!saveStateOfSearch&&saveLimit<limit&&saveLimit==100)
2988        solver->setIntParam(OsiMaxNumIterationHotStart,limit); 
2989    }
2990    // Say which one will be best
2991    int whichChoice=0;
2992    int bestChoice;
2993    if (iBestGot>=0)
2994      bestChoice=iBestGot;
2995    else
2996      bestChoice=iBestNot;
2997    assert (bestChoice>=0);
2998    // If we have hit max time don't do strong branching
2999    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
3000                        model->getDblParam(CbcModel::CbcMaximumSeconds));
3001    // also give up if we are looping round too much
3002    if (hitMaxTime||numberPassesLeft<=0) {
3003      int iObject = whichObject[bestChoice];
3004      OsiObject * object = model->modifiableObject(iObject);
3005      int preferredWay;
3006      object->infeasibility(&usefulInfo,preferredWay);
3007      CbcObject * obj =
3008        dynamic_cast <CbcObject *>(object) ;
3009      assert (obj);
3010      branch_=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3011      {
3012        CbcBranchingObject * branchObj =
3013          dynamic_cast <CbcBranchingObject *>(branch_) ;
3014        assert (branchObj);
3015        branchObj->way(preferredWay);
3016      }
3017      delete ws;
3018      ws=NULL;
3019      break;
3020    } else {
3021      // say do fast
3022      int easy=1;
3023      if (!skipAll)
3024        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
3025      int iDo;
3026#ifdef RANGING
3027      bool useRanging = model->allDynamic()&&!skipAll;
3028      if (useRanging) {
3029        double currentObjective = solver->getObjValue()*solver->getObjSense();
3030        double gap = cutoff-currentObjective;
3031        // relax a bit
3032        gap *= 1.0000001;
3033        gap = CoinMax(1.0e-5,gap);
3034        // off penalties if too much
3035        double needed = neededPenalties;
3036        needed *= numberRows;
3037        if (numberNodes) {
3038          if (needed>1.0e6) {
3039            neededPenalties=0;
3040          } else if (gap<1.0e5) {
3041            // maybe allow some not needed
3042            int extra = static_cast<int> ((1.0e6-needed)/numberRows);
3043            int nStored = numberObjects-optionalPenalties;
3044            extra = CoinMin(extra,nStored);
3045            for (int i=0;i<extra;i++) {
3046              objectMark[neededPenalties]=objectMark[optionalPenalties+i];
3047                which[neededPenalties++]=which[optionalPenalties+i];;
3048            }
3049          }
3050        }
3051        if (osiclp&&neededPenalties) {
3052          assert (!doneHotStart);
3053          xPen += neededPenalties;
3054          which--;
3055          which[0]=neededPenalties;
3056          osiclp->passInRanges(which);
3057          // Mark hot start and get ranges
3058          if (kPass) {
3059            // until can work out why solution can go funny
3060            int save = osiclp->specialOptions();
3061            osiclp->setSpecialOptions(save|256);
3062            solver->markHotStart();
3063            osiclp->setSpecialOptions(save);
3064          } else {
3065            solver->markHotStart();
3066          }
3067          doneHotStart=true;
3068          xMark++;
3069          kPass++;
3070          osiclp->passInRanges(NULL);
3071          const double * downCost=osiclp->upRange();
3072          const double * upCost=osiclp->downRange();
3073          bool problemFeasible=true;
3074          int numberFixed=0;
3075          for (int i=0;i<neededPenalties;i++) {
3076            int j = objectMark[i];
3077            int iObject = whichObject[j];
3078            OsiObject * object = model->modifiableObject(iObject);
3079            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3080              static_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3081            // Use this object's numberBeforeTrust
3082            int numberBeforeTrust = dynamicObject->numberBeforeTrust();
3083            int iSequence=dynamicObject->columnNumber();
3084            double value = saveSolution[iSequence];
3085            value -= floor(value);
3086            double upPenalty = CoinMin(upCost[i],1.0e110)*(1.0-value);
3087            double downPenalty = CoinMin(downCost[i],1.0e110)*value;
3088            int numberThisDown = dynamicObject->numberTimesDown();
3089            int numberThisUp = dynamicObject->numberTimesUp();
3090            if (!numberBeforeTrust) {
3091              // override
3092              downEstimate[iObject]=downPenalty;
3093              upEstimate[iObject]=upPenalty;
3094              double min1 = CoinMin(downEstimate[iObject],
3095                                    upEstimate[iObject]);
3096              double max1 = CoinMax(downEstimate[iObject],
3097                                    upEstimate[iObject]);
3098              min1 = 0.8*min1+0.2*max1;
3099              sort[j] = - min1;
3100            } else if (numberThisDown<numberBeforeTrust||
3101                       numberThisUp<numberBeforeTrust) {
3102              double invTrust = 1.0/static_cast<double> (numberBeforeTrust);
3103              if (numberThisDown<numberBeforeTrust) {
3104                double fraction = numberThisDown*invTrust;
3105                downEstimate[iObject] = fraction*downEstimate[iObject]+(1.0-fraction)*downPenalty;
3106              }
3107              if (numberThisUp<numberBeforeTrust) {
3108                double fraction = numberThisUp*invTrust;
3109                upEstimate[iObject] = fraction*upEstimate[iObject]+(1.0-fraction)*upPenalty;
3110              }
3111              double min1 = CoinMin(downEstimate[iObject],
3112                                    upEstimate[iObject]);
3113              double max1 = CoinMax(downEstimate[iObject],
3114                                    upEstimate[iObject]);
3115              min1 = 0.8*min1+0.2*max1;
3116              min1 *= 10.0;
3117              if (!(numberThisDown+numberThisUp))
3118                min1 *= 100.0;
3119              sort[j] = - min1;
3120            }
3121            if (CoinMax(downPenalty,upPenalty)>gap) {
3122              printf("gap %g object %d has down range %g, up %g\n",
3123                     gap,i,downPenalty,upPenalty);
3124              //sort[j] -= 1.0e50; // make more likely to be chosen
3125              int number;
3126              if (downPenalty>gap) {
3127                number = dynamicObject->numberTimesDown();
3128                if (upPenalty>gap) 
3129                  problemFeasible=false;
3130                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver,&usefulInfo,1);
3131                //branch->fix(solver,saveLower,saveUpper,1);
3132                delete branch;
3133              } else {
3134                number = dynamicObject->numberTimesUp();
3135                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver,&usefulInfo,1);
3136                //branch->fix(solver,saveLower,saveUpper,-1);
3137                delete branch;
3138              }
3139              if(number>=numberBeforeTrust)
3140                dynamicObject->setNumberBeforeTrust(number+1);
3141              numberFixed++;
3142            }
3143            if (!numberNodes)
3144              printf("%d pen down ps %g -> %g up ps %g -> %g\n",
3145                     iObject,downPenalty,downPenalty,upPenalty,upPenalty);
3146          }
3147          if(numberFixed&&problemFeasible) {
3148            assert(doneHotStart);
3149            solver->unmarkHotStart();
3150            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3151            solver->markHotStart();
3152            problemFeasible = solver->isProvenOptimal();
3153          }
3154          if(!problemFeasible) {
3155            fprintf(stdout,"both ways infeas on ranging - code needed\n");
3156            anyAction=-2;
3157            if (!choiceObject) {
3158              delete choice.possibleBranch;
3159              choice.possibleBranch=NULL;
3160            }
3161            //printf("Both infeasible for choice %d sequence %d\n",i,
3162            // model->object(choice.objectNumber)->columnNumber());
3163            // Delete the snapshot
3164            solver->unmarkHotStart();
3165            // back to normal
3166            solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3167            // restore basis
3168            solver->setWarmStart(ws);
3169            doneHotStart=false;
3170            delete ws;
3171            ws=NULL;
3172            break;
3173          }
3174        }
3175      }
3176#endif          /* RANGING */
3177      {
3178        int numberIterations = model->getIterationCount();
3179        //numberIterations += (model->numberExtraIterations()>>2);
3180        const int * strongInfo = model->strongInfo();
3181        //int numberDone = strongInfo[0]-strongInfo[3];
3182        int numberFixed = strongInfo[1]-strongInfo[4];
3183        int numberInfeasible = strongInfo[2]-strongInfo[5];
3184        assert (!strongInfo[3]);
3185        assert (!strongInfo[4]);
3186        assert (!strongInfo[5]);
3187        int numberStrongIterations = model->numberStrongIterations();
3188        int numberRows = solver->getNumRows();
3189        if (numberStrongIterations>numberIterations+CoinMin(100,10*numberRows)&&depth_>=4&&numberNodes>100) {
3190          if (20*numberInfeasible+4*numberFixed<numberNodes) {
3191            // Say never do
3192            skipAll=-1;
3193          }
3194        } 
3195      }
3196      // make sure best will be first
3197      if (iBestGot>=0)
3198        sort[iBestGot]=-COIN_DBL_MAX;
3199      // Actions 0 - exit for repeat, 1 resolve and try old choice,2 exit for continue
3200      if (anyAction)
3201        numberToDo=0; // skip as we will be trying again
3202      // Sort
3203      CoinSort_2(sort,sort+numberToDo,whichObject);
3204      // Change in objective opposite infeasible
3205      double worstFeasible=0.0;
3206      // Just first if strong off
3207      if (!numberStrong)
3208        numberToDo=CoinMin(numberToDo,1);
3209      if (searchStrategy==2)
3210        numberToDo=CoinMin(numberToDo,20);
3211      iDo=0;
3212      int saveLimit2;
3213      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit2);
3214      bool doQuickly = false; // numberToDo>2*numberStrong;
3215      if (searchStrategy==2)
3216        doQuickly=true;
3217      int numberTest=numberNotTrusted>0 ? numberStrong : (numberStrong+1)/2;
3218      if (searchStrategy==3) {
3219        // Previously decided we need strong
3220        doQuickly=false;
3221        numberTest = numberStrong;
3222      }
3223      // Try nearly always off
3224      if (skipAll>=0) {
3225        if (searchStrategy<2) {
3226          //if ((numberNodes%20)!=0) {
3227            if ((model->specialOptions()&8)==0) {
3228              numberTest=0;
3229              doQuickly=true;
3230            }
3231            //} else {
3232            //doQuickly=false;
3233            //numberTest=2*numberStrong;
3234            //skipAll=0;
3235            //}
3236        }
3237      } else {
3238        // Just take first
3239        doQuickly=true;
3240        numberTest=1;
3241      }
3242      int testDepth = (skipAll>=0) ? 8 : 4;
3243      if (depth_<testDepth&&numberStrong) {
3244        if (searchStrategy!=2) {
3245          doQuickly=false;
3246          int numberRows = solver->getNumRows();
3247          // whether to do this or not is important - think
3248          if (numberRows<300||numberRows+numberColumns<2500) {
3249            if (depth_<7)
3250              numberStrong = CoinMin(3*numberStrong,numberToDo);
3251            if (!depth_) 
3252              numberStrong=CoinMin(6*numberStrong,numberToDo);
3253          }
3254          numberTest=numberStrong;
3255          skipAll=0;
3256        }
3257      }
3258      // Do at least 5 strong
3259      if (numberColumns<1000&&(depth_<15||numberNodes<1000000))
3260        numberTest = CoinMax(numberTest,5);
3261      if ((model->specialOptions()&8)==0) {
3262        if (skipAll) {
3263          numberTest=0;
3264          doQuickly=true;
3265        }
3266      } else {
3267        // do 5 as strong is fixing
3268        numberTest = CoinMax(numberTest,5);
3269      }
3270      // see if switched off
3271      if (skipAll<0) {
3272        numberTest=0;
3273        doQuickly=true;
3274      }
3275      int realMaxHotIterations=999999;
3276      if (skipAll<0)
3277        numberToDo=1;
3278#ifdef DO_ALL_AT_ROOT
3279      if (!numberNodes) {
3280        printf("DOX %d test %d unsat %d - nobj %d\n",
3281               numberToDo,numberTest,numberUnsatisfied_,
3282               numberObjects);
3283        numberTest=numberToDo;
3284      }
3285#endif
3286      for ( iDo=0;iDo<numberToDo;iDo++) {
3287        int iObject = whichObject[iDo];
3288        OsiObject * object = model->modifiableObject(iObject);
3289        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3290          static_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3291        int iColumn = dynamicObject ? dynamicObject->columnNumber() : numberColumns+iObject;
3292        int preferredWay;
3293        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
3294        // may have become feasible
3295        if (!infeasibility)
3296          continue;
3297        CbcSimpleInteger * obj =
3298          dynamic_cast <CbcSimpleInteger *>(object) ;
3299        if (obj) {
3300          if (choiceObject) {
3301            obj->fillCreateBranch(choiceObject,&usefulInfo,preferredWay);
3302            choiceObject->setObject(dynamicObject);
3303          } else {
3304            choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3305          }
3306        } else {
3307          CbcObject * obj =
3308            dynamic_cast <CbcObject *>(object) ;
3309          assert (obj);
3310          choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3311        }
3312        // Save which object it was
3313        choice.objectNumber=iObject;
3314        choice.numIntInfeasUp = numberUnsatisfied_;
3315        choice.numIntInfeasDown = numberUnsatisfied_;
3316        choice.upMovement = upEstimate[iObject];
3317        choice.downMovement = downEstimate[iObject];
3318        assert (choice.upMovement>=0.0);
3319        assert (choice.downMovement>=0.0);
3320        choice.fix=0; // say not fixed
3321        // see if can skip strong branching
3322        int canSkip = choice.possibleBranch->fillStrongInfo(choice);
3323        if ((numberTest<=0||skipAll)) {
3324          if (iDo>20) {
3325#ifdef DO_ALL_AT_ROOT
3326            if (!numberNodes)
3327              printf("DOY test %d - iDo %d\n",
3328                     numberTest,iDo);
3329#endif
3330            if (!choiceObject) {
3331              delete choice.possibleBranch;
3332              choice.possibleBranch=NULL;
3333            }
3334            break; // give up anyway
3335          }
3336        }
3337        if (model->messageHandler()->logLevel()>3&&numberBeforeTrust&&dynamicObject) 
3338          dynamicObject->print(1,choice.possibleBranch->value());
3339        if (skipAll<0)
3340          canSkip=true;
3341        if (!canSkip) {
3342          if (!doneHotStart) {
3343            // Mark hot start
3344            doneHotStart=true;
3345            solver->markHotStart();
3346            xMark++;
3347          }
3348          numberTest--;
3349          // just do a few
3350          if (searchStrategy==2)
3351            solver->setIntParam(OsiMaxNumIterationHotStart,10); 
3352          double objectiveChange ;
3353          double newObjectiveValue=1.0e100;
3354          int j;
3355          // status is 0 finished, 1 infeasible and other
3356          int iStatus;
3357          /*
3358            Try the down direction first. (Specify the initial branching alternative as
3359            down with a call to way(-1). Each subsequent call to branch() performs the
3360            specified branch and advances the branch object state to the next branch
3361            alternative.)
3362          */
3363          choice.possibleBranch->way(-1) ;
3364          choice.possibleBranch->branch() ;
3365          solver->solveFromHotStart() ;
3366          bool needHotStartUpdate=false;
3367          numberStrongDone++;
3368          numberStrongIterations += solver->getIterationCount();
3369          /*
3370            We now have an estimate of objective degradation that we can use for strong
3371            branching. If we're over the cutoff, the variable is monotone up.
3372            If we actually made it to optimality, check for a solution, and if we have
3373            a good one, call setBestSolution to process it. Note that this may reduce the
3374            cutoff, so we check again to see if we can declare this variable monotone.
3375          */
3376          if (solver->isProvenOptimal())
3377            iStatus=0; // optimal
3378          else if (solver->isIterationLimitReached()
3379                   &&!solver->isDualObjectiveLimitReached())
3380            iStatus=2; // unknown
3381          else
3382            iStatus=1; // infeasible
3383          if (iStatus!=2&&solver->getIterationCount()>
3384              realMaxHotIterations)
3385            numberUnfinished++;
3386          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3387          choice.numItersDown = solver->getIterationCount();
3388          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3389          // Update branching information if wanted
3390          CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3391          if (cbcobj) {
3392            CbcObject * object = cbcobj->object();
3393            assert (object) ;
3394            CbcObjectUpdateData update = object->createUpdateInformation(solver,this,cbcobj);
3395            update.objectNumber_ = choice.objectNumber;
3396            model->addUpdateInformation(update);
3397          } else {
3398            decision->updateInformation( solver,this);
3399          }
3400          if (!iStatus) {
3401            choice.finishedDown = true ;
3402            if (newObjectiveValue>=cutoff) {
3403              objectiveChange = 1.0e100; // say infeasible
3404              numberStrongInfeasible++;
3405            } else {
3406              // See if integer solution
3407              if (model->feasibleSolution(choice.numIntInfeasDown,
3408                                          choice.numObjInfeasDown)
3409                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
3410                if (auxiliaryInfo->solutionAddsCuts()) {
3411                  needHotStartUpdate=true;
3412                  solver->unmarkHotStart();
3413                }
3414                model->setBestSolution(CBC_STRONGSOL,
3415                                       newObjectiveValue,
3416                                       solver->getColSolution()) ;
3417                if (needHotStartUpdate) {
3418                  model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3419                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3420                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3421                  model->feasibleSolution(choice.numIntInfeasDown,
3422                                          choice.numObjInfeasDown);
3423                }
3424                model->setLastHeuristic(NULL);
3425                model->incrementUsed(solver->getColSolution());
3426                cutoff =model->getCutoff();
3427                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3428                  objectiveChange = 1.0e100 ;
3429                  numberStrongInfeasible++;
3430                }
3431              }
3432            }
3433          } else if (iStatus==1) {
3434            objectiveChange = 1.0e100 ;
3435            numberStrongInfeasible++;
3436          } else {
3437            // Can't say much as we did not finish
3438            choice.finishedDown = false ;
3439            numberUnfinished++;
3440          }
3441          choice.downMovement = objectiveChange ;
3442         
3443          // restore bounds
3444          for ( j=0;j<numberColumns;j++) {
3445            if (saveLower[j] != lower[j])
3446              solver->setColLower(j,saveLower[j]);
3447            if (saveUpper[j] != upper[j])
3448              solver->setColUpper(j,saveUpper[j]);
3449          }
3450          if(needHotStartUpdate) {
3451            needHotStartUpdate = false;
3452            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3453            //we may again have an integer feasible solution
3454            int numberIntegerInfeasibilities;
3455            int numberObjectInfeasibilities;
3456            if (model->feasibleSolution(
3457                                        numberIntegerInfeasibilities,
3458                                        numberObjectInfeasibilities)) {
3459#ifdef BONMIN
3460              //In this case node has become integer feasible, let us exit the loop
3461              std::cout<<"Node has become integer feasible"<<std::endl;
3462              numberUnsatisfied_ = 0;
3463              break;
3464#endif
3465              double objValue = solver->getObjValue();
3466              model->setBestSolution(CBC_STRONGSOL,
3467                                     objValue,
3468                                     solver->getColSolution()) ;
3469              model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3470              cutoff =model->getCutoff();
3471            }
3472            solver->markHotStart();
3473          }
3474#ifdef DO_ALL_AT_ROOT
3475          if (!numberNodes)
3476            printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3477                   choice.objectNumber,iStatus,newObjectiveValue,choice.numItersDown,
3478                   choice.downMovement,choice.finishedDown,choice.numIntInfeasDown,
3479                   choice.numObjInfeasDown);
3480#endif
3481         
3482          // repeat the whole exercise, forcing the variable up
3483          choice.possibleBranch->branch();
3484          solver->solveFromHotStart() ;
3485          numberStrongDone++;
3486          numberStrongIterations += solver->getIterationCount();
3487          /*
3488            We now have an estimate of objective degradation that we can use for strong
3489            branching. If we're over the cutoff, the variable is monotone up.
3490            If we actually made it to optimality, check for a solution, and if we have
3491            a good one, call setBestSolution to process it. Note that this may reduce the
3492            cutoff, so we check again to see if we can declare this variable monotone.
3493          */
3494          if (solver->isProvenOptimal())
3495            iStatus=0; // optimal
3496          else if (solver->isIterationLimitReached()
3497                   &&!solver->isDualObjectiveLimitReached())
3498            iStatus=2; // unknown
3499          else
3500            iStatus=1; // infeasible
3501          if (iStatus!=2&&solver->getIterationCount()>
3502              realMaxHotIterations)
3503            numberUnfinished++;
3504          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3505          choice.numItersUp = solver->getIterationCount();
3506          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3507          // Update branching information if wanted
3508          cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3509          if (cbcobj) {
3510            CbcObject * object = cbcobj->object();
3511            assert (object) ;
3512            CbcObjectUpdateData update = object->createUpdateInformation(solver,this,cbcobj);
3513            update.objectNumber_ = choice.objectNumber;
3514            model->addUpdateInformation(update);
3515          } else {
3516            decision->updateInformation( solver,this);
3517          }
3518          if (!iStatus) {
3519            choice.finishedUp = true ;
3520            if (newObjectiveValue>=cutoff) {
3521              objectiveChange = 1.0e100; // say infeasible
3522              numberStrongInfeasible++;
3523            } else {
3524              // See if integer solution
3525              if (model->feasibleSolution(choice.numIntInfeasUp,
3526                                          choice.numObjInfeasUp)
3527                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
3528#ifdef BONMIN
3529                std::cout<<"Node has become integer feasible"<<std::endl;
3530                numberUnsatisfied_ = 0;
3531                break;
3532#endif
3533                if (auxiliaryInfo->solutionAddsCuts()) {
3534                  needHotStartUpdate=true;
3535                  solver->unmarkHotStart();
3536                }
3537                model->setBestSolution(CBC_STRONGSOL,
3538                                       newObjectiveValue,
3539                                       solver->getColSolution()) ;
3540                if (needHotStartUpdate) {
3541                  model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3542                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3543                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3544                  model->feasibleSolution(choice.numIntInfeasDown,
3545                                          choice.numObjInfeasDown);
3546                }
3547                model->setLastHeuristic(NULL);
3548                model->incrementUsed(solver->getColSolution());
3549                cutoff =model->getCutoff();
3550                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3551                  objectiveChange = 1.0e100 ;
3552                  numberStrongInfeasible++;
3553                }
3554              }
3555            }
3556          } else if (iStatus==1) {
3557            objectiveChange = 1.0e100 ;
3558            numberStrongInfeasible++;
3559          } else {
3560            // Can't say much as we did not finish
3561            choice.finishedUp = false ;
3562            numberUnfinished++;
3563          }
3564          choice.upMovement = objectiveChange ;
3565         
3566          // restore bounds
3567          for ( j=0;j<numberColumns;j++) {
3568            if (saveLower[j] != lower[j])
3569              solver->setColLower(j,saveLower[j]);
3570            if (saveUpper[j] != upper[j])
3571              solver->setColUpper(j,saveUpper[j]);
3572          }
3573          if(needHotStartUpdate) {
3574            needHotStartUpdate = false;
3575            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3576            //we may again have an integer feasible solution
3577            int numberIntegerInfeasibilities;
3578            int numberObjectInfeasibilities;
3579            if (model->feasibleSolution(
3580                                        numberIntegerInfeasibilities,
3581                                        numberObjectInfeasibilities)) {
3582              double objValue = solver->getObjValue();
3583              model->setBestSolution(CBC_STRONGSOL,
3584                                     objValue,
3585                                     solver->getColSolution()) ;
3586              model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3587              cutoff =model->getCutoff();
3588            }
3589            solver->markHotStart();
3590          }
3591         
3592#ifdef DO_ALL_AT_ROOT
3593          if (!numberNodes)
3594            printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3595                   choice.objectNumber,iStatus,newObjectiveValue,choice.numItersUp,
3596                   choice.upMovement,choice.finishedUp,choice.numIntInfeasUp,
3597                   choice.numObjInfeasUp);
3598#endif
3599        }
3600       
3601        solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit2); 
3602        /*
3603          End of evaluation for this candidate variable. Possibilities are:
3604          * Both sides below cutoff; this variable is a candidate for branching.
3605          * Both sides infeasible or above the objective cutoff: no further action
3606          here. Break from the evaluation loop and assume the node will be purged
3607          by the caller.
3608          * One side below cutoff: Install the branch (i.e., fix the variable). Break
3609          from the evaluation loop and assume the node will be reoptimised by the
3610          caller.
3611        */
3612        // reset
3613        choice.possibleBranch->resetNumberBranchesLeft();
3614        if (choice.upMovement<1.0e100) {
3615          if(choice.downMovement<1.0e100) {
3616            // In case solution coming in was odd
3617            choice.upMovement = CoinMax(0.0,choice.upMovement);
3618            choice.downMovement = CoinMax(0.0,choice.downMovement);
3619#if ZERO_ONE==2
3620            // branch on 0-1 first (temp)
3621            if (fabs(choice.possibleBranch->value())<1.0) {
3622              choice.upMovement *= ZERO_FAKE;
3623              choice.downMovement *= ZERO_FAKE;
3624            }
3625#endif
3626            // feasible - see which best
3627            if (!canSkip) {
3628              if (model->messageHandler()->logLevel()>3) 
3629                printf("sort %g downest %g upest %g ",sort[iDo],downEstimate[iObject],
3630                       upEstimate[iObject]);
3631              model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
3632                << iObject << iColumn
3633                <<choice.downMovement<<choice.numIntInfeasDown
3634                <<choice.upMovement<<choice.numIntInfeasUp
3635                <<choice.possibleBranch->value()
3636                <<CoinMessageEol;
3637            }
3638            int betterWay;
3639            {
3640              CbcBranchingObject * branchObj =
3641                dynamic_cast <CbcBranchingObject *>(branch_) ;
3642              if (branch_)
3643                assert (branchObj);
3644              betterWay = decision->betterBranch(choice.possibleBranch,
3645                                                 branchObj,
3646                                                 choice.upMovement,
3647                                                 choice.numIntInfeasUp ,
3648                                                 choice.downMovement,
3649                                                 choice.numIntInfeasDown );
3650            }
3651            if (betterWay) {
3652              // C) create branching object
3653              if (choiceObject) {
3654                delete branch_;
3655                branch_ = choice.possibleBranch->clone();
3656              } else {
3657                delete branch_;
3658                branch_ = choice.possibleBranch;
3659                choice.possibleBranch=NULL;
3660              }
3661              {
3662                CbcBranchingObject * branchObj =
3663                  dynamic_cast <CbcBranchingObject *>(branch_) ;
3664                assert (branchObj);
3665                //branchObj->way(preferredWay);
3666                branchObj->way(betterWay);
3667              }
3668              bestChoice = choice.objectNumber;
3669              whichChoice = iDo;
3670              if (numberStrong<=1) {
3671                delete ws;
3672                ws=NULL;
3673                break;
3674              }
3675            } else {
3676              if (!choiceObject) {
3677                delete choice.possibleBranch;
3678                choice.possibleBranch=NULL;
3679              }
3680              if (iDo>=2*numberStrong) {
3681                delete ws;
3682                ws=NULL;
3683                break;
3684              }
3685              if (!dynamicObject||dynamicObject->numberTimesUp()>1) {
3686                if (iDo-whichChoice>=numberStrong) {
3687                  if (!choiceObject) {
3688                    delete choice.possibleBranch;
3689                    choice.possibleBranch=NULL;
3690                  }
3691                  break; // give up
3692                }
3693              } else {
3694                if (iDo-whichChoice>=2*numberStrong) {
3695                  delete ws;
3696                  ws=NULL;
3697                  if (!choiceObject) {
3698                    delete choice.possibleBranch;
3699                    choice.possibleBranch=NULL;
3700                  }
3701                  break; // give up
3702                }
3703              }
3704            }
3705          } else {
3706            // up feasible, down infeasible
3707            anyAction=-1;
3708            worstFeasible = CoinMax(worstFeasible,choice.upMovement);
3709            model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
3710              << iObject << iColumn
3711              <<choice.downMovement<<choice.numIntInfeasDown
3712              <<choice.upMovement<<choice.numIntInfeasUp
3713              <<choice.possibleBranch->value()
3714              <<CoinMessageEol;
3715            //printf("Down infeasible for choice %d sequence %d\n",i,
3716            // model->object(choice.objectNumber)->columnNumber());
3717            choice.fix=1;
3718            numberToFix++;
3719            choice.possibleBranch->fix(solver,saveLower,saveUpper,1);
3720            if (!choiceObject) {
3721              choice.possibleBranch=NULL;
3722            } else {
3723              //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3724              choice.possibleBranch=choiceObject;
3725            }
3726            assert(doneHotStart);
3727            solver->unmarkHotStart();
3728            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3729            solver->markHotStart();
3730            // may be infeasible (if other way stopped on iterations)
3731            if (!solver->isProvenOptimal()) {
3732              // neither side feasible
3733              anyAction=-2;
3734              if (!choiceObject) {
3735                delete choice.possibleBranch;
3736                choice.possibleBranch=NULL;
3737              }
3738              //printf("Both infeasible for choice %d sequence %d\n",i,
3739              // model->object(choice.objectNumber)->columnNumber());
3740              delete ws;
3741              ws=NULL;
3742              break;
3743            }
3744          }
3745        } else {
3746          if(choice.downMovement<1.0e100) {
3747            // down feasible, up infeasible
3748            anyAction=-1;
3749            worstFeasible = CoinMax(worstFeasible,choice.downMovement);
3750            model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
3751              << iObject << iColumn
3752              <<choice.downMovement<<choice.numIntInfeasDown
3753              <<choice.upMovement<<choice.numIntInfeasUp
3754              <<choice.possibleBranch->value()
3755              <<CoinMessageEol;
3756            choice.fix=-1;
3757            numberToFix++;
3758            choice.possibleBranch->fix(solver,saveLower,saveUpper,-1);
3759            if (!choiceObject) {
3760              choice.possibleBranch=NULL;
3761            } else {
3762              //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3763              choice.possibleBranch=choiceObject;
3764            }
3765            assert(doneHotStart);
3766            solver->unmarkHotStart();
3767            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3768            solver->markHotStart();
3769            // may be infeasible (if other way stopped on iterations)
3770            if (!solver->isProvenOptimal()) {
3771              // neither side feasible
3772              anyAction=-2;
3773              if (!choiceObject) {
3774                delete choice.possibleBranch;
3775                choice.possibleBranch=NULL;
3776              }
3777              delete ws;
3778              ws=NULL;
3779              break;
3780            }
3781          } else {
3782            // neither side feasible
3783            anyAction=-2;
3784            if (!choiceObject) {
3785              delete choice.possibleBranch;
3786              choice.possibleBranch=NULL;
3787            }
3788            delete ws;
3789            ws=NULL;
3790            break;
3791          }
3792        }
3793        // Check max time
3794        hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
3795                       model->getDblParam(CbcModel::CbcMaximumSeconds));
3796        if (hitMaxTime) {
3797          // make sure rest are fast
3798          doQuickly=true;
3799          for ( int jDo=iDo+1;jDo<numberToDo;jDo++) {
3800            int iObject = whichObject[iDo];
3801            OsiObject * object = model->modifiableObject(iObject);
3802            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3803              dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3804            if (dynamicObject) 
3805              dynamicObject->setNumberBeforeTrust(0);
3806          }
3807          numberTest=0;
3808        }
3809        if (!choiceObject) {
3810          delete choice.possibleBranch;
3811        }
3812      }
3813      if (model->messageHandler()->logLevel()>3) { 
3814        if (anyAction==-2) {
3815          printf("infeasible\n");
3816        } else if(anyAction==-1) {
3817          printf("%d fixed AND choosing %d iDo %d iChosenWhen %d numberToDo %d\n",numberToFix,bestChoice,
3818                   iDo,whichChoice,numberToDo);
3819        } else {
3820          int iObject = whichObject[whichChoice];
3821          OsiObject * object = model->modifiableObject(iObject);
3822          CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3823            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3824          if (dynamicObject) {
3825            int iColumn = dynamicObject->columnNumber();
3826            printf("choosing %d (column %d) iChosenWhen %d numberToDo %d\n",bestChoice,
3827                   iColumn,whichChoice,numberToDo);
3828          }
3829        }
3830      }
3831      if (doneHotStart) {
3832        // Delete the snapshot
3833        solver->unmarkHotStart();
3834        // back to normal
3835        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3836        // restore basis
3837        solver->setWarmStart(ws);
3838      }
3839      solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
3840      // Unless infeasible we will carry on
3841      // But we could fix anyway
3842      if (numberToFix&&!hitMaxTime) {
3843        if (anyAction!=-2) {
3844          // apply and take off
3845          bool feasible=true;
3846          // can do quick optimality check
3847          int easy=2;
3848          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
3849          model->resolve(NULL,11,saveSolution,saveLower,saveUpper) ;
3850          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3851          feasible = solver->isProvenOptimal();
3852          if (feasible) {
3853            anyAction=0;
3854            // See if candidate still possible
3855            if (branch_) {
3856              const OsiObject * object = model->object(bestChoice);
3857              int preferredWay;
3858              double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
3859              if (!infeasibility) {
3860                // take out
3861                delete branch_;
3862                branch_=NULL;
3863              } else {
3864                CbcBranchingObject * branchObj =
3865                  dynamic_cast <CbcBranchingObject *>(branch_) ;
3866                assert (branchObj);
3867                branchObj->way(preferredWay);
3868              }
3869            }
3870          } else {
3871            anyAction=-2;
3872            finished=true;
3873          }
3874        }
3875      }
3876      // If  fixed then round again
3877      if (!branch_&&anyAction!=-2) {
3878        finished=false;
3879      }
3880      delete ws;
3881    }
3882  }
3883  // update number of strong iterations etc
3884  model->incrementStrongInfo(numberStrongDone,numberStrongIterations,
3885                             anyAction==-2 ? 0:numberToFix,anyAction==-2);
3886  if (model->searchStrategy()==-1) {
3887#ifndef COIN_DEVELOP
3888    if (solver->messageHandler()->logLevel()>1)
3889#endif
3890      printf("%d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3891             numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3892             numberNotTrusted);
3893    // decide what to do
3894    int strategy=1;
3895    if (((numberUnfinished*4>numberStrongDone&&
3896          numberStrongInfeasible*40<numberStrongDone)||
3897         numberStrongInfeasible<0)&&model->numberStrong()<10&&model->numberBeforeTrust()<=20&&model->numberObjects()>CoinMax(1000,solver->getNumRows())) {
3898      strategy=2;
3899#ifdef COIN_DEVELOP
3900      //if (model->logLevel()>1)
3901      printf("going to strategy 2\n");
3902#endif
3903      // Weaken
3904      model->setNumberStrong(2);
3905      model->setNumberBeforeTrust(1);
3906      model->synchronizeNumberBeforeTrust();
3907    }
3908    if (numberNodes)
3909      strategy=1;  // should only happen after hot start
3910    model->setSearchStrategy(strategy);
3911  } else if (numberStrongDone) {
3912    //printf("%d strongB, %d iters, %d inf, %d not finished, %d not trusted\n",
3913    //   numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3914    //   numberNotTrusted);
3915  }
3916  if (model->searchStrategy()==1&&numberNodes>500&&numberNodes<-510) {
3917#ifndef COIN_DEVELOP
3918    if (solver->messageHandler()->logLevel()>1)
3919#endif
3920      printf("after %d nodes - %d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3921             numberNodes,numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3922             numberNotTrusted);
3923    // decide what to do
3924    if (numberUnfinished*10>numberStrongDone+1||
3925        !numberStrongInfeasible) {
3926      printf("going to strategy 2\n");
3927      // Weaken
3928      model->setNumberStrong(2);
3929      model->setNumberBeforeTrust(1);
3930      model->synchronizeNumberBeforeTrust();
3931      model->setSearchStrategy(2);
3932    }
3933  }
3934  // Set guessed solution value
3935  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
3936 
3937  /*
3938    Cleanup, then we're finished
3939  */
3940  if (!model->branchingMethod())
3941    delete decision;
3942 
3943  delete choiceObject;
3944  delete [] sort;
3945  delete [] whichObject;
3946#ifdef RANGING
3947  delete [] objectMark;
3948#endif
3949  delete [] saveLower;
3950  delete [] saveUpper;
3951  delete [] upEstimate;
3952  delete [] downEstimate;
3953# ifdef COIN_HAS_CLP
3954  if (osiclp) {
3955    osiclp->setSpecialOptions(saveClpOptions);
3956  }
3957# endif
3958  // restore solution
3959  solver->setColSolution(saveSolution);
3960  model->reserveCurrentSolution(saveSolution);
3961  delete [] saveSolution;
3962  model->setStateOfSearch(saveStateOfSearch);
3963  model->setLogLevel(saveLogLevel);
3964  // delete extra regions
3965  if (usefulInfo.usefulRegion_) {
3966    delete [] usefulInfo.usefulRegion_;
3967    delete [] usefulInfo.indexRegion_;
3968    delete [] usefulInfo.pi_;
3969    usefulInfo.usefulRegion_ = NULL;
3970    usefulInfo.indexRegion_ = NULL;
3971    usefulInfo.pi_ = NULL;
3972  }
3973  return anyAction;
3974}
3975int CbcNode::analyze (CbcModel *model, double * results)
3976{
3977  int i;
3978  int numberIterationsAllowed = model->numberAnalyzeIterations();
3979  OsiSolverInterface * solver = model->solver();
3980  objectiveValue_ = solver->getObjSense()*solver->getObjValue();
3981  double cutoff =model->getCutoff();
3982  const double * lower = solver->getColLower();
3983  const double * upper = solver->getColUpper();
3984  const double * dj = solver->getReducedCost();
3985  int numberObjects = model->numberObjects();
3986  int numberColumns = model->getNumCols();
3987  // Initialize arrays
3988  int numberIntegers = model->numberIntegers();
3989  int * back = new int[numberColumns];
3990  const int * integerVariable = model->integerVariable();
3991  for (i=0;i<numberColumns;i++) 
3992    back[i]=-1;
3993  // What results is
3994  double * newLower = results;
3995  double * objLower = newLower+numberIntegers;
3996  double * newUpper = objLower+numberIntegers;
3997  double * objUpper = newUpper+numberIntegers;
3998  for (i=0;i<numberIntegers;i++) {
3999    int iColumn = integerVariable[i];
4000    back[iColumn]=i;
4001    newLower[i]=0.0;
4002    objLower[i]=-COIN_DBL_MAX;
4003    newUpper[i]=0.0;
4004    objUpper[i]=-COIN_DBL_MAX;
4005  }
4006  double * saveUpper = new double[numberColumns];
4007  double * saveLower = new double[numberColumns];
4008  int anyAction=0;
4009  // Save solution in case heuristics need good solution later
4010 
4011  double * saveSolution = new double[numberColumns];
4012  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
4013  model->reserveCurrentSolution(saveSolution);
4014  for (i=0;i<numberColumns;i++) {
4015    saveLower[i] = lower[i];
4016    saveUpper[i] = upper[i];
4017  }
4018  // Get arrays to sort
4019  double * sort = new double[numberObjects];
4020  int * whichObject = new int[numberObjects];
4021  int numberToFix=0;
4022  int numberToDo=0;
4023  double integerTolerance = 
4024    model->getDblParam(CbcModel::CbcIntegerTolerance);
4025  // point to useful information
4026  OsiBranchingInformation usefulInfo = model->usefulInformation();
4027  // and modify
4028  usefulInfo.depth_=depth_;
4029 
4030  // compute current state
4031  int numberObjectInfeasibilities; // just odd ones
4032  int numberIntegerInfeasibilities;
4033  model->feasibleSolution(
4034                          numberIntegerInfeasibilities,
4035                          numberObjectInfeasibilities);
4036# ifdef COIN_HAS_CLP
4037  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
4038  int saveClpOptions=0;
4039  bool fastIterations = (model->specialOptions()&8)!=0;
4040  if (osiclp&&fastIterations) {
4041    // for faster hot start
4042    saveClpOptions = osiclp->specialOptions();
4043    osiclp->setSpecialOptions(saveClpOptions|8192);
4044  }
4045# else
4046  bool fastIterations = false ;
4047# endif
4048  /*
4049    Scan for branching objects that indicate infeasibility. Choose candidates
4050    using priority as the first criteria, then integer infeasibility.
4051   
4052    The algorithm is to fill the array with a set of good candidates (by
4053    infeasibility) with priority bestPriority.  Finding a candidate with
4054    priority better (less) than bestPriority flushes the choice array. (This
4055    serves as initialization when the first candidate is found.)
4056   
4057  */
4058  numberToDo=0;
4059  for (i=0;i<numberObjects;i++) {
4060    OsiObject * object = model->modifiableObject(i);
4061    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4062      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4063    if(!dynamicObject)
4064      continue;
4065    int preferredWay;
4066    double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
4067    int iColumn = dynamicObject->columnNumber();
4068    if (saveUpper[iColumn]==saveLower[iColumn])
4069      continue;
4070    if (infeasibility)
4071      sort[numberToDo]=-1.0e10-infeasibility;
4072    else
4073      sort[numberToDo]=-fabs(dj[iColumn]);
4074    whichObject[numberToDo++]=i;
4075  }
4076  // Save basis
4077  CoinWarmStart * ws = solver->getWarmStart();
4078  int saveLimit;
4079  solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
4080  int targetIterations = CoinMax(500,numberIterationsAllowed/numberObjects);
4081  if (saveLimit<targetIterations)
4082    solver->setIntParam(OsiMaxNumIterationHotStart,targetIterations); 
4083  // Mark hot start
4084  solver->markHotStart();
4085  // Sort
4086  CoinSort_2(sort,sort+numberToDo,whichObject);
4087  double * currentSolution = model->currentSolution();
4088  double objMin = 1.0e50;
4089  double objMax = -1.0e50;
4090  bool needResolve=false;
4091  int iDo;
4092  for (iDo=0;iDo<numberToDo;iDo++) {
4093    CbcStrongInfo choice;
4094    int iObject = whichObject[iDo];
4095    OsiObject * object = model->modifiableObject(iObject);
4096    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4097      static_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4098    int iColumn = dynamicObject->columnNumber();
4099    int preferredWay;
4100    object->infeasibility(&usefulInfo,preferredWay);
4101    double value = currentSolution[iColumn];
4102    double nearest = floor(value+0.5);
4103    double lowerValue = floor(value);
4104    bool satisfied=false;
4105    if (fabs(value-nearest)<=integerTolerance||value<saveLower[iColumn]||value>saveUpper[iColumn]) {
4106      satisfied=true;
4107      double newValue;
4108      if (nearest<saveUpper[iColumn]) {
4109        newValue = nearest + 1.0001*integerTolerance;
4110        lowerValue = nearest;
4111      } else {
4112        newValue = nearest - 1.0001*integerTolerance;
4113        lowerValue = nearest-1;
4114      }
4115      currentSolution[iColumn]=newValue;
4116    }
4117    double upperValue = lowerValue+1.0;
4118    //CbcSimpleInteger * obj =
4119    //dynamic_cast <CbcSimpleInteger *>(object) ;
4120    //if (obj) {
4121    //choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
4122    //} else {
4123      CbcObject * obj =
4124        dynamic_cast <CbcObject *>(object) ;
4125      assert (obj);
4126      choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
4127      //}
4128    currentSolution[iColumn]=value;
4129    // Save which object it was
4130    choice.objectNumber=iObject;
4131    choice.numIntInfeasUp = numberUnsatisfied_;
4132    choice.numIntInfeasDown = numberUnsatisfied_;
4133    choice.downMovement = 0.0;
4134    choice.upMovement = 0.0;
4135    choice.numItersDown = 0;
4136    choice.numItersUp = 0;
4137    choice.fix=0; // say not fixed
4138    double objectiveChange ;
4139    double newObjectiveValue=1.0e100;
4140    int j;
4141    // status is 0 finished, 1 infeasible and other
4142    int iStatus;
4143    /*
4144      Try the down direction first. (Specify the initial branching alternative as
4145      down with a call to way(-1). Each subsequent call to branch() performs the
4146      specified branch and advances the branch object state to the next branch
4147      alternative.)
4148    */
4149    choice.possibleBranch->way(-1) ;
4150    choice.possibleBranch->branch() ;
4151    if (fabs(value-lowerValue)>integerTolerance) {
4152      solver->solveFromHotStart() ;
4153      /*
4154        We now have an estimate of objective degradation that we can use for strong
4155        branching. If we're over the cutoff, the variable is monotone up.
4156        If we actually made it to optimality, check for a solution, and if we have
4157        a good one, call setBestSolution to process it. Note that this may reduce the
4158        cutoff, so we check again to see if we can declare this variable monotone.
4159      */
4160      if (solver->isProvenOptimal())
4161        iStatus=0; // optimal
4162      else if (solver->isIterationLimitReached()
4163               &&!solver->isDualObjectiveLimitReached())
4164        iStatus=2; // unknown
4165      else
4166        iStatus=1; // infeasible
4167      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
4168      choice.numItersDown = solver->getIterationCount();
4169      numberIterationsAllowed -= choice.numItersDown;
4170      objectiveChange = newObjectiveValue  - objectiveValue_;
4171      if (!iStatus) {
4172        choice.finishedDown = true ;
4173        if (newObjectiveValue>=cutoff) {
4174          objectiveChange = 1.0e100; // say infeasible
4175        } else {
4176          // See if integer solution
4177          if (model->feasibleSolution(choice.numIntInfeasDown,
4178                                      choice.numObjInfeasDown)
4179              &&model->problemFeasibility()->feasible(model,-1)>=0) {
4180            model->setBestSolution(CBC_STRONGSOL,
4181                                   newObjectiveValue,
4182                                   solver->getColSolution()) ;
4183            model->setLastHeuristic(NULL);
4184            model->incrementUsed(solver->getColSolution());
4185            cutoff =model->getCutoff();
4186            if (newObjectiveValue >= cutoff)    //  *new* cutoff
4187              objectiveChange = 1.0e100 ;
4188          }
4189        }
4190      } else if (iStatus==1) {
4191        objectiveChange = 1.0e100 ;
4192      } else {
4193        // Can't say much as we did not finish
4194        choice.finishedDown = false ;
4195      }
4196      choice.downMovement = objectiveChange ;
4197    }
4198    // restore bounds
4199    for ( j=0;j<numberColumns;j++) {
4200      if (saveLower[j] != lower[j])
4201        solver->setColLower(j,saveLower[j]);
4202      if (saveUpper[j] != upper[j])
4203        solver->setColUpper(j,saveUpper[j]);
4204    }
4205    // repeat the whole exercise, forcing the variable up
4206    choice.possibleBranch->branch();
4207    if (fabs(value-upperValue)>integerTolerance) {
4208      solver->solveFromHotStart() ;
4209      /*
4210        We now have an estimate of objective degradation that we can use for strong
4211        branching. If we're over the cutoff, the variable is monotone up.
4212        If we actually made it to optimality, check for a solution, and if we have
4213        a good one, call setBestSolution to process it. Note that this may reduce the
4214        cutoff, so we check again to see if we can declare this variable monotone.
4215      */
4216      if (solver->isProvenOptimal())
4217        iStatus=0; // optimal
4218      else if (solver->isIterationLimitReached()
4219               &&!solver->isDualObjectiveLimitReached())
4220        iStatus=2; // unknown
4221      else
4222        iStatus=1; // infeasible
4223      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
4224      choice.numItersUp = solver->getIterationCount();
4225      numberIterationsAllowed -= choice.numItersUp;
4226      objectiveChange = newObjectiveValue  - objectiveValue_;
4227      if (!iStatus) {
4228        choice.finishedUp = true ;
4229        if (newObjectiveValue>=cutoff) {
4230          objectiveChange = 1.0e100; // say infeasible
4231        } else {
4232          // See if integer solution
4233          if (model->feasibleSolution(choice.numIntInfeasUp,
4234                                      choice.numObjInfeasUp)
4235              &&model->problemFeasibility()->feasible(model,-1)>=0) {
4236            model->setBestSolution(CBC_STRONGSOL,
4237                                   newObjectiveValue,
4238                                   solver->getColSolution()) ;
4239            model->setLastHeuristic(NULL);
4240            model->incrementUsed(solver->getColSolution());
4241            cutoff =model->getCutoff();
4242            if (newObjectiveValue >= cutoff)    //  *new* cutoff
4243              objectiveChange = 1.0e100 ;
4244          }
4245        }
4246      } else if (iStatus==1) {
4247        objectiveChange = 1.0e100 ;
4248      } else {
4249        // Can't say much as we did not finish
4250        choice.finishedUp = false ;
4251      }
4252      choice.upMovement = objectiveChange ;
4253     
4254      // restore bounds
4255      for ( j=0;j<numberColumns;j++) {
4256        if (saveLower[j] != lower[j])
4257          solver->setColLower(j,saveLower[j]);
4258        if (saveUpper[j] != upper[j])
4259          solver->setColUpper(j,saveUpper[j]);
4260      }
4261    }
4262    // If objective goes above certain amount we can set bound
4263    int jInt = back[iColumn];
4264    newLower[jInt]=upperValue;
4265    if (choice.finishedDown)
4266      objLower[jInt]=choice.downMovement+objectiveValue_;
4267    else
4268      objLower[jInt]=objectiveValue_;
4269    newUpper[jInt]=lowerValue;
4270    if (choice.finishedUp)
4271      objUpper[jInt]=choice.upMovement+objectiveValue_;
4272    else
4273      objUpper[jInt]=objectiveValue_;
4274    objMin = CoinMin(CoinMin(objLower[jInt],objUpper[jInt]),objMin);
4275    /*
4276      End of evaluation for this candidate variable. Possibilities are:
4277      * Both sides below cutoff; this variable is a candidate for branching.
4278      * Both sides infeasible or above the objective cutoff: no further action
4279      here. Break from the evaluation loop and assume the node will be purged
4280      by the caller.
4281      * One side below cutoff: Install the branch (i.e., fix the variable). Break
4282      from the evaluation loop and assume the node will be reoptimised by the
4283      caller.
4284    */
4285    if (choice.upMovement<1.0e100) {
4286      if(choice.downMovement<1.0e100) {
4287        objMax = CoinMax(CoinMax(objLower[jInt],objUpper[jInt]),objMax);
4288        // In case solution coming in was odd
4289        choice.upMovement = CoinMax(0.0,choice.upMovement);
4290        choice.downMovement = CoinMax(0.0,choice.downMovement);
4291        // feasible -
4292        model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
4293          << iObject << iColumn
4294          <<choice.downMovement<<choice.numIntInfeasDown
4295          <<choice.upMovement<<choice.numIntInfeasUp
4296          <<value
4297          <<CoinMessageEol;
4298      } else {
4299        // up feasible, down infeasible
4300        anyAction=-1;
4301        if (!satisfied)
4302          needResolve=true;
4303        choice.fix=1;
4304        numberToFix++;
4305        saveLower[iColumn]=upperValue;
4306        solver->setColLower(iColumn,upperValue);
4307      }
4308    } else {
4309      if(choice.downMovement<1.0e100) {
4310        // down feasible, up infeasible
4311        anyAction=-1;
4312        if (!satisfied)
4313          needResolve=true;
4314        choice.fix=-1;
4315        numberToFix++;
4316        saveUpper[iColumn]=lowerValue;
4317        solver->setColUpper(iColumn,lowerValue);
4318      } else {
4319        // neither side feasible
4320        anyAction=-2;
4321        printf("Both infeasible for choice %d sequence %d\n",i,
4322               model->object(choice.objectNumber)->columnNumber());
4323        delete ws;
4324        ws=NULL;
4325        //solver->writeMps("bad");
4326        numberToFix=-1;
4327        delete choice.possibleBranch;
4328        choice.possibleBranch=NULL;
4329        break;
4330      }
4331    }
4332    delete choice.possibleBranch;
4333    if (numberIterationsAllowed<=0)
4334      break;
4335    //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
4336    //     choice.downMovement,choice.upMovement,value);
4337  }
4338  printf("Best possible solution %g, can fix more if solution of %g found - looked at %d variables in %d iterations\n",
4339         objMin,objMax,iDo,model->numberAnalyzeIterations()-numberIterationsAllowed);
4340  model->setNumberAnalyzeIterations(numberIterationsAllowed);
4341  // Delete the snapshot
4342  solver->unmarkHotStart();
4343  // back to normal
4344  solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
4345  solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
4346  // restore basis
4347  solver->setWarmStart(ws);
4348  delete ws;
4349 
4350  delete [] sort;
4351  delete [] whichObject;
4352  delete [] saveLower;
4353  delete [] saveUpper;
4354  delete [] back;
4355  // restore solution
4356  solver->setColSolution(saveSolution);
4357# ifdef COIN_HAS_CLP
4358  if (osiclp) 
4359    osiclp->setSpecialOptions(saveClpOptions);
4360# endif
4361  model->reserveCurrentSolution(saveSolution);
4362  delete [] saveSolution;
4363  if (needResolve)
4364    solver->resolve();
4365  return numberToFix;
4366}
4367
4368
4369CbcNode::CbcNode(const CbcNode & rhs)
4370  : CoinTreeNode(rhs) 
4371{ 
4372#ifdef CHECK_NODE
4373  printf("CbcNode %x Constructor from rhs %x\n",this,&rhs);
4374#endif
4375  if (rhs.nodeInfo_)
4376    nodeInfo_ = rhs.nodeInfo_->clone();
4377  else
4378    nodeInfo_=NULL;
4379  objectiveValue_=rhs.objectiveValue_;
4380  guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4381  sumInfeasibilities_ = rhs.sumInfeasibilities_;
4382  if (rhs.branch_)
4383    branch_=rhs.branch_->clone();
4384  else
4385    branch_=NULL;
4386  depth_ = rhs.depth_;
4387  numberUnsatisfied_ = rhs.numberUnsatisfied_;
4388  nodeNumber_ = rhs.nodeNumber_;
4389  state_ = rhs.state_;
4390  if (nodeInfo_)
4391    assert ((state_&2)!=0);
4392  else
4393    assert ((state_&2)==0);
4394}
4395
4396CbcNode &
4397CbcNode::operator=(const CbcNode & rhs)
4398{
4399  if (this != &rhs) {
4400    delete nodeInfo_;
4401    if (rhs.nodeInfo_)
4402      nodeInfo_ = rhs.nodeInfo_->clone();
4403    else
4404      nodeInfo_ = NULL;
4405    objectiveValue_=rhs.objectiveValue_;
4406    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4407    sumInfeasibilities_ = rhs.sumInfeasibilities_;
4408    if (rhs.branch_)
4409      branch_=rhs.branch_->clone();
4410    else
4411      branch_=NULL,
4412        depth_ = rhs.depth_;
4413    numberUnsatisfied_ = rhs.numberUnsatisfied_;
4414    nodeNumber_ = rhs.nodeNumber_;
4415    state_ = rhs.state_;
4416    if (nodeInfo_)
4417      assert ((state_&2)!=0);
4418    else
4419      assert ((state_&2)==0);
4420  }
4421  return *this;
4422}
4423CbcNode::~CbcNode ()
4424{
4425#ifdef CHECK_NODE
4426  if (nodeInfo_) {
4427    printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
4428           this,nodeInfo_,nodeInfo_->numberPointingToThis());
4429    //assert(nodeInfo_->numberPointingToThis()>=0);
4430  } else {
4431    printf("CbcNode %x Destructor nodeInfo %x (?)\n",
4432           this,nodeInfo_);
4433  }
4434#endif
4435  if (nodeInfo_) {
4436    // was if (nodeInfo_&&(state_&2)!=0) {
4437    nodeInfo_->nullOwner();
4438    int numberToDelete=nodeInfo_->numberBranchesLeft();
4439    //    CbcNodeInfo * parent = nodeInfo_->parent();
4440    //assert (nodeInfo_->numberPointingToThis()>0);
4441    if (nodeInfo_->decrement(numberToDelete)==0||(state_&2)==0) {
4442      if ((state_&2)==0) 
4443        nodeInfo_->nullParent();
4444      delete nodeInfo_;
4445    } else {
4446      //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,nodeInfo_->parent());
4447      // anyway decrement parent
4448      //if (parent)
4449      ///parent->decrement(1);
4450    }
4451  }
4452  delete branch_;
4453}
4454// Decrement  active cut counts
4455void 
4456CbcNode::decrementCuts(int change)
4457{
4458  if (nodeInfo_)
4459    assert ((state_&2)!=0);
4460  else
4461    assert ((state_&2)==0);
4462  if(nodeInfo_) {
4463    nodeInfo_->decrementCuts(change);
4464  }
4465}
4466void 
4467CbcNode::decrementParentCuts(CbcModel * model, int change)
4468{
4469  if (nodeInfo_)
4470    assert ((state_&2)!=0);
4471  else
4472    assert ((state_&2)==0);
4473  if(nodeInfo_) {
4474    nodeInfo_->decrementParentCuts(model, change);
4475  }
4476}
4477
4478/*
4479  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
4480  in the attached nodeInfo_.
4481*/
4482void
4483CbcNode::initializeInfo()
4484{
4485  assert(nodeInfo_ && branch_) ;
4486  nodeInfo_->initializeInfo(branch_->numberBranches());
4487  assert ((state_&2)!=0);
4488  assert (nodeInfo_->numberBranchesLeft()==
4489          branch_->numberBranchesLeft());
4490}
4491// Nulls out node info
4492void 
4493CbcNode::nullNodeInfo()
4494{
4495  nodeInfo_=NULL;
4496  // say not active
4497  state_ &= ~2;
4498}
4499
4500int
4501CbcNode::branch(OsiSolverInterface * solver)
4502{
4503  double changeInGuessed;
4504  assert (nodeInfo_->numberBranchesLeft()==
4505          branch_->numberBranchesLeft());
4506  if (!solver)
4507    changeInGuessed=branch_->branch();
4508  else
4509    changeInGuessed=branch_->branch(solver);
4510  guessedObjectiveValue_+= changeInGuessed;
4511  //#define PRINTIT
4512#ifdef PRINTIT
4513  int numberLeft = nodeInfo_->numberBranchesLeft();
4514  CbcNodeInfo * parent = nodeInfo_->parent();
4515  int parentNodeNumber = -1;
4516  CbcBranchingObject * object1 = 
4517    dynamic_cast<CbcBranchingObject *>(branch_) ;
4518  //OsiObject * object = object1->
4519  //int sequence = object->columnNumber);
4520  int id=-1;
4521  double value=0.0;
4522  if (object1) {
4523    id = object1->variable();
4524    value = object1->value();
4525  }
4526  printf("id %d value %g objvalue %g\n",id,value,objectiveValue_);
4527  if (parent)
4528    parentNodeNumber = parent->nodeNumber();
4529  printf("Node number %d, %s, way %d, depth %d, parent node number %d\n",
4530         nodeInfo_->nodeNumber(),(numberLeft==2) ? "leftBranch" : "rightBranch",
4531         way(),depth_,parentNodeNumber);
4532  assert (parentNodeNumber!=nodeInfo_->nodeNumber());
4533#endif
4534  return nodeInfo_->branchedOn();
4535}
4536/* Active arm of the attached OsiBranchingObject.
4537   
4538   In the simplest instance, coded -1 for the down arm of the branch, +1 for
4539   the up arm. But see OsiBranchingObject::way()
4540   Use nodeInfo--.numberBranchesLeft_ to see how active
4541*/
4542int 
4543CbcNode::way() const
4544{
4545  if (branch_) {
4546    CbcBranchingObject * obj =
4547      dynamic_cast <CbcBranchingObject *>(branch_) ;
4548    if (obj) {
4549      return obj->way();
4550    } else {
4551      OsiTwoWayBranchingObject * obj2 =
4552        dynamic_cast <OsiTwoWayBranchingObject *>(branch_) ;
4553      assert (obj2);
4554      return obj2->way();
4555    }
4556  } else {
4557    return 0;
4558  }
4559}
4560/* Create a branching object for the node
4561   
4562   The routine scans the object list of the model and selects a set of
4563   unsatisfied objects as candidates for branching. The candidates are
4564   evaluated, and an appropriate branch object is installed.
4565   
4566   The numberPassesLeft is decremented to stop fixing one variable each time
4567   and going on and on (e.g. for stock cutting, air crew scheduling)
4568   
4569   If evaluation determines that an object is monotone or infeasible,
4570   the routine returns immediately. In the case of a monotone object,
4571   the branch object has already been called to modify the model.
4572   
4573   Return value:
4574   <ul>
4575   <li>  0: A branching object has been installed
4576   <li> -1: A monotone object was discovered
4577   <li> -2: An infeasible object was discovered
4578   </ul>
4579   Branch state:
4580   <ul>
4581   <li> -1: start
4582   <li> -1: A monotone object was discovered
4583   <li> -2: An infeasible object was discovered
4584   </ul>
4585*/
4586int 
4587CbcNode::chooseOsiBranch (CbcModel * model,
4588                          CbcNode * lastNode,
4589                          OsiBranchingInformation * usefulInfo,
4590                          int branchState)
4591{
4592  int returnStatus=0;
4593  if (lastNode)
4594    depth_ = lastNode->depth_+1;
4595  else
4596    depth_ = 0;
4597  OsiSolverInterface * solver = model->solver();
4598  objectiveValue_ = solver->getObjValue()*solver->getObjSense();
4599  usefulInfo->objectiveValue_ = objectiveValue_;
4600  usefulInfo->depth_ = depth_;
4601  const double * saveInfoSol = usefulInfo->solution_;
4602  double * saveSolution = new double[solver->getNumCols()];
4603  memcpy(saveSolution,solver->getColSolution(),solver->getNumCols()*sizeof(double));
4604  usefulInfo->solution_ = saveSolution;
4605  OsiChooseVariable * choose = model->branchingMethod()->chooseMethod();
4606  int numberUnsatisfied=-1;
4607  if (branchState<0) {
4608    // initialize
4609    // initialize sum of "infeasibilities"
4610    sumInfeasibilities_ = 0.0;
4611    numberUnsatisfied = choose->setupList(usefulInfo,true);
4612    numberUnsatisfied_ = numberUnsatisfied;
4613    branchState=0;
4614    if (numberUnsatisfied_<0) {
4615      // infeasible
4616      delete [] saveSolution;
4617      return -2;
4618    }
4619  }
4620  // unset best
4621  int best=-1;
4622  choose->setBestObjectIndex(-1);
4623  if (numberUnsatisfied) {
4624    if (branchState>0||!choose->numberOnList()) {
4625      // we need to return at once - don't do strong branching or anything
4626      if (choose->numberOnList()||!choose->numberStrong()) {
4627        best = choose->candidates()[0];
4628        choose->setBestObjectIndex(best);
4629      } else {
4630        // nothing on list - need to try again - keep any solution
4631        numberUnsatisfied = choose->setupList(usefulInfo, false);
4632        numberUnsatisfied_ = numberUnsatisfied;
4633        if (numberUnsatisfied) {
4634          best = choose->candidates()[0];
4635          choose->setBestObjectIndex(best);
4636        }
4637      }
4638    } else {
4639      // carry on with strong branching or whatever
4640      int returnCode = choose->chooseVariable(solver, usefulInfo,true);
4641      // update number of strong iterations etc
4642      model->incrementStrongInfo(choose->numberStrongDone(),choose->numberStrongIterations(),
4643                                 returnCode==-1 ? 0:choose->numberStrongFixed(),returnCode==-1);
4644      if (returnCode>1) {
4645        // has fixed some
4646        returnStatus=-1;
4647      } else if (returnCode==-1) {
4648        // infeasible
4649        returnStatus=-2;
4650      } else if (returnCode==0) {
4651        // normal
4652        returnStatus=0;
4653        numberUnsatisfied=1;
4654      } else {
4655        // ones on list satisfied - double check
4656        numberUnsatisfied = choose->setupList(usefulInfo, false);
4657        numberUnsatisfied_ = numberUnsatisfied;
4658        if (numberUnsatisfied) {
4659          best = choose->candidates()[0];
4660          choose->setBestObjectIndex(best);
4661        }
4662      }
4663    }
4664  } 
4665  delete branch_;
4666  branch_ = NULL;
4667  guessedObjectiveValue_ = COIN_DBL_MAX;//objectiveValue_; // for now
4668  if (!returnStatus) {
4669    if (numberUnsatisfied) {
4670      // create branching object
4671      const OsiObject * obj = model->solver()->object(choose->bestObjectIndex());
4672      //const OsiSolverInterface * solver = usefulInfo->solver_;
4673      branch_ = obj->createBranch(model->solver(),usefulInfo,obj->whichWay());
4674    }
4675  }
4676  usefulInfo->solution_=saveInfoSol;
4677  delete [] saveSolution;
4678  // may have got solution
4679  if (choose->goodSolution()
4680      &&model->problemFeasibility()->feasible(model,-1)>=0) {
4681    // yes
4682    double objValue = choose->goodObjectiveValue();
4683    model->setBestSolution(CBC_STRONGSOL,
4684                           objValue,
4685                           choose->goodSolution()) ;
4686    model->setLastHeuristic(NULL);
4687    model->incrementUsed(choose->goodSolution());
4688    choose->clearGoodSolution();
4689  }
4690  return returnStatus;
4691}
4692int 
4693CbcNode::chooseClpBranch (CbcModel * model,
4694                          CbcNode * lastNode)
4695{ 
4696  assert(lastNode);
4697  depth_ = lastNode->depth_+1;
4698  delete branch_;
4699  branch_=NULL;
4700  OsiSolverInterface * solver = model->solver();
4701  //double saveObjectiveValue = solver->getObjValue();
4702  //double objectiveValue = CoinMax(solver->getObjSense()*saveObjectiveValue,objectiveValue_);
4703  const double * lower = solver->getColLower();
4704  const double * upper = solver->getColUpper();
4705  // point to useful information
4706  OsiBranchingInformation usefulInfo = model->usefulInformation();
4707  // and modify
4708  usefulInfo.depth_=depth_;
4709  int i;
4710  //bool beforeSolution = model->getSolutionCount()==0;
4711  int numberObjects = model->numberObjects();
4712  int numberColumns = model->getNumCols();
4713  double * saveUpper = new double[numberColumns];
4714  double * saveLower = new double[numberColumns];
4715 
4716  // Save solution in case heuristics need good solution later
4717 
4718  double * saveSolution = new double[numberColumns];
4719  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
4720  model->reserveCurrentSolution(saveSolution);
4721  for (i=0;i<numberColumns;i++) {
4722    saveLower[i] = lower[i];
4723    saveUpper[i] = upper[i];
4724  }
4725  // Save basis
4726  CoinWarmStart * ws = solver->getWarmStart();
4727  numberUnsatisfied_ = 0;
4728  // initialize sum of "infeasibilities"
4729  sumInfeasibilities_ = 0.0;
4730  // Note looks as if off end (hidden one)
4731  OsiObject * object = model->modifiableObject(numberObjects);
4732  CbcGeneralDepth * thisOne = dynamic_cast <CbcGeneralDepth *> (object);
4733  assert (thisOne);
4734  OsiClpSolverInterface * clpSolver
4735    = dynamic_cast<OsiClpSolverInterface *> (solver);
4736  assert (clpSolver);
4737  ClpSimplex * simplex = clpSolver->getModelPtr();
4738  int preferredWay;
4739  double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
4740  if (thisOne->whichSolution()>=0) {
4741    ClpNode * nodeInfo = thisOne->nodeInfo(thisOne->whichSolution());
4742    nodeInfo->applyNode(simplex,2);
4743    int saveLogLevel = simplex->logLevel();
4744    simplex->setLogLevel(0);
4745    simplex->dual();
4746    simplex->setLogLevel(saveLogLevel);
4747    double cutoff = model->getCutoff();
4748    bool goodSolution=true;
4749    if (simplex->status()) {
4750      //simplex->writeMps("bad7.mps",2);
4751      if (nodeInfo->objectiveValue()>cutoff-1.0e-2)
4752        goodSolution=false;
4753      else
4754        assert (!simplex->status());
4755    }
4756    if (goodSolution) {
4757      double newObjectiveValue = solver->getObjSense()*solver->getObjValue();
4758      // See if integer solution
4759      int numInf;
4760      int numInf2;
4761      bool gotSol = model->feasibleSolution(numInf,numInf2);
4762      if (!gotSol) {
4763        printf("numinf %d\n",numInf);
4764        double * sol = simplex->primalColumnSolution();
4765        for (int i=0;i<numberColumns;i++) {
4766          if (simplex->isInteger(i)) {
4767            double value = floor(sol[i]+0.5);
4768            if (fabs(value-sol[i])>1.0e-7) {
4769              printf("%d value %g\n",i,sol[i]);
4770              if (fabs(value-sol[i])<1.0e-3) {
4771                sol[i]=value;
4772              }
4773            }
4774          }
4775        }
4776        simplex->writeMps("bad8.mps",2);
4777        bool gotSol = model->feasibleSolution(numInf,numInf2);
4778        if (!gotSol) 
4779          assert (gotSol);
4780      }
4781      model->setBestSolution(CBC_STRONGSOL,
4782                             newObjectiveValue,
4783                             solver->getColSolution()) ;
4784      model->setLastHeuristic(NULL);
4785      model->incrementUsed(solver->getColSolution());
4786    }
4787  }
4788  // restore bounds
4789  { for (int j=0;j<numberColumns;j++) {
4790      if (saveLower[j] != lower[j])
4791        solver->setColLower(j,saveLower[j]);
4792      if (saveUpper[j] != upper[j])
4793        solver->setColUpper(j,saveUpper[j]);
4794    }
4795  }
4796  // restore basis
4797  solver->setWarmStart(ws);
4798  delete ws;
4799  int anyAction;
4800  //#define CHECK_PATH
4801#ifdef CHECK_PATH
4802  extern int gotGoodNode_Z;
4803  if (gotGoodNode_Z>=0)
4804    printf("good node %d %g\n",gotGoodNode_Z,infeasibility);
4805#endif
4806  if (infeasibility>0.0) {
4807    if (infeasibility==COIN_DBL_MAX) {
4808      anyAction=-2; // infeasible
4809    } else {
4810      branch_=thisOne->createCbcBranch(solver,&usefulInfo,preferredWay);
4811      // Set to firat one (and change when re-pushing)
4812      CbcGeneralBranchingObject * branch = 
4813        dynamic_cast <CbcGeneralBranchingObject *> (branch_);
4814      branch->state(objectiveValue_,sumInfeasibilities_,
4815                    numberUnsatisfied_,0);
4816      branch->setNode(this);
4817      anyAction=0;
4818    }
4819  } else {
4820    anyAction=-1;
4821  }
4822#ifdef CHECK_PATH
4823  gotGoodNode_Z=-1;
4824#endif
4825  // Set guessed solution value
4826  guessedObjectiveValue_ = objectiveValue_+1.0e-5;
4827  delete [] saveLower;
4828  delete [] saveUpper;
4829 
4830  // restore solution
4831  solver->setColSolution(saveSolution);
4832  delete [] saveSolution;
4833  return anyAction;
4834}
4835/* Double checks in case node can change its mind!
4836   Returns objective value
4837   Can change objective etc */
4838double
4839CbcNode::checkIsCutoff(double cutoff)
4840{
4841  branch_->checkIsCutoff(cutoff);
4842  return objectiveValue_;
4843}
Note: See TracBrowser for help on using the repository browser.