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

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

fix odd accuracy problem when very large elements

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