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

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

Generated on Mon May 3 03:05:18 2010 by  doxygen 1.4.7