source: trunk/CbcNode.cpp @ 44

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

try and make more robust

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 50.7 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <string>
8//#define CBC_DEBUG 1
9//#define CHECK_CUT_COUNTS
10#include <cassert>
11#include <cfloat>
12#define CUTS
13#include "OsiSolverInterface.hpp"
14#include "CoinWarmStartBasis.hpp"
15#include "CoinTime.hpp"
16#include "CbcModel.hpp"
17#include "CbcNode.hpp"
18#include "CbcBranchActual.hpp"
19#include "OsiRowCut.hpp"
20#include "OsiRowCutDebugger.hpp"
21#include "OsiCuts.hpp"
22#include "CbcCountRowCut.hpp"
23#include "CbcMessage.hpp"
24#include "OsiClpSolverInterface.hpp"
25using namespace std;
26#include "CglCutGenerator.hpp"
27// Default Constructor
28CbcNodeInfo::CbcNodeInfo ()
29  :
30  numberPointingToThis_(0),
31  parent_(NULL),
32  owner_(NULL),
33  numberCuts_(0),
34  cuts_(NULL),
35  numberRows_(0),
36  numberBranchesLeft_(0)
37{
38#ifdef CHECK_NODE
39  printf("CbcNodeInfo %x Constructor\n",this);
40#endif
41}
42// Constructor given parent
43CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent)
44  :
45  numberPointingToThis_(2),
46  parent_(parent),
47  owner_(NULL),
48  numberCuts_(0),
49  cuts_(NULL),
50  numberRows_(0),
51  numberBranchesLeft_(2)
52{
53#ifdef CHECK_NODE
54  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
55#endif
56  if (parent_) {
57    numberRows_ = parent_->numberRows_+parent_->numberCuts_;
58  }
59}
60// Copy Constructor
61CbcNodeInfo::CbcNodeInfo (const CbcNodeInfo & rhs)
62  :
63  numberPointingToThis_(rhs.numberPointingToThis_),
64  parent_(rhs.parent_),
65  owner_(rhs.owner_),
66  numberCuts_(rhs.numberCuts_),
67  cuts_(NULL),
68  numberRows_(rhs.numberRows_),
69  numberBranchesLeft_(rhs.numberBranchesLeft_)
70{
71#ifdef CHECK_NODE
72  printf("CbcNodeInfo %x Copy constructor\n",this);
73#endif
74  if (numberCuts_) {
75    cuts_ = new CbcCountRowCut * [numberCuts_];
76    int n=0;
77    for (int i=0;i<numberCuts_;i++) {
78      CbcCountRowCut * thisCut = rhs.cuts_[i];
79      if (thisCut) {
80        // I think this is correct - new one should take priority
81        thisCut->setInfo(this,n);
82        thisCut->increment(numberBranchesLeft_); 
83        cuts_[n++] = thisCut;
84      }
85    }
86    numberCuts_=n;
87  }
88}
89// Constructor given parent and owner
90CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent, CbcNode * owner)
91  :
92  numberPointingToThis_(2),
93  parent_(parent),
94  owner_(owner),
95  numberCuts_(0),
96  cuts_(NULL),
97  numberRows_(0),
98  numberBranchesLeft_(2)
99{
100#ifdef CHECK_NODE
101  printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_);
102#endif
103  if (parent_) {
104    numberRows_ = parent_->numberRows_+parent_->numberCuts_;
105  }
106}
107
108/**
109  Take care to detach from the owning CbcNode and decrement the reference
110  count in the parent.  If this is the last nodeInfo object pointing to the
111  parent, make a recursive call to delete the parent.
112*/
113CbcNodeInfo::~CbcNodeInfo()
114{
115#ifdef CHECK_NODE
116  printf("CbcNodeInfo %x Destructor parent %x\n",this,parent_);
117#endif
118
119  assert(!numberPointingToThis_);
120  // But they may be some left (max nodes?)
121  for (int i=0;i<numberCuts_;i++) 
122    delete cuts_[i];
123
124  delete [] cuts_;
125  if (owner_) 
126    owner_->nullNodeInfo();
127  if (parent_) {
128    int numberLinks = parent_->decrement();
129    if (!numberLinks) delete parent_;
130  }
131}
132
133
134//#define ALLCUTS
135void
136CbcNodeInfo::decrementCuts(int change)
137{
138  int i;
139  // get rid of all remaining if negative
140  int changeThis;
141  if (change<0)
142    changeThis = numberBranchesLeft_;
143  else
144    changeThis = change;
145 // decrement cut counts
146  for (i=0;i<numberCuts_;i++) {
147    if (cuts_[i]) {
148      int number = cuts_[i]->decrement(changeThis);
149      if (!number) {
150        //printf("info %x del cut %d %x\n",this,i,cuts_[i]);
151        delete cuts_[i];
152        cuts_[i]=NULL;
153      }
154    }
155  }
156}
157void
158CbcNodeInfo::decrementParentCuts(int change)
159{
160  if (parent_) {
161    // get rid of all remaining if negative
162    int changeThis;
163    if (change<0)
164      changeThis = numberBranchesLeft_;
165    else
166      changeThis = change;
167    int i;
168    // Get over-estimate of space needed for basis
169    CoinWarmStartBasis dummy;
170    dummy.setSize(0,numberRows_+numberCuts_);
171    buildRowBasis(dummy);
172    /* everything is zero (i.e. free) so we can use to see
173       if latest basis */
174    CbcNodeInfo * thisInfo = parent_;
175    while (thisInfo) 
176      thisInfo = thisInfo->buildRowBasis(dummy);
177    // decrement cut counts
178    thisInfo = parent_;
179    int numberRows=numberRows_;
180    while (thisInfo) {
181      for (i=thisInfo->numberCuts_-1;i>=0;i--) {
182        CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
183#ifdef ALLCUTS
184        status = CoinWarmStartBasis::isFree;
185#endif
186        if (thisInfo->cuts_[i]) {
187          int number=1;
188          if (status!=CoinWarmStartBasis::basic) {
189            // tight - drop 1 or 2
190            if (change<0)
191              number = thisInfo->cuts_[i]->decrement(changeThis);
192            else
193              number = thisInfo->cuts_[i]->decrement(change);
194          }
195          if (!number) {
196            delete thisInfo->cuts_[i];
197            thisInfo->cuts_[i]=NULL;
198          }
199        }
200      }
201      thisInfo = thisInfo->parent_;
202    }
203  }
204}
205
206void
207CbcNodeInfo::incrementParentCuts(int change)
208{
209  if (parent_) {
210    int i;
211    // Get over-estimate of space needed for basis
212    CoinWarmStartBasis dummy;
213    dummy.setSize(0,numberRows_+numberCuts_);
214    /* everything is zero (i.e. free) so we can use to see
215       if latest basis */
216    buildRowBasis(dummy);
217    CbcNodeInfo * thisInfo = parent_;
218    while (thisInfo) 
219      thisInfo = thisInfo->buildRowBasis(dummy);
220    // increment cut counts
221    thisInfo = parent_;
222    int numberRows=numberRows_;
223    while (thisInfo) {
224      for (i=thisInfo->numberCuts_-1;i>=0;i--) {
225        CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
226#ifdef ALLCUTS
227        status = CoinWarmStartBasis::isFree;
228#endif
229        if (thisInfo->cuts_[i]&&status!=CoinWarmStartBasis::basic) {
230          thisInfo->cuts_[i]->increment(change);
231        }
232      }
233      thisInfo = thisInfo->parent_;
234    }
235  }
236}
237
238/*
239  Append cuts to the cuts_ array in a nodeInfo. The initial reference count
240  is set to numberToBranchOn, which will normally be the number of arms
241  defined for the CbcBranchingObject attached to the CbcNode that owns this
242  CbcNodeInfo.
243*/
244void
245CbcNodeInfo::addCuts (OsiCuts & cuts, int numberToBranchOn,
246                      int * whichGenerator)
247{
248  int numberCuts = cuts.sizeRowCuts();
249  if (numberCuts) {
250    int i;
251    if (!numberCuts_) {
252      cuts_ = new CbcCountRowCut * [numberCuts];
253    } else {
254      CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
255      memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *));
256      delete [] cuts_;
257      cuts_ = temp;
258    }
259    for (i=0;i<numberCuts;i++) {
260      CbcCountRowCut * thisCut = new CbcCountRowCut(*cuts.rowCutPtr(i),
261                                                    this,numberCuts_);
262      thisCut->increment(numberToBranchOn); 
263      cuts_[numberCuts_++] = thisCut;
264#ifdef CBC_DEBUG
265      int n=thisCut->row().getNumElements();
266#if CBC_DEBUG>1
267      printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
268             thisCut->ub());
269#endif
270      int j;
271#if CBC_DEBUG>1
272      const int * index = thisCut->row().getIndices();
273#endif
274      const double * element = thisCut->row().getElements();
275      for (j=0;j<n;j++) {
276#if CBC_DEBUG>1
277        printf(" (%d,%g)",index[j],element[j]);
278#endif
279        assert(fabs(element[j])>1.00e-12);
280      }
281      printf("\n");
282#endif
283    }
284  }
285}
286
287void
288CbcNodeInfo::addCuts(int numberCuts, CbcCountRowCut ** cut, 
289                     int numberToBranchOn)
290{
291  if (numberCuts) {
292    int i;
293    if (!numberCuts_) {
294      cuts_ = new CbcCountRowCut * [numberCuts];
295    } else {
296      CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
297      memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *));
298      delete [] cuts_;
299      cuts_ = temp;
300    }
301    for (i=0;i<numberCuts;i++) {
302      CbcCountRowCut * thisCut = cut[i];
303      thisCut->setInfo(this,numberCuts_);
304      //printf("info %x cut %d %x\n",this,i,thisCut);
305      thisCut->increment(numberToBranchOn); 
306      cuts_[numberCuts_++] = thisCut;
307#ifdef CBC_DEBUG
308      int n=thisCut->row().getNumElements();
309#if CBC_DEBUG>1
310      printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
311             thisCut->ub());
312#endif
313      int j;
314#if CBC_DEBUG>1
315      const int * index = thisCut->row().getIndices();
316#endif
317      const double * element = thisCut->row().getElements();
318      for (j=0;j<n;j++) {
319#if CBC_DEBUG>1
320        printf(" (%d,%g)",index[j],element[j]);
321#endif
322        assert(fabs(element[j])>1.00e-12);
323      }
324      printf("\n");
325#endif
326    }
327  }
328}
329
330// delete cuts
331void
332CbcNodeInfo::deleteCuts(int numberToDelete, CbcCountRowCut ** cuts)
333{
334  int i;
335  int j;
336  int last=-1;
337  for (i=0;i<numberToDelete;i++) {
338    CbcCountRowCut * next = cuts[i];
339    for (j=last+1;j<numberCuts_;j++) {
340      if (next==cuts_[j])
341        break;
342    }
343    if (j==numberCuts_) {
344      // start from beginning
345      for (j=0;j<last;j++) {
346        if (next==cuts_[j])
347          break;
348      }
349      assert(j<last);
350    }
351    last=j;
352    int number = cuts_[j]->decrement();
353    if (!number)
354      delete cuts_[j];
355    cuts_[j]=NULL;
356  }
357  j=0;
358  for (i=0;i<numberCuts_;i++) {
359    if (cuts_[i])
360      cuts_[j++]=cuts_[i];
361  }
362  numberCuts_ = j;
363}
364
365// delete cuts
366void
367CbcNodeInfo::deleteCuts(int numberToDelete, int * which)
368{
369  int i;
370  for (i=0;i<numberToDelete;i++) {
371    int iCut=which[i];
372    int number = cuts_[iCut]->decrement();
373    if (!number)
374      delete cuts_[iCut];
375    cuts_[iCut]=NULL;
376  }
377  int j=0;
378  for (i=0;i<numberCuts_;i++) {
379    if (cuts_[i])
380      cuts_[j++]=cuts_[i];
381  }
382  numberCuts_ = j;
383}
384
385// Really delete a cut
386void 
387CbcNodeInfo::deleteCut(int whichOne)
388{
389  assert(whichOne<numberCuts_);
390  cuts_[whichOne]=NULL;
391}
392
393CbcFullNodeInfo::CbcFullNodeInfo() :
394  CbcNodeInfo(),
395  basis_(),
396  numberIntegers_(0),
397  lower_(NULL),
398  upper_(NULL)
399{
400}
401CbcFullNodeInfo::CbcFullNodeInfo(CbcModel * model,
402                                 int numberRowsAtContinuous) :
403  CbcNodeInfo()
404{
405  OsiSolverInterface * solver = model->solver();
406  numberRows_ = numberRowsAtContinuous;
407  numberIntegers_ = model->numberIntegers();
408  int numberColumns = model->getNumCols();
409  lower_ = new double [numberColumns];
410  upper_ = new double [numberColumns];
411  const double * lower = solver->getColLower();
412  const double * upper = solver->getColUpper();
413  int i;
414
415  for (i=0;i<numberColumns;i++) {
416    lower_[i]=lower[i];
417    upper_[i]=upper[i];
418  }
419
420  basis_ =  dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
421}
422
423CbcFullNodeInfo::CbcFullNodeInfo(const CbcFullNodeInfo & rhs) :
424  CbcNodeInfo(rhs)
425{ 
426  basis_= dynamic_cast<CoinWarmStartBasis *>(rhs.basis_->clone()) ;
427  numberIntegers_=rhs.numberIntegers_;
428  lower_=NULL;
429  upper_=NULL;
430  if (rhs.lower_!=NULL) {
431    int numberColumns = basis_->getNumStructural();
432    lower_ = new double [numberColumns];
433    upper_ = new double [numberColumns];
434    assert (upper_!=NULL);
435    memcpy(lower_,rhs.lower_,numberColumns*sizeof(double));
436    memcpy(upper_,rhs.upper_,numberColumns*sizeof(double));
437  }
438}
439
440CbcNodeInfo * 
441CbcFullNodeInfo::clone() const
442{ 
443  return (new CbcFullNodeInfo(*this));
444}
445
446CbcFullNodeInfo::~CbcFullNodeInfo ()
447{
448  delete basis_ ;
449  delete [] lower_;
450  delete [] upper_;
451}
452
453/*
454   The basis supplied as a parameter is deleted and replaced with a new basis
455   appropriate for the node, and lower and upper bounds on variables are
456   reset according to the stored bounds arrays. Any cuts associated with this
457   node are added to the list in addCuts, but not actually added to the
458   constraint system in the model.
459
460   Why pass in a basis at all? The short answer is ``We need the parameter to
461   pass out a basis, so might as well use it to pass in the size.''
462   
463   A longer answer is that in practice we take a memory allocation hit up in
464   addCuts1 (the only place applyToModel is called) when we setSize() the
465   basis that's passed in. It's immediately tossed here in favour of a clone
466   of the basis attached to this nodeInfo. This can probably be fixed, given
467   a bit of thought.
468*/
469
470void CbcFullNodeInfo::applyToModel (CbcModel *model,
471                                    CoinWarmStartBasis *&basis,
472                                    CbcCountRowCut **addCuts,
473                                    int &currentNumberCuts) const 
474
475{ OsiSolverInterface *solver = model->solver() ;
476
477  // branch - do bounds
478  int i;
479  int numberColumns = model->getNumCols();
480  for (i=0;i<numberColumns;i++) {
481    solver->setColBounds(i,lower_[i],upper_[i]);
482  }
483  // move basis - but make sure size stays
484  int numberRows = basis->getNumArtificial();
485  delete basis ;
486  basis = dynamic_cast<CoinWarmStartBasis *>(basis_->clone()) ;
487  basis->resize(numberRows,numberColumns);
488  for (i=0;i<numberCuts_;i++) 
489    addCuts[currentNumberCuts+i]= cuts_[i];
490  currentNumberCuts += numberCuts_;
491  assert(!parent_);
492  return ;
493}
494
495/* Builds up row basis backwards (until original model).
496   Returns NULL or previous one to apply .
497   Depends on Free being 0 and impossible for cuts
498*/
499CbcNodeInfo * 
500CbcFullNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
501{
502  const unsigned int * saved = 
503    (const unsigned int *) basis_->getArtificialStatus();
504  unsigned int * now = 
505    (unsigned int *) basis.getArtificialStatus();
506  int number=basis_->getNumArtificial()>>4;;
507  int i;
508  for (i=0;i<number;i++) { 
509    if (!now[i])
510      now[i] = saved[i];
511  }
512  return NULL;
513}
514
515// Default constructor
516CbcPartialNodeInfo::CbcPartialNodeInfo()
517
518  : CbcNodeInfo(),
519    basisDiff_(NULL),
520    variables_(NULL),
521    newBounds_(NULL),
522    numberChangedBounds_(0)
523
524{ /* this space intentionally left blank */ }
525
526// Constructor from current state
527CbcPartialNodeInfo::CbcPartialNodeInfo (CbcNodeInfo *parent, CbcNode *owner,
528                                        int numberChangedBounds,
529                                        const int *variables,
530                                        const double *boundChanges,
531                                        const CoinWarmStartDiff *basisDiff)
532 : CbcNodeInfo(parent,owner)
533{
534  basisDiff_ = basisDiff->clone() ;
535
536  numberChangedBounds_ = numberChangedBounds;
537  variables_ = new int [numberChangedBounds_];
538  newBounds_ = new double [numberChangedBounds_];
539
540  int i ;
541  for (i=0;i<numberChangedBounds_;i++) {
542    variables_[i]=variables[i];
543    newBounds_[i]=boundChanges[i];
544  }
545}
546
547CbcPartialNodeInfo::CbcPartialNodeInfo (const CbcPartialNodeInfo & rhs)
548
549  : CbcNodeInfo(rhs.parent_)
550
551{ basisDiff_ = rhs.basisDiff_->clone() ;
552
553  numberChangedBounds_ = rhs.numberChangedBounds_;
554  variables_ = new int [numberChangedBounds_];
555  newBounds_ = new double [numberChangedBounds_];
556
557  int i ;
558  for (i=0;i<numberChangedBounds_;i++) {
559    variables_[i]=rhs.variables_[i];
560    newBounds_[i]=rhs.newBounds_[i];
561  }
562}
563
564CbcNodeInfo * 
565CbcPartialNodeInfo::clone() const
566{ 
567  return (new CbcPartialNodeInfo(*this));
568}
569
570
571CbcPartialNodeInfo::~CbcPartialNodeInfo ()
572{
573  delete basisDiff_ ;
574  delete [] variables_;
575  delete [] newBounds_;
576}
577
578
579/**
580   The basis supplied as a parameter is incrementally modified, and lower and
581   upper bounds on variables in the model are incrementally modified. Any
582   cuts associated with this node are added to the list in addCuts.
583*/
584
585void CbcPartialNodeInfo::applyToModel (CbcModel *model,
586                                       CoinWarmStartBasis *&basis,
587                                       CbcCountRowCut **addCuts,
588                                       int &currentNumberCuts) const 
589
590{ OsiSolverInterface *solver = model->solver();
591
592  basis->applyDiff(basisDiff_) ;
593
594  // branch - do bounds
595  int i;
596  for (i=0;i<numberChangedBounds_;i++) {
597    int variable = variables_[i];
598    if ((variable&0x80000000)==0) {
599      // lower bound changing
600      solver->setColLower(variable,newBounds_[i]);
601    } else {
602      // upper bound changing
603      solver->setColUpper(variable&0x7fffffff,newBounds_[i]);
604    }
605  }
606  for (i=0;i<numberCuts_;i++) 
607    addCuts[currentNumberCuts+i]= cuts_[i];
608  currentNumberCuts += numberCuts_;
609  return ;
610}
611
612/* Builds up row basis backwards (until original model).
613   Returns NULL or previous one to apply .
614   Depends on Free being 0 and impossible for cuts
615*/
616
617CbcNodeInfo * 
618CbcPartialNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
619
620{ basis.applyDiff(basisDiff_) ;
621
622  return parent_ ; }
623
624
625CbcNode::CbcNode() :
626  nodeInfo_(NULL),
627  objectiveValue_(1.0e100),
628  guessedObjectiveValue_(1.0e100),
629  branch_(NULL),
630  depth_(-1),
631  numberUnsatisfied_(0),
632  nodeNumber_(-1)
633{
634#ifdef CHECK_NODE
635  printf("CbcNode %x Constructor\n",this);
636#endif
637}
638
639CbcNode::CbcNode(CbcModel * model,
640                 CbcNode * lastNode) :
641  nodeInfo_(NULL),
642  objectiveValue_(1.0e100),
643  guessedObjectiveValue_(1.0e100),
644  branch_(NULL),
645  depth_(-1),
646  numberUnsatisfied_(0),
647  nodeNumber_(model->getNodeCount())
648{
649#ifdef CHECK_NODE
650  printf("CbcNode %x Constructor from model\n",this);
651#endif
652  OsiSolverInterface * solver = model->solver();
653  objectiveValue_ = solver->getObjSense()*solver->getObjValue();
654
655  if (lastNode)
656    lastNode->nodeInfo_->increment();
657}
658
659
660void
661CbcNode::createInfo (CbcModel *model,
662                     CbcNode *lastNode,
663                     const CoinWarmStartBasis *lastws,
664                     const double *lastLower, const double *lastUpper,
665                     int numberOldActiveCuts,int numberNewCuts)
666{ OsiSolverInterface * solver = model->solver();
667/*
668  The root --- no parent. Create full basis and bounds information.
669*/
670  if (!lastNode)
671  { nodeInfo_=new CbcFullNodeInfo(model,solver->getNumRows()); }
672/*
673  Not the root. Create an edit from the parent's basis & bound information.
674  This is not quite as straightforward as it seems. We need to reintroduce
675  cuts we may have dropped out of the basis, in the correct position, because
676  this whole process is strictly positional. Start by grabbing the current
677  basis.
678*/
679  else
680  { const CoinWarmStartBasis* ws =
681      dynamic_cast<const CoinWarmStartBasis*>(solver->getWarmStart());
682    assert(ws!=NULL); // make sure not volume
683    //int numberArtificials = lastws->getNumArtificial();
684    int numberColumns = solver->getNumCols();
685   
686    const double * lower = solver->getColLower();
687    const double * upper = solver->getColUpper();
688
689    int i;
690/*
691  Create a clone and resize it to hold all the structural constraints, plus
692  all the cuts: old cuts, both active and inactive (currentNumberCuts), and
693  new cuts (numberNewCuts).
694
695  TODO: You'd think that the set of constraints (logicals) in the expanded
696        basis should match the set represented in lastws. At least, that's
697        what I thought. But at the point I first looked hard at this bit of
698        code, it turned out that lastws was the stripped basis produced at
699        the end of addCuts(), rather than the raw basis handed back by
700        addCuts1(). The expanded basis here is equivalent to the raw basis of
701        addCuts1(). I said ``whoa, that's not good, I must have introduced a
702        bug'' and went back to John's code to see where I'd gone wrong.
703        And discovered the same `error' in his code.
704
705        After a bit of thought, my conclusion is that correctness is not
706        affected by whether lastws is the stripped or raw basis. The diffs
707        have no semantics --- just a set of changes that need to be made
708        to convert lastws into expanded. I think the only effect is that we
709        store a lot more diffs (everything in expanded that's not covered by
710        the stripped basis). But I need to give this more thought. There
711        may well be some subtle error cases.
712
713        In the mean time, I've twiddled addCuts() to set lastws to the raw
714        basis. Makes me (Lou) less nervous to compare apples to apples.
715*/
716    CoinWarmStartBasis *expanded = 
717      dynamic_cast<CoinWarmStartBasis *>(ws->clone()) ;
718    int numberRowsAtContinuous = model->numberRowsAtContinuous();
719    int iFull = numberRowsAtContinuous+model->currentNumberCuts()+
720      numberNewCuts;
721    //int numberArtificialsNow = iFull;
722    //int maxBasisLength = ((iFull+15)>>4)+((numberColumns+15)>>4);
723    //printf("l %d full %d\n",maxBasisLength,iFull);
724    expanded->resize(iFull,numberColumns);
725#ifdef FULL_DEBUG
726    printf("Before expansion: orig %d, old %d, new %d, current %d\n",
727           numberRowsAtContinuous,numberOldActiveCuts,numberNewCuts,
728           model->currentNumberCuts()) ;
729    ws->print();
730#endif
731/*
732  Now fill in the expanded basis. Any indices beyond nPartial must
733  be cuts created while processing this node --- they can be copied directly
734  into the expanded basis. From nPartial down, pull the status of active cuts
735  from ws, interleaving with a B entry for the deactivated (loose) cuts.
736*/
737    int numberDropped = model->currentNumberCuts()-numberOldActiveCuts;
738    int iCompact=iFull-numberDropped;
739    CbcCountRowCut ** cut = model->addedCuts();
740    int nPartial = model->currentNumberCuts()+numberRowsAtContinuous;
741    iFull--;
742    for (;iFull>=nPartial;iFull--) {
743      CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
744      //assert (status != CoinWarmStartBasis::basic); // may be permanent cut
745      expanded->setArtifStatus(iFull,status);
746    }
747    for (;iFull>=numberRowsAtContinuous;iFull--) {
748      if (cut[iFull-numberRowsAtContinuous]) {
749        CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
750        // If no cut generator being used then we may have basic variables
751        //if (model->getMaximumCutPasses()&&
752        //  status == CoinWarmStartBasis::basic)
753        //printf("cut basic\n");
754        expanded->setArtifStatus(iFull,status);
755      } else {
756        expanded->setArtifStatus(iFull,CoinWarmStartBasis::basic);
757      }
758    }
759#ifdef FULL_DEBUG
760    printf("Expanded basis\n");
761    expanded->print() ;
762    printf("Diffing against\n") ;
763    lastws->print() ;
764#endif   
765/*
766  Now that we have two bases in proper positional correspondence, creating
767  the actual diff is dead easy.
768*/
769
770    CoinWarmStartDiff *basisDiff = expanded->generateDiff(lastws) ;
771/*
772  Diff the bound vectors. It's assumed the number of structural variables is
773  not changing. Assuming that branching objects all involve integer variables,
774  we should see at least one bound change as a consequence of processing this
775  subproblem. Different types of branching objects could break this assertion.
776  Yes - cuts will break
777*/
778    double *boundChanges = new double [2*numberColumns] ;
779    int *variables = new int [2*numberColumns] ;
780    int numberChangedBounds=0;
781    for (i=0;i<numberColumns;i++) {
782      if (lower[i]!=lastLower[i]) {
783        variables[numberChangedBounds]=i;
784        boundChanges[numberChangedBounds++]=lower[i];
785      }
786      if (upper[i]!=lastUpper[i]) {
787        variables[numberChangedBounds]=i|0x80000000;
788        boundChanges[numberChangedBounds++]=upper[i];
789      }
790#ifdef CBC_DEBUG
791      if (lower[i]!=lastLower[i])
792        printf("lower on %d changed from %g to %g\n",
793               i,lastLower[i],lower[i]);
794      if (upper[i]!=lastUpper[i])
795        printf("upper on %d changed from %g to %g\n",
796               i,lastUpper[i],upper[i]);
797#endif
798    }
799#ifdef CBC_DEBUG
800    printf("%d changed bounds\n",numberChangedBounds) ;
801#endif
802    if (lastNode->branchingObject()->boundBranch())
803      assert (numberChangedBounds);
804/*
805  Hand the lot over to the CbcPartialNodeInfo constructor, then clean up and
806  return.
807*/
808    nodeInfo_ =
809      new CbcPartialNodeInfo(lastNode->nodeInfo_,this,numberChangedBounds,
810                             variables,boundChanges,basisDiff) ;
811    delete basisDiff ;
812    delete [] boundChanges;
813    delete [] variables;
814    delete expanded ;
815    delete ws;
816  }
817}
818
819/*
820  The routine scans through the object list of the model looking for objects
821  that indicate infeasibility. It tests each object using strong branching
822  and selects the one with the least objective degradation.  A corresponding
823  branching object is left attached to lastNode.
824
825  If strong branching is disabled, a candidate object is chosen essentially
826  at random (whatever object ends up in pos'n 0 of the candidate array).
827
828  If a branching candidate is found to be monotone, bounds are set to fix the
829  variable and the routine immediately returns (the caller is expected to
830  reoptimize).
831
832  If a branching candidate is found to result in infeasibility in both
833  directions, the routine immediately returns an indication of infeasibility.
834
835  Returns:  0   both branch directions are feasible
836           -1   branching variable is monotone
837           -2   infeasible
838
839  Original comments:
840    Here could go cuts etc etc
841    For now just fix on objective from strong branching.
842*/
843
844int CbcNode::chooseBranch (CbcModel *model, CbcNode *lastNode)
845
846{ if (lastNode)
847    depth_ = lastNode->depth_+1;
848  else
849    depth_ = 0;
850  delete branch_;
851  branch_=NULL;
852  OsiSolverInterface * solver = model->solver();
853  double saveObjectiveValue = solver->getObjValue();
854  double objectiveValue = solver->getObjSense()*saveObjectiveValue;
855  const double * lower = solver->getColLower();
856  const double * upper = solver->getColUpper();
857
858  int anyAction=0;
859  double integerTolerance = 
860    model->getDblParam(CbcModel::CbcIntegerTolerance);
861  int i;
862  bool beforeSolution = model->getSolutionCount()==0;
863  int numberStrong=model->numberStrong();
864  int numberObjects = model->numberObjects();
865  int maximumStrong = CoinMax(CoinMin(model->numberStrong(),numberObjects),1);
866  int numberColumns = model->getNumCols();
867  double * saveUpper = new double[numberColumns];
868  double * saveLower = new double[numberColumns];
869
870  // Save solution in case heuristics need good solution later
871
872  double * saveSolution = new double[numberColumns];
873  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
874
875/*
876  Get a branching decision object. Use the default decision criteria unless
877  the user has loaded a decision method into the model.
878*/
879  CbcBranchDecision *decision = model->branchingMethod();
880  if (!decision)
881    decision = new CbcBranchDefaultDecision();
882
883  typedef struct {
884    CbcBranchingObject * possibleBranch; // what a branch would do
885    double upMovement; // cost going up (and initial away from feasible)
886    double downMovement; // cost going down
887    int numIntInfeasUp ; // without odd ones
888    int numObjInfeasUp ; // just odd ones
889    bool finishedUp; // true if solver finished
890    int numItersUp ; // number of iterations in solver
891    int numIntInfeasDown ; // without odd ones
892    int numObjInfeasDown ; // just odd ones
893    bool finishedDown; // true if solver finished
894    int numItersDown; // number of iterations in solver
895    int objectNumber; // Which object it is
896  } Strong;
897  Strong * choice = new Strong[maximumStrong];
898  for (i=0;i<numberColumns;i++) {
899    saveLower[i] = lower[i];
900    saveUpper[i] = upper[i];
901  }
902  // Some objects may compute an estimate of best solution from here
903  double estimatedDegradation=0.0; 
904  int numberIntegerInfeasibilities=0; // without odd ones
905
906  // We may go round this loop twice (only if we think we have solution)
907  for (int iPass=0;iPass<2;iPass++) {
908
909    // compute current state
910    int numberObjectInfeasibilities; // just odd ones
911    model->feasibleSolution(
912                            numberIntegerInfeasibilities,
913                            numberObjectInfeasibilities);
914    // If forcePriority > 0 then we want best solution
915    const double * bestSolution = NULL;
916    int hotstartStrategy=model->getHotstartStrategy();
917    if (hotstartStrategy>0) {
918      bestSolution = model->bestSolution();
919    }
920   
921    // Some objects may compute an estimate of best solution from here
922    estimatedDegradation=0.0; 
923    numberUnsatisfied_ = 0;
924    int bestPriority=INT_MAX;
925    /*
926      Scan for branching objects that indicate infeasibility. Choose the best
927      maximumStrong candidates, using priority as the first criteria, then
928      integer infeasibility.
929     
930      The algorithm is to fill the choice array with a set of good candidates (by
931      infeasibility) with priority bestPriority.  Finding a candidate with
932      priority better (less) than bestPriority flushes the choice array. (This
933      serves as initialization when the first candidate is found.)
934     
935      A new candidate is added to choices only if its infeasibility exceeds the
936      current max infeasibility (mostAway). When a candidate is added, it
937      replaces the candidate with the smallest infeasibility (tracked by
938      iSmallest).
939    */
940    int iSmallest = 0;
941    double mostAway = integerTolerance;
942    for (i = 0 ; i < maximumStrong ; i++) choice[i].possibleBranch = NULL ;
943    numberStrong=0;
944    for (i=0;i<numberObjects;i++) {
945      const CbcObject * object = model->object(i);
946      int preferredWay;
947      double infeasibility = object->infeasibility(preferredWay);
948      int priorityLevel = model->priority(i);
949      if (bestSolution) {
950        // we are doing hot start
951        const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
952        if (thisOne) {
953          int iColumn = thisOne->modelSequence();
954          if (saveUpper[iColumn]>saveLower[iColumn]) {
955            double value = saveSolution[iColumn];
956            double targetValue = bestSolution[iColumn];
957            //double originalLower = thisOne->originalLower();
958            //double originalUpper = thisOne->originalUpper();
959            // switch off if not possible
960            if (targetValue>=saveLower[iColumn]&&targetValue<=saveUpper[iColumn]) {
961              /* priority outranks rest always if hotstartStrategy >1
962                 otherwise can be downgraded if at correct level.
963                 Infeasibility may be increased by targetValue to choose 1.0 values first.
964              */
965              if (fabs(value-targetValue)>integerTolerance) {
966                if (value>targetValue) {
967                  infeasibility += value;
968                  preferredWay=-1;
969                } else {
970                  infeasibility += targetValue;
971                  preferredWay=1;
972                }
973              } else if (hotstartStrategy>1) {
974                if (targetValue==saveLower[iColumn]) {
975                  infeasibility += integerTolerance+1.0e-12;
976                  preferredWay=-1;
977                } else if (targetValue==saveUpper[iColumn]) {
978                  infeasibility += integerTolerance+1.0e-12;
979                  preferredWay=1;
980                } else {
981                  infeasibility += integerTolerance+1.0e-12;
982                  preferredWay=1;
983                }
984              } else {
985                priorityLevel += 10000000;
986              }
987            } else {
988              // switch off if not possible
989              bestSolution=NULL;
990              model->setHotstartStrategy(0);
991            }
992          }
993        }
994      }
995      if (infeasibility>integerTolerance) {
996        // Increase estimated degradation to solution
997        estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate());
998        numberUnsatisfied_++;
999        // Better priority? Flush choices.
1000        if (priorityLevel<bestPriority) {
1001          int j;
1002          iSmallest=0;
1003          for (j=0;j<maximumStrong;j++) {
1004            choice[j].upMovement=0.0;
1005            delete choice[j].possibleBranch;
1006            choice[j].possibleBranch=NULL;
1007          }
1008          bestPriority = priorityLevel;
1009          mostAway=integerTolerance;
1010          numberStrong=0;
1011        } else if (priorityLevel>bestPriority) {
1012          continue;
1013        }
1014        // Check for suitability based on infeasibility.
1015        if (infeasibility>mostAway) {
1016          //add to list
1017          choice[iSmallest].upMovement=infeasibility;
1018          delete choice[iSmallest].possibleBranch;
1019          choice[iSmallest].possibleBranch=object->createBranch(preferredWay);
1020          numberStrong = CoinMax(numberStrong,iSmallest+1);
1021          // Save which object it was
1022          choice[iSmallest].objectNumber=i;
1023          int j;
1024          iSmallest=-1;
1025          mostAway = 1.0e50;
1026          for (j=0;j<maximumStrong;j++) {
1027            if (choice[j].upMovement<mostAway) {
1028              mostAway=choice[j].upMovement;
1029              iSmallest=j;
1030            }
1031          }
1032        }
1033      }
1034    }
1035    if (numberUnsatisfied_) {
1036      // some infeasibilities - go to next steps
1037      break;
1038    } else if (!iPass) {
1039      // looks like a solution - get paranoid
1040      bool roundAgain=false;
1041      // get basis
1042      CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
1043      if (!ws)
1044        break;
1045      for (i=0;i<numberColumns;i++) {
1046        double value = saveSolution[i];
1047        if (value<lower[i]) {
1048          saveSolution[i]=lower[i];
1049          roundAgain=true;
1050          ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
1051        } else if (value>upper[i]) {
1052          saveSolution[i]=upper[i];
1053          roundAgain=true;
1054          ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
1055        } 
1056      }
1057      if (roundAgain) {
1058        // restore basis
1059        solver->setWarmStart(ws);
1060        delete ws;
1061        solver->resolve();
1062        memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1063        if (!solver->isProvenOptimal()) {
1064          // infeasible
1065          anyAction=-2;
1066          break;
1067        }
1068      } else {
1069        delete ws;
1070        break;
1071      }
1072    }
1073  }
1074  /* Some solvers can do the strong branching calculations faster if
1075     they do them all at once.  At present only Clp does for ordinary
1076     integers but I think this coding would be easy to modify
1077  */
1078  bool allNormal=true; // to say if we can do fast strong branching
1079  // Say which one will be best
1080  int bestChoice=0;
1081  double worstInfeasibility=0.0;
1082  for (i=0;i<numberStrong;i++) {
1083    choice[i].numIntInfeasUp = numberUnsatisfied_;
1084    choice[i].numIntInfeasDown = numberUnsatisfied_;
1085    if (!dynamic_cast <const CbcSimpleInteger *> (model->object(choice[i].objectNumber)))
1086      allNormal=false; // Something odd so lets skip clever fast branching
1087    if ( !model->object(choice[i].objectNumber)->boundBranch())
1088      numberStrong=0; // switch off
1089    // Do best choice in case switched off
1090    if (choice[i].upMovement>worstInfeasibility) {
1091      worstInfeasibility=choice[i].upMovement;
1092      bestChoice=i;
1093    }
1094  }
1095/*
1096  Is strong branching enabled? If so, set up and do it. Otherwise, we'll
1097  fall through to simple branching.
1098
1099  Setup for strong branching involves saving the current basis (for restoration
1100  afterwards) and setting up for hot starts.
1101*/
1102  if (numberStrong&&model->numberStrong()) {
1103   
1104    bool solveAll=false; // set true to say look at all even if some fixed (experiment)
1105    // Save basis
1106    CoinWarmStart * ws = solver->getWarmStart();
1107    // save limit
1108    int saveLimit;
1109    solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
1110    if (beforeSolution)
1111      solver->setIntParam(OsiMaxNumIterationHotStart,10000); // go to end
1112
1113    /* If we are doing all strong branching in one go then we create new arrays
1114       to store information.  If clp NULL then doing old way.
1115       Going down -
1116       outputSolution[2*i] is final solution.
1117       outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown
1118       outputStuff[2*i+numberStrong] is number iterations
1119       On entry newUpper[i] is new upper bound, on exit obj change
1120       Going up -
1121       outputSolution[2*i+1] is final solution.
1122       outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
1123       outputStuff[2*i+1+numberStrong] is number iterations
1124       On entry newLower[i] is new lower bound, on exit obj change
1125    */
1126    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
1127    ClpSimplex * clp=NULL;
1128    double * newLower = NULL;
1129    double * newUpper = NULL;
1130    double ** outputSolution=NULL;
1131    int * outputStuff=NULL;
1132    int saveLogLevel=0;
1133    // For moment - until full testing - go back to normal way
1134    //allNormal=false;
1135    if (osiclp&& allNormal) {
1136      clp = osiclp->getModelPtr();
1137      saveLogLevel = clp->logLevel();
1138      clp->setLogLevel(0);
1139      int saveOptions = clp->specialOptions();
1140      int startFinishOptions;
1141      int specialOptions = osiclp->specialOptions();
1142      if((specialOptions&1)==0) {
1143        startFinishOptions=0;
1144        clp->setSpecialOptions(saveOptions|(64|1024));
1145      } else {
1146        startFinishOptions=1+2+4;
1147        //startFinishOptions=1+4; // for moment re-factorize
1148        if((specialOptions&4)==0) 
1149          clp->setSpecialOptions(saveOptions|(64|128|512|1024|4096));
1150        else
1151          clp->setSpecialOptions(saveOptions|(64|128|512|1024|2048|4096));
1152      }
1153      // User may want to clean up before strong branching
1154      if ((clp->specialOptions()&32)!=0) {
1155        clp->primal(1);
1156        if (clp->numberIterations())
1157          model->messageHandler()->message(CBC_ITERATE_STRONG,model->messages())
1158            << clp->numberIterations()
1159            <<CoinMessageEol;
1160      }
1161      int saveMaxIts = clp->maximumIterations();
1162      clp->setMaximumIterations(saveLimit);
1163      // Clp - do a different way
1164      newLower = new double[numberStrong];
1165      newUpper = new double[numberStrong];
1166      outputSolution = new double * [2*numberStrong];
1167      outputStuff = new int [4*numberStrong];
1168      int * which = new int[numberStrong];
1169      for (i=0;i<numberStrong;i++) {
1170        int iObject = choice[i].objectNumber;
1171        const CbcObject * object = model->object(iObject);
1172        const CbcSimpleInteger * simple = dynamic_cast <const CbcSimpleInteger *> (object);
1173        int iSequence = simple->modelSequence();
1174        newLower[i]= ceil(saveSolution[iSequence]);
1175        newUpper[i]= floor(saveSolution[iSequence]);
1176        which[i]=iSequence;
1177        outputSolution[2*i]= new double [numberColumns];
1178        outputSolution[2*i+1]= new double [numberColumns];
1179      }
1180      int returnCode=clp->strongBranching(numberStrong,which,
1181                                          newLower, newUpper,outputSolution,
1182                                          outputStuff,outputStuff+2*numberStrong,!solveAll,false,
1183                                          startFinishOptions);
1184      clp->setSpecialOptions(saveOptions); // restore
1185      clp->setMaximumIterations(saveMaxIts);
1186      if (returnCode==-2) {
1187        // bad factorization!!!
1188        // Doing normal way
1189        // Mark hot start
1190        solver->markHotStart();
1191        clp = NULL;
1192      }
1193      delete [] which;
1194    } else {
1195      // Doing normal way
1196      // Mark hot start
1197      solver->markHotStart();
1198    }
1199/*
1200  Open a loop to do the strong branching LPs. For each candidate variable,
1201  solve an LP with the variable forced down, then up. If a direction turns
1202  out to be infeasible or monotonic (i.e., over the dual objective cutoff),
1203  force the objective change to be big (1.0e100). If we determine the problem
1204  is infeasible, or find a monotone variable, escape the loop.
1205
1206  TODO: The `restore bounds' part might be better encapsulated as an
1207        unbranch() method. Branching objects more exotic than simple integers
1208        or cliques might not restrict themselves to variable bounds.
1209
1210  TODO: Virtuous solvers invalidate the current solution (or give bogus
1211        results :-) when the bounds are changed out from under them. So we
1212        need to do all the work associated with finding a new solution before
1213        restoring the bounds.
1214*/
1215    // If we have hit max time ignore fixing on one way
1216    bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
1217                        model->getDblParam(CbcModel::CbcMaximumSeconds));
1218    for (i = 0 ; i < numberStrong ; i++)
1219    { double objectiveChange ;
1220      double newObjectiveValue=1.0e100;
1221      // status is 0 finished, 1 infeasible and other
1222      int iStatus;
1223/*
1224  Try the down direction first. (Specify the initial branching alternative as
1225  down with a call to way(-1). Each subsequent call to branch() performs the
1226  specified branch and advances the branch object state to the next branch
1227  alternative.)
1228*/
1229      if (!clp) {
1230        choice[i].possibleBranch->way(-1) ;
1231        choice[i].possibleBranch->branch() ;
1232        solver->solveFromHotStart() ;
1233/*
1234  We now have an estimate of objective degradation that we can use for strong
1235  branching. If we're over the cutoff, the variable is monotone up.
1236  If we actually made it to optimality, check for a solution, and if we have
1237  a good one, call setBestSolution to process it. Note that this may reduce the
1238  cutoff, so we check again to see if we can declare this variable monotone.
1239*/
1240        if (solver->isProvenOptimal())
1241          iStatus=0; // optimal
1242        else if (solver->isIterationLimitReached()
1243                 &&!solver->isDualObjectiveLimitReached())
1244          iStatus=2; // unknown
1245        else
1246          iStatus=1; // infeasible
1247        newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1248        choice[i].numItersDown = solver->getIterationCount();
1249        objectiveChange = newObjectiveValue-objectiveValue ;
1250      } else {
1251        iStatus = outputStuff[2*i];
1252        choice[i].numItersDown = outputStuff[2*numberStrong+2*i];
1253        newObjectiveValue = objectiveValue+newUpper[i];
1254        solver->setColSolution(outputSolution[2*i]);
1255      }
1256      objectiveChange = newObjectiveValue  - objectiveValue;
1257      if (!iStatus) {
1258        choice[i].finishedDown = true ;
1259        if (newObjectiveValue>=model->getCutoff()) {
1260          objectiveChange = 1.0e100; // say infeasible
1261        } else {
1262          // See if integer solution
1263          if (model->feasibleSolution(choice[i].numIntInfeasDown,
1264                                      choice[i].numObjInfeasDown))
1265            { model->setBestSolution(CBC_STRONGSOL,
1266                                     newObjectiveValue,
1267                                     solver->getColSolution()) ;
1268            if (newObjectiveValue >= model->getCutoff())        //  *new* cutoff
1269              objectiveChange = 1.0e100 ;
1270            }
1271        }
1272      } else if (iStatus==1) {
1273        objectiveChange = 1.0e100 ;
1274      } else {
1275        // Can't say much as we did not finish
1276        choice[i].finishedDown = false ;
1277      }
1278      choice[i].downMovement = objectiveChange ;
1279
1280      // restore bounds
1281      if (!clp)
1282      { for (int j=0;j<numberColumns;j++) {
1283        if (saveLower[j] != lower[j])
1284          solver->setColLower(j,saveLower[j]);
1285        if (saveUpper[j] != upper[j])
1286          solver->setColUpper(j,saveUpper[j]);
1287        }
1288      }
1289      //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1290      //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersDown,
1291      //     choice[i].downMovement,choice[i].finishedDown,choice[i].numIntInfeasDown,
1292      //     choice[i].numObjInfeasDown);
1293
1294      // repeat the whole exercise, forcing the variable up
1295      if (!clp) {
1296        choice[i].possibleBranch->branch();
1297        solver->solveFromHotStart();
1298        if (solver->isProvenOptimal())
1299          iStatus=0; // optimal
1300        else if (solver->isIterationLimitReached()
1301                 &&!solver->isDualObjectiveLimitReached())
1302          iStatus=2; // unknown
1303        else
1304          iStatus=1; // infeasible
1305        newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1306        choice[i].numItersUp = solver->getIterationCount();
1307        objectiveChange = newObjectiveValue-objectiveValue ;
1308      } else {
1309        iStatus = outputStuff[2*i+1];
1310        choice[i].numItersUp = outputStuff[2*numberStrong+2*i+1];
1311        newObjectiveValue = objectiveValue+newLower[i];
1312        solver->setColSolution(outputSolution[2*i+1]);
1313      }
1314      objectiveChange = newObjectiveValue  - objectiveValue;
1315      if (!iStatus) {
1316        choice[i].finishedUp = true ;
1317        if (newObjectiveValue>=model->getCutoff()) {
1318          objectiveChange = 1.0e100; // say infeasible
1319        } else {
1320          // See if integer solution
1321          if (model->feasibleSolution(choice[i].numIntInfeasUp,
1322                                      choice[i].numObjInfeasUp))
1323            { model->setBestSolution(CBC_STRONGSOL,
1324                                     newObjectiveValue,
1325                                     solver->getColSolution()) ;
1326            if (newObjectiveValue >= model->getCutoff())        //  *new* cutoff
1327              objectiveChange = 1.0e100 ;
1328            }
1329        }
1330      } else if (iStatus==1) {
1331        objectiveChange = 1.0e100 ;
1332      } else {
1333        // Can't say much as we did not finish
1334        choice[i].finishedUp = false ;
1335      }
1336      choice[i].upMovement = objectiveChange ;
1337
1338      // restore bounds
1339      if (!clp)
1340      { for (int j=0;j<numberColumns;j++) {
1341        if (saveLower[j] != lower[j])
1342          solver->setColLower(j,saveLower[j]);
1343        if (saveUpper[j] != upper[j])
1344          solver->setColUpper(j,saveUpper[j]);
1345        }
1346      }
1347
1348      //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1349      //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersUp,
1350      //     choice[i].upMovement,choice[i].finishedUp,choice[i].numIntInfeasUp,
1351      //     choice[i].numObjInfeasUp);
1352
1353/*
1354  End of evaluation for this candidate variable. Possibilities are:
1355    * Both sides below cutoff; this variable is a candidate for branching.
1356    * Both sides infeasible or above the objective cutoff: no further action
1357      here. Break from the evaluation loop and assume the node will be purged
1358      by the caller.
1359    * One side below cutoff: Install the branch (i.e., fix the variable). Break
1360      from the evaluation loop and assume the node will be reoptimised by the
1361      caller.
1362*/
1363      if (choice[i].upMovement<1.0e100) {
1364        if(choice[i].downMovement<1.0e100) {
1365          // feasible - no action
1366        } else if (!hitMaxTime) {
1367          // up feasible, down infeasible
1368          anyAction=-1;
1369          if (!solveAll) {
1370            choice[i].possibleBranch->way(1);
1371            choice[i].possibleBranch->branch();
1372            break;
1373          }
1374        }
1375      } else {
1376        if(choice[i].downMovement<1.0e100) {
1377          if (!hitMaxTime) {
1378            // down feasible, up infeasible
1379            anyAction=-1;
1380            if (!solveAll) {
1381              choice[i].possibleBranch->way(-1);
1382              choice[i].possibleBranch->branch();
1383              break;
1384            }
1385          }
1386        } else {
1387          // neither side feasible
1388          anyAction=-2;
1389          break;
1390        }
1391      }
1392    }
1393
1394    if (!clp) {
1395      // Delete the snapshot
1396      solver->unmarkHotStart();
1397    } else {
1398      clp->setLogLevel(saveLogLevel);
1399      delete [] newLower;
1400      delete [] newUpper;
1401      delete [] outputStuff;
1402      int i;
1403      for (i=0;i<2*numberStrong;i++)
1404        delete [] outputSolution[i];
1405      delete [] outputSolution;
1406    }
1407    solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
1408    // restore basis
1409    solver->setWarmStart(ws);
1410    delete ws;
1411
1412/*
1413  anyAction >= 0 indicates that strong branching didn't produce any monotone
1414  variables. Sift through the candidates for the best one.
1415
1416  QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is
1417         local to this code block. Perhaps should be numberNodes_ from model?
1418         Unclear what this calculation is doing.
1419*/
1420    if (anyAction>=0) {
1421
1422      int numberNodes = model->getNodeCount();
1423      // get average cost per iteration and assume stopped ones
1424      // would stop after 50% more iterations at average cost??? !!! ???
1425      double averageCostPerIteration=0.0;
1426      double totalNumberIterations=1.0;
1427      int smallestNumberInfeasibilities=INT_MAX;
1428      for (i=0;i<numberStrong;i++) {
1429        totalNumberIterations += choice[i].numItersDown +
1430          choice[i].numItersUp ;
1431        averageCostPerIteration += choice[i].downMovement +
1432          choice[i].upMovement;
1433        smallestNumberInfeasibilities= 
1434          CoinMin(CoinMin(choice[i].numIntInfeasDown ,
1435                  choice[i].numIntInfeasUp ),
1436              smallestNumberInfeasibilities);
1437      }
1438      if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities)
1439        numberNodes=1000000; // switch off search for better solution
1440      numberNodes=1000000; // switch off anyway
1441      averageCostPerIteration /= totalNumberIterations;
1442      // all feasible - choose best bet
1443
1444#if 0
1445      for (i = 0 ; i < numberStrong ; i++)
1446      { int iColumn =
1447          model->integerVariable()[choice[i].possibleBranch->variable()] ;
1448         model->messageHandler()->message(CBC_STRONG,model->messages())
1449          << i << iColumn
1450          <<choice[i].downMovement<<choice[i].numIntInfeasDown
1451          <<choice[i].upMovement<<choice[i].numIntInfeasUp
1452          <<choice[i].possibleBranch->value()
1453          <<CoinMessageEol;
1454        int betterWay = decision->betterBranch(choice[i].possibleBranch,
1455                                              branch_,
1456                                              choice[i].upMovement,
1457                                              choice[i].numIntInfeasUp ,
1458                                              choice[i].downMovement,
1459                                              choice[i].numIntInfeasDown );
1460        if (betterWay) {
1461          delete branch_;
1462          // C) create branching object
1463          branch_ = choice[i].possibleBranch;
1464          choice[i].possibleBranch=NULL;
1465          branch_->way(betterWay);
1466        }
1467      }
1468#else
1469      // New method does all at once so it can be more sophisticated
1470      // in deciding how to balance actions.
1471      // But it does need arrays
1472      double * changeUp = new double [numberStrong];
1473      int * numberInfeasibilitiesUp = new int [numberStrong];
1474      double * changeDown = new double [numberStrong];
1475      int * numberInfeasibilitiesDown = new int [numberStrong];
1476      CbcBranchingObject ** objects = new CbcBranchingObject * [ numberStrong];
1477      for (i = 0 ; i < numberStrong ; i++) {
1478        int iColumn = choice[i].possibleBranch->variable() ;
1479        model->messageHandler()->message(CBC_STRONG,model->messages())
1480          << i << iColumn
1481          <<choice[i].downMovement<<choice[i].numIntInfeasDown
1482          <<choice[i].upMovement<<choice[i].numIntInfeasUp
1483          <<choice[i].possibleBranch->value()
1484          <<CoinMessageEol;
1485        changeUp[i]=choice[i].upMovement;
1486        numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
1487        changeDown[i]=choice[i].downMovement;
1488        numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
1489        objects[i] = choice[i].possibleBranch;
1490      }
1491      int whichObject = decision->bestBranch(objects,numberStrong,numberUnsatisfied_,
1492                                             changeUp,numberInfeasibilitiesUp,
1493                                             changeDown,numberInfeasibilitiesDown,
1494                                             objectiveValue);
1495      // move branching object and make sure it will not be deleted
1496      if (whichObject>=0) {
1497        branch_ = objects[whichObject];
1498        choice[whichObject].possibleBranch=NULL;
1499      }
1500      delete [] changeUp;
1501      delete [] numberInfeasibilitiesUp;
1502      delete [] changeDown;
1503      delete [] numberInfeasibilitiesDown;
1504      delete [] objects;
1505#endif
1506    }
1507  }
1508/*
1509  Simple branching. Probably just one, but we may have got here
1510  because of an odd branch e.g. a cut
1511*/
1512  else {
1513    // not strong
1514    // C) create branching object
1515    branch_ = choice[bestChoice].possibleBranch;
1516    choice[bestChoice].possibleBranch=NULL;
1517  }
1518  // Set guessed solution value
1519  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
1520/*
1521  Cleanup, then we're outta here.
1522*/
1523  if (!model->branchingMethod())
1524    delete decision;
1525   
1526  for (i=0;i<maximumStrong;i++)
1527    delete choice[i].possibleBranch;
1528  delete [] choice;
1529  delete [] saveLower;
1530  delete [] saveUpper;
1531 
1532  // restore solution
1533  solver->setColSolution(saveSolution);
1534  delete [] saveSolution;
1535  return anyAction;
1536}
1537
1538
1539CbcNode::CbcNode(const CbcNode & rhs) 
1540{ 
1541#ifdef CHECK_NODE
1542  printf("CbcNode %x Constructor from rhs %x\n",this,&rhs);
1543#endif
1544  if (rhs.nodeInfo_)
1545    nodeInfo_ = rhs.nodeInfo_->clone();
1546  else
1547    nodeInfo_=NULL;
1548  objectiveValue_=rhs.objectiveValue_;
1549  guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
1550  if (rhs.branch_)
1551    branch_=rhs.branch_->clone();
1552  else
1553    branch_=NULL;
1554  depth_ = rhs.depth_;
1555  numberUnsatisfied_ = rhs.numberUnsatisfied_;
1556  nodeNumber_=rhs.nodeNumber_;
1557}
1558
1559CbcNode &
1560CbcNode::operator=(const CbcNode & rhs)
1561{
1562  if (this != &rhs) {
1563    delete nodeInfo_;
1564    if (nodeInfo_)
1565      nodeInfo_ = rhs.nodeInfo_->clone();
1566    else
1567      nodeInfo_ = NULL;
1568    objectiveValue_=rhs.objectiveValue_;
1569    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
1570    if (rhs.branch_)
1571      branch_=rhs.branch_->clone();
1572    else
1573      branch_=NULL,
1574    depth_ = rhs.depth_;
1575    numberUnsatisfied_ = rhs.numberUnsatisfied_;
1576    nodeNumber_=rhs.nodeNumber_;
1577  }
1578  return *this;
1579}
1580
1581
1582CbcNode::~CbcNode ()
1583{
1584#ifdef CHECK_NODE
1585  if (nodeInfo_) 
1586    printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
1587         this,nodeInfo_,nodeInfo_->numberPointingToThis());
1588  else
1589    printf("CbcNode %x Destructor nodeInfo %x (?)\n",
1590         this,nodeInfo_);
1591#endif
1592  if (nodeInfo_) {
1593    nodeInfo_->nullOwner();
1594    int numberToDelete=nodeInfo_->numberBranchesLeft();
1595    //    CbcNodeInfo * parent = nodeInfo_->parent();
1596    if (nodeInfo_->decrement(numberToDelete)==0) {
1597      delete nodeInfo_;
1598    } else {
1599      //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,parent);
1600      // anyway decrement parent
1601      //if (parent)
1602      ///parent->decrement(1);
1603    }
1604  }
1605  delete branch_;
1606}
1607// Decrement  active cut counts
1608void 
1609CbcNode::decrementCuts(int change)
1610{
1611  if(nodeInfo_) {
1612    nodeInfo_->decrementCuts(change);
1613  }
1614}
1615void 
1616CbcNode::decrementParentCuts(int change)
1617{
1618  if(nodeInfo_) {
1619    nodeInfo_->decrementParentCuts(change);
1620  }
1621}
1622
1623/*
1624  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
1625  in the attached nodeInfo_.
1626*/
1627void
1628CbcNode::initializeInfo()
1629{
1630  assert(nodeInfo_ && branch_) ;
1631  nodeInfo_->initializeInfo(branch_->numberBranches());
1632}
1633// Nulls out node info
1634void 
1635CbcNode::nullNodeInfo()
1636{
1637  nodeInfo_=NULL;
1638}
1639
1640int
1641CbcNode::branch()
1642{
1643  double changeInGuessed=branch_->branch(true);
1644  guessedObjectiveValue_+= changeInGuessed;
1645  return nodeInfo_->branchedOn();
1646}
Note: See TracBrowser for help on using the repository browser.