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

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

Changes to try and make faster.
ALSO createBranch is now createCbcBranch. I don't think anyone uses derived classes from CbcBranchBase? but if they do they will get annoyed at me and have to change a few names. I am just trying to get rid of all but all warnings

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 161.1 KB
Line 
1/* $Id: CbcNode.cpp 1205 2009-07-30 15:19:11Z 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  if (numberObjects>model->numberIntegers()) {
2410    for (int i=model->numberIntegers();i<numberObjects;i++) {
2411      OsiObject * object = model->modifiableObject(i);
2412      CbcObject * obj = dynamic_cast <CbcObject *>(object) ;
2413      if (!obj || !obj->optionalObject()) {
2414        useOldWay=true;
2415        break;
2416      }
2417    }
2418  }
2419  if ((model->specialOptions()&128)!=0)
2420    useOldWay=false; // allow
2421  // For now return if not simple
2422  if (useOldWay)
2423    return -3;
2424  // point to useful information
2425  OsiBranchingInformation usefulInfo = model->usefulInformation();
2426  // and modify
2427  usefulInfo.depth_=depth_;
2428  if ((model->specialOptions()&128)!=0) {
2429    // SOS - shadow prices
2430    int numberRows = solver->getNumRows();
2431    const double * pi = usefulInfo.pi_;
2432    double sumPi=0.0;
2433    for (int i=0;i<numberRows;i++) 
2434      sumPi += fabs(pi[i]);
2435    sumPi /= static_cast<double> (numberRows);
2436    // and scale back
2437    sumPi *= 0.01;
2438    usefulInfo.defaultDual_ = sumPi; // switch on
2439    int numberColumns = solver->getNumCols();
2440    int size = CoinMax(numberColumns,2*numberRows);
2441    usefulInfo.usefulRegion_ = new double [size];
2442    CoinZeroN(usefulInfo.usefulRegion_,size);
2443    usefulInfo.indexRegion_ = new int [size];
2444    // pi may change
2445    usefulInfo.pi_=CoinCopyOfArray(usefulInfo.pi_,numberRows);
2446  }
2447  assert (auxiliaryInfo);
2448  double cutoff =model->getCutoff();
2449  const double * lower = solver->getColLower();
2450  const double * upper = solver->getColUpper();
2451  // See if user thinks infeasible
2452  int anyAction=model->problemFeasibility()->feasible(model,0);
2453  if (anyAction) {
2454    // will return -2 if infeasible , 0 if treat as integer
2455    return anyAction-1;
2456  }
2457  int i;
2458  int saveStateOfSearch = model->stateOfSearch()%10;
2459  int numberStrong=model->numberStrong();
2460  /* Ranging is switched off.
2461     The idea is that you can find out the effect of one iteration
2462     on each unsatisfied variable cheaply.  Then use this
2463     if you have not got much else to go on.
2464  */
2465  //#define RANGING
2466#ifdef RANGING
2467  // must have clp
2468#ifndef COIN_HAS_CLP
2469#  warning("Ranging switched off as not Clp");
2470#undef RANGING
2471#endif
2472  // Pass number
2473  int kPass=0;
2474  int numberRows = solver->getNumRows();
2475#endif
2476  int numberColumns = model->getNumCols();
2477  double * saveUpper = new double[numberColumns];
2478  double * saveLower = new double[numberColumns];
2479  for (i=0;i<numberColumns;i++) {
2480    saveLower[i] = lower[i];
2481    saveUpper[i] = upper[i];
2482  }
2483 
2484  // Save solution in case heuristics need good solution later
2485 
2486  double * saveSolution = new double[numberColumns];
2487  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
2488  model->reserveCurrentSolution(saveSolution);
2489  /*
2490    Get a branching decision object. Use the default dynamic decision criteria unless
2491    the user has loaded a decision method into the model.
2492  */
2493  CbcBranchDecision *decision = model->branchingMethod();
2494  if (!decision)
2495    decision = new CbcBranchDynamicDecision();
2496  int xMark=0;
2497  // Get arrays to sort
2498  double * sort = new double[numberObjects];
2499  int * whichObject = new int[numberObjects];
2500#ifdef RANGING
2501  int xPen=0;
2502  int * objectMark = new int[2*numberObjects+1];
2503#endif
2504  // Arrays with movements
2505  double * upEstimate = new double[numberObjects];
2506  double * downEstimate = new double[numberObjects];
2507  double estimatedDegradation=0.0; 
2508  int numberNodes=model->getNodeCount();
2509  int saveLogLevel = model->logLevel();
2510#if 0
2511  if ((numberNodes%500)==0) {
2512    model->setLogLevel(6);
2513    // Get average up and down costs
2514    double averageUp=0.0;
2515    double averageDown=0.0;
2516    int numberUp=0;
2517    int numberDown=0;
2518    int i;
2519    for ( i=0;i<numberObjects;i++) {
2520      OsiObject * object = model->modifiableObject(i);
2521      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2522        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2523      assert(dynamicObject);
2524      int  numberUp2=0;
2525      int numberDown2=0;
2526      double up=0.0;
2527      double down=0.0;
2528      if (dynamicObject->numberTimesUp()) {
2529        numberUp++;
2530        averageUp += dynamicObject->upDynamicPseudoCost();
2531        numberUp2 += dynamicObject->numberTimesUp();
2532        up = dynamicObject->upDynamicPseudoCost();
2533      }
2534      if (dynamicObject->numberTimesDown()) {
2535        numberDown++;
2536        averageDown += dynamicObject->downDynamicPseudoCost();
2537        numberDown2 += dynamicObject->numberTimesDown();
2538        down = dynamicObject->downDynamicPseudoCost();
2539      }
2540      if (numberUp2||numberDown2)
2541        printf("col %d - up %d times cost %g, - down %d times cost %g\n",
2542               dynamicObject->columnNumber(),numberUp2,up,numberDown2,down);
2543    }
2544    if (numberUp)
2545      averageUp /= static_cast<double> (numberUp);
2546    else
2547      averageUp=1.0;
2548    if (numberDown)
2549      averageDown /= static_cast<double> (numberDown);
2550    else
2551      averageDown=1.0;
2552    printf("total - up %d vars average %g, - down %d vars average %g\n",
2553           numberUp,averageUp,numberDown,averageDown);
2554  }
2555#endif
2556  int numberBeforeTrust = model->numberBeforeTrust();
2557  // May go round twice if strong branching fixes all local candidates
2558  bool finished=false;
2559  int numberToFix=0;
2560# ifdef COIN_HAS_CLP
2561  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
2562  int saveClpOptions=0;
2563  if (osiclp) {
2564    // for faster hot start
2565    saveClpOptions = osiclp->specialOptions();
2566    osiclp->setSpecialOptions(saveClpOptions|8192);
2567  }
2568# else
2569  OsiSolverInterface *osiclp = NULL ;
2570# endif
2571  const CglTreeProbingInfo * probingInfo = NULL; //model->probingInfo();
2572  // Old code left in with DEPRECATED_STRATEGY
2573  assert (model->searchStrategy()==-1||
2574          model->searchStrategy()==1||
2575          model->searchStrategy()==2);
2576#ifdef DEPRECATED_STRATEGY
2577  int saveSearchStrategy2 = model->searchStrategy();
2578#endif
2579  // Get average up and down costs
2580  double averageUp=0.0;
2581  double averageDown=0.0;
2582  {
2583    int numberUp=0;
2584    int numberDown=0;
2585    int i;
2586    for ( i=0;i<numberObjects;i++) {
2587      OsiObject * object = model->modifiableObject(i);
2588      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2589        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2590      if (dynamicObject) {
2591        if (dynamicObject->numberTimesUp()) {
2592          numberUp++;
2593          averageUp += dynamicObject->upDynamicPseudoCost();
2594        }
2595        if (dynamicObject->numberTimesDown()) {
2596          numberDown++;
2597          averageDown += dynamicObject->downDynamicPseudoCost();
2598        }
2599      }
2600    }
2601    if (numberUp) 
2602      averageUp /= static_cast<double> (numberUp);
2603    else
2604      averageUp=1.0;
2605    if (numberDown) 
2606      averageDown /= static_cast<double> (numberDown);
2607    else
2608      averageDown=1.0;
2609    for ( i=0;i<numberObjects;i++) {
2610      OsiObject * object = model->modifiableObject(i);
2611      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2612        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2613      if (dynamicObject) {
2614        if (!dynamicObject->numberTimesUp()) 
2615          dynamicObject->setUpDynamicPseudoCost(averageUp);
2616        if (!dynamicObject->numberTimesDown()) 
2617          dynamicObject->setDownDynamicPseudoCost(averageDown);
2618      }
2619    }
2620  }
2621#ifdef DEPRECATED_STRATEGY
2622  { // in for tabbing
2623  } else if (saveSearchStrategy2<1999) {
2624    // pseudo shadow prices
2625    model->pseudoShadow(NULL,NULL);
2626  } else if (saveSearchStrategy2<2999) {
2627    // leave old ones
2628  } else if (saveSearchStrategy2<3999) {
2629    // pseudo shadow prices at root
2630    if (!numberNodes)
2631      model->pseudoShadow(NULL,NULL);
2632  } else {
2633    abort();
2634  }
2635  if (saveSearchStrategy2>=0)
2636    saveSearchStrategy2 = saveSearchStrategy2 % 1000;
2637  if (saveSearchStrategy2==999)
2638    saveSearchStrategy2=-1;
2639  int saveSearchStrategy = saveSearchStrategy2<99 ? saveSearchStrategy2 : saveSearchStrategy2-100;
2640#endif //DEPRECATED_STRATEGY
2641  int numberNotTrusted=0;
2642  int numberStrongDone=0;
2643  int numberUnfinished=0;
2644  int numberStrongInfeasible=0;
2645  int numberStrongIterations=0;
2646  // so we can save lots of stuff
2647  CbcStrongInfo choice;
2648  CbcDynamicPseudoCostBranchingObject * choiceObject = NULL;
2649  if (model->allDynamic()) {
2650    CbcSimpleIntegerDynamicPseudoCost * object = NULL;
2651    choiceObject=new CbcDynamicPseudoCostBranchingObject(model,0,-1,0.5,object);
2652  }
2653  choice.possibleBranch=choiceObject;
2654  numberPassesLeft = CoinMax(numberPassesLeft,2);
2655  while(!finished) {
2656    numberPassesLeft--;
2657    finished=true;
2658    decision->initialize(model);
2659    // Some objects may compute an estimate of best solution from here
2660    estimatedDegradation=0.0; 
2661    numberToFix=0;
2662    int numberToDo=0;
2663    int iBestNot=-1;
2664    int iBestGot=-1;
2665    double best=0.0;
2666    numberNotTrusted=0;
2667    numberStrongDone=0;
2668    numberUnfinished=0;
2669    numberStrongInfeasible=0;
2670    numberStrongIterations=0;
2671#ifdef RANGING
2672    int * which = objectMark+numberObjects+1;
2673    int neededPenalties;
2674    int optionalPenalties;
2675#endif
2676    // We may go round this loop three times (only if we think we have solution)
2677    for (int iPass=0;iPass<3;iPass++) {
2678     
2679      // Some objects may compute an estimate of best solution from here
2680      estimatedDegradation=0.0; 
2681      numberUnsatisfied_ = 0;
2682      // initialize sum of "infeasibilities"
2683      sumInfeasibilities_ = 0.0;
2684      int bestPriority=COIN_INT_MAX;
2685      int number01 = 0;
2686      const cliqueEntry * entry = NULL;
2687      const int * toZero = NULL;
2688      const int * toOne = NULL;
2689      const int * backward = NULL;
2690      int numberUnsatisProbed=0;
2691      int numberUnsatisNotProbed=0; // 0-1
2692      if (probingInfo) {
2693        number01 = probingInfo->numberIntegers();
2694        entry = probingInfo->fixEntries();
2695        toZero = probingInfo->toZero();
2696        toOne = probingInfo->toOne();
2697        backward = probingInfo->backward();
2698        if (!toZero[number01]||number01<numberObjects||true) {
2699          // no info
2700          probingInfo=NULL;
2701        }
2702      }
2703      /*
2704        Scan for branching objects that indicate infeasibility. Choose candidates
2705        using priority as the first criteria, then integer infeasibility.
2706       
2707        The algorithm is to fill the array with a set of good candidates (by
2708        infeasibility) with priority bestPriority.  Finding a candidate with
2709        priority better (less) than bestPriority flushes the choice array. (This
2710        serves as initialization when the first candidate is found.)
2711       
2712      */
2713      numberToDo=0;
2714#ifdef RANGING
2715      neededPenalties=0;
2716      optionalPenalties=numberObjects;
2717#endif
2718      iBestNot=-1;
2719      double bestNot=0.0;
2720      iBestGot=-1;
2721      best=0.0;
2722      /* Problem type as set by user or found by analysis.  This will be extended
2723         0 - not known
2724         1 - Set partitioning <=
2725         2 - Set partitioning ==
2726         3 - Set covering
2727         4 - all +- 1 or all +1 and odd
2728      */
2729      int problemType = model->problemType();
2730      for (i=0;i<numberObjects;i++) {
2731        OsiObject * object = model->modifiableObject(i);
2732        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2733          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2734        int preferredWay;
2735        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
2736        int priorityLevel = object->priority();
2737#define ZERO_ONE 0
2738#define ZERO_FAKE 1.0e20;
2739#if ZERO_ONE==1
2740        // branch on 0-1 first (temp)
2741        if (fabs(saveSolution[dynamicObject->columnNumber()])<1.0)
2742          priorityLevel--;
2743#endif
2744#if ZERO_ONE==2
2745        if (fabs(saveSolution[dynamicObject->columnNumber()])<1.0)
2746          infeasibility *= ZERO_FAKE;
2747#endif
2748        if (infeasibility) {
2749          int iColumn = numberColumns+i;
2750          bool gotDown=false;
2751          int numberThisDown = 0;
2752          bool gotUp=false;
2753          int numberThisUp = 0;
2754          double downGuess=object->downEstimate();
2755          double upGuess=object->upEstimate();
2756          if (dynamicObject) {
2757            // Use this object's numberBeforeTrust
2758            int numberBeforeTrust = dynamicObject->numberBeforeTrust();
2759            iColumn = dynamicObject->columnNumber();
2760            gotDown=false;
2761            numberThisDown = dynamicObject->numberTimesDown();
2762            if (numberThisDown>=numberBeforeTrust)
2763              gotDown=true;
2764            gotUp=false;
2765            numberThisUp = dynamicObject->numberTimesUp();
2766            if (numberThisUp>=numberBeforeTrust)
2767              gotUp=true;
2768            if (!depth_&&false) {
2769              // try closest to 0.5
2770              double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2771              infeasibility = fabs(0.5-part);
2772            }
2773            if (problemType>0&&problemType<4&&false) {
2774              // try closest to 0.5
2775              double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2776              infeasibility = 0.5-fabs(0.5-part);
2777            }
2778            if (probingInfo) {
2779              int iSeq = backward[iColumn];
2780              assert (iSeq>=0);
2781              infeasibility = 1.0 + (toZero[iSeq+1]-toZero[iSeq])+ 
2782                5.0*CoinMin(toOne[iSeq]-toZero[iSeq],toZero[iSeq+1]-toOne[iSeq]);
2783              if (toZero[iSeq+1]>toZero[iSeq]) {
2784                numberUnsatisProbed++;
2785              } else {
2786                numberUnsatisNotProbed++;
2787              }
2788            }
2789          } else {
2790            // see if SOS
2791            CbcSOS * sosObject =
2792              dynamic_cast <CbcSOS *>(object) ;
2793            if (sosObject) {
2794              gotDown=false;
2795              numberThisDown = sosObject->numberTimesDown();
2796              if (numberThisDown>=numberBeforeTrust)
2797                gotDown=true;
2798              gotUp=false;
2799              numberThisUp = sosObject->numberTimesUp();
2800              if (numberThisUp>=numberBeforeTrust)
2801                gotUp=true;
2802            } else {
2803              gotDown=true;
2804              numberThisDown=999999;
2805              downGuess=1.0e20;
2806              gotUp=true;
2807              numberThisUp=999999;
2808              upGuess=1.0e20;
2809              numberPassesLeft=0;
2810            }
2811          }
2812          // Increase estimated degradation to solution
2813          estimatedDegradation += CoinMin(downGuess,upGuess);
2814          downEstimate[i]=downGuess;
2815          upEstimate[i]=upGuess;
2816          numberUnsatisfied_++;
2817          sumInfeasibilities_ += infeasibility;
2818          // Better priority? Flush choices.
2819          if (priorityLevel<bestPriority) {
2820            numberToDo=0;
2821            bestPriority = priorityLevel;
2822            iBestGot=-1;
2823            best=0.0;
2824            numberNotTrusted=0;
2825#ifdef RANGING
2826            neededPenalties=0;
2827            optionalPenalties=numberObjects;
2828#endif
2829          } else if (priorityLevel>bestPriority) {
2830            continue;
2831          }
2832          if (!gotUp||!gotDown)
2833            numberNotTrusted++;
2834          // Check for suitability based on infeasibility.
2835          if ((gotDown&&gotUp)&&numberStrong>0) {
2836            sort[numberToDo]=-infeasibility;
2837            if (infeasibility>best) {
2838              best=infeasibility;
2839              iBestGot=numberToDo;
2840            }
2841#ifdef RANGING
2842            if (dynamicObject) {
2843              objectMark[--optionalPenalties]=numberToDo;
2844              which[optionalPenalties]=iColumn;
2845            }
2846#endif
2847          } else {
2848#ifdef RANGING
2849            if (dynamicObject) {
2850              objectMark[neededPenalties]=numberToDo;
2851              which[neededPenalties++]=iColumn;
2852            }
2853#endif
2854            sort[numberToDo]=-10.0*infeasibility;
2855            if (!(numberThisUp+numberThisDown))
2856              sort[numberToDo] *= 100.0; // make even more likely
2857            if (iColumn<numberColumns) {
2858              double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2859              if (1.0-fabs(part-0.5)>bestNot) {
2860                iBestNot=numberToDo;
2861                bestNot = 1.0-fabs(part-0.5);
2862              }
2863            } else {
2864              // SOS
2865              if (-sort[numberToDo]>bestNot) {
2866                iBestNot=numberToDo;
2867                bestNot = -sort[numberToDo];
2868              }
2869            }
2870          }
2871          if (model->messageHandler()->logLevel()>3) { 
2872            printf("%d (%d) down %d %g up %d %g - infeas %g - sort %g solution %g\n",
2873                   i,iColumn,numberThisDown,object->downEstimate(),numberThisUp,object->upEstimate(),
2874                   infeasibility,sort[numberToDo],saveSolution[iColumn]);
2875          }
2876          whichObject[numberToDo++]=i;
2877        } else {
2878          // for debug
2879          downEstimate[i]=-1.0;
2880          upEstimate[i]=-1.0;
2881        }
2882      }
2883      if (numberUnsatisfied_) {
2884        if (probingInfo&&false)
2885          printf("nunsat %d, %d probed, %d other 0-1\n",numberUnsatisfied_,
2886                 numberUnsatisProbed,numberUnsatisNotProbed);
2887        // some infeasibilities - go to next steps
2888        break;
2889      } else if (!iPass) {
2890        // may just need resolve
2891        model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
2892        if (!solver->isProvenOptimal()) {
2893          // infeasible
2894          anyAction=-2;
2895          break;
2896        }
2897      } else if (iPass==1) {
2898        // looks like a solution - get paranoid
2899        bool roundAgain=false;
2900        // get basis
2901        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
2902        if (!ws)
2903          break;
2904        double tolerance;
2905        solver->getDblParam(OsiPrimalTolerance,tolerance);
2906        for (i=0;i<numberColumns;i++) {
2907          double value = saveSolution[i];
2908          if (value<lower[i]-tolerance) {
2909            saveSolution[i]=lower[i];
2910            roundAgain=true;
2911            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
2912          } else if (value>upper[i]+tolerance) {
2913            saveSolution[i]=upper[i];
2914            roundAgain=true;
2915            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
2916          } 
2917        }
2918        if (roundAgain) {
2919          // restore basis
2920          solver->setWarmStart(ws);
2921          solver->setColSolution(saveSolution);
2922          delete ws;
2923          bool takeHint;
2924          OsiHintStrength strength;
2925          solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
2926          solver->setHintParam(OsiDoDualInResolve,false,OsiHintDo) ;
2927          model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
2928          solver->setHintParam(OsiDoDualInResolve,takeHint,strength) ;
2929          if (!solver->isProvenOptimal()) {
2930            // infeasible
2931            anyAction=-2;
2932            break;
2933          }
2934        } else {
2935          delete ws;
2936          break;
2937        }
2938      }
2939    }
2940    if (anyAction==-2) {
2941      break;
2942    }
2943    // skip if solution
2944    if (!numberUnsatisfied_)
2945      break;
2946    int skipAll = (numberNotTrusted==0||numberToDo==1) ? 1 : 0;
2947    bool doneHotStart=false;
2948    //DEPRECATED_STRATEGYint searchStrategy = saveSearchStrategy>=0 ? (saveSearchStrategy%10) : -1;
2949    int searchStrategy = model->searchStrategy();
2950    // But adjust depending on ratio of iterations
2951    if (searchStrategy>0) {
2952      if (numberBeforeTrust>=5&&numberBeforeTrust<=10) {
2953        if (searchStrategy!=2) {
2954          assert (searchStrategy==1);
2955          if (depth_>5) {
2956            int numberIterations = model->getIterationCount();
2957            int numberStrongIterations = model->numberStrongIterations();
2958            if (numberStrongIterations>numberIterations+10000) {
2959              searchStrategy=2;
2960              skipAll=1;
2961            } else if (numberStrongIterations*4+1000<numberIterations) {
2962              searchStrategy=3;
2963              skipAll=0;
2964            }
2965          } else {
2966            searchStrategy=3;
2967            skipAll=0;
2968          }
2969        }
2970      }
2971    }
2972    // worth trying if too many times
2973    // Save basis
2974    CoinWarmStart * ws = NULL;
2975    // save limit
2976    int saveLimit=0;
2977    solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
2978    if (!numberPassesLeft)
2979      skipAll=1;
2980    if (!skipAll) {
2981      ws = solver->getWarmStart();
2982      int limit=100;
2983      if (!saveStateOfSearch&&saveLimit<limit&&saveLimit==100)
2984        solver->setIntParam(OsiMaxNumIterationHotStart,limit); 
2985    }
2986    // Say which one will be best
2987    int whichChoice=0;
2988    int bestChoice;
2989    if (iBestGot>=0)
2990      bestChoice=iBestGot;
2991    else
2992      bestChoice=iBestNot;
2993    assert (bestChoice>=0);
2994    // If we have hit max time don't do strong branching
2995    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
2996                        model->getDblParam(CbcModel::CbcMaximumSeconds));
2997    // also give up if we are looping round too much
2998    if (hitMaxTime||numberPassesLeft<=0) {
2999      int iObject = whichObject[bestChoice];
3000      OsiObject * object = model->modifiableObject(iObject);
3001      int preferredWay;
3002      object->infeasibility(&usefulInfo,preferredWay);
3003      CbcObject * obj =
3004        dynamic_cast <CbcObject *>(object) ;
3005      assert (obj);
3006      branch_=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3007      {
3008        CbcBranchingObject * branchObj =
3009          dynamic_cast <CbcBranchingObject *>(branch_) ;
3010        assert (branchObj);
3011        branchObj->way(preferredWay);
3012      }
3013      delete ws;
3014      ws=NULL;
3015      break;
3016    } else {
3017      // say do fast
3018      int easy=1;
3019      if (!skipAll)
3020        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
3021      int iDo;
3022#ifdef RANGING
3023      bool useRanging = model->allDynamic()&&!skipAll;
3024      if (useRanging) {
3025        double currentObjective = solver->getObjValue()*solver->getObjSense();
3026        double gap = cutoff-currentObjective;
3027        // relax a bit
3028        gap *= 1.0000001;
3029        gap = CoinMax(1.0e-5,gap);
3030        // off penalties if too much
3031        double needed = neededPenalties;
3032        needed *= numberRows;
3033        if (numberNodes) {
3034          if (needed>1.0e6) {
3035            neededPenalties=0;
3036          } else if (gap<1.0e5) {
3037            // maybe allow some not needed
3038            int extra = static_cast<int> ((1.0e6-needed)/numberRows);
3039            int nStored = numberObjects-optionalPenalties;
3040            extra = CoinMin(extra,nStored);
3041            for (int i=0;i<extra;i++) {
3042              objectMark[neededPenalties]=objectMark[optionalPenalties+i];
3043                which[neededPenalties++]=which[optionalPenalties+i];;
3044            }
3045          }
3046        }
3047        if (osiclp&&neededPenalties) {
3048          assert (!doneHotStart);
3049          xPen += neededPenalties;
3050          which--;
3051          which[0]=neededPenalties;
3052          osiclp->passInRanges(which);
3053          // Mark hot start and get ranges
3054          if (kPass) {
3055            // until can work out why solution can go funny
3056            int save = osiclp->specialOptions();
3057            osiclp->setSpecialOptions(save|256);
3058            solver->markHotStart();
3059            osiclp->setSpecialOptions(save);
3060          } else {
3061            solver->markHotStart();
3062          }
3063          doneHotStart=true;
3064          xMark++;
3065          kPass++;
3066          osiclp->passInRanges(NULL);
3067          const double * downCost=osiclp->upRange();
3068          const double * upCost=osiclp->downRange();
3069          bool problemFeasible=true;
3070          int numberFixed=0;
3071          for (int i=0;i<neededPenalties;i++) {
3072            int j = objectMark[i];
3073            int iObject = whichObject[j];
3074            OsiObject * object = model->modifiableObject(iObject);
3075            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3076              static_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3077            // Use this object's numberBeforeTrust
3078            int numberBeforeTrust = dynamicObject->numberBeforeTrust();
3079            int iSequence=dynamicObject->columnNumber();
3080            double value = saveSolution[iSequence];
3081            value -= floor(value);
3082            double upPenalty = CoinMin(upCost[i],1.0e110)*(1.0-value);
3083            double downPenalty = CoinMin(downCost[i],1.0e110)*value;
3084            int numberThisDown = dynamicObject->numberTimesDown();
3085            int numberThisUp = dynamicObject->numberTimesUp();
3086            if (!numberBeforeTrust) {
3087              // override
3088              downEstimate[iObject]=downPenalty;
3089              upEstimate[iObject]=upPenalty;
3090              double min1 = CoinMin(downEstimate[iObject],
3091                                    upEstimate[iObject]);
3092              double max1 = CoinMax(downEstimate[iObject],
3093                                    upEstimate[iObject]);
3094              min1 = 0.8*min1+0.2*max1;
3095              sort[j] = - min1;
3096            } else if (numberThisDown<numberBeforeTrust||
3097                       numberThisUp<numberBeforeTrust) {
3098              double invTrust = 1.0/static_cast<double> (numberBeforeTrust);
3099              if (numberThisDown<numberBeforeTrust) {
3100                double fraction = numberThisDown*invTrust;
3101                downEstimate[iObject] = fraction*downEstimate[iObject]+(1.0-fraction)*downPenalty;
3102              }
3103              if (numberThisUp<numberBeforeTrust) {
3104                double fraction = numberThisUp*invTrust;
3105                upEstimate[iObject] = fraction*upEstimate[iObject]+(1.0-fraction)*upPenalty;
3106              }
3107              double min1 = CoinMin(downEstimate[iObject],
3108                                    upEstimate[iObject]);
3109              double max1 = CoinMax(downEstimate[iObject],
3110                                    upEstimate[iObject]);
3111              min1 = 0.8*min1+0.2*max1;
3112              min1 *= 10.0;
3113              if (!(numberThisDown+numberThisUp))
3114                min1 *= 100.0;
3115              sort[j] = - min1;
3116            }
3117            if (CoinMax(downPenalty,upPenalty)>gap) {
3118              printf("gap %g object %d has down range %g, up %g\n",
3119                     gap,i,downPenalty,upPenalty);
3120              //sort[j] -= 1.0e50; // make more likely to be chosen
3121              int number;
3122              if (downPenalty>gap) {
3123                number = dynamicObject->numberTimesDown();
3124                if (upPenalty>gap) 
3125                  problemFeasible=false;
3126                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver,&usefulInfo,1);
3127                //branch->fix(solver,saveLower,saveUpper,1);
3128                delete branch;
3129              } else {
3130                number = dynamicObject->numberTimesUp();
3131                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver,&usefulInfo,1);
3132                //branch->fix(solver,saveLower,saveUpper,-1);
3133                delete branch;
3134              }
3135              if(number>=numberBeforeTrust)
3136                dynamicObject->setNumberBeforeTrust(number+1);
3137              numberFixed++;
3138            }
3139            if (!numberNodes)
3140              printf("%d pen down ps %g -> %g up ps %g -> %g\n",
3141                     iObject,downPenalty,downPenalty,upPenalty,upPenalty);
3142          }
3143          if(numberFixed&&problemFeasible) {
3144            assert(doneHotStart);
3145            solver->unmarkHotStart();
3146            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3147            solver->markHotStart();
3148            problemFeasible = solver->isProvenOptimal();
3149          }
3150          if(!problemFeasible) {
3151            fprintf(stdout,"both ways infeas on ranging - code needed\n");
3152            anyAction=-2;
3153            if (!choiceObject) {
3154              delete choice.possibleBranch;
3155              choice.possibleBranch=NULL;
3156            }
3157            //printf("Both infeasible for choice %d sequence %d\n",i,
3158            // model->object(choice.objectNumber)->columnNumber());
3159            // Delete the snapshot
3160            solver->unmarkHotStart();
3161            // back to normal
3162            solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3163            // restore basis
3164            solver->setWarmStart(ws);
3165            doneHotStart=false;
3166            delete ws;
3167            ws=NULL;
3168            break;
3169          }
3170        }
3171      }
3172#endif          /* RANGING */
3173      {
3174        int numberIterations = model->getIterationCount();
3175        //numberIterations += (model->numberExtraIterations()>>2);
3176        const int * strongInfo = model->strongInfo();
3177        //int numberDone = strongInfo[0]-strongInfo[3];
3178        int numberFixed = strongInfo[1]-strongInfo[4];
3179        int numberInfeasible = strongInfo[2]-strongInfo[5];
3180        assert (!strongInfo[3]);
3181        assert (!strongInfo[4]);
3182        assert (!strongInfo[5]);
3183        int numberStrongIterations = model->numberStrongIterations();
3184        int numberRows = solver->getNumRows();
3185        if (numberStrongIterations>numberIterations+CoinMin(100,10*numberRows)&&depth_>=4&&numberNodes>100) {
3186          if (20*numberInfeasible+4*numberFixed<numberNodes) {
3187            // Say never do
3188            skipAll=-1;
3189          }
3190        } 
3191      }
3192      // make sure best will be first
3193      if (iBestGot>=0)
3194        sort[iBestGot]=-COIN_DBL_MAX;
3195      // Actions 0 - exit for repeat, 1 resolve and try old choice,2 exit for continue
3196      if (anyAction)
3197        numberToDo=0; // skip as we will be trying again
3198      // Sort
3199      CoinSort_2(sort,sort+numberToDo,whichObject);
3200      // Change in objective opposite infeasible
3201      double worstFeasible=0.0;
3202      // Just first if strong off
3203      if (!numberStrong)
3204        numberToDo=CoinMin(numberToDo,1);
3205      if (searchStrategy==2)
3206        numberToDo=CoinMin(numberToDo,20);
3207      iDo=0;
3208      int saveLimit2;
3209      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit2);
3210      bool doQuickly = false; // numberToDo>2*numberStrong;
3211      if (searchStrategy==2)
3212        doQuickly=true;
3213      int numberTest=numberNotTrusted>0 ? numberStrong : (numberStrong+1)/2;
3214      if (searchStrategy==3) {
3215        // Previously decided we need strong
3216        doQuickly=false;
3217        numberTest = numberStrong;
3218      }
3219      // Try nearly always off
3220      if (skipAll>=0) {
3221        if (searchStrategy<2) {
3222          //if ((numberNodes%20)!=0) {
3223            if ((model->specialOptions()&8)==0) {
3224              numberTest=0;
3225              doQuickly=true;
3226            }
3227            //} else {
3228            //doQuickly=false;
3229            //numberTest=2*numberStrong;
3230            //skipAll=0;
3231            //}
3232        }
3233      } else {
3234        // Just take first
3235        doQuickly=true;
3236        numberTest=1;
3237      }
3238      int testDepth = (skipAll>=0) ? 8 : 4;
3239      if (depth_<testDepth&&numberStrong) {
3240        if (searchStrategy!=2) {
3241          doQuickly=false;
3242          int numberRows = solver->getNumRows();
3243          // whether to do this or not is important - think
3244          if (numberRows<300||numberRows+numberColumns<2500) {
3245            if (depth_<7)
3246              numberStrong = CoinMin(3*numberStrong,numberToDo);
3247            if (!depth_) 
3248              numberStrong=CoinMin(6*numberStrong,numberToDo);
3249          }
3250          numberTest=numberStrong;
3251          skipAll=0;
3252        }
3253      }
3254      // Do at least 5 strong
3255      if (numberColumns<1000&&(depth_<15||numberNodes<1000000))
3256        numberTest = CoinMax(numberTest,5);
3257      if ((model->specialOptions()&8)==0) {
3258        if (skipAll) {
3259          numberTest=0;
3260          doQuickly=true;
3261        }
3262      } else {
3263        // do 5 as strong is fixing
3264        numberTest = CoinMax(numberTest,5);
3265      }
3266      // see if switched off
3267      if (skipAll<0) {
3268        numberTest=0;
3269        doQuickly=true;
3270      }
3271      int realMaxHotIterations=999999;
3272      if (skipAll<0)
3273        numberToDo=1;
3274#ifdef DO_ALL_AT_ROOT
3275      if (!numberNodes) {
3276        printf("DOX %d test %d unsat %d - nobj %d\n",
3277               numberToDo,numberTest,numberUnsatisfied_,
3278               numberObjects);
3279        numberTest=numberToDo;
3280      }
3281#endif
3282      for ( iDo=0;iDo<numberToDo;iDo++) {
3283        int iObject = whichObject[iDo];
3284        OsiObject * object = model->modifiableObject(iObject);
3285        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3286          static_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3287        int iColumn = dynamicObject ? dynamicObject->columnNumber() : numberColumns+iObject;
3288        int preferredWay;
3289        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
3290        // may have become feasible
3291        if (!infeasibility)
3292          continue;
3293        CbcSimpleInteger * obj =
3294          dynamic_cast <CbcSimpleInteger *>(object) ;
3295        if (obj) {
3296          if (choiceObject) {
3297            obj->fillCreateBranch(choiceObject,&usefulInfo,preferredWay);
3298            choiceObject->setObject(dynamicObject);
3299          } else {
3300            choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3301          }
3302        } else {
3303          CbcObject * obj =
3304            dynamic_cast <CbcObject *>(object) ;
3305          assert (obj);
3306          choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3307        }
3308        // Save which object it was
3309        choice.objectNumber=iObject;
3310        choice.numIntInfeasUp = numberUnsatisfied_;
3311        choice.numIntInfeasDown = numberUnsatisfied_;
3312        choice.upMovement = upEstimate[iObject];
3313        choice.downMovement = downEstimate[iObject];
3314        assert (choice.upMovement>=0.0);
3315        assert (choice.downMovement>=0.0);
3316        choice.fix=0; // say not fixed
3317        // see if can skip strong branching
3318        int canSkip = choice.possibleBranch->fillStrongInfo(choice);
3319        if ((numberTest<=0||skipAll)) {
3320          if (iDo>20) {
3321#ifdef DO_ALL_AT_ROOT
3322            if (!numberNodes)
3323              printf("DOY test %d - iDo %d\n",
3324                     numberTest,iDo);
3325#endif
3326            if (!choiceObject) {
3327              delete choice.possibleBranch;
3328              choice.possibleBranch=NULL;
3329            }
3330            break; // give up anyway
3331          }
3332        }
3333        if (model->messageHandler()->logLevel()>3&&numberBeforeTrust&&dynamicObject) 
3334          dynamicObject->print(1,choice.possibleBranch->value());
3335        if (skipAll<0)
3336          canSkip=true;
3337        if (!canSkip) {
3338          if (!doneHotStart) {
3339            // Mark hot start
3340            doneHotStart=true;
3341            solver->markHotStart();
3342            xMark++;
3343          }
3344          numberTest--;
3345          // just do a few
3346          if (searchStrategy==2)
3347            solver->setIntParam(OsiMaxNumIterationHotStart,10); 
3348          double objectiveChange ;
3349          double newObjectiveValue=1.0e100;
3350          int j;
3351          // status is 0 finished, 1 infeasible and other
3352          int iStatus;
3353          /*
3354            Try the down direction first. (Specify the initial branching alternative as
3355            down with a call to way(-1). Each subsequent call to branch() performs the
3356            specified branch and advances the branch object state to the next branch
3357            alternative.)
3358          */
3359          choice.possibleBranch->way(-1) ;
3360          choice.possibleBranch->branch() ;
3361          solver->solveFromHotStart() ;
3362          bool needHotStartUpdate=false;
3363          numberStrongDone++;
3364          numberStrongIterations += solver->getIterationCount();
3365          /*
3366            We now have an estimate of objective degradation that we can use for strong
3367            branching. If we're over the cutoff, the variable is monotone up.
3368            If we actually made it to optimality, check for a solution, and if we have
3369            a good one, call setBestSolution to process it. Note that this may reduce the
3370            cutoff, so we check again to see if we can declare this variable monotone.
3371          */
3372          if (solver->isProvenOptimal())
3373            iStatus=0; // optimal
3374          else if (solver->isIterationLimitReached()
3375                   &&!solver->isDualObjectiveLimitReached())
3376            iStatus=2; // unknown
3377          else
3378            iStatus=1; // infeasible
3379          if (iStatus!=2&&solver->getIterationCount()>
3380              realMaxHotIterations)
3381            numberUnfinished++;
3382          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3383          choice.numItersDown = solver->getIterationCount();
3384          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3385          // Update branching information if wanted
3386          CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3387          if (cbcobj) {
3388            CbcObject * object = cbcobj->object();
3389            assert (object) ;
3390            CbcObjectUpdateData update = object->createUpdateInformation(solver,this,cbcobj);
3391            update.objectNumber_ = choice.objectNumber;
3392            model->addUpdateInformation(update);
3393          } else {
3394            decision->updateInformation( solver,this);
3395          }
3396          if (!iStatus) {
3397            choice.finishedDown = true ;
3398            if (newObjectiveValue>=cutoff) {
3399              objectiveChange = 1.0e100; // say infeasible
3400              numberStrongInfeasible++;
3401            } else {
3402              // See if integer solution
3403              if (model->feasibleSolution(choice.numIntInfeasDown,
3404                                          choice.numObjInfeasDown)
3405                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
3406                if (auxiliaryInfo->solutionAddsCuts()) {
3407                  needHotStartUpdate=true;
3408                  solver->unmarkHotStart();
3409                }
3410                model->setBestSolution(CBC_STRONGSOL,
3411                                       newObjectiveValue,
3412                                       solver->getColSolution()) ;
3413                if (needHotStartUpdate) {
3414                  model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3415                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3416                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3417                  model->feasibleSolution(choice.numIntInfeasDown,
3418                                          choice.numObjInfeasDown);
3419                }
3420                model->setLastHeuristic(NULL);
3421                model->incrementUsed(solver->getColSolution());
3422                cutoff =model->getCutoff();
3423                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3424                  objectiveChange = 1.0e100 ;
3425                  numberStrongInfeasible++;
3426                }
3427              }
3428            }
3429          } else if (iStatus==1) {
3430            objectiveChange = 1.0e100 ;
3431            numberStrongInfeasible++;
3432          } else {
3433            // Can't say much as we did not finish
3434            choice.finishedDown = false ;
3435            numberUnfinished++;
3436          }
3437          choice.downMovement = objectiveChange ;
3438         
3439          // restore bounds
3440          for ( j=0;j<numberColumns;j++) {
3441            if (saveLower[j] != lower[j])
3442              solver->setColLower(j,saveLower[j]);
3443            if (saveUpper[j] != upper[j])
3444              solver->setColUpper(j,saveUpper[j]);
3445          }
3446          if(needHotStartUpdate) {
3447            needHotStartUpdate = false;
3448            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3449            //we may again have an integer feasible solution
3450            int numberIntegerInfeasibilities;
3451            int numberObjectInfeasibilities;
3452            if (model->feasibleSolution(
3453                                        numberIntegerInfeasibilities,
3454                                        numberObjectInfeasibilities)) {
3455#ifdef BONMIN
3456              //In this case node has become integer feasible, let us exit the loop
3457              std::cout<<"Node has become integer feasible"<<std::endl;
3458              numberUnsatisfied_ = 0;
3459              break;
3460#endif
3461              double objValue = solver->getObjValue();
3462              model->setBestSolution(CBC_STRONGSOL,
3463                                     objValue,
3464                                     solver->getColSolution()) ;
3465              model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3466              cutoff =model->getCutoff();
3467            }
3468            solver->markHotStart();
3469          }
3470#ifdef DO_ALL_AT_ROOT
3471          if (!numberNodes)
3472            printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3473                   choice.objectNumber,iStatus,newObjectiveValue,choice.numItersDown,
3474                   choice.downMovement,choice.finishedDown,choice.numIntInfeasDown,
3475                   choice.numObjInfeasDown);
3476#endif
3477         
3478          // repeat the whole exercise, forcing the variable up
3479          choice.possibleBranch->branch();
3480          solver->solveFromHotStart() ;
3481          numberStrongDone++;
3482          numberStrongIterations += solver->getIterationCount();
3483          /*
3484            We now have an estimate of objective degradation that we can use for strong
3485            branching. If we're over the cutoff, the variable is monotone up.
3486            If we actually made it to optimality, check for a solution, and if we have
3487            a good one, call setBestSolution to process it. Note that this may reduce the
3488            cutoff, so we check again to see if we can declare this variable monotone.
3489          */
3490          if (solver->isProvenOptimal())
3491            iStatus=0; // optimal
3492          else if (solver->isIterationLimitReached()
3493                   &&!solver->isDualObjectiveLimitReached())
3494            iStatus=2; // unknown
3495          else
3496            iStatus=1; // infeasible
3497          if (iStatus!=2&&solver->getIterationCount()>
3498              realMaxHotIterations)
3499            numberUnfinished++;
3500          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3501          choice.numItersUp = solver->getIterationCount();
3502          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3503          // Update branching information if wanted
3504          cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3505          if (cbcobj) {
3506            CbcObject * object = cbcobj->object();
3507            assert (object) ;
3508            CbcObjectUpdateData update = object->createUpdateInformation(solver,this,cbcobj);
3509            update.objectNumber_ = choice.objectNumber;
3510            model->addUpdateInformation(update);
3511          } else {
3512            decision->updateInformation( solver,this);
3513          }
3514          if (!iStatus) {
3515            choice.finishedUp = true ;
3516            if (newObjectiveValue>=cutoff) {
3517              objectiveChange = 1.0e100; // say infeasible
3518              numberStrongInfeasible++;
3519            } else {
3520              // See if integer solution
3521              if (model->feasibleSolution(choice.numIntInfeasUp,
3522                                          choice.numObjInfeasUp)
3523                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
3524#ifdef BONMIN
3525                std::cout<<"Node has become integer feasible"<<std::endl;
3526                numberUnsatisfied_ = 0;
3527                break;
3528#endif
3529                if (auxiliaryInfo->solutionAddsCuts()) {
3530                  needHotStartUpdate=true;
3531                  solver->unmarkHotStart();
3532                }
3533                model->setBestSolution(CBC_STRONGSOL,
3534                                       newObjectiveValue,
3535                                       solver->getColSolution()) ;
3536                if (needHotStartUpdate) {
3537                  model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3538                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3539                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
3540                  model->feasibleSolution(choice.numIntInfeasDown,
3541                                          choice.numObjInfeasDown);
3542                }
3543                model->setLastHeuristic(NULL);
3544                model->incrementUsed(solver->getColSolution());
3545                cutoff =model->getCutoff();
3546                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3547                  objectiveChange = 1.0e100 ;
3548                  numberStrongInfeasible++;
3549                }
3550              }
3551            }
3552          } else if (iStatus==1) {
3553            objectiveChange = 1.0e100 ;
3554            numberStrongInfeasible++;
3555          } else {
3556            // Can't say much as we did not finish
3557            choice.finishedUp = false ;
3558            numberUnfinished++;
3559          }
3560          choice.upMovement = objectiveChange ;
3561         
3562          // restore bounds
3563          for ( j=0;j<numberColumns;j++) {
3564            if (saveLower[j] != lower[j])
3565              solver->setColLower(j,saveLower[j]);
3566            if (saveUpper[j] != upper[j])
3567              solver->setColUpper(j,saveUpper[j]);
3568          }
3569          if(needHotStartUpdate) {
3570            needHotStartUpdate = false;
3571            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3572            //we may again have an integer feasible solution
3573            int numberIntegerInfeasibilities;
3574            int numberObjectInfeasibilities;
3575            if (model->feasibleSolution(
3576                                        numberIntegerInfeasibilities,
3577                                        numberObjectInfeasibilities)) {
3578              double objValue = solver->getObjValue();
3579              model->setBestSolution(CBC_STRONGSOL,
3580                                     objValue,
3581                                     solver->getColSolution()) ;
3582              model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3583              cutoff =model->getCutoff();
3584            }
3585            solver->markHotStart();
3586          }
3587         
3588#ifdef DO_ALL_AT_ROOT
3589          if (!numberNodes)
3590            printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3591                   choice.objectNumber,iStatus,newObjectiveValue,choice.numItersUp,
3592                   choice.upMovement,choice.finishedUp,choice.numIntInfeasUp,
3593                   choice.numObjInfeasUp);
3594#endif
3595        }
3596       
3597        solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit2); 
3598        /*
3599          End of evaluation for this candidate variable. Possibilities are:
3600          * Both sides below cutoff; this variable is a candidate for branching.
3601          * Both sides infeasible or above the objective cutoff: no further action
3602          here. Break from the evaluation loop and assume the node will be purged
3603          by the caller.
3604          * One side below cutoff: Install the branch (i.e., fix the variable). Break
3605          from the evaluation loop and assume the node will be reoptimised by the
3606          caller.
3607        */
3608        // reset
3609        choice.possibleBranch->resetNumberBranchesLeft();
3610        if (choice.upMovement<1.0e100) {
3611          if(choice.downMovement<1.0e100) {
3612            // In case solution coming in was odd
3613            choice.upMovement = CoinMax(0.0,choice.upMovement);
3614            choice.downMovement = CoinMax(0.0,choice.downMovement);
3615#if ZERO_ONE==2
3616            // branch on 0-1 first (temp)
3617            if (fabs(choice.possibleBranch->value())<1.0) {
3618              choice.upMovement *= ZERO_FAKE;
3619              choice.downMovement *= ZERO_FAKE;
3620            }
3621#endif
3622            // feasible - see which best
3623            if (!canSkip) {
3624              if (model->messageHandler()->logLevel()>3) 
3625                printf("sort %g downest %g upest %g ",sort[iDo],downEstimate[iObject],
3626                       upEstimate[iObject]);
3627              model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
3628                << iObject << iColumn
3629                <<choice.downMovement<<choice.numIntInfeasDown
3630                <<choice.upMovement<<choice.numIntInfeasUp
3631                <<choice.possibleBranch->value()
3632                <<CoinMessageEol;
3633            }
3634            int betterWay;
3635            {
3636              CbcBranchingObject * branchObj =
3637                dynamic_cast <CbcBranchingObject *>(branch_) ;
3638              if (branch_)
3639                assert (branchObj);
3640              betterWay = decision->betterBranch(choice.possibleBranch,
3641                                                 branchObj,
3642                                                 choice.upMovement,
3643                                                 choice.numIntInfeasUp ,
3644                                                 choice.downMovement,
3645                                                 choice.numIntInfeasDown );
3646            }
3647            if (betterWay) {
3648              // C) create branching object
3649              if (choiceObject) {
3650                delete branch_;
3651                branch_ = choice.possibleBranch->clone();
3652              } else {
3653                delete branch_;
3654                branch_ = choice.possibleBranch;
3655                choice.possibleBranch=NULL;
3656              }
3657              {
3658                CbcBranchingObject * branchObj =
3659                  dynamic_cast <CbcBranchingObject *>(branch_) ;
3660                assert (branchObj);
3661                //branchObj->way(preferredWay);
3662                branchObj->way(betterWay);
3663              }
3664              bestChoice = choice.objectNumber;
3665              whichChoice = iDo;
3666              if (numberStrong<=1) {
3667                delete ws;
3668                ws=NULL;
3669                break;
3670              }
3671            } else {
3672              if (!choiceObject) {
3673                delete choice.possibleBranch;
3674                choice.possibleBranch=NULL;
3675              }
3676              if (iDo>=2*numberStrong) {
3677                delete ws;
3678                ws=NULL;
3679                break;
3680              }
3681              if (!dynamicObject||dynamicObject->numberTimesUp()>1) {
3682                if (iDo-whichChoice>=numberStrong) {
3683                  if (!choiceObject) {
3684                    delete choice.possibleBranch;
3685                    choice.possibleBranch=NULL;
3686                  }
3687                  break; // give up
3688                }
3689              } else {
3690                if (iDo-whichChoice>=2*numberStrong) {
3691                  delete ws;
3692                  ws=NULL;
3693                  if (!choiceObject) {
3694                    delete choice.possibleBranch;
3695                    choice.possibleBranch=NULL;
3696                  }
3697                  break; // give up
3698                }
3699              }
3700            }
3701          } else {
3702            // up feasible, down infeasible
3703            anyAction=-1;
3704            worstFeasible = CoinMax(worstFeasible,choice.upMovement);
3705            model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
3706              << iObject << iColumn
3707              <<choice.downMovement<<choice.numIntInfeasDown
3708              <<choice.upMovement<<choice.numIntInfeasUp
3709              <<choice.possibleBranch->value()
3710              <<CoinMessageEol;
3711            //printf("Down infeasible for choice %d sequence %d\n",i,
3712            // model->object(choice.objectNumber)->columnNumber());
3713            choice.fix=1;
3714            numberToFix++;
3715            choice.possibleBranch->fix(solver,saveLower,saveUpper,1);
3716            if (!choiceObject) {
3717              choice.possibleBranch=NULL;
3718            } else {
3719              //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3720              choice.possibleBranch=choiceObject;
3721            }
3722            assert(doneHotStart);
3723            solver->unmarkHotStart();
3724            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3725            solver->markHotStart();
3726            // may be infeasible (if other way stopped on iterations)
3727            if (!solver->isProvenOptimal()) {
3728              // neither side feasible
3729              anyAction=-2;
3730              if (!choiceObject) {
3731                delete choice.possibleBranch;
3732                choice.possibleBranch=NULL;
3733              }
3734              //printf("Both infeasible for choice %d sequence %d\n",i,
3735              // model->object(choice.objectNumber)->columnNumber());
3736              delete ws;
3737              ws=NULL;
3738              break;
3739            }
3740          }
3741        } else {
3742          if(choice.downMovement<1.0e100) {
3743            // down feasible, up infeasible
3744            anyAction=-1;
3745            worstFeasible = CoinMax(worstFeasible,choice.downMovement);
3746            model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
3747              << iObject << iColumn
3748              <<choice.downMovement<<choice.numIntInfeasDown
3749              <<choice.upMovement<<choice.numIntInfeasUp
3750              <<choice.possibleBranch->value()
3751              <<CoinMessageEol;
3752            choice.fix=-1;
3753            numberToFix++;
3754            choice.possibleBranch->fix(solver,saveLower,saveUpper,-1);
3755            if (!choiceObject) {
3756              choice.possibleBranch=NULL;
3757            } else {
3758              //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3759              choice.possibleBranch=choiceObject;
3760            }
3761            assert(doneHotStart);
3762            solver->unmarkHotStart();
3763            model->resolve(NULL,11,saveSolution,saveLower,saveUpper);
3764            solver->markHotStart();
3765            // may be infeasible (if other way stopped on iterations)
3766            if (!solver->isProvenOptimal()) {
3767              // neither side feasible
3768              anyAction=-2;
3769              if (!choiceObject) {
3770                delete choice.possibleBranch;
3771                choice.possibleBranch=NULL;
3772              }
3773              delete ws;
3774              ws=NULL;
3775              break;
3776            }
3777          } else {
3778            // neither side feasible
3779            anyAction=-2;
3780            if (!choiceObject) {
3781              delete choice.possibleBranch;
3782              choice.possibleBranch=NULL;
3783            }
3784            delete ws;
3785            ws=NULL;
3786            break;
3787          }
3788        }
3789        // Check max time
3790        hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
3791                       model->getDblParam(CbcModel::CbcMaximumSeconds));
3792        if (hitMaxTime) {
3793          // make sure rest are fast
3794          doQuickly=true;
3795          for ( int jDo=iDo+1;jDo<numberToDo;jDo++) {
3796            int iObject = whichObject[iDo];
3797            OsiObject * object = model->modifiableObject(iObject);
3798            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3799              dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3800            if (dynamicObject) 
3801              dynamicObject->setNumberBeforeTrust(0);
3802          }
3803          numberTest=0;
3804        }
3805        if (!choiceObject) {
3806          delete choice.possibleBranch;
3807        }
3808      }
3809      if (model->messageHandler()->logLevel()>3) { 
3810        if (anyAction==-2) {
3811          printf("infeasible\n");
3812        } else if(anyAction==-1) {
3813          printf("%d fixed AND choosing %d iDo %d iChosenWhen %d numberToDo %d\n",numberToFix,bestChoice,
3814                   iDo,whichChoice,numberToDo);
3815        } else {
3816          int iObject = whichObject[whichChoice];
3817          OsiObject * object = model->modifiableObject(iObject);
3818          CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3819            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3820          if (dynamicObject) {
3821            int iColumn = dynamicObject->columnNumber();
3822            printf("choosing %d (column %d) iChosenWhen %d numberToDo %d\n",bestChoice,
3823                   iColumn,whichChoice,numberToDo);
3824          }
3825        }
3826      }
3827      if (doneHotStart) {
3828        // Delete the snapshot
3829        solver->unmarkHotStart();
3830        // back to normal
3831        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3832        // restore basis
3833        solver->setWarmStart(ws);
3834      }
3835      solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
3836      // Unless infeasible we will carry on
3837      // But we could fix anyway
3838      if (numberToFix&&!hitMaxTime) {
3839        if (anyAction!=-2) {
3840          // apply and take off
3841          bool feasible=true;
3842          // can do quick optimality check
3843          int easy=2;
3844          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
3845          model->resolve(NULL,11,saveSolution,saveLower,saveUpper) ;
3846          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3847          feasible = solver->isProvenOptimal();
3848          if (feasible) {
3849            anyAction=0;
3850            // See if candidate still possible
3851            if (branch_) {
3852              const OsiObject * object = model->object(bestChoice);
3853              int preferredWay;
3854              double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
3855              if (!infeasibility) {
3856                // take out
3857                delete branch_;
3858                branch_=NULL;
3859              } else {
3860                CbcBranchingObject * branchObj =
3861                  dynamic_cast <CbcBranchingObject *>(branch_) ;
3862                assert (branchObj);
3863                branchObj->way(preferredWay);
3864              }
3865            }
3866          } else {
3867            anyAction=-2;
3868            finished=true;
3869          }
3870        }
3871      }
3872      // If  fixed then round again
3873      if (!branch_&&anyAction!=-2) {
3874        finished=false;
3875      }
3876      delete ws;
3877    }
3878  }
3879  // update number of strong iterations etc
3880  model->incrementStrongInfo(numberStrongDone,numberStrongIterations,
3881                             anyAction==-2 ? 0:numberToFix,anyAction==-2);
3882  if (model->searchStrategy()==-1) {
3883#ifndef COIN_DEVELOP
3884    if (solver->messageHandler()->logLevel()>1)
3885#endif
3886      printf("%d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3887             numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3888             numberNotTrusted);
3889    // decide what to do
3890    int strategy=1;
3891    if (((numberUnfinished*4>numberStrongDone&&
3892          numberStrongInfeasible*40<numberStrongDone)||
3893         numberStrongInfeasible<0)&&model->numberStrong()<10&&model->numberBeforeTrust()<=20&&model->numberObjects()>CoinMax(1000,solver->getNumRows())) {
3894      strategy=2;
3895#ifdef COIN_DEVELOP
3896      //if (model->logLevel()>1)
3897      printf("going to strategy 2\n");
3898#endif
3899      // Weaken
3900      model->setNumberStrong(2);
3901      model->setNumberBeforeTrust(1);
3902      model->synchronizeNumberBeforeTrust();
3903    }
3904    if (numberNodes)
3905      strategy=1;  // should only happen after hot start
3906    model->setSearchStrategy(strategy);
3907  } else if (numberStrongDone) {
3908    //printf("%d strongB, %d iters, %d inf, %d not finished, %d not trusted\n",
3909    //   numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3910    //   numberNotTrusted);
3911  }
3912  if (model->searchStrategy()==1&&numberNodes>500&&numberNodes<-510) {
3913#ifndef COIN_DEVELOP
3914    if (solver->messageHandler()->logLevel()>1)
3915#endif
3916      printf("after %d nodes - %d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3917             numberNodes,numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3918             numberNotTrusted);
3919    // decide what to do
3920    if (numberUnfinished*10>numberStrongDone+1||
3921        !numberStrongInfeasible) {
3922      printf("going to strategy 2\n");
3923      // Weaken
3924      model->setNumberStrong(2);
3925      model->setNumberBeforeTrust(1);
3926      model->synchronizeNumberBeforeTrust();
3927      model->setSearchStrategy(2);
3928    }
3929  }
3930  // Set guessed solution value
3931  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
3932 
3933  /*
3934    Cleanup, then we're finished
3935  */
3936  if (!model->branchingMethod())
3937    delete decision;
3938 
3939  delete choiceObject;
3940  delete [] sort;
3941  delete [] whichObject;
3942#ifdef RANGING
3943  delete [] objectMark;
3944#endif
3945  delete [] saveLower;
3946  delete [] saveUpper;
3947  delete [] upEstimate;
3948  delete [] downEstimate;
3949# ifdef COIN_HAS_CLP
3950  if (osiclp) {
3951    osiclp->setSpecialOptions(saveClpOptions);
3952  }
3953# endif
3954  // restore solution
3955  solver->setColSolution(saveSolution);
3956  model->reserveCurrentSolution(saveSolution);
3957  delete [] saveSolution;
3958  model->setStateOfSearch(saveStateOfSearch);
3959  model->setLogLevel(saveLogLevel);
3960  // delete extra regions
3961  if (usefulInfo.usefulRegion_) {
3962    delete [] usefulInfo.usefulRegion_;
3963    delete [] usefulInfo.indexRegion_;
3964    delete [] usefulInfo.pi_;
3965    usefulInfo.usefulRegion_ = NULL;
3966    usefulInfo.indexRegion_ = NULL;
3967    usefulInfo.pi_ = NULL;
3968  }
3969  return anyAction;
3970}
3971int CbcNode::analyze (CbcModel *model, double * results)
3972{
3973  int i;
3974  int numberIterationsAllowed = model->numberAnalyzeIterations();
3975  OsiSolverInterface * solver = model->solver();
3976  objectiveValue_ = solver->getObjSense()*solver->getObjValue();
3977  double cutoff =model->getCutoff();
3978  const double * lower = solver->getColLower();
3979  const double * upper = solver->getColUpper();
3980  const double * dj = solver->getReducedCost();
3981  int numberObjects = model->numberObjects();
3982  int numberColumns = model->getNumCols();
3983  // Initialize arrays
3984  int numberIntegers = model->numberIntegers();
3985  int * back = new int[numberColumns];
3986  const int * integerVariable = model->integerVariable();
3987  for (i=0;i<numberColumns;i++) 
3988    back[i]=-1;
3989  // What results is
3990  double * newLower = results;
3991  double * objLower = newLower+numberIntegers;
3992  double * newUpper = objLower+numberIntegers;
3993  double * objUpper = newUpper+numberIntegers;
3994  for (i=0;i<numberIntegers;i++) {
3995    int iColumn = integerVariable[i];
3996    back[iColumn]=i;
3997    newLower[i]=0.0;
3998    objLower[i]=-COIN_DBL_MAX;
3999    newUpper[i]=0.0;
4000    objUpper[i]=-COIN_DBL_MAX;
4001  }
4002  double * saveUpper = new double[numberColumns];
4003  double * saveLower = new double[numberColumns];
4004  int anyAction=0;
4005  // Save solution in case heuristics need good solution later
4006 
4007  double * saveSolution = new double[numberColumns];
4008  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
4009  model->reserveCurrentSolution(saveSolution);
4010  for (i=0;i<numberColumns;i++) {
4011    saveLower[i] = lower[i];
4012    saveUpper[i] = upper[i];
4013  }
4014  // Get arrays to sort
4015  double * sort = new double[numberObjects];
4016  int * whichObject = new int[numberObjects];
4017  int numberToFix=0;
4018  int numberToDo=0;
4019  double integerTolerance = 
4020    model->getDblParam(CbcModel::CbcIntegerTolerance);
4021  // point to useful information
4022  OsiBranchingInformation usefulInfo = model->usefulInformation();
4023  // and modify
4024  usefulInfo.depth_=depth_;
4025 
4026  // compute current state
4027  int numberObjectInfeasibilities; // just odd ones
4028  int numberIntegerInfeasibilities;
4029  model->feasibleSolution(
4030                          numberIntegerInfeasibilities,
4031                          numberObjectInfeasibilities);
4032# ifdef COIN_HAS_CLP
4033  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
4034  int saveClpOptions=0;
4035  bool fastIterations = (model->specialOptions()&8)!=0;
4036  if (osiclp&&fastIterations) {
4037    // for faster hot start
4038    saveClpOptions = osiclp->specialOptions();
4039    osiclp->setSpecialOptions(saveClpOptions|8192);
4040  }
4041# else
4042  bool fastIterations = false ;
4043# endif
4044  /*
4045    Scan for branching objects that indicate infeasibility. Choose candidates
4046    using priority as the first criteria, then integer infeasibility.
4047   
4048    The algorithm is to fill the array with a set of good candidates (by
4049    infeasibility) with priority bestPriority.  Finding a candidate with
4050    priority better (less) than bestPriority flushes the choice array. (This
4051    serves as initialization when the first candidate is found.)
4052   
4053  */
4054  numberToDo=0;
4055  for (i=0;i<numberObjects;i++) {
4056    OsiObject * object = model->modifiableObject(i);
4057    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4058      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4059    if(!dynamicObject)
4060      continue;
4061    int preferredWay;
4062    double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
4063    int iColumn = dynamicObject->columnNumber();
4064    if (saveUpper[iColumn]==saveLower[iColumn])
4065      continue;
4066    if (infeasibility)
4067      sort[numberToDo]=-1.0e10-infeasibility;
4068    else
4069      sort[numberToDo]=-fabs(dj[iColumn]);
4070    whichObject[numberToDo++]=i;
4071  }
4072  // Save basis
4073  CoinWarmStart * ws = solver->getWarmStart();
4074  int saveLimit;
4075  solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
4076  int targetIterations = CoinMax(500,numberIterationsAllowed/numberObjects);
4077  if (saveLimit<targetIterations)
4078    solver->setIntParam(OsiMaxNumIterationHotStart,targetIterations); 
4079  // Mark hot start
4080  solver->markHotStart();
4081  // Sort
4082  CoinSort_2(sort,sort+numberToDo,whichObject);
4083  double * currentSolution = model->currentSolution();
4084  double objMin = 1.0e50;
4085  double objMax = -1.0e50;
4086  bool needResolve=false;
4087  int iDo;
4088  for (iDo=0;iDo<numberToDo;iDo++) {
4089    CbcStrongInfo choice;
4090    int iObject = whichObject[iDo];
4091    OsiObject * object = model->modifiableObject(iObject);
4092    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4093      static_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4094    int iColumn = dynamicObject->columnNumber();
4095    int preferredWay;
4096    object->infeasibility(&usefulInfo,preferredWay);
4097    double value = currentSolution[iColumn];
4098    double nearest = floor(value+0.5);
4099    double lowerValue = floor(value);
4100    bool satisfied=false;
4101    if (fabs(value-nearest)<=integerTolerance||value<saveLower[iColumn]||value>saveUpper[iColumn]) {
4102      satisfied=true;
4103      double newValue;
4104      if (nearest<saveUpper[iColumn]) {
4105        newValue = nearest + 1.0001*integerTolerance;
4106        lowerValue = nearest;
4107      } else {
4108        newValue = nearest - 1.0001*integerTolerance;
4109        lowerValue = nearest-1;
4110      }
4111      currentSolution[iColumn]=newValue;
4112    }
4113    double upperValue = lowerValue+1.0;
4114    //CbcSimpleInteger * obj =
4115    //dynamic_cast <CbcSimpleInteger *>(object) ;
4116    //if (obj) {
4117    //choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
4118    //} else {
4119      CbcObject * obj =
4120        dynamic_cast <CbcObject *>(object) ;
4121      assert (obj);
4122      choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
4123      //}
4124    currentSolution[iColumn]=value;
4125    // Save which object it was
4126    choice.objectNumber=iObject;
4127    choice.numIntInfeasUp = numberUnsatisfied_;
4128    choice.numIntInfeasDown = numberUnsatisfied_;
4129    choice.downMovement = 0.0;
4130    choice.upMovement = 0.0;
4131    choice.numItersDown = 0;
4132    choice.numItersUp = 0;
4133    choice.fix=0; // say not fixed
4134    double objectiveChange ;
4135    double newObjectiveValue=1.0e100;
4136    int j;
4137    // status is 0 finished, 1 infeasible and other
4138    int iStatus;
4139    /*
4140      Try the down direction first. (Specify the initial branching alternative as
4141      down with a call to way(-1). Each subsequent call to branch() performs the
4142      specified branch and advances the branch object state to the next branch
4143      alternative.)
4144    */
4145    choice.possibleBranch->way(-1) ;
4146    choice.possibleBranch->branch() ;
4147    if (fabs(value-lowerValue)>integerTolerance) {
4148      solver->solveFromHotStart() ;
4149      /*
4150        We now have an estimate of objective degradation that we can use for strong
4151        branching. If we're over the cutoff, the variable is monotone up.
4152        If we actually made it to optimality, check for a solution, and if we have
4153        a good one, call setBestSolution to process it. Note that this may reduce the
4154        cutoff, so we check again to see if we can declare this variable monotone.
4155      */
4156      if (solver->isProvenOptimal())
4157        iStatus=0; // optimal
4158      else if (solver->isIterationLimitReached()
4159               &&!solver->isDualObjectiveLimitReached())
4160        iStatus=2; // unknown
4161      else
4162        iStatus=1; // infeasible
4163      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
4164      choice.numItersDown = solver->getIterationCount();
4165      numberIterationsAllowed -= choice.numItersDown;
4166      objectiveChange = newObjectiveValue  - objectiveValue_;
4167      if (!iStatus) {
4168        choice.finishedDown = true ;
4169        if (newObjectiveValue>=cutoff) {
4170          objectiveChange = 1.0e100; // say infeasible
4171        } else {
4172          // See if integer solution
4173          if (model->feasibleSolution(choice.numIntInfeasDown,
4174                                      choice.numObjInfeasDown)
4175              &&model->problemFeasibility()->feasible(model,-1)>=0) {
4176            model->setBestSolution(CBC_STRONGSOL,
4177                                   newObjectiveValue,
4178                                   solver->getColSolution()) ;
4179            model->setLastHeuristic(NULL);
4180            model->incrementUsed(solver->getColSolution());
4181            cutoff =model->getCutoff();
4182            if (newObjectiveValue >= cutoff)    //  *new* cutoff
4183              objectiveChange = 1.0e100 ;
4184          }
4185        }
4186      } else if (iStatus==1) {
4187        objectiveChange = 1.0e100 ;
4188      } else {
4189        // Can't say much as we did not finish
4190        choice.finishedDown = false ;
4191      }
4192      choice.downMovement = objectiveChange ;
4193    }
4194    // restore bounds
4195    for ( j=0;j<numberColumns;j++) {
4196      if (saveLower[j] != lower[j])
4197        solver->setColLower(j,saveLower[j]);
4198      if (saveUpper[j] != upper[j])
4199        solver->setColUpper(j,saveUpper[j]);
4200    }
4201    // repeat the whole exercise, forcing the variable up
4202    choice.possibleBranch->branch();
4203    if (fabs(value-upperValue)>integerTolerance) {
4204      solver->solveFromHotStart() ;
4205      /*
4206        We now have an estimate of objective degradation that we can use for strong
4207        branching. If we're over the cutoff, the variable is monotone up.
4208        If we actually made it to optimality, check for a solution, and if we have
4209        a good one, call setBestSolution to process it. Note that this may reduce the
4210        cutoff, so we check again to see if we can declare this variable monotone.
4211      */
4212      if (solver->isProvenOptimal())
4213        iStatus=0; // optimal
4214      else if (solver->isIterationLimitReached()
4215               &&!solver->isDualObjectiveLimitReached())
4216        iStatus=2; // unknown
4217      else
4218        iStatus=1; // infeasible
4219      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
4220      choice.numItersUp = solver->getIterationCount();
4221      numberIterationsAllowed -= choice.numItersUp;
4222      objectiveChange = newObjectiveValue  - objectiveValue_;
4223      if (!iStatus) {
4224        choice.finishedUp = true ;
4225        if (newObjectiveValue>=cutoff) {
4226          objectiveChange = 1.0e100; // say infeasible
4227        } else {
4228          // See if integer solution
4229          if (model->feasibleSolution(choice.numIntInfeasUp,
4230                                      choice.numObjInfeasUp)
4231              &&model->problemFeasibility()->feasible(model,-1)>=0) {
4232            model->setBestSolution(CBC_STRONGSOL,
4233                                   newObjectiveValue,
4234                                   solver->getColSolution()) ;
4235            model->setLastHeuristic(NULL);
4236            model->incrementUsed(solver->getColSolution());
4237            cutoff =model->getCutoff();
4238            if (newObjectiveValue >= cutoff)    //  *new* cutoff
4239              objectiveChange = 1.0e100 ;
4240          }
4241        }
4242      } else if (iStatus==1) {
4243        objectiveChange = 1.0e100 ;
4244      } else {
4245        // Can't say much as we did not finish
4246        choice.finishedUp = false ;
4247      }
4248      choice.upMovement = objectiveChange ;
4249     
4250      // restore bounds
4251      for ( j=0;j<numberColumns;j++) {
4252        if (saveLower[j] != lower[j])
4253          solver->setColLower(j,saveLower[j]);
4254        if (saveUpper[j] != upper[j])
4255          solver->setColUpper(j,saveUpper[j]);
4256      }
4257    }
4258    // If objective goes above certain amount we can set bound
4259    int jInt = back[iColumn];
4260    newLower[jInt]=upperValue;
4261    if (choice.finishedDown)
4262      objLower[jInt]=choice.downMovement+objectiveValue_;
4263    else
4264      objLower[jInt]=objectiveValue_;
4265    newUpper[jInt]=lowerValue;
4266    if (choice.finishedUp)
4267      objUpper[jInt]=choice.upMovement+objectiveValue_;
4268    else
4269      objUpper[jInt]=objectiveValue_;
4270    objMin = CoinMin(CoinMin(objLower[jInt],objUpper[jInt]),objMin);
4271    /*
4272      End of evaluation for this candidate variable. Possibilities are:
4273      * Both sides below cutoff; this variable is a candidate for branching.
4274      * Both sides infeasible or above the objective cutoff: no further action
4275      here. Break from the evaluation loop and assume the node will be purged
4276      by the caller.
4277      * One side below cutoff: Install the branch (i.e., fix the variable). Break
4278      from the evaluation loop and assume the node will be reoptimised by the
4279      caller.
4280    */
4281    if (choice.upMovement<1.0e100) {
4282      if(choice.downMovement<1.0e100) {
4283        objMax = CoinMax(CoinMax(objLower[jInt],objUpper[jInt]),objMax);
4284        // In case solution coming in was odd
4285        choice.upMovement = CoinMax(0.0,choice.upMovement);
4286        choice.downMovement = CoinMax(0.0,choice.downMovement);
4287        // feasible -
4288        model->messageHandler()->message(CBC_STRONG,*model->messagesPointer())
4289          << iObject << iColumn
4290          <<choice.downMovement<<choice.numIntInfeasDown
4291          <<choice.upMovement<<choice.numIntInfeasUp
4292          <<value
4293          <<CoinMessageEol;
4294      } else {
4295        // up feasible, down infeasible
4296        anyAction=-1;
4297        if (!satisfied)
4298          needResolve=true;
4299        choice.fix=1;
4300        numberToFix++;
4301        saveLower[iColumn]=upperValue;
4302        solver->setColLower(iColumn,upperValue);
4303      }
4304    } else {
4305      if(choice.downMovement<1.0e100) {
4306        // down feasible, up infeasible
4307        anyAction=-1;
4308        if (!satisfied)
4309          needResolve=true;
4310        choice.fix=-1;
4311        numberToFix++;
4312        saveUpper[iColumn]=lowerValue;
4313        solver->setColUpper(iColumn,lowerValue);
4314      } else {
4315        // neither side feasible
4316        anyAction=-2;
4317        printf("Both infeasible for choice %d sequence %d\n",i,
4318               model->object(choice.objectNumber)->columnNumber());
4319        delete ws;
4320        ws=NULL;
4321        //solver->writeMps("bad");
4322        numberToFix=-1;
4323        delete choice.possibleBranch;
4324        choice.possibleBranch=NULL;
4325        break;
4326      }
4327    }
4328    delete choice.possibleBranch;
4329    if (numberIterationsAllowed<=0)
4330      break;
4331    //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
4332    //     choice.downMovement,choice.upMovement,value);
4333  }
4334  printf("Best possible solution %g, can fix more if solution of %g found - looked at %d variables in %d iterations\n",
4335         objMin,objMax,iDo,model->numberAnalyzeIterations()-numberIterationsAllowed);
4336  model->setNumberAnalyzeIterations(numberIterationsAllowed);
4337  // Delete the snapshot
4338  solver->unmarkHotStart();
4339  // back to normal
4340  solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
4341  solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
4342  // restore basis
4343  solver->setWarmStart(ws);
4344  delete ws;
4345 
4346  delete [] sort;
4347  delete [] whichObject;
4348  delete [] saveLower;
4349  delete [] saveUpper;
4350  delete [] back;
4351  // restore solution
4352  solver->setColSolution(saveSolution);
4353# ifdef COIN_HAS_CLP
4354  if (osiclp) 
4355    osiclp->setSpecialOptions(saveClpOptions);
4356# endif
4357  model->reserveCurrentSolution(saveSolution);
4358  delete [] saveSolution;
4359  if (needResolve)
4360    solver->resolve();
4361  return numberToFix;
4362}
4363
4364
4365CbcNode::CbcNode(const CbcNode & rhs)
4366  : CoinTreeNode(rhs) 
4367{ 
4368#ifdef CHECK_NODE
4369  printf("CbcNode %x Constructor from rhs %x\n",this,&rhs);
4370#endif
4371  if (rhs.nodeInfo_)
4372    nodeInfo_ = rhs.nodeInfo_->clone();
4373  else
4374    nodeInfo_=NULL;
4375  objectiveValue_=rhs.objectiveValue_;
4376  guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4377  sumInfeasibilities_ = rhs.sumInfeasibilities_;
4378  if (rhs.branch_)
4379    branch_=rhs.branch_->clone();
4380  else
4381    branch_=NULL;
4382  depth_ = rhs.depth_;
4383  numberUnsatisfied_ = rhs.numberUnsatisfied_;
4384  nodeNumber_ = rhs.nodeNumber_;
4385  state_ = rhs.state_;
4386  if (nodeInfo_)
4387    assert ((state_&2)!=0);
4388  else
4389    assert ((state_&2)==0);
4390}
4391
4392CbcNode &
4393CbcNode::operator=(const CbcNode & rhs)
4394{
4395  if (this != &rhs) {
4396    delete nodeInfo_;
4397    if (rhs.nodeInfo_)
4398      nodeInfo_ = rhs.nodeInfo_->clone();
4399    else
4400      nodeInfo_ = NULL;
4401    objectiveValue_=rhs.objectiveValue_;
4402    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4403    sumInfeasibilities_ = rhs.sumInfeasibilities_;
4404    if (rhs.branch_)
4405      branch_=rhs.branch_->clone();
4406    else
4407      branch_=NULL,
4408        depth_ = rhs.depth_;
4409    numberUnsatisfied_ = rhs.numberUnsatisfied_;
4410    nodeNumber_ = rhs.nodeNumber_;
4411    state_ = rhs.state_;
4412    if (nodeInfo_)
4413      assert ((state_&2)!=0);
4414    else
4415      assert ((state_&2)==0);
4416  }
4417  return *this;
4418}
4419CbcNode::~CbcNode ()
4420{
4421#ifdef CHECK_NODE
4422  if (nodeInfo_) {
4423    printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
4424           this,nodeInfo_,nodeInfo_->numberPointingToThis());
4425    //assert(nodeInfo_->numberPointingToThis()>=0);
4426  } else {
4427    printf("CbcNode %x Destructor nodeInfo %x (?)\n",
4428           this,nodeInfo_);
4429  }
4430#endif
4431  if (nodeInfo_) {
4432    // was if (nodeInfo_&&(state_&2)!=0) {
4433    nodeInfo_->nullOwner();
4434    int numberToDelete=nodeInfo_->numberBranchesLeft();
4435    //    CbcNodeInfo * parent = nodeInfo_->parent();
4436    //assert (nodeInfo_->numberPointingToThis()>0);
4437    if (nodeInfo_->decrement(numberToDelete)==0||(state_&2)==0) {
4438      if ((state_&2)==0) 
4439        nodeInfo_->nullParent();
4440      delete nodeInfo_;
4441    } else {
4442      //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,nodeInfo_->parent());
4443      // anyway decrement parent
4444      //if (parent)
4445      ///parent->decrement(1);
4446    }
4447  }
4448  delete branch_;
4449}
4450// Decrement  active cut counts
4451void 
4452CbcNode::decrementCuts(int change)
4453{
4454  if (nodeInfo_)
4455    assert ((state_&2)!=0);
4456  else
4457    assert ((state_&2)==0);
4458  if(nodeInfo_) {
4459    nodeInfo_->decrementCuts(change);
4460  }
4461}
4462void 
4463CbcNode::decrementParentCuts(CbcModel * model, int change)
4464{
4465  if (nodeInfo_)
4466    assert ((state_&2)!=0);
4467  else
4468    assert ((state_&2)==0);
4469  if(nodeInfo_) {
4470    nodeInfo_->decrementParentCuts(model, change);
4471  }
4472}
4473
4474/*
4475  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
4476  in the attached nodeInfo_.
4477*/
4478void
4479CbcNode::initializeInfo()
4480{
4481  assert(nodeInfo_ && branch_) ;
4482  nodeInfo_->initializeInfo(branch_->numberBranches());
4483  assert ((state_&2)!=0);
4484  assert (nodeInfo_->numberBranchesLeft()==
4485          branch_->numberBranchesLeft());
4486}
4487// Nulls out node info
4488void 
4489CbcNode::nullNodeInfo()
4490{
4491  nodeInfo_=NULL;
4492  // say not active
4493  state_ &= ~2;
4494}
4495
4496int
4497CbcNode::branch(OsiSolverInterface * solver)
4498{
4499  double changeInGuessed;
4500  assert (nodeInfo_->numberBranchesLeft()==
4501          branch_->numberBranchesLeft());
4502  if (!solver)
4503    changeInGuessed=branch_->branch();
4504  else
4505    changeInGuessed=branch_->branch(solver);
4506  guessedObjectiveValue_+= changeInGuessed;
4507  //#define PRINTIT
4508#ifdef PRINTIT
4509  int numberLeft = nodeInfo_->numberBranchesLeft();
4510  CbcNodeInfo * parent = nodeInfo_->parent();
4511  int parentNodeNumber = -1;
4512  CbcBranchingObject * object1 = 
4513    dynamic_cast<CbcBranchingObject *>(branch_) ;
4514  //OsiObject * object = object1->
4515  //int sequence = object->columnNumber);
4516  int id=-1;
4517  double value=0.0;
4518  if (object1) {
4519    id = object1->variable();
4520    value = object1->value();
4521  }
4522  printf("id %d value %g objvalue %g\n",id,value,objectiveValue_);
4523  if (parent)
4524    parentNodeNumber = parent->nodeNumber();
4525  printf("Node number %d, %s, way %d, depth %d, parent node number %d\n",
4526         nodeInfo_->nodeNumber(),(numberLeft==2) ? "leftBranch" : "rightBranch",
4527         way(),depth_,parentNodeNumber);
4528  assert (parentNodeNumber!=nodeInfo_->nodeNumber());
4529#endif
4530  return nodeInfo_->branchedOn();
4531}
4532/* Active arm of the attached OsiBranchingObject.
4533   
4534   In the simplest instance, coded -1 for the down arm of the branch, +1 for
4535   the up arm. But see OsiBranchingObject::way()
4536   Use nodeInfo--.numberBranchesLeft_ to see how active
4537*/
4538int 
4539CbcNode::way() const
4540{
4541  if (branch_) {
4542    CbcBranchingObject * obj =
4543      dynamic_cast <CbcBranchingObject *>(branch_) ;
4544    if (obj) {
4545      return obj->way();
4546    } else {
4547      OsiTwoWayBranchingObject * obj2 =
4548        dynamic_cast <OsiTwoWayBranchingObject *>(branch_) ;
4549      assert (obj2);
4550      return obj2->way();
4551    }
4552  } else {
4553    return 0;
4554  }
4555}
4556/* Create a branching object for the node
4557   
4558   The routine scans the object list of the model and selects a set of
4559   unsatisfied objects as candidates for branching. The candidates are
4560   evaluated, and an appropriate branch object is installed.
4561   
4562   The numberPassesLeft is decremented to stop fixing one variable each time
4563   and going on and on (e.g. for stock cutting, air crew scheduling)
4564   
4565   If evaluation determines that an object is monotone or infeasible,
4566   the routine returns immediately. In the case of a monotone object,
4567   the branch object has already been called to modify the model.
4568   
4569   Return value:
4570   <ul>
4571   <li>  0: A branching object has been installed
4572   <li> -1: A monotone object was discovered
4573   <li> -2: An infeasible object was discovered
4574   </ul>
4575   Branch state:
4576   <ul>
4577   <li> -1: start
4578   <li> -1: A monotone object was discovered
4579   <li> -2: An infeasible object was discovered
4580   </ul>
4581*/
4582int 
4583CbcNode::chooseOsiBranch (CbcModel * model,
4584                          CbcNode * lastNode,
4585                          OsiBranchingInformation * usefulInfo,
4586                          int branchState)
4587{
4588  int returnStatus=0;
4589  if (lastNode)
4590    depth_ = lastNode->depth_+1;
4591  else
4592    depth_ = 0;
4593  OsiSolverInterface * solver = model->solver();
4594  objectiveValue_ = solver->getObjValue()*solver->getObjSense();
4595  usefulInfo->objectiveValue_ = objectiveValue_;
4596  usefulInfo->depth_ = depth_;
4597  const double * saveInfoSol = usefulInfo->solution_;
4598  double * saveSolution = new double[solver->getNumCols()];
4599  memcpy(saveSolution,solver->getColSolution(),solver->getNumCols()*sizeof(double));
4600  usefulInfo->solution_ = saveSolution;
4601  OsiChooseVariable * choose = model->branchingMethod()->chooseMethod();
4602  int numberUnsatisfied=-1;
4603  if (branchState<0) {
4604    // initialize
4605    // initialize sum of "infeasibilities"
4606    sumInfeasibilities_ = 0.0;
4607    numberUnsatisfied = choose->setupList(usefulInfo,true);
4608    numberUnsatisfied_ = numberUnsatisfied;
4609    branchState=0;
4610    if (numberUnsatisfied_<0) {
4611      // infeasible
4612      delete [] saveSolution;
4613      return -2;
4614    }
4615  }
4616  // unset best
4617  int best=-1;
4618  choose->setBestObjectIndex(-1);
4619  if (numberUnsatisfied) {
4620    if (branchState>0||!choose->numberOnList()) {
4621      // we need to return at once - don't do strong branching or anything
4622      if (choose->numberOnList()||!choose->numberStrong()) {
4623        best = choose->candidates()[0];
4624        choose->setBestObjectIndex(best);
4625      } else {
4626        // nothing on list - need to try again - keep any solution
4627        numberUnsatisfied = choose->setupList(usefulInfo, false);
4628        numberUnsatisfied_ = numberUnsatisfied;
4629        if (numberUnsatisfied) {
4630          best = choose->candidates()[0];
4631          choose->setBestObjectIndex(best);
4632        }
4633      }
4634    } else {
4635      // carry on with strong branching or whatever
4636      int returnCode = choose->chooseVariable(solver, usefulInfo,true);
4637      // update number of strong iterations etc
4638      model->incrementStrongInfo(choose->numberStrongDone(),choose->numberStrongIterations(),
4639                                 returnCode==-1 ? 0:choose->numberStrongFixed(),returnCode==-1);
4640      if (returnCode>1) {
4641        // has fixed some
4642        returnStatus=-1;
4643      } else if (returnCode==-1) {
4644        // infeasible
4645        returnStatus=-2;
4646      } else if (returnCode==0) {
4647        // normal
4648        returnStatus=0;
4649        numberUnsatisfied=1;
4650      } else {
4651        // ones on list satisfied - double check
4652        numberUnsatisfied = choose->setupList(usefulInfo, false);
4653        numberUnsatisfied_ = numberUnsatisfied;
4654        if (numberUnsatisfied) {
4655          best = choose->candidates()[0];
4656          choose->setBestObjectIndex(best);
4657        }
4658      }
4659    }
4660  } 
4661  delete branch_;
4662  branch_ = NULL;
4663  guessedObjectiveValue_ = COIN_DBL_MAX;//objectiveValue_; // for now
4664  if (!returnStatus) {
4665    if (numberUnsatisfied) {
4666      // create branching object
4667      const OsiObject * obj = model->solver()->object(choose->bestObjectIndex());
4668      //const OsiSolverInterface * solver = usefulInfo->solver_;
4669      branch_ = obj->createBranch(model->solver(),usefulInfo,obj->whichWay());
4670    }
4671  }
4672  usefulInfo->solution_=saveInfoSol;
4673  delete [] saveSolution;
4674  // may have got solution
4675  if (choose->goodSolution()
4676      &&model->problemFeasibility()->feasible(model,-1)>=0) {
4677    // yes
4678    double objValue = choose->goodObjectiveValue();
4679    model->setBestSolution(CBC_STRONGSOL,
4680                           objValue,
4681                           choose->goodSolution()) ;
4682    model->setLastHeuristic(NULL);
4683    model->incrementUsed(choose->goodSolution());
4684    choose->clearGoodSolution();
4685  }
4686  return returnStatus;
4687}
4688int 
4689CbcNode::chooseClpBranch (CbcModel * model,
4690                          CbcNode * lastNode)
4691{ 
4692  assert(lastNode);
4693  depth_ = lastNode->depth_+1;
4694  delete branch_;
4695  branch_=NULL;
4696  OsiSolverInterface * solver = model->solver();
4697  //double saveObjectiveValue = solver->getObjValue();
4698  //double objectiveValue = CoinMax(solver->getObjSense()*saveObjectiveValue,objectiveValue_);
4699  const double * lower = solver->getColLower();
4700  const double * upper = solver->getColUpper();
4701  // point to useful information
4702  OsiBranchingInformation usefulInfo = model->usefulInformation();
4703  // and modify
4704  usefulInfo.depth_=depth_;
4705  int i;
4706  //bool beforeSolution = model->getSolutionCount()==0;
4707  int numberObjects = model->numberObjects();
4708  int numberColumns = model->getNumCols();
4709  double * saveUpper = new double[numberColumns];
4710  double * saveLower = new double[numberColumns];
4711 
4712  // Save solution in case heuristics need good solution later
4713 
4714  double * saveSolution = new double[numberColumns];
4715  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
4716  model->reserveCurrentSolution(saveSolution);
4717  for (i=0;i<numberColumns;i++) {
4718    saveLower[i] = lower[i];
4719    saveUpper[i] = upper[i];
4720  }
4721  // Save basis
4722  CoinWarmStart * ws = solver->getWarmStart();
4723  numberUnsatisfied_ = 0;
4724  // initialize sum of "infeasibilities"
4725  sumInfeasibilities_ = 0.0;
4726  // Note looks as if off end (hidden one)
4727  OsiObject * object = model->modifiableObject(numberObjects);
4728  CbcGeneralDepth * thisOne = dynamic_cast <CbcGeneralDepth *> (object);
4729  assert (thisOne);
4730  OsiClpSolverInterface * clpSolver
4731    = dynamic_cast<OsiClpSolverInterface *> (solver);
4732  assert (clpSolver);
4733  ClpSimplex * simplex = clpSolver->getModelPtr();
4734  int preferredWay;
4735  double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
4736  if (thisOne->whichSolution()>=0) {
4737    ClpNode * nodeInfo = thisOne->nodeInfo(thisOne->whichSolution());
4738    nodeInfo->applyNode(simplex,2);
4739    int saveLogLevel = simplex->logLevel();
4740    simplex->setLogLevel(0);
4741    simplex->dual();
4742    simplex->setLogLevel(saveLogLevel);
4743    double cutoff = model->getCutoff();
4744    bool goodSolution=true;
4745    if (simplex->status()) {
4746      //simplex->writeMps("bad7.mps",2);
4747      if (nodeInfo->objectiveValue()>cutoff-1.0e-2)
4748        goodSolution=false;
4749      else
4750        assert (!simplex->status());
4751    }
4752    if (goodSolution) {
4753      double newObjectiveValue = solver->getObjSense()*solver->getObjValue();
4754      // See if integer solution
4755      int numInf;
4756      int numInf2;
4757      bool gotSol = model->feasibleSolution(numInf,numInf2);
4758      if (!gotSol) {
4759        printf("numinf %d\n",numInf);
4760        double * sol = simplex->primalColumnSolution();
4761        for (int i=0;i<numberColumns;i++) {
4762          if (simplex->isInteger(i)) {
4763            double value = floor(sol[i]+0.5);
4764            if (fabs(value-sol[i])>1.0e-7) {
4765              printf("%d value %g\n",i,sol[i]);
4766              if (fabs(value-sol[i])<1.0e-3) {
4767                sol[i]=value;
4768              }
4769            }
4770          }
4771        }
4772        simplex->writeMps("bad8.mps",2);
4773        bool gotSol = model->feasibleSolution(numInf,numInf2);
4774        if (!gotSol) 
4775          assert (gotSol);
4776      }
4777      model->setBestSolution(CBC_STRONGSOL,
4778                             newObjectiveValue,
4779                             solver->getColSolution()) ;
4780      model->setLastHeuristic(NULL);
4781      model->incrementUsed(solver->getColSolution());
4782    }
4783  }
4784  // restore bounds
4785  { for (int j=0;j<numberColumns;j++) {
4786      if (saveLower[j] != lower[j])
4787        solver->setColLower(j,saveLower[j]);
4788      if (saveUpper[j] != upper[j])
4789        solver->setColUpper(j,saveUpper[j]);
4790    }
4791  }
4792  // restore basis
4793  solver->setWarmStart(ws);
4794  delete ws;
4795  int anyAction;
4796  //#define CHECK_PATH
4797#ifdef CHECK_PATH
4798  extern int gotGoodNode_Z;
4799  if (gotGoodNode_Z>=0)
4800    printf("good node %d %g\n",gotGoodNode_Z,infeasibility);
4801#endif
4802  if (infeasibility>0.0) {
4803    if (infeasibility==COIN_DBL_MAX) {
4804      anyAction=-2; // infeasible
4805    } else {
4806      branch_=thisOne->createCbcBranch(solver,&usefulInfo,preferredWay);
4807      // Set to firat one (and change when re-pushing)
4808      CbcGeneralBranchingObject * branch = 
4809        dynamic_cast <CbcGeneralBranchingObject *> (branch_);
4810      branch->state(objectiveValue_,sumInfeasibilities_,
4811                    numberUnsatisfied_,0);
4812      branch->setNode(this);
4813      anyAction=0;
4814    }
4815  } else {
4816    anyAction=-1;
4817  }
4818#ifdef CHECK_PATH
4819  gotGoodNode_Z=-1;
4820#endif
4821  // Set guessed solution value
4822  guessedObjectiveValue_ = objectiveValue_+1.0e-5;
4823  delete [] saveLower;
4824  delete [] saveUpper;
4825 
4826  // restore solution
4827  solver->setColSolution(saveSolution);
4828  delete [] saveSolution;
4829  return anyAction;
4830}
4831/* Double checks in case node can change its mind!
4832   Returns objective value
4833   Can change objective etc */
4834double
4835CbcNode::checkIsCutoff(double cutoff)
4836{
4837  branch_->checkIsCutoff(cutoff);
4838  return objectiveValue_;
4839}
Note: See TracBrowser for help on using the repository browser.