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

Last change on this file since 862 was 854, checked in by forrest, 12 years ago

try and make a bit faster

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