source: trunk/CbcNode.cpp @ 201

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

mostly NDEBUG

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