source: trunk/CbcNode.cpp @ 118

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

always back to correct level

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 52.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,int numberPassesLeft)
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  int anyAction=0;
858  double integerTolerance = 
859    model->getDblParam(CbcModel::CbcIntegerTolerance);
860  int i;
861  bool beforeSolution = model->getSolutionCount()==0;
862  int numberStrong=model->numberStrong();
863  int numberObjects = model->numberObjects();
864  bool checkFeasibility = numberObjects>model->numberIntegers();
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 = 1.0e-100;
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) {
996        // Increase estimated degradation to solution
997        estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate());
998        numberUnsatisfied_++;
999        // Better priority? Flush choices.
1000        if (priorityLevel<bestPriority) {
1001          int j;
1002          iSmallest=0;
1003          for (j=0;j<maximumStrong;j++) {
1004            choice[j].upMovement=0.0;
1005            delete choice[j].possibleBranch;
1006            choice[j].possibleBranch=NULL;
1007          }
1008          bestPriority = priorityLevel;
1009          mostAway=1.0e-100;
1010          numberStrong=0;
1011        } else if (priorityLevel>bestPriority) {
1012          continue;
1013        }
1014        // Check for suitability based on infeasibility.
1015        if (infeasibility>mostAway) {
1016          //add to list
1017          choice[iSmallest].upMovement=infeasibility;
1018          delete choice[iSmallest].possibleBranch;
1019          choice[iSmallest].possibleBranch=object->createBranch(preferredWay);
1020          numberStrong = CoinMax(numberStrong,iSmallest+1);
1021          // Save which object it was
1022          choice[iSmallest].objectNumber=i;
1023          int j;
1024          iSmallest=-1;
1025          mostAway = 1.0e50;
1026          for (j=0;j<maximumStrong;j++) {
1027            if (choice[j].upMovement<mostAway) {
1028              mostAway=choice[j].upMovement;
1029              iSmallest=j;
1030            }
1031          }
1032        }
1033      }
1034    }
1035    if (numberUnsatisfied_) {
1036      // some infeasibilities - go to next steps
1037      break;
1038    } else if (!iPass) {
1039      // looks like a solution - get paranoid
1040      bool roundAgain=false;
1041      // get basis
1042      CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
1043      if (!ws)
1044        break;
1045      for (i=0;i<numberColumns;i++) {
1046        double value = saveSolution[i];
1047        if (value<lower[i]) {
1048          saveSolution[i]=lower[i];
1049          roundAgain=true;
1050          ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
1051        } else if (value>upper[i]) {
1052          saveSolution[i]=upper[i];
1053          roundAgain=true;
1054          ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
1055        } 
1056      }
1057      if (roundAgain) {
1058        // restore basis
1059        solver->setWarmStart(ws);
1060        delete ws;
1061        solver->resolve();
1062        memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1063        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  // If we have hit max time don't do strong branching
1096  bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
1097                        model->getDblParam(CbcModel::CbcMaximumSeconds));
1098  // also give up if we are looping round too much
1099  if (hitMaxTime||numberPassesLeft<=0)
1100    numberStrong=0;
1101/*
1102  Is strong branching enabled? If so, set up and do it. Otherwise, we'll
1103  fall through to simple branching.
1104
1105  Setup for strong branching involves saving the current basis (for restoration
1106  afterwards) and setting up for hot starts.
1107*/
1108  if (numberStrong&&model->numberStrong()) {
1109   
1110    bool solveAll=false; // set true to say look at all even if some fixed (experiment)
1111    // worth trying if too many times
1112    // Save basis
1113    CoinWarmStart * ws = solver->getWarmStart();
1114    // save limit
1115    int saveLimit;
1116    solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
1117    if (beforeSolution)
1118      solver->setIntParam(OsiMaxNumIterationHotStart,10000); // go to end
1119
1120    /* If we are doing all strong branching in one go then we create new arrays
1121       to store information.  If clp NULL then doing old way.
1122       Going down -
1123       outputSolution[2*i] is final solution.
1124       outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown
1125       outputStuff[2*i+numberStrong] is number iterations
1126       On entry newUpper[i] is new upper bound, on exit obj change
1127       Going up -
1128       outputSolution[2*i+1] is final solution.
1129       outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
1130       outputStuff[2*i+1+numberStrong] is number iterations
1131       On entry newLower[i] is new lower bound, on exit obj change
1132    */
1133    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
1134    ClpSimplex * clp=NULL;
1135    double * newLower = NULL;
1136    double * newUpper = NULL;
1137    double ** outputSolution=NULL;
1138    int * outputStuff=NULL;
1139    int saveLogLevel=0;
1140    // Go back to normal way if user wants it
1141    if (osiclp&&(osiclp->specialOptions()&16)!=0)
1142      allNormal=false;
1143    if (osiclp&& allNormal) {
1144      clp = osiclp->getModelPtr();
1145      saveLogLevel = clp->logLevel();
1146      clp->setLogLevel(0);
1147      int saveOptions = clp->specialOptions();
1148      int startFinishOptions;
1149      int specialOptions = osiclp->specialOptions();
1150      if((specialOptions&1)==0) {
1151        startFinishOptions=0;
1152        clp->setSpecialOptions(saveOptions|(64|1024));
1153      } else {
1154        startFinishOptions=1+2+4;
1155        //startFinishOptions=1+4; // for moment re-factorize
1156        if((specialOptions&4)==0) 
1157          clp->setSpecialOptions(saveOptions|(64|128|512|1024|4096));
1158        else
1159          clp->setSpecialOptions(saveOptions|(64|128|512|1024|2048|4096));
1160      }
1161      // User may want to clean up before strong branching
1162      if ((clp->specialOptions()&32)!=0) {
1163        clp->primal(1);
1164        if (clp->numberIterations())
1165          model->messageHandler()->message(CBC_ITERATE_STRONG,model->messages())
1166            << clp->numberIterations()
1167            <<CoinMessageEol;
1168      }
1169      int saveMaxIts = clp->maximumIterations();
1170      clp->setMaximumIterations(saveLimit);
1171      // Clp - do a different way
1172      newLower = new double[numberStrong];
1173      newUpper = new double[numberStrong];
1174      outputSolution = new double * [2*numberStrong];
1175      outputStuff = new int [4*numberStrong];
1176      int * which = new int[numberStrong];
1177      for (i=0;i<numberStrong;i++) {
1178        int iObject = choice[i].objectNumber;
1179        const CbcObject * object = model->object(iObject);
1180        const CbcSimpleInteger * simple = dynamic_cast <const CbcSimpleInteger *> (object);
1181        int iSequence = simple->modelSequence();
1182        newLower[i]= ceil(saveSolution[iSequence]);
1183        newUpper[i]= floor(saveSolution[iSequence]);
1184        which[i]=iSequence;
1185        outputSolution[2*i]= new double [numberColumns];
1186        outputSolution[2*i+1]= new double [numberColumns];
1187      }
1188      int returnCode=clp->strongBranching(numberStrong,which,
1189                                          newLower, newUpper,outputSolution,
1190                                          outputStuff,outputStuff+2*numberStrong,!solveAll,false,
1191                                          startFinishOptions);
1192      clp->setSpecialOptions(saveOptions); // restore
1193      clp->setMaximumIterations(saveMaxIts);
1194      if (returnCode==-2) {
1195        // bad factorization!!!
1196        // Doing normal way
1197        // Mark hot start
1198        solver->markHotStart();
1199        clp->setLogLevel(saveLogLevel);
1200        clp = NULL;
1201      }
1202      delete [] which;
1203    } else {
1204      // Doing normal way
1205      // Mark hot start
1206      solver->markHotStart();
1207    }
1208/*
1209  Open a loop to do the strong branching LPs. For each candidate variable,
1210  solve an LP with the variable forced down, then up. If a direction turns
1211  out to be infeasible or monotonic (i.e., over the dual objective cutoff),
1212  force the objective change to be big (1.0e100). If we determine the problem
1213  is infeasible, or find a monotone variable, escape the loop.
1214
1215  TODO: The `restore bounds' part might be better encapsulated as an
1216        unbranch() method. Branching objects more exotic than simple integers
1217        or cliques might not restrict themselves to variable bounds.
1218
1219  TODO: Virtuous solvers invalidate the current solution (or give bogus
1220        results :-) when the bounds are changed out from under them. So we
1221        need to do all the work associated with finding a new solution before
1222        restoring the bounds.
1223*/
1224    for (i = 0 ; i < numberStrong ; i++)
1225    { double objectiveChange ;
1226      double newObjectiveValue=1.0e100;
1227      // status is 0 finished, 1 infeasible and other
1228      int iStatus;
1229/*
1230  Try the down direction first. (Specify the initial branching alternative as
1231  down with a call to way(-1). Each subsequent call to branch() performs the
1232  specified branch and advances the branch object state to the next branch
1233  alternative.)
1234*/
1235      if (!clp) {
1236        choice[i].possibleBranch->way(-1) ;
1237        choice[i].possibleBranch->branch() ;
1238        bool feasible=true;
1239        if (checkFeasibility) {
1240          // check branching did not make infeasible
1241          int iColumn;
1242          int numberColumns = solver->getNumCols();
1243          const double * columnLower = solver->getColLower();
1244          const double * columnUpper = solver->getColUpper();
1245          for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1246            if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1247              feasible=false;
1248          }
1249        }
1250        if (feasible) {
1251          solver->solveFromHotStart() ;
1252/*
1253  We now have an estimate of objective degradation that we can use for strong
1254  branching. If we're over the cutoff, the variable is monotone up.
1255  If we actually made it to optimality, check for a solution, and if we have
1256  a good one, call setBestSolution to process it. Note that this may reduce the
1257  cutoff, so we check again to see if we can declare this variable monotone.
1258*/
1259          if (solver->isProvenOptimal())
1260            iStatus=0; // optimal
1261          else if (solver->isIterationLimitReached()
1262                   &&!solver->isDualObjectiveLimitReached())
1263            iStatus=2; // unknown
1264          else
1265            iStatus=1; // infeasible
1266          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1267          choice[i].numItersDown = solver->getIterationCount();
1268        } else {
1269          iStatus=1; // infeasible
1270          newObjectiveValue = 1.0e100;
1271          choice[i].numItersDown = 0;
1272        }
1273        objectiveChange = newObjectiveValue-objectiveValue ;
1274      } else {
1275        iStatus = outputStuff[2*i];
1276        choice[i].numItersDown = outputStuff[2*numberStrong+2*i];
1277        newObjectiveValue = objectiveValue+newUpper[i];
1278        solver->setColSolution(outputSolution[2*i]);
1279      }
1280      objectiveChange = newObjectiveValue  - objectiveValue;
1281      if (!iStatus) {
1282        choice[i].finishedDown = true ;
1283        if (newObjectiveValue>=model->getCutoff()) {
1284          objectiveChange = 1.0e100; // say infeasible
1285        } else {
1286          // See if integer solution
1287          if (model->feasibleSolution(choice[i].numIntInfeasDown,
1288                                      choice[i].numObjInfeasDown))
1289            { model->setBestSolution(CBC_STRONGSOL,
1290                                     newObjectiveValue,
1291                                     solver->getColSolution()) ;
1292            if (newObjectiveValue >= model->getCutoff())        //  *new* cutoff
1293              objectiveChange = 1.0e100 ;
1294            }
1295        }
1296      } else if (iStatus==1) {
1297        objectiveChange = 1.0e100 ;
1298      } else {
1299        // Can't say much as we did not finish
1300        choice[i].finishedDown = false ;
1301      }
1302      choice[i].downMovement = objectiveChange ;
1303
1304      // restore bounds
1305      if (!clp)
1306      { for (int j=0;j<numberColumns;j++) {
1307        if (saveLower[j] != lower[j])
1308          solver->setColLower(j,saveLower[j]);
1309        if (saveUpper[j] != upper[j])
1310          solver->setColUpper(j,saveUpper[j]);
1311        }
1312      }
1313      //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1314      //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersDown,
1315      //     choice[i].downMovement,choice[i].finishedDown,choice[i].numIntInfeasDown,
1316      //     choice[i].numObjInfeasDown);
1317
1318      // repeat the whole exercise, forcing the variable up
1319      if (!clp) {
1320        choice[i].possibleBranch->branch();
1321        bool feasible=true;
1322        if (checkFeasibility) {
1323          // check branching did not make infeasible
1324          int iColumn;
1325          int numberColumns = solver->getNumCols();
1326          const double * columnLower = solver->getColLower();
1327          const double * columnUpper = solver->getColUpper();
1328          for (iColumn= 0;iColumn<numberColumns;iColumn++) {
1329            if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
1330              feasible=false;
1331          }
1332        }
1333        if (feasible) {
1334          solver->solveFromHotStart() ;
1335/*
1336  We now have an estimate of objective degradation that we can use for strong
1337  branching. If we're over the cutoff, the variable is monotone up.
1338  If we actually made it to optimality, check for a solution, and if we have
1339  a good one, call setBestSolution to process it. Note that this may reduce the
1340  cutoff, so we check again to see if we can declare this variable monotone.
1341*/
1342          if (solver->isProvenOptimal())
1343            iStatus=0; // optimal
1344          else if (solver->isIterationLimitReached()
1345                   &&!solver->isDualObjectiveLimitReached())
1346            iStatus=2; // unknown
1347          else
1348            iStatus=1; // infeasible
1349          newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1350          choice[i].numItersUp = solver->getIterationCount();
1351        } else {
1352          iStatus=1; // infeasible
1353          newObjectiveValue = 1.0e100;
1354          choice[i].numItersDown = 0;
1355        }
1356        objectiveChange = newObjectiveValue-objectiveValue ;
1357      } else {
1358        iStatus = outputStuff[2*i+1];
1359        choice[i].numItersUp = outputStuff[2*numberStrong+2*i+1];
1360        newObjectiveValue = objectiveValue+newLower[i];
1361        solver->setColSolution(outputSolution[2*i+1]);
1362      }
1363      objectiveChange = newObjectiveValue  - objectiveValue;
1364      if (!iStatus) {
1365        choice[i].finishedUp = true ;
1366        if (newObjectiveValue>=model->getCutoff()) {
1367          objectiveChange = 1.0e100; // say infeasible
1368        } else {
1369          // See if integer solution
1370          if (model->feasibleSolution(choice[i].numIntInfeasUp,
1371                                      choice[i].numObjInfeasUp))
1372            { model->setBestSolution(CBC_STRONGSOL,
1373                                     newObjectiveValue,
1374                                     solver->getColSolution()) ;
1375            if (newObjectiveValue >= model->getCutoff())        //  *new* cutoff
1376              objectiveChange = 1.0e100 ;
1377            }
1378        }
1379      } else if (iStatus==1) {
1380        objectiveChange = 1.0e100 ;
1381      } else {
1382        // Can't say much as we did not finish
1383        choice[i].finishedUp = false ;
1384      }
1385      choice[i].upMovement = objectiveChange ;
1386
1387      // restore bounds
1388      if (!clp)
1389      { for (int j=0;j<numberColumns;j++) {
1390        if (saveLower[j] != lower[j])
1391          solver->setColLower(j,saveLower[j]);
1392        if (saveUpper[j] != upper[j])
1393          solver->setColUpper(j,saveUpper[j]);
1394        }
1395      }
1396
1397      //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1398      //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersUp,
1399      //     choice[i].upMovement,choice[i].finishedUp,choice[i].numIntInfeasUp,
1400      //     choice[i].numObjInfeasUp);
1401
1402/*
1403  End of evaluation for this candidate variable. Possibilities are:
1404    * Both sides below cutoff; this variable is a candidate for branching.
1405    * Both sides infeasible or above the objective cutoff: no further action
1406      here. Break from the evaluation loop and assume the node will be purged
1407      by the caller.
1408    * One side below cutoff: Install the branch (i.e., fix the variable). Break
1409      from the evaluation loop and assume the node will be reoptimised by the
1410      caller.
1411*/
1412      if (choice[i].upMovement<1.0e100) {
1413        if(choice[i].downMovement<1.0e100) {
1414          // feasible - no action
1415        } else {
1416          // up feasible, down infeasible
1417          anyAction=-1;
1418          if (!solveAll) {
1419            choice[i].possibleBranch->way(1);
1420            choice[i].possibleBranch->branch();
1421            break;
1422          }
1423        }
1424      } else {
1425        if(choice[i].downMovement<1.0e100) {
1426          // down feasible, up infeasible
1427          anyAction=-1;
1428          if (!solveAll) {
1429            choice[i].possibleBranch->way(-1);
1430            choice[i].possibleBranch->branch();
1431            break;
1432          }
1433        } else {
1434          // neither side feasible
1435          anyAction=-2;
1436          break;
1437        }
1438      }
1439    }
1440
1441    if (!clp) {
1442      // Delete the snapshot
1443      solver->unmarkHotStart();
1444    } else {
1445      clp->setLogLevel(saveLogLevel);
1446      delete [] newLower;
1447      delete [] newUpper;
1448      delete [] outputStuff;
1449      int i;
1450      for (i=0;i<2*numberStrong;i++)
1451        delete [] outputSolution[i];
1452      delete [] outputSolution;
1453    }
1454    solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
1455    // restore basis
1456    solver->setWarmStart(ws);
1457    delete ws;
1458
1459/*
1460  anyAction >= 0 indicates that strong branching didn't produce any monotone
1461  variables. Sift through the candidates for the best one.
1462
1463  QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is
1464         local to this code block. Perhaps should be numberNodes_ from model?
1465         Unclear what this calculation is doing.
1466*/
1467    if (anyAction>=0) {
1468
1469      int numberNodes = model->getNodeCount();
1470      // get average cost per iteration and assume stopped ones
1471      // would stop after 50% more iterations at average cost??? !!! ???
1472      double averageCostPerIteration=0.0;
1473      double totalNumberIterations=1.0;
1474      int smallestNumberInfeasibilities=INT_MAX;
1475      for (i=0;i<numberStrong;i++) {
1476        totalNumberIterations += choice[i].numItersDown +
1477          choice[i].numItersUp ;
1478        averageCostPerIteration += choice[i].downMovement +
1479          choice[i].upMovement;
1480        smallestNumberInfeasibilities= 
1481          CoinMin(CoinMin(choice[i].numIntInfeasDown ,
1482                  choice[i].numIntInfeasUp ),
1483              smallestNumberInfeasibilities);
1484      }
1485      if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities)
1486        numberNodes=1000000; // switch off search for better solution
1487      numberNodes=1000000; // switch off anyway
1488      averageCostPerIteration /= totalNumberIterations;
1489      // all feasible - choose best bet
1490
1491#if 0
1492      for (i = 0 ; i < numberStrong ; i++)
1493      { int iColumn =
1494          model->integerVariable()[choice[i].possibleBranch->variable()] ;
1495         model->messageHandler()->message(CBC_STRONG,model->messages())
1496          << i << iColumn
1497          <<choice[i].downMovement<<choice[i].numIntInfeasDown
1498          <<choice[i].upMovement<<choice[i].numIntInfeasUp
1499          <<choice[i].possibleBranch->value()
1500          <<CoinMessageEol;
1501        int betterWay = decision->betterBranch(choice[i].possibleBranch,
1502                                              branch_,
1503                                              choice[i].upMovement,
1504                                              choice[i].numIntInfeasUp ,
1505                                              choice[i].downMovement,
1506                                              choice[i].numIntInfeasDown );
1507        if (betterWay) {
1508          delete branch_;
1509          // C) create branching object
1510          branch_ = choice[i].possibleBranch;
1511          choice[i].possibleBranch=NULL;
1512          branch_->way(betterWay);
1513        }
1514      }
1515#else
1516      // New method does all at once so it can be more sophisticated
1517      // in deciding how to balance actions.
1518      // But it does need arrays
1519      double * changeUp = new double [numberStrong];
1520      int * numberInfeasibilitiesUp = new int [numberStrong];
1521      double * changeDown = new double [numberStrong];
1522      int * numberInfeasibilitiesDown = new int [numberStrong];
1523      CbcBranchingObject ** objects = new CbcBranchingObject * [ numberStrong];
1524      for (i = 0 ; i < numberStrong ; i++) {
1525        int iColumn = choice[i].possibleBranch->variable() ;
1526        model->messageHandler()->message(CBC_STRONG,model->messages())
1527          << i << iColumn
1528          <<choice[i].downMovement<<choice[i].numIntInfeasDown
1529          <<choice[i].upMovement<<choice[i].numIntInfeasUp
1530          <<choice[i].possibleBranch->value()
1531          <<CoinMessageEol;
1532        changeUp[i]=choice[i].upMovement;
1533        numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
1534        changeDown[i]=choice[i].downMovement;
1535        numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
1536        objects[i] = choice[i].possibleBranch;
1537      }
1538      int whichObject = decision->bestBranch(objects,numberStrong,numberUnsatisfied_,
1539                                             changeUp,numberInfeasibilitiesUp,
1540                                             changeDown,numberInfeasibilitiesDown,
1541                                             objectiveValue);
1542      // move branching object and make sure it will not be deleted
1543      if (whichObject>=0) {
1544        branch_ = objects[whichObject];
1545        choice[whichObject].possibleBranch=NULL;
1546      }
1547      delete [] changeUp;
1548      delete [] numberInfeasibilitiesUp;
1549      delete [] changeDown;
1550      delete [] numberInfeasibilitiesDown;
1551      delete [] objects;
1552#endif
1553    }
1554  }
1555/*
1556  Simple branching. Probably just one, but we may have got here
1557  because of an odd branch e.g. a cut
1558*/
1559  else {
1560    // not strong
1561    // C) create branching object
1562    branch_ = choice[bestChoice].possibleBranch;
1563    choice[bestChoice].possibleBranch=NULL;
1564  }
1565  // Set guessed solution value
1566  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
1567/*
1568  Cleanup, then we're outta here.
1569*/
1570  if (!model->branchingMethod())
1571    delete decision;
1572   
1573  for (i=0;i<maximumStrong;i++)
1574    delete choice[i].possibleBranch;
1575  delete [] choice;
1576  delete [] saveLower;
1577  delete [] saveUpper;
1578 
1579  // restore solution
1580  solver->setColSolution(saveSolution);
1581  delete [] saveSolution;
1582  return anyAction;
1583}
1584
1585
1586CbcNode::CbcNode(const CbcNode & rhs) 
1587{ 
1588#ifdef CHECK_NODE
1589  printf("CbcNode %x Constructor from rhs %x\n",this,&rhs);
1590#endif
1591  if (rhs.nodeInfo_)
1592    nodeInfo_ = rhs.nodeInfo_->clone();
1593  else
1594    nodeInfo_=NULL;
1595  objectiveValue_=rhs.objectiveValue_;
1596  guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
1597  if (rhs.branch_)
1598    branch_=rhs.branch_->clone();
1599  else
1600    branch_=NULL;
1601  depth_ = rhs.depth_;
1602  numberUnsatisfied_ = rhs.numberUnsatisfied_;
1603  nodeNumber_=rhs.nodeNumber_;
1604}
1605
1606CbcNode &
1607CbcNode::operator=(const CbcNode & rhs)
1608{
1609  if (this != &rhs) {
1610    delete nodeInfo_;
1611    if (nodeInfo_)
1612      nodeInfo_ = rhs.nodeInfo_->clone();
1613    else
1614      nodeInfo_ = NULL;
1615    objectiveValue_=rhs.objectiveValue_;
1616    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
1617    if (rhs.branch_)
1618      branch_=rhs.branch_->clone();
1619    else
1620      branch_=NULL,
1621    depth_ = rhs.depth_;
1622    numberUnsatisfied_ = rhs.numberUnsatisfied_;
1623    nodeNumber_=rhs.nodeNumber_;
1624  }
1625  return *this;
1626}
1627
1628
1629CbcNode::~CbcNode ()
1630{
1631#ifdef CHECK_NODE
1632  if (nodeInfo_) 
1633    printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
1634         this,nodeInfo_,nodeInfo_->numberPointingToThis());
1635  else
1636    printf("CbcNode %x Destructor nodeInfo %x (?)\n",
1637         this,nodeInfo_);
1638#endif
1639  if (nodeInfo_) {
1640    nodeInfo_->nullOwner();
1641    int numberToDelete=nodeInfo_->numberBranchesLeft();
1642    //    CbcNodeInfo * parent = nodeInfo_->parent();
1643    if (nodeInfo_->decrement(numberToDelete)==0) {
1644      delete nodeInfo_;
1645    } else {
1646      //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,parent);
1647      // anyway decrement parent
1648      //if (parent)
1649      ///parent->decrement(1);
1650    }
1651  }
1652  delete branch_;
1653}
1654// Decrement  active cut counts
1655void 
1656CbcNode::decrementCuts(int change)
1657{
1658  if(nodeInfo_) {
1659    nodeInfo_->decrementCuts(change);
1660  }
1661}
1662void 
1663CbcNode::decrementParentCuts(int change)
1664{
1665  if(nodeInfo_) {
1666    nodeInfo_->decrementParentCuts(change);
1667  }
1668}
1669
1670/*
1671  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
1672  in the attached nodeInfo_.
1673*/
1674void
1675CbcNode::initializeInfo()
1676{
1677  assert(nodeInfo_ && branch_) ;
1678  nodeInfo_->initializeInfo(branch_->numberBranches());
1679}
1680// Nulls out node info
1681void 
1682CbcNode::nullNodeInfo()
1683{
1684  nodeInfo_=NULL;
1685}
1686
1687int
1688CbcNode::branch()
1689{
1690  double changeInGuessed=branch_->branch(true);
1691  guessedObjectiveValue_+= changeInGuessed;
1692  return nodeInfo_->branchedOn();
1693}
Note: See TracBrowser for help on using the repository browser.