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

Go to the documentation of this file.
00001 /* $Id: fake_tightening.cpp 141 2009-06-03 04:19:19Z pbelotti $ */
00002 /*
00003  * Name:    fake_tightening.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: fake single bounds in variables to exclude parts of the solution space 
00006  *
00007  * (C) Carnegie-Mellon University, 2007. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "CouenneProblem.hpp"
00012 #include "BonBabInfos.hpp"
00013 
00014 //#define DEBUG
00015 
00016 #define MAX_ITER  10 // max # fake tightening (inner) iterations 
00017 #define AGGR_MUL  2 // the larger,  the more conservative. Must be > 0
00018 #define AGGR_DIV  2 // the smaller, the more conservative. Must be > 1
00019 
00020 // golden ratio, used to find the ideal bound
00021 const CouNumber phi = 0.5 * (1. + sqrt (5.));
00022 
00023 // create fictitious bounds to tighten current interval
00024 CouNumber fictBounds (char direction,
00025                       CouNumber  x,
00026                       CouNumber  lb,   
00027                       CouNumber  ub) {
00028 
00029   if   (lb < -COUENNE_INFINITY / 10) {
00030     if (ub >  COUENNE_INFINITY / 10) { // ]-inf,+inf[
00031 
00032       if (fabs (x) < COUENNE_EPS) return (direction ? AGGR_MUL : - AGGR_MUL);
00033       else                        return (direction ? AGGR_MUL : - AGGR_MUL) * fabs (x);
00034 
00035     } else { // ]-inf,u]
00036 
00037       if      (x < -COUENNE_EPS) return (direction ? CoinMin (0., (x+ub)/2) : AGGR_MUL * x);
00038       else if (x >  COUENNE_EPS) return (direction ? (x + (ub-x)/AGGR_DIV) : 0);
00039       else                       return (direction ? (ub/AGGR_DIV) : -AGGR_MUL);
00040     }
00041   }
00042   else {
00043     if (ub >  COUENNE_INFINITY / 10) { // [l,+inf[
00044 
00045       if      (x < -COUENNE_EPS) return (direction ? 0 : (x - (x-lb) / AGGR_DIV));
00046       else if (x >  COUENNE_EPS) return (direction ? (AGGR_MUL * x) : CoinMax (0.,(x+lb)/2));
00047       else                       return (direction ? AGGR_MUL : lb/AGGR_DIV);
00048     } else // [l,u]
00049       return (direction ? (x + (ub-x) / AGGR_DIV) : x - (x-lb) / AGGR_DIV);
00050   }
00051 }
00052 
00053 
00054 // Single fake tightening. Return
00055 //
00056 // -1   if infeasible
00057 //  0   if no improvement
00058 // +1   if improved
00059 int CouenneProblem::
00060 fake_tighten (char direction,  // 0: left, 1: right
00061               int index,       // index of the variable tested
00062               const double *X, // point round which tightening is done
00063               CouNumber *olb,  // cur. lower bound
00064               CouNumber *oub,  //      upper
00065               t_chg_bounds *chg_bds,
00066               t_chg_bounds *f_chg) const {
00067   int
00068     ncols    = nVars (),
00069     objind   = Obj (0) -> Body  () -> Index ();
00070 
00071   assert (objind >= 0);
00072 
00073   bool 
00074     tightened = false,
00075     intvar    = variables_ [index] -> isInteger ();
00076 
00077   CouNumber 
00078     xcur      = X [index],
00079     inner     = xcur,                                                 // point closest to current x
00080     outer     = (direction ? oub : olb) [index],                      // point closest to bound
00081     fb        = fictBounds (direction, xcur, Lb (index), Ub (index)); // starting point
00082 
00083   // This is a one-dimensional optimization problem between inner and
00084   // outer, on a monotone function of which we can compute the value
00085   // (with relative expense) but not the derivative.
00086 
00087 #ifdef DEBUG
00088   CouNumber curdb     = Lb (objind);// : Ub (objind);  // current dual bound
00089   printf ("  x_%d.  x = %10g, lb = %g, cutoff = %g-----------------\n", index,xcur,curdb,getCutOff());
00090 #endif
00091 
00092   /*if (index == objind)
00093     printf ("  x_%d [%g,%g].  x = %10g, break at %g, cutoff = %g-----------------\n", 
00094     index, Lb (index), Ub (index), xcur, fb, getCutOff());*/
00095 
00096   for (int iter = 0; iter < MAX_ITER; iter++) {
00097 
00098     if (intvar) {
00099 
00100       if (!direction) {inner = floor (inner); outer = ceil  (outer);}
00101       else            {inner = ceil  (inner); outer = floor (outer);}
00102 
00103       if ( direction && (inner > outer) ||
00104           !direction && (inner < outer)) {
00105 
00106         // apply bound
00107         if (direction) {oub[index] = Ub (index) = fb; chg_bds[index].setUpper(t_chg_bounds::CHANGED);}
00108         else           {olb[index] = Lb (index) = fb; chg_bds[index].setLower(t_chg_bounds::CHANGED);}
00109 
00110         tightened = true;
00111 
00112         // restore initial bound
00113         CoinCopyN (chg_bds, ncols, f_chg);
00114         CoinCopyN (olb, ncols, Lb ());
00115         CoinCopyN (oub, ncols, Ub ());
00116 
00117         break;
00118       }
00119 
00120       if (direction  && ((fb < inner) || (fb > outer)) ||
00121           !direction && ((fb > inner) || (fb < outer)))
00122         fb = 0.5 * (inner + outer);
00123     }
00124 
00125     if (direction) {
00126       Lb (index) = intvar ? ceil (fb) : fb; 
00127       f_chg [index].setLower (t_chg_bounds::CHANGED);
00128     } else {
00129       Ub (index) = intvar ? floor (fb) : fb; 
00130       f_chg [index].setUpper (t_chg_bounds::CHANGED);
00131     }
00132 
00133     //    (direction ? lb_ : ub_) [index] = fb; 
00134 
00135 #ifdef DEBUG
00136     char c1 = direction ? '-' : '>', c2 = direction ? '<' : '-';
00137     printf ("    #%3d: [%+10g -%c %+10g %c- %+10g] /\\/\\ ",iter,olb[index],c1,fb,c2, oub [index]);
00138     printf (" [%10g,%10g]<%g,%g>=> ",Lb (index),Ub (index),CoinMin(inner,outer),CoinMax(inner,outer));
00139 #endif
00140 
00141     bool
00142       feasible  = btCore (f_chg),             // true if feasible with fake bound
00143       betterbds = Lb (objind) > getCutOff (); // true if over cutoff
00144 
00145 #ifdef DEBUG
00146     printf(" [%10g,%10g] lb = %g {fea=%d,btr=%d} ",
00147            Lb (index), Ub (index), Lb (objind),feasible,betterbds);
00148 #endif
00149 
00150     if (feasible && !betterbds) {
00151 
00152       // case 1: too tight, move inner out
00153       inner = fb;
00154 
00155       // restore initial bound
00156       CoinCopyN (chg_bds, ncols, f_chg);
00157       CoinCopyN (olb, ncols, Lb ());
00158       CoinCopyN (oub, ncols, Ub ());
00159 
00160     } else {
00161 
00162       // case 2: tightening valid, apply and move outer in
00163 
00164 #ifdef DEBUG
00165       printf (" --> %cbound [x_%d]: %g --> %g",direction?'U':'L',index,(direction?oub:olb)[index],fb);
00166       if (optimum_ && 
00167           ((!direction &&
00168             (optimum_ [index] >= olb [index]) && 
00169             (optimum_ [index] <= Lb (index) - COUENNE_EPS)) ||
00170            (direction &&
00171             (optimum_ [index] <= oub [index]) && 
00172             (optimum_ [index] >= Ub (index) + COUENNE_EPS)))) {
00173         printf ("fake tightening cuts out optimum: x%d=%g in [%g,%g] but not in [%g,%g]\n",
00174                 index, olb [index], oub [index], Lb (index), Ub (index));
00175       }
00176 #endif
00177 
00178       /*bool do_not_tighten = false;
00179 
00180       // check if cut best known solution
00181       if (optimum_) {
00182         if (direction) {
00183           if ((oub [index] > optimum_ [index]) && 
00184               (fb          < optimum_ [index])) {
00185             printf ("aggressive bt cuts optimum ub %d: %g < %g < %g\n", 
00186                     index, fb, optimum_ [index], oub [index]);
00187             do_not_tighten = true;
00188           }
00189         } else {
00190           if ((olb [index] < optimum_ [index]) && 
00191               (fb          > optimum_ [index])) {
00192             printf ("aggressive bt cuts optimum lb %d: %g < %g < %g\n", 
00193                     index, olb [index], optimum_ [index], fb);
00194             do_not_tighten = true;
00195           }
00196         }
00197         }*/
00198 
00199       //if (!do_not_tighten) {
00200 
00201         // apply bound
00202       if (direction) {
00203         oub [index] = Ub (index) = intvar ? floor (fb) : fb; 
00204         chg_bds [index]. setUpper (t_chg_bounds::CHANGED);
00205       }
00206       else {
00207         olb [index] = Lb (index) = intvar ? ceil (fb) : fb; 
00208         chg_bds [index].setLower (t_chg_bounds::CHANGED);
00209       }
00210 
00211       outer = fb; // we have at least a tightened bound, save it 
00212 
00213       tightened = true;
00214         //}
00215 
00216       // restore initial bound
00217       CoinCopyN (chg_bds, ncols, f_chg);
00218       CoinCopyN (olb, ncols, Lb ());
00219       CoinCopyN (oub, ncols, Ub ());
00220 
00221       //#if BR_TEST_LOG < 0 // for fair testing
00222       // check tightened problem for feasibility
00223       if (!(btCore (chg_bds))) {
00224 #ifdef DEBUG
00225         printf ("\n    pruned by aggressive BT\n");
00226 #endif
00227         return -1;
00228       }
00229       //#endif
00230     }
00231 
00232     // TODO: compute golden section
00233     //fb = (inner + outer) / 2;
00234 
00235     //fb = fictBounds (direction, fb, CoinMin (inner, outer), CoinMax (inner, outer));
00236 
00237     CouNumber diff = fabs (inner-outer);
00238 
00239     if (diff == 0.) break;
00240 
00241     if (diff > 1) {
00242 
00243       CouNumber L = log (diff) / log (10.);
00244 
00245       if (direction) fb = inner + exp (log (10.) * L/2);
00246       else           fb = inner - exp (log (10.) * L/2);
00247 
00248     } else fb = (inner+outer)/2;
00249 
00250     //    if () fb = (          inner + (phi-1) * outer) / phi;
00251     //    else  fb = ((phi-1) * inner +           outer) / phi;
00252 
00253     //  if (!feasible)       
00254     //    fb = fictBounds (direction, xcur, 
00255     //       direction ? lb [index] : outer,
00256     //       direction ? outer      : ub [index]);
00257 
00258 #ifdef DEBUG
00259     printf ("\n");
00260 #endif
00261   }
00262 
00263   Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING, "\n");
00264   if (tightened) 
00265     Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_BOUNDTIGHTENING, 
00266                     "  [x%2d] pruned %s [%g, %g] -- lb = %g cutoff = %g\n", 
00267                     index,direction?"right":"left",
00268                     olb[index],oub[index], Lb (objind), getCutOff ());
00269 
00270   return tightened ? 1 : 0;
00271 }

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