00001
00002
00003 #if defined(_MSC_VER)
00004
00005 # pragma warning(disable:4786)
00006 #endif
00007 #include <string>
00008
00009
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
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
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
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
00102 void
00103 SbbNodeInfo::decrementCuts(int change)
00104 {
00105 int i;
00106
00107 int changeThis;
00108 if (change<0)
00109 changeThis = numberBranchesLeft_;
00110 else
00111 changeThis = change;
00112
00113 for (i=0;i<numberCuts_;i++) {
00114 if (cuts_[i]) {
00115 int number = cuts_[i]->decrement(changeThis);
00116 if (!number) {
00117
00118 delete cuts_[i];
00119 cuts_[i]=NULL;
00120 }
00121 }
00122 }
00123 }
00124 void
00125 SbbNodeInfo::decrementParentCuts(int change)
00126 {
00127 if (parent_) {
00128
00129 int changeThis;
00130 if (change<0)
00131 changeThis = numberBranchesLeft_;
00132 else
00133 changeThis = change;
00134 int i;
00135
00136 CoinWarmStartBasis dummy;
00137 dummy.setSize(0,numberRows_+numberCuts_);
00138 buildRowBasis(dummy);
00139
00140
00141 SbbNodeInfo * thisInfo = parent_;
00142 while (thisInfo)
00143 thisInfo = thisInfo->buildRowBasis(dummy);
00144
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
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
00179 CoinWarmStartBasis dummy;
00180 dummy.setSize(0,numberRows_+numberCuts_);
00181
00182
00183 buildRowBasis(dummy);
00184 SbbNodeInfo * thisInfo = parent_;
00185 while (thisInfo)
00186 thisInfo = thisInfo->buildRowBasis(dummy);
00187
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
00207
00208
00209
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
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
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
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
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
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
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437 void SbbFullNodeInfo::applyToModel (SbbModel *model,
00438 CoinWarmStartBasis *&basis,
00439 SbbCountRowCut **addCuts,
00440 int ¤tNumberCuts) const
00441
00442 { OsiSolverInterface *solver = model->solver() ;
00443
00444
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
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
00463
00464
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
00483 SbbPartialNodeInfo::SbbPartialNodeInfo()
00484
00485 : SbbNodeInfo(),
00486 basisDiff_(NULL),
00487 variables_(NULL),
00488 newBounds_(NULL),
00489 numberChangedBounds_(0)
00490
00491 { }
00492
00493
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 ¤tNumberCuts) const
00556
00557 { OsiSolverInterface *solver = model->solver();
00558
00559 basis->applyDiff(basisDiff_) ;
00560
00561
00562 int i;
00563 for (i=0;i<numberChangedBounds_;i++) {
00564 int variable = variables_[i];
00565 if ((variable&0x80000000)==0) {
00566
00567 solver->setColLower(variable,newBounds_[i]);
00568 } else {
00569
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
00580
00581
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
00634
00635 if (!lastNode)
00636 { nodeInfo_=new SbbFullNodeInfo(model,solver->getNumRows()); }
00637
00638
00639
00640
00641
00642
00643
00644 else
00645 { const CoinWarmStartBasis* ws =
00646 dynamic_cast<const CoinWarmStartBasis*>(solver->getWarmStart());
00647 assert(ws!=NULL);
00648
00649 int numberColumns = solver->getNumCols();
00650
00651 const double * lower = solver->getColLower();
00652 const double * upper = solver->getColUpper();
00653
00654 int i;
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680 CoinWarmStartBasis *expanded =
00681 dynamic_cast<CoinWarmStartBasis *>(ws->clone()) ;
00682 int numberRowsAtContinuous = model->numberRowsAtContinuous();
00683 int iFull = numberRowsAtContinuous+model->currentNumberCuts()+
00684 numberNewCuts;
00685
00686
00687
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
00697
00698
00699
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
00728
00729
00730
00731 CoinWarmStartDiff *basisDiff = expanded->generateDiff(lastws) ;
00732
00733
00734
00735
00736
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
00765
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
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
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
00830
00831 double * saveSolution = new double[numberColumns];
00832 memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
00833
00834
00835
00836
00837
00838 SbbBranchDecision *decision = model->branchingMethod();
00839 if (!decision)
00840 decision = new SbbBranchDefaultDecision();
00841
00842 typedef struct {
00843 SbbBranchingObject * possibleBranch;
00844 double upMovement;
00845 double downMovement;
00846 int numIntInfeasUp ;
00847 int numObjInfeasUp ;
00848 bool finishedUp;
00849 int numItersUp ;
00850 int numIntInfeasDown ;
00851 int numObjInfeasDown ;
00852 bool finishedDown;
00853 int numItersDown;
00854 int objectNumber;
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
00863 int numberIntegerInfeasibilities;
00864 int numberObjectInfeasibilities;
00865 model->feasibleSolution(
00866 numberIntegerInfeasibilities,
00867 numberObjectInfeasibilities);
00868
00869
00870
00871
00872 assert(model->getForcePriority()<0);
00873
00874 double saveOtherWay=0.0;
00875 numberUnsatisfied_ = 0;
00876 int bestPriority=INT_MAX;
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
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
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
00920 if (infeasibility>mostAway) {
00921
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
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
00943
00944
00945
00946 bool allNormal=true;
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;
00952 }
00953
00954
00955
00956
00957
00958
00959
00960 if (model->numberStrong()) {
00961
00962 bool solveAll=false;
00963
00964 CoinWarmStart * ws = solver->getWarmStart();
00965
00966 int saveLimit;
00967 solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
00968 if (beforeSolution)
00969 solver->setIntParam(OsiMaxNumIterationHotStart,10000);
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
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
00992
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
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
01023
01024 solver->markHotStart();
01025 }
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042 for (i = 0 ; i < numberStrong ; i++)
01043 { double objectiveChange ;
01044 double newObjectiveValue=1.0e100;
01045
01046 int iStatus;
01047
01048
01049
01050
01051
01052
01053 if (!clp) {
01054 choice[i].possibleBranch->way(-1) ;
01055 choice[i].possibleBranch->branch() ;
01056 solver->solveFromHotStart() ;
01057
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
01067
01068
01069
01070
01071
01072 if (solver->isProvenOptimal())
01073 iStatus=0;
01074 else if (solver->isIterationLimitReached()
01075 &&!solver->isDualObjectiveLimitReached())
01076 iStatus=2;
01077 else
01078 iStatus=1;
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;
01093 } else {
01094
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())
01101 objectiveChange = 1.0e100 ;
01102 }
01103 }
01104 } else if (iStatus==1) {
01105 objectiveChange = 1.0e100 ;
01106 } else {
01107
01108 choice[i].finishedDown = false ;
01109 }
01110 choice[i].downMovement = objectiveChange ;
01111
01112
01113
01114
01115
01116
01117 if (!clp) {
01118 choice[i].possibleBranch->branch();
01119 solver->solveFromHotStart();
01120
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;
01129 else if (solver->isIterationLimitReached()
01130 &&!solver->isDualObjectiveLimitReached())
01131 iStatus=2;
01132 else
01133 iStatus=1;
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;
01148 } else {
01149
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())
01156 objectiveChange = 1.0e100 ;
01157 }
01158 }
01159 } else if (iStatus==1) {
01160 objectiveChange = 1.0e100 ;
01161 } else {
01162
01163 choice[i].finishedUp = false ;
01164 }
01165 choice[i].upMovement = objectiveChange ;
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181 if (choice[i].upMovement<1.0e100) {
01182 if(choice[i].downMovement<1.0e100) {
01183
01184 } else {
01185
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
01196 anyAction=-1;
01197 if (!solveAll) {
01198 choice[i].possibleBranch->way(-1);
01199 choice[i].possibleBranch->branch();
01200 break;
01201 }
01202 } else {
01203
01204 anyAction=-2;
01205 break;
01206 }
01207 }
01208 }
01209
01210 if (!clp) {
01211
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
01225 solver->setWarmStart(ws);
01226 delete ws;
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236 if (anyAction>=0) {
01237
01238 int numberNodes = model->getNodeCount();
01239
01240
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;
01256 numberNodes=1000000;
01257 averageCostPerIteration /= totalNumberIterations;
01258
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
01279 branch_ = choice[i].possibleBranch;
01280 choice[i].possibleBranch=NULL;
01281 branch_->way(betterWay);
01282 }
01283 }
01284 #else
01285
01286
01287
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
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
01327
01328
01329
01330
01331 else {
01332
01333
01334 branch_ = choice[0].possibleBranch;
01335 choice[0].possibleBranch=NULL;
01336
01337 }
01338
01339
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
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
01412 if (nodeInfo_->decrement(numberToDelete)==0) {
01413 delete nodeInfo_;
01414 } else {
01415
01416
01417
01419 }
01420 }
01421 delete branch_;
01422 }
01423
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
01441
01442
01443 void
01444 SbbNode::initializeInfo()
01445 {
01446 assert(nodeInfo_ && branch_) ;
01447 nodeInfo_->initializeInfo(branch_->numberBranches());
01448 }
01449
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 }