00001
00002
00003
00004 #include <climits>
00005 #include "CoinPragma.hpp"
00006 #include "BonNWayChoose.hpp"
00007 #include "BonNWayObject.hpp"
00008 #include "CoinTime.hpp"
00009
00010 #ifndef NDEBUG
00011 #define ASSERTED_CAST static_cast
00012 #else
00013 #define ASSERTED_CAST dynamic_cast
00014 #endif
00015 namespace Bonmin
00016 {
00017
00018
00019 BonNWayChoose::BonNWayChoose(BabSetupBase &b, const OsiSolverInterface* solver):
00020 OsiChooseVariable(solver),
00021 br_depth_(0),
00022 bounds_(),
00023 unit_changes_(),
00024 num_ps_costs_(),
00025 num_eval_(),
00026 geo_means_(0)
00027 {
00028 Ipopt::SmartPtr<Ipopt::OptionsList> options = b.options();
00029 options->GetNumericValue("time_limit", time_limit_, b.prefix());
00030 options->GetNumericValue("cutoff_multiplier", cutoff_multiplier_, b.prefix());
00031 options->GetNumericValue("pseudocost_trust_value", pseudocost_trust_value_, b.prefix());
00032 options->GetIntegerValue("strong_branch_depth", br_depth_, b.prefix());
00033 options->GetIntegerValue("nway_branch_log_level", log_, b.prefix());
00034 options->GetEnumValue("do_fixings", do_fixings_, b.prefix());
00035 options->GetEnumValue("use_geo_means", geo_means_, b.prefix());
00037 int numberObjects = solver_->numberObjects();
00038 std::cout<<"Number objects "<<numberObjects<<std::endl;
00039 start_time_ = CoinCpuTime();
00040 OsiObject ** object = solver->objects();
00041 for (int i=0;i<numberObjects;i++) {
00042 BonNWayObject * nway = dynamic_cast<BonNWayObject *>(object[i]);
00043 if(!nway) continue;
00044 start_nway_ = i;
00045 break;
00046 }
00047 numberObjects -= start_nway_;
00048 }
00049
00050 BonNWayChoose::BonNWayChoose(const BonNWayChoose & rhs) :
00051 OsiChooseVariable(rhs),
00052 br_depth_(rhs.br_depth_),
00053 do_fixings_(rhs.do_fixings_),
00054 cutoff_multiplier_(rhs.cutoff_multiplier_),
00055 pseudocost_trust_value_(rhs.pseudocost_trust_value_),
00056 time_limit_(rhs.time_limit_),
00057 start_time_(rhs.start_time_),
00058 start_nway_(rhs.start_nway_),
00059 log_(rhs.log_),
00060 bounds_(rhs.bounds_),
00061 unit_changes_(rhs.unit_changes_),
00062 num_ps_costs_(rhs.num_ps_costs_),
00063 num_eval_(rhs.num_eval_),
00064 geo_means_(rhs.geo_means_)
00065 {
00066 }
00067
00068 BonNWayChoose &
00069 BonNWayChoose::operator=(const BonNWayChoose & rhs)
00070 {
00071 if (this != &rhs) {
00072 br_depth_ = rhs.br_depth_;
00073 do_fixings_ = rhs.do_fixings_;
00074 cutoff_multiplier_ = rhs.cutoff_multiplier_;
00075 pseudocost_trust_value_ = rhs.pseudocost_trust_value_;
00076 time_limit_ = rhs.time_limit_;
00077 start_time_ = rhs.start_time_;
00078 log_ = rhs.log_;
00079 start_nway_ = rhs.start_nway_;
00080 OsiChooseVariable::operator=(rhs);
00081 bounds_ = rhs.bounds_;
00082 unit_changes_ = rhs.unit_changes_;
00083 num_ps_costs_ = rhs.num_ps_costs_;
00084 num_eval_ = rhs.num_eval_;
00085 geo_means_ = rhs.geo_means_;
00086 }
00087 return *this;
00088 }
00089
00090 OsiChooseVariable *
00091 BonNWayChoose::clone() const
00092 {
00093 return new BonNWayChoose(*this);
00094 }
00095
00096 BonNWayChoose::~BonNWayChoose ()
00097 {
00098 }
00099
00100 void
00101 BonNWayChoose::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)
00102 {
00103 roptions->SetRegisteringCategory("NWay Strong branching setup", RegisteredOptions::BonminCategory);
00104 roptions->AddLowerBoundedIntegerOption("nway_branch_log_level",
00105 "Log level for the branching on nways",
00106 0,1,
00107 "");
00108
00109 roptions->AddLowerBoundedIntegerOption("strong_branch_depth",
00110 "To which level do we perform strong-branching",
00111 0,0,
00112 "");
00113
00114 roptions->AddLowerBoundedNumberOption("cutoff_multiplier",
00115 "multiplier applied to cutoff_ for computing pseudo-cost of infeasible sub-problems",
00116 1.,0,3.,
00117 "");
00118
00119 roptions->AddLowerBoundedNumberOption("pseudocost_trust_value",
00120 "Trust pseudo cost of best nway object if it is above this value",
00121 0.,0,0,
00122 "");
00123
00124 roptions->AddStringOption2("use_geo_means", "Use geometrical means to average pseudo-costs",
00125 "yes",
00126 "no", "Use artihmetical means",
00127 "yes", "Use geometrical means","");
00128
00129 roptions->AddStringOption4("do_fixings",
00130 "Do we fix variables in strong branching?",
00131 "all",
00132 "none", "Don't do any.",
00133 "in-tree", "Fix only variables in the tree",
00134 "strong-branching", "Fix variable in strong branching only",
00135 "all", "Fix whenever possible",
00136 "");
00137
00138
00139 }
00140
00141 double
00142 BonNWayChoose::compute_usefulness(int objectIndex, const OsiBranchingInformation * info) const
00143 {
00144 int nwayIndex = objectIndex - start_nway_;
00145
00146 BonNWayObject * nway = ASSERTED_CAST<BonNWayObject *>(info->solver_->objects()[objectIndex]);
00147 assert(nway);
00148 size_t n = nway->numberMembers();
00149 const int * vars = nway->members();
00150
00151 std::vector<double> unit_changes(unit_changes_[nwayIndex]);
00152 if(geo_means_){
00153 for(size_t k = 0 ; k < unit_changes.size() ; k++) unit_changes[k] = (num_ps_costs_[nwayIndex][k]) ? pow(unit_changes[k], 1./(double) num_ps_costs_[nwayIndex][k]): unit_changes[k];
00154 }
00155 else {
00156 for(size_t k = 0 ; k < unit_changes.size() ; k++) unit_changes[k] = (num_ps_costs_[nwayIndex][k]) ? unit_changes[k]/(double) num_ps_costs_[nwayIndex][k]: unit_changes[k];
00157 }
00158
00159
00160 double r_val = compute_usefulness(info, n, vars, bounds_[nwayIndex], unit_changes);
00161 return r_val;
00162 }
00163
00164 double
00165 BonNWayChoose::compute_usefulness(const OsiBranchingInformation * info,
00166 size_t n, const int * vars,const std::vector<double> &bounds, const std::vector<double> &unit_changes) const
00167 {
00168 const double * solution = info->solution_;
00169 double obj_val = info->objectiveValue_;
00170 const double * lower = info->lower_;
00171 const double * upper = info->upper_;
00172 double integerTolerance = info->integerTolerance_;
00173 double cutoff = info->cutoff_*cutoff_multiplier_;
00174 double r_val = (geo_means_) ? 1 : 0;
00175
00176 for(size_t i = 0 ; i < n ; i++){
00177 int iCol = vars[i];
00178 if(fabs(lower[iCol] - upper[iCol]) < integerTolerance) {
00179
00180 continue;
00181 }
00182 assert(lower[iCol] < upper[iCol]);
00183 double residual = upper[iCol] - solution[iCol];
00184 double score = std::min(cutoff - obj_val, std::max(residual*unit_changes[i],bounds[i] - obj_val));
00185 if(geo_means_)
00186 r_val*=score;
00187 else
00188 r_val += score;
00189 }
00190 return r_val;
00191 }
00192
00193 int
00194 BonNWayChoose::setupList ( OsiBranchingInformation *info, bool initialize)
00195 {
00196 assert(initialize);
00197 if (initialize) {
00198 status_=-2;
00199 delete [] goodSolution_;
00200 bestObjectIndex_=-1;
00201 goodSolution_ = NULL;
00202 goodObjectiveValue_ = COIN_DBL_MAX;
00203 }
00204
00205
00206 numberOnList_=0;
00207 numberUnsatisfied_=0;
00208 int numberObjects = info->solver_->numberObjects() - start_nway_;
00209 double cutoff = info->cutoff_;
00210 if(info->depth_ == 0){
00211 bounds_.resize(0);
00212 unit_changes_.resize(0);
00213 bounds_.resize(numberObjects, std::vector<double>(numberObjects,cutoff));
00214 unit_changes_.resize(numberObjects, std::vector<double>(numberObjects,0));
00215 num_eval_.resize(numberObjects, 0);
00216 num_ps_costs_.resize(numberObjects, std::vector<int>(numberObjects,0));
00217 }
00218 else {
00219 assert(unit_changes_.size() == numberObjects);
00220 assert(bounds_.size() == unit_changes_.size());
00221 }
00222
00223 int maximumStrong = numberObjects;
00224
00225 double check = -COIN_DBL_MAX;
00226 int checkIndex=0;
00227 int bestPriority=COIN_INT_MAX;
00228 int putOther = numberObjects;
00229 int i;
00230 for (i=0;i<numberObjects;i++) {
00231 list_[i]=-1;
00232 useful_[i]=0.0;
00233 }
00234
00235 OsiObject ** object = info->solver_->objects();
00236 object += start_nway_;
00237
00238
00239 bool feasible = true;
00240 for ( i=0;i<numberObjects;i++) {
00241 int way;
00242 double value = object[i]->infeasibility(info,way);
00243 if (value>0.0) {
00244 numberUnsatisfied_++;
00245 int priorityLevel = object[i]->priority();
00246
00247 if (priorityLevel<bestPriority) {
00248 for (int j=maximumStrong-1;j>=0;j--) {
00249 if (list_[j]>=0) {
00250 int iObject = list_[j];
00251 list_[j]=-1;
00252 useful_[j]=0.0;
00253 list_[--putOther]=iObject;
00254 }
00255 }
00256 bestPriority = priorityLevel;
00257 check=-COIN_DBL_MAX;
00258 checkIndex=0;
00259 }
00260 if (priorityLevel==bestPriority) {
00261
00262
00263 if(info->depth_ != 0){
00264 value = compute_usefulness(i + start_nway_, info);
00265 }
00266
00267 if (value>check) {
00268
00269 int iObject = list_[checkIndex];
00270 if (iObject>=0) {
00271 assert (list_[putOther-1]<0);
00272 list_[--putOther]=iObject;
00273 }
00274 list_[checkIndex]= i + start_nway_;
00275 assert (checkIndex<putOther);
00276 useful_[checkIndex]=value;
00277
00278 check=COIN_DBL_MAX;
00279 maximumStrong = CoinMin(maximumStrong,putOther);
00280 for (int j=0;j<maximumStrong;j++) {
00281 if (list_[j]>=0) {
00282 if (useful_[j]<check) {
00283 check=useful_[j];
00284 checkIndex=j;
00285 }
00286 }
00287 else {
00288 check=0.0;
00289 checkIndex = j;
00290 break;
00291 }
00292 }
00293 }
00294 else {
00295
00296 assert (list_[putOther-1]<0);
00297 list_[--putOther]=i + start_nway_;
00298 maximumStrong = CoinMin(maximumStrong,putOther);
00299 }
00300 }
00301 else {
00302
00303
00304 assert (list_[putOther-1]<0);
00305 list_[--putOther]=i + start_nway_;
00306 maximumStrong = CoinMin(maximumStrong,putOther);
00307 }
00308 }
00309 }
00310
00311
00312 numberOnList_=0;
00313 if (feasible) {
00314 maximumStrong = CoinMin(maximumStrong,putOther);
00315 for (i=0;i<maximumStrong;i++) {
00316 if (list_[i]>=0) {
00317 list_[numberOnList_]=list_[i];
00318 useful_[numberOnList_++]=-useful_[i];
00319 }
00320 }
00321 if (numberOnList_) {
00322
00323 CoinSort_2(useful_,useful_+numberOnList_,list_);
00324
00325 i = numberOnList_;
00326 for (;putOther<numberObjects;putOther++)
00327 list_[i++]=list_[putOther];
00328 assert (i==numberUnsatisfied_);
00329 if (!numberStrong_)
00330 numberOnList_=0;
00331 }
00332 }
00333 else {
00334
00335 numberUnsatisfied_=-1;
00336 }
00337
00338 info->defaultDual_ = -1.0;
00339 delete [] info->usefulRegion_;
00340 delete [] info->indexRegion_;
00341
00342 return numberUnsatisfied_;
00343 }
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357 int
00358 BonNWayChoose::chooseVariable(
00359 OsiSolverInterface * solver,
00360 OsiBranchingInformation *info,
00361 bool fixVariables)
00362 {
00363 if(!numberUnsatisfied_) return 1;
00364
00365
00366 double cutoff = info->cutoff_;
00367 double obj_val = info->objectiveValue_;
00368 if(info->depth_ > br_depth_ && ( (- useful_[0] > pseudocost_trust_value_ )|| num_eval_[list_[0] - start_nway_] >= 18) ){
00369 const double * lower = info->lower_;
00370 const double * upper = info->upper_;
00371
00372
00373
00374 int n_fixed = 0;
00375 if(do_fixings_ > 1){
00376 for(int i = 0 ; i < numberUnsatisfied_ ; i++){
00377 int iObject = list_[i];
00378 int nwayIndex = iObject - start_nway_;
00379 const BonNWayObject * nway = ASSERTED_CAST<const BonNWayObject *>(solver->object(iObject));
00380
00381 size_t n = nway->numberMembers();
00382 const int * vars = nway->members();
00383 for(size_t j = 0 ; j < n ; j++){
00384 int iCol = vars[j];
00385 if(upper[iCol] < lower[iCol] + 0.5) continue;
00386 if(bounds_[nwayIndex][j] > cutoff){
00387 solver->setColUpper(iCol, lower[iCol]);
00388 n_fixed ++;
00389 }
00390 }
00391 }
00392 if(n_fixed && log_ > 1)
00393 printf("NWAY: Fixed %i variables\n", n_fixed);
00394 }
00395
00396 assert(bounds_.size() == unit_changes_.size());
00397 assert(unit_changes_.size() == info->solver_->numberObjects() - start_nway_);
00398
00399 bestObjectIndex_ = list_[0];
00400 bestWhichWay_ = 1;
00401 OsiObject * obj = solver->objects()[bestObjectIndex_];
00402 obj->setWhichWay(bestWhichWay_);
00403
00404 if(log_ > 1){
00405 printf("level %i branch on %i bound %g usefullness %g.\n",
00406 info->depth_, bestObjectIndex_ - start_nway_, obj_val, - useful_[0]);
00407 }
00408 if(n_fixed) return 2;
00409 return 0;
00410 }
00411
00412
00413 if(log_ > 0)
00414 printf("Restarting strong branching loop....\n\n");
00415
00416 numberStrongIterations_ = 0;
00417 numberStrongDone_ = 0;
00418 int numberLeft = numberOnList_;
00419 int returnCode=0;
00420 bestObjectIndex_ = -1;
00421 bestWhichWay_ = -1;
00422 firstForcedObjectIndex_ = -1;
00423 firstForcedWhichWay_ =-1;
00424 double best_score = -COIN_DBL_MAX;
00425 int bestPriority=0;
00426
00427 int n = solver->getNumCols();
00428 std::vector<double> saveLower(n);
00429 std::vector<double> saveUpper(n);
00430 std::copy(info->lower_, info->lower_ + n, saveLower.begin());
00431 std::copy(info->upper_, info->upper_ + n, saveUpper.begin());
00432
00433
00434 solver->markHotStart();
00435 for (int i=0;i<numberLeft;i++) {
00436 int iObject = list_[i];
00437 const int objectPriority = solver->object(iObject)->priority();
00438 if (objectPriority >= bestPriority){
00439 bestPriority = objectPriority;
00440 }
00441 else break;
00442 double score;
00443 int r_val = doStrongBranching(solver, info, iObject, saveLower.data(),
00444 saveUpper.data(), score);
00445 if(r_val == -1) {
00446 if(log_ > 0)
00447 std::cout<<"This is Infeasible"<<std::endl;
00448 returnCode = -1;
00449 break;
00450 }
00451
00452 if(do_fixings_ > 1 && r_val == 1 && info->depth_ == 0) returnCode=2;
00453
00454 if(log_ > 0)
00455 printf("Usefullness from strong branching on %i : %g\n", iObject - start_nway_, score);
00456 if(score > best_score){
00457 best_score = score;
00458 bestObjectIndex_ = iObject;
00459 bestWhichWay_ = 0;
00460 }
00461 if (r_val==3) {
00462 returnCode = 3;
00463
00464 if(bestObjectIndex_ < 0){
00465 bestObjectIndex_ = list_[0];
00466 bestWhichWay_ = 0;
00467 }
00468 break;
00469 }
00470 }
00471 solver->unmarkHotStart();
00472 return returnCode;
00473 }
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484 int
00485 BonNWayChoose::doStrongBranching( OsiSolverInterface * solver,
00486 OsiBranchingInformation *info,
00487 int objectIndex,
00488 double * saveLower,
00489 double * saveUpper, double & score)
00490 {
00491 int nwayIndex = objectIndex - start_nway_;
00492 const double * lower = info->lower_;
00493 const double * upper = info->upper_;
00494
00495 int numberColumns = solver->getNumCols();
00496 double timeStart = CoinCpuTime();
00497
00498 int numberObjects = info->solver_->numberObjects();
00499 const BonNWayObject * nway = ASSERTED_CAST<const BonNWayObject *>(solver->object(objectIndex));
00500 assert(nway);
00501 BonNWayBranchingObject * branch = ASSERTED_CAST<BonNWayBranchingObject *>(nway->createBranch(solver, info, 1));
00502
00503 int branches_left = branch->numberBranchesLeft();
00504 int number_branches = branch->numberBranchesLeft();
00505 int n_can_be_fixed = 0;
00506
00507 double big_val = cutoff_multiplier_*info->cutoff_;
00508 if(big_val > 1e10){ big_val = 10*info->objectiveValue_;}
00509 big_val += fabs(big_val)*1e-5;
00510 std::vector<double> unit_changes(numberObjects - start_nway_, -DBL_MAX);
00511
00512 while(branches_left){
00513
00514 branch->branch(solver);
00515 int v_br = branch->var_branched_on();
00516 int s_br = branch->seq_branched_on();
00517 double residual = upper[v_br] - info->solution_[v_br];
00518 solver->solveFromHotStart() ;
00519 numberStrongIterations_ += solver->getIterationCount();
00520 numberStrongDone_++;
00521
00522 double obj_val = solver->getObjValue();
00523
00524 if(solver->isProvenPrimalInfeasible() ||
00525 (solver->isProvenOptimal() && obj_val > info->cutoff_)){
00526 if(info->depth_ == 0){
00527 bounds_[nwayIndex][s_br] = big_val;
00528 }
00529 unit_changes[s_br] = (big_val - info->objectiveValue_)/residual;
00530 if(do_fixings_ > 1){
00531 n_can_be_fixed++;
00532 if(log_ > 0)
00533 printf("Fixing variable %i to 0 the cutoff is %g\n", v_br, big_val);
00534 saveUpper[v_br] = saveLower[v_br];
00535 }
00536 }
00537 else{
00538 if(info->depth_ == 0){
00539 bounds_[nwayIndex][s_br] = obj_val;
00540 }
00541 unit_changes[s_br] = (obj_val - info->objectiveValue_)/residual;
00542 }
00543
00544
00545 for (int j=0;j<numberColumns;j++) {
00546 if (saveLower[j] != lower[j])
00547 solver->setColLower(j,saveLower[j]);
00548 if (saveUpper[j] != upper[j])
00549 solver->setColUpper(j,saveUpper[j]);
00550 }
00551 branches_left = branch->numberBranchesLeft();
00552 }
00553
00554 score = compute_usefulness(info, nway->numberMembers(), nway->members(), bounds_[nwayIndex], unit_changes);
00555 if(info->depth_ == 0){
00556 if(do_fixings_ == 1 || do_fixings_ == 3)
00557 nway->set_bounds(bounds_[nwayIndex]);
00558 for(size_t k = 0 ; k < unit_changes.size() ; k++){
00559 num_ps_costs_[nwayIndex][k]=1;
00560 }
00561 unit_changes_[nwayIndex] = unit_changes;
00562 num_eval_[nwayIndex] = 1;
00563 }
00564 else if (n_can_be_fixed < number_branches -1){
00565 num_eval_[nwayIndex]++;
00566 for(size_t k = 0 ; k < unit_changes.size() ; k++){
00567 if(unit_changes[k] > 0.){
00568 if(geo_means_)
00569 unit_changes_[nwayIndex][k] *= unit_changes[k];
00570 else
00571 unit_changes_[nwayIndex][k] += unit_changes[k];
00572 num_ps_costs_[nwayIndex][k]++;
00573 }
00574 }
00575 }
00576 if(n_can_be_fixed == number_branches){
00577 return -1;
00578 }
00579 if(n_can_be_fixed){
00580 return 1;
00581 }
00582 bool hitMaxTime = ( CoinCpuTime()-timeStart > info->timeRemaining_)
00583 || ( CoinCpuTime() - start_time_ > time_limit_);
00584 if (hitMaxTime) {
00585 return 3;
00586 }
00587 return 0;
00588 }
00589
00590 }
00591