00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "assert.h"
00012
00013 #include <stdio.h>
00014
00015 #include "CouenneExprDiv.hpp"
00016 #include "CouenneExprConst.hpp"
00017 #include "CouenneExprClone.hpp"
00018 #include "CouenneExprMul.hpp"
00019 #include "CouenneExprOpp.hpp"
00020 #include "CouenneExprPow.hpp"
00021 #include "CouenneExprInv.hpp"
00022 #include "CouenneExprSub.hpp"
00023 #include "CouenneExprBDiv.hpp"
00024
00025 #include "CouennePrecisions.hpp"
00026
00027 #include "CoinFinite.hpp"
00028
00029 using namespace Couenne;
00030
00031
00032
00033 expression *exprDiv::simplify () {
00034
00035 exprOp:: simplify ();
00036 expression *ret = NULL;
00037
00038 if ((*arglist_) -> Type () == CONST) {
00039
00040 CouNumber c0 = (*arglist_) -> Value ();
00041
00042 if (arglist_ [1] -> Type () == CONST) {
00043
00044 CouNumber c1 = arglist_ [1] -> Value ();
00045
00046 if (fabs (c1) == 0.) {
00047 printf ("Couenne: Warning, division by zero -- "); print (); printf ("\n");
00048 }
00049 else {
00050
00051
00052 ret = new exprConst (c0 / c1);
00053
00054 }
00055 }
00056 else {
00057 if (fabs (c0) == 0.)
00058 return new exprConst (0.);
00059
00060
00061
00062
00063
00064 if (fabs (arglist_ [0] -> Value () - 1.) == 0.) {
00065
00066
00067 ret = new exprInv (arglist_ [1]);
00068 arglist_ [1] = NULL;
00069 }
00070 else {
00071
00072 ret = new exprMul (arglist_ [0], new exprInv (arglist_ [1]));
00073 arglist_ [0] = arglist_ [1] = NULL;
00074 }
00075
00076
00077
00078 }
00079 }
00080 else
00081
00082 if (arglist_ [1] -> Type () == CONST) {
00083
00084
00085 ret = new exprMul (new exprConst (1. / (arglist_ [1] -> Value ())), arglist_ [0]);
00086
00087
00088 arglist_ [0] = NULL;
00089
00090 }
00091
00092 return ret;
00093 }
00094
00095
00096
00097
00098
00099 expression *exprDiv::differentiate (int index) {
00100
00101 bool
00102 diffNum = arglist_ [0] -> dependsOn (index),
00103 diffDen = arglist_ [1] -> dependsOn (index);
00104
00105 if (diffNum) {
00106
00107 if (diffDen) {
00108
00109
00110
00111 return new exprDiv (new exprSub (new exprMul (arglist_ [1] -> differentiate (index),
00112 arglist_ [0] -> clone ()),
00113 new exprMul (arglist_ [1] -> clone (),
00114 arglist_ [0] -> differentiate (index))),
00115 new exprPow (arglist_ [1] -> clone (), new exprConst (2.)));
00116
00117 } else {
00118
00119 return new exprDiv (arglist_ [0] -> differentiate (index),
00120 arglist_ [1] -> clone ());
00121 }
00122
00123 } else {
00124
00125 if (diffDen) {
00126
00127 return new exprOpp (new exprDiv (new exprMul (arglist_ [0] -> clone (),
00128 arglist_ [1] -> differentiate (index)),
00129 new exprPow (arglist_ [1] -> clone (),
00130 new exprConst (2.))));
00131
00132 } else
00133 return new exprConst (0.);
00134 }
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 }
00155
00156
00157
00158
00159 void exprDiv::getBounds (expression *&lb, expression *&ub) {
00160
00161 expression **almin = new expression * [4];
00162 expression **almax = new expression * [4];
00163
00164 arglist_ [0] -> getBounds (almin [0], almin [1]);
00165 arglist_ [1] -> getBounds (almin [2], almin [3]);
00166
00167 almax [0] = new exprClone (almin [0]);
00168 almax [1] = new exprClone (almin [1]);
00169 almax [2] = new exprClone (almin [2]);
00170 almax [3] = new exprClone (almin [3]);
00171
00172 lb = new exprLBDiv (almin, 4);
00173 ub = new exprUBDiv (almax, 4);
00174 }
00175
00176
00177
00178
00179 void exprDiv::getBounds (CouNumber &lb, CouNumber &ub) {
00180
00181
00182
00183 CouNumber ln, un, ld, ud;
00184
00185 arglist_ [0] -> getBounds (ln, un);
00186 arglist_ [1] -> getBounds (ld, ud);
00187
00188 if (ld > 0)
00189 if (ln > 0) lb = safeDiv (ln,ud,-1);
00190 else lb = safeDiv (ln,ld,-1);
00191 else {
00192 if (ud > 0) lb = - COUENNE_INFINITY;
00193 else if (un > 0) lb = safeDiv (un,ud,-1);
00194 else lb = safeDiv (un,ld,-1);
00195 }
00196
00197
00198
00199 if (ld > 0)
00200 if (un < 0) ub = safeDiv (un,ud,1);
00201 else ub = safeDiv (un,ld,1);
00202 else {
00203 if (ud > 0) ub = + COUENNE_INFINITY;
00204 else if (ln < 0) ub = safeDiv (ln,ud,1);
00205 else ub = safeDiv (ln,ld,1);
00206 }
00207 }
00208
00210 bool exprDiv::isInteger () {
00211
00212
00213
00214
00215
00216
00217
00218 CouNumber dl, du, nl, nu;
00219
00220 arglist_ [1] -> getBounds (dl, du);
00221 arglist_ [0] -> getBounds (nl, nu);
00222
00223
00224
00225
00226
00227 bool
00228 denzero = (fabs (dl) < COUENNE_EPS),
00229 numconst = (fabs (nl - nu) < COUENNE_EPS);
00230
00231 if ((fabs (nl) < COUENNE_EPS) &&
00232 numconst &&
00233 !denzero)
00234
00235 return true;
00236
00237
00238
00239 if (fabs (dl - du) < COUENNE_EPS) {
00240
00241 if (fabs (fabs (dl) - 1) < COUENNE_EPS)
00242 return arglist_ [0] -> isInteger ();
00243
00244 if (denzero)
00245 return false;
00246
00247 if (numconst) {
00248
00249 CouNumber quot = nl / dl;
00250
00251 if (fabs (COUENNE_round (quot) - quot) < COUENNE_EPS)
00252 return true;
00253 }
00254 }
00255
00256 return false;
00257 }
00258
00259
00261 void exprDiv::closestFeasible (expression *varind,
00262 expression *vardep,
00263 CouNumber &left,
00264 CouNumber &right) const {
00265
00266 expression *varoth = arglist_ [0];
00267
00268 bool numerator = false;
00269
00270 if (varoth -> Index () == varind -> Index ()) {
00271 varoth = arglist_ [1];
00272 numerator = true;
00273 } else assert (arglist_ [1] -> Index () == varind -> Index ());
00274
00275 CouNumber
00276 x = (*varind) (),
00277 y = (*vardep) (),
00278 c = (*varoth) ();
00279
00280 if (numerator)
00281
00282 if (c < 0.)
00283 if (c*y > x) {assert (c*y > right); right = c*y;}
00284 else {assert (c*y < left); left = c*y;}
00285 else if (c > 0.)
00286 if (c*y < x) {assert (c*y < left); left = c*y;}
00287 else {assert (c*y > right); right = c*y;}
00288 else left = - (right = COIN_DBL_MAX);
00289
00290 else
00291
00292 if (y < 0.)
00293 if (x*y > c) {assert (c/y > right); right = c/y;}
00294 else {assert (c/y < left); left = c/y;}
00295 else if (y > 0.)
00296 if (x*y > c) {assert (c/y < left); left = c/y;}
00297 else {assert (c/y > right); right = c/y;}
00298 else left = - (right = COIN_DBL_MAX);
00299 }
00300
00301
00303 CouNumber exprDiv::gradientNorm (const double *x) {
00304
00305 int
00306 ind0 = arglist_ [0] -> Index (),
00307 ind1 = arglist_ [1] -> Index ();
00308
00309 CouNumber
00310 x0 = (ind0 < 0) ? fabs (arglist_ [0] -> Value ()) : fabs (x [ind0]),
00311 x1 = (ind1 < 0) ? fabs (arglist_ [1] -> Value ()) : fabs (x [ind1]),
00312 x1sq = x1 * x1;
00313
00314 if (x1sq < 1/COUENNE_INFINITY) {
00315 x1sq = 1/COUENNE_INFINITY;
00316 if (x1 < 1/COUENNE_INFINITY)
00317 x1 = 1/COUENNE_INFINITY;
00318 }
00319
00320 if (ind0 < 0)
00321 if (ind1 < 0) return 0.;
00322 else return fabs (x0/(x1sq));
00323 else
00324 if (ind1 < 0) return 1. / x1;
00325 else return sqrt (1. / x1sq + x0*x0 / (x1sq * x1sq));
00326 }