source: trunk/CbcNode.cpp @ 216

Last change on this file since 216 was 216, checked in by forrest, 15 years ago

add branching stuff

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