source: trunk/CbcNode.cpp @ 271

Last change on this file since 271 was 271, checked in by forrest, 14 years ago

for Pierre (and to annoy Lou:-)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 129.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#include <string>
8//#define CBC_DEBUG 1
9//#define CHECK_CUT_COUNTS
10//#define CHECK_NODE
11#define CBC_WEAK_STRONG
12#include <cassert>
13#include <cfloat>
14#define CUTS
15#include "OsiSolverInterface.hpp"
16#include "OsiAuxInfo.hpp"
17#include "OsiSolverBranch.hpp"
18#include "CoinWarmStartBasis.hpp"
19#include "CoinTime.hpp"
20#include "CbcModel.hpp"
21#include "CbcNode.hpp"
22#include "CbcStatistics.hpp"
23#include "CbcStrategy.hpp"
24#include "CbcBranchActual.hpp"
25#include "CbcBranchDynamic.hpp"
26#include "OsiRowCut.hpp"
27#include "OsiRowCutDebugger.hpp"
28#include "OsiCuts.hpp"
29#include "CbcCountRowCut.hpp"
30#include "CbcFeasibilityBase.hpp"
31#include "CbcMessage.hpp"
32#include "OsiClpSolverInterface.hpp"
33#include "ClpSimplexOther.hpp"
34using namespace std;
35#include "CglCutGenerator.hpp"
36// Default Constructor
37CbcNodeInfo::CbcNodeInfo ()
38  :
39  numberPointingToThis_(0),
40  parent_(NULL),
41  owner_(NULL),
42  numberCuts_(0),
43  nodeNumber_(0),
44  cuts_(NULL),
45  numberRows_(0),
46  numberBranchesLeft_(0)
47{
48#ifdef CHECK_NODE
49  printf("CbcNodeInfo %x Constructor\n",this);
50#endif
51}
52// Constructor given parent
53CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent)
54  :
55  numberPointingToThis_(2),
56  parent_(parent),
57  owner_(NULL),
58  numberCuts_(0),
59  nodeNumber_(0),
60  cuts_(NULL),
61  numberRows_(0),
62  numberBranchesLeft_(2)
63{
64#ifdef CHECK_NODE
65  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
66#endif
67  if (parent_) {
68    numberRows_ = parent_->numberRows_+parent_->numberCuts_;
69  }
70}
71// Copy Constructor
72CbcNodeInfo::CbcNodeInfo (const CbcNodeInfo & rhs)
73  :
74  numberPointingToThis_(rhs.numberPointingToThis_),
75  parent_(rhs.parent_),
76  owner_(rhs.owner_),
77  numberCuts_(rhs.numberCuts_),
78  nodeNumber_(rhs.nodeNumber_),
79  cuts_(NULL),
80  numberRows_(rhs.numberRows_),
81  numberBranchesLeft_(rhs.numberBranchesLeft_)
82{
83#ifdef CHECK_NODE
84  printf("CbcNodeInfo %x Copy constructor\n",this);
85#endif
86  if (numberCuts_) {
87    cuts_ = new CbcCountRowCut * [numberCuts_];
88    int n=0;
89    for (int i=0;i<numberCuts_;i++) {
90      CbcCountRowCut * thisCut = rhs.cuts_[i];
91      if (thisCut) {
92        // I think this is correct - new one should take priority
93        thisCut->setInfo(this,n);
94        thisCut->increment(numberBranchesLeft_); 
95        cuts_[n++] = thisCut;
96      }
97    }
98    numberCuts_=n;
99  }
100}
101// Constructor given parent and owner
102CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent, CbcNode * owner)
103  :
104  numberPointingToThis_(2),
105  parent_(parent),
106  owner_(owner),
107  numberCuts_(0),
108  nodeNumber_(0),
109  cuts_(NULL),
110  numberRows_(0),
111  numberBranchesLeft_(2)
112{
113#ifdef CHECK_NODE
114  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
115#endif
116  if (parent_) {
117    numberRows_ = parent_->numberRows_+parent_->numberCuts_;
118  }
119}
120
121/**
122  Take care to detach from the owning CbcNode and decrement the reference
123  count in the parent.  If this is the last nodeInfo object pointing to the
124  parent, make a recursive call to delete the parent.
125*/
126CbcNodeInfo::~CbcNodeInfo()
127{
128#ifdef CHECK_NODE
129  printf("CbcNodeInfo %x Destructor parent %x\n",this,parent_);
130#endif
131
132  assert(!numberPointingToThis_);
133  // But they may be some left (max nodes?)
134  for (int i=0;i<numberCuts_;i++) 
135    delete cuts_[i];
136
137  delete [] cuts_;
138  if (owner_) 
139    owner_->nullNodeInfo();
140  if (parent_) {
141    int numberLinks = parent_->decrement();
142    if (!numberLinks) delete parent_;
143  }
144}
145
146
147//#define ALLCUTS
148void
149CbcNodeInfo::decrementCuts(int change)
150{
151  int i;
152  // get rid of all remaining if negative
153  int changeThis;
154  if (change<0)
155    changeThis = numberBranchesLeft_;
156  else
157    changeThis = change;
158 // decrement cut counts
159  for (i=0;i<numberCuts_;i++) {
160    if (cuts_[i]) {
161      int number = cuts_[i]->decrement(changeThis);
162      if (!number) {
163        //printf("info %x del cut %d %x\n",this,i,cuts_[i]);
164        delete cuts_[i];
165        cuts_[i]=NULL;
166      }
167    }
168  }
169}
170void
171CbcNodeInfo::decrementParentCuts(int change)
172{
173  if (parent_) {
174    // get rid of all remaining if negative
175    int changeThis;
176    if (change<0)
177      changeThis = numberBranchesLeft_;
178    else
179      changeThis = change;
180    int i;
181    // Get over-estimate of space needed for basis
182    CoinWarmStartBasis dummy;
183    dummy.setSize(0,numberRows_+numberCuts_);
184    buildRowBasis(dummy);
185    /* everything is zero (i.e. free) so we can use to see
186       if latest basis */
187    CbcNodeInfo * thisInfo = parent_;
188    while (thisInfo) 
189      thisInfo = thisInfo->buildRowBasis(dummy);
190    // decrement cut counts
191    thisInfo = parent_;
192    int numberRows=numberRows_;
193    while (thisInfo) {
194      for (i=thisInfo->numberCuts_-1;i>=0;i--) {
195        CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
196#ifdef ALLCUTS
197        status = CoinWarmStartBasis::isFree;
198#endif
199        if (thisInfo->cuts_[i]) {
200          int number=1;
201          if (status!=CoinWarmStartBasis::basic) {
202            // tight - drop 1 or 2
203            if (change<0)
204              number = thisInfo->cuts_[i]->decrement(changeThis);
205            else
206              number = thisInfo->cuts_[i]->decrement(change);
207          }
208          if (!number) {
209            delete thisInfo->cuts_[i];
210            thisInfo->cuts_[i]=NULL;
211          }
212        }
213      }
214      thisInfo = thisInfo->parent_;
215    }
216  }
217}
218
219void
220CbcNodeInfo::incrementParentCuts(int change)
221{
222  if (parent_) {
223    int i;
224    // Get over-estimate of space needed for basis
225    CoinWarmStartBasis dummy;
226    dummy.setSize(0,numberRows_+numberCuts_);
227    /* everything is zero (i.e. free) so we can use to see
228       if latest basis */
229    buildRowBasis(dummy);
230    CbcNodeInfo * thisInfo = parent_;
231    while (thisInfo) 
232      thisInfo = thisInfo->buildRowBasis(dummy);
233    // increment cut counts
234    thisInfo = parent_;
235    int numberRows=numberRows_;
236    while (thisInfo) {
237      for (i=thisInfo->numberCuts_-1;i>=0;i--) {
238        CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
239#ifdef ALLCUTS
240        status = CoinWarmStartBasis::isFree;
241#endif
242        if (thisInfo->cuts_[i]&&status!=CoinWarmStartBasis::basic) {
243          thisInfo->cuts_[i]->increment(change);
244        }
245      }
246      thisInfo = thisInfo->parent_;
247    }
248  }
249}
250
251/*
252  Append cuts to the cuts_ array in a nodeInfo. The initial reference count
253  is set to numberToBranchOn, which will normally be the number of arms
254  defined for the CbcBranchingObject attached to the CbcNode that owns this
255  CbcNodeInfo.
256*/
257void
258CbcNodeInfo::addCuts (OsiCuts & cuts, int numberToBranchOn,
259                      int * whichGenerator)
260{
261  int numberCuts = cuts.sizeRowCuts();
262  if (numberCuts) {
263    int i;
264    if (!numberCuts_) {
265      cuts_ = new CbcCountRowCut * [numberCuts];
266    } else {
267      CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
268      memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *));
269      delete [] cuts_;
270      cuts_ = temp;
271    }
272    for (i=0;i<numberCuts;i++) {
273      CbcCountRowCut * thisCut = new CbcCountRowCut(*cuts.rowCutPtr(i),
274                                                    this,numberCuts_);
275      thisCut->increment(numberToBranchOn); 
276      cuts_[numberCuts_++] = thisCut;
277#ifdef CBC_DEBUG
278      int n=thisCut->row().getNumElements();
279#if CBC_DEBUG>1
280      printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
281             thisCut->ub());
282#endif
283      int j;
284#if CBC_DEBUG>1
285      const int * index = thisCut->row().getIndices();
286#endif
287      const double * element = thisCut->row().getElements();
288      for (j=0;j<n;j++) {
289#if CBC_DEBUG>1
290        printf(" (%d,%g)",index[j],element[j]);
291#endif
292        assert(fabs(element[j])>1.00e-12);
293      }
294      printf("\n");
295#endif
296    }
297  }
298}
299
300void
301CbcNodeInfo::addCuts(int numberCuts, CbcCountRowCut ** cut, 
302                     int numberToBranchOn)
303{
304  if (numberCuts) {
305    int i;
306    if (!numberCuts_) {
307      cuts_ = new CbcCountRowCut * [numberCuts];
308    } else {
309      CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
310      memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *));
311      delete [] cuts_;
312      cuts_ = temp;
313    }
314    for (i=0;i<numberCuts;i++) {
315      CbcCountRowCut * thisCut = cut[i];
316      thisCut->setInfo(this,numberCuts_);
317      //printf("info %x cut %d %x\n",this,i,thisCut);
318      thisCut->increment(numberToBranchOn); 
319      cuts_[numberCuts_++] = thisCut;
320#ifdef CBC_DEBUG
321      int n=thisCut->row().getNumElements();
322#if CBC_DEBUG>1
323      printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
324             thisCut->ub());
325#endif
326      int j;
327#if CBC_DEBUG>1
328      const int * index = thisCut->row().getIndices();
329#endif
330      const double * element = thisCut->row().getElements();
331      for (j=0;j<n;j++) {
332#if CBC_DEBUG>1
333        printf(" (%d,%g)",index[j],element[j]);
334#endif
335        assert(fabs(element[j])>1.00e-12);
336      }
337      printf("\n");
338#endif
339    }
340  }
341}
342
343// delete cuts
344void
345CbcNodeInfo::deleteCuts(int numberToDelete, CbcCountRowCut ** cuts)
346{
347  int i;
348  int j;
349  int last=-1;
350  for (i=0;i<numberToDelete;i++) {
351    CbcCountRowCut * next = cuts[i];
352    for (j=last+1;j<numberCuts_;j++) {
353      if (next==cuts_[j])
354        break;
355    }
356    if (j==numberCuts_) {
357      // start from beginning
358      for (j=0;j<last;j++) {
359        if (next==cuts_[j])
360          break;
361      }
362      assert(j<last);
363    }
364    last=j;
365    int number = cuts_[j]->decrement();
366    if (!number)
367      delete cuts_[j];
368    cuts_[j]=NULL;
369  }
370  j=0;
371  for (i=0;i<numberCuts_;i++) {
372    if (cuts_[i])
373      cuts_[j++]=cuts_[i];
374  }
375  numberCuts_ = j;
376}
377
378// delete cuts
379void
380CbcNodeInfo::deleteCuts(int numberToDelete, int * which)
381{
382  int i;
383  for (i=0;i<numberToDelete;i++) {
384    int iCut=which[i];
385    int number = cuts_[iCut]->decrement();
386    if (!number)
387      delete cuts_[iCut];
388    cuts_[iCut]=NULL;
389  }
390  int j=0;
391  for (i=0;i<numberCuts_;i++) {
392    if (cuts_[i])
393      cuts_[j++]=cuts_[i];
394  }
395  numberCuts_ = j;
396}
397
398// Really delete a cut
399void 
400CbcNodeInfo::deleteCut(int whichOne)
401{
402  assert(whichOne<numberCuts_);
403  cuts_[whichOne]=NULL;
404}
405
406CbcFullNodeInfo::CbcFullNodeInfo() :
407  CbcNodeInfo(),
408  basis_(),
409  numberIntegers_(0),
410  lower_(NULL),
411  upper_(NULL)
412{
413}
414CbcFullNodeInfo::CbcFullNodeInfo(CbcModel * model,
415                                 int numberRowsAtContinuous) :
416  CbcNodeInfo()
417{
418  OsiSolverInterface * solver = model->solver();
419  numberRows_ = numberRowsAtContinuous;
420  numberIntegers_ = model->numberIntegers();
421  int numberColumns = model->getNumCols();
422  lower_ = new double [numberColumns];
423  upper_ = new double [numberColumns];
424  const double * lower = solver->getColLower();
425  const double * upper = solver->getColUpper();
426  int i;
427
428  for (i=0;i<numberColumns;i++) {
429    lower_[i]=lower[i];
430    upper_[i]=upper[i];
431  }
432
433  basis_ =  dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
434}
435
436CbcFullNodeInfo::CbcFullNodeInfo(const CbcFullNodeInfo & rhs) :
437  CbcNodeInfo(rhs)
438{ 
439  basis_= dynamic_cast<CoinWarmStartBasis *>(rhs.basis_->clone()) ;
440  numberIntegers_=rhs.numberIntegers_;
441  lower_=NULL;
442  upper_=NULL;
443  if (rhs.lower_!=NULL) {
444    int numberColumns = basis_->getNumStructural();
445    lower_ = new double [numberColumns];
446    upper_ = new double [numberColumns];
447    assert (upper_!=NULL);
448    memcpy(lower_,rhs.lower_,numberColumns*sizeof(double));
449    memcpy(upper_,rhs.upper_,numberColumns*sizeof(double));
450  }
451}
452
453CbcNodeInfo * 
454CbcFullNodeInfo::clone() const
455{ 
456  return (new CbcFullNodeInfo(*this));
457}
458
459CbcFullNodeInfo::~CbcFullNodeInfo ()
460{
461  delete basis_ ;
462  delete [] lower_;
463  delete [] upper_;
464}
465
466/*
467   The basis supplied as a parameter is deleted and replaced with a new basis
468   appropriate for the node, and lower and upper bounds on variables are
469   reset according to the stored bounds arrays. Any cuts associated with this
470   node are added to the list in addCuts, but not actually added to the
471   constraint system in the model.
472
473   Why pass in a basis at all? The short answer is ``We need the parameter to
474   pass out a basis, so might as well use it to pass in the size.''
475   
476   A longer answer is that in practice we take a memory allocation hit up in
477   addCuts1 (the only place applyToModel is called) when we setSize() the
478   basis that's passed in. It's immediately tossed here in favour of a clone
479   of the basis attached to this nodeInfo. This can probably be fixed, given
480   a bit of thought.
481*/
482
483void CbcFullNodeInfo::applyToModel (CbcModel *model,
484                                    CoinWarmStartBasis *&basis,
485                                    CbcCountRowCut **addCuts,
486                                    int &currentNumberCuts) const 
487
488{ OsiSolverInterface *solver = model->solver() ;
489
490  // branch - do bounds
491  int i;
492  solver->setColLower(lower_);
493  solver->setColUpper(upper_);
494  int numberColumns = model->getNumCols();
495  // move basis - but make sure size stays
496  // for bon-min - should not be needed int numberRows = model->getNumRows();
497  int numberRows=basis->getNumArtificial();
498  delete basis ;
499  if (basis_) {
500    basis = dynamic_cast<CoinWarmStartBasis *>(basis_->clone()) ;
501    basis->resize(numberRows,numberColumns);
502  } else {
503    // We have a solver without a basis
504    basis=NULL;
505  }
506  for (i=0;i<numberCuts_;i++) 
507    addCuts[currentNumberCuts+i]= cuts_[i];
508  currentNumberCuts += numberCuts_;
509  assert(!parent_);
510  return ;
511}
512
513/* Builds up row basis backwards (until original model).
514   Returns NULL or previous one to apply .
515   Depends on Free being 0 and impossible for cuts
516*/
517CbcNodeInfo * 
518CbcFullNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
519{
520  const unsigned int * saved = 
521    (const unsigned int *) basis_->getArtificialStatus();
522  unsigned int * now = 
523    (unsigned int *) basis.getArtificialStatus();
524  int number=basis_->getNumArtificial()>>4;;
525  int i;
526  for (i=0;i<number;i++) { 
527    if (!now[i])
528      now[i] = saved[i];
529  }
530  return NULL;
531}
532
533// Default constructor
534CbcPartialNodeInfo::CbcPartialNodeInfo()
535
536  : CbcNodeInfo(),
537    basisDiff_(NULL),
538    variables_(NULL),
539    newBounds_(NULL),
540    numberChangedBounds_(0)
541
542{ /* this space intentionally left blank */ }
543
544// Constructor from current state
545CbcPartialNodeInfo::CbcPartialNodeInfo (CbcNodeInfo *parent, CbcNode *owner,
546                                        int numberChangedBounds,
547                                        const int *variables,
548                                        const double *boundChanges,
549                                        const CoinWarmStartDiff *basisDiff)
550 : CbcNodeInfo(parent,owner)
551{
552  basisDiff_ = basisDiff->clone() ;
553
554  numberChangedBounds_ = numberChangedBounds;
555  variables_ = new int [numberChangedBounds_];
556  newBounds_ = new double [numberChangedBounds_];
557
558  int i ;
559  for (i=0;i<numberChangedBounds_;i++) {
560    variables_[i]=variables[i];
561    newBounds_[i]=boundChanges[i];
562  }
563}
564
565CbcPartialNodeInfo::CbcPartialNodeInfo (const CbcPartialNodeInfo & rhs)
566
567  : CbcNodeInfo(rhs.parent_)
568
569{ basisDiff_ = rhs.basisDiff_->clone() ;
570
571  numberChangedBounds_ = rhs.numberChangedBounds_;
572  variables_ = new int [numberChangedBounds_];
573  newBounds_ = new double [numberChangedBounds_];
574
575  int i ;
576  for (i=0;i<numberChangedBounds_;i++) {
577    variables_[i]=rhs.variables_[i];
578    newBounds_[i]=rhs.newBounds_[i];
579  }
580}
581
582CbcNodeInfo * 
583CbcPartialNodeInfo::clone() const
584{ 
585  return (new CbcPartialNodeInfo(*this));
586}
587
588
589CbcPartialNodeInfo::~CbcPartialNodeInfo ()
590{
591  delete basisDiff_ ;
592  delete [] variables_;
593  delete [] newBounds_;
594}
595
596
597/**
598   The basis supplied as a parameter is incrementally modified, and lower and
599   upper bounds on variables in the model are incrementally modified. Any
600   cuts associated with this node are added to the list in addCuts.
601*/
602
603void CbcPartialNodeInfo::applyToModel (CbcModel *model,
604                                       CoinWarmStartBasis *&basis,
605                                       CbcCountRowCut **addCuts,
606                                       int &currentNumberCuts) const 
607
608{ OsiSolverInterface *solver = model->solver();
609
610  basis->applyDiff(basisDiff_) ;
611
612  // branch - do bounds
613  int i;
614  for (i=0;i<numberChangedBounds_;i++) {
615    int variable = variables_[i];
616    if ((variable&0x80000000)==0) {
617      // lower bound changing
618      solver->setColLower(variable,newBounds_[i]);
619    } else {
620      // upper bound changing
621      solver->setColUpper(variable&0x7fffffff,newBounds_[i]);
622    }
623  }
624  for (i=0;i<numberCuts_;i++) {
625    addCuts[currentNumberCuts+i]= cuts_[i];
626    if (cuts_[i]&&model->messageHandler()->logLevel()>2) {
627      cuts_[i]->print();
628    }
629  }
630   
631  currentNumberCuts += numberCuts_;
632  return ;
633}
634
635/* Builds up row basis backwards (until original model).
636   Returns NULL or previous one to apply .
637   Depends on Free being 0 and impossible for cuts
638*/
639
640CbcNodeInfo * 
641CbcPartialNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
642
643{ basis.applyDiff(basisDiff_) ;
644
645  return parent_ ; }
646
647
648CbcNode::CbcNode() :
649  nodeInfo_(NULL),
650  objectiveValue_(1.0e100),
651  guessedObjectiveValue_(1.0e100),
652  branch_(NULL),
653  depth_(-1),
654  numberUnsatisfied_(0)
655{
656#ifdef CHECK_NODE
657  printf("CbcNode %x Constructor\n",this);
658#endif
659}
660
661CbcNode::CbcNode(CbcModel * model,
662                 CbcNode * lastNode) :
663  nodeInfo_(NULL),
664  objectiveValue_(1.0e100),
665  guessedObjectiveValue_(1.0e100),
666  branch_(NULL),
667  depth_(-1),
668  numberUnsatisfied_(0)
669{
670#ifdef CHECK_NODE
671  printf("CbcNode %x Constructor from model\n",this);
672#endif
673  model->setObjectiveValue(this,lastNode);
674
675  if (lastNode)
676    lastNode->nodeInfo_->increment();
677}
678
679
680void
681CbcNode::createInfo (CbcModel *model,
682                     CbcNode *lastNode,
683                     const CoinWarmStartBasis *lastws,
684                     const double *lastLower, const double *lastUpper,
685                     int numberOldActiveCuts,int numberNewCuts)
686{ OsiSolverInterface * solver = model->solver();
687 CbcStrategy * strategy = model->strategy();
688/*
689  The root --- no parent. Create full basis and bounds information.
690*/
691  if (!lastNode)
692  { 
693    if (!strategy)
694      nodeInfo_=new CbcFullNodeInfo(model,solver->getNumRows());
695    else
696      nodeInfo_ = strategy->fullNodeInfo(model,solver->getNumRows());
697  }
698/*
699  Not the root. Create an edit from the parent's basis & bound information.
700  This is not quite as straightforward as it seems. We need to reintroduce
701  cuts we may have dropped out of the basis, in the correct position, because
702  this whole process is strictly positional. Start by grabbing the current
703  basis.
704*/
705  else
706  { const CoinWarmStartBasis* ws =
707      dynamic_cast<const CoinWarmStartBasis*>(solver->getWarmStart());
708    assert(ws!=NULL); // make sure not volume
709    //int numberArtificials = lastws->getNumArtificial();
710    int numberColumns = solver->getNumCols();
711   
712    const double * lower = solver->getColLower();
713    const double * upper = solver->getColUpper();
714
715    int i;
716/*
717  Create a clone and resize it to hold all the structural constraints, plus
718  all the cuts: old cuts, both active and inactive (currentNumberCuts), and
719  new cuts (numberNewCuts).
720
721  TODO: You'd think that the set of constraints (logicals) in the expanded
722        basis should match the set represented in lastws. At least, that's
723        what I thought. But at the point I first looked hard at this bit of
724        code, it turned out that lastws was the stripped basis produced at
725        the end of addCuts(), rather than the raw basis handed back by
726        addCuts1(). The expanded basis here is equivalent to the raw basis of
727        addCuts1(). I said ``whoa, that's not good, I must have introduced a
728        bug'' and went back to John's code to see where I'd gone wrong.
729        And discovered the same `error' in his code.
730
731        After a bit of thought, my conclusion is that correctness is not
732        affected by whether lastws is the stripped or raw basis. The diffs
733        have no semantics --- just a set of changes that need to be made
734        to convert lastws into expanded. I think the only effect is that we
735        store a lot more diffs (everything in expanded that's not covered by
736        the stripped basis). But I need to give this more thought. There
737        may well be some subtle error cases.
738
739        In the mean time, I've twiddled addCuts() to set lastws to the raw
740        basis. Makes me (Lou) less nervous to compare apples to apples.
741*/
742    CoinWarmStartBasis *expanded = 
743      dynamic_cast<CoinWarmStartBasis *>(ws->clone()) ;
744    int numberRowsAtContinuous = model->numberRowsAtContinuous();
745    int iFull = numberRowsAtContinuous+model->currentNumberCuts()+
746      numberNewCuts;
747    //int numberArtificialsNow = iFull;
748    //int maxBasisLength = ((iFull+15)>>4)+((numberColumns+15)>>4);
749    //printf("l %d full %d\n",maxBasisLength,iFull);
750    if (expanded) 
751      expanded->resize(iFull,numberColumns);
752#ifdef FULL_DEBUG
753    printf("Before expansion: orig %d, old %d, new %d, current %d\n",
754           numberRowsAtContinuous,numberOldActiveCuts,numberNewCuts,
755           model->currentNumberCuts()) ;
756    ws->print();
757#endif
758/*
759  Now fill in the expanded basis. Any indices beyond nPartial must
760  be cuts created while processing this node --- they can be copied directly
761  into the expanded basis. From nPartial down, pull the status of active cuts
762  from ws, interleaving with a B entry for the deactivated (loose) cuts.
763*/
764    int numberDropped = model->currentNumberCuts()-numberOldActiveCuts;
765    int iCompact=iFull-numberDropped;
766    CbcCountRowCut ** cut = model->addedCuts();
767    int nPartial = model->currentNumberCuts()+numberRowsAtContinuous;
768    iFull--;
769    for (;iFull>=nPartial;iFull--) {
770      CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
771      //assert (status != CoinWarmStartBasis::basic); // may be permanent cut
772      expanded->setArtifStatus(iFull,status);
773    }
774    for (;iFull>=numberRowsAtContinuous;iFull--) {
775      if (cut[iFull-numberRowsAtContinuous]) {
776        CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
777        // If no cut generator being used then we may have basic variables
778        //if (model->getMaximumCutPasses()&&
779        //  status == CoinWarmStartBasis::basic)
780        //printf("cut basic\n");
781        expanded->setArtifStatus(iFull,status);
782      } else {
783        expanded->setArtifStatus(iFull,CoinWarmStartBasis::basic);
784      }
785    }
786#ifdef FULL_DEBUG
787    printf("Expanded basis\n");
788    expanded->print() ;
789    printf("Diffing against\n") ;
790    lastws->print() ;
791#endif   
792/*
793  Now that we have two bases in proper positional correspondence, creating
794  the actual diff is dead easy.
795*/
796
797    CoinWarmStartDiff *basisDiff = expanded->generateDiff(lastws) ;
798/*
799  Diff the bound vectors. It's assumed the number of structural variables is
800  not changing. Assuming that branching objects all involve integer variables,
801  we should see at least one bound change as a consequence of processing this
802  subproblem. Different types of branching objects could break this assertion.
803  Not true at all - we have not applied current branch - JJF.
804*/
805    double *boundChanges = new double [2*numberColumns] ;
806    int *variables = new int [2*numberColumns] ;
807    int numberChangedBounds=0;
808    for (i=0;i<numberColumns;i++) {
809      if (lower[i]!=lastLower[i]) {
810        variables[numberChangedBounds]=i;
811        boundChanges[numberChangedBounds++]=lower[i];
812      }
813      if (upper[i]!=lastUpper[i]) {
814        variables[numberChangedBounds]=i|0x80000000;
815        boundChanges[numberChangedBounds++]=upper[i];
816      }
817#ifdef CBC_DEBUG
818      if (lower[i]!=lastLower[i])
819        printf("lower on %d changed from %g to %g\n",
820               i,lastLower[i],lower[i]);
821      if (upper[i]!=lastUpper[i])
822        printf("upper on %d changed from %g to %g\n",
823               i,lastUpper[i],upper[i]);
824#endif
825    }
826#ifdef CBC_DEBUG
827    printf("%d changed bounds\n",numberChangedBounds) ;
828#endif
829    //if (lastNode->branchingObject()->boundBranch())
830    //assert (numberChangedBounds);
831/*
832  Hand the lot over to the CbcPartialNodeInfo constructor, then clean up and
833  return.
834*/
835    if (!strategy)
836      nodeInfo_ =
837        new CbcPartialNodeInfo(lastNode->nodeInfo_,this,numberChangedBounds,
838                               variables,boundChanges,basisDiff) ;
839    else
840      nodeInfo_ = strategy->partialNodeInfo(model, lastNode->nodeInfo_,this,numberChangedBounds,
841                               variables,boundChanges,basisDiff) ;
842    delete basisDiff ;
843    delete [] boundChanges;
844    delete [] variables;
845    delete expanded ;
846    delete ws;
847  }
848  // Set node number
849  nodeInfo_->setNodeNumber(model->getNodeCount2());
850}
851
852/*
853  The routine scans through the object list of the model looking for objects
854  that indicate infeasibility. It tests each object using strong branching
855  and selects the one with the least objective degradation.  A corresponding
856  branching object is left attached to lastNode.
857
858  If strong branching is disabled, a candidate object is chosen essentially
859  at random (whatever object ends up in pos'n 0 of the candidate array).
860
861  If a branching candidate is found to be monotone, bounds are set to fix the
862  variable and the routine immediately returns (the caller is expected to
863  reoptimize).
864
865  If a branching candidate is found to result in infeasibility in both
866  directions, the routine immediately returns an indication of infeasibility.
867
868  Returns:  0   both branch directions are feasible
869           -1   branching variable is monotone
870           -2   infeasible
871
872  Original comments:
873    Here could go cuts etc etc
874    For now just fix on objective from strong branching.
875*/
876
877int CbcNode::chooseBranch (CbcModel *model, CbcNode *lastNode,int numberPassesLeft)
878
879{ if (lastNode)
880    depth_ = lastNode->depth_+1;
881  else
882    depth_ = 0;
883  delete branch_;
884  branch_=NULL;
885  OsiSolverInterface * solver = model->solver();
886  double saveObjectiveValue = solver->getObjValue();
887  double objectiveValue = CoinMax(solver->getObjSense()*saveObjectiveValue,objectiveValue_);
888  const double * lower = solver->getColLower();
889  const double * upper = solver->getColUpper();
890  // See what user thinks
891  int anyAction=model->problemFeasibility()->feasible(model,0);
892  if (anyAction) {
893    // will return -2 if infeasible , 0 if treat as integer
894    return anyAction-1;
895  }
896  double integerTolerance = 
897    model->getDblParam(CbcModel::CbcIntegerTolerance);
898  int i;
899  bool beforeSolution = model->getSolutionCount()==0;
900  int numberStrong=model->numberStrong();
901  // switch off strong if hotstart
902  if (model->hotstartSolution())
903    numberStrong=0;
904  int numberStrongDone=0;
905  int numberUnfinished=0;
906  int numberStrongInfeasible=0;
907  int numberStrongIterations=0;
908  int saveNumberStrong=numberStrong;
909  int numberObjects = model->numberObjects();
910  bool checkFeasibility = numberObjects>model->numberIntegers();
911  int maximumStrong = CoinMax(CoinMin(model->numberStrong(),numberObjects),1);
912  int numberColumns = model->getNumCols();
913  double * saveUpper = new double[numberColumns];
914  double * saveLower = new double[numberColumns];
915
916  // Save solution in case heuristics need good solution later
917 
918  double * saveSolution = new double[numberColumns];
919  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
920  model->reserveCurrentSolution(saveSolution);
921  /*
922    Get a branching decision object. Use the default decision criteria unless
923    the user has loaded a decision method into the model.
924  */
925  CbcBranchDecision *decision = model->branchingMethod();
926  if (!decision)
927    decision = new CbcBranchDefaultDecision();
928
929  CbcStrongInfo * choice = new CbcStrongInfo[maximumStrong];
930  for (i=0;i<numberColumns;i++) {
931    saveLower[i] = lower[i];
932    saveUpper[i] = upper[i];
933  }
934  // May go round twice if strong branching fixes all local candidates
935  bool finished=false;
936  double estimatedDegradation=0.0; 
937  while(!finished) {
938    finished=true;
939    // Some objects may compute an estimate of best solution from here
940    estimatedDegradation=0.0; 
941    //int numberIntegerInfeasibilities=0; // without odd ones
942    numberStrongDone=0;
943    numberUnfinished=0;
944    numberStrongInfeasible=0;
945    numberStrongIterations=0;
946   
947    // We may go round this loop twice (only if we think we have solution)
948    for (int iPass=0;iPass<2;iPass++) {
949     
950      // compute current state
951      //int numberObjectInfeasibilities; // just odd ones
952      //model->feasibleSolution(
953      //                      numberIntegerInfeasibilities,
954      //                      numberObjectInfeasibilities);
955      const double * hotstartSolution = model->hotstartSolution();
956      const int * hotstartPriorities = model->hotstartPriorities();
957     
958      // Some objects may compute an estimate of best solution from here
959      estimatedDegradation=0.0; 
960      numberUnsatisfied_ = 0;
961      int bestPriority=INT_MAX;
962      /*
963        Scan for branching objects that indicate infeasibility. Choose the best
964        maximumStrong candidates, using priority as the first criteria, then
965        integer infeasibility.
966       
967        The algorithm is to fill the choice array with a set of good candidates (by
968        infeasibility) with priority bestPriority.  Finding a candidate with
969        priority better (less) than bestPriority flushes the choice array. (This
970        serves as initialization when the first candidate is found.)
971       
972        A new candidate is added to choices only if its infeasibility exceeds the
973        current max infeasibility (mostAway). When a candidate is added, it
974        replaces the candidate with the smallest infeasibility (tracked by
975        iSmallest).
976      */
977      int iSmallest = 0;
978      double mostAway = 1.0e-100;
979      for (i = 0 ; i < maximumStrong ; i++)
980        choice[i].possibleBranch = NULL ;
981      numberStrong=0;
982      for (i=0;i<numberObjects;i++) {
983        CbcObject * object = model->modifiableObject(i);
984        int preferredWay;
985        double infeasibility = object->infeasibility(preferredWay);
986        int priorityLevel = object->priority();
987        if (hotstartSolution) {
988          // we are doing hot start
989          const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
990          if (thisOne) {
991            int iColumn = thisOne->modelSequence();
992            double targetValue = hotstartSolution[iColumn];
993            if (saveUpper[iColumn]>saveLower[iColumn]) {
994              double value = saveSolution[iColumn];
995              if (hotstartPriorities)
996                priorityLevel=hotstartPriorities[iColumn]; 
997              //double originalLower = thisOne->originalLower();
998              //double originalUpper = thisOne->originalUpper();
999              // switch off if not possible
1000              if (targetValue>=saveLower[iColumn]&&targetValue<=saveUpper[iColumn]) {
1001                /* priority outranks rest always if negative
1002                   otherwise can be downgraded if at correct level.
1003                   Infeasibility may be increased to choose 1.0 values first.
1004                   choose one near wanted value
1005                */
1006                if (fabs(value-targetValue)>integerTolerance) {
1007                  infeasibility = 1.0-fabs(value-targetValue);
1008                  if (targetValue==1.0)
1009                    infeasibility += 1.0;
1010                  if (value>targetValue) {
1011                    preferredWay=-1;
1012                  } else {
1013                    preferredWay=1;
1014                  }
1015                  priorityLevel = CoinAbs(priorityLevel);
1016                } else if (priorityLevel<0) {
1017                  priorityLevel = CoinAbs(priorityLevel);
1018                  if (targetValue==saveLower[iColumn]) {
1019                    infeasibility = integerTolerance+1.0e-12;
1020                    preferredWay=-1;
1021                  } else if (targetValue==saveUpper[iColumn]) {
1022                    infeasibility = integerTolerance+1.0e-12;
1023                    preferredWay=1;
1024                  } else {
1025                    // can't
1026                    priorityLevel += 10000000;
1027                  }
1028                } else {
1029                  priorityLevel += 10000000;
1030                }
1031              } else {
1032                // switch off if not possible
1033                hotstartSolution=NULL;
1034                model->setHotstartSolution(NULL,NULL);
1035              }
1036            } else if (targetValue<saveLower[iColumn]||targetValue>saveUpper[iColumn]) {
1037              // switch off as not possible
1038              hotstartSolution=NULL;
1039              model->setHotstartSolution(NULL,NULL);
1040            }
1041          } else {
1042            priorityLevel += 10000000;
1043          }
1044        }
1045        if (infeasibility) {
1046          // Increase estimated degradation to solution
1047          estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate());
1048          numberUnsatisfied_++;
1049          // Better priority? Flush choices.
1050          if (priorityLevel<bestPriority) {
1051            int j;
1052            iSmallest=0;
1053            for (j=0;j<maximumStrong;j++) {
1054              choice[j].upMovement=0.0;
1055              delete choice[j].possibleBranch;
1056              choice[j].possibleBranch=NULL;
1057            }
1058            bestPriority = priorityLevel;
1059            mostAway=1.0e-100;
1060            numberStrong=0;
1061          } else if (priorityLevel>bestPriority) {
1062            continue;
1063          }
1064          // Check for suitability based on infeasibility.
1065          if (infeasibility>mostAway) {
1066            //add to list
1067            choice[iSmallest].upMovement=infeasibility;
1068            delete choice[iSmallest].possibleBranch;
1069            choice[iSmallest].possibleBranch=object->createBranch(preferredWay);
1070            numberStrong = CoinMax(numberStrong,iSmallest+1);
1071            // Save which object it was
1072            choice[iSmallest].objectNumber=i;
1073            int j;
1074            iSmallest=-1;
1075            mostAway = 1.0e50;
1076            for (j=0;j<maximumStrong;j++) {
1077              if (choice[j].upMovement<mostAway) {
1078                mostAway=choice[j].upMovement;
1079                iSmallest=j;
1080              }
1081            }
1082          }
1083        }
1084      }
1085      if (numberUnsatisfied_) {
1086        // some infeasibilities - go to next steps
1087        break;
1088      } else if (!iPass) {
1089        // looks like a solution - get paranoid
1090        bool roundAgain=false;
1091        // get basis
1092        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
1093        if (!ws)
1094          break;
1095        for (i=0;i<numberColumns;i++) {
1096          double value = saveSolution[i];
1097          if (value<lower[i]) {
1098            saveSolution[i]=lower[i];
1099            roundAgain=true;
1100            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
1101          } else if (value>upper[i]) {
1102            saveSolution[i]=upper[i];
1103            roundAgain=true;
1104            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
1105          } 
1106        }
1107        if (roundAgain&&saveNumberStrong) {
1108          // restore basis
1109          solver->setWarmStart(ws);
1110          delete ws;
1111          solver->resolve();
1112          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1113          model->reserveCurrentSolution(saveSolution);
1114          if (!solver->isProvenOptimal()) {
1115            // infeasible
1116            anyAction=-2;
1117            break;
1118          }
1119        } else {
1120          delete ws;
1121          break;
1122        }
1123      }
1124    }
1125    /* Some solvers can do the strong branching calculations faster if
1126       they do them all at once.  At present only Clp does for ordinary
1127       integers but I think this coding would be easy to modify
1128    */
1129    bool allNormal=true; // to say if we can do fast strong branching
1130    // Say which one will be best
1131    int bestChoice=0;
1132    double worstInfeasibility=0.0;
1133    for (i=0;i<numberStrong;i++) {
1134      choice[i].numIntInfeasUp = numberUnsatisfied_;
1135      choice[i].numIntInfeasDown = numberUnsatisfied_;
1136      choice[i].fix=0; // say not fixed
1137      if (!dynamic_cast <const CbcSimpleInteger *> (model->object(choice[i].objectNumber)))
1138        allNormal=false; // Something odd so lets skip clever fast branching
1139      if ( !model->object(choice[i].objectNumber)->boundBranch())
1140        numberStrong=0; // switch off
1141      if ( choice[i].possibleBranch->numberBranches()>2)
1142        numberStrong=0; // switch off
1143      // Do best choice in case switched off
1144      if (choice[i].upMovement>worstInfeasibility) {
1145        worstInfeasibility=choice[i].upMovement;
1146        bestChoice=i;
1147      }
1148    }
1149    // If we have hit max time don't do strong branching
1150    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
1151                        model->getDblParam(CbcModel::CbcMaximumSeconds));
1152    // also give up if we are looping round too much
1153    if (hitMaxTime||numberPassesLeft<=0)
1154      numberStrong=0;
1155    /*
1156      Is strong branching enabled? If so, set up and do it. Otherwise, we'll
1157      fall through to simple branching.
1158     
1159      Setup for strong branching involves saving the current basis (for restoration
1160      afterwards) and setting up for hot starts.
1161    */
1162    if (numberStrong&&saveNumberStrong) {
1163     
1164      bool solveAll=false; // set true to say look at all even if some fixed (experiment)
1165      solveAll=true;
1166      // worth trying if too many times
1167      // Save basis
1168      CoinWarmStart * ws = solver->getWarmStart();
1169      // save limit
1170      int saveLimit;
1171      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
1172      if (beforeSolution&&saveLimit<100)
1173        solver->setIntParam(OsiMaxNumIterationHotStart,100); // go to end
1174     
1175      /* If we are doing all strong branching in one go then we create new arrays
1176         to store information.  If clp NULL then doing old way.
1177         Going down -
1178         outputSolution[2*i] is final solution.
1179         outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown
1180         outputStuff[2*i+numberStrong] is number iterations
1181         On entry newUpper[i] is new upper bound, on exit obj change
1182         Going up -
1183         outputSolution[2*i+1] is final solution.
1184         outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
1185         outputStuff[2*i+1+numberStrong] is number iterations
1186       On entry newLower[i] is new lower bound, on exit obj change
1187      */
1188      OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
1189      ClpSimplex * clp=NULL;
1190      double * newLower = NULL;
1191      double * newUpper = NULL;
1192      double ** outputSolution=NULL;
1193      int * outputStuff=NULL;
1194      // Go back to normal way if user wants it
1195      if (osiclp&&(osiclp->specialOptions()&16)!=0&&osiclp->specialOptions()>0)
1196        allNormal=false;
1197      if (osiclp&&!allNormal) {
1198        // say do fast
1199        int easy=1;
1200        osiclp->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
1201      }
1202      if (osiclp&& allNormal) {
1203        clp = osiclp->getModelPtr();
1204        // Clp - do a different way
1205        newLower = new double[numberStrong];
1206        newUpper = new double[numberStrong];
1207        outputSolution = new double * [2*numberStrong];
1208        outputStuff = new int [4*numberStrong];
1209        int * which = new int[numberStrong];
1210        int startFinishOptions;
1211        int specialOptions = osiclp->specialOptions();
1212        int clpOptions = clp->specialOptions();
1213        int returnCode=0;
1214#define CRUNCH
1215#ifdef CRUNCH
1216        // Crunch down problem
1217        int numberRows = clp->numberRows();
1218        // Use dual region
1219        double * rhs = clp->dualRowSolution();
1220        int * whichRow = new int[3*numberRows];
1221        int * whichColumn = new int[2*numberColumns];
1222        int nBound;
1223        ClpSimplex * small = ((ClpSimplexOther *) clp)->crunch(rhs,whichRow,whichColumn,nBound,true);
1224        if (!small) {
1225          anyAction=-2;
1226          //printf("XXXX Inf by inspection\n");
1227          delete [] whichColumn;
1228          whichColumn=NULL;
1229          delete [] whichRow;
1230          whichRow=NULL;
1231          break;
1232        } else {
1233          clp = small;
1234        }
1235#else
1236        int saveLogLevel = clp->logLevel();
1237        int saveMaxIts = clp->maximumIterations();
1238#endif
1239        clp->setLogLevel(0);
1240        if((specialOptions&1)==0) {
1241          startFinishOptions=0;
1242          clp->setSpecialOptions(clpOptions|(64|1024));
1243        } else {
1244          startFinishOptions=1+2+4;
1245          //startFinishOptions=1+4; // for moment re-factorize
1246          if((specialOptions&4)==0) 
1247            clp->setSpecialOptions(clpOptions|(64|128|512|1024|4096));
1248          else
1249            clp->setSpecialOptions(clpOptions|(64|128|512|1024|2048|4096));
1250        }
1251        // User may want to clean up before strong branching
1252        if ((clp->specialOptions()&32)!=0) {
1253          clp->primal(1);
1254          if (clp->numberIterations())
1255            model->messageHandler()->message(CBC_ITERATE_STRONG,model->messages())
1256              << clp->numberIterations()
1257              <<CoinMessageEol;
1258        }
1259        clp->setMaximumIterations(saveLimit);
1260#ifdef CRUNCH
1261        int * backColumn = whichColumn+numberColumns;
1262#endif
1263        for (i=0;i<numberStrong;i++) {
1264          int iObject = choice[i].objectNumber;
1265          const CbcObject * object = model->object(iObject);
1266          const CbcSimpleInteger * simple = dynamic_cast <const CbcSimpleInteger *> (object);
1267          int iSequence = simple->modelSequence();
1268          newLower[i]= ceil(saveSolution[iSequence]);
1269          newUpper[i]= floor(saveSolution[iSequence]);
1270#ifdef CRUNCH
1271          iSequence = backColumn[iSequence];
1272          assert (iSequence>=0);
1273#endif
1274          which[i]=iSequence;
1275          outputSolution[2*i]= new double [numberColumns];
1276          outputSolution[2*i+1]= new double [numberColumns];
1277        }
1278        //clp->writeMps("bad");
1279        returnCode=clp->strongBranching(numberStrong,which,
1280                                            newLower, newUpper,outputSolution,
1281                                            outputStuff,outputStuff+2*numberStrong,!solveAll,false,
1282                                            startFinishOptions);
1283#ifndef CRUNCH
1284        clp->setSpecialOptions(clpOptions); // restore
1285        clp->setMaximumIterations(saveMaxIts);
1286        clp->setLogLevel(saveLogLevel);
1287#endif
1288        if (returnCode==-2) {
1289          // bad factorization!!!
1290          // Doing normal way
1291          // Mark hot start
1292          solver->markHotStart();
1293          clp = NULL;
1294        } else {
1295#ifdef CRUNCH
1296          // extract solution
1297          //bool checkSol=true;
1298          for (i=0;i<numberStrong;i++) {
1299            int iObject = choice[i].objectNumber;
1300            const CbcObject * object = model->object(iObject);
1301            const CbcSimpleInteger * simple = dynamic_cast <const CbcSimpleInteger *> (object);
1302            int iSequence = simple->modelSequence();
1303            which[i]=iSequence;
1304            double * sol = outputSolution[2*i];
1305            double * sol2 = outputSolution[2*i+1];
1306            //bool x=true;
1307            //bool x2=true;
1308            for (int iColumn=numberColumns-1;iColumn>=0;iColumn--) {
1309              int jColumn = backColumn[iColumn];
1310              if (jColumn>=0) {
1311                sol[iColumn]=sol[jColumn];
1312                sol2[iColumn]=sol2[jColumn];
1313              } else {
1314                sol[iColumn]=saveSolution[iColumn];
1315                sol2[iColumn]=saveSolution[iColumn];
1316              }
1317            }
1318          }
1319#endif
1320        }
1321#ifdef CRUNCH
1322        delete [] whichColumn;
1323        delete [] whichRow;
1324        delete small;
1325#endif
1326        delete [] which;
1327      } else {
1328        // Doing normal way
1329        // Mark hot start
1330        solver->markHotStart();
1331      }
1332      /*
1333        Open a loop to do the strong branching LPs. For each candidate variable,
1334        solve an LP with the variable forced down, then up. If a direction turns
1335        out to be infeasible or monotonic (i.e., over the dual objective cutoff),
1336        force the objective change to be big (1.0e100). If we determine the problem
1337        is infeasible, or find a monotone variable, escape the loop.
1338       
1339        TODO: The `restore bounds' part might be better encapsulated as an
1340        unbranch() method. Branching objects more exotic than simple integers
1341        or cliques might not restrict themselves to variable bounds.
1342
1343        TODO: Virtuous solvers invalidate the current solution (or give bogus
1344        results :-) when the bounds are changed out from under them. So we
1345        need to do all the work associated with finding a new solution before
1346        restoring the bounds.
1347      */
1348      for (i = 0 ; i < numberStrong ; i++)
1349        { double objectiveChange ;
1350        double newObjectiveValue=1.0e100;
1351        // status is 0 finished, 1 infeasible and other
1352        int iStatus;
1353        /*
1354          Try the down direction first. (Specify the initial branching alternative as
1355          down with a call to way(-1). Each subsequent call to branch() performs the
1356          specified branch and advances the branch object state to the next branch
1357          alternative.)
1358        */
1359        if (!clp) {
1360          choice[i].possibleBranch->way(-1) ;
1361          choice[i].possibleBranch->branch() ;
1362          bool feasible=true;
1363          if (checkFeasibility) {
1364            // check branching did not make infeasible
1365            int iColumn;
1366            int numberColumns = solver->getNumCols();
1367            const double * columnLower = solver->getColLower();
1368            const double * columnUpper = solver->getColUpper();
1369            for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1370              if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1371                feasible=false;
1372            }
1373          }
1374          if (feasible) {
1375            solver->solveFromHotStart() ;
1376            numberStrongDone++;
1377            numberStrongIterations += solver->getIterationCount();
1378            /*
1379              We now have an estimate of objective degradation that we can use for strong
1380              branching. If we're over the cutoff, the variable is monotone up.
1381              If we actually made it to optimality, check for a solution, and if we have
1382              a good one, call setBestSolution to process it. Note that this may reduce the
1383              cutoff, so we check again to see if we can declare this variable monotone.
1384            */
1385            if (solver->isProvenOptimal())
1386              iStatus=0; // optimal
1387            else if (solver->isIterationLimitReached()
1388                     &&!solver->isDualObjectiveLimitReached())
1389              iStatus=2; // unknown
1390            else
1391              iStatus=1; // infeasible
1392            newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1393            choice[i].numItersDown = solver->getIterationCount();
1394          } else {
1395            iStatus=1; // infeasible
1396            newObjectiveValue = 1.0e100;
1397            choice[i].numItersDown = 0;
1398          }
1399        } else {
1400          iStatus = outputStuff[2*i];
1401          choice[i].numItersDown = outputStuff[2*numberStrong+2*i];
1402          numberStrongDone++;
1403          numberStrongIterations += choice[i].numItersDown;
1404          newObjectiveValue = objectiveValue+newUpper[i];
1405          solver->setColSolution(outputSolution[2*i]);
1406        }
1407        objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
1408        if (!iStatus) {
1409          choice[i].finishedDown = true ;
1410          if (newObjectiveValue>=model->getCutoff()) {
1411            objectiveChange = 1.0e100; // say infeasible
1412            numberStrongInfeasible++;
1413          } else {
1414            // See if integer solution
1415            if (model->feasibleSolution(choice[i].numIntInfeasDown,
1416                                        choice[i].numObjInfeasDown)
1417                &&model->problemFeasibility()->feasible(model,-1)>=0) {
1418              model->setBestSolution(CBC_STRONGSOL,
1419                                     newObjectiveValue,
1420                                     solver->getColSolution()) ;
1421              // only needed for odd solvers
1422              newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1423              objectiveChange = CoinMax(newObjectiveValue-objectiveValue_,0.0) ;
1424              model->setLastHeuristic(NULL);
1425              model->incrementUsed(solver->getColSolution());
1426              if (newObjectiveValue >= model->getCutoff()) {    //  *new* cutoff
1427                objectiveChange = 1.0e100 ;
1428                numberStrongInfeasible++;
1429              }
1430            }
1431          }
1432        } else if (iStatus==1) {
1433          objectiveChange = 1.0e100 ;
1434          numberStrongInfeasible++;
1435        } else {
1436          // Can't say much as we did not finish
1437          choice[i].finishedDown = false ;
1438          numberUnfinished++;
1439        }
1440        choice[i].downMovement = objectiveChange ;
1441       
1442        // restore bounds
1443        if (!clp)
1444          { for (int j=0;j<numberColumns;j++) {
1445            if (saveLower[j] != lower[j])
1446              solver->setColLower(j,saveLower[j]);
1447            if (saveUpper[j] != upper[j])
1448              solver->setColUpper(j,saveUpper[j]);
1449          }
1450          }
1451        //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1452        //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersDown,
1453        //     choice[i].downMovement,choice[i].finishedDown,choice[i].numIntInfeasDown,
1454        //     choice[i].numObjInfeasDown);
1455       
1456        // repeat the whole exercise, forcing the variable up
1457        if (!clp) {
1458          bool feasible=true;
1459          // If odd branching then maybe just one possibility
1460          if(choice[i].possibleBranch->numberBranchesLeft()>0) {
1461            choice[i].possibleBranch->branch();
1462            if (checkFeasibility) {
1463              // check branching did not make infeasible
1464              int iColumn;
1465              int numberColumns = solver->getNumCols();
1466              const double * columnLower = solver->getColLower();
1467              const double * columnUpper = solver->getColUpper();
1468              for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1469                if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1470                  feasible=false;
1471              }
1472            }
1473          } else {
1474            // second branch infeasible
1475            feasible=false;
1476          }
1477          if (feasible) {
1478            solver->solveFromHotStart() ;
1479            numberStrongDone++;
1480            numberStrongIterations += solver->getIterationCount();
1481            /*
1482              We now have an estimate of objective degradation that we can use for strong
1483              branching. If we're over the cutoff, the variable is monotone up.
1484              If we actually made it to optimality, check for a solution, and if we have
1485              a good one, call setBestSolution to process it. Note that this may reduce the
1486              cutoff, so we check again to see if we can declare this variable monotone.
1487            */
1488            if (solver->isProvenOptimal())
1489              iStatus=0; // optimal
1490            else if (solver->isIterationLimitReached()
1491                     &&!solver->isDualObjectiveLimitReached())
1492              iStatus=2; // unknown
1493            else
1494              iStatus=1; // infeasible
1495            newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1496            choice[i].numItersUp = solver->getIterationCount();
1497          } else {
1498            iStatus=1; // infeasible
1499            newObjectiveValue = 1.0e100;
1500            choice[i].numItersDown = 0;
1501          }
1502        } else {
1503          iStatus = outputStuff[2*i+1];
1504          choice[i].numItersUp = outputStuff[2*numberStrong+2*i+1];
1505          numberStrongDone++;
1506          numberStrongIterations += choice[i].numItersUp;
1507          newObjectiveValue = objectiveValue+newLower[i];
1508          solver->setColSolution(outputSolution[2*i+1]);
1509        }
1510        objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
1511        if (!iStatus) {
1512          choice[i].finishedUp = true ;
1513          if (newObjectiveValue>=model->getCutoff()) {
1514            objectiveChange = 1.0e100; // say infeasible
1515            numberStrongInfeasible++;
1516          } else {
1517            // See if integer solution
1518            if (model->feasibleSolution(choice[i].numIntInfeasUp,
1519                                        choice[i].numObjInfeasUp)
1520                &&model->problemFeasibility()->feasible(model,-1)>=0) {
1521              model->setBestSolution(CBC_STRONGSOL,
1522                                     newObjectiveValue,
1523                                     solver->getColSolution()) ;
1524              // only needed for odd solvers
1525              newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1526              objectiveChange = CoinMax(newObjectiveValue-objectiveValue_,0.0) ;
1527              model->setLastHeuristic(NULL);
1528              model->incrementUsed(solver->getColSolution());
1529              if (newObjectiveValue >= model->getCutoff()) {    //  *new* cutoff
1530                objectiveChange = 1.0e100 ;
1531                numberStrongInfeasible++;
1532              }
1533            }
1534          }
1535        } else if (iStatus==1) {
1536          objectiveChange = 1.0e100 ;
1537          numberStrongInfeasible++;
1538        } else {
1539          // Can't say much as we did not finish
1540          choice[i].finishedUp = false ;
1541          numberUnfinished++;
1542        }
1543        choice[i].upMovement = objectiveChange ;
1544       
1545        // restore bounds
1546        if (!clp)
1547          { for (int j=0;j<numberColumns;j++) {
1548            if (saveLower[j] != lower[j])
1549              solver->setColLower(j,saveLower[j]);
1550            if (saveUpper[j] != upper[j])
1551              solver->setColUpper(j,saveUpper[j]);
1552          }
1553          }
1554       
1555        //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1556        //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersUp,
1557        //     choice[i].upMovement,choice[i].finishedUp,choice[i].numIntInfeasUp,
1558        //     choice[i].numObjInfeasUp);
1559       
1560        /*
1561          End of evaluation for this candidate variable. Possibilities are:
1562          * Both sides below cutoff; this variable is a candidate for branching.
1563          * Both sides infeasible or above the objective cutoff: no further action
1564          here. Break from the evaluation loop and assume the node will be purged
1565          by the caller.
1566          * One side below cutoff: Install the branch (i.e., fix the variable). Break
1567          from the evaluation loop and assume the node will be reoptimised by the
1568          caller.
1569        */
1570        // reset
1571        choice[i].possibleBranch->resetNumberBranchesLeft();
1572        if (choice[i].upMovement<1.0e100) {
1573          if(choice[i].downMovement<1.0e100) {
1574            // feasible - no action
1575          } else {
1576            // up feasible, down infeasible
1577            anyAction=-1;
1578            //printf("Down infeasible for choice %d sequence %d\n",i,
1579            // model->object(choice[i].objectNumber)->columnNumber());
1580            if (!solveAll) {
1581              choice[i].possibleBranch->way(1);
1582              choice[i].possibleBranch->branch();
1583              break;
1584            } else {
1585              choice[i].fix=1;
1586            }
1587          }
1588        } else {
1589          if(choice[i].downMovement<1.0e100) {
1590            // down feasible, up infeasible
1591            anyAction=-1;
1592            //printf("Up infeasible for choice %d sequence %d\n",i,
1593            // model->object(choice[i].objectNumber)->columnNumber());
1594            if (!solveAll) {
1595              choice[i].possibleBranch->way(-1);
1596              choice[i].possibleBranch->branch();
1597              break;
1598            } else {
1599              choice[i].fix=-1;
1600            }
1601          } else {
1602            // neither side feasible
1603            anyAction=-2;
1604            //printf("Both infeasible for choice %d sequence %d\n",i,
1605            // model->object(choice[i].objectNumber)->columnNumber());
1606            break;
1607          }
1608        }
1609        bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
1610                            model->getDblParam(CbcModel::CbcMaximumSeconds));
1611        if (hitMaxTime) {
1612          numberStrong=i+1;
1613          break;
1614        }
1615        }
1616      if (!clp) {
1617        // Delete the snapshot
1618        solver->unmarkHotStart();
1619      } else {
1620        delete [] newLower;
1621        delete [] newUpper;
1622        delete [] outputStuff;
1623        int i;
1624        for (i=0;i<2*numberStrong;i++)
1625          delete [] outputSolution[i];
1626        delete [] outputSolution;
1627      }
1628      solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
1629      // restore basis
1630      solver->setWarmStart(ws);
1631      // Unless infeasible we will carry on
1632      // But we could fix anyway
1633      if (anyAction==-1&&solveAll) {
1634        // apply and take off
1635        for (i = 0 ; i < numberStrong ; i++) {
1636          if (choice[i].fix) {
1637            choice[i].possibleBranch->way(choice[i].fix) ;
1638            choice[i].possibleBranch->branch() ;
1639          }
1640        }
1641        bool feasible=true;
1642        if (checkFeasibility) {
1643          // check branching did not make infeasible
1644          int iColumn;
1645          int numberColumns = solver->getNumCols();
1646          const double * columnLower = solver->getColLower();
1647          const double * columnUpper = solver->getColUpper();
1648          for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1649            if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1650              feasible=false;
1651          }
1652        }
1653        if (feasible) {
1654          // can do quick optimality check
1655          int easy=2;
1656          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
1657          solver->resolve() ;
1658          solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
1659          feasible = solver->isProvenOptimal();
1660        }
1661        if (feasible) {
1662          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1663          model->reserveCurrentSolution(saveSolution);
1664          memcpy(saveLower,solver->getColLower(),numberColumns*sizeof(double));
1665          memcpy(saveUpper,solver->getColUpper(),numberColumns*sizeof(double));
1666          // Clean up all candidates whih are fixed
1667          int numberLeft=0;
1668          for (i = 0 ; i < numberStrong ; i++) {
1669            CbcStrongInfo thisChoice = choice[i];
1670            choice[i].possibleBranch=NULL;
1671            const CbcObject * object = model->object(thisChoice.objectNumber);
1672            int preferredWay;
1673            double infeasibility = object->infeasibility(preferredWay);
1674            if (!infeasibility) {
1675              // take out
1676              delete thisChoice.possibleBranch;
1677            } else {
1678              choice[numberLeft++]=thisChoice;
1679            }
1680          }
1681          numberStrong=numberLeft;
1682          for (;i<maximumStrong;i++) {
1683            delete choice[i].possibleBranch;
1684            choice[i].possibleBranch=NULL;
1685          }
1686          // If all fixed then round again
1687          if (!numberLeft) {
1688            finished=false;
1689            numberStrong=0;
1690            saveNumberStrong=0;
1691            maximumStrong=1;
1692          } else {
1693            anyAction=0;
1694          }
1695          // If these two uncommented then different action
1696          anyAction=-1;
1697          finished=true;
1698          //printf("some fixed but continuing %d left\n",numberLeft);
1699        } else {
1700          anyAction=-2; // say infeasible
1701        }
1702      }
1703      delete ws;
1704      int numberNodes = model->getNodeCount();
1705      // update number of strong iterations etc
1706      model->incrementStrongInfo(numberStrongDone,numberStrongIterations,
1707                                 anyAction==-2 ? 0:numberStrongInfeasible,anyAction==-2);
1708     
1709      /*
1710        anyAction >= 0 indicates that strong branching didn't produce any monotone
1711        variables. Sift through the candidates for the best one.
1712       
1713        QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is
1714        local to this code block. Perhaps should be numberNodes_ from model?
1715        Unclear what this calculation is doing.
1716      */
1717      if (anyAction>=0) {
1718       
1719        // get average cost per iteration and assume stopped ones
1720        // would stop after 50% more iterations at average cost??? !!! ???
1721        double averageCostPerIteration=0.0;
1722        double totalNumberIterations=1.0;
1723        int smallestNumberInfeasibilities=INT_MAX;
1724        for (i=0;i<numberStrong;i++) {
1725          totalNumberIterations += choice[i].numItersDown +
1726            choice[i].numItersUp ;
1727          averageCostPerIteration += choice[i].downMovement +
1728            choice[i].upMovement;
1729          smallestNumberInfeasibilities= 
1730            CoinMin(CoinMin(choice[i].numIntInfeasDown ,
1731                            choice[i].numIntInfeasUp ),
1732                    smallestNumberInfeasibilities);
1733        }
1734        //if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities)
1735        //numberNodes=1000000; // switch off search for better solution
1736        numberNodes=1000000; // switch off anyway
1737        averageCostPerIteration /= totalNumberIterations;
1738        // all feasible - choose best bet
1739       
1740        // New method does all at once so it can be more sophisticated
1741        // in deciding how to balance actions.
1742        // But it does need arrays
1743        double * changeUp = new double [numberStrong];
1744        int * numberInfeasibilitiesUp = new int [numberStrong];
1745        double * changeDown = new double [numberStrong];
1746        int * numberInfeasibilitiesDown = new int [numberStrong];
1747        CbcBranchingObject ** objects = new CbcBranchingObject * [ numberStrong];
1748        for (i = 0 ; i < numberStrong ; i++) {
1749          int iColumn = choice[i].possibleBranch->variable() ;
1750          model->messageHandler()->message(CBC_STRONG,model->messages())
1751            << i << iColumn
1752            <<choice[i].downMovement<<choice[i].numIntInfeasDown
1753            <<choice[i].upMovement<<choice[i].numIntInfeasUp
1754            <<choice[i].possibleBranch->value()
1755            <<CoinMessageEol;
1756          changeUp[i]=choice[i].upMovement;
1757          numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
1758          changeDown[i]=choice[i].downMovement;
1759          numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
1760          objects[i] = choice[i].possibleBranch;
1761        }
1762        int whichObject = decision->bestBranch(objects,numberStrong,numberUnsatisfied_,
1763                                               changeUp,numberInfeasibilitiesUp,
1764                                               changeDown,numberInfeasibilitiesDown,
1765                                               objectiveValue_);
1766        // move branching object and make sure it will not be deleted
1767        if (whichObject>=0) {
1768          branch_ = objects[whichObject];
1769          if (model->messageHandler()->logLevel()>3) 
1770            printf("Choosing column %d\n",choice[whichObject].possibleBranch->variable()) ;
1771          choice[whichObject].possibleBranch=NULL;
1772        }
1773        delete [] changeUp;
1774        delete [] numberInfeasibilitiesUp;
1775        delete [] changeDown;
1776        delete [] numberInfeasibilitiesDown;
1777        delete [] objects;
1778      }
1779      if (osiclp&&!allNormal) {
1780        // back to normal
1781        osiclp->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
1782      }
1783    }
1784    /*
1785      Simple branching. Probably just one, but we may have got here
1786      because of an odd branch e.g. a cut
1787    */
1788    else {
1789      // not strong
1790      // C) create branching object
1791      branch_ = choice[bestChoice].possibleBranch;
1792      choice[bestChoice].possibleBranch=NULL;
1793    }
1794  }
1795  // Set guessed solution value
1796  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
1797/*
1798  Cleanup, then we're outta here.
1799*/
1800  if (!model->branchingMethod())
1801    delete decision;
1802   
1803  for (i=0;i<maximumStrong;i++)
1804    delete choice[i].possibleBranch;
1805  delete [] choice;
1806  delete [] saveLower;
1807  delete [] saveUpper;
1808 
1809  // restore solution
1810  solver->setColSolution(saveSolution);
1811  delete [] saveSolution;
1812  return anyAction;
1813}
1814
1815/*
1816  Version for dynamic pseudo costs.
1817
1818  **** For now just return if anything odd
1819  later allow even if odd
1820
1821  The routine scans through the object list of the model looking for objects
1822  that indicate infeasibility. It tests each object using strong branching
1823  and selects the one with the least objective degradation.  A corresponding
1824  branching object is left attached to lastNode.
1825  This version gives preference in evaluation to variables which
1826  have not been evaluated many times.  It also uses numberStrong
1827  to say give up if last few tries have not changed incumbent.
1828  See Achterberg, Koch and Martin.
1829
1830  If strong branching is disabled, a candidate object is chosen essentially
1831  at random (whatever object ends up in pos'n 0 of the candidate array).
1832
1833  If a branching candidate is found to be monotone, bounds are set to fix the
1834  variable and the routine immediately returns (the caller is expected to
1835  reoptimize).
1836
1837  If a branching candidate is found to result in infeasibility in both
1838  directions, the routine immediately returns an indication of infeasibility.
1839
1840  Returns:  0   both branch directions are feasible
1841           -1   branching variable is monotone
1842           -2   infeasible
1843           -3   Use another method
1844
1845           For now just fix on objective from strong branching.
1846*/
1847
1848int CbcNode::chooseDynamicBranch (CbcModel *model, CbcNode *lastNode,
1849                                  OsiSolverBranch * & branches,int numberPassesLeft)
1850
1851{ if (lastNode)
1852    depth_ = lastNode->depth_+1;
1853  else
1854    depth_ = 0;
1855  delete branch_;
1856  branch_=NULL;
1857  OsiSolverInterface * solver = model->solver();
1858  // get information on solver type
1859  const OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
1860  const OsiBabSolver * auxiliaryInfo = dynamic_cast<const OsiBabSolver *> (auxInfo);
1861  //assert(objectiveValue_ == solver->getObjSense()*solver->getObjValue());
1862  double cutoff =model->getCutoff();
1863  double distanceToCutoff=cutoff-objectiveValue_;
1864  const double * lower = solver->getColLower();
1865  const double * upper = solver->getColUpper();
1866  // See what user thinks
1867  int anyAction=model->problemFeasibility()->feasible(model,0);
1868  if (anyAction) {
1869    // will return -2 if infeasible , 0 if treat as integer
1870    return anyAction-1;
1871  }
1872  int i;
1873  int saveStateOfSearch = model->stateOfSearch();
1874  int numberStrong=model->numberStrong();
1875  if (!auxiliaryInfo->warmStart())
1876    numberStrong=0;
1877  // But make more likely to get out after some times
1878  int changeStrategy=numberStrong;
1879  double changeFactor=1.0;
1880  // Use minimum of this and one stored in objects
1881  //int numberBeforeTrust = model->numberBeforeTrust();
1882  int numberObjects = model->numberObjects();
1883  bool checkFeasibility = numberObjects>model->numberIntegers();
1884  // For now return if not simple
1885  if (checkFeasibility)
1886    return -3;
1887  // Return if doing hot start (in BAB sense)
1888  if (model->hotstartSolution()) 
1889    return -3;
1890#define RANGING
1891#ifdef RANGING
1892  // Pass number
1893  int kPass=0;
1894  int numberRows = solver->getNumRows();
1895#endif
1896  int numberColumns = model->getNumCols();
1897  double * saveUpper = new double[numberColumns];
1898  double * saveLower = new double[numberColumns];
1899
1900  // Save solution in case heuristics need good solution later
1901 
1902  double * saveSolution = new double[numberColumns];
1903  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1904  model->reserveCurrentSolution(saveSolution);
1905  /*
1906    Get a branching decision object. Use the default dynamic decision criteria unless
1907    the user has loaded a decision method into the model.
1908  */
1909  CbcBranchDecision *decision = model->branchingMethod();
1910  if (!decision)
1911    decision = new CbcBranchDynamicDecision();
1912  int numberMini=0;
1913  int xPen=0;
1914  int xMark=0;
1915  for (i=0;i<numberColumns;i++) {
1916    saveLower[i] = lower[i];
1917    saveUpper[i] = upper[i];
1918  }
1919  // Get arrays to sort
1920  double * sort = new double[numberObjects];
1921  int * whichObject = new int[numberObjects];
1922  int * objectMark = new int[2*numberObjects+1];
1923  // Arrays with movements
1924  double * upEstimate = new double[numberObjects];
1925  double * downEstimate = new double[numberObjects];
1926  CbcStrongInfo * fixObject = new CbcStrongInfo[numberObjects];
1927  double estimatedDegradation=0.0; 
1928  int numberNodes=model->getNodeCount();
1929  int numberBeforeTrust = model->numberBeforeTrust();
1930  int numberPenalties = model->numberPenalties();
1931  if (numberBeforeTrust>=1000000) {
1932    numberBeforeTrust = numberBeforeTrust % 1000000;
1933    numberPenalties=0;
1934  } else if (numberBeforeTrust<0) {
1935    numberPenalties=numberColumns;
1936    numberBeforeTrust=0;
1937  }
1938  // May go round twice if strong branching fixes all local candidates
1939  bool finished=false;
1940  int numberToFix=0;
1941  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
1942  int saveClpOptions=0;
1943  if (osiclp) {
1944    // for faster hot start
1945    saveClpOptions = osiclp->specialOptions();
1946    osiclp->setSpecialOptions(saveClpOptions|1024);
1947  }
1948  int saveSearchStrategy2 = model->searchStrategy();
1949  if (saveSearchStrategy2<999) {
1950    // Get average up and down costs
1951    double averageUp=0.0;
1952    double averageDown=0.0;
1953    {
1954      int numberUp=0;
1955      int numberDown=0;
1956      int i;
1957      for ( i=0;i<numberObjects;i++) {
1958        CbcObject * object = model->modifiableObject(i);
1959        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
1960          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
1961        assert(dynamicObject);
1962        if (dynamicObject->numberTimesUp()) {
1963          numberUp++;
1964          averageUp += dynamicObject->upDynamicPseudoCost();
1965        }
1966        if (dynamicObject->numberTimesDown()) {
1967          numberDown++;
1968          averageDown += dynamicObject->downDynamicPseudoCost();
1969        }
1970      }
1971      if (numberUp) 
1972        averageUp /= (double) numberUp;
1973      else
1974        averageUp=1.0;
1975      if (numberDown) 
1976        averageDown /= (double) numberDown;
1977      else
1978        averageDown=1.0;
1979      for ( i=0;i<numberObjects;i++) {
1980        CbcObject * object = model->modifiableObject(i);
1981        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
1982          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
1983        assert(dynamicObject);
1984        if (!dynamicObject->numberTimesUp()) 
1985          dynamicObject->setUpDynamicPseudoCost(averageUp);
1986      if (!dynamicObject->numberTimesDown()) 
1987        dynamicObject->setDownDynamicPseudoCost(averageDown);
1988      }
1989    }
1990  } else if (saveSearchStrategy2<1999) {
1991    // pseudo shadow prices
1992    model->pseudoShadow(NULL,NULL);
1993  } else if (saveSearchStrategy2<2999) {
1994    // leave old ones
1995  } else if (saveSearchStrategy2<3999) {
1996    // pseudo shadow prices at root
1997    if (!numberNodes)
1998      model->pseudoShadow(NULL,NULL);
1999  } else {
2000    abort();
2001  }
2002  if (saveSearchStrategy2>=0)
2003    saveSearchStrategy2 = saveSearchStrategy2 % 1000;
2004  if (saveSearchStrategy2==999)
2005    saveSearchStrategy2=-1;
2006  int px[4]={-1,-1,-1,-1};
2007  int saveSearchStrategy = saveSearchStrategy2<99 ? saveSearchStrategy2 : saveSearchStrategy2-100;
2008  bool newWay = saveSearchStrategy2>98;
2009  int numberNotTrusted=0;
2010  int numberStrongDone;
2011  int numberUnfinished;
2012  int numberStrongInfeasible;
2013  int numberStrongIterations;
2014  while(!finished) {
2015    finished=true;
2016    decision->initialize(model);
2017    // Some objects may compute an estimate of best solution from here
2018    estimatedDegradation=0.0; 
2019    numberToFix=0;
2020    int numberIntegerInfeasibilities=0; // without odd ones
2021    int numberToDo=0;
2022    int iBestNot=-1;
2023    int iBestGot=-1;
2024    double best=0.0;
2025    numberNotTrusted=0;
2026    numberStrongDone=0;
2027    numberUnfinished=0;
2028    numberStrongInfeasible=0;
2029    numberStrongIterations=0;
2030    int * which = objectMark+numberObjects+1;
2031    int neededPenalties;
2032    // We may go round this loop three times (only if we think we have solution)
2033    for (int iPass=0;iPass<3;iPass++) {
2034     
2035      // compute current state
2036      int numberObjectInfeasibilities; // just odd ones
2037      model->feasibleSolution(
2038                              numberIntegerInfeasibilities,
2039                              numberObjectInfeasibilities);
2040     
2041      // Some objects may compute an estimate of best solution from here
2042      estimatedDegradation=0.0; 
2043      numberUnsatisfied_ = 0;
2044      int bestPriority=INT_MAX;
2045      /*
2046        Scan for branching objects that indicate infeasibility. Choose candidates
2047        using priority as the first criteria, then integer infeasibility.
2048       
2049        The algorithm is to fill the array with a set of good candidates (by
2050        infeasibility) with priority bestPriority.  Finding a candidate with
2051        priority better (less) than bestPriority flushes the choice array. (This
2052        serves as initialization when the first candidate is found.)
2053       
2054      */
2055      numberToDo=0;
2056      neededPenalties=0;
2057      iBestNot=-1;
2058      double bestNot=0.0;
2059      iBestGot=-1;
2060      best=0.0;
2061#define PRINT_STUFF -1
2062      for (i=0;i<numberObjects;i++) {
2063        CbcObject * object = model->modifiableObject(i);
2064        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2065          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2066        assert(dynamicObject);
2067        int preferredWay;
2068        double infeasibility = object->infeasibility(preferredWay);
2069        int priorityLevel = object->priority();
2070#define ZERO_ONE 0
2071#define ZERO_FAKE 1.0e20;
2072#if ZERO_ONE==1
2073        // branch on 0-1 first (temp)
2074        if (fabs(saveSolution[dynamicObject->columnNumber()])<1.0)
2075          priorityLevel--;
2076#endif
2077#if ZERO_ONE==2
2078        if (fabs(saveSolution[dynamicObject->columnNumber()])<1.0)
2079          infeasibility *= ZERO_FAKE;
2080#endif
2081        if (infeasibility) {
2082          int iColumn = dynamicObject->columnNumber();
2083          double gap = saveUpper[iColumn]-saveLower[iColumn];
2084          // Give precedence to ones with gap of 1.0
2085          assert(gap>0.0);
2086          infeasibility /= CoinMin(gap,100.0);
2087          if (!depth_&&false) {
2088            // try closest to 0.5
2089            double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2090            infeasibility = fabs(0.5-part);
2091          }
2092          bool gotDown=false;
2093          int numberThisDown = dynamicObject->numberTimesDown();
2094          if (numberThisDown>=numberBeforeTrust)
2095            gotDown=true;
2096          bool gotUp=false;
2097          int numberThisUp = dynamicObject->numberTimesUp();
2098          if (numberThisUp>=numberBeforeTrust)
2099            gotUp=true;
2100          if ((numberNodes%PRINT_STUFF)==0&&PRINT_STUFF>0)
2101            printf("%d down %d %g up %d %g - infeas %g\n",
2102                   i,numberThisDown,object->downEstimate(),numberThisUp,object->upEstimate(),
2103                   infeasibility);
2104          // Increase estimated degradation to solution
2105          estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate());
2106          downEstimate[i]=object->downEstimate();
2107          upEstimate[i]=object->upEstimate();
2108          numberUnsatisfied_++;
2109          // Better priority? Flush choices.
2110          if (priorityLevel<bestPriority) {
2111            numberToDo=0;
2112            bestPriority = priorityLevel;
2113            iBestGot=-1;
2114            best=0.0;
2115            numberNotTrusted=0;
2116          } else if (priorityLevel>bestPriority) {
2117            continue;
2118          }
2119          if (!gotUp||!gotDown)
2120            numberNotTrusted++;
2121          // Check for suitability based on infeasibility.
2122          if ((gotDown&&gotUp)&&numberStrong>0) {
2123            sort[numberToDo]=-infeasibility;
2124            if (infeasibility>best) {
2125              best=infeasibility;
2126              iBestGot=numberToDo;
2127            }
2128          } else {
2129            objectMark[neededPenalties]=numberToDo;
2130            which[neededPenalties++]=dynamicObject->columnNumber();
2131            int iColumn = dynamicObject->columnNumber();
2132            double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
2133            sort[numberToDo]=-10.0*infeasibility;
2134            if (!(numberThisUp+numberThisDown))
2135              sort[numberToDo] *= 100.0; // make even more likely
2136            if (1.0-fabs(part-0.5)>bestNot) {
2137              iBestNot=numberToDo;
2138              bestNot = 1.0-fabs(part-0.5);
2139            }
2140          }
2141          whichObject[numberToDo++]=i;
2142        } else {
2143          // for debug
2144          downEstimate[i]=-1.0;
2145          upEstimate[i]=-1.0;
2146        }
2147      }
2148      if (numberUnsatisfied_) {
2149        // some infeasibilities - go to next steps
2150        break;
2151      } else if (!iPass) {
2152        // may just need resolve
2153        solver->resolve();
2154        memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
2155        model->reserveCurrentSolution(saveSolution);
2156        if (!solver->isProvenOptimal()) {
2157          // infeasible
2158          anyAction=-2;
2159          break;
2160        }
2161      } else if (iPass==1) {
2162        // looks like a solution - get paranoid
2163        bool roundAgain=false;
2164        // get basis
2165        CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
2166        if (!ws)
2167          break;
2168        double tolerance;
2169        solver->getDblParam(OsiPrimalTolerance,tolerance);
2170        for (i=0;i<numberColumns;i++) {
2171          double value = saveSolution[i];
2172          if (value<lower[i]-tolerance) {
2173            saveSolution[i]=lower[i];
2174            roundAgain=true;
2175            ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
2176          } else if (value>upper[i]+tolerance) {
2177            saveSolution[i]=upper[i];
2178            roundAgain=true;
2179            ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
2180          } 
2181        }
2182        if (roundAgain) {
2183          // restore basis
2184          solver->setWarmStart(ws);
2185          solver->setColSolution(saveSolution);
2186          delete ws;
2187          bool takeHint;
2188          OsiHintStrength strength;
2189          solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
2190          solver->setHintParam(OsiDoDualInResolve,false,OsiHintDo) ;
2191          solver->resolve();
2192          solver->setHintParam(OsiDoDualInResolve,takeHint,strength) ;
2193          memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
2194          model->reserveCurrentSolution(saveSolution);
2195          if (!solver->isProvenOptimal()) {
2196            // infeasible
2197            anyAction=-2;
2198            break;
2199          }
2200        } else {
2201          delete ws;
2202          break;
2203        }
2204      }
2205    }
2206    if (anyAction==-2) {
2207      break;
2208    }
2209    bool solveAll=false; // set true to say look at all even if some fixed (experiment)
2210    solveAll=true;
2211    // skip if solution
2212    if (!numberUnsatisfied_)
2213      break;
2214    //bool skipAll = (numberBeforeTrust>20&&numberNodes>20000&&numberNotTrusted==0);
2215    bool skipAll = numberNotTrusted==0;
2216    bool doneHotStart=false;
2217#ifndef CBC_WEAK_STRONG
2218    if ((numberNodes%20)==0||(model->specialOptions()&8)!=0)
2219      skipAll=false;
2220#endif
2221    int searchStrategy = saveSearchStrategy>=0 ? (saveSearchStrategy%10) : -1;
2222    if (!newWay) {
2223    // 10 up always use %10, 20 up as 10 and allow penalties
2224    // But adjust depending on ratio of iterations
2225    if (searchStrategy>0&&saveSearchStrategy<10) {
2226      if (numberBeforeTrust>=5&&numberBeforeTrust<=10) {
2227        if (searchStrategy!=2) {
2228          int numberIterations = model->getIterationCount();
2229          int numberStrongIterations = model->numberStrongIterations();
2230          if (numberStrongIterations>numberIterations+10000) {
2231            searchStrategy=2;
2232            //skipAll=true;
2233          } else if (numberStrongIterations*4+1000<numberIterations||depth_<5) {
2234            searchStrategy=3;
2235            skipAll=false;
2236          }
2237        } else {
2238          skipAll=true;
2239        }
2240      }
2241    }
2242    } else {
2243    // But adjust depending on ratio of iterations
2244    if (saveSearchStrategy<0) {
2245      // unset
2246      if ((numberNodes%20)==0||(model->specialOptions()&8)!=0) {
2247        // Do numberStrong
2248        searchStrategy=3;
2249      } else if (depth_<5) {
2250        // Do numberStrong
2251        searchStrategy=2;
2252      } else {
2253        int numberIterations = model->getIterationCount();
2254        int numberStrongIterations = model->numberStrongIterations();
2255        int numberRows = solver->getNumRows();
2256        if (numberStrongIterations>numberIterations+CoinMin(10000,10*numberRows)) {
2257          // off
2258          searchStrategy=0;
2259        } else if (numberStrongIterations*4+1000<numberIterations) {
2260          // Do numberStrong if not trusted
2261          searchStrategy=2;
2262        } else {
2263          searchStrategy=1;
2264        }
2265      }
2266    }
2267    if (searchStrategy<3&&(!numberNotTrusted||!searchStrategy))
2268      skipAll=true;
2269    else
2270      skipAll=false;
2271    }
2272    // worth trying if too many times
2273    // Save basis
2274    CoinWarmStart * ws = NULL;
2275    // save limit
2276    int saveLimit=0;
2277    solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
2278    if (!skipAll) {
2279      ws = solver->getWarmStart();
2280      int limit=100;
2281#if 0
2282      int averageBranchIterations = model->getIterationCount()/(model->getNodeCount()+1);
2283      if (numberNodes)
2284        limit = CoinMin(CoinMax(limit,2*averageBranchIterations),500);
2285      else
2286        limit = 500;
2287#endif
2288      if ((!saveStateOfSearch||searchStrategy>3)&&saveLimit<limit&&saveLimit==100)
2289        solver->setIntParam(OsiMaxNumIterationHotStart,limit); 
2290    }
2291    // Say which one will be best
2292    int whichChoice=0;
2293    int bestChoice;
2294    if (iBestGot>=0)
2295      bestChoice=iBestGot;
2296    else
2297      bestChoice=iBestNot;
2298    assert (bestChoice>=0);
2299    // If we have hit max time don't do strong branching
2300    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
2301                        model->getDblParam(CbcModel::CbcMaximumSeconds));
2302    // also give up if we are looping round too much
2303    if (hitMaxTime||numberPassesLeft<=0||(!numberNotTrusted&&false)) {
2304      int iObject = whichObject[bestChoice];
2305      CbcObject * object = model->modifiableObject(iObject);
2306      int preferredWay;
2307      object->infeasibility(preferredWay);
2308      branch_=object->createBranch(preferredWay);
2309      branch_->way(preferredWay);
2310      delete ws;
2311      ws=NULL;
2312      break;
2313    } else {
2314      // say do fast
2315      int easy=1;
2316      if (!skipAll)
2317        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
2318      int iDo;
2319#ifdef RANGING
2320      if ((skipAll&&numberBeforeTrust&&saveSearchStrategy<20)||saveSearchStrategy<10)
2321        numberPenalties=0;
2322      {
2323        // off penalties if too much
2324        double needed = neededPenalties;
2325        needed *= numberRows;
2326        if (needed>1.0e6&&numberNodes&&saveSearchStrategy<20) {
2327          numberPenalties=0;
2328          neededPenalties=0;
2329        }
2330      }
2331      if (osiclp&&numberPenalties&&neededPenalties) {
2332        xPen += neededPenalties;
2333        which--;
2334        which[0]=neededPenalties;
2335        osiclp->passInRanges(which);
2336        // Mark hot start and get ranges
2337        if (kPass) {
2338          // until can work out why solution can go funny
2339          int save = osiclp->specialOptions();
2340          osiclp->setSpecialOptions(save|256);
2341          solver->markHotStart();
2342          osiclp->setSpecialOptions(save);
2343        } else {
2344          solver->markHotStart();
2345        }
2346        assert (auxiliaryInfo->warmStart());
2347        doneHotStart=true;
2348        xMark++;
2349        kPass++;
2350        osiclp->passInRanges(NULL);
2351        const double * downCost=osiclp->upRange();
2352        const double * upCost=osiclp->downRange();
2353        //printf("numberTodo %d needed %d numberpenalties %d\n",numberToDo,neededPenalties,numberPenalties);
2354        double invTrust = 1.0/((double) numberBeforeTrust);
2355        for (int i=0;i<neededPenalties;i++) {
2356          int j = objectMark[i];
2357          int iObject = whichObject[j];
2358          CbcObject * object = model->modifiableObject(iObject);
2359          CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2360            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2361          int iSequence=dynamicObject->columnNumber();
2362          double value = saveSolution[iSequence];
2363          value -= floor(value);
2364          double upPenalty = CoinMin(upCost[i],1.0e110)*(1.0-value);
2365          double downPenalty = CoinMin(downCost[i],1.0e110)*value;
2366          if (!numberBeforeTrust) {
2367            // override
2368            downEstimate[iObject]=downPenalty;
2369            upEstimate[iObject]=upPenalty;
2370          } else {
2371            int numberThisDown = dynamicObject->numberTimesDown();
2372            if (numberThisDown<numberBeforeTrust) {
2373              double fraction = ((double) numberThisDown)*invTrust;
2374              downEstimate[iObject] = fraction*downEstimate[iObject]+(1.0-fraction)*downPenalty;
2375            }
2376            int numberThisUp = dynamicObject->numberTimesUp();
2377            if (numberThisUp<numberBeforeTrust) {
2378              double fraction = ((double) numberThisUp)*invTrust;
2379              upEstimate[iObject] = fraction*upEstimate[iObject]+(1.0-fraction)*upPenalty;
2380            }
2381          }
2382          sort[j] = - CoinMin(downEstimate[iObject],upEstimate[iObject]);
2383#ifdef CBC_WEAK_STRONG
2384          sort[j] -= 1.0e10; // make more likely to be chosen
2385#endif
2386          //if ((numberNodes%PRINT_STUFF)==0&&PRINT_STUFF>0)
2387          if (!numberNodes)
2388            printf("%d pen down ps %g -> %g up ps %g -> %g\n",
2389                   iObject,downCost[i],downPenalty,upCost[i],upPenalty);
2390        }
2391      } else {
2392        if (!skipAll) {
2393          // Mark hot start
2394          solver->markHotStart();
2395          doneHotStart=true;
2396          assert (auxiliaryInfo->warmStart());
2397          xMark++;
2398          //if (solver->isProvenPrimalInfeasible())
2399          //printf("**** Hot start says node infeasible\n");
2400        }
2401        // make sure best will be first
2402        if (iBestGot>=0)
2403          sort[iBestGot]=-1.0e120;
2404      }
2405#else
2406      if (!skipAll) {
2407        // Mark hot start
2408        doneHotStart=true;
2409        assert (auxiliaryInfo->warmStart());
2410        solver->markHotStart();
2411        xMark++;
2412      }
2413      // make sure best will be first
2414      if (iBestGot>=0)
2415        sort[iBestGot]=-COIN_DBL_MAX;
2416#endif
2417      // Actions 0 - exit for repeat, 1 resolve and try old choice,2 exit for continue
2418#define ACTION 0
2419#if ACTION<2
2420      if (anyAction)
2421        numberToDo=0; // skip as we will be trying again
2422#endif
2423      // Sort
2424      CoinSort_2(sort,sort+numberToDo,whichObject);
2425      // Change in objective opposite infeasible
2426      double worstFeasible=0.0;
2427      // Just first if strong off
2428      if (!numberStrong)
2429        numberToDo=CoinMin(numberToDo,1);
2430      iDo=0;
2431      int saveLimit2;
2432      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit2);
2433      bool doQuickly = false; // numberToDo>2*numberStrong;
2434      //printf("todo %d, strong %d\n",numberToDo,numberStrong);
2435      int numberTest=numberNotTrusted>0 ? numberStrong : (numberStrong+1)/2;
2436      int numberTest2 = 2*numberStrong;
2437      if (!newWay) {
2438      if (searchStrategy==3) {
2439        // Previously decided we need strong
2440        doQuickly=false;
2441        numberTest = numberStrong;
2442        //numberTest2 = 1000000;
2443      }
2444      if (searchStrategy<0||searchStrategy==1)
2445        //numberTest2 = 1000000;
2446#if 0
2447      if (numberBeforeTrust>20&&(numberNodes>20000||(numberNodes>200&&numberNotTrusted==0))) {
2448        if ((numberNodes%20)!=0) {
2449          numberTest=0;
2450          doQuickly=true;
2451        }
2452      }
2453#else
2454      // Try nearly always off
2455      if (searchStrategy<2) {
2456        if ((numberNodes%20)!=0) {
2457          if ((model->specialOptions()&8)==0) {
2458            numberTest=0;
2459            doQuickly=true;
2460          }
2461        } else {
2462          doQuickly=false;
2463          numberTest=2*numberStrong;
2464        }
2465      } else if (searchStrategy!=3) {
2466        doQuickly=true;
2467        numberTest=numberStrong;
2468      }
2469#endif
2470      if (depth_<10&&numberStrong) {
2471        if (searchStrategy!=2) {
2472          doQuickly=false;
2473          if (depth_<7)
2474            numberStrong *=3;
2475          if (!depth_) 
2476            numberStrong=CoinMin(6*numberStrong,numberToDo);
2477          numberTest=numberStrong;
2478        }
2479        model->setStateOfSearch(2); // use min min
2480      }
2481      // could adjust using average iterations per branch
2482      // double average = ((double)model->getIterationCount())/
2483      //((double) model->getNodeCount()+1.0);
2484      // if too many and big then just do 10 its
2485      if (!skipAll&&saveStateOfSearch) {
2486        //if (numberNotTrusted>3*numberStrong&&numberRows>250&&numberColumns>1000&&saveLimit==100)
2487          // off solver->setIntParam(OsiMaxNumIterationHotStart,10);
2488      }
2489      // make negative for test
2490      distanceToCutoff = - distanceToCutoff;
2491      if (numberObjects>-100) {
2492        // larger
2493        distanceToCutoff *= 100.0;
2494      }
2495        distanceToCutoff = -COIN_DBL_MAX;
2496      // Do at least 5 strong
2497      if (numberColumns<1000&&(depth_<15||numberNodes<1000000))
2498        numberTest = CoinMax(numberTest,5);
2499      if ((model->specialOptions()&8)==0) {
2500        if (skipAll) {
2501          numberTest=0;
2502          doQuickly=true;
2503        }
2504      } else {
2505        // do 5 as strong is fixing
2506        numberTest = CoinMax(numberTest,5);
2507      }
2508      } else {
2509      int numberTest=numberNotTrusted>0 ? numberStrong : (numberStrong+1)/2;
2510      int numberTest2 = 2*numberStrong;
2511      if (searchStrategy>=3) {
2512        // Previously decided we need strong
2513        doQuickly=false;
2514        if (depth_<7)
2515          numberStrong *=3;
2516        if (!depth_) 
2517          numberStrong=CoinMin(6*numberStrong,numberToDo);
2518        numberTest = numberStrong;
2519        numberTest2 *= 2;
2520      } else if (searchStrategy==2||(searchStrategy==1&&depth_<6)) {
2521        numberStrong *=2;
2522        if (!depth_) 
2523          numberStrong=CoinMin(2*numberStrong,numberToDo);
2524        numberTest = numberStrong;
2525      } else if (searchStrategy==1&&numberNotTrusted) {
2526        numberTest = numberStrong;
2527      } else {
2528        numberTest=0;
2529        skipAll=true;
2530      }
2531      distanceToCutoff=model->getCutoff()-objectiveValue_;
2532      // make negative for test
2533      distanceToCutoff = - distanceToCutoff;
2534      if (numberObjects>-100) {
2535        // larger
2536        distanceToCutoff *= 100.0;
2537      }
2538      distanceToCutoff = -COIN_DBL_MAX;
2539      if (skipAll) {
2540        numberTest=0;
2541        doQuickly=true;
2542      }
2543      }
2544      px[0]=numberTest;
2545      px[1]=numberTest2;
2546      px[2]= doQuickly ? 1 : -1;
2547      px[3]=numberStrong;
2548      //printf("skipAll %c doQuickly %c numberTest %d numberTest2 %d numberNot %d\n",
2549      //     skipAll ? 'Y' : 'N',doQuickly ? 'Y' : 'N',numberTest,numberTest2,numberNotTrusted);
2550      // See if we want mini tree
2551      bool wantMiniTree=false;
2552      if (model->sizeMiniTree()&&depth_>7&&saveStateOfSearch>0)
2553        wantMiniTree=true;
2554      numberMini=0;
2555      for ( iDo=0;iDo<numberToDo;iDo++) {
2556        CbcStrongInfo choice;
2557        int iObject = whichObject[iDo];
2558        CbcObject * object = model->modifiableObject(iObject);
2559        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2560          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2561        int iColumn = dynamicObject->columnNumber();
2562        int preferredWay;
2563        object->infeasibility(preferredWay);
2564        choice.possibleBranch=object->createBranch(preferredWay);
2565        // Save which object it was
2566        choice.objectNumber=iObject;
2567        choice.numIntInfeasUp = numberUnsatisfied_;
2568        choice.numIntInfeasDown = numberUnsatisfied_;
2569        choice.upMovement = upEstimate[iObject];
2570        choice.downMovement = downEstimate[iObject];
2571        assert (choice.upMovement>=0.0);
2572        assert (choice.downMovement>=0.0);
2573        choice.fix=0; // say not fixed
2574        // see if can skip strong branching
2575        int canSkip = choice.possibleBranch->fillStrongInfo(choice);
2576        if (!newWay) {
2577        if (!doQuickly||(numberTest>0&&searchStrategy!=2))
2578          canSkip=0;
2579        } else {
2580        if (skipAll)
2581          canSkip=1;
2582        else if (numberTest>0&&searchStrategy>=3)
2583          canSkip=0;
2584        }
2585        if (!numberBeforeTrust) {
2586          canSkip=1;
2587        }
2588        if (sort[iDo]<distanceToCutoff)
2589          canSkip=0;
2590        if (((numberTest2<=0&&numberTest<=0)||skipAll)&&sort[iDo]>distanceToCutoff) {
2591          canSkip=1; // always skip
2592          if (iDo>20)
2593            break; // give up anyway
2594        }
2595        if (model->messageHandler()->logLevel()>3&&numberBeforeTrust) 
2596          dynamicObject->print(1,choice.possibleBranch->value());
2597        // was if (!canSkip)
2598        if (newWay)
2599        numberTest2--;
2600        if (!canSkip) {
2601          numberTest--;
2602          if (!newWay)
2603          numberTest2--;
2604          // just do a few
2605          //if (canSkip)
2606          //solver->setIntParam(OsiMaxNumIterationHotStart,10);
2607          double objectiveChange ;
2608          double newObjectiveValue=1.0e100;
2609          int j;
2610          // status is 0 finished, 1 infeasible and other
2611          int iStatus;
2612          /*
2613            Try the down direction first. (Specify the initial branching alternative as
2614            down with a call to way(-1). Each subsequent call to branch() performs the
2615            specified branch and advances the branch object state to the next branch
2616            alternative.)
2617          */
2618          choice.possibleBranch->way(-1) ;
2619          decision->saveBranchingObject( choice.possibleBranch);
2620          choice.possibleBranch->branch() ;
2621          solver->solveFromHotStart() ;
2622          bool needHotStartUpdate=false;
2623          numberStrongDone++;
2624          numberStrongIterations += solver->getIterationCount();
2625          /*
2626            We now have an estimate of objective degradation that we can use for strong
2627            branching. If we're over the cutoff, the variable is monotone up.
2628            If we actually made it to optimality, check for a solution, and if we have
2629            a good one, call setBestSolution to process it. Note that this may reduce the
2630            cutoff, so we check again to see if we can declare this variable monotone.
2631          */
2632          if (solver->isProvenOptimal())
2633            iStatus=0; // optimal
2634          else if (solver->isIterationLimitReached()
2635                   &&!solver->isDualObjectiveLimitReached())
2636            iStatus=2; // unknown
2637          else
2638            iStatus=1; // infeasible
2639          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2640          choice.numItersDown = solver->getIterationCount();
2641          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
2642          decision->updateInformation( solver,this);
2643          if (!iStatus) {
2644            choice.finishedDown = true ;
2645            if (newObjectiveValue>=cutoff) {
2646              objectiveChange = 1.0e100; // say infeasible
2647              numberStrongInfeasible++;
2648            } else {
2649              // See if integer solution
2650              if (model->feasibleSolution(choice.numIntInfeasDown,
2651                                          choice.numObjInfeasDown)
2652                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
2653                if (auxiliaryInfo->solutionAddsCuts()) {
2654                  needHotStartUpdate=true;
2655                  solver->unmarkHotStart();
2656                }
2657                model->setBestSolution(CBC_STRONGSOL,
2658                                       newObjectiveValue,
2659                                       solver->getColSolution()) ;
2660                if (needHotStartUpdate) {
2661                  solver->resolve();
2662                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2663                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
2664                  model->feasibleSolution(choice.numIntInfeasDown,
2665                                          choice.numObjInfeasDown);
2666                }
2667                model->setLastHeuristic(NULL);
2668                model->incrementUsed(solver->getColSolution());
2669                cutoff =model->getCutoff();
2670                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
2671                  objectiveChange = 1.0e100 ;
2672                  numberStrongInfeasible++;
2673                }
2674              }
2675            }
2676          } else if (iStatus==1) {
2677            objectiveChange = 1.0e100 ;
2678            numberStrongInfeasible++;
2679          } else {
2680            // Can't say much as we did not finish
2681            choice.finishedDown = false ;
2682            numberUnfinished++;
2683          }
2684          choice.downMovement = objectiveChange ;
2685         
2686          // restore bounds
2687          for ( j=0;j<numberColumns;j++) {
2688            if (saveLower[j] != lower[j])
2689              solver->setColLower(j,saveLower[j]);
2690            if (saveUpper[j] != upper[j])
2691              solver->setColUpper(j,saveUpper[j]);
2692          }
2693          if(needHotStartUpdate) {
2694            needHotStartUpdate = false;
2695            solver->resolve();
2696            //we may again have an integer feasible solution
2697            int numberIntegerInfeasibilities;
2698            int numberObjectInfeasibilities;
2699            if (model->feasibleSolution(
2700                                        numberIntegerInfeasibilities,
2701                                        numberObjectInfeasibilities)) {
2702              double objValue = solver->getObjValue();
2703              model->setBestSolution(CBC_STRONGSOL,
2704                                     objValue,
2705                                     solver->getColSolution()) ;
2706              solver->resolve();
2707              cutoff =model->getCutoff();
2708            }
2709            solver->markHotStart();
2710          }
2711          //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2712          //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2713          //     choice.objectNumber,iStatus,newObjectiveValue,choice.numItersDown,
2714          //     choice.downMovement,choice.finishedDown,choice.numIntInfeasDown,
2715          //     choice.numObjInfeasDown);
2716         
2717          // repeat the whole exercise, forcing the variable up
2718          decision->saveBranchingObject( choice.possibleBranch);
2719          choice.possibleBranch->branch();
2720          solver->solveFromHotStart() ;
2721          numberStrongDone++;
2722          numberStrongIterations += solver->getIterationCount();
2723          /*
2724            We now have an estimate of objective degradation that we can use for strong
2725            branching. If we're over the cutoff, the variable is monotone up.
2726            If we actually made it to optimality, check for a solution, and if we have
2727            a good one, call setBestSolution to process it. Note that this may reduce the
2728            cutoff, so we check again to see if we can declare this variable monotone.
2729          */
2730          if (solver->isProvenOptimal())
2731            iStatus=0; // optimal
2732          else if (solver->isIterationLimitReached()
2733                   &&!solver->isDualObjectiveLimitReached())
2734            iStatus=2; // unknown
2735          else
2736            iStatus=1; // infeasible
2737          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2738          choice.numItersUp = solver->getIterationCount();
2739          objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
2740          decision->updateInformation( solver,this);
2741          if (!iStatus) {
2742            choice.finishedUp = true ;
2743            if (newObjectiveValue>=cutoff) {
2744              objectiveChange = 1.0e100; // say infeasible
2745              numberStrongInfeasible++;
2746            } else {
2747              // See if integer solution
2748              if (model->feasibleSolution(choice.numIntInfeasUp,
2749                                          choice.numObjInfeasUp)
2750                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
2751                if (auxiliaryInfo->solutionAddsCuts()) {
2752                  needHotStartUpdate=true;
2753                  solver->unmarkHotStart();
2754                }
2755                model->setBestSolution(CBC_STRONGSOL,
2756                                       newObjectiveValue,
2757                                       solver->getColSolution()) ;
2758                if (needHotStartUpdate) {
2759                  solver->resolve();
2760                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
2761                  objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_,0.0);
2762                  model->feasibleSolution(choice.numIntInfeasDown,
2763                                          choice.numObjInfeasDown);
2764                }
2765                model->setLastHeuristic(NULL);
2766                model->incrementUsed(solver->getColSolution());
2767                cutoff =model->getCutoff();
2768                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
2769                  objectiveChange = 1.0e100 ;
2770                  numberStrongInfeasible++;
2771                }
2772              }
2773            }
2774          } else if (iStatus==1) {
2775            objectiveChange = 1.0e100 ;
2776            numberStrongInfeasible++;
2777          } else {
2778            // Can't say much as we did not finish
2779            choice.finishedUp = false ;
2780            numberUnfinished++;
2781          }
2782          choice.upMovement = objectiveChange ;
2783         
2784          // restore bounds
2785          for ( j=0;j<numberColumns;j++) {
2786            if (saveLower[j] != lower[j])
2787              solver->setColLower(j,saveLower[j]);
2788            if (saveUpper[j] != upper[j])
2789              solver->setColUpper(j,saveUpper[j]);
2790          }
2791          if(needHotStartUpdate) {
2792            needHotStartUpdate = false;
2793            solver->resolve();
2794            //we may again have an integer feasible solution
2795            int numberIntegerInfeasibilities;
2796            int numberObjectInfeasibilities;
2797            if (model->feasibleSolution(
2798                                        numberIntegerInfeasibilities,
2799                                        numberObjectInfeasibilities)) {
2800              double objValue = solver->getObjValue();
2801              model->setBestSolution(CBC_STRONGSOL,
2802                                     objValue,
2803                                     solver->getColSolution()) ;
2804              solver->resolve();
2805              cutoff =model->getCutoff();
2806            }
2807            solver->markHotStart();
2808          }
2809         
2810          //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2811          //     choice.objectNumber,iStatus,newObjectiveValue,choice.numItersUp,
2812          //     choice.upMovement,choice.finishedUp,choice.numIntInfeasUp,
2813          //     choice.numObjInfeasUp);
2814        }
2815   
2816        solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit2); 
2817        /*
2818          End of evaluation for this candidate variable. Possibilities are:
2819          * Both sides below cutoff; this variable is a candidate for branching.
2820          * Both sides infeasible or above the objective cutoff: no further action
2821          here. Break from the evaluation loop and assume the node will be purged
2822          by the caller.
2823          * One side below cutoff: Install the branch (i.e., fix the variable). Break
2824          from the evaluation loop and assume the node will be reoptimised by the
2825          caller.
2826        */
2827        // reset
2828        choice.possibleBranch->resetNumberBranchesLeft();
2829        if (choice.upMovement<1.0e100) {
2830          if(choice.downMovement<1.0e100) {
2831            // In case solution coming in was odd
2832            choice.upMovement = CoinMax(0.0,choice.upMovement);
2833            choice.downMovement = CoinMax(0.0,choice.downMovement);
2834#if ZERO_ONE==2
2835            // branch on 0-1 first (temp)
2836            if (fabs(choice.possibleBranch->value())<1.0) {
2837              choice.upMovement *= ZERO_FAKE;
2838              choice.downMovement *= ZERO_FAKE;
2839            }
2840#endif
2841            // feasible - see which best
2842            if (!canSkip) {
2843              if (iColumn==-46) {
2844                printf("sort %g downest %g upest %g ",sort[iDo],downEstimate[iObject],
2845                     upEstimate[iObject]);
2846                printf("downMove %g upMove %g value %g current pseudo %g %g\n",
2847                       choice.downMovement,choice.upMovement,choice.possibleBranch->value(),
2848                       dynamicObject->downDynamicPseudoCost(),dynamicObject->upDynamicPseudoCost());
2849              }
2850              if (model->messageHandler()->logLevel()>3) 
2851                printf("sort %g downest %g upest %g ",sort[iDo],downEstimate[iObject],
2852                     upEstimate[iObject]);
2853              model->messageHandler()->message(CBC_STRONG,model->messages())
2854                << iObject << iColumn
2855                <<choice.downMovement<<choice.numIntInfeasDown
2856                <<choice.upMovement<<choice.numIntInfeasUp
2857                <<choice.possibleBranch->value()
2858                <<CoinMessageEol;
2859            }
2860            //if (!stateOfSearch)
2861            //choice.numIntInfeasDown=99999; // temp fudge
2862            if (wantMiniTree)
2863              decision->setBestCriterion(-1.0);
2864            double bestCriterion = -1.0;
2865            double gap = saveUpper[iColumn]-saveLower[iColumn];
2866            // Give precedence to ones with gap of 1.0
2867            assert(gap>0.0);
2868            double factor = changeFactor/CoinMin(gap,100.0);
2869            int betterWay = decision->betterBranch(choice.possibleBranch,
2870                                                   branch_,
2871                                                   choice.upMovement*factor,
2872                                                   choice.numIntInfeasUp ,
2873                                                   choice.downMovement*factor,
2874                                                   choice.numIntInfeasDown );
2875            if (wantMiniTree) {
2876              double criterion = decision->getBestCriterion();
2877              sort[numberMini]=-criterion;
2878              whichObject[numberMini++]=whichObject[iDo];
2879              assert (betterWay);
2880              if (criterion>bestCriterion) 
2881                bestCriterion=criterion;
2882              else
2883                betterWay=0;
2884            }
2885            if (iDo>=changeStrategy) {
2886              // make less likely
2887              changeStrategy+=numberStrong;
2888              changeFactor *= 0.9;
2889            }
2890            if (betterWay) {
2891              delete branch_;
2892              // C) create branching object
2893              branch_ = choice.possibleBranch;
2894              choice.possibleBranch=NULL;
2895              branch_->way(betterWay);
2896              bestChoice = choice.objectNumber;
2897              whichChoice = iDo;
2898              if (numberStrong<=1) {
2899                delete ws;
2900                ws=NULL;
2901                break;
2902              }
2903            } else {
2904              delete choice.possibleBranch;
2905              if (iDo>=2*numberStrong) {
2906                delete ws;
2907                ws=NULL;
2908                break;
2909              }
2910              if (!dynamicObject||dynamicObject->numberTimesUp()>1) {
2911                if (iDo-whichChoice>=numberStrong)
2912                  break; // give up
2913              } else {
2914                if (iDo-whichChoice>=2*numberStrong) {
2915                  delete ws;
2916                  ws=NULL;
2917                  break; // give up
2918                }
2919              }
2920            }
2921          } else {
2922            // up feasible, down infeasible
2923            anyAction=-1;
2924            worstFeasible = CoinMax(worstFeasible,choice.upMovement);
2925            //printf("Down infeasible for choice %d sequence %d\n",i,
2926            // model->object(choice.objectNumber)->columnNumber());
2927            if (!solveAll) {
2928              choice.possibleBranch->way(1);
2929              choice.possibleBranch->branch();
2930              delete choice.possibleBranch;
2931              delete ws;
2932              ws=NULL;
2933              break;
2934            } else {
2935              choice.fix=1;
2936              fixObject[numberToFix++]=choice;
2937#define FIXNOW
2938#ifdef FIXNOW
2939              double value = ceil(saveSolution[iColumn]);
2940              saveLower[iColumn]=value;
2941              solver->setColLower(iColumn,value);
2942              assert(doneHotStart);
2943              solver->unmarkHotStart();
2944              solver->markHotStart();
2945#endif
2946            }
2947          }
2948        } else {
2949          if(choice.downMovement<1.0e100) {
2950            // down feasible, up infeasible
2951            anyAction=-1;
2952            worstFeasible = CoinMax(worstFeasible,choice.downMovement);
2953            //printf("Up infeasible for choice %d sequence %d\n",i,
2954            // model->object(choice.objectNumber)->columnNumber());
2955            if (!solveAll) {
2956              choice.possibleBranch->way(-1);
2957              choice.possibleBranch->branch();
2958              delete choice.possibleBranch;
2959              delete ws;
2960              ws=NULL;
2961              break;
2962            } else {
2963              choice.fix=-1;
2964              fixObject[numberToFix++]=choice;
2965#ifdef FIXNOW
2966              double value = floor(saveSolution[iColumn]);
2967              saveUpper[iColumn]=value;
2968              solver->setColUpper(iColumn,value);
2969              assert(doneHotStart);
2970              solver->unmarkHotStart();
2971              solver->markHotStart();
2972#endif
2973            }
2974          } else {
2975            // neither side feasible
2976            anyAction=-2;
2977            delete choice.possibleBranch;
2978            //printf("Both infeasible for choice %d sequence %d\n",i,
2979            // model->object(choice.objectNumber)->columnNumber());
2980            delete ws;
2981            ws=NULL;
2982            break;
2983          }
2984        }
2985        // Check max time
2986        hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
2987                       model->getDblParam(CbcModel::CbcMaximumSeconds));
2988        if (hitMaxTime) {
2989          // make sure rest are fast
2990          doQuickly=true;
2991          for ( int jDo=iDo+1;jDo<numberToDo;jDo++) {
2992            int iObject = whichObject[iDo];
2993            CbcObject * object = model->modifiableObject(iObject);
2994            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2995              dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2996            dynamicObject->setNumberBeforeTrust(0);
2997          }
2998          numberTest=0;
2999          distanceToCutoff=-COIN_DBL_MAX;
3000        }
3001      }
3002      double averageChange = model->sumChangeObjective()/((double) model->getNodeCount());
3003      if (depth_<10||worstFeasible>0.2*averageChange) 
3004        solveAll=false;
3005      if (model->messageHandler()->logLevel()>3||false) { 
3006        if (anyAction==-2)
3007          printf("infeasible\n");
3008        else if(anyAction==-1)
3009          if (!solveAll)
3010            printf("%d fixed\n",numberToFix);
3011          else
3012            printf("%d fixed AND choosing %d iDo %d iChosenWhen %d numberToDo %d\n",numberToFix,bestChoice,
3013                   iDo,whichChoice,numberToDo);
3014        else
3015          printf("choosing %d iDo %d iChosenWhen %d numberToDo %d\n",bestChoice,
3016                 iDo,whichChoice,numberToDo);
3017      }
3018      if (doneHotStart) {
3019        // Delete the snapshot
3020        solver->unmarkHotStart();
3021        // back to normal
3022        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3023        // restore basis
3024        solver->setWarmStart(ws);
3025      }
3026      solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
3027      // Unless infeasible we will carry on
3028      // But we could fix anyway
3029      if (numberToFix&&!hitMaxTime) {
3030        if (anyAction==-2) {
3031          // take off
3032          for (i = 0 ; i < numberToFix ; i++) {
3033            delete fixObject[i].possibleBranch;
3034          }
3035        } else {
3036          // apply and take off
3037          for (i = 0 ; i < numberToFix ; i++) {
3038#ifndef FIXNOW
3039            fixObject[i].possibleBranch->way(fixObject[i].fix) ;
3040            fixObject[i].possibleBranch->branch() ;
3041#endif
3042            delete fixObject[i].possibleBranch;
3043          }
3044          bool feasible=true;
3045#if ACTION <2
3046          if (solveAll) {
3047            // can do quick optimality check
3048            int easy=2;
3049            solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
3050            solver->resolve() ;
3051            solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3052            feasible = solver->isProvenOptimal();
3053            if (feasible) {
3054              anyAction=0;
3055              numberMini=0;
3056              memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
3057              model->reserveCurrentSolution(saveSolution);
3058              memcpy(saveLower,solver->getColLower(),numberColumns*sizeof(double));
3059              memcpy(saveUpper,solver->getColUpper(),numberColumns*sizeof(double));
3060              model->setPointers(solver);
3061              // See if candidate still possible
3062              if (branch_) {
3063                const CbcObject * object = model->object(bestChoice);
3064                int preferredWay;
3065                double infeasibility = object->infeasibility(preferredWay);
3066                if (!infeasibility) {
3067                  // take out
3068                  delete branch_;
3069                  branch_=NULL;
3070                } else {
3071                  branch_->way(preferredWay);
3072                }
3073              }
3074            } else {
3075              anyAction=-2;
3076              finished=true;
3077            }
3078          }
3079#endif
3080          // If  fixed then round again
3081          if (!branch_&&anyAction!=-2) {
3082            finished=false;
3083          }
3084          // If these in then different action
3085#if ACTION == 1
3086          if (!anyAction)
3087            anyAction=-1;
3088          finished=true;
3089#endif
3090        }
3091      }
3092      delete ws;
3093    }
3094  }
3095  if (model->messageHandler()->logLevel()>2) 
3096    printf("%d strong, %d iters, %d pen, %d mark, %d fixed, action %d nnott %d nt %d, %d dq %s ns %d\n",
3097         numberStrongDone,numberStrongIterations,xPen,xMark,
3098           numberToFix,anyAction,numberNotTrusted,px[0],px[1],px[2]>0 ? "y" : "n",px[3]);
3099  // update number of strong iterations etc
3100  model->incrementStrongInfo(numberStrongDone,numberStrongIterations,
3101                             anyAction==-2 ? 0:numberToFix,anyAction==-2);
3102  if (!newWay) {
3103  if (((model->searchStrategy()+1)%1000)==0) {
3104    if (solver->messageHandler()->logLevel()>1)
3105      printf("%d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3106             numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3107             numberNotTrusted);
3108    // decide what to do
3109    int strategy=1;
3110    if (numberUnfinished*4>numberStrongDone&&numberStrongInfeasible*10<numberStrongDone) {
3111      strategy=2;
3112      if (model->logLevel()>1)
3113        printf("going to strategy 2\n");
3114    }
3115    if (numberNodes)
3116      strategy=1;  // should only happen after hot start
3117    model->setSearchStrategy(strategy);
3118  }
3119  }
3120  //if (numberToFix&&depth_<5)
3121  //printf("%d fixed by strong at depth %d\n",numberToFix,depth_);
3122  // Set guessed solution value
3123  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
3124 
3125  // Get collection of branches if mini tree wanted
3126  if (anyAction==0&&numberMini&&numberMini>1) {
3127    // Sort
3128    CoinSort_2(sort,sort+numberMini,whichObject);
3129    delete branch_;
3130    branch_=NULL;
3131    numberMini = CoinMin(numberMini,model->sizeMiniTree());
3132    anyAction=numberMini;
3133    branches = new OsiSolverBranch[numberMini];
3134    for (int iDo=0;iDo<numberMini;iDo++) {
3135      int iObject = whichObject[iDo];
3136      CbcObject * object = model->modifiableObject(iObject);
3137      OsiSolverBranch * oneBranch = object->solverBranch();
3138      branches[iDo]=*oneBranch;
3139      delete oneBranch;
3140    }
3141  }
3142/*
3143  Cleanup, then we're finished
3144*/
3145  if (!model->branchingMethod())
3146    delete decision;
3147   
3148  delete [] fixObject;
3149  delete [] sort;
3150  delete [] whichObject;
3151  delete [] objectMark;
3152  delete [] saveLower;
3153  delete [] saveUpper;
3154  delete [] upEstimate;
3155  delete [] downEstimate;
3156  if (osiclp) 
3157    osiclp->setSpecialOptions(saveClpOptions);
3158 
3159  // restore solution
3160  solver->setColSolution(saveSolution);
3161  model->reserveCurrentSolution(saveSolution);
3162  delete [] saveSolution;
3163  model->setStateOfSearch(saveStateOfSearch);
3164  return anyAction;
3165}
3166int CbcNode::analyze (CbcModel *model, double * results)
3167{
3168  int i;
3169  int numberIterationsAllowed = model->numberAnalyzeIterations();
3170  OsiSolverInterface * solver = model->solver();
3171  objectiveValue_ = solver->getObjSense()*solver->getObjValue();
3172  double cutoff =model->getCutoff();
3173  const double * lower = solver->getColLower();
3174  const double * upper = solver->getColUpper();
3175  const double * dj = solver->getReducedCost();
3176  int numberObjects = model->numberObjects();
3177  int numberColumns = model->getNumCols();
3178  // Initialize arrays
3179  int numberIntegers = model->numberIntegers();
3180  int * back = new int[numberColumns];
3181  const int * integerVariable = model->integerVariable();
3182  for (i=0;i<numberColumns;i++) 
3183    back[i]=-1;
3184  // What results is
3185  double * newLower = results;
3186  double * objLower = newLower+numberIntegers;
3187  double * newUpper = objLower+numberIntegers;
3188  double * objUpper = newUpper+numberIntegers;
3189  for (i=0;i<numberIntegers;i++) {
3190    int iColumn = integerVariable[i];
3191    back[iColumn]=i;
3192    newLower[i]=0.0;
3193    objLower[i]=-COIN_DBL_MAX;
3194    newUpper[i]=0.0;
3195    objUpper[i]=-COIN_DBL_MAX;
3196  }
3197  double * saveUpper = new double[numberColumns];
3198  double * saveLower = new double[numberColumns];
3199  int anyAction=0;
3200  // Save solution in case heuristics need good solution later
3201 
3202  double * saveSolution = new double[numberColumns];
3203  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
3204  model->reserveCurrentSolution(saveSolution);
3205  for (i=0;i<numberColumns;i++) {
3206    saveLower[i] = lower[i];
3207    saveUpper[i] = upper[i];
3208  }
3209  // Get arrays to sort
3210  double * sort = new double[numberObjects];
3211  int * whichObject = new int[numberObjects];
3212  int numberToFix=0;
3213  int numberToDo=0;
3214     
3215  // compute current state
3216  int numberObjectInfeasibilities; // just odd ones
3217  int numberIntegerInfeasibilities;
3218  model->feasibleSolution(
3219                          numberIntegerInfeasibilities,
3220                          numberObjectInfeasibilities);
3221     
3222  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
3223  int saveClpOptions=0;
3224  bool fastIterations = (model->specialOptions()&8)!=0;
3225  if (osiclp&&fastIterations) {
3226    // for faster hot start
3227    saveClpOptions = osiclp->specialOptions();
3228    osiclp->setSpecialOptions(saveClpOptions|1024);
3229  }
3230  /*
3231    Scan for branching objects that indicate infeasibility. Choose candidates
3232    using priority as the first criteria, then integer infeasibility.
3233   
3234    The algorithm is to fill the array with a set of good candidates (by
3235    infeasibility) with priority bestPriority.  Finding a candidate with
3236    priority better (less) than bestPriority flushes the choice array. (This
3237    serves as initialization when the first candidate is found.)
3238   
3239  */
3240  numberToDo=0;
3241  for (i=0;i<numberObjects;i++) {
3242    CbcObject * object = model->modifiableObject(i);
3243    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3244      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3245    if(!dynamicObject)
3246      continue;
3247    int preferredWay;
3248    double infeasibility = object->infeasibility(preferredWay);
3249    int iColumn = dynamicObject->columnNumber();
3250    if (saveUpper[iColumn]==saveLower[iColumn])
3251      continue;
3252    if (infeasibility)
3253      sort[numberToDo]=-1.0e10-infeasibility;
3254    else
3255      sort[numberToDo]=-fabs(dj[iColumn]);
3256    whichObject[numberToDo++]=i;
3257  }
3258  // Save basis
3259  CoinWarmStart * ws = solver->getWarmStart();
3260  int saveLimit;
3261  solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
3262  int targetIterations = CoinMax(500,numberIterationsAllowed/numberObjects);
3263  if (saveLimit<targetIterations)
3264    solver->setIntParam(OsiMaxNumIterationHotStart,targetIterations); 
3265  // Mark hot start
3266  solver->markHotStart();
3267  // Sort
3268  CoinSort_2(sort,sort+numberToDo,whichObject);
3269  //double distanceToCutoff=model->getCutoff()-objectiveValue_;
3270  double * currentSolution = model->currentSolution();
3271  double integerTolerance = 
3272    model->getDblParam(CbcModel::CbcIntegerTolerance);
3273  double objMin = 1.0e50;
3274  double objMax = -1.0e50;
3275  bool needResolve=false;
3276  int iDo;
3277  for (iDo=0;iDo<numberToDo;iDo++) {
3278    CbcStrongInfo choice;
3279    int iObject = whichObject[iDo];
3280    CbcObject * object = model->modifiableObject(iObject);
3281    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3282      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3283    int iColumn = dynamicObject->columnNumber();
3284    int preferredWay;
3285    object->infeasibility(preferredWay);
3286    double value = currentSolution[iColumn];
3287    double nearest = floor(value+0.5);
3288    double lowerValue = floor(value);
3289    bool satisfied=false;
3290    if (fabs(value-nearest)<=integerTolerance||value<saveLower[iColumn]||value>saveUpper[iColumn]) {
3291      satisfied=true;
3292      double newValue;
3293      if (nearest<saveUpper[iColumn]) {
3294        newValue = nearest + 1.0001*integerTolerance;
3295        lowerValue = nearest;
3296      } else {
3297        newValue = nearest - 1.0001*integerTolerance;
3298        lowerValue = nearest-1;
3299      }
3300      currentSolution[iColumn]=newValue;
3301    }
3302    double upperValue = lowerValue+1.0;
3303    choice.possibleBranch=object->createBranch(preferredWay);
3304    currentSolution[iColumn]=value;
3305    // Save which object it was
3306    choice.objectNumber=iObject;
3307    choice.numIntInfeasUp = numberUnsatisfied_;
3308    choice.numIntInfeasDown = numberUnsatisfied_;
3309    choice.downMovement = 0.0;
3310    choice.upMovement = 0.0;
3311    choice.numItersDown = 0;
3312    choice.numItersUp = 0;
3313    choice.fix=0; // say not fixed
3314    double objectiveChange ;
3315    double newObjectiveValue=1.0e100;
3316    int j;
3317    // status is 0 finished, 1 infeasible and other
3318    int iStatus;
3319    /*
3320      Try the down direction first. (Specify the initial branching alternative as
3321      down with a call to way(-1). Each subsequent call to branch() performs the
3322      specified branch and advances the branch object state to the next branch
3323      alternative.)
3324    */
3325    choice.possibleBranch->way(-1) ;
3326    choice.possibleBranch->branch() ;
3327    if (fabs(value-lowerValue)>integerTolerance) {
3328      solver->solveFromHotStart() ;
3329      /*
3330        We now have an estimate of objective degradation that we can use for strong
3331        branching. If we're over the cutoff, the variable is monotone up.
3332        If we actually made it to optimality, check for a solution, and if we have
3333        a good one, call setBestSolution to process it. Note that this may reduce the
3334        cutoff, so we check again to see if we can declare this variable monotone.
3335      */
3336      if (solver->isProvenOptimal())
3337        iStatus=0; // optimal
3338      else if (solver->isIterationLimitReached()
3339               &&!solver->isDualObjectiveLimitReached())
3340        iStatus=2; // unknown
3341      else
3342        iStatus=1; // infeasible
3343      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3344      choice.numItersDown = solver->getIterationCount();
3345      numberIterationsAllowed -= choice.numItersDown;
3346      objectiveChange = newObjectiveValue  - objectiveValue_;
3347      if (!iStatus) {
3348        choice.finishedDown = true ;
3349        if (newObjectiveValue>=cutoff) {
3350          objectiveChange = 1.0e100; // say infeasible
3351        } else {
3352          // See if integer solution
3353          if (model->feasibleSolution(choice.numIntInfeasDown,
3354                                      choice.numObjInfeasDown)
3355              &&model->problemFeasibility()->feasible(model,-1)>=0) {
3356            model->setBestSolution(CBC_STRONGSOL,
3357                                   newObjectiveValue,
3358                                   solver->getColSolution()) ;
3359            model->setLastHeuristic(NULL);
3360            model->incrementUsed(solver->getColSolution());
3361            cutoff =model->getCutoff();
3362            if (newObjectiveValue >= cutoff)    //  *new* cutoff
3363              objectiveChange = 1.0e100 ;
3364          }
3365        }
3366      } else if (iStatus==1) {
3367        objectiveChange = 1.0e100 ;
3368      } else {
3369        // Can't say much as we did not finish
3370        choice.finishedDown = false ;
3371      }
3372      choice.downMovement = objectiveChange ;
3373    }
3374    // restore bounds
3375    for ( j=0;j<numberColumns;j++) {
3376      if (saveLower[j] != lower[j])
3377        solver->setColLower(j,saveLower[j]);
3378      if (saveUpper[j] != upper[j])
3379        solver->setColUpper(j,saveUpper[j]);
3380    }
3381    // repeat the whole exercise, forcing the variable up
3382    choice.possibleBranch->branch();
3383    if (fabs(value-upperValue)>integerTolerance) {
3384      solver->solveFromHotStart() ;
3385      /*
3386        We now have an estimate of objective degradation that we can use for strong
3387        branching. If we're over the cutoff, the variable is monotone up.
3388        If we actually made it to optimality, check for a solution, and if we have
3389        a good one, call setBestSolution to process it. Note that this may reduce the
3390        cutoff, so we check again to see if we can declare this variable monotone.
3391      */
3392      if (solver->isProvenOptimal())
3393        iStatus=0; // optimal
3394      else if (solver->isIterationLimitReached()
3395               &&!solver->isDualObjectiveLimitReached())
3396        iStatus=2; // unknown
3397      else
3398        iStatus=1; // infeasible
3399      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
3400      choice.numItersUp = solver->getIterationCount();
3401      numberIterationsAllowed -= choice.numItersUp;
3402      objectiveChange = newObjectiveValue  - objectiveValue_;
3403      if (!iStatus) {
3404        choice.finishedUp = true ;
3405        if (newObjectiveValue>=cutoff) {
3406          objectiveChange = 1.0e100; // say infeasible
3407        } else {
3408          // See if integer solution
3409          if (model->feasibleSolution(choice.numIntInfeasUp,
3410                                      choice.numObjInfeasUp)
3411              &&model->problemFeasibility()->feasible(model,-1)>=0) {
3412            model->setBestSolution(CBC_STRONGSOL,
3413                                   newObjectiveValue,
3414                                   solver->getColSolution()) ;
3415            model->setLastHeuristic(NULL);
3416            model->incrementUsed(solver->getColSolution());
3417            cutoff =model->getCutoff();
3418            if (newObjectiveValue >= cutoff)    //  *new* cutoff
3419              objectiveChange = 1.0e100 ;
3420          }
3421        }
3422      } else if (iStatus==1) {
3423        objectiveChange = 1.0e100 ;
3424      } else {
3425        // Can't say much as we did not finish
3426        choice.finishedUp = false ;
3427      }
3428      choice.upMovement = objectiveChange ;
3429     
3430      // restore bounds
3431      for ( j=0;j<numberColumns;j++) {
3432        if (saveLower[j] != lower[j])
3433          solver->setColLower(j,saveLower[j]);
3434        if (saveUpper[j] != upper[j])
3435          solver->setColUpper(j,saveUpper[j]);
3436      }
3437    }
3438    // If objective goes above certain amount we can set bound
3439    int jInt = back[iColumn];
3440    newLower[jInt]=upperValue;
3441    if (choice.finishedDown||!fastIterations)
3442      objLower[jInt]=choice.downMovement+objectiveValue_;
3443    else
3444      objLower[jInt]=objectiveValue_;
3445    newUpper[jInt]=lowerValue;
3446    if (choice.finishedUp||!fastIterations)
3447      objUpper[jInt]=choice.upMovement+objectiveValue_;
3448    else
3449      objUpper[jInt]=objectiveValue_;
3450    objMin = CoinMin(CoinMin(objLower[jInt],objUpper[jInt]),objMin);
3451    /*
3452      End of evaluation for this candidate variable. Possibilities are:
3453      * Both sides below cutoff; this variable is a candidate for branching.
3454      * Both sides infeasible or above the objective cutoff: no further action
3455      here. Break from the evaluation loop and assume the node will be purged
3456      by the caller.
3457      * One side below cutoff: Install the branch (i.e., fix the variable). Break
3458      from the evaluation loop and assume the node will be reoptimised by the
3459      caller.
3460    */
3461    if (choice.upMovement<1.0e100) {
3462      if(choice.downMovement<1.0e100) {
3463        objMax = CoinMax(CoinMax(objLower[jInt],objUpper[jInt]),objMax);
3464        // In case solution coming in was odd
3465        choice.upMovement = CoinMax(0.0,choice.upMovement);
3466        choice.downMovement = CoinMax(0.0,choice.downMovement);
3467        // feasible -
3468        model->messageHandler()->message(CBC_STRONG,model->messages())
3469          << iObject << iColumn
3470          <<choice.downMovement<<choice.numIntInfeasDown
3471          <<choice.upMovement<<choice.numIntInfeasUp
3472          <<value
3473          <<CoinMessageEol;
3474      } else {
3475        // up feasible, down infeasible
3476        anyAction=-1;
3477        if (!satisfied)
3478          needResolve=true;
3479        choice.fix=1;
3480        numberToFix++;
3481        saveLower[iColumn]=upperValue;
3482        solver->setColLower(iColumn,upperValue);
3483      }
3484    } else {
3485      if(choice.downMovement<1.0e100) {
3486        // down feasible, up infeasible
3487        anyAction=-1;
3488        if (!satisfied)
3489          needResolve=true;
3490        choice.fix=-1;
3491        numberToFix++;
3492        saveUpper[iColumn]=lowerValue;
3493        solver->setColUpper(iColumn,lowerValue);
3494      } else {
3495        // neither side feasible
3496        anyAction=-2;
3497        printf("Both infeasible for choice %d sequence %d\n",i,
3498               model->object(choice.objectNumber)->columnNumber());
3499        delete ws;
3500        ws=NULL;
3501        solver->writeMps("bad");
3502        numberToFix=-1;
3503        break;
3504      }
3505    }
3506    delete choice.possibleBranch;
3507    if (numberIterationsAllowed<=0)
3508      break;
3509    //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
3510    //     choice.downMovement,choice.upMovement,value);
3511  }
3512  printf("Best possible solution %g, can fix more if solution of %g found - looked at %d variables in %d iterations\n",
3513         objMin,objMax,iDo,model->numberAnalyzeIterations()-numberIterationsAllowed);
3514  model->setNumberAnalyzeIterations(numberIterationsAllowed);
3515  // Delete the snapshot
3516  solver->unmarkHotStart();
3517  // back to normal
3518  solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3519  solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
3520  // restore basis
3521  solver->setWarmStart(ws);
3522  delete ws;
3523   
3524  delete [] sort;
3525  delete [] whichObject;
3526  delete [] saveLower;
3527  delete [] saveUpper;
3528  delete [] back;
3529  // restore solution
3530  solver->setColSolution(saveSolution);
3531  if (osiclp) 
3532    osiclp->setSpecialOptions(saveClpOptions);
3533  model->reserveCurrentSolution(saveSolution);
3534  delete [] saveSolution;
3535  if (needResolve)
3536    solver->resolve();
3537  return numberToFix;
3538}
3539
3540
3541CbcNode::CbcNode(const CbcNode & rhs) 
3542{ 
3543#ifdef CHECK_NODE
3544  printf("CbcNode %x Constructor from rhs %x\n",this,&rhs);
3545#endif
3546  if (rhs.nodeInfo_)
3547    nodeInfo_ = rhs.nodeInfo_->clone();
3548  else
3549    nodeInfo_=NULL;
3550  objectiveValue_=rhs.objectiveValue_;
3551  guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
3552  if (rhs.branch_)
3553    branch_=rhs.branch_->clone();
3554  else
3555    branch_=NULL;
3556  depth_ = rhs.depth_;
3557  numberUnsatisfied_ = rhs.numberUnsatisfied_;
3558}
3559
3560CbcNode &
3561CbcNode::operator=(const CbcNode & rhs)
3562{
3563  if (this != &rhs) {
3564    delete nodeInfo_;
3565    if (nodeInfo_)
3566      nodeInfo_ = rhs.nodeInfo_->clone();
3567    else
3568      nodeInfo_ = NULL;
3569    objectiveValue_=rhs.objectiveValue_;
3570    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
3571    if (rhs.branch_)
3572      branch_=rhs.branch_->clone();
3573    else
3574      branch_=NULL,
3575    depth_ = rhs.depth_;
3576    numberUnsatisfied_ = rhs.numberUnsatisfied_;
3577  }
3578  return *this;
3579}
3580
3581
3582CbcNode::~CbcNode ()
3583{
3584#ifdef CHECK_NODE
3585  if (nodeInfo_) 
3586    printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
3587         this,nodeInfo_,nodeInfo_->numberPointingToThis());
3588  else
3589    printf("CbcNode %x Destructor nodeInfo %x (?)\n",
3590         this,nodeInfo_);
3591#endif
3592  if (nodeInfo_) {
3593    nodeInfo_->nullOwner();
3594    int numberToDelete=nodeInfo_->numberBranchesLeft();
3595    //    CbcNodeInfo * parent = nodeInfo_->parent();
3596    if (nodeInfo_->decrement(numberToDelete)==0) {
3597      delete nodeInfo_;
3598    } else {
3599      //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,parent);
3600      // anyway decrement parent
3601      //if (parent)
3602      ///parent->decrement(1);
3603    }
3604  }
3605  delete branch_;
3606}
3607// Decrement  active cut counts
3608void 
3609CbcNode::decrementCuts(int change)
3610{
3611  if(nodeInfo_) {
3612    nodeInfo_->decrementCuts(change);
3613  }
3614}
3615void 
3616CbcNode::decrementParentCuts(int change)
3617{
3618  if(nodeInfo_) {
3619    nodeInfo_->decrementParentCuts(change);
3620  }
3621}
3622
3623/*
3624  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
3625  in the attached nodeInfo_.
3626*/
3627void
3628CbcNode::initializeInfo()
3629{
3630  assert(nodeInfo_ && branch_) ;
3631  nodeInfo_->initializeInfo(branch_->numberBranches());
3632}
3633// Nulls out node info
3634void 
3635CbcNode::nullNodeInfo()
3636{
3637  nodeInfo_=NULL;
3638}
3639
3640int
3641CbcNode::branch()
3642{
3643  double changeInGuessed=branch_->branch(true);
3644  guessedObjectiveValue_+= changeInGuessed;
3645  //#define PRINTIT
3646#ifdef PRINTIT
3647  int numberLeft = nodeInfo_->numberBranchesLeft();
3648  CbcNodeInfo * parent = nodeInfo_->parent();
3649  int parentNodeNumber = -1;
3650  //CbcBranchingObject * object1 = branch_->object_;
3651  //CbcObject * object = object1->
3652  //int sequence = object->columnNumber);
3653  int id=-1;
3654  double value=0.0;
3655  if (branch_) {
3656    id = branch_->variable();
3657    value = branch_->value();
3658  }
3659  printf("id %d value %g objvalue %g\n",id,value,objectiveValue_);
3660  if (parent)
3661    parentNodeNumber = parent->nodeNumber();
3662  printf("Node number %d, %s, way %d, depth %d, parent node number %d\n",
3663         nodeInfo_->nodeNumber(),(numberLeft==2) ? "leftBranch" : "rightBranch",
3664         way(),depth_,parentNodeNumber);
3665#endif
3666  return nodeInfo_->branchedOn();
3667}
Note: See TracBrowser for help on using the repository browser.