00001
00002
00003
00004 #include <climits>
00005 #include "CoinPragma.hpp"
00006 #include "BonChooseVariable.hpp"
00007 #include "CoinTime.hpp"
00008 #include "IpBlas.hpp"
00009 #include "BonMsgUtils.hpp"
00010
00011
00012 #include "CbcModel.hpp"
00013
00014 namespace Bonmin
00015 {
00016
00017 BonChooseVariable::Messages::Messages():
00018 CoinMessages((int) BON_CHOOSE_MESSAGES_DUMMY_END)
00019 {
00020 strcpy(source_,"BON");
00021 ADD_MSG(PS_COST_HISTORY,std_m,6,"%3d up %3d %.8e down %3d %.8e");
00022 ADD_MSG(PS_COST_MULT,std_m, 6, "upMultiplier = %e downMultiplier = %e");
00023 ADD_MSG(PS_COST_ESTIMATES, std_m, 6, "%3d value = %e upEstimate = %e downEstimate = %e infeas = %e value2 = %e");
00024 ADD_MSG(CANDIDATE_LIST,std_m,6,
00025 "list_[%5d] = %5d, usefull_[%5d] = %.16e %.16e");
00026 ADD_MSG(CANDIDATE_LIST2, std_m, 6,
00027 "list_[%3d] = %3d useful_[%3d] = %e");
00028 ADD_MSG(CANDIDATE_LIST3, std_m, 6,
00029 "list2[%3d] = %3d useful2[%3d] = %e");
00030 ADD_MSG(SB_START, std_m,5,
00031 " Starting strong branching. Obj. val = %g\n");
00032 ADD_MSG(SB_HEADER, std_m,5,
00033 " Var Value DownStat DownChange UpStat UpChange");
00034 ADD_MSG(SB_RES, std_m, 5,
00035 " %3d %3d %.6e %6s %.6e %6s %.6e");
00036 ADD_MSG(BRANCH_VAR, std_m, 4, "Branched on variable %i, bestWhichWay: %i");
00037 ADD_MSG(CHOSEN_VAR, std_m, 4," Choosing %d");
00038 ADD_MSG(UPDATE_PS_COST, std_m, 4,"update %3d %3d %e %e %3d");
00039 }
00040 const std::string BonChooseVariable::CNAME = "BonChooseVariable";
00041
00042 BonChooseVariable::BonChooseVariable(BabSetupBase &b, const OsiSolverInterface* solver):
00043 OsiChooseVariable(solver),
00044 results_(),
00045 cbc_model_(NULL),
00046 only_pseudo_when_trusted_(false),
00047 pseudoCosts_()
00048 {
00049 jnlst_ = b.journalist();
00050 Ipopt::SmartPtr<Ipopt::OptionsList> options = b.options();
00051
00052 handler_ = new CoinMessageHandler;
00053
00054 options->GetIntegerValue("bb_log_level", bb_log_level_, b.prefix());
00055 handler_->setLogLevel(bb_log_level_);
00056 options->GetNumericValue("time_limit", time_limit_, b.prefix());
00057 options->GetNumericValue("setup_pseudo_frac", setup_pseudo_frac_, b.prefix());
00058 options->GetNumericValue("maxmin_crit_no_sol", maxmin_crit_no_sol_, b.prefix());
00059 options->GetNumericValue("maxmin_crit_have_sol", maxmin_crit_have_sol_, b.prefix());
00060 options->GetEnumValue("trust_strong_branching_for_pseudo_cost",trustStrongForPseudoCosts_ , b.prefix());
00061 int sortCrit;
00062 options->GetEnumValue("candidate_sort_criterion", sortCrit, b.prefix());
00063 #ifndef OLD_USEFULLNESS
00064 sortCrit_ = (CandidateSortCriterion) sortCrit;
00065 #endif
00066
00067 int numberObjects = solver_->numberObjects();
00068
00069 pseudoCosts_.initialize(numberObjects);
00070 int numberBeforeTrusted = b.getIntParameter(BabSetupBase::MinReliability);
00071 pseudoCosts_.setNumberBeforeTrusted(numberBeforeTrusted);
00072
00073 setNumberStrong(b.getIntParameter(BabSetupBase::NumberStrong));
00074
00076 if (!options->GetIntegerValue("number_before_trust_list", numberBeforeTrustedList_, b.prefix())) {
00077
00078 numberBeforeTrustedList_ = numberBeforeTrusted;
00079 }
00080 options->GetIntegerValue("number_strong_branch_root", numberStrongRoot_, b.prefix());
00081 options->GetIntegerValue("min_number_strong_branch", minNumberStrongBranch_, b.prefix());
00082 options->GetIntegerValue("number_look_ahead", numberLookAhead_, b.prefix());
00083
00084 start_time_ = CoinCpuTime();
00085 }
00086
00087 BonChooseVariable::BonChooseVariable(const BonChooseVariable & rhs) :
00088 OsiChooseVariable(rhs),
00089 results_(rhs.results_),
00090 time_limit_(rhs.time_limit_),
00091 start_time_(CoinCpuTime()),
00092 cbc_model_(rhs.cbc_model_),
00093 only_pseudo_when_trusted_(rhs.only_pseudo_when_trusted_),
00094 maxmin_crit_no_sol_(rhs.maxmin_crit_no_sol_),
00095 maxmin_crit_have_sol_(rhs.maxmin_crit_have_sol_),
00096 setup_pseudo_frac_(rhs.setup_pseudo_frac_),
00097 numberBeforeTrustedList_(rhs.numberBeforeTrustedList_),
00098 numberStrongRoot_(rhs.numberStrongRoot_),
00099 #ifndef OLD_USEFULLNESS
00100 sortCrit_(rhs.sortCrit_),
00101 #endif
00102 numberLookAhead_(rhs.numberLookAhead_),
00103 minNumberStrongBranch_(rhs.minNumberStrongBranch_),
00104 pseudoCosts_(rhs.pseudoCosts_),
00105 trustStrongForPseudoCosts_(rhs.trustStrongForPseudoCosts_)
00106 {
00107 jnlst_ = rhs.jnlst_;
00108 handler_ = rhs.handler_->clone();
00109 bb_log_level_ = rhs.bb_log_level_;
00110 DBG_ASSERT(IsValid(jnlst_));
00111 }
00112
00113 BonChooseVariable &
00114 BonChooseVariable::operator=(const BonChooseVariable & rhs)
00115 {
00116 if (this != &rhs) {
00117 OsiChooseVariable::operator=(rhs);
00118 delete handler_;
00119 handler_ = rhs.handler_->clone();
00120 jnlst_ = rhs.jnlst_;
00121 bb_log_level_ = rhs.bb_log_level_;
00122 cbc_model_ = rhs.cbc_model_;
00123 only_pseudo_when_trusted_ = rhs.only_pseudo_when_trusted_;
00124 maxmin_crit_no_sol_ = rhs.maxmin_crit_no_sol_;
00125 maxmin_crit_have_sol_ = rhs.maxmin_crit_have_sol_;
00126 setup_pseudo_frac_ = rhs.setup_pseudo_frac_;
00127 numberBeforeTrustedList_ = rhs.numberBeforeTrustedList_;
00128 numberStrongRoot_ = rhs.numberStrongRoot_;
00129 #ifndef OLD_USEFULLNESS
00130 sortCrit_ = rhs.sortCrit_;
00131 #endif
00132 minNumberStrongBranch_ = rhs.minNumberStrongBranch_;
00133 pseudoCosts_ = rhs.pseudoCosts_;
00134 trustStrongForPseudoCosts_ = rhs.trustStrongForPseudoCosts_;
00135 numberLookAhead_ = rhs.numberLookAhead_;
00136 results_ = rhs.results_;
00137 }
00138 return *this;
00139 }
00140
00141 OsiChooseVariable *
00142 BonChooseVariable::clone() const
00143 {
00144 return new BonChooseVariable(*this);
00145 }
00146
00147 BonChooseVariable::~BonChooseVariable ()
00148 {
00149 delete handler_;
00150 }
00151
00152 void
00153 BonChooseVariable::registerOptions(
00154 Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)
00155 {
00156 roptions->SetRegisteringCategory("Strong branching setup", RegisteredOptions::BonminCategory);
00157 roptions->AddStringOption4("candidate_sort_criterion",
00158 "Choice of the criterion to choose candidates in strong-branching",
00159 "best-ps-cost",
00160 "best-ps-cost","Sort by decreasing pseudo-cost",
00161 "worst-ps-cost", "Sort by increasing pseudo-cost",
00162 "most-fractional", "Sort by decreasing integer infeasibility",
00163 "least-fractional", "Sort by increasing integer infeasibility","");
00164
00165 roptions->setOptionExtraInfo("candidate_sort_criterion",63);
00166
00167 roptions->AddBoundedNumberOption("setup_pseudo_frac", "Proportion of strong branching list that has to be taken from most-integer-infeasible list.",
00168 0., false, 1., false, 0.5);
00169 roptions->setOptionExtraInfo("setup_pseudo_frac",63);
00170 roptions->AddBoundedNumberOption("maxmin_crit_no_sol", "Weight towards minimum in of lower and upper branching estimates when no solution has been found yet.",
00171 0., false, 1., false, 0.7);
00172 roptions->setOptionExtraInfo("maxmin_crit_no_sol",63);
00173 roptions->AddBoundedNumberOption("maxmin_crit_have_sol", "Weight towards minimum in of lower and upper branching estimates when a solution has been found.",
00174 0., false, 1., false, 0.1);
00175 roptions->setOptionExtraInfo("maxmin_crit_have_sol",63);
00176 roptions->AddLowerBoundedIntegerOption("number_before_trust_list",
00177 "Set the number of branches on a variable before its pseudo costs are to be believed during setup of strong branching candidate list.",
00178 -1, 0, "The default value is that of \"number_before_trust\"");
00179 roptions->setOptionExtraInfo("number_before_trust_list",63);
00180 roptions->AddLowerBoundedIntegerOption("number_strong_branch_root",
00181 "Maximum number of variables considered for strong branching in root node.",
00182 0, COIN_INT_MAX, "");
00183 roptions->setOptionExtraInfo("number_strong_branch_root",63);
00184
00185 roptions->AddLowerBoundedIntegerOption("min_number_strong_branch", "Sets minimum number of variables for strong branching (overriding trust)",
00186 0, 0,"");
00187 roptions->setOptionExtraInfo("min_number_strong_branch",63);
00188 roptions->AddStringOption2("trust_strong_branching_for_pseudo_cost",
00189 "Whether or not to trust strong branching results for updating pseudo costs.",
00190 "yes",
00191 "no","",
00192 "yes","",
00193 ""
00194 );
00195 roptions->setOptionExtraInfo("trust_strong_branching_for_pseudo_cost", 63);
00196
00197 roptions->AddLowerBoundedIntegerOption("number_look_ahead", "Sets limit of look-ahead strong-branching trials",
00198 0, 0,"");
00199 roptions->setOptionExtraInfo("number_look_ahead", 31);
00200 }
00201
00202
00203 void
00204 BonChooseVariable::computeMultipliers(double& upMult, double& downMult) const
00205 {
00206 const double* upTotalChange = pseudoCosts_.upTotalChange();
00207 const double* downTotalChange = pseudoCosts_.downTotalChange();
00208 const int* upNumber = pseudoCosts_.upNumber();
00209 const int* downNumber = pseudoCosts_.downNumber();
00210 double sumUp=0.0;
00211 double numberUp=0.0;
00212 double sumDown=0.0;
00213 double numberDown=0.0;
00214 for (int i=pseudoCosts_.numberObjects() - 1; i >= 0; --i) {
00215 sumUp += upTotalChange[i];
00216 numberUp += upNumber[i];
00217 sumDown += downTotalChange[i];
00218 numberDown += downNumber[i];
00219 message(PS_COST_HISTORY)
00220 <<i<< upNumber[i]<< upTotalChange[i]
00221 << downNumber[i]<< downTotalChange[i]<<CoinMessageEol;
00222 }
00223 upMult=(1.0+sumUp)/(1.0+numberUp);
00224 downMult=(1.0+sumDown)/(1.0+numberDown);
00225
00226 message(PS_COST_MULT)
00227 <<upMult<< downMult<<CoinMessageEol;
00228 }
00229
00230 double
00231 BonChooseVariable::computeUsefulness(const double MAXMIN_CRITERION,
00232 const double upMult, const double downMult,
00233 const double value,
00234 const OsiObject* object, int i,
00235 double& value2) const
00236 {
00237 #ifdef OLD_USEFULLNESS
00238 double sumUp = pseudoCosts_.upTotalChange()[i]+1.0e-30;
00239 int numberUp = pseudoCosts_.upNumber()[i];
00240 double sumDown = pseudoCosts_.downTotalChange()[i]+1.0e-30;
00241 int numberDown = pseudoCosts_.downNumber()[i];
00242 double upEst = object->upEstimate();
00243 double downEst = object->downEstimate();
00244 upEst = numberUp ? ((upEst*sumUp)/numberUp) : (upEst*upMult);
00245
00246
00247 downEst = numberDown ? ((downEst*sumDown)/numberDown) : (downEst*downMult);
00248
00249
00250 double useful = ( MAXMIN_CRITERION*CoinMin(upEst,downEst) +
00251 (1.0-MAXMIN_CRITERION)*CoinMax(upEst,downEst) );
00252 value2 = -COIN_DBL_MAX;
00253 if (numberUp < numberBeforeTrustedList_ ||
00254 numberDown < numberBeforeTrustedList_) {
00255 value2 = value;
00256 }
00257 #else
00258
00259 int numberUp = pseudoCosts_.upNumber()[i];
00260 int numberDown = pseudoCosts_.downNumber()[i];
00261 if (sortCrit_ >= DecrPs && sortCrit_ <= IncrPs) {
00262 double sumUp = pseudoCosts_.upTotalChange()[i]+1.0e-30;
00263 double sumDown = pseudoCosts_.downTotalChange()[i]+1.0e-30;
00264 double upEst = object->upEstimate();
00265 double downEst = object->downEstimate();
00266 upEst = numberUp ? ((upEst*sumUp)/numberUp) : (upEst*upMult);
00267
00268
00269 downEst = numberDown ? ((downEst*sumDown)/numberDown) : (downEst*downMult);
00270
00271
00272 double useful = ( MAXMIN_CRITERION*CoinMin(upEst,downEst) +
00273 (1.0-MAXMIN_CRITERION)*CoinMax(upEst,downEst) );
00274 if (numberUp == 0 || numberDown == 0) {
00275 if (value == 0.) useful = 0;
00276 else if (MAXMIN_CRITERION >= 0.5)
00277 useful = CoinMin(upEst, downEst);
00278 else {
00279 useful = CoinMax(upEst, downEst);
00280 }
00281 }
00282 value2 = -COIN_DBL_MAX;
00283 if (numberUp < numberBeforeTrustedList_ ||
00284 numberDown < numberBeforeTrustedList_) {
00285 value2 = value;
00286 }
00287 #endif
00288 message(PS_COST_ESTIMATES)
00289 <<i<< useful<< upEst<< downEst<< value<< value2<< CoinMessageEol;
00290 return useful;
00291 }
00292 #ifndef OLD_USEFULLNESS
00293 else if (sortCrit_ >= DecrInfeas && sortCrit_ <= IncrInfeas) {
00294 double usefull = value;
00295 value2 = value;
00296 return usefull;
00297 }
00298 else {
00299 throw CoinError("BonChooseVariable", "computeUsefullnee",
00300 "Unknown criterion for soring candidates");
00301 return COIN_DBL_MAX;
00302 }
00303 }
00304 #endif
00305
00306 int
00307 BonChooseVariable::setupList ( OsiBranchingInformation *info, bool initialize)
00308 {
00309 if (numberBeforeTrustedList_ < 0) {
00310 number_not_trusted_ = 1;
00311 return OsiChooseVariable::setupList(info, initialize);
00312 }
00313 if (initialize) {
00314 status_=-2;
00315 delete [] goodSolution_;
00316 bestObjectIndex_=-1;
00317 numberStrongDone_=0;
00318 numberStrongIterations_ = 0;
00319 numberStrongFixed_ = 0;
00320 goodSolution_ = NULL;
00321 goodObjectiveValue_ = COIN_DBL_MAX;
00322 number_not_trusted_=0;
00323 }
00324 else {
00325 throw CoinError(CNAME,"setupList","Should not be called with initialize==false");
00326 }
00327 numberOnList_=0;
00328 numberUnsatisfied_=0;
00329 int numberObjects = solver_->numberObjects();
00330 assert (numberObjects);
00331 if (numberObjects>pseudoCosts_.numberObjects()) {
00332
00333
00334
00335
00336
00337 int saveNumberBeforeTrusted = pseudoCosts_.numberBeforeTrusted();
00338 pseudoCosts_.initialize(numberObjects);
00339 pseudoCosts_.setNumberBeforeTrusted(saveNumberBeforeTrusted);
00340 }
00341 double check = -COIN_DBL_MAX;
00342 int checkIndex=0;
00343 int bestPriority=COIN_INT_MAX;
00344 int maximumStrong= CoinMin(CoinMax(numberStrong_,numberStrongRoot_),
00345 numberObjects) ;
00346 int putOther = numberObjects;
00347 int i;
00348 for (i=0;i<numberObjects;i++) {
00349 list_[i]=-1;
00350 useful_[i]=0.0;
00351 }
00352
00353 int* list2 = NULL;
00354 double* useful2 = NULL;
00355 double check2 = -COIN_DBL_MAX;
00356 int checkIndex2=0;
00357 int max_most_fra = setup_pseudo_frac_ > 0. ? (int)floor(setup_pseudo_frac_*(double)maximumStrong): 0;
00358 if (setup_pseudo_frac_ > 0.) {
00359 max_most_fra = CoinMax(1, max_most_fra);
00360 }
00361 if (max_most_fra) {
00362 list2 = new int[max_most_fra];
00363 useful2 = new double[max_most_fra];
00364 for (i=0;i<max_most_fra;i++) {
00365 list2[i]=-1;
00366 useful2[i]=0.0;
00367 }
00368 }
00369
00370 OsiObject ** object = info->solver_->objects();
00371 double upMultiplier, downMultiplier;
00372 computeMultipliers(upMultiplier, downMultiplier);
00373
00374
00375 bool feasible = true;
00376 const double MAXMIN_CRITERION = maxminCrit(info);
00377 for ( i=0;i<numberObjects;i++) {
00378 int way;
00379 double value = object[i]->infeasibility(info,way);
00380 if (value>0.0) {
00381 numberUnsatisfied_++;
00382 if (value>=1e50) {
00383
00384 feasible=false;
00385 break;
00386 }
00387 int priorityLevel = object[i]->priority();
00388
00389 if (priorityLevel<bestPriority) {
00390 for (int j=maximumStrong-1;j>=0;j--) {
00391 if (list_[j]>=0) {
00392 int iObject = list_[j];
00393 list_[j]=-1;
00394 useful_[j]=0.0;
00395 list_[--putOther]=iObject;
00396 }
00397 }
00398 maximumStrong = CoinMin(maximumStrong,putOther);
00399 bestPriority = priorityLevel;
00400 check=-COIN_DBL_MAX;
00401 checkIndex=0;
00402 check2=-COIN_DBL_MAX;
00403 checkIndex2=0;
00404 number_not_trusted_=0;
00405 if (max_most_fra>0) {
00406 for (int j=0;j<max_most_fra;j++) {
00407 list2[j]=-1;
00408 useful2[j]=0.0;
00409 }
00410 }
00411 }
00412 if (priorityLevel==bestPriority) {
00413
00414 double value2;
00415 value = computeUsefulness(MAXMIN_CRITERION,
00416 upMultiplier, downMultiplier, value,
00417 object[i], i, value2);
00418 if (value>check) {
00419
00420 int iObject = list_[checkIndex];
00421 if (iObject>=0) {
00422 assert (list_[putOther-1]<0);
00423 list_[--putOther]=iObject;
00424 }
00425 list_[checkIndex]=i;
00426 assert (checkIndex<putOther);
00427 useful_[checkIndex]=value;
00428
00429 check=COIN_DBL_MAX;
00430 maximumStrong = CoinMin(maximumStrong,putOther);
00431 for (int j=0;j<maximumStrong;j++) {
00432 if (list_[j]>=0) {
00433 if (useful_[j]<check) {
00434 check=useful_[j];
00435 checkIndex=j;
00436 }
00437 }
00438 else {
00439 check=0.0;
00440 checkIndex = j;
00441 break;
00442 }
00443 }
00444 }
00445 else {
00446
00447 assert (list_[putOther-1]<0);
00448 list_[--putOther]=i;
00449 maximumStrong = CoinMin(maximumStrong,putOther);
00450 }
00451 if (max_most_fra > 0 && value2>check2) {
00452
00453 number_not_trusted_++;
00454 list2[checkIndex2]=i;
00455 useful2[checkIndex2]=value2;
00456
00457 check2=COIN_DBL_MAX;
00458 for (int j=0;j<max_most_fra;j++) {
00459 if (list2[j]>=0) {
00460 if (useful2[j]<check2) {
00461 check2=useful2[j];
00462 checkIndex2=j;
00463 }
00464 }
00465 else {
00466 check2=0.0;
00467 checkIndex2 = j;
00468 break;
00469 }
00470 }
00471 }
00472 }
00473 else {
00474
00475
00476 assert (list_[putOther-1]<0);
00477 list_[--putOther]=i;
00478 maximumStrong = CoinMin(maximumStrong,putOther);
00479 }
00480 }
00481 }
00482 #if 0
00483 for (int i=0; i<maximumStrong; i++) {
00484 int way;
00485 message(CANDIDATE_LIST)<<i
00486 <<list_[i]<<i<<useful_[i]<<object[list_[i]]->infeasibility(info,way)
00487 <<CoinMessageEol;
00488 }
00489 #endif
00490
00491 numberOnList_=0;
00492 if (feasible) {
00493 maximumStrong = CoinMin(maximumStrong,putOther);
00494 for (i=0;i<maximumStrong;i++) {
00495 if (list_[i]>=0) {
00496 #ifdef OLD_USEFULLNESS
00497 list_[numberOnList_]=list_[i];
00498 useful_[numberOnList_++]=-useful_[i];
00499
00500 #else
00501 list_[numberOnList_]=list_[i];
00502 if ((sortCrit_ & 1) == 0) {
00503 useful_[numberOnList_++]=-useful_[i];
00504 }
00505 else useful_[numberOnList_++] = useful_[i];
00506 #endif
00507 message(CANDIDATE_LIST2)<<numberOnList_-1
00508 <<list_[numberOnList_-1]<<numberOnList_-1<<useful_[numberOnList_-1]
00509 <<CoinMessageEol;
00510 }
00511 }
00512 if (numberOnList_) {
00513 int tmp_on_list = 0;
00514 if (max_most_fra > 0 && numberOnList_ >= maximumStrong) {
00515
00516
00517 number_not_trusted_=0;
00518 for (i=0;i<max_most_fra;i++) {
00519 if (list2[i]>=0) {
00520 list2[number_not_trusted_] = list2[i];
00521 useful2[number_not_trusted_++] = useful2[i];
00522 message(CANDIDATE_LIST3)<<number_not_trusted_-1
00523 <<list2[number_not_trusted_-1]<<number_not_trusted_-1
00524 <<useful2[number_not_trusted_-1]<<CoinMessageEol;
00525 }
00526 }
00527 if (number_not_trusted_) {
00528 CoinSort_2(list_,list_+numberOnList_,useful_);
00529 CoinSort_2(list2,list2+number_not_trusted_,useful2);
00530 int i1=0;
00531 int i2=0;
00532 for (i=0; i<numberObjects; i++) {
00533 bool found1 = (list_[i1]==i);
00534 bool found2 = (list2[i2]==i);
00535 if (found1 && found2) {
00536 useful_[i1] = -1e150*(1.+useful2[i2]);
00537 list2[i2] = -1;
00538 }
00539 if (found1) i1++;
00540 if (found2) i2++;
00541 if (i2==max_most_fra) break;
00542 }
00543 for (i=0; i<number_not_trusted_; i++) {
00544 if (list2[i] >= 0) {
00545 list_[numberOnList_+tmp_on_list] = list2[i];
00546 useful_[numberOnList_+tmp_on_list] = -1e150*(1.+useful2[i]);
00547 tmp_on_list++;
00548 }
00549 }
00550 }
00551 }
00552
00553 CoinSort_2(useful_,useful_+numberOnList_+tmp_on_list,list_);
00554
00555 i = numberOnList_;
00556 for (;putOther<numberObjects;putOther++)
00557 list_[i++]=list_[putOther];
00558 assert (i==numberUnsatisfied_);
00559 if (!CoinMax(numberStrong_,numberStrongRoot_))
00560 numberOnList_=0;
00561 }
00562 }
00563 else {
00564
00565 numberUnsatisfied_=-1;
00566 }
00567
00568 info->defaultDual_ = -1.0;
00569 delete [] info->usefulRegion_;
00570 delete [] info->indexRegion_;
00571 delete [] list2;
00572 delete [] useful2;
00573 int way;
00574 if (bb_log_level_>3) {
00575
00576 for (int i=0; i<numberOnList_; i++){
00577 message(CANDIDATE_LIST)<<i<< list_[i]<< i<< useful_[i]
00578 <<object[list_[i]]->infeasibility(info,way)
00579 <<CoinMessageEol;
00580 }
00581 }
00582 return numberUnsatisfied_;
00583
00584 }
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598 int
00599 BonChooseVariable::chooseVariable(
00600 OsiSolverInterface * solver,
00601 OsiBranchingInformation *info,
00602 bool fixVariables)
00603 {
00604
00605
00606
00607 bool isRoot = isRootNode(info);
00608 int numberStrong = numberStrong_;
00609 if (isRoot) {
00610 numberStrong = CoinMax(numberStrong_, numberStrongRoot_);
00611 }
00612 std::vector<double> save_sol;
00613 if (bb_log_level_>=3) {
00614 save_sol.resize(info->numberColumns_);
00615 std::copy(info->solution_, info->solution_ + info->numberColumns_ , save_sol.begin());
00616 }
00617 if (numberUnsatisfied_) {
00618 const double* upTotalChange = pseudoCosts_.upTotalChange();
00619 const double* downTotalChange = pseudoCosts_.downTotalChange();
00620 const int* upNumber = pseudoCosts_.upNumber();
00621 const int* downNumber = pseudoCosts_.downNumber();
00622 int numberBeforeTrusted = pseudoCosts_.numberBeforeTrusted();
00623 int numberLeft = CoinMin(numberStrong -numberStrongDone_,numberUnsatisfied_);
00624 results_.clear();
00625 int returnCode=0;
00626 bestObjectIndex_ = -1;
00627 bestWhichWay_ = -1;
00628 firstForcedObjectIndex_ = -1;
00629 firstForcedWhichWay_ =-1;
00630 double bestTrusted=-COIN_DBL_MAX;
00631 for (int i=0;i<numberLeft;i++) {
00632 int iObject = list_[i];
00633 if (numberBeforeTrusted == 0||
00634 i < minNumberStrongBranch_ ||
00635 (
00636 only_pseudo_when_trusted_ && number_not_trusted_>0 ) ||
00637 ( !isRoot && (upNumber[iObject]<numberBeforeTrusted ||
00638 downNumber[iObject]<numberBeforeTrusted ))||
00639 ( isRoot && (!upNumber[iObject] && !downNumber[iObject])) ) {
00640
00641 results_.push_back(HotInfo(solver, info,
00642 solver->objects(), iObject));
00643 }
00644 else {
00645 const OsiObject * obj = solver->object(iObject);
00646 double upEstimate = (upTotalChange[iObject]*obj->upEstimate())/upNumber[iObject];
00647 double downEstimate = (downTotalChange[iObject]*obj->downEstimate())/downNumber[iObject];
00648 double MAXMIN_CRITERION = maxminCrit(info);
00649 double value = MAXMIN_CRITERION*CoinMin(upEstimate,downEstimate) + (1.0-MAXMIN_CRITERION)*CoinMax(upEstimate,downEstimate);
00650 if (value > bestTrusted) {
00651 bestObjectIndex_=iObject;
00652 bestWhichWay_ = upEstimate>downEstimate ? 0 : 1;
00653 bestTrusted = value;
00654 }
00655 }
00656 }
00657 int numberFixed=0;
00658 if (results_.size() > 0) {
00659 returnCode = doStrongBranching(solver, info, (int)results_.size(), 1);
00660 if (bb_log_level_>=3) {
00661 OsiObject ** obj = solver->objects();
00662 const char* stat_msg[] = {"NOTDON", "FEAS", "INFEAS", "NOFINI"};
00663 message(SB_START)<<info->objectiveValue_<<CoinMessageEol;
00664 message(SB_HEADER)<<CoinMessageEol;
00665 for (unsigned int i = 0; i< results_.size(); i++) {
00666 double up_change = results_[i].upChange();
00667 double down_change = results_[i].downChange();
00668 int up_status = results_[i].upStatus();
00669 int down_status = results_[i].downStatus();
00670 int icol = obj[results_[i].whichObject()]->columnNumber();
00671 double val = save_sol[icol];
00672 message(SB_RES)<<(int) i<<icol<<val<<stat_msg[down_status+1]<<down_change
00673 <<stat_msg[up_status+1]<< up_change<< CoinMessageEol;
00674 }
00675 }
00676 if (returnCode>=0&&returnCode<=2) {
00677 if (returnCode) {
00678 returnCode=4;
00679 if (bestObjectIndex_>=0)
00680 returnCode=3;
00681 }
00682 for (unsigned int i=0;i < results_.size();i++) {
00683 int iObject = results_[i].whichObject();
00684 double upEstimate;
00685 if (results_[i].downStatus()== 2 || results_[i].upStatus()==2) {
00686
00687 }
00688 if (results_[i].upStatus()!=1) {
00689 assert (results_[i].upStatus()>=0);
00690 upEstimate = results_[i].upChange();
00691 }
00692 else {
00693
00694 if (info->cutoff_<1.0e50)
00695 upEstimate = 2.0*(info->cutoff_-info->objectiveValue_);
00696 else
00697 upEstimate = 2.0*fabs(info->objectiveValue_);
00698 if (firstForcedObjectIndex_ <0) {
00699
00700 firstForcedObjectIndex_ = iObject;
00701 firstForcedWhichWay_ =0;
00702 }
00703 numberFixed++;
00704 if (fixVariables) {
00705 const OsiObject * obj = solver->object(iObject);
00706 OsiBranchingObject * branch = obj->createBranch(solver,info,0);
00707 branch->branch(solver);
00708 delete branch;
00709 }
00710 }
00711 double downEstimate;
00712 if (results_[i].downStatus()!=1) {
00713 assert (results_[i].downStatus()>=0);
00714 downEstimate = results_[i].downChange();
00715 }
00716 else {
00717
00718 if (info->cutoff_<1.0e50)
00719 downEstimate = 2.0*(info->cutoff_-info->objectiveValue_);
00720 else
00721 downEstimate = 2.0*fabs(info->objectiveValue_);
00722 if (firstForcedObjectIndex_ <0) {
00723 firstForcedObjectIndex_ = iObject;
00724 firstForcedWhichWay_ =1;
00725 }
00726 numberFixed++;
00727 if (fixVariables) {
00728 const OsiObject * obj = solver->object(iObject);
00729 OsiBranchingObject * branch = obj->createBranch(solver,info,1);
00730 branch->branch(solver);
00731 delete branch;
00732 }
00733 }
00734 double MAXMIN_CRITERION = maxminCrit(info);
00735 double value = MAXMIN_CRITERION*CoinMin(upEstimate,downEstimate) + (1.0-MAXMIN_CRITERION)*CoinMax(upEstimate,downEstimate);
00736 if (value>bestTrusted) {
00737 bestTrusted = value;
00738 bestObjectIndex_ = iObject;
00739 bestWhichWay_ = upEstimate>downEstimate ? 0 : 1;
00740
00741 const OsiObject * obj = solver->object(iObject);
00742 if (obj->preferredWay()>=0&&obj->infeasibility())
00743 bestWhichWay_ = obj->preferredWay();
00744 if (returnCode)
00745 returnCode=2;
00746 }
00747 }
00748 }
00749 else if (returnCode==3) {
00750
00751 bestObjectIndex_ = list_[0];
00752 bestWhichWay_ = 0;
00753 returnCode=0;
00754 }
00755 }
00756 else {
00757 bestObjectIndex_=list_[0];
00758 }
00759 if ( bestObjectIndex_ >=0 ) {
00760 OsiObject * obj = solver->objects()[bestObjectIndex_];
00761 obj->setWhichWay(bestWhichWay_);
00762 message(BRANCH_VAR)<<obj->columnNumber()<< bestWhichWay_
00763 <<CoinMessageEol;
00764 }
00765 message(CHOSEN_VAR)<<bestObjectIndex_<<CoinMessageEol;
00766 if (numberFixed==numberUnsatisfied_&&numberFixed)
00767 returnCode=4;
00768 return returnCode;
00769 }
00770 else {
00771 return 1;
00772 }
00773 }
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786 int
00787 BonChooseVariable::doStrongBranching( OsiSolverInterface * solver,
00788 OsiBranchingInformation *info,
00789 int numberToDo, int returnCriterion)
00790 {
00791
00792 double bestLookAhead_ = -COIN_DBL_MAX;
00793 int trialsSinceBest_ = 0;
00794 bool isRoot = isRootNode(info);
00795
00796 double * saveLower = NULL;
00797 double * saveUpper = NULL;
00798 int numberColumns = solver->getNumCols();
00799 solver->markHotStart();
00800 const double * lower = info->lower_;
00801 const double * upper = info->upper_;
00802 saveLower = CoinCopyOfArray(info->lower_,numberColumns);
00803 saveUpper = CoinCopyOfArray(info->upper_,numberColumns);
00804 int returnCode=0;
00805 double timeStart = CoinCpuTime();
00806 int iDo = 0;
00807 for (;iDo<numberToDo;iDo++) {
00808 HotInfo * result = results_() + iDo;
00809
00810 OsiBranchingObject * branch = result->branchingObject();
00811 assert (branch->numberBranches()==2);
00812
00813
00814
00815
00816
00817 OsiSolverInterface * thisSolver = solver;
00818 if (branch->boundBranch()) {
00819
00820 branch->branch(solver);
00821
00822 solver->solveFromHotStart() ;
00823
00824
00825 } else {
00826
00827 thisSolver = solver->clone();
00828 branch->branch(thisSolver);
00829
00830 int limit;
00831 thisSolver->getIntParam(OsiMaxNumIterationHotStart,limit);
00832 thisSolver->setIntParam(OsiMaxNumIteration,limit);
00833 thisSolver->resolve();
00834 }
00835
00836
00837 int status0 = result->updateInformation(thisSolver,info,this);
00838 if (status0==3) {
00839
00840 if (trustStrongForSolution_) {
00841 info->cutoff_ = goodObjectiveValue_;
00842 status0=0;
00843 }
00844 }
00845 if(solver->getRowCutDebugger() && status0 == 1 ){
00846 OsiTMINLPInterface * tminlp_solver = dynamic_cast<OsiTMINLPInterface *> (solver);
00847 throw tminlp_solver->newUnsolvedError(1, tminlp_solver->problem(), "SB");
00848 }
00849 numberStrongIterations_ += thisSolver->getIterationCount();
00850 if (solver!=thisSolver)
00851 delete thisSolver;
00852
00853 for (int j=0;j<numberColumns;j++) {
00854 if (saveLower[j] != lower[j])
00855 solver->setColLower(j,saveLower[j]);
00856 if (saveUpper[j] != upper[j])
00857 solver->setColUpper(j,saveUpper[j]);
00858 }
00859
00860
00861
00862 thisSolver = solver;
00863 if (branch->boundBranch()) {
00864
00865 branch->branch(solver);
00866
00867 fflush(stdout);
00868 solver->solveFromHotStart() ;
00869 } else {
00870
00871 thisSolver = solver->clone();
00872 branch->branch(thisSolver);
00873
00874 int limit;
00875 thisSolver->getIntParam(OsiMaxNumIterationHotStart,limit);
00876 thisSolver->setIntParam(OsiMaxNumIteration,limit);
00877
00878
00879 thisSolver->resolve();
00880 }
00881
00882
00883 int status1 = result->updateInformation(thisSolver,info,this);
00884 numberStrongDone_++;
00885 if (status1==3) {
00886
00887 if (trustStrongForSolution_) {
00888 info->cutoff_ = goodObjectiveValue_;
00889 status1=0;
00890 }
00891 }
00892 if(solver->getRowCutDebugger() && status1 == 1){
00893 OsiTMINLPInterface * tminlp_solver = dynamic_cast<OsiTMINLPInterface *> (solver);
00894 throw tminlp_solver->newUnsolvedError(1, tminlp_solver->problem(), "SB");
00895 }
00896 numberStrongIterations_ += thisSolver->getIterationCount();
00897 if (solver!=thisSolver)
00898 delete thisSolver;
00899
00900 for (int j=0;j<numberColumns;j++) {
00901 if (saveLower[j] != lower[j])
00902 solver->setColLower(j,saveLower[j]);
00903 if (saveUpper[j] != upper[j])
00904 solver->setColUpper(j,saveUpper[j]);
00905 }
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916 if (status0==1&&status1==1) {
00917
00918 returnCode=-1;
00919
00920 } else if (status0==1||status1==1) {
00921 numberStrongFixed_++;
00922 returnCode=1;
00923 }
00924 bool hitMaxTime = ( CoinCpuTime()-timeStart > info->timeRemaining_)
00925 || ( CoinCpuTime() - start_time_ > time_limit_);
00926 if (hitMaxTime) {
00927 returnCode=3;
00928 break;
00929 }
00930
00931 if (!isRoot && numberLookAhead_) {
00932 assert(status0==0 && status1==0);
00933 double upEstimate = result->upChange();
00934 double downEstimate = result->downChange();
00935 double MAXMIN_CRITERION = maxminCrit(info);
00936 double value = MAXMIN_CRITERION*CoinMin(upEstimate,downEstimate) + (1.0-MAXMIN_CRITERION)*CoinMax(upEstimate,downEstimate);
00937 if (value > bestLookAhead_) {
00938 bestLookAhead_ = value;
00939 trialsSinceBest_ = 0;
00940 }
00941 else {
00942 trialsSinceBest_++;
00943 if (trialsSinceBest_ >= numberLookAhead_) {
00944 break;
00945 }
00946 }
00947 }
00948 }
00949 if(iDo < numberToDo) iDo++;
00950 assert(iDo <= (int) results_.size());
00951 results_.resize(iDo);
00952 delete [] saveLower;
00953 delete [] saveUpper;
00954
00955 solver->unmarkHotStart();
00956 return returnCode;
00957 }
00958
00959 bool BonChooseVariable::isRootNode(const OsiBranchingInformation *info) const
00960 {
00961 return info->depth_ == 0;
00962 }
00963
00964 double
00965 BonChooseVariable::maxminCrit(const OsiBranchingInformation *info) const
00966 {
00967 double retval = maxmin_crit_no_sol_;
00968 if (cbc_model_) {
00969
00970 const int stateOfSearch = cbc_model_->stateOfSearch();
00971 const int depth = info->depth_;
00972 if (stateOfSearch>1 && depth>10) {
00973 retval = maxmin_crit_have_sol_;
00974 }
00975 }
00976 return retval;
00977 }
00978
00979
00980 void
00981 BonChooseVariable::updateInformation(const OsiBranchingInformation *info,
00982 int branch, OsiHotInfo * hotInfo)
00983 {
00984 if(!trustStrongForPseudoCosts_) return;
00985 int index = hotInfo->whichObject();
00986 assert (index<solver_->numberObjects());
00987 const OsiObject * object = info->solver_->object(index);
00988 assert (object->upEstimate()>0.0&&object->downEstimate()>0.0);
00989 assert (branch<2);
00990 double* upTotalChange = pseudoCosts_.upTotalChange();
00991 double* downTotalChange = pseudoCosts_.downTotalChange();
00992 int* upNumber = pseudoCosts_.upNumber();
00993 int* downNumber = pseudoCosts_.downNumber();
00994 if (branch) {
00995
00996
00997
00998 if (hotInfo->upStatus()==0) {
00999 assert (hotInfo->upStatus()>=0);
01000 upTotalChange[index] += hotInfo->upChange()/object->upEstimate();
01001 upNumber[index]++;
01002 }
01003 else if (hotInfo->upStatus()==1) {
01004
01005 upNumber[index]++;
01006 if (info->cutoff_<1.0e50)
01007 upTotalChange[index] += 2.0*(info->cutoff_-info->objectiveValue_)/object->upEstimate();
01008 else
01009 upTotalChange[index] += 2.0*fabs(info->objectiveValue_)/object->upEstimate();
01010 }
01011 }
01012 else {
01013 if (hotInfo->downStatus()==0) {
01014 assert (hotInfo->downStatus()>=0);
01015 downTotalChange[index] += hotInfo->downChange()/object->downEstimate();
01016 downNumber[index]++;
01017 }
01018 else if (hotInfo->downStatus()==1) {
01019 downNumber[index]++;
01020
01021 if (info->cutoff_<1.0e50)
01022 downTotalChange[index] += 2.0*(info->cutoff_-info->objectiveValue_)/object->downEstimate();
01023 else
01024 downTotalChange[index] += 2.0*fabs(info->objectiveValue_)/object->downEstimate();
01025 }
01026 }
01027 }
01028
01029
01030 void
01031 BonChooseVariable::updateInformation( int index, int branch,
01032 double changeInObjective, double changeInValue,
01033 int status)
01034 {
01035 if(cbc_model_ == NULL) return;
01036 assert (index<solver_->numberObjects());
01037 assert (branch<2);
01038
01039 if(fabs(changeInValue) < 1e-6) return;
01040
01041 double* upTotalChange = pseudoCosts_.upTotalChange();
01042 double* downTotalChange = pseudoCosts_.downTotalChange();
01043 int* upNumber = pseudoCosts_.upNumber();
01044 int* downNumber = pseudoCosts_.downNumber();
01045 message(UPDATE_PS_COST)<<index<< branch
01046 <<changeInObjective<<changeInValue<<status
01047 <<CoinMessageEol;
01048
01049 if (branch) {
01050 if (status!=1) {
01051 assert (status>=0);
01052 upTotalChange[index] += changeInObjective/changeInValue;
01053 upNumber[index]++;
01054 } else {
01055
01056 upNumber[index]++;
01057 assert(cbc_model_);
01058 double cutoff = cbc_model_->getCutoff();
01059 double objectiveValue = cbc_model_->getCurrentObjValue();
01060 if (cutoff<1.0e50)
01061 upTotalChange[index] += 2.0*(cutoff-objectiveValue)/changeInValue;
01062 else
01063 upTotalChange[index] += 2.0*fabs(objectiveValue)/changeInValue;
01064 }
01065 } else {
01066 if (status!=1) {
01067 assert (status>=0);
01068 downTotalChange[index] += changeInObjective/changeInValue;
01069 downNumber[index]++;
01070 } else {
01071 assert(cbc_model_);
01072
01073 downNumber[index]++;
01074 double cutoff = cbc_model_->getCutoff();
01075 double objectiveValue = cbc_model_->getCurrentObjValue();
01076 if (cutoff<1.0e50)
01077 downTotalChange[index] += 2.0*(cutoff-objectiveValue)/changeInValue;
01078 else
01079 downTotalChange[index] += 2.0*fabs(objectiveValue)/changeInValue;
01080 }
01081 }
01082 }
01083
01084
01085 HotInfo::HotInfo(): OsiHotInfo(), infeasibilities_(){
01086 }
01087
01088 HotInfo::HotInfo( OsiSolverInterface * solver,
01089 const OsiBranchingInformation *info,
01090 const OsiObject * const * objects,
01091 int whichObject):
01092 OsiHotInfo(solver, info, objects, whichObject),
01093 infeasibilities_(){
01094 infeasibilities_.resize(branchingObject_->numberBranches());
01095 }
01096
01097 HotInfo::HotInfo(const HotInfo& other): OsiHotInfo(other),
01098 infeasibilities_(other.infeasibilities_){
01099 }
01100
01101 OsiHotInfo *
01102 HotInfo::clone() const {
01103 return new HotInfo(*this);
01104 }
01105
01106 HotInfo &
01107 HotInfo::operator=(const HotInfo &rhs){
01108 if(this != &rhs){
01109 OsiHotInfo::operator=(rhs);
01110 infeasibilities_ = rhs.infeasibilities_;
01111 }
01112 return (*this);
01113 }
01114
01115 HotInfo::~HotInfo(){
01116 }
01117
01118 int
01119 HotInfo::updateInformation(const OsiSolverInterface * solver,
01120 const OsiBranchingInformation * info,
01121 OsiChooseVariable * choose){
01122
01123 int iBranch = branchingObject_->branchIndex()-1;
01124 double & infeasibility = infeasibilities_[iBranch] = 0.;
01125
01126 OsiObject ** objects = solver->objects();
01127 int numObject = solver->numberObjects();
01128 for(int i = 0 ; i < numObject ; i++){
01129 infeasibility += objects[i]->checkInfeasibility(info);
01130 }
01131 int status = OsiHotInfo::updateInformation(solver, info, choose);
01132 #if 1
01133 if(!solver->isProvenPrimalInfeasible() && !solver->isProvenOptimal()){
01134 status = 2;
01135 statuses_[iBranch] = status;
01136 }
01137 else if(solver->isProvenPrimalInfeasible() && fabs(solver->getObjValue()) < 1e-06) {
01138 assert(solver->messageHandler() != NULL);
01139 *solver->messageHandler() << "Very small infeasibility: " << solver->getObjValue() << CoinMessageEol;
01140 status = 2;
01141 statuses_[iBranch] = status;
01142 }
01143 #endif
01144 return status;
01145 }
01146
01147 }
01148