/home/coin/SVN-release/OS-2.1.0/Couenne/src/expression/operators/bounds/exprBQuad.cpp

Go to the documentation of this file.
00001 /* $Id: exprBQuad.cpp 141 2009-06-03 04:19:19Z pbelotti $ */
00002 /*
00003  * Name:    exprBQuad.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: method to compute value of an expr?BQuad
00006  *
00007  * (C) Carnegie-Mellon University, 2006. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "CoinHelperFunctions.hpp"
00012 #include "exprBQuad.hpp"
00013 
00014 //#define DEBUG
00015 
00016 CouNumber exprQuad::computeQBound (int sign) {
00017 
00018   // 1) loose: disaggregated bound
00019   // 2) tighter: aggregate coefficient per variable
00020   // 3) tightest: solve (convex) QP on available alpha-convexification
00021 
00022   // 1) loose bound
00023   //
00024   // w = a0 + a'x + x'Qx means that its lower bound is
00025   //
00026   // w_l = a0 + a'_+ x_l + a'_- x_u + \sum_
00027 
00028 
00029   // 2) tighter bound -- TODO
00030   //
00031   // Compute lower (if sign == -1) or upper (sign == +1) bound of an
00032   // exprQuad based on the information obtained through
00033   // alpha-convexification, if any, or as follows:
00034   //
00035   // w = a0 + a'x + x'Qx = 
00036   //   = a0 + sum{i} [(a_i + sum {j>=i} q_ij * x_j) * x_i] =
00037   //   = a0 + sum{i} [                          z_i * x_i]
00038   // 
00039   // Thus, some bound on z_i can be computed and a bound on the whole
00040   // expression should be better than what can be obtained by summing
00041   // all bounds separately.
00042   //
00043   // Notice that the above computation is fast and may be better than
00044   // the convexification after some updates in the variable bounds
00045   // without updating the convexification. Notice also that the
00046   // direction can also be vertical, not only horizontal.
00047 
00048   CouNumber bound = c0_, term;
00049 
00050   // derive linear part (obtain constant)
00051   for (lincoeff::iterator el = lcoeff_.begin (); el != lcoeff_.end (); ++el) {
00052 
00053     CouNumber 
00054       coe = el -> second, term = 0.,
00055       lb  = el -> first -> lb (),
00056       ub  = el -> first -> ub ();
00057 
00058     if ((coe < 0.) && (sign < 0) || 
00059         (coe > 0.) && (sign > 0)) 
00060       {if    ((term=ub) >  COUENNE_INFINITY) return (sign < 0)? -COUENNE_INFINITY : COUENNE_INFINITY;}
00061     else {if ((term=lb) < -COUENNE_INFINITY) return (sign < 0)? -COUENNE_INFINITY : COUENNE_INFINITY;}
00062 
00063     bound += coe * term;
00064   }
00065 
00066 #ifdef DEBUG
00067   printf ("quadBound --- linear, %cb = %g\n", (sign < 0) ? 'l' : 'u', bound);
00068 #endif
00069 
00070   // derive quadratic part (obtain linear part)
00071   for (sparseQ::iterator row = matrix_.begin (); row != matrix_.end (); ++row) {
00072 
00073     int xind = row -> first -> Index ();
00074 
00075       CouNumber 
00076         lbi = row -> first -> lb (),
00077         ubi = row -> first -> ub ();
00078 
00079     for (sparseQcol::iterator col = row -> second.begin (); col != row -> second.end (); ++col) {
00080 
00081       int yind = col -> first -> Index ();
00082 
00083       CouNumber coe = col -> second;
00084 
00085       if (xind == yind) { // term of the form q_ii x_i^2
00086 
00087         if ((coe > 0.) && (sign < 0) ||
00088             (coe < 0.) && (sign > 0)) 
00089           term = (ubi < 0) ? (ubi * ubi) : (lbi > 0) ? (lbi * lbi) : 0.; //min{xi^2: xi in [lbi,ubi]
00090         else 
00091           if ((term = CoinMax (lbi*lbi, ubi*ubi)) > COUENNE_INFINITY) 
00092             return (sign < 0) ? -COUENNE_INFINITY : COUENNE_INFINITY;
00093 
00094         term *= coe;
00095 
00096 #ifdef DEBUG
00097         printf ("Qii %d %g %g -> %g\n", xind, coe, term, bound + term);
00098 #endif
00099       } else {
00100 
00101         coe *= 2;
00102 
00103         CouNumber
00104           lbj = col -> first -> lb (),
00105           ubj = col -> first -> ub (),
00106           b1 = coe * lbi * lbj, 
00107           b2 = coe * lbi * ubj,
00108           b3 = coe * ubi * lbj, 
00109           b4 = coe * ubi * ubj;
00110 
00111         if (fabs (lbi) == 0) b1 = b2 = 0;
00112         if (fabs (lbj) == 0) b1 = b3 = 0;
00113         if (fabs (ubi) == 0) b3 = b4 = 0;
00114         if (fabs (ubj) == 0) b2 = b4 = 0;
00115 
00116         if (sign < 0) {
00117           if ((term = CoinMin (CoinMin (b1, b2), CoinMin (b3, b4))) < -COUENNE_INFINITY) 
00118             return -COUENNE_INFINITY; 
00119         } else
00120           if ((term = CoinMax (CoinMax (b1, b2), CoinMax (b3, b4))) >  COUENNE_INFINITY)
00121             return  COUENNE_INFINITY;
00122 
00123 #ifdef DEBUG
00124         printf ("Qij %d %d %g %g -> %g\n", xind, yind, coe, term, bound + term);
00125 #endif
00126       }
00127 
00128       //      if ((i!=j) || (lbi >= 0) || (ubi <= 0))
00129       bound += term;
00130     }
00131   }
00132 
00133   return bound;
00134 }

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