source: stable/2.2/Cbc/src/CbcNode.cpp @ 1060

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

out printf

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