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

Go to the documentation of this file.
00001 /* $Id: impliedBounds-exprPow.cpp 560 2011-04-17 10:01:15Z stefan $
00002  *
00003  * Name:    impliedBounds-exprPow.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: implied bounds for power operators
00006  *
00007  * (C) Carnegie-Mellon University, 2008-10.
00008  * This file is licensed under the Eclipse Public License (EPL)
00009  */
00010 
00011 #include "CouenneExprPow.hpp"
00012 #include "CouenneExpression.hpp"
00013 #include "CouenneConfig.h"
00014 #include "CoinFinite.hpp"
00015 #include "CoinHelperFunctions.hpp"
00016 
00017 using namespace Couenne;
00018 
00021 
00022 void invPowImplBounds (int, int, CouNumber *, CouNumber *, CouNumber, bool &, bool &, enum expression::auxSign);
00023 
00024 
00027 
00028 bool exprPow::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg, enum expression::auxSign sign) {
00029 
00030   //int xi = arglist_ [0] -> Index ();
00031   //if (xi>=0) printf ("in implBound-pow: %g,%g\n", l [xi], u [xi]);
00032 
00033   bool resL, resU = resL = false;
00034 
00035   if (arglist_ [0] -> Type () == CONST)   // base is constant or zero, nothing to do
00036     return false;
00037 
00038   assert (arglist_ [1] -> Type () == CONST);
00039 
00040   int index = arglist_ [0] -> Index ();
00041 
00042   CouNumber k = arglist_ [1] -> Value (); // exponent
00043 
00044   if ((fabs (k) < COUENNE_EPS) || 
00045       (fabs (k) > COUENNE_INFINITY)) // a null or infinite k is of little use
00046     return false;
00047 
00048   int intk; // integer (or integer inverse of) exponent
00049 
00050   bool 
00051     isint    =           (            fabs (k    - (intk = COUENNE_round (k)))    < COUENNE_EPS),  // k   integer
00052     isinvint = !isint && (k != 0. && (fabs (1./k - (intk = COUENNE_round (1./k))) < COUENNE_EPS)); // 1/k integer
00053 
00054   CouNumber 
00055     wl = ((sign == expression::AUX_GEQ) ? -COIN_DBL_MAX : l [wind]), // lower w
00056     wu = ((sign == expression::AUX_LEQ) ?  COIN_DBL_MAX : u [wind]); // upper w
00057 
00058   if ((isint || isinvint) && (intk % 2)) { 
00059 
00060     // k or 1/k integer and odd, or non-integer --> it is a monotone
00061     // increasing function, apart when k negative
00062 
00063     if (k > 0.) { // simple, just follow bounds
00064 
00065       if (wl > - COUENNE_INFINITY) resL = updateBound (-1, l + index, safe_pow (wl, 1./k)); 
00066       if (wu <   COUENNE_INFINITY) resU = updateBound (+1, u + index, safe_pow (wu, 1./k));
00067 
00068     } else // slightly more complicated, resort to same method as in exprInv
00069       invPowImplBounds (wind, index, l, u, 1./k, resL, resU, sign);
00070   } 
00071   else 
00072     if (isint) { // x^k, k integer and even --> non monotone
00073 
00074       CouNumber bound = (k<0) ? wl : wu;
00075 
00076       // |x| <= b^(1/k), where b is wl or wu depending on k negative
00077       // or positive, respectively
00078 
00079       if (bound > COUENNE_EPS) {
00080 
00081         if (fabs (bound) < COUENNE_INFINITY) {
00082           resL = updateBound (-1, l + index, - safe_pow (bound, 1./k));
00083           resU = updateBound (+1, u + index,   safe_pow (bound, 1./k));
00084         } /*else {
00085           resL = updateBound (-1, l + index, - COUENNE_INFINITY);
00086           resU = updateBound (+1, u + index,   COUENNE_INFINITY);
00087           }*/
00088       }
00089 
00090       // invert check, if bounds on x do not contain 0 we may improve them
00091 
00092       bound = (k>0) ? wl : wu;
00093 
00094       CouNumber xl = l [index], 
00095                 xu = u [index],
00096                 xb = safe_pow (bound, 1./k);
00097 
00098       if      (xl > - xb + COUENNE_EPS) resL = updateBound (-1, l + index,   xb) || resL;
00099       else if (xu <   xb - COUENNE_EPS) resU = updateBound ( 1, u + index, - xb) || resU;
00100 
00101     } else { 
00102 
00103       // Two cases:
00104       // 1) x^k, k=(1/h), h integer and even, or 
00105       // 2) x^k, neither k nor 1/k integer
00106 
00107       CouNumber lb = wl, ub = wu;
00108 
00109       if (k < 0) { // swap bounds as they swap on the curve x^k when 
00110         lb = wu;
00111         ub = wl;
00112       }
00113 
00114       if (lb > 0. || k > 0.) resL = updateBound (-1, l + index, safe_pow (lb, 1./k));
00115 
00116       if ((fabs (ub) < COUENNE_INFINITY) && (ub > 0 || k > 0.)) 
00117         resU = updateBound (+1, u + index, safe_pow (ub, 1./k));
00118       //else                  resU = updateBound (+1, u + index, COUENNE_INFINITY);
00119     }
00120 
00121   if (resL) chg [index].setLower(t_chg_bounds::CHANGED);
00122   if (resU) chg [index].setUpper(t_chg_bounds::CHANGED);
00123 
00124   bool xInt = arglist_ [0] -> isInteger ();
00125 
00126   if ((resL || resU) && xInt) {
00127 
00128     // careful with what "integer" means when a bound is 1e-8 (see minlp/deb[789].nl)
00129     if (resL && (fabs (l [index]) > COUENNE_EPS)) l [index] = ceil  (l [index] - COUENNE_EPS);
00130     if (resU && (fabs (u [index]) > COUENNE_EPS)) u [index] = floor (u [index] + COUENNE_EPS);
00131   }
00132 
00133   return (resL || resU);
00134 }

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