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

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

compiler warnings, deterministic parallel and stability

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