/home/coin/SVN-release/OS-2.4.1/Couenne/src/expression/operators/exprPow.cpp

Go to the documentation of this file.
00001 /* $Id: exprPow.cpp 698 2011-06-20 13:36:43Z pbelotti $
00002  *
00003  * Name:    exprPow.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: definition of powers
00006  *
00007  * (C) Carnegie-Mellon University, 2006-11.
00008  * This file is licensed under the Eclipse Public License (EPL)
00009  */
00010 
00011 #include <math.h>
00012 #include <assert.h>
00013 
00014 #include "CouennePrecisions.hpp"
00015 #include "CouenneExprPow.hpp"
00016 #include "CouenneExprSum.hpp"
00017 #include "CouenneExprMul.hpp"
00018 #include "CouenneExprDiv.hpp"
00019 #include "CouenneExprLog.hpp"
00020 #include "CouenneExprConst.hpp"
00021 #include "CouenneProblem.hpp"
00022 
00023 #include "CouenneConfig.h"
00024 #include "CoinHelperFunctions.hpp"
00025 #include "CoinFinite.hpp"
00026 
00027 using namespace Couenne;
00028 
00030 
00031 expression *exprPow::simplify () {
00032 
00033   exprOp:: simplify ();
00034 
00035   if ((*arglist_) -> Type () == CONST) { // expr = c1 ^ g(x)
00036 
00037     CouNumber c0 = (*arglist_) -> Value ();
00038 
00039     if (arglist_ [1] -> Type () == CONST) { // expr = c1 ^ c2
00040 
00041       CouNumber c1 = arglist_ [1] -> Value ();
00042 
00043       delete arglist_ [0]; 
00044       delete arglist_ [1];
00045 
00046       arglist_ [0] = arglist_ [1] = NULL;
00047 
00048       return new exprConst (pow (c0, c1));
00049     }
00050     else 
00051       if (fabs (c0) <= COUENNE_EPS_SIMPL) 
00052         return new exprConst (0.);
00053   }
00054   else // only need to check if g(x) == 0
00055 
00056     if (arglist_ [1] -> Type () == CONST) {
00057 
00058       CouNumber expon = arglist_ [1] -> Value ();
00059 
00060       if (fabs (expon) <= COUENNE_EPS_SIMPL) // expr = x ^ 0 = 1
00061         return new exprConst (1.);
00062 
00063       else if (fabs (expon - 1.) <= COUENNE_EPS_SIMPL) { // expr = x ^ 1 = x
00064 
00065         delete arglist_ [1];
00066         expression *ret = arglist_ [0];
00067         arglist_ [0] = arglist_ [1] = NULL;
00068         return ret;
00069       }
00070 
00071       else if (fabs (expon + 1.) <= COUENNE_EPS_SIMPL) { // expr = x ^ -1 = 1/x
00072 
00073         delete arglist_ [1];
00074         expression *ret = new exprInv (arglist_ [0]);
00075         arglist_ [0] = arglist_ [1] = NULL;
00076         return ret;
00077       }
00078 
00079       //
00080       // x^k = x for x binary. Too bad we don't know bounds yet, so the code below will give segfault
00081 
00082       //       // is it an integer variable with bounds [-1,0] or [0,1]
00083       //       else if ((arglist_ [0] -> Type () == VAR) && (arglist_ [0] -> isDefinedInteger ())) {
00084 
00085       //        CouNumber lb, ub;
00086       //        arglist_ [0] -> getBounds (lb, ub);
00087 
00088       //        if ((fabs (lb)      < COUENNE_EPS) &&
00089       //            (fabs (ub - 1.) < COUENNE_EPS)) {  // {0,1}
00090 
00091       //          delete arglist_ [1];
00092       //          expression *ret = arglist_ [0];
00093       //          arglist_ [0] = arglist_ [1] = NULL;
00094       //          return ret;
00095 
00096       //        } else if ((fabs (lb + 1.) < COUENNE_EPS) &&
00097       //                   (fabs (ub)      < COUENNE_EPS)) { // {-1,0}
00098 
00099       //          delete arglist_ [1];
00100       //          expression *ret = new exprOpp (arglist_ [0]);
00101       //          arglist_ [0] = arglist_ [1] = NULL;
00102       //          return ret;
00103       //        }
00104       //       }
00105 
00106     }
00107 
00108   return NULL;
00109 }
00110 
00111 
00113 
00114 expression *exprPow::differentiate (int index) {
00115 
00116   if (!(arglist_ [0] -> dependsOn (index))  &&
00117       !(arglist_ [1] -> dependsOn (index)))
00118     return new exprConst (0.);
00119 
00120   if (arglist_ [0] -> Type () == CONST) { // k^f(x), k constant
00121 
00122     CouNumber base = arglist_ [0] -> Value ();
00123 
00124     if (base == 0.)
00125       return new exprConst (0.);
00126 
00127     return new exprMul (new exprConst (log (base)),
00128                         new exprMul (new exprPow (new exprConst (base), 
00129                                                   arglist_ [1] -> clone ()),
00130                                      arglist_ [1] -> differentiate (index)));
00131 
00132   } else if (arglist_ [1] -> Type () == CONST) { // f(x)^k, k constant
00133 
00134     CouNumber exponent = arglist_ [1] -> Value ();
00135 
00136     return new exprMul (new exprConst (exponent),
00137                         new exprMul (new exprPow (arglist_ [0] -> clone (),
00138                                                   new exprConst (exponent - 1.)),
00139                                      arglist_ [0] -> differentiate (index)));
00140   }
00141 
00142   // all other cases: f(x)^g(x)
00143 
00144   expression **alm  = new expression * [2];
00145   expression **alp  = new expression * [2];
00146   expression **als  = new expression * [2];
00147   expression **alm1 = new expression * [2];
00148   expression **alm2 = new expression * [2];
00149   expression **ald  = new expression * [2];
00150 
00151   alp [0] = new exprClone (arglist_ [0]);
00152   alp [1] = new exprClone (arglist_ [1]);
00153 
00154   alm [0] = new exprPow (alp, 2);
00155 
00156   alm1 [0] = arglist_ [1] -> differentiate (index);
00157   alm1 [1] = new exprLog (new exprClone (arglist_ [0]));
00158 
00159   als [0] = new exprMul (alm1, 2);
00160 
00161   ald [0] = new exprClone (arglist_ [1]);
00162   ald [1] = new exprClone (arglist_ [0]);
00163 
00164   alm2 [0] = new exprDiv (ald, 2);
00165   alm2 [1] = arglist_ [0] -> differentiate (index);
00166 
00167   als [1] = new exprMul (alm2, 2);
00168 
00169   alm [1] = new exprSum (als, 2);
00170 
00171   return new exprMul (alm, 2);
00172 }
00173 
00174 
00182 
00183 int exprPow::Linearity () {
00184 
00185   if (arglist_ [0] -> Type () == CONST) {
00186 
00187     if (arglist_ [1] -> Type () == CONST) return CONSTANT;
00188     else                                  return NONLINEAR;
00189   }
00190   else {
00191 
00192     double exponent = arglist_ [1] -> Value ();
00193 
00194     if (fabs (exponent - COUENNE_round (exponent)) > COUENNE_EPS)
00195       return NONLINEAR;
00196 
00197     if (arglist_ [1] -> Type () == CONST) { 
00198 
00199       int expInt = (int) COUENNE_round (exponent);
00200 
00201       if (arglist_ [0] -> Linearity () == LINEAR) {
00202 
00203         switch (expInt) {
00204 
00205         case 0:  return CONSTANT;
00206         case 1:  return LINEAR;
00207         case 2:  return QUADRATIC;
00208 
00209         default: return NONLINEAR;
00210         }
00211       }
00212       else 
00213         if (arglist_ [0] -> Linearity () == QUADRATIC) 
00214           switch (expInt) {
00215 
00216           case 0:  return CONSTANT;
00217           case 1:  return QUADRATIC;
00218 
00219           default: return NONLINEAR;
00220           }
00221         else return NONLINEAR;
00222     }
00223     else return NONLINEAR;
00224   }
00225 }
00226 
00227 
00229 bool exprPow::isInteger () {
00230 
00231   // base
00232 
00233   if (!(arglist_ [0] -> isInteger ())) { 
00234 
00235     // base not integer: check if constant and integer
00236     CouNumber lb, ub;
00237     arglist_ [0] -> getBounds (lb, ub);
00238 
00239     if ((fabs (lb - ub) > COUENNE_EPS) ||
00240         !::isInteger (lb))
00241       return false;
00242   }
00243 
00244   // exponent
00245 
00246   if (!(arglist_ [1] -> isInteger ())) { 
00247 
00248     // exponent not defined integer: check if constant and at integer
00249     // value (and positive, or base negative integer)
00250 
00251     CouNumber lb, ub;
00252     arglist_ [1] -> getBounds (lb, ub);
00253 
00254     if ((fabs (lb - ub) > COUENNE_EPS) ||
00255         !::isInteger (lb))
00256       return false;
00257 
00258     if (lb < 0) { // exponent negative, check again base
00259 
00260       arglist_ [0] -> getBounds (lb, ub);
00261 
00262       if ((fabs (lb - ub) > COUENNE_EPS) ||
00263           (fabs (lb) < COUENNE_EPS) ||
00264           !::isInteger (1. / lb))
00265         return false;
00266     }
00267   } else {
00268 
00269     // if base integer and exponent integer, must check that exponent
00270     // is nonnegative
00271 
00272     CouNumber lb, ub;
00273     arglist_ [1] -> getBounds (lb, ub);
00274 
00275     if (lb < .5)
00276       return false;
00277   }
00278 
00279   return true;
00280 }
00281 
00282 
00284 void exprPow::closestFeasible (expression *varind,
00285                                expression *vardep, 
00286                                CouNumber &left,
00287                                CouNumber &right) const {
00288   CouNumber
00289     x  = (*varind) (),
00290     y  = (*vardep) (),
00291     k  = arglist_ [1] -> Value (),
00292     xk = safe_pow (x, k),
00293     yk = safe_pow (y, 1./k);
00294 
00295   int intk = 0;
00296 
00297   bool isInt    =            fabs (k    - (double) (intk = COUENNE_round (k)))    < COUENNE_EPS,
00298        isInvInt = !isInt && (fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS);
00299 
00300   // three cases: 
00301   // 1) k or  1/k odd,        => have either left or right
00302   // 2) k or  1/k even,       => may have both
00303   // 3) k and 1/k fractional  => have either left or right
00304 
00305   if (isInt || isInvInt)
00306 
00307     if (intk % 2) // case 1
00308 
00309       if (k > 0) 
00310         ((y < xk) ? left : right) = yk; // easy, x^k is continuous
00311 
00312       else
00313 
00314         if      (y < 0.)          // third, fourth orthant
00315           if (y < xk) right = yk; // in convex region y < 1/x within third orthant
00316           else        left  = yk; // remaining non-convex area
00317 
00318         else                      // first, second orthant
00319           if (y > xk) left  = yk; // in convex region y > 1/x within first orthant
00320           else        right = yk; // remaining non-convex area
00321 
00322     else // case 2
00323 
00324       if (y <= 0.) // third, fourth orthant => no solution
00325         left = - (right = COIN_DBL_MAX);
00326 
00327       else
00328 
00329         if (k > 0) 
00330 
00331           if (k < 1) // roots, have x >= 0
00332 
00333             if (x > yk) left  = yk;
00334             else        right = yk;
00335 
00336           else
00337 
00338             if (x > yk)       left  =  yk;
00339             else if (x < -yk) right = -yk;
00340             else              left  = - (right = yk);
00341 
00342         else // k negative
00343           if (y < xk) // between asymptotes
00344             left = - (right = yk);
00345           else  // in one of the two convex areas
00346             if (x > 0) left  =  yk;
00347             else       right = -yk;
00348 
00349   else // case 3: assume x bounded from below by 0
00350 
00351     if (k > 0) ((y < xk) ? left : right) = yk;
00352     else       ((y > xk) ? left : right) = yk;
00353 }
00354 
00355 
00357 CouNumber exprPow::gradientNorm (const double *x) {
00358 
00359   int ind0 = arglist_ [0] -> Index ();
00360   CouNumber exponent = arglist_ [1] -> Value ();
00361   return (ind0 < 0) ? 0. : fabs (exponent * safe_pow (x [ind0], exponent - 1));
00362 }
00363 
00364 
00367 bool exprPow::isCuttable (CouenneProblem *problem, int index) const {
00368 
00369   CouNumber exponent = arglist_ [1] -> Value ();
00370 
00371   bool
00372     isInt    = ::isInteger (exponent),
00373     isInvInt = (exponent != 0.) && ::isInteger (1. / exponent);
00374 
00375   int intExp = (isInt ? COUENNE_round (exponent) : (isInvInt ? COUENNE_round (1. / exponent) : 0));
00376 
00377   if (exponent > 0.) {
00378 
00379     if (isInt || isInvInt) {
00380 
00381       if (intExp % 2) return false; // exponent odd or 1/odd
00382 
00383       CouNumber 
00384         x = problem -> X (arglist_ [0] -> Index ()),
00385         y = problem -> X (index);
00386 
00387       if (isInt) return (y <= safe_pow (x, exponent)); // below convex curve ==> cuttable
00388 
00389       return (y >= safe_pow (x, exponent)); // above concave k-th root curve ==> cuttable
00390     } else {
00391 
00392       // non-integer exponent
00393       CouNumber 
00394         x = problem -> X (arglist_ [0] -> Index ()),
00395         y = problem -> X (index);
00396 
00397       return (((exponent <= 1.) && (y >= safe_pow (x, exponent))) ||
00398               ((exponent >= 1.) && (y <= safe_pow (x, exponent))));
00399     }      
00400   } else {
00401 
00402     // non-integer exponent
00403     CouNumber 
00404       x  = problem -> X (arglist_ [0] -> Index ()),
00405       y  = problem -> X (index),
00406       lb = problem -> Lb (index),
00407       ub = problem -> Ub (index);
00408 
00409     if (isInt || isInvInt)
00410 
00411       if (!(intExp % 2)) return (((lb > 0) || (ub < 0)) && (y * safe_pow (fabs (x), -exponent) <= 1.));
00412       else               return (((lb > 0) || (ub < 0)) && (y * safe_pow (x,        -exponent) <= 1.));
00413     else                 return                            (y * safe_pow (x,        -exponent) <= 1.);
00414   }
00415 }

Generated on Thu Nov 10 03:05:44 2011 by  doxygen 1.4.7