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

Last change on this file since 912 was 912, checked in by ladanyi, 13 years ago

Incorporated changes from branches/heur

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