Main Page | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | Related Pages

SbbNode.cpp

00001 // Copyright (C) 2002, International Business Machines
00002 // Corporation and others.  All Rights Reserved.
00003 #if defined(_MSC_VER)
00004 // Turn off compiler warning about long names
00005 #  pragma warning(disable:4786)
00006 #endif
00007 #include <string>
00008 //#define SBB_DEBUG 1
00009 //#define CHECK_CUT_COUNTS
00010 #include <cassert>
00011 #include <cfloat>
00012 #define CUTS
00013 #include "OsiSolverInterface.hpp"
00014 #include "CoinWarmStartBasis.hpp"
00015 #include "SbbModel.hpp"
00016 #include "SbbNode.hpp"
00017 #include "SbbBranchActual.hpp"
00018 #include "OsiRowCut.hpp"
00019 #include "OsiRowCutDebugger.hpp"
00020 #include "OsiCuts.hpp"
00021 #include "SbbCountRowCut.hpp"
00022 #include "SbbMessage.hpp"
00023 #include "OsiClpSolverInterface.hpp"
00024 using namespace std;
00025 #include "CglCutGenerator.hpp"
00026 // Default Constructor 
00027 SbbNodeInfo::SbbNodeInfo ()
00028   :
00029   numberPointingToThis_(0),
00030   parent_(NULL),
00031   owner_(NULL),
00032   numberCuts_(0),
00033   cuts_(NULL),
00034   numberRows_(0),
00035   numberBranchesLeft_(0)
00036 {
00037 #ifdef CHECK_NODE
00038   printf("SbbNodeInfo %x Constructor\n",this);
00039 #endif
00040 }
00041 // Constructor given parent
00042 SbbNodeInfo::SbbNodeInfo (SbbNodeInfo * parent)
00043   :
00044   numberPointingToThis_(2),
00045   parent_(parent),
00046   owner_(NULL),
00047   numberCuts_(0),
00048   cuts_(NULL),
00049   numberRows_(0),
00050   numberBranchesLeft_(2)
00051 {
00052 #ifdef CHECK_NODE
00053   printf("SbbNodeInfo %x Constructor from parent %x\n",this,parent_);
00054 #endif
00055   if (parent_) {
00056     numberRows_ = parent_->numberRows_+parent_->numberCuts_;
00057   }
00058 }
00059 // Constructor given parent and owner
00060 SbbNodeInfo::SbbNodeInfo (SbbNodeInfo * parent, SbbNode * owner)
00061   :
00062   numberPointingToThis_(2),
00063   parent_(parent),
00064   owner_(owner),
00065   numberCuts_(0),
00066   cuts_(NULL),
00067   numberRows_(0),
00068   numberBranchesLeft_(2)
00069 {
00070 #ifdef CHECK_NODE
00071   printf("SbbNodeInfo %x Constructor from parent %x\n",this,parent_);
00072 #endif
00073   if (parent_) {
00074     numberRows_ = parent_->numberRows_+parent_->numberCuts_;
00075   }
00076 }
00077 
00083 SbbNodeInfo::~SbbNodeInfo()
00084 {
00085 #ifdef CHECK_NODE
00086   printf("SbbNodeInfo %x Destructor parent %x\n",this,parent_);
00087 #endif
00088 
00089   assert(!numberPointingToThis_);
00090 
00091   delete [] cuts_;
00092   if (owner_) 
00093     owner_->nullNodeInfo();
00094   if (parent_) {
00095     int numberLinks = parent_->decrement();
00096     if (!numberLinks) delete parent_;
00097   }
00098 }
00099 
00100 
00101 //#define ALLCUTS
00102 void
00103 SbbNodeInfo::decrementCuts(int change)
00104 {
00105   int i;
00106   // get rid of all remaining if negative
00107   int changeThis;
00108   if (change<0)
00109     changeThis = numberBranchesLeft_;
00110   else
00111     changeThis = change;
00112  // decrement cut counts
00113   for (i=0;i<numberCuts_;i++) {
00114     if (cuts_[i]) {
00115       int number = cuts_[i]->decrement(changeThis);
00116       if (!number) {
00117         //printf("info %x del cut %d %x\n",this,i,cuts_[i]);
00118         delete cuts_[i];
00119         cuts_[i]=NULL;
00120       }
00121     }
00122   }
00123 }
00124 void
00125 SbbNodeInfo::decrementParentCuts(int change)
00126 {
00127   if (parent_) {
00128     // get rid of all remaining if negative
00129     int changeThis;
00130     if (change<0)
00131       changeThis = numberBranchesLeft_;
00132     else
00133       changeThis = change;
00134     int i;
00135     // Get over-estimate of space needed for basis
00136     CoinWarmStartBasis dummy;
00137     dummy.setSize(0,numberRows_+numberCuts_);
00138     buildRowBasis(dummy);
00139     /* everything is zero (i.e. free) so we can use to see
00140        if latest basis */
00141     SbbNodeInfo * thisInfo = parent_;
00142     while (thisInfo) 
00143       thisInfo = thisInfo->buildRowBasis(dummy);
00144     // decrement cut counts
00145     thisInfo = parent_;
00146     int numberRows=numberRows_;
00147     while (thisInfo) {
00148       for (i=thisInfo->numberCuts_-1;i>=0;i--) {
00149         CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
00150 #ifdef ALLCUTS
00151         status = CoinWarmStartBasis::isFree;
00152 #endif
00153         if (thisInfo->cuts_[i]) {
00154           int number=1;
00155           if (status!=CoinWarmStartBasis::basic) {
00156             // tight - drop 1 or 2
00157             if (change<0)
00158               number = thisInfo->cuts_[i]->decrement(changeThis);
00159             else
00160               number = thisInfo->cuts_[i]->decrement(change);
00161           }
00162           if (!number) {
00163             delete thisInfo->cuts_[i];
00164             thisInfo->cuts_[i]=NULL;
00165           }
00166         }
00167       }
00168       thisInfo = thisInfo->parent_;
00169     }
00170   }
00171 }
00172 
00173 void
00174 SbbNodeInfo::incrementParentCuts(int change)
00175 {
00176   if (parent_) {
00177     int i;
00178     // Get over-estimate of space needed for basis
00179     CoinWarmStartBasis dummy;
00180     dummy.setSize(0,numberRows_+numberCuts_);
00181     /* everything is zero (i.e. free) so we can use to see
00182        if latest basis */
00183     buildRowBasis(dummy);
00184     SbbNodeInfo * thisInfo = parent_;
00185     while (thisInfo) 
00186       thisInfo = thisInfo->buildRowBasis(dummy);
00187     // increment cut counts
00188     thisInfo = parent_;
00189     int numberRows=numberRows_;
00190     while (thisInfo) {
00191       for (i=thisInfo->numberCuts_-1;i>=0;i--) {
00192         CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
00193 #ifdef ALLCUTS
00194         status = CoinWarmStartBasis::isFree;
00195 #endif
00196         if (thisInfo->cuts_[i]&&status!=CoinWarmStartBasis::basic) {
00197           thisInfo->cuts_[i]->increment(change);
00198         }
00199       }
00200       thisInfo = thisInfo->parent_;
00201     }
00202   }
00203 }
00204 
00205 /*
00206   Append cuts to the cuts_ array in a nodeInfo. The initial reference count
00207   is set to numberToBranchOn, which will normally be the number of arms
00208   defined for the SbbBranchingObject attached to the SbbNode that owns this
00209   SbbNodeInfo.
00210 */
00211 void
00212 SbbNodeInfo::addCuts (OsiCuts & cuts, int numberToBranchOn,
00213                       int * whichGenerator)
00214 {
00215   int numberCuts = cuts.sizeRowCuts();
00216   if (numberCuts) {
00217     int i;
00218     if (!numberCuts_) {
00219       cuts_ = new SbbCountRowCut * [numberCuts];
00220     } else {
00221       SbbCountRowCut ** temp = new SbbCountRowCut * [numberCuts+numberCuts_];
00222       memcpy(temp,cuts_,numberCuts_*sizeof(SbbCountRowCut *));
00223       delete [] cuts_;
00224       cuts_ = temp;
00225     }
00226     for (i=0;i<numberCuts;i++) {
00227       SbbCountRowCut * thisCut = new SbbCountRowCut(*cuts.rowCutPtr(i),
00228                                                     this,numberCuts_);
00229       thisCut->increment(numberToBranchOn); 
00230       cuts_[numberCuts_++] = thisCut;
00231 #ifdef SBB_DEBUG
00232       int n=thisCut->row().getNumElements();
00233 #if SBB_DEBUG>1
00234       printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
00235              thisCut->ub());
00236 #endif
00237       int j;
00238 #if SBB_DEBUG>1
00239       const int * index = thisCut->row().getIndices();
00240 #endif
00241       const double * element = thisCut->row().getElements();
00242       for (j=0;j<n;j++) {
00243 #if SBB_DEBUG>1
00244         printf(" (%d,%g)",index[j],element[j]);
00245 #endif
00246         assert(fabs(element[j])>1.00e-12);
00247       }
00248       printf("\n");
00249 #endif
00250     }
00251   }
00252 }
00253 
00254 void
00255 SbbNodeInfo::addCuts(int numberCuts, SbbCountRowCut ** cut, 
00256                      int numberToBranchOn)
00257 {
00258   if (numberCuts) {
00259     int i;
00260     if (!numberCuts_) {
00261       cuts_ = new SbbCountRowCut * [numberCuts];
00262     } else {
00263       SbbCountRowCut ** temp = new SbbCountRowCut * [numberCuts+numberCuts_];
00264       memcpy(temp,cuts_,numberCuts_*sizeof(SbbCountRowCut *));
00265       delete [] cuts_;
00266       cuts_ = temp;
00267     }
00268     for (i=0;i<numberCuts;i++) {
00269       SbbCountRowCut * thisCut = cut[i];
00270       thisCut->setInfo(this,numberCuts_);
00271       //printf("info %x cut %d %x\n",this,i,thisCut);
00272       thisCut->increment(numberToBranchOn); 
00273       cuts_[numberCuts_++] = thisCut;
00274 #ifdef SBB_DEBUG
00275       int n=thisCut->row().getNumElements();
00276 #if SBB_DEBUG>1
00277       printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
00278              thisCut->ub());
00279 #endif
00280       int j;
00281 #if SBB_DEBUG>1
00282       const int * index = thisCut->row().getIndices();
00283 #endif
00284       const double * element = thisCut->row().getElements();
00285       for (j=0;j<n;j++) {
00286 #if SBB_DEBUG>1
00287         printf(" (%d,%g)",index[j],element[j]);
00288 #endif
00289         assert(fabs(element[j])>1.00e-12);
00290       }
00291       printf("\n");
00292 #endif
00293     }
00294   }
00295 }
00296 
00297 // delete cuts
00298 void
00299 SbbNodeInfo::deleteCuts(int numberToDelete, SbbCountRowCut ** cuts)
00300 {
00301   int i;
00302   int j;
00303   int last=-1;
00304   for (i=0;i<numberToDelete;i++) {
00305     SbbCountRowCut * next = cuts[i];
00306     for (j=last+1;j<numberCuts_;j++) {
00307       if (next==cuts_[j])
00308         break;
00309     }
00310     if (j==numberCuts_) {
00311       // start from beginning
00312       for (j=0;j<last;j++) {
00313         if (next==cuts_[j])
00314           break;
00315       }
00316       assert(j<last);
00317     }
00318     last=j;
00319     int number = cuts_[j]->decrement();
00320     if (!number)
00321       delete cuts_[j];
00322     cuts_[j]=NULL;
00323   }
00324   j=0;
00325   for (i=0;i<numberCuts_;i++) {
00326     if (cuts_[i])
00327       cuts_[j++]=cuts_[i];
00328   }
00329   numberCuts_ = j;
00330 }
00331 
00332 // delete cuts
00333 void
00334 SbbNodeInfo::deleteCuts(int numberToDelete, int * which)
00335 {
00336   int i;
00337   for (i=0;i<numberToDelete;i++) {
00338     int iCut=which[i];
00339     int number = cuts_[iCut]->decrement();
00340     if (!number)
00341       delete cuts_[iCut];
00342     cuts_[iCut]=NULL;
00343   }
00344   int j=0;
00345   for (i=0;i<numberCuts_;i++) {
00346     if (cuts_[i])
00347       cuts_[j++]=cuts_[i];
00348   }
00349   numberCuts_ = j;
00350 }
00351 
00352 // Really delete a cut
00353 void 
00354 SbbNodeInfo::deleteCut(int whichOne)
00355 {
00356   assert(whichOne<numberCuts_);
00357   cuts_[whichOne]=NULL;
00358 }
00359 
00360 SbbFullNodeInfo::SbbFullNodeInfo() :
00361   SbbNodeInfo(),
00362   basis_(),
00363   numberIntegers_(0),
00364   lower_(NULL),
00365   upper_(NULL)
00366 {
00367 }
00368 SbbFullNodeInfo::SbbFullNodeInfo(SbbModel * model,
00369                                  int numberRowsAtContinuous) :
00370   SbbNodeInfo()
00371 {
00372   OsiSolverInterface * solver = model->solver();
00373   numberRows_ = numberRowsAtContinuous;
00374   numberIntegers_ = model->numberIntegers();
00375   int numberColumns = model->getNumCols();
00376   lower_ = new double [numberColumns];
00377   upper_ = new double [numberColumns];
00378   const double * lower = solver->getColLower();
00379   const double * upper = solver->getColUpper();
00380   int i;
00381 
00382   for (i=0;i<numberColumns;i++) {
00383     lower_[i]=lower[i];
00384     upper_[i]=upper[i];
00385   }
00386 
00387   basis_ =  dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
00388 }
00389 
00390 SbbFullNodeInfo::SbbFullNodeInfo(const SbbFullNodeInfo & rhs) :
00391   SbbNodeInfo(rhs)
00392 {  
00393   basis_= dynamic_cast<CoinWarmStartBasis *>(rhs.basis_->clone()) ;
00394   numberIntegers_=rhs.numberIntegers_;
00395   lower_=NULL;
00396   upper_=NULL;
00397   if (rhs.lower_!=NULL) {
00398     int numberColumns = basis_->getNumStructural();
00399     lower_ = new double [numberColumns];
00400     upper_ = new double [numberColumns];
00401     assert (upper_!=NULL);
00402     memcpy(lower_,rhs.lower_,numberColumns*sizeof(double));
00403     memcpy(upper_,rhs.upper_,numberColumns*sizeof(double));
00404   }
00405 }
00406 
00407 SbbNodeInfo * 
00408 SbbFullNodeInfo::clone() const
00409 { 
00410   return (new SbbFullNodeInfo(*this));
00411 }
00412 
00413 SbbFullNodeInfo::~SbbFullNodeInfo ()
00414 {
00415   delete basis_ ;
00416   delete [] lower_;
00417   delete [] upper_;
00418 }
00419 
00420 /*
00421    The basis supplied as a parameter is deleted and replaced with a new basis
00422    appropriate for the node, and lower and upper bounds on variables are
00423    reset according to the stored bounds arrays. Any cuts associated with this
00424    node are added to the list in addCuts, but not actually added to the
00425    constraint system in the model.
00426 
00427    Why pass in a basis at all? The short answer is ``We need the parameter to
00428    pass out a basis, so might as well use it to pass in the size.''
00429    
00430    A longer answer is that in practice we take a memory allocation hit up in
00431    addCuts1 (the only place applyToModel is called) when we setSize() the
00432    basis that's passed in. It's immediately tossed here in favour of a clone
00433    of the basis attached to this nodeInfo. This can probably be fixed, given
00434    a bit of thought.
00435 */
00436 
00437 void SbbFullNodeInfo::applyToModel (SbbModel *model,
00438                                     CoinWarmStartBasis *&basis,
00439                                     SbbCountRowCut **addCuts,
00440                                     int &currentNumberCuts) const 
00441 
00442 { OsiSolverInterface *solver = model->solver() ;
00443 
00444   // branch - do bounds
00445   int i;
00446   int numberColumns = model->getNumCols();
00447   for (i=0;i<numberColumns;i++) {
00448     solver->setColBounds(i,lower_[i],upper_[i]);
00449   }
00450   // move basis - but make sure size stays
00451   int numberRows = basis->getNumArtificial();
00452   delete basis ;
00453   basis = dynamic_cast<CoinWarmStartBasis *>(basis_->clone()) ;
00454   basis->resize(numberRows,numberColumns);
00455   for (i=0;i<numberCuts_;i++) 
00456     addCuts[currentNumberCuts+i]= cuts_[i];
00457   currentNumberCuts += numberCuts_;
00458   assert(!parent_);
00459   return ;
00460 }
00461 
00462 /* Builds up row basis backwards (until original model).
00463    Returns NULL or previous one to apply .
00464    Depends on Free being 0 and impossible for cuts
00465 */
00466 SbbNodeInfo * 
00467 SbbFullNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
00468 {
00469   const unsigned int * saved = 
00470     (const unsigned int *) basis_->getArtificialStatus();
00471   unsigned int * now = 
00472     (unsigned int *) basis.getArtificialStatus();
00473   int number=basis_->getNumArtificial()>>4;;
00474   int i;
00475   for (i=0;i<number;i++) { 
00476     if (!now[i])
00477       now[i] = saved[i];
00478   }
00479   return NULL;
00480 }
00481 
00482 // Default constructor
00483 SbbPartialNodeInfo::SbbPartialNodeInfo()
00484 
00485   : SbbNodeInfo(),
00486     basisDiff_(NULL),
00487     variables_(NULL),
00488     newBounds_(NULL),
00489     numberChangedBounds_(0)
00490 
00491 { /* this space intentionally left blank */ }
00492 
00493 // Constructor from current state 
00494 SbbPartialNodeInfo::SbbPartialNodeInfo (SbbNodeInfo *parent, SbbNode *owner,
00495                                         int numberChangedBounds,
00496                                         const int *variables,
00497                                         const double *boundChanges,
00498                                         const CoinWarmStartDiff *basisDiff)
00499  : SbbNodeInfo(parent,owner)
00500 {
00501   basisDiff_ = basisDiff->clone() ;
00502 
00503   numberChangedBounds_ = numberChangedBounds;
00504   variables_ = new int [numberChangedBounds_];
00505   newBounds_ = new double [numberChangedBounds_];
00506 
00507   int i ;
00508   for (i=0;i<numberChangedBounds_;i++) {
00509     variables_[i]=variables[i];
00510     newBounds_[i]=boundChanges[i];
00511   }
00512 }
00513 
00514 SbbPartialNodeInfo::SbbPartialNodeInfo (const SbbPartialNodeInfo & rhs)
00515 
00516   : SbbNodeInfo(rhs.parent_)
00517 
00518 { basisDiff_ = rhs.basisDiff_->clone() ;
00519 
00520   numberChangedBounds_ = rhs.numberChangedBounds_;
00521   variables_ = new int [numberChangedBounds_];
00522   newBounds_ = new double [numberChangedBounds_];
00523 
00524   int i ;
00525   for (i=0;i<numberChangedBounds_;i++) {
00526     variables_[i]=rhs.variables_[i];
00527     newBounds_[i]=rhs.newBounds_[i];
00528   }
00529 }
00530 
00531 SbbNodeInfo * 
00532 SbbPartialNodeInfo::clone() const
00533 { 
00534   return (new SbbPartialNodeInfo(*this));
00535 }
00536 
00537 
00538 SbbPartialNodeInfo::~SbbPartialNodeInfo ()
00539 {
00540   delete basisDiff_ ;
00541   delete [] variables_;
00542   delete [] newBounds_;
00543 }
00544 
00545 
00552 void SbbPartialNodeInfo::applyToModel (SbbModel *model,
00553                                        CoinWarmStartBasis *&basis,
00554                                        SbbCountRowCut **addCuts,
00555                                        int &currentNumberCuts) const 
00556 
00557 { OsiSolverInterface *solver = model->solver();
00558 
00559   basis->applyDiff(basisDiff_) ;
00560 
00561   // branch - do bounds
00562   int i;
00563   for (i=0;i<numberChangedBounds_;i++) {
00564     int variable = variables_[i];
00565     if ((variable&0x80000000)==0) {
00566       // lower bound changing
00567       solver->setColLower(variable,newBounds_[i]);
00568     } else {
00569       // upper bound changing
00570       solver->setColUpper(variable&0x7fffffff,newBounds_[i]);
00571     }
00572   }
00573   for (i=0;i<numberCuts_;i++) 
00574     addCuts[currentNumberCuts+i]= cuts_[i];
00575   currentNumberCuts += numberCuts_;
00576   return ;
00577 }
00578 
00579 /* Builds up row basis backwards (until original model).
00580    Returns NULL or previous one to apply .
00581    Depends on Free being 0 and impossible for cuts
00582 */
00583 
00584 SbbNodeInfo * 
00585 SbbPartialNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const 
00586 
00587 { basis.applyDiff(basisDiff_) ;
00588 
00589   return parent_ ; }
00590 
00591 
00592 SbbNode::SbbNode() :
00593   nodeInfo_(NULL),
00594   objectiveValue_(1.0e100),
00595   guessedObjectiveValue_(1.0e100),
00596   branch_(NULL),
00597   depth_(-1),
00598   numberUnsatisfied_(0)
00599 {
00600 #ifdef CHECK_NODE
00601   printf("SbbNode %x Constructor\n",this);
00602 #endif
00603 }
00604 
00605 SbbNode::SbbNode(SbbModel * model,
00606                  SbbNode * lastNode) :
00607   nodeInfo_(NULL),
00608   objectiveValue_(1.0e100),
00609   guessedObjectiveValue_(1.0e100),
00610   branch_(NULL),
00611   depth_(-1),
00612   numberUnsatisfied_(0)
00613 {
00614 #ifdef CHECK_NODE
00615   printf("SbbNode %x Constructor from model\n",this);
00616 #endif
00617   OsiSolverInterface * solver = model->solver();
00618   objectiveValue_ = solver->getObjSense()*solver->getObjValue();
00619 
00620   if (lastNode)
00621     lastNode->nodeInfo_->increment();
00622 }
00623 
00624 
00625 void
00626 SbbNode::createInfo (SbbModel *model,
00627                      SbbNode *lastNode,
00628                      const CoinWarmStartBasis *lastws,
00629                      const double *lastLower, const double *lastUpper,
00630                      int numberOldActiveCuts,int numberNewCuts)
00631 { OsiSolverInterface * solver = model->solver();
00632 /*
00633   The root --- no parent. Create full basis and bounds information.
00634 */
00635   if (!lastNode)
00636   { nodeInfo_=new SbbFullNodeInfo(model,solver->getNumRows()); }
00637 /*
00638   Not the root. Create an edit from the parent's basis & bound information.
00639   This is not quite as straightforward as it seems. We need to reintroduce
00640   cuts we may have dropped out of the basis, in the correct position, because
00641   this whole process is strictly positional. Start by grabbing the current
00642   basis.
00643 */
00644   else
00645   { const CoinWarmStartBasis* ws =
00646       dynamic_cast<const CoinWarmStartBasis*>(solver->getWarmStart());
00647     assert(ws!=NULL); // make sure not volume
00648     //int numberArtificials = lastws->getNumArtificial();
00649     int numberColumns = solver->getNumCols();
00650     
00651     const double * lower = solver->getColLower();
00652     const double * upper = solver->getColUpper();
00653 
00654     int i;
00655 /*
00656   Create a clone and resize it to hold all the structural constraints, plus
00657   all the cuts: old cuts, both active and inactive (currentNumberCuts), and
00658   new cuts (numberNewCuts).
00659 
00660   TODO: You'd think that the set of constraints (logicals) in the expanded
00661         basis should match the set represented in lastws. At least, that's
00662         what I thought. But it turns out that lastws is actually the stripped
00663         basis produced at the end of addCuts(), rather than the raw basis
00664         handed back by addCuts1(). The expanded basis here is equivalent to
00665         the raw basis of addCuts1(). I said ``whoa, that's not good'' and went
00666         back to John's code to see where I'd gone wrong. And discovered the
00667         same `error' in his code.
00668 
00669         After a bit of thought, my conclusion is that correctness is not
00670         affected by whether lastws is the stripped or raw basis. The diffs
00671         have no semantics --- just a set of changes that need to be made
00672         to convert lastws into expanded. I think the only effect is that we
00673         store a lot more diffs (everything in expanded that's not covered by
00674         the stripped basis). But I need to give this more thought. There
00675         may well be some subtle error cases.
00676 
00677         In the mean time, I've twiddled addCuts() to set lastws to the raw
00678         basis. Makes me less nervous to compare apples to apples.
00679 */
00680     CoinWarmStartBasis *expanded = 
00681       dynamic_cast<CoinWarmStartBasis *>(ws->clone()) ;
00682     int numberRowsAtContinuous = model->numberRowsAtContinuous();
00683     int iFull = numberRowsAtContinuous+model->currentNumberCuts()+
00684       numberNewCuts;
00685     //int numberArtificialsNow = iFull;
00686     //int maxBasisLength = ((iFull+15)>>4)+((numberColumns+15)>>4);
00687     //printf("l %d full %d\n",maxBasisLength,iFull);
00688     expanded->resize(iFull,numberColumns);
00689 #ifdef FULL_DEBUG
00690     printf("Before expansion: orig %d, old %d, new %d, current %d\n",
00691            numberRowsAtContinuous,numberOldActiveCuts,numberNewCuts,
00692            model->currentNumberCuts()) ;
00693     ws->print();
00694 #endif
00695 /*
00696   Now fill in the expanded basis. Anything indices beyond nPartial must
00697   be cuts created while processing this node --- they can be copied directly
00698   into the expanded basis. From nPartial down, pull the status of active cuts
00699   from ws, interleaving with a B entry for the deactivated (loose) cuts.
00700 */
00701     int numberDropped = model->currentNumberCuts()-numberOldActiveCuts;
00702     int iCompact=iFull-numberDropped;
00703     SbbCountRowCut ** cut = model->addedCuts();
00704     int nPartial = model->currentNumberCuts()+numberRowsAtContinuous;
00705     iFull--;
00706     for (;iFull>=nPartial;iFull--) {
00707       CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
00708       assert (status != CoinWarmStartBasis::basic);
00709       expanded->setArtifStatus(iFull,status);
00710     }
00711     for (;iFull>=numberRowsAtContinuous;iFull--) {
00712       if (cut[iFull-numberRowsAtContinuous]) {
00713         CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
00714         assert (status != CoinWarmStartBasis::basic);
00715         expanded->setArtifStatus(iFull,status);
00716       } else {
00717         expanded->setArtifStatus(iFull,CoinWarmStartBasis::basic);
00718       }
00719     }
00720 #ifdef FULL_DEBUG
00721     printf("Expanded basis\n");
00722     expanded->print() ;
00723     printf("Diffing against\n") ;
00724     lastws->print() ;
00725 #endif    
00726 /*
00727   Now that we have two bases in proper positional correspondence, creating
00728   the actual diff is dead easy.
00729 */
00730 
00731     CoinWarmStartDiff *basisDiff = expanded->generateDiff(lastws) ;
00732 /*
00733   Diff the bound vectors. It's assumed the number of structural variables is
00734   not changing. Assuming that branching objects all involve integer variables,
00735   we should see at least one bound change as a consequence of processing this
00736   subproblem. Different types of branching objects could break this assertion.
00737 */
00738     double *boundChanges = new double [2*numberColumns] ;
00739     int *variables = new int [2*numberColumns] ;
00740     int numberChangedBounds=0;
00741     for (i=0;i<numberColumns;i++) {
00742       if (lower[i]!=lastLower[i]) {
00743         variables[numberChangedBounds]=i;
00744         boundChanges[numberChangedBounds++]=lower[i];
00745       }
00746       if (upper[i]!=lastUpper[i]) {
00747         variables[numberChangedBounds]=i|0x80000000;
00748         boundChanges[numberChangedBounds++]=upper[i];
00749       }
00750 #ifdef SBB_DEBUG
00751       if (lower[i]!=lastLower[i])
00752         printf("lower on %d changed from %g to %g\n",
00753                i,lastLower[i],lower[i]);
00754       if (upper[i]!=lastUpper[i])
00755         printf("upper on %d changed from %g to %g\n",
00756                i,lastUpper[i],upper[i]);
00757 #endif
00758     }
00759 #ifdef SBB_DEBUG
00760     printf("%d changed bounds\n",numberChangedBounds) ;
00761 #endif
00762     assert (numberChangedBounds);
00763 /*
00764   Hand the lot over to the SbbPartialNodeInfo constructor, then clean up and
00765   return.
00766 */
00767     nodeInfo_ =
00768       new SbbPartialNodeInfo(lastNode->nodeInfo_,this,numberChangedBounds,
00769                              variables,boundChanges,basisDiff) ;
00770     delete basisDiff ;
00771     delete [] boundChanges;
00772     delete [] variables;
00773     delete expanded ;
00774     delete ws;
00775   }
00776 }
00777 
00778 /*
00779   The routine scans through the object list of the model looking for objects
00780   that indicate infeasibility. It tests each object using strong branching
00781   and selects the one with the least objective degradation.  A corresponding
00782   branching object is left attached to lastNode.
00783 
00784   If strong branching is disabled, a candidate object is chosen essentially
00785   at random (whatever object ends up in pos'n 0 of the candidate array).
00786 
00787   If a branching candidate is found to be monotone, bounds are set to fix the
00788   variable and the routine immediately returns (the caller is expected to
00789   reoptimize).
00790 
00791   If a branching candidate is found to result in infeasibility in both
00792   directions, the routine immediately returns an indication of infeasibility.
00793 
00794   Returns:  0   both branch directions are feasible
00795            -1   branching variable is monotone
00796            -2   infeasible
00797 
00798   Original comments:
00799     Here could go cuts etc etc
00800     For now just fix on objective from strong branching.
00801 */
00802 
00803 int SbbNode::chooseBranch (SbbModel *model, SbbNode *lastNode)
00804 
00805 { if (lastNode)
00806     depth_ = lastNode->depth_+1;
00807   else
00808     depth_ = 0;
00809   delete branch_;
00810   branch_=NULL;
00811   OsiSolverInterface * solver = model->solver();
00812   double saveObjectiveValue = solver->getObjValue();
00813   double objectiveValue = solver->getObjSense()*saveObjectiveValue;
00814   const double * lower = solver->getColLower();
00815   const double * upper = solver->getColUpper();
00816 
00817   int anyAction=0;
00818   double integerTolerance = 
00819     model->getDblParam(SbbModel::SbbIntegerTolerance);
00820   int i;
00821   bool beforeSolution = model->getSolutionCount()==0;
00822   int numberStrong=model->numberStrong();
00823   int numberObjects = model->numberObjects();
00824   int maximumStrong = max(min(model->numberStrong(),numberObjects),1);
00825   int numberColumns = model->getNumCols();
00826   double * saveUpper = new double[numberColumns];
00827   double * saveLower = new double[numberColumns];
00828 
00829   // Save solution in case heuristics need good solution later
00830 
00831   double * saveSolution = new double[numberColumns];
00832   memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
00833 
00834 /*
00835   Get a branching decision object. Use the default decision criteria unless
00836   the user has loaded a decision method into the model.
00837 */
00838   SbbBranchDecision *decision = model->branchingMethod();
00839   if (!decision)
00840     decision = new SbbBranchDefaultDecision();
00841 
00842   typedef struct {
00843     SbbBranchingObject * possibleBranch; // what a branch would do
00844     double upMovement; // cost going up (and initial away from feasible)
00845     double downMovement; // cost going down
00846     int numIntInfeasUp ; // without odd ones
00847     int numObjInfeasUp ; // just odd ones
00848     bool finishedUp; // true if solver finished
00849     int numItersUp ; // number of iterations in solver
00850     int numIntInfeasDown ; // without odd ones
00851     int numObjInfeasDown ; // just odd ones
00852     bool finishedDown; // true if solver finished
00853     int numItersDown; // number of iterations in solver
00854     int objectNumber; // Which object it is
00855   } Strong;
00856   Strong * choice = new Strong[maximumStrong];
00857   for (i=0;i<numberColumns;i++) {
00858     saveLower[i] = lower[i];
00859     saveUpper[i] = upper[i];
00860   }
00861 
00862   // compute current state
00863   int numberIntegerInfeasibilities; // without odd ones
00864   int numberObjectInfeasibilities; // just odd ones
00865   model->feasibleSolution(
00866                            numberIntegerInfeasibilities,
00867                            numberObjectInfeasibilities);
00868 /*
00869   Original comment ``Put in coding for forcePriority for hot starts.'' Looks
00870   like a `to do' to me, given the assert.   -- lh, 2003.08.
00871 */
00872   assert(model->getForcePriority()<0);
00873 
00874   double saveOtherWay=0.0; // just for non-strong branching
00875   numberUnsatisfied_ = 0;
00876   int bestPriority=INT_MAX;
00877 /*
00878   Scan for branching objects that indicate infeasibility. Choose the best
00879   maximumStrong candidates, using priority as the first criteria, then
00880   integer infeasibility.
00881 
00882   The algorithm is to fill the choice array with a set of good candidates (by
00883   infeasibility) with priority bestPriority.  Finding a candidate with
00884   priority better (less) than bestPriority flushes the choice array. (This
00885   serves as initialization when the first candidate is found.)
00886 
00887   A new candidate is added to choices only if its infeasibility exceeds the
00888   current max infeasibility (mostAway). When a candidate is added, it
00889   replaces the candidate with the smallest infeasibility (tracked by
00890   iSmallest).
00891 */
00892   int iSmallest ;
00893   double mostAway ;
00894   for (i = 0 ; i < maximumStrong ; i++) choice[i].possibleBranch = NULL ;
00895   numberStrong=0;
00896   for (i=0;i<numberObjects;i++) {
00897     const SbbObject * object = model->object(i);
00898     int preferredWay;
00899     double otherWay;
00900     double infeasibility = object->infeasibility(preferredWay,otherWay);
00901     if (infeasibility>integerTolerance) {
00902       numberUnsatisfied_++;
00903       int priorityLevel = model->priority(i);
00904       // Better priority? Flush choices.
00905       if (priorityLevel<bestPriority) {
00906         int j;
00907         iSmallest=0;
00908         for (j=0;j<maximumStrong;j++) {
00909           choice[j].upMovement=0.0;
00910           delete choice[j].possibleBranch;
00911           choice[j].possibleBranch=NULL;
00912         }
00913         bestPriority = priorityLevel;
00914         mostAway=integerTolerance;
00915         numberStrong=0;
00916       } else if (priorityLevel>bestPriority) {
00917         continue;
00918       }
00919       // Check for suitability based on infeasibility.
00920       if (infeasibility>mostAway) {
00921         //add to list
00922         choice[iSmallest].upMovement=infeasibility;
00923         choice[iSmallest].downMovement=otherWay;
00924         saveOtherWay=otherWay;
00925         delete choice[iSmallest].possibleBranch;
00926         choice[iSmallest].possibleBranch=object->createBranch(preferredWay);
00927         numberStrong = max(numberStrong,iSmallest+1);
00928         // Save which object it was
00929         choice[iSmallest].objectNumber=i;
00930         int j;
00931         iSmallest=-1;
00932         mostAway = 1.0e50;
00933         for (j=0;j<maximumStrong;j++) {
00934           if (choice[j].upMovement<mostAway) {
00935             mostAway=choice[j].upMovement;
00936             iSmallest=j;
00937           }
00938         }
00939       }
00940     }
00941   }
00942   /* Some solvers can do the strong branching calculations faster if
00943      they do them all at once.  At present only Clp does for ordinary
00944      integers but I think this coding would be easy to modify 
00945   */
00946   bool allNormal=true; // to say if we can do fast strong branching
00947   for (i=0;i<numberStrong;i++) {
00948     choice[i].numIntInfeasUp = numberUnsatisfied_;
00949     choice[i].numIntInfeasDown = numberUnsatisfied_;
00950     if (!dynamic_cast <const SbbSimpleInteger *> (model->object(choice[i].objectNumber)))
00951       allNormal=false; // Something odd so lets skip clever fast branching
00952   }
00953 /*
00954   Is strong branching enabled? If so, set up and do it. Otherwise, we'll
00955   fall through to simple branching.
00956 
00957   Setup for strong branching involves saving the current basis (for restoration
00958   afterwards) and setting up for hot starts.
00959 */
00960   if (model->numberStrong()) {
00961     
00962     bool solveAll=false; // set true to say look at all even if some fixed (experiment)
00963     // Save basis
00964     CoinWarmStart * ws = solver->getWarmStart();
00965     // save limit
00966     int saveLimit;
00967     solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
00968     if (beforeSolution)
00969       solver->setIntParam(OsiMaxNumIterationHotStart,10000); // go to end
00970 
00971     /* If we are doing all strong branching in one go then we create new arrays
00972        to store information.  If clp NULL then doing old way.
00973        Going down -
00974        outputSolution[2*i] is final solution.
00975        outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown
00976        outputStuff[2*i+numberStrong] is number iterations
00977        On entry newUpper[i] is new upper bound, on exit obj change
00978        Going up -
00979        outputSolution[2*i+1] is final solution.
00980        outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
00981        outputStuff[2*i+1+numberStrong] is number iterations
00982        On entry newLower[i] is new lower bound, on exit obj change
00983     */
00984     OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
00985     ClpSimplex * clp=NULL;
00986     double * newLower = NULL;
00987     double * newUpper = NULL;
00988     double ** outputSolution=NULL;
00989     int * outputStuff=NULL;
00990     int saveLogLevel;
00991     // For moment - until full testing - go back to normal way
00992     //allNormal=false;
00993     if (osiclp&& allNormal) {
00994       clp = osiclp->getModelPtr();
00995       saveLogLevel = clp->logLevel();
00996       int saveMaxIts = clp->maximumIterations();
00997       clp->setMaximumIterations(saveLimit);
00998       clp->setLogLevel(0);
00999       // Clp - do a different way
01000       newLower = new double[numberStrong];
01001       newUpper = new double[numberStrong];
01002       outputSolution = new double * [2*numberStrong];
01003       outputStuff = new int [4*numberStrong];
01004       int * which = new int[numberStrong];
01005       for (i=0;i<numberStrong;i++) {
01006         int iObject = choice[i].objectNumber;
01007         const SbbObject * object = model->object(iObject);
01008         const SbbSimpleInteger * simple = dynamic_cast <const SbbSimpleInteger *> (object);
01009         int iSequence = simple->modelSequence();
01010         newLower[i]= ceil(saveSolution[iSequence]);
01011         newUpper[i]= floor(saveSolution[iSequence]);
01012         which[i]=iSequence;
01013         outputSolution[2*i]= new double [numberColumns];
01014         outputSolution[2*i+1]= new double [numberColumns];
01015       }
01016       clp->strongBranching(numberStrong,which,
01017                            newLower, newUpper,outputSolution,
01018                            outputStuff,outputStuff+2*numberStrong,!solveAll,false);
01019       clp->setMaximumIterations(saveMaxIts);
01020       delete [] which;
01021     } else {
01022       // Doing normal way
01023       // Mark hot start
01024       solver->markHotStart();
01025     }
01026 /*
01027   Open a loop to do the strong branching LPs. For each candidate variable,
01028   solve an LP with the variable forced down, then up. If a direction turns
01029   out to be infeasible or monotonic (i.e., over the dual objective cutoff),
01030   force the objective change to be big (1.0e100). If we determine the problem
01031   is infeasible, or find a monotone variable, escape the loop.
01032 
01033   TODO: The `restore bounds' part might be better encapsulated as an
01034         unbranch() method. Branching objects more exotic than simple integers
01035         or cliques might not restrict themselves to variable bounds.
01036 
01037   TODO: Virtuous solvers invalidate the current solution (or give bogus
01038         results :-) when the bounds are changed out from under them. So we
01039         need to do all the work associated with finding a new solution before
01040         restoring the bounds.
01041 */
01042     for (i = 0 ; i < numberStrong ; i++)
01043     { double objectiveChange ;
01044       double newObjectiveValue=1.0e100;
01045       // status is 0 finished, 1 infeasible and other
01046       int iStatus;
01047 /*
01048   Try the down direction first. (Specify the initial branching alternative as
01049   down with a call to way(-1). Each subsequent call to branch() performs the
01050   specified branch and advances the branch object state to the next branch
01051   alternative.)
01052 */
01053       if (!clp) {
01054         choice[i].possibleBranch->way(-1) ;
01055         choice[i].possibleBranch->branch() ;
01056         solver->solveFromHotStart() ;
01057         // restore bounds
01058         for (int j=0;j<numberColumns;j++) {
01059           if (saveLower[j] != lower[j])
01060             solver->setColLower(j,saveLower[j]);
01061           if (saveUpper[j] != upper[j])
01062             solver->setColUpper(j,saveUpper[j]);
01063         }
01064       
01065 /*
01066   We now have an estimate of objective degradation that we can use for strong
01067   branching. If we're over the cutoff, the variable is monotone up.
01068   If we actually made it to optimality, check for a solution, and if we have
01069   a good one, call setBestSolution to process it. Note that this may reduce the
01070   cutoff, so we check again to see if we can declare this variable monotone.
01071 */
01072         if (solver->isProvenOptimal())
01073           iStatus=0; // optimal
01074         else if (solver->isIterationLimitReached()
01075                  &&!solver->isDualObjectiveLimitReached())
01076           iStatus=2; // unknown 
01077         else
01078           iStatus=1; // infeasible
01079         newObjectiveValue = solver->getObjSense()*solver->getObjValue();
01080         choice[i].numItersDown = solver->getIterationCount();
01081         objectiveChange = newObjectiveValue-objectiveValue ;
01082       } else {
01083         iStatus = outputStuff[2*i];
01084         choice[i].numItersDown = outputStuff[2*numberStrong+2*i];
01085         newObjectiveValue = objectiveValue+newUpper[i];
01086         solver->setColSolution(outputSolution[2*i]);
01087       }
01088       objectiveChange = newObjectiveValue  - objectiveValue;
01089       if (!iStatus) {
01090         choice[i].finishedDown = true ;
01091         if (newObjectiveValue>model->getCutoff()) {
01092           objectiveChange = 1.0e100; // say infeasible
01093         } else {
01094           // See if integer solution
01095           if (model->feasibleSolution(choice[i].numIntInfeasDown,
01096                                       choice[i].numObjInfeasDown))
01097             { model->setBestSolution(SBB_STRONGSOL,
01098                                      newObjectiveValue,
01099                                      solver->getColSolution()) ;
01100             if (newObjectiveValue > model->getCutoff()) //  *new* cutoff
01101               objectiveChange = 1.0e100 ;
01102             }
01103         }
01104       } else if (iStatus==1) {
01105         objectiveChange = 1.0e100 ;
01106       } else {
01107         // Can't say much as we did not finish
01108         choice[i].finishedDown = false ;
01109       }
01110       choice[i].downMovement = objectiveChange ;
01111       //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
01112       //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersDown,
01113       //     choice[i].downMovement,choice[i].finishedDown,choice[i].numIntInfeasDown,
01114       //     choice[i].numObjInfeasDown);
01115 
01116       // repeat the whole exercise, forcing the variable up
01117       if (!clp) {
01118         choice[i].possibleBranch->branch();
01119         solver->solveFromHotStart();
01120         // restore bounds
01121         for (int j=0;j<numberColumns;j++) {
01122           if (saveLower[j] != lower[j])
01123             solver->setColLower(j,saveLower[j]);
01124           if (saveUpper[j] != upper[j])
01125             solver->setColUpper(j,saveUpper[j]);
01126         }
01127         if (solver->isProvenOptimal())
01128           iStatus=0; // optimal
01129         else if (solver->isIterationLimitReached()
01130                  &&!solver->isDualObjectiveLimitReached())
01131           iStatus=2; // unknown 
01132         else
01133           iStatus=1; // infeasible
01134         newObjectiveValue = solver->getObjSense()*solver->getObjValue();
01135         choice[i].numItersUp = solver->getIterationCount();
01136         objectiveChange = newObjectiveValue-objectiveValue ;
01137       } else {
01138         iStatus = outputStuff[2*i+1];
01139         choice[i].numItersUp = outputStuff[2*numberStrong+2*i+1];
01140         newObjectiveValue = objectiveValue+newLower[i];
01141         solver->setColSolution(outputSolution[2*i+1]);
01142       }
01143       objectiveChange = newObjectiveValue  - objectiveValue;
01144       if (!iStatus) {
01145         choice[i].finishedUp = true ;
01146         if (newObjectiveValue>model->getCutoff()) {
01147           objectiveChange = 1.0e100; // say infeasible
01148         } else {
01149           // See if integer solution
01150           if (model->feasibleSolution(choice[i].numIntInfeasUp,
01151                                       choice[i].numObjInfeasUp))
01152             { model->setBestSolution(SBB_STRONGSOL,
01153                                      newObjectiveValue,
01154                                      solver->getColSolution()) ;
01155             if (newObjectiveValue > model->getCutoff()) //  *new* cutoff
01156               objectiveChange = 1.0e100 ;
01157             }
01158         }
01159       } else if (iStatus==1) {
01160         objectiveChange = 1.0e100 ;
01161       } else {
01162         // Can't say much as we did not finish
01163         choice[i].finishedUp = false ;
01164       }
01165       choice[i].upMovement = objectiveChange ;
01166       //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
01167       //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersUp,
01168       //     choice[i].upMovement,choice[i].finishedUp,choice[i].numIntInfeasUp,
01169       //     choice[i].numObjInfeasUp);
01170 
01171 /*
01172   End of evaluation for this candidate variable. Possibilities are:
01173     * Both sides below cutoff; this variable is a candidate for branching.
01174     * Both sides infeasible or above the objective cutoff: no further action
01175       here. Break from the evaluation loop and assume the node will be purged
01176       by the caller.
01177     * One side below cutoff: Install the branch (i.e., fix the variable). Break
01178       from the evaluation loop and assume the node will be reoptimised by the
01179       caller.
01180 */
01181       if (choice[i].upMovement<1.0e100) {
01182         if(choice[i].downMovement<1.0e100) {
01183           // feasible - no action
01184         } else {
01185           // up feasible, down infeasible
01186           anyAction=-1;
01187           if (!solveAll) {
01188             choice[i].possibleBranch->way(1);
01189             choice[i].possibleBranch->branch();
01190             break;
01191           }
01192         }
01193       } else {
01194         if(choice[i].downMovement<1.0e100) {
01195           // down feasible, up infeasible
01196           anyAction=-1;
01197           if (!solveAll) {
01198             choice[i].possibleBranch->way(-1);
01199             choice[i].possibleBranch->branch();
01200             break;
01201           }
01202         } else {
01203           // neither side feasible
01204           anyAction=-2;
01205           break;
01206         }
01207       }
01208     }
01209 
01210     if (!clp) {
01211       // Delete the snapshot
01212       solver->unmarkHotStart();
01213     } else {
01214       clp->setLogLevel(saveLogLevel);
01215       delete [] newLower;
01216       delete [] newUpper;
01217       delete [] outputStuff;
01218       int i;
01219       for (i=0;i<2*numberStrong;i++)
01220         delete [] outputSolution[i];
01221       delete [] outputSolution;
01222     }
01223     solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
01224     // restore basis
01225     solver->setWarmStart(ws);
01226     delete ws;
01227 
01228 /*
01229   anyAction >= 0 indicates that strong branching didn't produce any monotone
01230   variables. Sift through the candidates for the best one.
01231 
01232   QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is
01233          local to this code block. Perhaps should be numberNodes_ from model?
01234          Unclear what this calculation is doing.
01235 */
01236     if (anyAction>=0) {
01237 
01238       int numberNodes = model->getNodeCount();
01239       // get average cost per iteration and assume stopped ones
01240       // would stop after 50% more iterations at average cost??? !!! ???
01241       double averageCostPerIteration=0.0;
01242       double totalNumberIterations=1.0;
01243       int smallestNumberInfeasibilities=INT_MAX;
01244       for (i=0;i<numberStrong;i++) {
01245         totalNumberIterations += choice[i].numItersDown +
01246           choice[i].numItersUp ;
01247         averageCostPerIteration += choice[i].downMovement +
01248           choice[i].upMovement;
01249         smallestNumberInfeasibilities= 
01250           min(min(choice[i].numIntInfeasDown ,
01251                   choice[i].numIntInfeasUp ),
01252               smallestNumberInfeasibilities);
01253       }
01254       if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities)
01255         numberNodes=1000000; // switch off search for better solution
01256       numberNodes=1000000; // switch off anyway
01257       averageCostPerIteration /= totalNumberIterations;
01258       // all feasible - choose best bet
01259 
01260 #if 0
01261       for (i = 0 ; i < numberStrong ; i++)
01262       { int iColumn =
01263           model->integerVariable()[choice[i].possibleBranch->variable()] ;
01264          model->messageHandler()->message(SBB_STRONG,model->messages())
01265           << i << iColumn
01266           <<choice[i].downMovement<<choice[i].numIntInfeasDown 
01267           <<choice[i].upMovement<<choice[i].numIntInfeasUp 
01268           <<choice[i].possibleBranch->value()
01269           <<CoinMessageEol;
01270         int betterWay = decision->betterBranch(choice[i].possibleBranch,
01271                                               branch_,
01272                                               choice[i].upMovement,
01273                                               choice[i].numIntInfeasUp ,
01274                                               choice[i].downMovement,
01275                                               choice[i].numIntInfeasDown );
01276         if (betterWay) {
01277           delete branch_;
01278           // C) create branching object
01279           branch_ = choice[i].possibleBranch;
01280           choice[i].possibleBranch=NULL;
01281           branch_->way(betterWay);
01282         }
01283       }
01284 #else
01285       // New method does all at once so it can be more sophisticated
01286       // in deciding how to balance actions.
01287       // But it does need arrays
01288       double * changeUp = new double [numberStrong];
01289       int * numberInfeasibilitiesUp = new int [numberStrong];
01290       double * changeDown = new double [numberStrong];
01291       int * numberInfeasibilitiesDown = new int [numberStrong];
01292       SbbBranchingObject ** objects = new SbbBranchingObject * [ numberStrong];
01293       for (i = 0 ; i < numberStrong ; i++) {
01294         int iColumn =
01295           model->integerVariable()[choice[i].possibleBranch->variable()] ;
01296         model->messageHandler()->message(SBB_STRONG,model->messages())
01297           << i << iColumn
01298           <<choice[i].downMovement<<choice[i].numIntInfeasDown 
01299           <<choice[i].upMovement<<choice[i].numIntInfeasUp 
01300           <<choice[i].possibleBranch->value()
01301           <<CoinMessageEol;
01302         changeUp[i]=choice[i].upMovement;
01303         numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
01304         changeDown[i]=choice[i].downMovement;
01305         numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
01306         objects[i] = choice[i].possibleBranch;
01307       }
01308       int whichObject = decision->bestBranch(objects,numberStrong,numberUnsatisfied_,
01309                                              changeUp,numberInfeasibilitiesUp,
01310                                              changeDown,numberInfeasibilitiesDown,
01311                                              objectiveValue);
01312       // move branching object and make sure it will not be deleted
01313       if (whichObject>=0) {
01314         branch_ = objects[whichObject];
01315         choice[whichObject].possibleBranch=NULL;
01316       }
01317       delete [] changeUp;
01318       delete [] numberInfeasibilitiesUp;
01319       delete [] changeDown;
01320       delete [] numberInfeasibilitiesDown;
01321       delete [] objects;
01322 #endif 
01323     }
01324   }
01325 /*
01326   Simple branching. Why are we choosing choice[0]? Seems arbitrary. Why not
01327   the final candidate? (It will have the greatest infeasibility.)
01328 
01329   (jjf) you can choose the first or last as if we get here choice has length 1!
01330 */
01331   else {
01332     // not strong
01333     // C) create branching object
01334     branch_ = choice[0].possibleBranch;
01335     choice[0].possibleBranch=NULL;
01336     // we could use otherWay - think a bit
01337   }
01338 /*
01339   Cleanup, then we're outta here.
01340 */
01341   if (!model->branchingMethod())
01342     delete decision;
01343     
01344   for (i=0;i<maximumStrong;i++)
01345     delete choice[i].possibleBranch;
01346   delete [] choice;
01347   delete [] saveLower;
01348   delete [] saveUpper;
01349 
01350   // restore solution
01351   solver->setColSolution(saveSolution);
01352   delete [] saveSolution;
01353   return anyAction;
01354 }
01355 
01356 
01357 SbbNode::SbbNode(const SbbNode & rhs) 
01358 {  
01359 #ifdef CHECK_NODE
01360   printf("SbbNode %x Constructor from rhs %x\n",this,&rhs);
01361 #endif
01362   if (rhs.nodeInfo_)
01363     nodeInfo_ = rhs.nodeInfo_->clone();
01364   else
01365     nodeInfo_=NULL;
01366   objectiveValue_=rhs.objectiveValue_;
01367   guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
01368   if (rhs.branch_)
01369     branch_=rhs.branch_->clone();
01370   else
01371     branch_=NULL,
01372   depth_ = rhs.depth_;
01373   numberUnsatisfied_ = rhs.numberUnsatisfied_;
01374 }
01375 
01376 SbbNode &
01377 SbbNode::operator=(const SbbNode & rhs)
01378 {
01379   if (this != &rhs) {
01380     delete nodeInfo_;
01381     if (nodeInfo_)
01382       nodeInfo_ = rhs.nodeInfo_->clone();
01383     else
01384       nodeInfo_ = NULL;
01385     objectiveValue_=rhs.objectiveValue_;
01386     guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
01387     if (rhs.branch_)
01388       branch_=rhs.branch_->clone();
01389     else
01390       branch_=NULL,
01391     depth_ = rhs.depth_;
01392     numberUnsatisfied_ = rhs.numberUnsatisfied_;
01393   }
01394   return *this;
01395 }
01396 
01397 
01398 SbbNode::~SbbNode ()
01399 {
01400 #ifdef CHECK_NODE
01401   if (nodeInfo_) 
01402     printf("SbbNode %x Destructor nodeInfo %x (%d)\n",
01403          this,nodeInfo_,nodeInfo_->numberPointingToThis());
01404   else
01405     printf("SbbNode %x Destructor nodeInfo %x (?)\n",
01406          this,nodeInfo_);
01407 #endif
01408   if (nodeInfo_) {
01409     nodeInfo_->nullOwner();
01410     int numberToDelete=nodeInfo_->numberBranchesLeft();
01411     //    SbbNodeInfo * parent = nodeInfo_->parent();
01412     if (nodeInfo_->decrement(numberToDelete)==0) {
01413       delete nodeInfo_;
01414     } else {
01415       //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,parent);
01416       // anyway decrement parent
01417       //if (parent)
01419     }
01420   }
01421   delete branch_;
01422 }
01423 // Decrement  active cut counts 
01424 void 
01425 SbbNode::decrementCuts(int change)
01426 {
01427   if(nodeInfo_) {
01428     nodeInfo_->decrementCuts(change);
01429   }
01430 }
01431 void 
01432 SbbNode::decrementParentCuts(int change)
01433 {
01434   if(nodeInfo_) {
01435     nodeInfo_->decrementParentCuts(change);
01436   }
01437 }
01438 
01439 /*
01440   Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
01441   in the attached nodeInfo_.
01442 */
01443 void
01444 SbbNode::initializeInfo()
01445 {
01446   assert(nodeInfo_ && branch_) ;
01447   nodeInfo_->initializeInfo(branch_->numberBranches());
01448 }
01449 // Nulls out node info
01450 void 
01451 SbbNode::nullNodeInfo()
01452 {
01453   nodeInfo_=NULL;
01454 }
01455 
01456 int
01457 SbbNode::branch()
01458 {
01459   branch_->branch();
01460   return nodeInfo_->branchedOn();
01461 }

Generated on Wed Dec 3 14:36:21 2003 for Sbb by doxygen 1.3.5