/home/coin/SVN-release/OS-2.0.0/Couenne/src/convex/operators/conv-exprPow-getBounds.cpp

Go to the documentation of this file.
00001 /* $Id: conv-exprPow-getBounds.cpp 139 2009-06-03 03:56:30Z pbelotti $ 
00002  *
00003  * Name:    conv-exprPow-getBounds.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: method to get lower and upper bounds of a power x^y
00006  *
00007  * (C) Carnegie-Mellon University, 2006-09.
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include <math.h>
00012 
00013 #include "CouenneTypes.hpp"
00014 #include "exprPow.hpp"
00015 #include "exprConst.hpp"
00016 #include "exprClone.hpp"
00017 #include "exprMax.hpp"
00018 #include "exprMin.hpp"
00019 #include "exprOpp.hpp"
00020 #include "CouennePrecisions.hpp"
00021 #include "CouenneProblem.hpp"
00022 
00023 
00024 // compute expressions for lower and upper bounds of a power x^y,
00025 // based on the lower/upper bounds of x and y
00026 
00027 void exprPow::getBounds (expression *&lb, expression *&ub) {
00028 
00029   // We have a standardized expression of the form w = x^y, where x or
00030   // y could be constant. Let us study each case separately.
00031 
00032   assert (arglist_ [0] -> Type () != CONST);
00033 
00034   // x is not constant, so it has (possibly different) lower and
00035   // upper bounds. The expression is x^b, b constant (the case x^y
00036   // has been decomposed by simplify() into exp(y log x).
00037 
00038   expression *lbbase, *ubbase;
00039   arglist_ [0] -> getBounds (lbbase, ubbase);
00040 
00041   //    printf ("ubbase = "); ubbase -> print (std::cout); printf ("\n"); 
00042 
00043   if (arglist_ [1] -> Type () == CONST) { 
00044 
00045     // expression = x^b, b!=0. There are four cases:
00046     // 
00047     // 1) b   is integer and odd  (cube, x^5, etc)
00048     // 2) b   is integer and even (square, x^8, etc)
00049     // 3) 1/b is integer and odd  (cube root, x^(1/7), etc)
00050     // 4) 1/b is integer and even (square root, x^(1/4), etc)
00051     // 5) none of the above
00052     //
00053     // For all of these, need to check if the exponent is negative...
00054 
00055     CouNumber expon = arglist_ [1] -> Value ();
00056     int rndexp;
00057 
00058     bool isInt    =  fabs (expon - (rndexp = COUENNE_round (expon))) < COUENNE_EPS,
00059       isInvInt = !isInt &&
00060       ((fabs (expon) > COUENNE_EPS) &&
00061        (fabs (1/expon - (rndexp = COUENNE_round (1/expon))) < COUENNE_EPS));
00062 
00063     if ((isInt || isInvInt) && (rndexp % 2) && (rndexp > 0)) { 
00064 
00065       // the exponent is integer (or inverse integer), odd and
00066       // positive, hence the function is monotone non decreasing
00067 
00068       lb = new exprPow (lbbase, new exprConst (expon));
00069       ub = new exprPow (ubbase, new exprConst (expon));
00070     } 
00071     else {
00072 
00073       // the exponent is either negative, integer even, or fractional
00074 
00075       expression **all = new expression * [6];
00076 
00077       all [0] = new exprOpp   (lbbase);
00078       all [2] = new exprConst (0.);
00079       all [4] = ubbase;
00080 
00081       if (expon > 0) 
00082         all    [1] = new exprPow (new exprClone (lbbase), new exprConst (expon));
00083       else all [1] = new exprPow (new exprClone (ubbase), new exprConst (expon));
00084 
00085       // all [3] is lower bound when lbbase <= 0 <= ubbase
00086 
00087       if (expon > COUENNE_EPS) all [3] = new exprConst (0.);
00088       else if (isInt || isInvInt) {
00089         if (rndexp % 2) 
00090           all [3] = new exprConst (-COUENNE_INFINITY);
00091         else all [3] = new exprMin (new exprClone (all [1]),
00092                                     new exprPow (new exprClone (lbbase), 
00093                                                  new exprConst (expon)));
00094       }
00095       else all [3] = new exprClone (all [1]);
00096 
00097       // all [5] is the lower bound value when lbbase <= ubbase <= 0
00098 
00099       if (expon > COUENNE_EPS) {
00100         if (isInt && !(rndexp % 2))
00101           all [5] = new exprPow (new exprClone (ubbase), new exprConst (expon));
00102         else all [5] = new exprConst (0.);
00103       }
00104       else {
00105         if (isInt || isInvInt) {
00106           if (rndexp % 2)
00107             all    [5] = new exprPow (new exprClone (ubbase), new exprConst (expon));
00108           else all [5] = new exprPow (new exprClone (lbbase), new exprConst (expon));
00109         }
00110         else all [5] = new exprConst (0.);
00111       }
00112 
00113       lb = new exprMin (all, 6);
00114 
00115       // And now the upper bound ///////////////////////////////////
00116 
00117       if (expon > 0) {
00118 
00119         // special case: upper bound depends to variable bounds only:
00120         // $max {lb^k, ub^k}$
00121 
00122         ub = new exprMax (new exprPow (new exprClone (lbbase), new exprConst (expon)),
00123                           new exprPow (new exprClone (ubbase), new exprConst (expon)));
00124 
00125       } else { // from this point on, expon < 0
00126 
00127         expression **alu = new expression * [6];
00128 
00129         alu [0] = new exprClone (all [0]);
00130         alu [2] = new exprConst (0.);
00131         alu [4] = new exprClone (ubbase);
00132 
00133         //if ((isInt || isInvInt) && !(rndexp % 2))
00134         //alu    [1] = new exprPow (new exprClone (ubbase), new exprConst (expon));
00135         //else 
00136 
00137         // if negative exponent and base has nonnegative lower bound,
00138         // the upper bound can only be lb^k
00139         alu [1] = new exprPow (new exprClone (lbbase), new exprConst (expon));
00140 
00141         // alu [3] is upper bound when lbbase <= 0 <= ubbase
00142 
00143         //if (expon < - COUENNE_EPS) 
00144         alu [3] = new exprConst (COUENNE_INFINITY);
00145         //else if (isInt && !(rndexp % 2))
00146         //alu [3] = new exprPow (new exprMax (new exprClone (lbbase), new exprClone (ubbase)),
00147         //new exprConst (expon));
00148         //else alu [3] = new exprPow (new exprClone (ubbase), new exprConst (expon));
00149 
00150         // alu [5] is the upper bound value when lbbase <= ubbase <= 0
00151 
00152         /*if (expon > COUENNE_EPS) {
00153 
00154           if (isInt && !(rndexp % 2)) 
00155             alu [5] = new exprPow (new exprClone(ubbase), new exprConst(expon));
00156           else alu [5] = new exprConst (0.);
00157         }
00158         else {*/
00159         if (isInt || isInvInt) 
00160           alu [5] = new exprPow (new exprClone (ubbase), new exprConst (expon));
00161         else alu [5] = new exprConst (COUENNE_INFINITY);
00162           //}
00163 
00164         ub = new exprMin (alu, 6);
00165       }
00166     }
00167   }
00168   else // should NOT happen, exponent is not constant...
00169     printf ("exprPow::getBounds(): Warning, exponent not constant\n");
00170 
00171   /*CouNumber l, u;
00172   arglist_ [0] -> getBounds (l,u);
00173 
00174   printf ("pow::bound: ["); 
00175   lb -> print (); printf ("=%g, ", (*lb) ());
00176   ub -> print (); printf ("=%g [%g,%g]\n", (*ub) (), l, u);*/
00177 }
00178 
00179 
00180 // get value of lower and upper bound for the expression
00181 void exprPow::getBounds (CouNumber &lb, CouNumber &ub) {
00182 
00183   CouNumber lba, uba, k = (*(arglist_ [1])) ();
00184   arglist_ [0] -> getBounds (lba, uba);
00185   int intk;
00186 
00187   bool 
00188     isInt    =           fabs (k    - (double) (intk = COUENNE_round (k)))    < COUENNE_EPS,
00189     isInvInt = !isInt && fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS;
00190 
00191   if (!isInt && (!isInvInt || !(intk % 2))) {
00192 
00193     // if exponent is fractional or 1/even, the base better be nonnegative
00194 
00195     if (lba < 0.) lba = 0.;
00196     if (uba < 0.) uba = 0.;
00197   }
00198 
00199   if (isInt && !(intk % 2) && (k > 0)) { // x^{2h}
00200 
00201     if (uba < 0) {
00202       lb = safe_pow (-uba, k);
00203       ub = safe_pow (-lba, k);
00204     } else if (lba > 0) {
00205       lb = safe_pow (lba, k);
00206       ub = safe_pow (uba, k);
00207     } else {
00208       lb = 0;
00209       ub = safe_pow (CoinMax (-lba, uba), k);
00210     }
00211 
00212   } else if (k > 0) { // monotone increasing: x^{2h+1} with h integer, or x^h with h real
00213 
00214     lb = safe_pow (lba, k);
00215     ub = safe_pow (uba, k);
00216 
00217   } else if (isInt && !(intk % 2)) { // x^{-2h} or x^{-1/2h} with h integer
00218 
00219     if (uba < 0) {
00220       lb = safe_pow (-lba, k);
00221       ub = safe_pow (-uba, k);
00222     } else if (lba > 0) {
00223       lb = safe_pow (uba, k);
00224       ub = safe_pow (lba, k);
00225     } else {
00226       lb = safe_pow (CoinMax (-lba, uba), k);
00227       ub = COUENNE_INFINITY;
00228     }
00229 
00230   } else { // x^k, k<0
00231 
00232     if (uba < 0) {
00233       lb = safe_pow (uba, k);
00234       ub = safe_pow (lba, k);
00235     } else if (lba > 0) {
00236       lb = safe_pow (uba, k);
00237       ub = safe_pow (lba, k);
00238     } else {
00239       lb = -COIN_DBL_MAX; // !!! may not be reached
00240       ub =  COIN_DBL_MAX;
00241     }
00242   }
00243 }

Generated on Mon Aug 3 03:02:19 2009 by  doxygen 1.4.7