00001
00002
00003
00004
00005
00006
00007
00008
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) {
00036
00037 CouNumber c0 = (*arglist_) -> Value ();
00038
00039 if (arglist_ [1] -> Type () == CONST) {
00040
00041 CouNumber c1 = arglist_ [1] -> Value ();
00042
00043
00044
00045
00046
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
00055
00056 if (arglist_ [1] -> Type () == CONST) {
00057
00058 CouNumber expon = arglist_ [1] -> Value ();
00059
00060 if (fabs (expon) <= COUENNE_EPS_SIMPL)
00061 return new exprConst (1.);
00062
00063 else if (fabs (expon - 1.) <= COUENNE_EPS_SIMPL) {
00064
00065
00066 expression *ret = arglist_ [0];
00067 arglist_ [0] = NULL;
00068
00069 return ret;
00070 }
00071
00072 else if (fabs (expon + 1.) <= COUENNE_EPS_SIMPL) {
00073
00074
00075 expression *ret = new exprInv (arglist_ [0]);
00076 arglist_ [0] = NULL;
00077 return ret;
00078 }
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109 }
00110
00111 return NULL;
00112 }
00113
00114
00116
00117 expression *exprPow::differentiate (int index) {
00118
00119 if (!(arglist_ [0] -> dependsOn (index)) &&
00120 !(arglist_ [1] -> dependsOn (index)))
00121 return new exprConst (0.);
00122
00123 if (arglist_ [0] -> Type () == CONST) {
00124
00125 CouNumber base = arglist_ [0] -> Value ();
00126
00127 if (base == 0.)
00128 return new exprConst (0.);
00129
00130 return new exprMul (new exprConst (log (base)),
00131 new exprMul (new exprPow (new exprConst (base),
00132 arglist_ [1] -> clone ()),
00133 arglist_ [1] -> differentiate (index)));
00134
00135 } else if (arglist_ [1] -> Type () == CONST) {
00136
00137 CouNumber exponent = arglist_ [1] -> Value ();
00138
00139 return new exprMul (new exprConst (exponent),
00140 new exprMul (new exprPow (arglist_ [0] -> clone (),
00141 new exprConst (exponent - 1.)),
00142 arglist_ [0] -> differentiate (index)));
00143 }
00144
00145
00146
00147 expression **alm = new expression * [2];
00148 expression **alp = new expression * [2];
00149 expression **als = new expression * [2];
00150 expression **alm1 = new expression * [2];
00151 expression **alm2 = new expression * [2];
00152 expression **ald = new expression * [2];
00153
00154 alp [0] = new exprClone (arglist_ [0]);
00155 alp [1] = new exprClone (arglist_ [1]);
00156
00157 alm [0] = new exprPow (alp, 2);
00158
00159 alm1 [0] = arglist_ [1] -> differentiate (index);
00160 alm1 [1] = new exprLog (new exprClone (arglist_ [0]));
00161
00162 als [0] = new exprMul (alm1, 2);
00163
00164 ald [0] = new exprClone (arglist_ [1]);
00165 ald [1] = new exprClone (arglist_ [0]);
00166
00167 alm2 [0] = new exprDiv (ald, 2);
00168 alm2 [1] = arglist_ [0] -> differentiate (index);
00169
00170 als [1] = new exprMul (alm2, 2);
00171
00172 alm [1] = new exprSum (als, 2);
00173
00174 return new exprMul (alm, 2);
00175 }
00176
00177
00185
00186 int exprPow::Linearity () {
00187
00188 if (arglist_ [0] -> Type () == CONST) {
00189
00190 if (arglist_ [1] -> Type () == CONST) return CONSTANT;
00191 else return NONLINEAR;
00192 }
00193 else {
00194
00195 double exponent = arglist_ [1] -> Value ();
00196
00197 if (fabs (exponent - COUENNE_round (exponent)) > COUENNE_EPS)
00198 return NONLINEAR;
00199
00200 if (arglist_ [1] -> Type () == CONST) {
00201
00202 int expInt = (int) COUENNE_round (exponent);
00203
00204 if (arglist_ [0] -> Linearity () == LINEAR) {
00205
00206 switch (expInt) {
00207
00208 case 0: return CONSTANT;
00209 case 1: return LINEAR;
00210 case 2: return QUADRATIC;
00211
00212 default: return NONLINEAR;
00213 }
00214 }
00215 else
00216 if (arglist_ [0] -> Linearity () == QUADRATIC)
00217 switch (expInt) {
00218
00219 case 0: return CONSTANT;
00220 case 1: return QUADRATIC;
00221
00222 default: return NONLINEAR;
00223 }
00224 else return NONLINEAR;
00225 }
00226 else return NONLINEAR;
00227 }
00228 }
00229
00230
00232 bool exprPow::isInteger () {
00233
00234
00235
00236 if (!(arglist_ [0] -> isInteger ())) {
00237
00238
00239 CouNumber lb, ub;
00240 arglist_ [0] -> getBounds (lb, ub);
00241
00242 if ((fabs (lb - ub) > COUENNE_EPS) ||
00243 !::isInteger (lb))
00244 return false;
00245 }
00246
00247
00248
00249 if (!(arglist_ [1] -> isInteger ())) {
00250
00251
00252
00253
00254 CouNumber lb, ub;
00255 arglist_ [1] -> getBounds (lb, ub);
00256
00257 if ((fabs (lb - ub) > COUENNE_EPS) ||
00258 !::isInteger (lb))
00259 return false;
00260
00261 if (lb < 0) {
00262
00263 arglist_ [0] -> getBounds (lb, ub);
00264
00265 if ((fabs (lb - ub) > COUENNE_EPS) ||
00266 (fabs (lb) < COUENNE_EPS) ||
00267 !::isInteger (1. / lb))
00268 return false;
00269 }
00270 } else {
00271
00272
00273
00274
00275 CouNumber lb, ub;
00276 arglist_ [1] -> getBounds (lb, ub);
00277
00278 if (lb < .5)
00279 return false;
00280 }
00281
00282 return true;
00283 }
00284
00285
00287 void exprPow::closestFeasible (expression *varind,
00288 expression *vardep,
00289 CouNumber &left,
00290 CouNumber &right) const {
00291 CouNumber
00292 x = (*varind) (),
00293 y = (*vardep) (),
00294 k = arglist_ [1] -> Value (),
00295 xk = safe_pow (x, k),
00296 yk = safe_pow (y, 1./k);
00297
00298 int intk = 0;
00299
00300 bool isInt = fabs (k - (double) (intk = COUENNE_round (k))) < COUENNE_EPS,
00301 isInvInt = !isInt && (fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS);
00302
00303
00304
00305
00306
00307
00308 if (isInt || isInvInt)
00309
00310 if (intk % 2)
00311
00312 if (k > 0)
00313 ((y < xk) ? left : right) = yk;
00314
00315 else
00316
00317 if (y < 0.)
00318 if (y < xk) right = yk;
00319 else left = yk;
00320
00321 else
00322 if (y > xk) left = yk;
00323 else right = yk;
00324
00325 else
00326
00327 if (y <= 0.)
00328 left = - (right = COIN_DBL_MAX);
00329
00330 else
00331
00332 if (k > 0)
00333
00334 if (k < 1)
00335
00336 if (x > yk) left = yk;
00337 else right = yk;
00338
00339 else
00340
00341 if (x > yk) left = yk;
00342 else if (x < -yk) right = -yk;
00343 else left = - (right = yk);
00344
00345 else
00346 if (y < xk)
00347 left = - (right = yk);
00348 else
00349 if (x > 0) left = yk;
00350 else right = -yk;
00351
00352 else
00353
00354 if (k > 0) ((y < xk) ? left : right) = yk;
00355 else ((y > xk) ? left : right) = yk;
00356 }
00357
00358
00360 CouNumber exprPow::gradientNorm (const double *x) {
00361
00362 int ind0 = arglist_ [0] -> Index ();
00363 CouNumber exponent = arglist_ [1] -> Value ();
00364 return (ind0 < 0) ? 0. : fabs (exponent * safe_pow (x [ind0], exponent - 1));
00365 }
00366
00367
00370 bool exprPow::isCuttable (CouenneProblem *problem, int index) const {
00371
00372 CouNumber exponent = arglist_ [1] -> Value ();
00373
00374 bool
00375 isInt = ::isInteger (exponent),
00376 isInvInt = (exponent != 0.) && ::isInteger (1. / exponent);
00377
00378 int intExp = (isInt ? COUENNE_round (exponent) : (isInvInt ? COUENNE_round (1. / exponent) : 0));
00379
00380 if (exponent > 0.) {
00381
00382 if (isInt || isInvInt) {
00383
00384 if (intExp % 2) return false;
00385
00386 CouNumber
00387 x = problem -> X (arglist_ [0] -> Index ()),
00388 y = problem -> X (index);
00389
00390 if (isInt) return (y <= safe_pow (x, exponent));
00391
00392 return (y >= safe_pow (x, exponent));
00393 } else {
00394
00395
00396 CouNumber
00397 x = problem -> X (arglist_ [0] -> Index ()),
00398 y = problem -> X (index);
00399
00400 return (((exponent <= 1.) && (y >= safe_pow (x, exponent))) ||
00401 ((exponent >= 1.) && (y <= safe_pow (x, exponent))));
00402 }
00403 } else {
00404
00405
00406 CouNumber
00407 x = problem -> X (arglist_ [0] -> Index ()),
00408 y = problem -> X (index),
00409 lb = problem -> Lb (index),
00410 ub = problem -> Ub (index);
00411
00412 if (isInt || isInvInt)
00413
00414 if (!(intExp % 2)) return (((lb > 0) || (ub < 0)) && (y * safe_pow (fabs (x), -exponent) <= 1.));
00415 else return (((lb > 0) || (ub < 0)) && (y * safe_pow (x, -exponent) <= 1.));
00416 else return (y * safe_pow (x, -exponent) <= 1.);
00417 }
00418 }