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

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

Generated on Tue Mar 30 03:04:36 2010 by  doxygen 1.4.7