source: trunk/CbcNode.cpp @ 129

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

major changes

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