source: trunk/CbcNode.cpp @ 48

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

skip messages etc

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 50.8 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  int maximumStrong = CoinMax(CoinMin(model->numberStrong(),numberObjects),1);
865  int numberColumns = model->getNumCols();
866  double * saveUpper = new double[numberColumns];
867  double * saveLower = new double[numberColumns];
868
869  // Save solution in case heuristics need good solution later
870
871  double * saveSolution = new double[numberColumns];
872  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
873
874/*
875  Get a branching decision object. Use the default decision criteria unless
876  the user has loaded a decision method into the model.
877*/
878  CbcBranchDecision *decision = model->branchingMethod();
879  if (!decision)
880    decision = new CbcBranchDefaultDecision();
881
882  typedef struct {
883    CbcBranchingObject * possibleBranch; // what a branch would do
884    double upMovement; // cost going up (and initial away from feasible)
885    double downMovement; // cost going down
886    int numIntInfeasUp ; // without odd ones
887    int numObjInfeasUp ; // just odd ones
888    bool finishedUp; // true if solver finished
889    int numItersUp ; // number of iterations in solver
890    int numIntInfeasDown ; // without odd ones
891    int numObjInfeasDown ; // just odd ones
892    bool finishedDown; // true if solver finished
893    int numItersDown; // number of iterations in solver
894    int objectNumber; // Which object it is
895  } Strong;
896  Strong * choice = new Strong[maximumStrong];
897  for (i=0;i<numberColumns;i++) {
898    saveLower[i] = lower[i];
899    saveUpper[i] = upper[i];
900  }
901  // Some objects may compute an estimate of best solution from here
902  double estimatedDegradation=0.0; 
903  int numberIntegerInfeasibilities=0; // without odd ones
904
905  // We may go round this loop twice (only if we think we have solution)
906  for (int iPass=0;iPass<2;iPass++) {
907
908    // compute current state
909    int numberObjectInfeasibilities; // just odd ones
910    model->feasibleSolution(
911                            numberIntegerInfeasibilities,
912                            numberObjectInfeasibilities);
913    // If forcePriority > 0 then we want best solution
914    const double * bestSolution = NULL;
915    int hotstartStrategy=model->getHotstartStrategy();
916    if (hotstartStrategy>0) {
917      bestSolution = model->bestSolution();
918    }
919   
920    // Some objects may compute an estimate of best solution from here
921    estimatedDegradation=0.0; 
922    numberUnsatisfied_ = 0;
923    int bestPriority=INT_MAX;
924    /*
925      Scan for branching objects that indicate infeasibility. Choose the best
926      maximumStrong candidates, using priority as the first criteria, then
927      integer infeasibility.
928     
929      The algorithm is to fill the choice array with a set of good candidates (by
930      infeasibility) with priority bestPriority.  Finding a candidate with
931      priority better (less) than bestPriority flushes the choice array. (This
932      serves as initialization when the first candidate is found.)
933     
934      A new candidate is added to choices only if its infeasibility exceeds the
935      current max infeasibility (mostAway). When a candidate is added, it
936      replaces the candidate with the smallest infeasibility (tracked by
937      iSmallest).
938    */
939    int iSmallest = 0;
940    double mostAway = integerTolerance;
941    for (i = 0 ; i < maximumStrong ; i++) choice[i].possibleBranch = NULL ;
942    numberStrong=0;
943    for (i=0;i<numberObjects;i++) {
944      const CbcObject * object = model->object(i);
945      int preferredWay;
946      double infeasibility = object->infeasibility(preferredWay);
947      int priorityLevel = model->priority(i);
948      if (bestSolution) {
949        // we are doing hot start
950        const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
951        if (thisOne) {
952          int iColumn = thisOne->modelSequence();
953          if (saveUpper[iColumn]>saveLower[iColumn]) {
954            double value = saveSolution[iColumn];
955            double targetValue = bestSolution[iColumn];
956            //double originalLower = thisOne->originalLower();
957            //double originalUpper = thisOne->originalUpper();
958            // switch off if not possible
959            if (targetValue>=saveLower[iColumn]&&targetValue<=saveUpper[iColumn]) {
960              /* priority outranks rest always if hotstartStrategy >1
961                 otherwise can be downgraded if at correct level.
962                 Infeasibility may be increased by targetValue to choose 1.0 values first.
963              */
964              if (fabs(value-targetValue)>integerTolerance) {
965                if (value>targetValue) {
966                  infeasibility += value;
967                  preferredWay=-1;
968                } else {
969                  infeasibility += targetValue;
970                  preferredWay=1;
971                }
972              } else if (hotstartStrategy>1) {
973                if (targetValue==saveLower[iColumn]) {
974                  infeasibility += integerTolerance+1.0e-12;
975                  preferredWay=-1;
976                } else if (targetValue==saveUpper[iColumn]) {
977                  infeasibility += integerTolerance+1.0e-12;
978                  preferredWay=1;
979                } else {
980                  infeasibility += integerTolerance+1.0e-12;
981                  preferredWay=1;
982                }
983              } else {
984                priorityLevel += 10000000;
985              }
986            } else {
987              // switch off if not possible
988              bestSolution=NULL;
989              model->setHotstartStrategy(0);
990            }
991          }
992        }
993      }
994      if (infeasibility>integerTolerance) {
995        // Increase estimated degradation to solution
996        estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate());
997        numberUnsatisfied_++;
998        // Better priority? Flush choices.
999        if (priorityLevel<bestPriority) {
1000          int j;
1001          iSmallest=0;
1002          for (j=0;j<maximumStrong;j++) {
1003            choice[j].upMovement=0.0;
1004            delete choice[j].possibleBranch;
1005            choice[j].possibleBranch=NULL;
1006          }
1007          bestPriority = priorityLevel;
1008          mostAway=integerTolerance;
1009          numberStrong=0;
1010        } else if (priorityLevel>bestPriority) {
1011          continue;
1012        }
1013        // Check for suitability based on infeasibility.
1014        if (infeasibility>mostAway) {
1015          //add to list
1016          choice[iSmallest].upMovement=infeasibility;
1017          delete choice[iSmallest].possibleBranch;
1018          choice[iSmallest].possibleBranch=object->createBranch(preferredWay);
1019          numberStrong = CoinMax(numberStrong,iSmallest+1);
1020          // Save which object it was
1021          choice[iSmallest].objectNumber=i;
1022          int j;
1023          iSmallest=-1;
1024          mostAway = 1.0e50;
1025          for (j=0;j<maximumStrong;j++) {
1026            if (choice[j].upMovement<mostAway) {
1027              mostAway=choice[j].upMovement;
1028              iSmallest=j;
1029            }
1030          }
1031        }
1032      }
1033    }
1034    if (numberUnsatisfied_) {
1035      // some infeasibilities - go to next steps
1036      break;
1037    } else if (!iPass) {
1038      // looks like a solution - get paranoid
1039      bool roundAgain=false;
1040      // get basis
1041      CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
1042      if (!ws)
1043        break;
1044      for (i=0;i<numberColumns;i++) {
1045        double value = saveSolution[i];
1046        if (value<lower[i]) {
1047          saveSolution[i]=lower[i];
1048          roundAgain=true;
1049          ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
1050        } else if (value>upper[i]) {
1051          saveSolution[i]=upper[i];
1052          roundAgain=true;
1053          ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
1054        } 
1055      }
1056      if (roundAgain) {
1057        // restore basis
1058        solver->setWarmStart(ws);
1059        delete ws;
1060        solver->resolve();
1061        memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
1062        if (!solver->isProvenOptimal()) {
1063          // infeasible
1064          anyAction=-2;
1065          break;
1066        }
1067      } else {
1068        delete ws;
1069        break;
1070      }
1071    }
1072  }
1073  /* Some solvers can do the strong branching calculations faster if
1074     they do them all at once.  At present only Clp does for ordinary
1075     integers but I think this coding would be easy to modify
1076  */
1077  bool allNormal=true; // to say if we can do fast strong branching
1078  // Say which one will be best
1079  int bestChoice=0;
1080  double worstInfeasibility=0.0;
1081  for (i=0;i<numberStrong;i++) {
1082    choice[i].numIntInfeasUp = numberUnsatisfied_;
1083    choice[i].numIntInfeasDown = numberUnsatisfied_;
1084    if (!dynamic_cast <const CbcSimpleInteger *> (model->object(choice[i].objectNumber)))
1085      allNormal=false; // Something odd so lets skip clever fast branching
1086    if ( !model->object(choice[i].objectNumber)->boundBranch())
1087      numberStrong=0; // switch off
1088    // Do best choice in case switched off
1089    if (choice[i].upMovement>worstInfeasibility) {
1090      worstInfeasibility=choice[i].upMovement;
1091      bestChoice=i;
1092    }
1093  }
1094  // If we have hit max time don't do strong branching
1095  bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > 
1096                        model->getDblParam(CbcModel::CbcMaximumSeconds));
1097  // also give up if we are looping round too much
1098  if (hitMaxTime||numberPassesLeft<=0)
1099    numberStrong=0;
1100/*
1101  Is strong branching enabled? If so, set up and do it. Otherwise, we'll
1102  fall through to simple branching.
1103
1104  Setup for strong branching involves saving the current basis (for restoration
1105  afterwards) and setting up for hot starts.
1106*/
1107  if (numberStrong&&model->numberStrong()) {
1108   
1109    bool solveAll=false; // set true to say look at all even if some fixed (experiment)
1110    // worth trying if too many times
1111    // Save basis
1112    CoinWarmStart * ws = solver->getWarmStart();
1113    // save limit
1114    int saveLimit;
1115    solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
1116    if (beforeSolution)
1117      solver->setIntParam(OsiMaxNumIterationHotStart,10000); // go to end
1118
1119    /* If we are doing all strong branching in one go then we create new arrays
1120       to store information.  If clp NULL then doing old way.
1121       Going down -
1122       outputSolution[2*i] is final solution.
1123       outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown
1124       outputStuff[2*i+numberStrong] is number iterations
1125       On entry newUpper[i] is new upper bound, on exit obj change
1126       Going up -
1127       outputSolution[2*i+1] is final solution.
1128       outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
1129       outputStuff[2*i+1+numberStrong] is number iterations
1130       On entry newLower[i] is new lower bound, on exit obj change
1131    */
1132    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
1133    ClpSimplex * clp=NULL;
1134    double * newLower = NULL;
1135    double * newUpper = NULL;
1136    double ** outputSolution=NULL;
1137    int * outputStuff=NULL;
1138    int saveLogLevel=0;
1139    // For moment - until full testing - go back to normal way
1140    //allNormal=false;
1141    if (osiclp&& allNormal) {
1142      clp = osiclp->getModelPtr();
1143      saveLogLevel = clp->logLevel();
1144      clp->setLogLevel(0);
1145      int saveOptions = clp->specialOptions();
1146      int startFinishOptions;
1147      int specialOptions = osiclp->specialOptions();
1148      if((specialOptions&1)==0) {
1149        startFinishOptions=0;
1150        clp->setSpecialOptions(saveOptions|(64|1024));
1151      } else {
1152        startFinishOptions=1+2+4;
1153        //startFinishOptions=1+4; // for moment re-factorize
1154        if((specialOptions&4)==0) 
1155          clp->setSpecialOptions(saveOptions|(64|128|512|1024|4096));
1156        else
1157          clp->setSpecialOptions(saveOptions|(64|128|512|1024|2048|4096));
1158      }
1159      // User may want to clean up before strong branching
1160      if ((clp->specialOptions()&32)!=0) {
1161        clp->primal(1);
1162        if (clp->numberIterations())
1163          model->messageHandler()->message(CBC_ITERATE_STRONG,model->messages())
1164            << clp->numberIterations()
1165            <<CoinMessageEol;
1166      }
1167      int saveMaxIts = clp->maximumIterations();
1168      clp->setMaximumIterations(saveLimit);
1169      // Clp - do a different way
1170      newLower = new double[numberStrong];
1171      newUpper = new double[numberStrong];
1172      outputSolution = new double * [2*numberStrong];
1173      outputStuff = new int [4*numberStrong];
1174      int * which = new int[numberStrong];
1175      for (i=0;i<numberStrong;i++) {
1176        int iObject = choice[i].objectNumber;
1177        const CbcObject * object = model->object(iObject);
1178        const CbcSimpleInteger * simple = dynamic_cast <const CbcSimpleInteger *> (object);
1179        int iSequence = simple->modelSequence();
1180        newLower[i]= ceil(saveSolution[iSequence]);
1181        newUpper[i]= floor(saveSolution[iSequence]);
1182        which[i]=iSequence;
1183        outputSolution[2*i]= new double [numberColumns];
1184        outputSolution[2*i+1]= new double [numberColumns];
1185      }
1186      int returnCode=clp->strongBranching(numberStrong,which,
1187                                          newLower, newUpper,outputSolution,
1188                                          outputStuff,outputStuff+2*numberStrong,!solveAll,false,
1189                                          startFinishOptions);
1190      clp->setSpecialOptions(saveOptions); // restore
1191      clp->setMaximumIterations(saveMaxIts);
1192      if (returnCode==-2) {
1193        // bad factorization!!!
1194        // Doing normal way
1195        // Mark hot start
1196        solver->markHotStart();
1197        clp = NULL;
1198      }
1199      delete [] which;
1200    } else {
1201      // Doing normal way
1202      // Mark hot start
1203      solver->markHotStart();
1204    }
1205/*
1206  Open a loop to do the strong branching LPs. For each candidate variable,
1207  solve an LP with the variable forced down, then up. If a direction turns
1208  out to be infeasible or monotonic (i.e., over the dual objective cutoff),
1209  force the objective change to be big (1.0e100). If we determine the problem
1210  is infeasible, or find a monotone variable, escape the loop.
1211
1212  TODO: The `restore bounds' part might be better encapsulated as an
1213        unbranch() method. Branching objects more exotic than simple integers
1214        or cliques might not restrict themselves to variable bounds.
1215
1216  TODO: Virtuous solvers invalidate the current solution (or give bogus
1217        results :-) when the bounds are changed out from under them. So we
1218        need to do all the work associated with finding a new solution before
1219        restoring the bounds.
1220*/
1221    for (i = 0 ; i < numberStrong ; i++)
1222    { double objectiveChange ;
1223      double newObjectiveValue=1.0e100;
1224      // status is 0 finished, 1 infeasible and other
1225      int iStatus;
1226/*
1227  Try the down direction first. (Specify the initial branching alternative as
1228  down with a call to way(-1). Each subsequent call to branch() performs the
1229  specified branch and advances the branch object state to the next branch
1230  alternative.)
1231*/
1232      if (!clp) {
1233        choice[i].possibleBranch->way(-1) ;
1234        choice[i].possibleBranch->branch() ;
1235        solver->solveFromHotStart() ;
1236/*
1237  We now have an estimate of objective degradation that we can use for strong
1238  branching. If we're over the cutoff, the variable is monotone up.
1239  If we actually made it to optimality, check for a solution, and if we have
1240  a good one, call setBestSolution to process it. Note that this may reduce the
1241  cutoff, so we check again to see if we can declare this variable monotone.
1242*/
1243        if (solver->isProvenOptimal())
1244          iStatus=0; // optimal
1245        else if (solver->isIterationLimitReached()
1246                 &&!solver->isDualObjectiveLimitReached())
1247          iStatus=2; // unknown
1248        else
1249          iStatus=1; // infeasible
1250        newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1251        choice[i].numItersDown = solver->getIterationCount();
1252        objectiveChange = newObjectiveValue-objectiveValue ;
1253      } else {
1254        iStatus = outputStuff[2*i];
1255        choice[i].numItersDown = outputStuff[2*numberStrong+2*i];
1256        newObjectiveValue = objectiveValue+newUpper[i];
1257        solver->setColSolution(outputSolution[2*i]);
1258      }
1259      objectiveChange = newObjectiveValue  - objectiveValue;
1260      if (!iStatus) {
1261        choice[i].finishedDown = true ;
1262        if (newObjectiveValue>=model->getCutoff()) {
1263          objectiveChange = 1.0e100; // say infeasible
1264        } else {
1265          // See if integer solution
1266          if (model->feasibleSolution(choice[i].numIntInfeasDown,
1267                                      choice[i].numObjInfeasDown))
1268            { model->setBestSolution(CBC_STRONGSOL,
1269                                     newObjectiveValue,
1270                                     solver->getColSolution()) ;
1271            if (newObjectiveValue >= model->getCutoff())        //  *new* cutoff
1272              objectiveChange = 1.0e100 ;
1273            }
1274        }
1275      } else if (iStatus==1) {
1276        objectiveChange = 1.0e100 ;
1277      } else {
1278        // Can't say much as we did not finish
1279        choice[i].finishedDown = false ;
1280      }
1281      choice[i].downMovement = objectiveChange ;
1282
1283      // restore bounds
1284      if (!clp)
1285      { for (int j=0;j<numberColumns;j++) {
1286        if (saveLower[j] != lower[j])
1287          solver->setColLower(j,saveLower[j]);
1288        if (saveUpper[j] != upper[j])
1289          solver->setColUpper(j,saveUpper[j]);
1290        }
1291      }
1292      //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1293      //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersDown,
1294      //     choice[i].downMovement,choice[i].finishedDown,choice[i].numIntInfeasDown,
1295      //     choice[i].numObjInfeasDown);
1296
1297      // repeat the whole exercise, forcing the variable up
1298      if (!clp) {
1299        choice[i].possibleBranch->branch();
1300        solver->solveFromHotStart();
1301        if (solver->isProvenOptimal())
1302          iStatus=0; // optimal
1303        else if (solver->isIterationLimitReached()
1304                 &&!solver->isDualObjectiveLimitReached())
1305          iStatus=2; // unknown
1306        else
1307          iStatus=1; // infeasible
1308        newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1309        choice[i].numItersUp = solver->getIterationCount();
1310        objectiveChange = newObjectiveValue-objectiveValue ;
1311      } else {
1312        iStatus = outputStuff[2*i+1];
1313        choice[i].numItersUp = outputStuff[2*numberStrong+2*i+1];
1314        newObjectiveValue = objectiveValue+newLower[i];
1315        solver->setColSolution(outputSolution[2*i+1]);
1316      }
1317      objectiveChange = newObjectiveValue  - objectiveValue;
1318      if (!iStatus) {
1319        choice[i].finishedUp = true ;
1320        if (newObjectiveValue>=model->getCutoff()) {
1321          objectiveChange = 1.0e100; // say infeasible
1322        } else {
1323          // See if integer solution
1324          if (model->feasibleSolution(choice[i].numIntInfeasUp,
1325                                      choice[i].numObjInfeasUp))
1326            { model->setBestSolution(CBC_STRONGSOL,
1327                                     newObjectiveValue,
1328                                     solver->getColSolution()) ;
1329            if (newObjectiveValue >= model->getCutoff())        //  *new* cutoff
1330              objectiveChange = 1.0e100 ;
1331            }
1332        }
1333      } else if (iStatus==1) {
1334        objectiveChange = 1.0e100 ;
1335      } else {
1336        // Can't say much as we did not finish
1337        choice[i].finishedUp = false ;
1338      }
1339      choice[i].upMovement = objectiveChange ;
1340
1341      // restore bounds
1342      if (!clp)
1343      { for (int j=0;j<numberColumns;j++) {
1344        if (saveLower[j] != lower[j])
1345          solver->setColLower(j,saveLower[j]);
1346        if (saveUpper[j] != upper[j])
1347          solver->setColUpper(j,saveUpper[j]);
1348        }
1349      }
1350
1351      //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
1352      //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersUp,
1353      //     choice[i].upMovement,choice[i].finishedUp,choice[i].numIntInfeasUp,
1354      //     choice[i].numObjInfeasUp);
1355
1356/*
1357  End of evaluation for this candidate variable. Possibilities are:
1358    * Both sides below cutoff; this variable is a candidate for branching.
1359    * Both sides infeasible or above the objective cutoff: no further action
1360      here. Break from the evaluation loop and assume the node will be purged
1361      by the caller.
1362    * One side below cutoff: Install the branch (i.e., fix the variable). Break
1363      from the evaluation loop and assume the node will be reoptimised by the
1364      caller.
1365*/
1366      if (choice[i].upMovement<1.0e100) {
1367        if(choice[i].downMovement<1.0e100) {
1368          // feasible - no action
1369        } else {
1370          // up feasible, down infeasible
1371          anyAction=-1;
1372          if (!solveAll) {
1373            choice[i].possibleBranch->way(1);
1374            choice[i].possibleBranch->branch();
1375            break;
1376          }
1377        }
1378      } else {
1379        if(choice[i].downMovement<1.0e100) {
1380          // down feasible, up infeasible
1381          anyAction=-1;
1382          if (!solveAll) {
1383            choice[i].possibleBranch->way(-1);
1384            choice[i].possibleBranch->branch();
1385            break;
1386          }
1387        } else {
1388          // neither side feasible
1389          anyAction=-2;
1390          break;
1391        }
1392      }
1393    }
1394
1395    if (!clp) {
1396      // Delete the snapshot
1397      solver->unmarkHotStart();
1398    } else {
1399      clp->setLogLevel(saveLogLevel);
1400      delete [] newLower;
1401      delete [] newUpper;
1402      delete [] outputStuff;
1403      int i;
1404      for (i=0;i<2*numberStrong;i++)
1405        delete [] outputSolution[i];
1406      delete [] outputSolution;
1407    }
1408    solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
1409    // restore basis
1410    solver->setWarmStart(ws);
1411    delete ws;
1412
1413/*
1414  anyAction >= 0 indicates that strong branching didn't produce any monotone
1415  variables. Sift through the candidates for the best one.
1416
1417  QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is
1418         local to this code block. Perhaps should be numberNodes_ from model?
1419         Unclear what this calculation is doing.
1420*/
1421    if (anyAction>=0) {
1422
1423      int numberNodes = model->getNodeCount();
1424      // get average cost per iteration and assume stopped ones
1425      // would stop after 50% more iterations at average cost??? !!! ???
1426      double averageCostPerIteration=0.0;
1427      double totalNumberIterations=1.0;
1428      int smallestNumberInfeasibilities=INT_MAX;
1429      for (i=0;i<numberStrong;i++) {
1430        totalNumberIterations += choice[i].numItersDown +
1431          choice[i].numItersUp ;
1432        averageCostPerIteration += choice[i].downMovement +
1433          choice[i].upMovement;
1434        smallestNumberInfeasibilities= 
1435          CoinMin(CoinMin(choice[i].numIntInfeasDown ,
1436                  choice[i].numIntInfeasUp ),
1437              smallestNumberInfeasibilities);
1438      }
1439      if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities)
1440        numberNodes=1000000; // switch off search for better solution
1441      numberNodes=1000000; // switch off anyway
1442      averageCostPerIteration /= totalNumberIterations;
1443      // all feasible - choose best bet
1444
1445#if 0
1446      for (i = 0 ; i < numberStrong ; i++)
1447      { int iColumn =
1448          model->integerVariable()[choice[i].possibleBranch->variable()] ;
1449         model->messageHandler()->message(CBC_STRONG,model->messages())
1450          << i << iColumn
1451          <<choice[i].downMovement<<choice[i].numIntInfeasDown
1452          <<choice[i].upMovement<<choice[i].numIntInfeasUp
1453          <<choice[i].possibleBranch->value()
1454          <<CoinMessageEol;
1455        int betterWay = decision->betterBranch(choice[i].possibleBranch,
1456                                              branch_,
1457                                              choice[i].upMovement,
1458                                              choice[i].numIntInfeasUp ,
1459                                              choice[i].downMovement,
1460                                              choice[i].numIntInfeasDown );
1461        if (betterWay) {
1462          delete branch_;
1463          // C) create branching object
1464          branch_ = choice[i].possibleBranch;
1465          choice[i].possibleBranch=NULL;
1466          branch_->way(betterWay);
1467        }
1468      }
1469#else
1470      // New method does all at once so it can be more sophisticated
1471      // in deciding how to balance actions.
1472      // But it does need arrays
1473      double * changeUp = new double [numberStrong];
1474      int * numberInfeasibilitiesUp = new int [numberStrong];
1475      double * changeDown = new double [numberStrong];
1476      int * numberInfeasibilitiesDown = new int [numberStrong];
1477      CbcBranchingObject ** objects = new CbcBranchingObject * [ numberStrong];
1478      for (i = 0 ; i < numberStrong ; i++) {
1479        int iColumn = choice[i].possibleBranch->variable() ;
1480        model->messageHandler()->message(CBC_STRONG,model->messages())
1481          << i << iColumn
1482          <<choice[i].downMovement<<choice[i].numIntInfeasDown
1483          <<choice[i].upMovement<<choice[i].numIntInfeasUp
1484          <<choice[i].possibleBranch->value()
1485          <<CoinMessageEol;
1486        changeUp[i]=choice[i].upMovement;
1487        numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
1488        changeDown[i]=choice[i].downMovement;
1489        numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
1490        objects[i] = choice[i].possibleBranch;
1491      }
1492      int whichObject = decision->bestBranch(objects,numberStrong,numberUnsatisfied_,
1493                                             changeUp,numberInfeasibilitiesUp,
1494                                             changeDown,numberInfeasibilitiesDown,
1495                                             objectiveValue);
1496      // move branching object and make sure it will not be deleted
1497      if (whichObject>=0) {
1498        branch_ = objects[whichObject];
1499        choice[whichObject].possibleBranch=NULL;
1500      }
1501      delete [] changeUp;
1502      delete [] numberInfeasibilitiesUp;
1503      delete [] changeDown;
1504      delete [] numberInfeasibilitiesDown;
1505      delete [] objects;
1506#endif
1507    }
1508  }
1509/*
1510  Simple branching. Probably just one, but we may have got here
1511  because of an odd branch e.g. a cut
1512*/
1513  else {
1514    // not strong
1515    // C) create branching object
1516    branch_ = choice[bestChoice].possibleBranch;
1517    choice[bestChoice].possibleBranch=NULL;
1518  }
1519  // Set guessed solution value
1520  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
1521/*
1522  Cleanup, then we're outta here.
1523*/
1524  if (!model->branchingMethod())
1525    delete decision;
1526   
1527  for (i=0;i<maximumStrong;i++)
1528    delete choice[i].possibleBranch;
1529  delete [] choice;
1530  delete [] saveLower;
1531  delete [] saveUpper;
1532 
1533  // restore solution
1534  solver->setColSolution(saveSolution);
1535  delete [] saveSolution;
1536  return anyAction;
1537}
1538
1539
1540CbcNode::CbcNode(const CbcNode & rhs) 
1541{ 
1542#ifdef CHECK_NODE
1543  printf("CbcNode %x Constructor from rhs %x\n",this,&rhs);
1544#endif
1545  if (rhs.nodeInfo_)
1546    nodeInfo_ = rhs.nodeInfo_->clone();
1547  else
1548    nodeInfo_=NULL;
1549  objectiveValue_=rhs.objectiveValue_;
1550  guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
1551  if (rhs.branch_)
1552    branch_=rhs.branch_->clone();
1553  else
1554    branch_=NULL;
1555  depth_ = rhs.depth_;
1556  numberUnsatisfied_ = rhs.numberUnsatisfied_;
1557  nodeNumber_=rhs.nodeNumber_;
1558}
1559
1560CbcNode &
1561CbcNode::operator=(const CbcNode & rhs)
1562{
1563  if (this != &rhs) {
1564    delete nodeInfo_;
1565    if (nodeInfo_)
1566      nodeInfo_ = rhs.nodeInfo_->clone();
1567    else
1568      nodeInfo_ = NULL;
1569    objectiveValue_=rhs.objectiveValue_;
1570    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
1571    if (rhs.branch_)
1572      branch_=rhs.branch_->clone();
1573    else
1574      branch_=NULL,
1575    depth_ = rhs.depth_;
1576    numberUnsatisfied_ = rhs.numberUnsatisfied_;
1577    nodeNumber_=rhs.nodeNumber_;
1578  }
1579  return *this;
1580}
1581
1582
1583CbcNode::~CbcNode ()
1584{
1585#ifdef CHECK_NODE
1586  if (nodeInfo_) 
1587    printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
1588         this,nodeInfo_,nodeInfo_->numberPointingToThis());
1589  else
1590    printf("CbcNode %x Destructor nodeInfo %x (?)\n",
1591         this,nodeInfo_);
1592#endif
1593  if (nodeInfo_) {
1594    nodeInfo_->nullOwner();
1595    int numberToDelete=nodeInfo_->numberBranchesLeft();
1596    //    CbcNodeInfo * parent = nodeInfo_->parent();
1597    if (nodeInfo_->decrement(numberToDelete)==0) {
1598      delete nodeInfo_;
1599    } else {
1600      //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,parent);
1601      // anyway decrement parent
1602      //if (parent)
1603      ///parent->decrement(1);
1604    }
1605  }
1606  delete branch_;
1607}
1608// Decrement  active cut counts
1609void 
1610CbcNode::decrementCuts(int change)
1611{
1612  if(nodeInfo_) {
1613    nodeInfo_->decrementCuts(change);
1614  }
1615}
1616void 
1617CbcNode::decrementParentCuts(int change)
1618{
1619  if(nodeInfo_) {
1620    nodeInfo_->decrementParentCuts(change);
1621  }
1622}
1623
1624/*
1625  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
1626  in the attached nodeInfo_.
1627*/
1628void
1629CbcNode::initializeInfo()
1630{
1631  assert(nodeInfo_ && branch_) ;
1632  nodeInfo_->initializeInfo(branch_->numberBranches());
1633}
1634// Nulls out node info
1635void 
1636CbcNode::nullNodeInfo()
1637{
1638  nodeInfo_=NULL;
1639}
1640
1641int
1642CbcNode::branch()
1643{
1644  double changeInGuessed=branch_->branch(true);
1645  guessedObjectiveValue_+= changeInGuessed;
1646  return nodeInfo_->branchedOn();
1647}
Note: See TracBrowser for help on using the repository browser.