/home/coin/SVN-release/OS-2.4.2/Couenne/src/bound_tightening/operators/impliedBounds-exprSum.cpp

Go to the documentation of this file.
00001 /* $Id: impliedBounds-exprSum.cpp 736 2011-07-04 02:27:17Z pbelotti $
00002  *
00003  * Name:    impliedBounds-exprSum.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: implied bound enforcing for exprSum and exprGroup
00006  *
00007  * (C) Carnegie-Mellon University, 2006-10.
00008  * This file is licensed under the Eclipse Public License (EPL)
00009  */
00010 
00011 #include "CouenneExprSum.hpp"
00012 #include "CouenneExprGroup.hpp"
00013 #include "CouenneConfig.h"
00014 #include "CoinHelperFunctions.hpp"
00015 #include "CoinFinite.hpp"
00016 
00017 using namespace Couenne;
00018 
00020 static CouNumber scanBounds (int, int, int *, CouNumber *, CouNumber *, int *);
00021 
00022 
00025 
00026 bool exprSum::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg, enum auxSign sign) {
00027 
00052   // compute number of pos/neg coefficients and sum up constants
00053 
00054   CouNumber
00055     a0 = 0.,   // constant term in the sum
00056     wl = sign == expression::AUX_GEQ ? -COIN_DBL_MAX : l [wind],
00057     wu = sign == expression::AUX_LEQ ?  COIN_DBL_MAX : u [wind];
00058 
00059   // quick check: if both are infinite, nothing is implied...
00060 
00061   if ((wl < -COUENNE_INFINITY) &&
00062       (wu >  COUENNE_INFINITY))
00063     return false;
00064 
00065   int 
00066     nterms = nargs_, // # nonconstant terms
00067     nlin   = 0;      // # linear terms
00068 
00069   exprGroup *eg = NULL;
00070 
00071   bool isExprGroup = (code () == COU_EXPRGROUP);
00072 
00073   if (isExprGroup) { // if exprGroup, compute no. linear terms
00074 
00075     eg = dynamic_cast <exprGroup *> (this);
00076 
00077     a0   += eg -> getc0 ();
00078     nlin += eg -> lcoeff (). size ();
00079   }
00080 
00081   nterms += nlin;
00082 
00083   // Coefficients and indices of the positive and the negative
00084   // non-constant parts of the sum (at most nlin are negative, as all
00085   // "nonlinear" terms have coefficient 1)
00086 
00087   CouNumber *C1 = (CouNumber *) malloc (nterms * sizeof (CouNumber)),
00088             *C2 = (CouNumber *) malloc (nlin   * sizeof (CouNumber));
00089   int       *I1 = (int       *) malloc (nterms * sizeof (int)),
00090             *I2 = (int       *) malloc (nlin   * sizeof (int));
00091 
00092   int ipos, ineg = ipos = 0; // #pos and #neg terms
00093 
00094   std::set <int> intSet; // contains indices of integer variables
00095 
00096   // classify terms as positive or constant for the exprSum
00097 
00098   for (int i = nargs_; i--;) {
00099 
00100     int index = arglist_ [i] -> Index ();
00101 
00102     if (index < 0) // assume a non variable is a constant
00103       a0 += arglist_ [i] -> Value ();
00104     else {
00105       I1 [ipos]   = index;
00106       C1 [ipos++] = 1.;
00107 
00108       // add entry to integer variable
00109       if (arglist_ [i] -> isDefinedInteger ())
00110         intSet.insert (index);
00111     }
00112   }
00113 
00114   // classify terms as pos/neg or constant for the remaining of exprGroup
00115 
00116   if (isExprGroup) {
00117 
00118     exprGroup::lincoeff &lcoe = eg -> lcoeff ();
00119 
00120     for (register exprGroup::lincoeff::iterator el = lcoe.begin ();
00121          el != lcoe.end (); ++el) {
00122 
00123       CouNumber coe = el -> second;
00124       int       ind = el -> first -> Index ();
00125 
00126       if      (coe >  0.) {I1 [ipos] = ind; C1 [ipos++] = coe;}
00127       else if (coe < -0.) {I2 [ineg] = ind; C2 [ineg++] = coe;}
00128 
00129       if (el -> first -> isDefinedInteger ())
00130         intSet.insert (ind);
00131     }
00132   }
00133 
00134   // realloc to save some memory
00135 
00136   /*printf ("  a0 = %g, w in [%g, %g]\n", a0, wl, wu);
00137   printf ("  pos: "); for (int i=0; i<ipos; i++) printf ("%g x%d [%g,%g] ", C1 [i], I1 [i], l [I1 [i]], u [I1 [i]]); printf ("\n");
00138   printf ("  neg: "); for (int i=0; i<ineg; i++) printf ("%g x%d [%g,%g] ", C2 [i], I2 [i], l [I2 [i]], u [I2 [i]]); printf ("\n");*/
00139 
00140   // now we have correct  I1, I2, C1, C2, ipos, ineg, and a0
00141 
00142   // indices of the variable in I1 or I2 with infinite lower or upper
00143   // bound. If more than one is found, it is set to -2
00144 
00145   int infLo1 = -1, infLo2 = -1,
00146       infUp1 = -1, infUp2 = -1;
00147 
00148   // upper bound of the sum, considering lower/upper bounds of the
00149   // variables but negliging the infinite ones:
00150   //
00151   // lower = $a0 + \sum_{i in I_1: a_i <  \infinity} a_i l_i
00152   //             + \sum_{i in I_2: a_i > -\infinity} a_i u_i$
00153   //
00154   // upper = $a0 + \sum_{i in I_1: a_i <  \infinity} a_i u_i
00155   //             + \sum_{i in I_2: a_i > -\infinity} a_i l_i$
00156 
00157   CouNumber
00158 
00159     lower = a0 + scanBounds (ipos, -1, I1, C1, l, &infLo1)
00160                + scanBounds (ineg, +1, I2, C2, u, &infUp2),
00161 
00162     upper = a0 + scanBounds (ipos, +1, I1, C1, u, &infUp1)
00163                + scanBounds (ineg, -1, I2, C2, l, &infLo2);
00164 
00165   // Now compute lower bound for all or for some of the variables:
00166   // There is a bound for all variables only if both infUp1 and infLo2
00167   // are -1, otherwise there is a bound only for one variable if one
00168   // is -1 and the other is nonnegative.
00169 
00170   bool tighter = false;
00171 
00172   // check if there is room for improvement
00173 
00174   {
00175     CouNumber 
00176       slackL = lower - wl, // each is positive if no implication 
00177       slackU = wu - upper; // 
00178 
00179     // if lower < wl or upper > wu, some bounds can be tightened.
00180     // 
00181     // otherwise, there is no implication, but if lower
00182     //
00183     // steal some work to bound propagation... 
00184 
00185     if ((slackL > 0.) &&
00186         (infLo1 == -1) && 
00187         (infUp2 == -1)) {   // no implication on lower
00188 
00189       // propagate lower bound to w
00190       if ((sign != expression::AUX_LEQ) && (updateBound (-1, l + wind, lower))) {
00191         chg [wind].setLower(t_chg_bounds::CHANGED);
00192         tighter = true;
00193         if (intSet.find (wind)!= intSet.end ()) 
00194           l [wind] = ceil (l [wind] - COUENNE_EPS);
00195       }
00196 
00197       if ((slackU > 0.) &&
00198           (infLo2 == -1) && 
00199           (infUp1 == -1)) { // no implication on upper
00200 
00201         // propagate upper bound to w
00202         if ((sign != expression::AUX_GEQ) && (updateBound (+1, u + wind, upper))) {
00203           chg [wind].setUpper(t_chg_bounds::CHANGED);
00204           tighter = true;
00205           if (intSet.find (wind)!= intSet.end ()) 
00206             u [wind] = floor (u [wind] + COUENNE_EPS);
00207         }
00208 
00209         free (I1); free (I2);
00210         free (C1); free (C2);
00211 
00212         return false; // both bounds were weak, no implication possible
00213       }
00214     }
00215     else if ((slackU > 0.) && 
00216              (infLo2 == -1) && 
00217              (infUp1 == -1)) {
00218 
00219       // propagate upper bound to w
00220       if ((sign != expression::AUX_GEQ) && (updateBound (+1, u + wind, upper))) {
00221         tighter = true;
00222         chg [wind].setUpper(t_chg_bounds::CHANGED);
00223         if (intSet.find (wind)!= intSet.end ()) 
00224           u [wind] = floor (u [wind] + COUENNE_EPS);
00225       }
00226     }
00227   }
00228 
00229   // Subtle... make two copies of lower and upper to avoid updating
00230   // bounds with some previously updated bounds
00231 
00232   // first of all, find maximum index in I1 and I2
00233 
00234   int maxind = -1;
00235 
00236   for (register int i=ipos; i--; I1++) if (*I1 > maxind) maxind = *I1;
00237   for (register int i=ineg; i--; I2++) if (*I2 > maxind) maxind = *I2;
00238 
00239   I1 -= ipos;
00240   I2 -= ineg;
00241 
00242   CouNumber *lc = (CouNumber *) malloc (++maxind * sizeof (CouNumber));
00243   CouNumber *uc = (CouNumber *) malloc (maxind   * sizeof (CouNumber));
00244 
00245   CoinCopyN (l, maxind, lc);
00246   CoinCopyN (u, maxind, uc);
00247 
00248   // Update lowers in I1 and uppers in I2
00249 
00250   if ((infLo1 == -1) && (infUp2 == -1) && (wu < COUENNE_INFINITY / 1e10)) { 
00251     // All finite bounds. All var. bounds can be tightened.
00252 
00253     // tighten upper bound of variables in I1
00254     for (register int i=ipos; i--;) {
00255       int ind = I1 [i];
00256       if ((tighter = (updateBound (+1, u + ind, (wu - lower) / C1 [i] + lc [ind]) || tighter))) {
00257         chg [ind].setUpper(t_chg_bounds::CHANGED);
00258         if (intSet.find (ind)!= intSet.end ()) 
00259           u [ind] = floor (u [ind] + COUENNE_EPS);
00260       }
00261     }
00262 
00263     // tighten lower bound of variables in I2
00264     for (register int i=ineg; i--;) {
00265       int ind = I2 [i];
00266       if ((tighter = (updateBound (-1, l + ind, (wu - lower) / C2 [i] + uc [ind]) || tighter))) {
00267         chg [ind].setLower(t_chg_bounds::CHANGED);
00268         if (intSet.find (ind)!= intSet.end ()) 
00269           l [ind] = ceil (l [ind] - COUENNE_EPS);
00270       }
00271     }
00272   } else
00273 
00274     if ((infLo1 >= 0) && (infUp2 == -1)) {    // There is one infinite lower bound in I1
00275       int ind = I1 [infLo1];
00276       if ((tighter = (updateBound (+1, u + ind, (wu - lower) / C1 [infLo1]) || tighter))) {
00277         chg [ind].setUpper(t_chg_bounds::CHANGED);
00278         if (intSet.find (ind)!= intSet.end ()) 
00279           u [ind] = floor (u [ind] + COUENNE_EPS);
00280       }
00281     }
00282     else 
00283       if ((infLo1 == -1) && (infUp2 >= 0)) {  // There is one infinite upper bound in I2
00284         int ind = I2 [infUp2];
00285         if ((tighter = (updateBound (-1, l + ind, (wu - lower) / C2 [infUp2]) || tighter))) {
00286           chg [ind].setLower(t_chg_bounds::CHANGED);
00287           if (intSet.find (ind)!= intSet.end ()) 
00288             l [ind] = ceil (l [ind] - COUENNE_EPS);
00289         }
00290       }
00291 
00292   // Update uppers in I1 and lowers in I2
00293 
00294   if ((infUp1 == -1) && (infLo2 == -1) && (wl > -COUENNE_INFINITY / 1e10)) { 
00295     // All finite bounds. All var. bounds can be tightened.
00296 
00297     for (register int i=ipos; i--;) {
00298       int ind = I1 [i];
00299       if ((tighter = (updateBound (-1, l + ind, (wl - upper) / C1 [i] + uc [ind]) || tighter))) {
00300         chg [ind].setLower(t_chg_bounds::CHANGED); 
00301         if (intSet.find (ind) != intSet.end ()) 
00302           l [ind] = ceil (l [ind] - COUENNE_EPS);
00303       }
00304     }
00305 
00306     for (register int i=ineg; i--;) {
00307       int ind = I2 [i];
00308       if ((tighter = (updateBound (+1, u + ind, (wl - upper) / C2 [i] + lc [ind]) || tighter))) {
00309         chg [ind].setUpper(t_chg_bounds::CHANGED);
00310         if (intSet.find (ind) != intSet.end ()) 
00311           u [ind] = floor (u [ind] + COUENNE_EPS);
00312       }
00313     }
00314   } else 
00315 
00316     if ((infUp1 >= 0) && (infLo2 == -1)) { // There is one infinite lower bound in I2
00317       int ind = I1 [infUp1];
00318       if ((tighter = (updateBound (-1, l + ind, (wl - upper) / C1 [infUp1]) || tighter))) {
00319         chg [ind].setLower(t_chg_bounds::CHANGED);
00320         if (intSet.find (ind) != intSet.end ()) 
00321           l [ind] = ceil (l [ind] - COUENNE_EPS);
00322       }
00323     }
00324     else 
00325       if ((infUp1 == -1) && (infLo2 >= 0)) {  // There is one infinite upper bound in I1
00326         int ind = I2 [infLo2];
00327         if ((tighter = (updateBound (+1, u + ind, (wl - upper) / C2 [infLo2]) || tighter))) {
00328           chg [ind].setUpper(t_chg_bounds::CHANGED);
00329           if (intSet.find (ind) != intSet.end ()) 
00330             u [ind] = floor (u [ind] + COUENNE_EPS);
00331         }
00332       }
00333 
00334   // ...phew!
00335 
00336   free (I1); free (I2);
00337   free (C1); free (C2);
00338   free (lc); free (uc);
00339 
00340   return tighter;
00341 }
00342 
00343 
00345 
00346 static CouNumber scanBounds (int        num,      
00347                              int        sign,     
00348                              int       *indices,  
00349                              CouNumber *coeff,    
00350                              CouNumber *bounds,   
00351                              int       *infnum) { 
00352 
00353   CouNumber bound = 0.;
00354 
00355   for (register int i = num; i--;) {
00356 
00357     CouNumber bd = bounds [indices [i]];
00358 
00359     // be sensitive here, check for bounds a little within the finite realm
00360 
00361     if (((sign > 0) ? bd : -bd) > COUENNE_INFINITY / 1e10) {
00362 
00363       bounds [indices [i]] = (sign > 0) ? COIN_DBL_MAX : -COIN_DBL_MAX;
00364       //bounds [indices [i]] = (sign > 0) ? 1e300 : -1e300;
00365 
00366       // this variable has an infinite bound, mark it
00367       if      (*infnum == -1) *infnum =  i; // first variable with infinite bound, so far
00368       else if (*infnum >=  0) *infnum = -2; // oops... more than one found, no finite bound
00369     }
00370     else bound += coeff [i] * bd; // no infinity detected, sum a weighted, finite bound
00371   }
00372 
00373   return bound;
00374 }

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