/home/coin/SVN-release/OS-2.4.0/Couenne/src/bound_tightening/operators/impliedBounds-exprQuad.cpp

Go to the documentation of this file.
00001 /* $Id: impliedBounds-exprQuad.cpp 560 2011-04-17 10:01:15Z stefan $
00002  *
00003  * Name:    impliedBounds-exprQuad.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: inferring bounds on independent variables of an exprQuad
00006  *          given bounds on the auxiliary variable
00007  *
00008  * (C) Carnegie-Mellon University, 2006-10.
00009  * This file is licensed under the Eclipse Public License (EPL)
00010  */
00011 
00012 #include <set>
00013 
00014 #include "CouenneExprQuad.hpp"
00015 #include "CouenneConfig.h"
00016 #include "CoinFinite.hpp"
00017 #include "CoinHelperFunctions.hpp"
00018 
00019 using namespace Couenne;
00020 
00021 // struct compExpr;
00022 
00028 struct compVar {
00029   inline bool operator () (exprVar* e0, exprVar* e1) const
00030   {return (e0 -> Index () < e1 -> Index ());}
00031 };
00032 
00033 
00036 
00037 bool exprQuad::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg, enum auxSign sign) {
00038 
00039   //return false; // !!!
00040 
00041   // three implied bound computations, of increasing complexity. The
00042   // first reduces this sum to a sum of aux and tries, after
00043   // tightening bounds of each aux, to infer new bounds
00044 
00045   CouNumber 
00046     wl = sign == expression::AUX_GEQ ? -COIN_DBL_MAX : l [wind],
00047     wu = sign == expression::AUX_LEQ ?  COIN_DBL_MAX : u [wind];
00048 
00049 #ifdef DEBUG
00050   printf ("################ implied bounds: [%g,%g], ", wl, wu);
00051   print (); printf ("\n");
00052 #endif
00053 
00054   // Nevermind about the nonlinear part stored in arglist_...
00055 
00056   //std::set <int> indexSet;
00057   std::set <exprVar *, compVar> indexSet;
00058 
00059   // insert indices of linear part
00060   for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el)
00061     indexSet.insert (el -> first);
00062 
00063   // insert indices of quadratic part
00064   for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00065 
00066     indexSet.insert (row -> first);
00067 
00068     for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col)
00069       indexSet.insert (col -> first);
00070   }
00071 
00072   // CAUTION: this relies on the first version of bound expressions
00073   // for quadratic form, i.e. the sum of bounds of independent terms
00074 
00075   CouNumber qMin, qMax;
00076 
00077   int indInfLo = -1, indInfUp = -1;
00078 
00079   // get bounds for nonlinear part of the sum
00080   expression *qlb, *qub;
00081   exprSum::getBounds (qlb, qub);
00082   qMin = (*qlb) (); delete qlb;
00083   qMax = (*qub) (); delete qub;
00084   //exprSum::getBounds (qMin, qMax);
00085 
00086   if (qMin < -COUENNE_INFINITY) indInfLo = -2;
00087   if (qMax >  COUENNE_INFINITY) indInfUp = -2;
00088 
00089   if ((indInfLo == -2) && 
00090       (indInfUp == -2)) 
00091     return false;
00092   
00093 #ifdef DEBUG
00094   printf ("1st phase... inf=(%d,%d) q=[%g,%g].\n", indInfLo, indInfUp, qMin, qMax);
00095   for (std::set <int>:: iterator iter = indexSet.begin ();
00096        iter != indexSet.end (); ++iter) 
00097     printf ("%4d [%+6g %+6g]\n", *iter, l [iter -> Index ()], u [iter -> Index ()]);
00098 #endif
00099 
00100   // compute bound on expression using only finite variable bounds
00101 
00102   computeQuadFiniteBound (qMin, qMax, l, u, indInfLo, indInfUp);
00103 
00104   // similar to impliedBounds-exprGroup. indInf* are -1 if no
00105   // variables are unbounded, i>0 if variable x_i is unbounded and -2
00106   // if at least two are, thus making implied bounds useless.
00107 
00108   qMin += c0_;
00109   qMax += c0_;
00110 
00111   if (((indInfLo == -2) && (indInfUp == -2)) || // at least two variables are unbounded
00112       ((indInfLo == -1) && (indInfUp == -1) &&  // or, none is (and qMin, qMax are exact) 
00113        (qMin > wl) && (qMax < wu))) // but explicit bounds are loose
00114     return false;
00115 
00116 #ifdef DEBUG
00117   printf ("2nd phase... inf=(%d,%d) q=[%g,%g].\n", indInfLo, indInfUp, qMin, qMax);
00118 #endif
00119 
00121   //
00122   // now fill in b_i (constant term in both linear coefficient vectors)
00123 
00124   // prepare data structure for scanning all variables
00125   int
00126     minindex = (*(indexSet.begin  ())) -> Index (), // minimum index
00127     maxindex = (*(indexSet.rbegin ())) -> Index (), // maximum index
00128     nvars = maxindex - minindex + 1;           // upper bound on # variables involved
00129 
00130   CouNumber 
00131     *linCoeMin = new CouNumber [nvars], // min coeff of var x_i
00132     *linCoeMax = new CouNumber [nvars], // max coeff of var x_i
00133     *qii       = new CouNumber [nvars], // quadratic coeff
00134     *bCutLb    = new CouNumber [nvars], // amount to be drawn from the overall lower bound
00135     *bCutUb    = new CouNumber [nvars]; // amount to be drawn from the overall upper bound
00136 
00137   // fill all with 0
00138   CoinFillN (linCoeMin, nvars, 0.);
00139   CoinFillN (linCoeMax, nvars, 0.);
00140   CoinFillN (qii,       nvars, 0.);
00141   CoinFillN (bCutLb,    nvars, 0.);
00142   CoinFillN (bCutUb,    nvars, 0.);
00143 
00144   // assume all coefficients are finite
00145   for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el) {
00146     //  for (int i=0; i<nlterms_; i++) {
00147 
00148     int ind = el -> first -> Index ();//index_ [i];
00149 
00150     CouNumber
00151       coe = el -> second, //coeff_ [i], 
00152       li  = l [ind], 
00153       ui  = u [ind];
00154 
00155     ind -= minindex;
00156 
00157     linCoeMin [ind] += coe;
00158     linCoeMax [ind] += coe;
00159 
00160     if (coe > 0) { // contribution of linear term to bound depends on its coefficient
00161       if (li > -COUENNE_INFINITY) bCutLb [ind] += coe * li;
00162       if (ui <  COUENNE_INFINITY) bCutUb [ind] += coe * ui;
00163     } else {
00164       if (ui <  COUENNE_INFINITY) bCutLb [ind] += coe * ui;
00165       if (li > -COUENNE_INFINITY) bCutUb [ind] += coe * li;
00166     }
00167   }
00168 
00169 #ifdef DEBUG
00170   printf ("linear filling (%d,%d): -----------------------\n", minindex, maxindex);
00171   for (std::set <int>:: iterator iter = indexSet.begin ();
00172        iter != indexSet.end (); ++iter) 
00173     printf ("%4d [%+6g %+6g] [%+6g %+6g]\n", iter -> Index (), 
00174             linCoeMin [iter -> Index () - minindex], linCoeMax [iter -> Index () - minindex],
00175             bCutLb    [iter -> Index () - minindex], bCutUb    [iter -> Index () - minindex]);
00176 #endif
00177 
00178   // fill in remaining linear coefficients and quadratic ones
00179   for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00180 
00181     for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
00182 
00183       int
00184         qi = row -> first -> Index (),
00185         qj = col -> first -> Index (); //qindexJ_ [i];
00186 
00187       CouNumber coe = col -> second, //qcoeff_ [i],
00188         li = l [qi], lj = l [qj], 
00189         ui = u [qi], uj = u [qj];
00190 
00191       if (qi == qj) { // quadratic term
00192 
00193         qi -= minindex;
00194 
00195         qii [qi] = coe; // quadratic term
00196 
00197         CouNumber
00198           maxbUb = CoinMax (fabs (li), fabs (ui)),
00199           maxbLb = (li >= 0) ? (li) : (ui <= 0) ? (ui) : 0;
00200 
00201         if (maxbUb > COUENNE_INFINITY) maxbUb = 0;
00202 
00203         maxbUb *= maxbUb * coe;
00204         maxbLb *= maxbLb * coe;
00205 
00206         if (coe > 0) {
00207           bCutUb [qi] += maxbUb;
00208           bCutLb [qi] += maxbLb;
00209         } else {
00210           bCutUb [qi] += maxbLb;
00211           bCutLb [qi] += maxbUb;
00212         }
00213       } else { // product term
00214 
00215         coe *= 2;
00216 
00217         CouNumber *b1, *b2;
00218 
00219         if (coe > 0) {b1 = l; b2 = u;} 
00220         else         {b1 = u; b2 = l;}
00221 
00222         b1 += minindex;
00223         b2 += minindex;
00224 
00225         qi -= minindex;
00226         qj -= minindex;
00227 
00228         linCoeMin [qi] += coe * b1 [qj];
00229         linCoeMin [qj] += coe * b1 [qi];
00230 
00231         linCoeMax [qi] += coe * b2 [qj];
00232         linCoeMax [qj] += coe * b2 [qi];
00233 
00234         CouNumber
00235           addLo = CoinMin (CoinMin (li*lj, ui*uj),
00236                            CoinMin (ui*lj, li*uj)),
00237           addUp = CoinMax (CoinMax (li*lj, ui*uj), 
00238                            CoinMax (ui*lj, li*uj));
00239 
00240         if (addLo < -COUENNE_INFINITY) addLo = 0;
00241         if (addUp >  COUENNE_INFINITY) addUp = 0;
00242 
00243         addLo *= coe;
00244         addUp *= coe;
00245 
00246         if (coe > 0) {
00247           bCutLb [qi] += addLo; bCutUb [qi] += addUp;
00248           bCutLb [qj] += addLo; bCutUb [qj] += addUp;
00249         } else {
00250           bCutLb [qi] += addUp; bCutUb [qi] += addLo;
00251           bCutLb [qj] += addUp; bCutUb [qj] += addLo;
00252         }
00253       }
00254     }
00255   }
00256 
00257 #ifdef DEBUG
00258   printf ("quad filling: -----------------------\n");
00259   for (std::set <int>:: iterator iter = indexSet.begin ();
00260        iter != indexSet.end (); ++iter) 
00261     printf ("%4d [%+6g %+6g] [%+6g %+6g]\n", *iter, 
00262             linCoeMin [iter -> Index () - minindex], linCoeMax [iter -> Index () - minindex],
00263             bCutLb    [iter -> Index () - minindex], bCutUb    [iter -> Index () - minindex]);
00264 #endif
00265 
00266   // Done filling vectors /////////////////////////////////////////////////////////////
00267   // Now improve each independent variable in the set indexSet
00268 
00269   bool one_updated = false;
00270 
00271   for (std::set <exprVar *, compVar>:: iterator iter = indexSet.begin ();
00272        iter != indexSet.end (); ++iter) {
00273 
00274     bool
00275       updatedL = false,
00276       updatedU = false;
00277 
00278     int ind  = (*iter) -> Index (), 
00279         indn = ind - minindex;
00280 
00281     CouNumber 
00282       al = linCoeMin [indn],
00283       au = linCoeMax [indn],
00284       q  = qii       [indn];
00285 
00286 #ifdef DEBUG
00287     CouNumber
00288       ol = l [ind], 
00289       ou = u [ind];
00290 #endif
00291 
00292     if (fabs (q) < COUENNE_EPS) { // case 1: qii is zero, term is "linear"
00293 
00294       if ((al > 0) || (au < 0)) { // otherwise, min and max lin coe
00295                                   // contain zero, not much to do...
00296 
00297         //          l         <= b_i x_i + c <= u
00298         // <===>    l - c_MAX <= b_i x_i     <= u - c_MIN
00299         //
00300         // c_MAX = qMax - bCutUb [indn]
00301         // c_MIN = qMin - bCutLb [indn]
00302 
00303         CouNumber 
00304           l_b = wl - qMax + bCutUb [indn],
00305           u_b = wu - qMin + bCutLb [indn];
00306 
00307         if (al > 0) { // care about min -- see outer if, it means 0 < al < au
00308 
00309           if ((indInfUp == -1) || (indInfUp == ind))
00310             updatedL = updateBound (-1, l + ind, (l_b) / ((l_b < 0) ? al : au)) || updatedL;
00311           if ((indInfLo == -1) || (indInfLo == ind)) 
00312             updatedU = updateBound (+1, u + ind, (u_b) / ((u_b < 0) ? au : al)) || updatedU;
00313 
00314 #ifdef DEBUG
00315           if (l [ind] > ol) printf ("0. l%d: %g --> %g\n", ind, ol, l [ind]);
00316           if (u [ind] < ou) printf ("0. u%d: %g --> %g\n", ind, ou, u [ind]);
00317 #endif
00318         } else { // only look at max, as al < au < 0
00319 
00320           if ((indInfLo == -1) || (indInfLo == ind)) 
00321             updatedL = updateBound (-1, l + ind, (u_b) / ((u_b < 0) ? al : au)) || updatedL;
00322           if ((indInfUp == -1) || (indInfUp == ind)) 
00323             updatedU = updateBound (+1, u + ind, (l_b) / ((l_b < 0) ? au : al)) || updatedU;
00324 
00325 #ifdef DEBUG
00326           if (l [ind] > ol) printf ("1. l%d: %g --> %g\n", ind, ol, l [ind]);
00327           if (u [ind] < ou) printf ("1. u%d: %g --> %g\n", ind, ou, u [ind]);
00328 #endif
00329         }
00330       }
00331 
00332     } else if (q > 0) {
00333 
00334       // skip if constant term is unbounded from below
00335 
00336       if ((indInfLo != -1) && 
00337           (indInfLo != ind)) 
00338         continue;
00339 
00340       // case 2: qii is positive, the parabola is facing upwards and
00341       // we only need to look at w_U = wu
00342 
00343       // there are two parabolas, with linear coefficient linCoeMin
00344       // and linCoeMax, respectively. If both cut the line $w=w_U$
00345       // then take the minimum of the cut points as new lower bound
00346       // and similarly for the new upper bound.
00347 
00348       // Fortunately there are just two values of linCoe which we have
00349       // to look at, as the parabola is of the form $q_{ii} x_i^2 +
00350       // \hat b_i x_i + c \le w_u$ and the solution contains $\sqrt
00351       // {b_i^2 - 4q_{ii}(c-u)}$, whose maximum is attained with b_i
00352       // maximum or minimum, that is, at linCoeMax and linCoeMin.
00353 
00354       CouNumber 
00355         deltaSecond = 4 * q * (qMin - bCutLb [indn] - wu),
00356         deltaUp     = au*au - deltaSecond,
00357         deltaLo     = al*al - deltaSecond;
00358 
00359       // First case, both parabolas cut the line
00360 
00361       if ((deltaUp >= 0) && 
00362           (deltaLo >= 0)) {
00363 
00364         updatedL = updateBound (-1, l + ind, (- au - sqrt (deltaUp)) / (2*q)) || updatedL;
00365         updatedU = updateBound (+1, u + ind, (- al + sqrt (deltaLo)) / (2*q)) || updatedU;
00366 
00367 #ifdef DEBUG
00368         if (l [ind] > ol) printf ("2. l%d: %g --> %g\n", ind, ol, l [ind]);
00369         if (u [ind] < ou) printf ("2. u%d: %g --> %g\n", ind, ou, u [ind]);
00370 #endif
00371 
00372       } else if (deltaUp >= 0) { // only  left-leaning parabola does
00373 
00374         updatedL = updateBound (-1, l + ind, (- au - sqrt (deltaUp)) / (2*q)) || updatedL;
00375         updatedU = updateBound (+1, u + ind, (- au + sqrt (deltaUp)) / (2*q)) || updatedU;
00376         //      addCoeffCut ();
00377 
00378 #ifdef DEBUG
00379         if (l [ind] > ol) printf ("3. l%d: %g --> %g\n", ind, ol, l [ind]);
00380         if (u [ind] < ou) printf ("3. u%d: %g --> %g\n", ind, ou, u [ind]);
00381 #endif
00382 
00383       } else if (deltaLo >= 0) { // only right-leaning parabola does
00384 
00385         updatedL = updateBound (-1, l + ind, (- al - sqrt (deltaLo)) / (2*q)) || updatedL;
00386         updatedU = updateBound (+1, u + ind, (- al + sqrt (deltaLo)) / (2*q)) || updatedU;
00387         //      addCoeffCut ();
00388 
00389 #ifdef DEBUG
00390         if (l [ind] > ol) printf ("4. l%d: %g --> %g\n", ind, ol, l [ind]);
00391         if (u [ind] < ou) printf ("4. u%d: %g --> %g\n", ind, ou, u [ind]);
00392 #endif
00393 
00394       } else { // none of them does, the problem is infeasible
00395 
00396         updatedL = updateBound (-1, l+ind, +1) || updatedL;
00397         updatedU = updateBound (+1, u+ind, -1) || updatedU;
00398 
00399 #ifdef DEBUG
00400         printf ("5. infeasible!\n");
00401 #endif
00402       }
00403 
00404       // if only one parabola cuts the line, take its lower and upper
00405       // bounds as the new bounds
00406 
00407       // TODO: in the 2nd and 3rd cases, add constraint $\hat b \ge
00408       // \sqrt {4q_{ii}(c-u)}$ or the opposite, to ensure that
00409       // parabola is at least tangent to line w=w_U
00410 
00411       // Also, lower bound can be of help if x_i has a bound between
00412       // the two intersections of the parabola with line w=w_L
00413 
00414     } else {
00415 
00416       // case 3: qii is negative, the parabola is facing downwards and
00417       // we only need to look at wl
00418 
00419       // skip if constant term is unbounded from above
00420       if ((indInfUp != -1) && 
00421           (indInfUp != ind)) 
00422         continue;
00423 
00424       CouNumber 
00425         deltaSecond = 4 * q * (qMax - bCutUb [indn] - wl),
00426         deltaUp     = au*au - deltaSecond,
00427         deltaLo     = al*al - deltaSecond;
00428 
00429       // First case, both parabolas cut the line
00430 
00431       if ((deltaUp >= 0) && 
00432           (deltaLo >= 0)) {
00433 
00434         updatedL = updateBound (-1, l + ind, (al - sqrt (deltaLo)) / (-2*q)) || updatedL;
00435         updatedU = updateBound (+1, u + ind, (au + sqrt (deltaUp)) / (-2*q)) || updatedU;
00436 
00437 #ifdef DEBUG
00438         if (l [ind] > ol) printf ("6. l%d: %g --> %g\n", ind, ol, l [ind]);
00439         if (u [ind] < ou) printf ("6. u%d: %g --> %g\n", ind, ou, u [ind]);
00440 #endif
00441       } else if (deltaUp >= 0) { // only  left-leaning parabola does
00442 
00443         updatedL = updateBound (-1, l + ind, (au - sqrt (deltaUp)) / (-2*q)) || updatedU;
00444         updatedU = updateBound (+1, u + ind, (au + sqrt (deltaUp)) / (-2*q)) || updatedL;
00445         //      addCoeffCut ();
00446 
00447 #ifdef DEBUG
00448         if (l [ind] > ol) printf ("7. l%d: %g --> %g\n", ind, ol, l [ind]);
00449         if (u [ind] < ou) printf ("7. u%d: %g --> %g\n", ind, ou, u [ind]);
00450 #endif
00451 
00452       } else if (deltaLo >= 0) { // only right-leaning parabola does
00453 
00454         updatedL = updateBound (-1, l + ind, (al - sqrt (deltaLo)) / (-2*q)) || updatedL;
00455         updatedU = updateBound (+1, u + ind, (al + sqrt (deltaLo)) / (-2*q)) || updatedU;
00456         //      addCoeffCut ();
00457 #ifdef DEBUG
00458         if (l [ind] > ol) printf ("8. l%d: %g --> %g\n", ind, ol, l [ind]);
00459         if (u [ind] < ou) printf ("8. u%d: %g --> %g\n", ind, ou, u [ind]);
00460 #endif
00461 
00462       } else { // none of them does, the problem is infeasible
00463 
00464         updatedL = updateBound (-1, l+ind, +1) || updatedL;
00465         updatedU = updateBound (+1, u+ind, -1) || updatedU;
00466 
00467 #ifdef DEBUG
00468         printf ("9. infeasible\n");
00469 #endif
00470       }
00471     }
00472 
00473     // 
00474 
00475     if (updatedL) {
00476       one_updated = true;
00477       chg [ind].setLower(t_chg_bounds::CHANGED);
00478       if ((*iter) -> isInteger ())
00479         l [ind] = ceil (l [ind] - COUENNE_EPS);
00480     }
00481 
00482     if (updatedU) {
00483       one_updated = true;
00484       chg [ind].setUpper(t_chg_bounds::CHANGED);
00485       if ((*iter) -> isInteger ())
00486         u [ind] = floor (u [ind] + COUENNE_EPS);
00487     }
00488   }
00489 
00490   delete [] linCoeMin;
00491   delete [] linCoeMax;
00492   delete [] qii;
00493   delete [] bCutLb;
00494   delete [] bCutUb;
00495 
00496   return one_updated;
00497 }

Generated on Thu Sep 22 03:05:55 2011 by  doxygen 1.4.7