/home/coin/SVN-release/OS-2.0.1/Couenne/src/bound_tightening/obbt_iter.cpp

Go to the documentation of this file.
00001 /* $Id: obbt_iter.cpp 141 2009-06-03 04:19:19Z pbelotti $ */
00002 /*
00003  * Name:    obbt.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: Optimality-Based Bound Tightening
00006  *
00007  * (C) Carnegie-Mellon University, 2006-08.
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "CglCutGenerator.hpp"
00012 #include "CouenneCutGenerator.hpp"
00013 #include "CouenneProblem.hpp"
00014 #include "CouenneSolverInterface.hpp"
00015 
00016 #define OBBT_EPS 1e-3
00017 
00018 // TODO: seems like Clp doesn't like large bounds and crashes on
00019 // explicit bounds around 1e200 or so. For now simply use fictitious
00020 // bounds around 1e14. Fix.
00021 
00023 static bool obbt_updateBound (CouenneSolverInterface *csi, 
00024                               int sense,                   
00025                               CouNumber &bound,            
00026                               bool isint) {                
00027 
00028   //csi -> deleteScaleFactors ();
00029   csi -> setDblParam (OsiDualObjectiveLimit, COIN_DBL_MAX); 
00030   csi -> setDblParam (OsiPrimalObjectiveLimit, (sense==1) ? bound : -bound);
00031   csi -> setObjSense (1); // always minimize, just change the sign of the variable
00032 
00034 
00035   //csi -> resolve_nobt (); // this is a time-expensive part, be considerate...
00036   csi -> resolve (); // this is a time-expensive part, be considerate...
00037 
00039 
00040   if (csi -> isProvenOptimal ()) {
00041 
00042     double opt = csi -> getObjValue ();
00043 
00044     if (sense > 0) 
00045          {if (opt       >bound+OBBT_EPS) {bound=(isint ? ceil (opt-COUENNE_EPS) : opt); return true;}}
00046     else {if ((opt=-opt)<bound-OBBT_EPS) {bound=(isint ? floor(opt+COUENNE_EPS) : opt); return true;}}
00047   }
00048 
00049   return false;
00050 }
00051 
00052 
00054 
00055 int CouenneProblem::
00056 obbt_iter (CouenneSolverInterface *csi, 
00057            t_chg_bounds *chg_bds, 
00058            const CoinWarmStart *warmstart, 
00059            Bonmin::BabInfo *babInfo,
00060            double *objcoe,
00061            int sense, 
00062            int index) const {
00063 
00064   // TODO: do NOT apply OBBT if this is a variable of the form
00065   // w2=c*w1, as it suffices to multiply result. More in general, do
00066   // not apply if w2 is a unary monotone function of w1. Even more in
00067   // general, if w2 is a unary function of w1, apply bound propagation
00068   // from w1 to w2 and mark it as exact (depending on whether it is
00069   // non-decreasing or non-increasing
00070 
00071   //  static int iter = 0;
00072 
00073   std::set <int> deplist;
00074   int deplistsize;
00075 
00076   bool issimple = false;
00077 
00078   exprVar *var = Var (index);
00079 
00080   int
00081     objind  = Obj (0) -> Body () -> Index (),
00082     ncols   = csi -> getNumCols (),
00083     nImprov = 0;
00084 
00085   if ((var -> Type  () == AUX) &&
00086       ((deplistsize = var -> Image () -> DepList (deplist, STOP_AT_AUX)) <= 1)) {
00087 
00088     if (!deplistsize) { // funny, the expression is constant...
00089 
00090       CouNumber value = (*(var -> Image ())) ();
00091 
00092       if (csi -> getColLower () [index] < value - COUENNE_EPS) {
00093         csi -> setColLower (index, value); 
00094         chg_bds    [index].setLowerBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
00095       }
00096       else chg_bds [index].setLowerBits(t_chg_bounds::EXACT);
00097 
00098       if (csi -> getColUpper () [index] > value + COUENNE_EPS) {
00099         csi -> setColUpper (index, value); 
00100         chg_bds    [index].setUpperBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
00101       }
00102       else chg_bds [index].setUpperBits(t_chg_bounds::EXACT);
00103 
00104       issimple = true;
00105 
00106     } else { // the expression only depends on one variable, meaning
00107              // that bound propagation should be sufficient
00108 
00109       int indInd = *(deplist.begin ());
00110 
00111       //      expression *image = var -> Image ();
00112       // TODO: write code for monotone functions...
00113 
00114       if // independent variable is exactly bounded in both ways
00115         ((chg_bds [indInd].lower() & t_chg_bounds::EXACT) && 
00116          (chg_bds [indInd].upper() & t_chg_bounds::EXACT) ||
00117          // or if this expression is of the form w=cx+d, that is, it
00118          // depends on one variable only and it is linear
00119          (var -> Image () -> Linearity () <= LINEAR)) {
00120 
00121         issimple = true;
00122 
00123         CouNumber lb, ub;
00124         var -> Image () -> getBounds (lb, ub);
00125 
00126         if (csi -> getColLower () [index] < lb - COUENNE_EPS) {
00127           csi -> setColLower (index, lb); 
00128           chg_bds      [index].setLowerBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
00129         } else chg_bds [index].setLowerBits(t_chg_bounds::EXACT);
00130 
00131         if (csi -> getColUpper () [index] > ub + COUENNE_EPS) {
00132           csi -> setColUpper (index, ub); 
00133           chg_bds      [index].setUpperBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
00134         } else chg_bds [index].setUpperBits(t_chg_bounds::EXACT);
00135       }
00136     }
00137   }
00138 
00139   // only improve bounds if
00140   if (!issimple &&
00141       ((Var (index) -> Type () == VAR) ||             // it is an original variable 
00142        (Var (index) -> Multiplicity () > 0)) &&       // or its multiplicity is at least 1
00143       (Lb (index) < Ub (index) - COUENNE_EPS) && // in any case, bounds are not equal
00144 
00145       ((index != objind) // this is not the objective
00146        // or it is, so we use it for re-solving // TODO: check!
00147        || ((sense ==  1) && !(chg_bds [index].lower() & t_chg_bounds::EXACT))
00148        )) {
00149        //((sense==-1) && (psense == MAXIMIZE) && !(chg_bds [index].upper() & t_chg_bounds::EXACT)))) {
00150 
00151     bool isInt = (Var (index) -> isInteger ());
00152 
00153     objcoe [index] = sense;
00154 
00155     csi -> setObjective (objcoe);
00156     csi -> setObjSense (1); // minimization
00157 
00158     // TODO: Use something else!
00159 #if 0
00160     for (int iv=0; iv<csi->getNumCols (); iv++) {
00161       if (fabs (csi -> getColLower () [iv]) > 1e7) csi -> setColLower (iv, -1e14);
00162       if (fabs (csi -> getColUpper () [iv]) > 1e7) csi -> setColUpper (iv,  1e14);
00163     }
00164 #endif
00165 
00166     CouNumber &bound = 
00167       (sense == 1) ? 
00168       (Lb (index)) : 
00169       (Ub (index)); 
00170 
00171     // m{in,ax}imize xi on csi
00172 
00173     /*
00174     if (Jnlst()->ProduceOutput(J_MOREVECTOR, J_BOUNDTIGHTENING)) {
00175       Jnlst()->Printf(J_MOREVECTOR, J_BOUNDTIGHTENING,
00176                       "m%simizing x%d [%g,%g] %c= %g\n",
00177             (sense==1) ? "in" : "ax", index, Lb (index), Ub (index),
00178             (sense==1) ? '>'  : '<',  bound); fflush (stdout);
00179       if (Jnlst()->ProduceOutput(J_MOREMATRIX, J_BOUNDTIGHTENING)) {
00180         char fname [20];
00181         sprintf (fname, "m%s_w%03d_%03d", (sense == 1) ? "in" : "ax", index, iter);
00182         //Jnlst()->Printf(J_MOREVECTOR, J_BOUNDTIGHTENING,"writing %s\n", fname);
00183         csi -> writeLp (fname);
00184       }
00185     }
00186     */
00187 
00188     csi -> setWarmStart (warmstart);
00189     //csi -> continuousModel () -> setPerturbation (50);
00190 
00191     /* From ClpSimplex.cpp:
00192        
00193        If you are re-using the same matrix again and again then the
00194        setup time to do scaling may be significant.  Also you may not
00195        want to initialize all values or return all values (especially
00196        if infeasible).  While an auxiliary model exists it will be
00197        faster.  If options -1 then model is switched off.  Otherwise
00198        switched on with following options.
00199 
00200        1 - rhs is constant
00201        2 - bounds are constant
00202        4 - objective is constant
00203        8 - solution in by basis and no djs etc in
00204        16 - no duals out (but reduced costs)
00205        32 - no output if infeasible
00206     */
00207 
00208     //csi -> continuousModel () -> auxiliaryModel (1|8|16|32);
00209 
00210     //Jnlst () -> Printf (J_MATRIX, J_BOUNDTIGHTENING,
00211     //"obbt___ index = %d [sen=%d,bd=%g,int=%d]\n", 
00212     //index, sense, bound, isInt);
00213 
00214     if (obbt_updateBound (csi, sense, bound, isInt)) {
00215 
00216       if (bestSol ()) {
00217         if (sense == 1) {
00218           if ((Lb (index) < bestSol () [index]) && 
00219               (bound       > COUENNE_EPS + bestSol () [index]))
00220             Jnlst()->Printf(J_STRONGWARNING, J_BOUNDTIGHTENING,
00221                             "#### OBBT error on x%d: lb = %g, opt = %g, new lb = %g\n", 
00222                             index, Lb (index), bestSol () [index], bound);
00223         } else {
00224           if ((Ub (index) > bestSol () [index]) && 
00225               (bound       < -COUENNE_EPS + bestSol () [index]))
00226             Jnlst()->Printf(J_STRONGWARNING, J_BOUNDTIGHTENING,
00227                             "#### OBBT error on x%d: ub = %g, opt = %g, new ub = %g\n", 
00228                             index, Ub (index), bestSol () [index], bound);
00229         }
00230       }
00231 
00232       // more conservative, only change (and set CHANGED) if improve
00233 
00234       if (sense==1)
00235         if (csi -> getColLower () [index] < bound - COUENNE_EPS) {
00236           Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,"l_%d: %g --> %g\n", 
00237                           index, csi -> getColLower () [index], bound);
00238           csi -> setColLower (index, bound); 
00239           chg_bds      [index].setLowerBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
00240         } else chg_bds [index].setLowerBits(t_chg_bounds::EXACT);
00241       else
00242         if (csi -> getColUpper () [index] > bound + COUENNE_EPS) {
00243           Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,"u_%d: %g --> %g\n", 
00244                           index, csi -> getColUpper () [index], bound);
00245           csi -> setColUpper (index, bound); 
00246           chg_bds      [index].setUpperBits(t_chg_bounds::CHANGED | t_chg_bounds::EXACT);
00247         } else chg_bds [index].setUpperBits(t_chg_bounds::EXACT);
00248 
00249       /*
00250       if (sense==1) {csi -> setColLower (index, bound); chg_bds [index].lower |= CHANGED | EXACT;}
00251       else          {csi -> setColUpper (index, bound); chg_bds [index].upper |= CHANGED | EXACT;}
00252       */
00253 
00254       // check value and bounds of other variables
00255 
00256       const double *sol = csi -> getColSolution ();
00257 
00258       for (int j=0; j<ncols; j++) 
00259         if ((j!=index) && (j!=objind)) {
00260 
00261           if (sol [j] <= Lb (j) + COUENNE_EPS) chg_bds [j].setLowerBits(t_chg_bounds::EXACT);
00262           if (sol [j] >= Ub (j) - COUENNE_EPS) chg_bds [j].setUpperBits(t_chg_bounds::EXACT);
00263         }
00264 
00265 #if 0
00266       // re-check considering reduced costs (more expensive)
00267 
00268       CouNumber *redcost = NULL;
00269 
00270       // first, compute reduced cost when c = c - e_i, where e_i is
00271       // a vector with all zero except a one in position i. This
00272       // serves as a base to compute modified reduced cost below.
00273 
00274       for (int j=0; j<ncols; j++) 
00275         if ((j!=index) && (j!=objind)) {
00276 
00277           // fake a change in the objective function and compute
00278           // reduced cost. If resulting vector is all positive
00279           // (negative), this solution is also optimal for the
00280           // minimization (maximization) of x_j and the corresponding
00281           // chg_bds[j].lower (.upper) can be set to EXACT.
00282 
00283           if (!(chg_bds [j].lower & EXACT)) {
00284           }
00285 
00286           if (!(chg_bds [j].upper & EXACT)) {
00287           }
00288         }
00289 #endif  
00290 
00291       // re-apply bound tightening -- here WITHOUT reduced cost
00292       // (first argument =NULL is pointer to solverInterface) as csi
00293       // is not our problem
00294 
00295       Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
00296                       "  OBBT: x_%d: [%g, %g]\n", index, 
00297                       csi -> getColLower () [index], 
00298                       csi -> getColUpper () [index]);
00299 
00300       if (doFBBT_ && !(boundTightening (chg_bds, babInfo))) {
00301         Jnlst()->Printf(J_DETAILED, J_BOUNDTIGHTENING,
00302                         "node is infeasible after post-OBBT tightening\n");
00303         return -1; // tell caller this is infeasible
00304       }
00305 
00306       nImprov++;
00307     }
00308 
00309     // if we solved the problem on the objective function's
00310     // auxiliary variable (that is, we re-solved the extended
00311     // problem), it is worth updating the current point (it will be
00312     // used later to generate new cuts).
00313 
00314     // TODO: is it, really? And shouldn't we check the opt sense too?
00315     /*
00316     if ((objind == index) && (csi -> isProvenOptimal ()) && (sense == 1))
00317       update (csi -> getColSolution (), NULL, NULL);
00318     */
00319     // restore null obj fun
00320     objcoe [index] = 0;
00321   }
00322 
00323   if (nImprov && jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_BOUNDTIGHTENING)) {
00324     Jnlst () -> Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING, "OBBT: tightened ", nImprov);
00325     if (jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_BOUNDTIGHTENING)) Var (index) -> print ();
00326     Jnlst () -> Printf (J_ITERSUMMARY, J_BOUNDTIGHTENING, "\n");
00327   }
00328 
00329   return nImprov;
00330 }

Generated on Thu Oct 8 03:02:55 2009 by  doxygen 1.4.7