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

Last change on this file since 838 was 838, checked in by forrest, 13 years ago

for deterministic parallel

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