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

Go to the documentation of this file.
00001 /* $Id: conv-exprPow.cpp 151 2009-06-16 03:31:06Z pbelotti $
00002  *
00003  * Name:    conv-exprPow.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: methods to convexify an expression x^k, k constant
00006  *
00007  * (C) Carnegie-Mellon University, 2006. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include <math.h>
00012 #ifndef M_E
00013 # define M_E  2.7182818284590452354
00014 #endif
00015 
00016 #include "CouenneTypes.hpp"
00017 #include "rootQ.hpp"
00018 #include "exprPow.hpp"
00019 #include "exprExp.hpp"
00020 #include "exprConst.hpp"
00021 #include "exprClone.hpp"
00022 #include "exprMul.hpp"
00023 #include "exprSum.hpp"
00024 #include "exprLog.hpp"
00025 #include "CouennePrecisions.hpp"
00026 #include "CouenneProblem.hpp"
00027 #include "CouenneCutGenerator.hpp"
00028 
00029 
00030 std::map <int, CouNumber> Qroot::Qmap;
00031 
00032 // Create standard formulation of this expression
00033 
00034 exprAux *exprPow::standardize (CouenneProblem *p, bool addAux) {
00035 
00036   expression *ret;
00037 
00038   if (arglist_ [0] -> Type () == CONST) { // expression is a^y, reduce
00039                                           // to exp (x * log a)
00040     exprOp::standardize (p);
00041 
00042     CouNumber base = arglist_ [0] -> Value ();
00043 
00044     if (fabs (base - M_E) < COUENNE_EPS) // is base e = 2.7182818...
00045       ret = new exprExp (new exprClone (arglist_ [1]));
00046     else // no? convert k^x to e^(x log (k))
00047       ret = new exprExp (new exprClone 
00048                          (p -> addAuxiliary (new exprMul (new exprClone (arglist_ [1]), 
00049                                                           new exprConst (log (base))))));
00050   } else if (arglist_ [1] -> Type () != CONST) { //  x^y, convert to exp (y*log(x));
00051 
00052     exprOp::standardize (p);
00053 
00054     ret = new exprExp (new exprClone (p -> addAuxiliary 
00055                                       (new exprMul 
00056                                        (new exprClone (arglist_ [1]), 
00057                                         new exprClone 
00058                                         (p -> addAuxiliary 
00059                                          (new exprLog (new exprClone (arglist_ [0]))))))));
00060 
00061   } else { // expression is x^k, return as it is
00062 
00063     // TODO: check if it's of the form ||x||_k, as this admits a
00064     // better (lower) convexification -- replace exprOp::standardize 
00065     exprOp::standardize (p);
00066 
00067     // if binary -- NO! Bounds not known yet
00068     /*if (arglist_ [0] -> isInteger () &&  // integer
00069         (fabs (p -> Lb (arglist_ [0] -> Index ()))      < COUENNE_EPS) && // >= 0
00070         (fabs (p -> Ub (arglist_ [0] -> Index ()) - 1.) < COUENNE_EPS))   // <= 1
00071       {
00072 
00073         printf ("------------------ returning binary! [%g,%g] ",
00074                 p -> Lb (arglist_ [0] -> Index ()),
00075                 p -> Ub (arglist_ [0] -> Index ()));
00076         print ();
00077         printf ("\n");
00078 
00079       return (addAux ? (p -> addAuxiliary (arglist_ [0])) :  // return same variable
00080               new exprAux (arglist_ [0], p -> domain ()));
00081       }
00082       else*/  // otherwise return normal power
00083     {
00084 
00085       /*printf ("------------------ returning normal [%g,%g] ",
00086               arglist_ [0] -> ),
00087               p -> Ub (arglist_ [0] -> Index ()));
00088       print ();
00089       printf ("\n");*/
00090 
00091       return (addAux ? (p -> addAuxiliary (this)) : new exprAux (this, p -> domain ()));
00092     }
00093   }
00094 
00095   return (addAux ? (p -> addAuxiliary (ret)) : new exprAux (ret, p -> domain ()));
00096 }
00097 
00098 
00099 // generate convexification cut for constraint w = x^k
00100 
00101 void exprPow::generateCuts (expression *aux, const OsiSolverInterface &si, 
00102                             OsiCuts &cs, const CouenneCutGenerator *cg,
00103                             t_chg_bounds *chg, int wind, 
00104                             CouNumber lbw, CouNumber ubw) {
00105 
00106   // after standardization, all such expressions are of the form x^k,
00107   // with k constant
00108 
00109   CouNumber k = arglist_ [1] -> Value ();
00110 
00111   // get bounds of base
00112 
00113   expression *xe = arglist_ [0];
00114 
00115   CouNumber l, u;
00116   xe -> getBounds (l, u);
00117 
00118   int w_ind = aux -> Index ();
00119   int x_ind = xe  -> Index ();
00120 
00121   bool cL = !chg || (chg [x_ind].lower() != t_chg_bounds::UNCHANGED) || cg -> isFirst ();
00122   bool cR = !chg || (chg [x_ind].upper() != t_chg_bounds::UNCHANGED) || cg -> isFirst ();
00123 
00124   CouNumber 
00125     w = (*aux) (), 
00126     x = (*xe)  ();
00127 
00128   // if xl and xu are too close, approximate it as a line: sum the
00129   // segment through the two extreme points (l,l^k) and (u,u^k), and
00130   // the tangent at the midpoint ((l+u)/2, ((l+u)/2)^k)
00131 
00132   if (fabs (u-l) < COUENNE_EPS) {
00133 
00134     CouNumber avg     = 0.5 * (l+u),
00135               avg_k_1 = safe_pow (avg, k-1),
00136               lk      = safe_pow (l,   k),
00137               uk      = safe_pow (u,   k);
00138 
00139     if (cL || cR) 
00140       cg -> createCut (cs, u*lk - l*uk + avg * avg_k_1 * (1-k), 0,
00141                        w_ind, u - l + 1, x_ind, lk-uk - k * avg_k_1);
00142     return;
00143   }
00144 
00145   // classify power
00146 
00147   int intk = 0;
00148 
00149   if (k < - COUENNE_INFINITY) { // w=x^{-inf} means w=0
00150     if (cL || cR) cg -> createCut (cs, 0., 0, w_ind, 1.);
00151     return;
00152   }
00153 
00154   if (k > COUENNE_INFINITY) // w=x^{inf} means not much...
00155     return;
00156 
00157   if (fabs (k) < COUENNE_EPS) { // w = x^0 means w=1
00158     if (cL || cR) cg -> createCut (cs, 1., 0, w_ind, 1.);
00159     return;
00160   }
00161 
00162   bool isInt    =            fabs (k    - (double) (intk = COUENNE_round (k)))    < COUENNE_EPS,
00163        isInvInt = !isInt && (fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS);
00164 
00165   // two macro-cases: 
00166 
00167   if (   (isInt || isInvInt)
00168       && (intk % 2) 
00169       && (k >   COUENNE_EPS) 
00170       && (l < - COUENNE_EPS) 
00171       && (u >   COUENNE_EPS)) {
00172 
00173     // 1) k (or its inverse) is positive, integer, and odd, and 0 is
00174     //    an internal point of the interval [l,u].
00175 
00176     // this case is somewhat simpler than the second, although if the
00177     // bound interval includes the origin we have to resort to
00178     // numerical procedures to find the (unique) root of a polynomial
00179     // Q(x) (see Liberti and Pantelides, 2003).
00180 
00181     CouNumber q = 1;
00182 
00183     if ((l<0) && (u>0)) {
00184 
00185       Qroot qmap;
00186       q = qmap (intk);
00187     }
00188 
00189     int sign;
00190 
00191     if (isInvInt) {
00192       if (cg -> isFirst ()) {
00193         w = (l>0) ? 1 : (u<0) ? -1 : 0;
00194         x = 0;
00195       }
00196       q = safe_pow (q, k);
00197       sign = -1;
00198     }
00199     else {
00200       if (cg -> isFirst ()) {
00201         x = (l>0) ? l : (u<0) ? u : 0;
00202         w = 0;
00203       }
00204       sign = 1;
00205     }
00206 
00207     // don't want big coefficients -- check only when k>1
00208     CouNumber powThres = (k<=1) ? COUENNE_INFINITY : pow (COU_MAX_COEFF, 1./k);
00209 
00210     // lower envelope
00211     if (l > -powThres) {
00212       if (l>0) addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k,   l, u, sign); // 0<l<u, tangents only
00213       else if (u > q * l) { // upper x is after "turning point", add lower envelope
00214         addPowEnvelope        (cg, cs, w_ind, x_ind, x, w, k, q*l, u, sign);
00215         cg      -> addSegment     (cs, w_ind, x_ind, l, safe_pow (l,k), q*l, safe_pow (q*l,k), sign);
00216       } else cg -> addSegment     (cs, w_ind, x_ind, l, safe_pow (l,k), u,   safe_pow (u,  k), sign);
00217     }
00218 
00219     // upper envelope
00220     if (u < powThres) {
00221       if (u<0) addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l,   u, -sign);  // l<u<0, tangents only
00222       else if (l < q * u) { // lower x is before "turning point", add upper envelope
00223         addPowEnvelope        (cg, cs, w_ind, x_ind, x, w, k, l, q*u, -sign);
00224         cg      -> addSegment     (cs, w_ind, x_ind, q*u, safe_pow (q*u,k), u, safe_pow (u,k), -sign);
00225       } else cg -> addSegment     (cs, w_ind, x_ind, l,   safe_pow (l,k),   u, safe_pow (u,k), -sign);
00226     }
00227   }
00228   else {
00229 
00230     // 2) all other cases.
00231 
00232     // if k is real or inv(k) is even, then lift l to max (0,l) but if
00233     // u is also negative, there is no convexification -- this
00234     // function is only defined on non-negative numbers
00235 
00236     if (!isInt 
00237         && (!isInvInt || !(intk % 2))
00238         && (l < - COUENNE_EPS) 
00239         && (u < (l=0)))        // CAUTION! l updated if negative (k real)
00240       return;
00241 
00242     // if k is negative and 0 is an internal point of [l,u], no
00243     // convexification is possible -- just add a segment between two
00244     // tails of the asymptotes.
00245 
00246     if ((k < 0) && 
00247         (l < - COUENNE_EPS) && 
00248         (u >   COUENNE_EPS)) {
00249 
00250       if (!(intk % 2))
00251         cg -> addSegment (cs, w_ind, arglist_ [0] -> Index (), 
00252                           l, safe_pow (l,k), u, safe_pow (u,k), +1);
00253       return;
00254     }
00255 
00256     // Between l and u we have a convex/concave function that needs to
00257     // be enveloped. Standard segment and tangent cuts can be applied.
00258 
00259     // create upper envelope (segment)
00260 
00261     int sign = 1; // sign based on k
00262 
00263     // invert sign if 
00264     if (   ((l < - COUENNE_EPS) && (intk % 2) && (k < -COUENNE_EPS)) // k<0 odd, l<0
00265         || ((u < - COUENNE_EPS) && (intk % 2) && (k >  COUENNE_EPS)) // k>0 odd, u<0
00266         || (fabs (k-0.5) < 0.5 - COUENNE_EPS))                       // k in [0,1]
00267       sign = -1;
00268 
00269     CouNumber powThres = CoinMin (COUENNE_INFINITY, 
00270                                   pow (COU_MAX_COEFF, 1./k)), // don't want big coefficients
00271               powStep  = 1;
00272 
00273     // lower envelope for k negative even
00274     if ((k <  COUENNE_EPS) &&
00275         isInt && !(intk % 2) &&
00276         (l < -COUENNE_EPS) &&       // bounds do not contain 0
00277         (u >  COUENNE_EPS) &&
00278         (l > - powThres) &&         // and are finite
00279         (u <   powThres)) 
00280 
00281       cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l, k), u, safe_pow (u, k), 1);
00282 
00283 
00284     // upper envelope
00285     if ((  (k > COUENNE_EPS)        // when k negative, add only if
00286         || (l > COUENNE_EPS)        // bounds do not contain 0
00287         || (u < - COUENNE_EPS)) &&
00288         (l > - powThres) &&         // and are finite
00289         (u <   powThres) &&
00290         (fabs (l+u) > COUENNE_EPS)) // bounds are not opposite (otherwise it's a variable bound)
00291 
00292       cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l, k), u, safe_pow (u, k), -sign);
00293 
00294     // similarly, pay attention not to add infinite slopes
00295 
00296     if (cg -> isFirst()) {
00297       x = (k<0) ? ((u<0) ? u : (l>0) ? l : 0) : 0;
00298       w = 0;
00299     }
00300 
00301     if (k > COUENNE_EPS) {
00302 
00303       if (u >   powThres) u = CoinMax (x,l) + powStep;
00304       if (l < - powThres) l = CoinMin (x,u) - powStep;
00305     }
00306     else {
00307 
00308       if (fabs (l) < COUENNE_EPS) l =  1. / powThres; // l --> 0+
00309       if (fabs (u) < COUENNE_EPS) u = -1. / powThres; // u --> 0-
00310     }
00311 
00312     addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, u, sign);
00313   }
00314 }

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