/home/coin/SVN-release/OS-2.2.0/Couenne/src/expression/operators/exprDiv.cpp

Go to the documentation of this file.
00001 /* $Id: exprDiv.cpp 154 2009-06-16 18:52:53Z pbelotti $ */
00002 /*
00003  * Name:    exprDiv.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: definition of divisions
00006  *
00007  * (C) Carnegie-Mellon University, 2006-08. 
00008  * This file is licensed under the Common Public License (CPL)
00009  */
00010 
00011 #include "exprDiv.hpp"
00012 #include "exprConst.hpp"
00013 #include "exprClone.hpp"
00014 #include "exprMul.hpp"
00015 #include "exprInv.hpp"
00016 #include "exprSub.hpp"
00017 #include "exprBDiv.hpp"
00018 
00019 #include "CouennePrecisions.hpp"
00020 
00021 
00022 // simplify division
00023 
00024 expression *exprDiv::simplify () {
00025 
00026   exprOp:: simplify ();
00027 
00028   if ((*arglist_) -> Type () == CONST) { // expr = a / y 
00029 
00030     CouNumber c0 = (*arglist_) -> Value ();
00031 
00032     if (arglist_ [1] -> Type () == CONST) { // expr = a / b
00033 
00034       CouNumber c1 = arglist_ [1] -> Value ();
00035 
00036       delete arglist_ [0]; 
00037       delete arglist_ [1];
00038       arglist_ [0] = arglist_ [1] = NULL;
00039 
00040       return new exprConst (c0 / c1);
00041     }
00042     else {
00043       if (fabs (c0) < COUENNE_EPS_SIMPL) // expr = 0/y
00044         return new exprConst (0.);
00045 
00046       // otherwise, expression = k/y, return k*inv(y)
00047 
00048       expression *ret;
00049 
00050       if (fabs (arglist_ [0] -> Value ()-1) < COUENNE_EPS) {
00051         delete *arglist_;
00052         *arglist_ = NULL;
00053         ret = new exprInv (arglist_ [1]);
00054       }
00055       else ret = new exprMul (arglist_ [0], new exprInv (arglist_ [1]));
00056 
00057       arglist_ = NULL;
00058       return ret;
00059     }
00060   }
00061   else // only need to check if f2 == 0
00062 
00063     if (arglist_ [1] -> Type () == CONST) { // expression = x/h,
00064                                             // transform into (1/h)*x
00065 
00066       expression *ret = new exprMul (new exprConst (1. / (arglist_ [1] -> Value ())), 
00067                                      arglist_ [0]);
00068       delete arglist_ [1];
00069       arglist_ = NULL;
00070       return ret;
00071     }
00072 
00073   return NULL;
00074 }
00075 
00076 
00077 // differentiate quotient of expressions
00078 
00079 expression *exprDiv::differentiate (int index) {
00080 
00081   if (!(arglist_ [0] -> dependsOn (index))  &&
00082       !(arglist_ [1] -> dependsOn (index)))
00083     return new exprConst (0.);
00084 
00085   expression **alm  = new expression * [2];
00086   expression **als  = new expression * [2];
00087   expression **alm2 = new expression * [3];
00088 
00089   exprInv *invg = new exprInv (new exprClone (arglist_ [1]));
00090 
00091   alm [0] = invg; // evaluated before alm2 [2]
00092 
00093   alm2 [0] = new exprClone (arglist_ [0]);
00094   alm2 [1] = arglist_ [1] -> differentiate (index);
00095   alm2 [2] = new exprClone (invg);
00096 
00097   als [0] = arglist_ [0] -> differentiate (index);
00098   als [1] = new exprMul (alm2, 3);
00099 
00100   alm [1] = new exprSub (als, 2);
00101 
00102   return new exprMul (alm, 2);
00103 }
00104 
00105 
00106 // get lower/upper bounds as a function of the arguments' lower/upper
00107 // bounds
00108 void exprDiv::getBounds (expression *&lb, expression *&ub) {
00109 
00110   expression **almin = new expression * [4];
00111   expression **almax = new expression * [4];
00112 
00113   arglist_ [0] -> getBounds (almin [0], almin [1]);
00114   arglist_ [1] -> getBounds (almin [2], almin [3]);
00115 
00116   almax [0] = new exprClone (almin [0]);
00117   almax [1] = new exprClone (almin [1]);
00118   almax [2] = new exprClone (almin [2]);
00119   almax [3] = new exprClone (almin [3]);
00120 
00121   lb = new exprLBDiv (almin, 4);
00122   ub = new exprUBDiv (almax, 4);
00123 }
00124 
00125 
00126 // get lower/upper bounds as a function of the arguments' lower/upper
00127 // bounds
00128 void exprDiv::getBounds (CouNumber &lb, CouNumber &ub) {
00129 
00130   // lower
00131 
00132   CouNumber ln, un, ld, ud;
00133 
00134   arglist_ [0] -> getBounds (ln, un);
00135   arglist_ [1] -> getBounds (ld, ud);
00136 
00137   if (ld > 0)                                      // (?,?,+,+)
00138     if   (ln > 0)    lb = safeDiv (ln,ud,-1);      // (+,+,+,+) --> ln/ud
00139     else             lb = safeDiv (ln,ld,-1);      // (-,?,+,+) --> ln/ld
00140   else { // ld <= 0
00141     if      (ud > 0) lb = - COUENNE_INFINITY;      // (?,?,-,+) --> unbounded
00142     else if (un > 0) lb = safeDiv (un,ud,-1);      // (?,+,-,-) --> un/ud
00143     else             lb = safeDiv (un,ld,-1);      // (-,-,-,-) --> un/ld
00144   }
00145 
00146   // upper
00147 
00148   if (ld > 0)                                     // (ln,un,ld,ud)     lb 
00149     if   (un < 0) ub = safeDiv (un,ud,1);         // (-,-,+,+) --> un/ud
00150     else          ub = safeDiv (un,ld,1);         // (?,+,+,+) --> un/ld
00151   else { // ld <= 0
00152     if      (ud > 0) ub = + COUENNE_INFINITY;     // (?,?,-,+) --> unbounded
00153     else if (ln < 0) ub = safeDiv (ln,ud,1);      // (-,?,-,-) --> ln/ud
00154     else             ub = safeDiv (ln,ld,1);      // (+,+,-,-) --> ln/ld
00155   }
00156 }
00157 
00159 bool exprDiv::isInteger () {
00160 
00161   // only check if arguments (specifically, the denominator) are, *at
00162   // this point in the algorithm*, constant -- due to branching rules,
00163   // for instance. If so, check if the corresponding evaluated
00164   // expression is integer. Otherwise, check if denominator is +1 or
00165   // -1.
00166 
00167   CouNumber dl, du, nl, nu;
00168 
00169   arglist_ [1] -> getBounds (dl, du);
00170   arglist_ [0] -> getBounds (nl, nu);
00171 
00172   //register CouNumber 
00173   //num = (*nl) (), 
00174   //den = (*dl) ();
00175 
00176   bool
00177     denzero  = (fabs (dl)      < COUENNE_EPS),
00178     numconst = (fabs (nl - nu) < COUENNE_EPS);
00179 
00180   if ((fabs (nl) < COUENNE_EPS)  && // numerator is zero
00181       numconst                   && // constant
00182       !denzero)                     // and denominator is nonzero
00183 
00184     return true;
00185 
00186   // otherwise...
00187 
00188   if (fabs (dl - du) < COUENNE_EPS) { // denominator is constant
00189 
00190     if (fabs (fabs (dl) - 1) < COUENNE_EPS) // it is +1 or -1, check numerator
00191       return arglist_ [0] -> isInteger ();
00192 
00193     if (denzero) // it is zero, better leave...
00194       return false;
00195 
00196     if (numconst) { // numerator is constant, too
00197 
00198       CouNumber quot = nl / dl;
00199 
00200       if (fabs (COUENNE_round (quot) - quot) < COUENNE_EPS)
00201         return true;
00202     }
00203   }
00204 
00205   return false;
00206 }
00207 
00208 
00210 void exprDiv::closestFeasible (expression *varind,
00211                                expression *vardep, 
00212                                CouNumber &left,
00213                                CouNumber &right) const {
00214 
00215   expression *varoth = arglist_ [0]; // assume y = c/x
00216 
00217   bool numerator = false;
00218 
00219   if (varoth -> Index () == varind -> Index ()) { // actually y = x/c
00220     varoth = arglist_ [1];
00221     numerator = true;
00222   } else assert (arglist_ [1] -> Index () == varind -> Index ()); // right to assume y = c/x
00223 
00224   CouNumber 
00225     x = (*varind) (),
00226     y = (*vardep) (),
00227     c = (*varoth) ();
00228 
00229   if (numerator) // checking y = x/c
00230 
00231     if (c < 0.)
00232       if (c*y > x) {assert (c*y > right); right = c*y;}
00233       else         {assert (c*y < left);  left  = c*y;}
00234     else if (c > 0.)
00235       if (c*y < x) {assert (c*y < left);  left  = c*y;}
00236       else         {assert (c*y > right); right = c*y;}
00237     else left = - (right = COIN_DBL_MAX);
00238 
00239   else           // checking y = c/x
00240 
00241     if      (y < 0.)
00242       if (x*y > c) {assert (c/y > right); right = c/y;} // convex area in third orthant
00243       else         {assert (c/y < left);  left  = c/y;} // remaining of third+fourth orthant
00244     else if (y > 0.) 
00245       if (x*y > c) {assert (c/y < left);  left  = c/y;} // convex area in first orthant
00246       else         {assert (c/y > right); right = c/y;} // remaining of first+second orthant
00247     else left = - (right = COIN_DBL_MAX);
00248 }
00249 
00250 
00252 CouNumber exprDiv::gradientNorm (const double *x) {
00253 
00254   int 
00255     ind0 = arglist_ [0] -> Index (),
00256     ind1 = arglist_ [1] -> Index ();
00257 
00258   CouNumber
00259     x0 = (ind0 < 0) ? fabs (arglist_ [0] -> Value ()) : fabs (x [ind0]),
00260     x1 = (ind1 < 0) ? fabs (arglist_ [1] -> Value ()) : fabs (x [ind1]),
00261     x1sq = x1 * x1;
00262 
00263   if (x1sq < 1/COUENNE_INFINITY) {
00264     x1sq = 1/COUENNE_INFINITY;
00265     if (x1 < 1/COUENNE_INFINITY) // implied
00266       x1 = 1/COUENNE_INFINITY;
00267   }
00268 
00269   if (ind0 < 0)
00270     if (ind1 < 0) return 0.;                // c/d
00271     else          return fabs (x0/(x1sq)); // c/y
00272   else 
00273     if (ind1 < 0) return 1. / x1;                                // x/d
00274     else          return sqrt (1. / x1sq + x0*x0 / (x1sq * x1sq)); // x/y
00275 }

Generated on Thu Aug 5 03:02:57 2010 by  doxygen 1.4.7