source: trunk/CbcNode.cpp @ 154

Last change on this file since 154 was 154, checked in by forrest, 16 years ago

back to old node count

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