source: trunk/CbcNode.cpp @ 165

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

CbcFeasibility?

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