/home/coin/SVN-release/OS-2.4.2/Couenne/src/branch/doStrongBranching.cpp

Go to the documentation of this file.
00001 /* $Id: doStrongBranching.cpp 726 2011-06-30 19:42:37Z fmargot $
00002  *
00003  * Name:    doStrongBranching.cpp
00004  * Authors: Andreas Waechter, IBM Corp.
00005  *          Pietro Belotti, CMU
00006  * Purpose: actual strong branching method
00007  *
00008  * (C) Carnegie-Mellon University, 2008-10.
00009  * This file is licensed under the Eclipse Public License (EPL)
00010  */
00011 
00012 #include "CoinTime.hpp"
00013 #include "BonChooseVariable.hpp"
00014 #include "CouenneChooseStrong.hpp"
00015 #include "CouenneProblem.hpp"
00016 #include "CouenneObject.hpp"
00017 #include "CouenneBranchingObject.hpp"
00018 
00019 //#define TRACE_STRONG
00020 //#define TRACE_STRONG2
00021 
00022 using namespace Ipopt;
00023 using namespace Couenne;
00024 
00027 bool BranchingFBBT (CouenneProblem *problem,
00028                     OsiObject *Object,
00029                     OsiSolverInterface *solver);
00030 
00033 double distance (const double *p1, const double *p2, int size, double k=2.) {
00034 
00035   double 
00036     result = 0.,
00037     element;
00038 
00039   if (k == 2.) // a bit faster, probably
00040 
00041     while (size--) {
00042       element = *p1++ - *p2++;
00043       result += element * element;
00044     }
00045 
00046   else
00047 
00048     while (size--) {
00049       element = *p1++ - *p2++;
00050       result += pow (element, k);
00051     }
00052 
00053   return pow (result, 1. / k);
00054 }
00055 
00056 
00069   int CouenneChooseStrong::doStrongBranching (OsiSolverInterface *solver, 
00070                                               OsiBranchingInformation *info,
00071                                               int numberToDo, int returnCriterion) {
00072 
00073 #ifdef TRACE_STRONG2
00074     int pv = -1;
00075     if(problem_->doPrint_) {
00076       if(pv > -1) {
00077         printf("doSB: beg: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
00078                pv, solver->getColSolution()[pv], solver->getColLower()[pv], solver->getColUpper()[pv]);
00079         printf("doSB: info: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
00080                pv, info->solution_[pv], info->lower_[pv], info->upper_[pv]);
00081         printf("doSB: problem: lb: %10.4f  ub: %10.4f\n",
00082                problem_->Lb(pv), problem_->Ub(pv));
00083       }
00084     }
00085 #endif
00086 
00087     jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING, 
00088                       "\n-\n------- CCS: trying %d objects:\n", numberToDo);
00089 
00090     //solver -> doingResolve () = false; // turns off setCutoff and restoreUnused
00091 
00092     int numberColumns = solver -> getNumCols ();
00093 
00094     solver -> markHotStart (); // save current LP point
00095 
00096     // save initial bounds
00097     const double
00098       *initLower = info -> lower_,
00099       *initUpper = info -> upper_;
00100 
00101     // save intersection of bounds obtained by branching on all objects
00102     double 
00103       *saveLower = CoinCopyOfArray (info -> lower_, numberColumns),
00104       *saveUpper = CoinCopyOfArray (info -> upper_, numberColumns),
00105 
00106       *Lower0 = NULL,
00107       *Upper0 = NULL,
00108 
00109       // save union of bounds on both branches of one object;
00110       // reset to (current) saveLower when branching on a new object
00111       *unionLower  = new double [numberColumns],
00112       *unionUpper  = new double [numberColumns],
00113 
00114       *lpSol     = NULL, 
00115        timeStart = CoinCpuTime ();
00116 
00117     if (jnlst_ -> ProduceOutput (J_DETAILED, J_BRANCHING)) {
00118       Lower0 = CoinCopyOfArray (info -> lower_, numberColumns); // delete afterwards
00119       Upper0 = CoinCopyOfArray (info -> upper_, numberColumns);
00120     }
00121 
00122     // LP solution for distance
00123     if (pseudoUpdateLP_) 
00124       lpSol = CoinCopyOfArray (info -> solution_, numberColumns);
00125 
00126     // provide Couenne problem with point/bounds contained in info
00127     // problem_ -> domain () -> push
00128     //   (problem_ -> nVars (),
00129     //    info -> solution_,
00130     //    info -> lower_,
00131     //    info -> upper_);
00132 
00133     //Bonmin::HotInfo * results = results_ ();
00134 
00135     int returnCode = 0, iDo = 0;
00136 
00137     for (iDo = 0; iDo < numberToDo; iDo++) {
00138 
00139       Bonmin::HotInfo * result = results_ () + iDo; // retrieve i-th object to test
00140 
00141       OsiObject *Object = solver_ -> objects () [result -> whichObject ()];
00142 
00143       // TODO: apply isCuttable()     
00144 
00145       // TODO: set a cutoff for dual bound in dual simplex
00146       //       do the same for primal based on SB's alpha
00147 
00148       // For now just 2 way
00149       OsiBranchingObject * branch = result -> branchingObject ();
00150       assert (branch->numberBranches()==2);
00151 
00152       CouenneBranchingObject *cb = dynamic_cast <CouenneBranchingObject *> (branch);
00153 
00154       if (cb) cb -> setSimulate (true);
00155 
00156       int 
00157         status0 = -1, 
00158         status1 = -1;
00159 
00161 
00162       /* Try the first direction.  Each subsequent call to branch()
00163          performs the specified branch and advances the branch object
00164          state to the next branch alternative. */
00165 
00166       bool isInf0 = false;
00167       bool isInf1 = false;
00168 
00169       double indUb = 0, indLb = 0;
00170       CouenneObject *CouObj = dynamic_cast <CouenneObject *> (Object);
00171       OsiSimpleInteger *simpl = dynamic_cast <OsiSimpleInteger *>(solver_->objects()[result->whichObject ()]);
00172       // if OsiSimpleInteger Object with branching point outside 
00173       // current solver bound interval, one branch must
00174       // be set as infeasible, otherwise bounds are enlarged 
00175       // in one branch
00176       int objVarIndex = -1;
00177       if(CouObj) {
00178         objVarIndex = CouObj->Reference()->Index();
00179       }
00180       else {
00181         objVarIndex = Object->columnNumber();
00182       }
00183       if(simpl) { 
00184         if(objVarIndex >= 0) {
00185           indUb = solver->getColUpper()[objVarIndex];
00186           indLb = solver->getColLower()[objVarIndex];
00187           if(info->solution_[objVarIndex] < indLb) {
00188             isInf0 = true;
00189           }
00190           if(info->solution_[objVarIndex] > indUb) {
00191             isInf1 = true;
00192           }
00193         }
00194       }
00195 
00196       status0 = simulateBranch (Object, info, branch, solver, result, -1);
00197       if(isInf0) {
00198         status0 = 1; // branch was known to be infeasible
00199         result->setDownStatus(1);
00200       }
00201 
00202       // save current bounds as tightened by the down branch; will be           
00203       // used below to update global bounding box in solver                     
00204       // if status0 == 1 unionLower will be ignored below                         
00205       CoinCopyN (solver->getColLower(), numberColumns, unionLower);
00206       CoinCopyN (solver->getColUpper(), numberColumns, unionUpper);
00207 
00208       // Restore pre-left-branch bounds in solver and problem
00209       for (int j=0; j<numberColumns; j++) {
00210         if(problem_->Lb(j) > unionLower[j]) {
00211           unionLower[j] = problem_->Lb(j);
00212         }
00213         if(problem_->Ub(j) < unionLower[j]) {
00214           unionLower[j] = problem_->Ub(j);
00215         }
00216 
00217         solver->setColLower(j, saveLower [j]);
00218         solver->setColUpper (j, saveUpper [j]);
00219         problem_ -> Lb (j) = saveLower [j];
00220         problem_ -> Ub (j) = saveUpper [j];
00221       }
00222 
00223       /* second direction */
00224 
00225       status1 = simulateBranch (Object, info, branch, solver, result, +1);
00226       if(isInf1) {
00227         status1 = 1; // branch was known to be infeasible
00228         result->setUpStatus(1);
00229       }
00230 
00231 #ifdef TRACE_STRONG
00232       if(problem_->doPrint_) {
00233         printf("Strong on object %d: status0: %d  status1: %d\n", 
00234                result->whichObject(), status0, status1);
00235       }
00236 #endif
00237 
00239 
00240       jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING, "-------\n");
00241 
00242       if (cb) 
00243         cb -> setSimulate (false);
00244 
00246 
00247       bool tightened = false;
00248 
00249       t_chg_bounds *chg_bds = new t_chg_bounds [numberColumns];
00250 
00251       const double *sLb = solver->getColLower();
00252       const double *sUb = solver->getColUpper();
00253 
00254       if(status1 != 1) { // feasible
00255         if(status0 != 1) { // feasible; take union of both branches
00256           for (int j=0; j<numberColumns; j++) {
00257             double maxLb = (problem_->Lb(j) < sLb[j] ? sLb[j] : problem_->Lb(j));
00258             double minUb = (problem_->Ub(j) < sUb[j] ? problem_->Ub(j) : sUb[j]);
00259             problem_->Lb(j) = (unionLower[j] > maxLb ? maxLb : unionLower [j]);
00260             problem_->Ub(j) = (unionUpper[j] < minUb ? minUb : unionUpper [j]);
00261           }
00262         }
00263         else { // keep current bounds; best of problem_ and solver
00264           for (int j=0; j<numberColumns; j++) {
00265             problem_->Lb(j) = (problem_->Lb(j) < sLb[j] ? sLb[j] : problem_->Lb(j));
00266             problem_->Ub(j) = (problem_->Ub(j) < sUb[j] ? problem_->Ub(j) : sUb[j]);
00267           }
00268         }
00269       }
00270       else { // branch 1 infeasible
00271         if(status0 != 1) { // feasible; otherwise both branches are infeasible  
00272                            // keep current inconsistant bounds in solver        
00273           for (int j=0; j<numberColumns; j++) {                                 
00274             problem_->Lb (j) = unionLower [j];                                  
00275             problem_->Ub (j) = unionUpper [j];                                  
00276           }                                                                     
00277         }                                                                       
00278       }                                                                         
00279 
00280       if((status0 == 1) && (status1 == 1)) {
00281         tightened = false;
00282 
00283         // make sure that bounds in solver proves problem is
00284         // infeasible
00285         double lbVar0 = solver->getColLower()[0];
00286         if(lbVar0 < 1) {
00287           solver->setColLower(0, 1);
00288           solver->setColUpper(0, 0);
00289         }
00290         else {
00291           solver->setColUpper(0, lbVar0-1);
00292         }
00293       }
00294       else {
00295         for (int j=0; j<numberColumns; j++) {
00296           if (problem_ -> Lb (j) > initLower [j] + COUENNE_EPS) {
00297             chg_bds [j].setLower (t_chg_bounds::CHANGED);
00298             tightened = true;
00299           }
00300           
00301           if (problem_ -> Ub (j) < initUpper [j] - COUENNE_EPS) {
00302             chg_bds [j].setUpper (t_chg_bounds::CHANGED);
00303             tightened = true;
00304           }
00305         }
00306       }
00307       if (tightened &&                     // have tighter bounds
00308           (problem_ -> doFBBT ()) &&       // selected FBBT
00309           !(problem_ -> btCore (chg_bds))) // tighten again on root
00310 
00311         status0 = status1 = 1;             // if returns false, problem is infeasible
00312 
00313       delete [] chg_bds;
00314 
00315 
00316       if((status0 != 1) || (status1 != 1)) {
00317 
00318         // set new bounding box as the possibly tightened one (a subset
00319         // of the initial)
00320         for (int j=0; j<numberColumns; j++) {
00321           solver -> setColLower (j, saveLower [j] = problem_ -> Lb (j));
00322           solver -> setColUpper (j, saveUpper [j] = problem_ -> Ub (j));
00323         }
00324       }
00325 
00326       /*
00327         End of evaluation for this candidate object. Possibilities are:
00328 
00329         * Both sides below cutoff; this variable is a candidate for
00330           branching.
00331 
00332         * Both sides infeasible or above the objective cutoff: no
00333           further action here. Break from the evaluation loop and
00334           assume the node will be purged by the caller.
00335 
00336         * One side feasible and below cutoff: Install the branch
00337           (i.e., fix the variable). Possibly break from the evaluation
00338           loop and assume the node will be reoptimised by the caller.
00339       */
00340 
00341       if (status0 == 1 && 
00342           status1 == 1) { // infeasible
00343         returnCode=-1;
00344         break; // exit loop
00345       } else if (status0==1 || status1==1) {
00346         numberStrongFixed_++;
00347         if (!returnCriterion) {
00348           returnCode=1;
00349         } else {
00350           returnCode=2;
00351           break;
00352         }
00353       }
00354 
00355       bool hitMaxTime = ( CoinCpuTime()-timeStart > info->timeRemaining_);
00356       if (hitMaxTime) {
00357         returnCode=3;
00358         break;
00359       }
00360     } // end loop /***********************************/
00361 
00362     if (jnlst_ -> ProduceOutput (J_DETAILED, J_BRANCHING)) {
00363       printf ("strong branching: tightened bounds. ");
00364       // create union of bounding box from both branching directions
00365       for (int j=0; j<numberColumns; j++) {
00366       
00367         if (problem_ -> Lb (j) > Lower0 [j]) printf ("l%d (%g-->%g) ", j,Lower0[j], problem_->Lb (j));
00368         if (problem_ -> Ub (j) < Upper0 [j]) printf ("u%d (%g-->%g) ", j,Upper0[j], problem_->Ub (j));
00369       }
00370 
00371       delete [] Lower0;
00372       delete [] Upper0;
00373     }
00374 
00375     //problem_ -> domain () -> pop (); // discard current point/bounds from problem
00376 
00377     delete [] lpSol;
00378 
00379     jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING, "----------------------done\n\n\n");
00380 
00381 #ifdef TRACE_STRONG2
00382     if(problem_->doPrint_) {
00383       if(pv > -1) {
00384         printf("doSB: beg: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
00385                pv, solver->getColSolution()[pv], solver->getColLower()[pv], solver->getColUpper()[pv]);
00386         printf("doSB: info: x[%d]: %10.4f  lb: %10.4f  ub: %10.4f\n",
00387                pv, info->solution_[pv], info->lower_[pv], info->upper_[pv]);
00388         printf("doSB: problem: lb: %10.4f  ub: %10.4f\n",
00389                problem_->Lb(pv), problem_->Ub(pv));
00390       }
00391     }
00392 #endif
00393 
00394     if (iDo < numberToDo) iDo++; // exited due to infeasibility
00395     assert (iDo <= (int) results_.size());
00396     results_.resize (iDo);
00397 
00398     delete [] unionLower;
00399     delete [] unionUpper;
00400 
00401     delete [] saveLower;
00402     delete [] saveUpper;
00403 
00404     solver -> unmarkHotStart ();     // Delete the snapshot
00405 
00406     //solver -> doingResolve () = true;
00407     branchtime_ += CoinCpuTime () - timeStart;
00408 
00409     jnlst_ -> Printf (J_DETAILED, J_BRANCHING, "Done doStrongBranching\n");
00410 
00411     return returnCode;
00412   }
00413 
00414 // Do one side of strong branching
00415 int CouenneChooseStrong::simulateBranch (OsiObject *Object,
00416                                          OsiBranchingInformation *info,
00417                                          OsiBranchingObject *branch,
00418                                          OsiSolverInterface *solver,
00419                                          Bonmin::HotInfo * result,
00420                                          int direction) {
00421 
00422   bool boundBranch = branch -> boundBranch ();
00423 
00424   int status = -1;
00425 
00426   OsiSolverInterface *thisSolver = 
00427     boundBranch ? solver : solver -> clone ();
00428 
00429   CouenneObject *CouObj = dynamic_cast <CouenneObject *> (Object);
00430 
00431   if ((branch -> branch (thisSolver) > COUENNE_INFINITY) || // branch is infeasible
00432       // Bound tightening if not a CouenneObject -- explicit since
00433       // FBBT is done at ::branch() for CouenneObjects
00434       (!CouObj && !BranchingFBBT (problem_, Object, thisSolver))) {
00435 
00436     status = 1;
00437 
00438     if (direction < 0) result -> setDownStatus (1);
00439     else               result -> setUpStatus   (1);
00440 
00441   } else {
00442 
00443     if (boundBranch) // branching rule is a variable bound, can use hotstart
00444 
00445       thisSolver -> solveFromHotStart ();
00446 
00447     else { // branching rule is more complicated, need a resolve
00448 
00449       int limit;
00450       thisSolver -> getIntParam (OsiMaxNumIterationHotStart, limit);
00451       thisSolver -> setIntParam (OsiMaxNumIteration,         limit); 
00452 
00453       thisSolver -> resolve ();
00454     }
00455 
00456     if (pseudoUpdateLP_ && CouObj && thisSolver -> isProvenOptimal ()) {
00457       CouNumber dist = distance (info -> solution_, thisSolver -> getColSolution (), 
00458                                  problem_ -> nVars ());
00459       if (dist > COUENNE_EPS)
00460         CouObj -> setEstimate (dist, direction < 0 ? 0 : 1);
00461     }
00462   }
00463 
00464   // Can check if we got solution
00465   // status is 0 finished, 1 infeasible and 2 unfinished and 3 is solution
00466 
00467   // only update information if this branch is feasible
00468   if (status < 0)
00469     status = result -> updateInformation (thisSolver, info, this);
00470 
00471   numberStrongIterations_ += thisSolver -> getIterationCount ();
00472 
00473   if ((status == 3) && (trustStrongForSolution_)) {
00474     // new solution already saved
00475     info -> cutoff_ = goodObjectiveValue_;
00476     //problem_ -> setCutOff (goodObjectiveValue_);
00477     status = 0;
00478   }
00479 
00480   if (solver != thisSolver)
00481     delete thisSolver;
00482 
00483   return status;
00484 }
00485 
00488 bool BranchingFBBT (CouenneProblem *problem,
00489                     OsiObject *Object,
00490                     OsiSolverInterface *solver) {
00491 
00492   bool feasible = true;
00493 
00494   if (problem -> doFBBT ()) {
00495 
00496     problem -> domain () -> push (solver);
00497 
00498     int 
00499       indVar = Object  -> columnNumber (),
00500       nvars  = problem -> nVars ();
00501 
00502     // do not perform this if object is not a variable object
00503 
00504     if (indVar >= 0) {
00505 
00506       // Tell Couenne one bound has changed
00507 
00508       t_chg_bounds *chg_bds = new t_chg_bounds [nvars];
00509       chg_bds [indVar].setUpper (t_chg_bounds::CHANGED);
00510       chg_bds [indVar].setLower (t_chg_bounds::CHANGED);
00511 
00512       problem -> installCutOff ();
00513 
00514       if ((feasible = problem -> btCore (chg_bds))) {
00515 
00516         const double
00517           *lb  = solver -> getColLower (),
00518           *ub  = solver -> getColUpper (),
00519           *nLB = problem -> Lb (),
00520           *nUB = problem -> Ub ();
00521 
00522         for (int i=0; i<nvars; i++) {
00523           if (nLB [i] > lb [i]) solver -> setColLower (i, nLB [i]);
00524           if (nUB [i] < ub [i]) solver -> setColUpper (i, nUB [i]);
00525         }
00526       }
00527 
00528       delete [] chg_bds;
00529     }
00530 
00531     problem -> domain () -> pop ();
00532   }
00533 
00534   return feasible;
00535 }

Generated on Wed Nov 30 03:03:58 2011 by  doxygen 1.4.7