00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <math.h>
00012 #ifndef M_E
00013 # define M_E 2.7182818284590452354
00014 #endif
00015
00016 #include "CouenneCutGenerator.hpp"
00017
00018 #include "CouenneTypes.hpp"
00019 #include "CouenneRootQ.hpp"
00020 #include "CouenneExprAux.hpp"
00021 #include "CouenneExprPow.hpp"
00022 #include "CouenneExprExp.hpp"
00023 #include "CouenneExprConst.hpp"
00024 #include "CouenneExprClone.hpp"
00025 #include "CouenneExprMul.hpp"
00026 #include "CouenneExprSum.hpp"
00027 #include "CouenneExprLog.hpp"
00028 #include "CouennePrecisions.hpp"
00029 #include "CouenneProblem.hpp"
00030
00031 using namespace Couenne;
00032
00033 std::map <int, CouNumber> Qroot::Qmap;
00034
00035
00036
00037 exprAux *exprPow::standardize (CouenneProblem *p, bool addAux) {
00038
00039 expression *ret;
00040
00041 if (arglist_ [0] -> Type () == CONST) {
00042
00043 exprOp::standardize (p);
00044
00045 CouNumber base = arglist_ [0] -> Value ();
00046
00047 if (fabs (base - M_E) < COUENNE_EPS_SIMPL)
00048 ret = new exprExp (new exprClone (arglist_ [1]));
00049 else
00050 ret = new exprExp (new exprClone
00051 (p -> addAuxiliary (new exprMul (new exprClone (arglist_ [1]),
00052 new exprConst (log (base))))));
00053 } else if (arglist_ [1] -> Type () != CONST) {
00054
00055 exprOp::standardize (p);
00056
00057 ret = new exprExp (new exprClone (p -> addAuxiliary
00058 (new exprMul
00059 (new exprClone (arglist_ [1]),
00060 new exprClone
00061 (p -> addAuxiliary
00062 (new exprLog (new exprClone (arglist_ [0]))))))));
00063
00064 } else {
00065
00066
00067
00068 exprOp::standardize (p);
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 {
00087
00088
00089
00090
00091
00092
00093
00094 return (addAux ? (p -> addAuxiliary (this)) : new exprAux (this, p -> domain ()));
00095 }
00096 }
00097
00098 return (addAux ? (p -> addAuxiliary (ret)) : new exprAux (ret, p -> domain ()));
00099 }
00100
00101
00102
00103
00104 void exprPow::generateCuts (expression *aux,
00105 OsiCuts &cs, const CouenneCutGenerator *cg,
00106 t_chg_bounds *chg, int wind,
00107 CouNumber lbw, CouNumber ubw) {
00108
00109
00110
00111
00112 CouNumber k = arglist_ [1] -> Value ();
00113
00114
00115
00116 expression *xe = arglist_ [0];
00117
00118 CouNumber l, u;
00119 xe -> getBounds (l, u);
00120
00121 int w_ind = aux -> Index ();
00122 int x_ind = xe -> Index ();
00123
00124 bool cL = !chg || (chg [x_ind].lower() != t_chg_bounds::UNCHANGED) || cg -> isFirst ();
00125 bool cR = !chg || (chg [x_ind].upper() != t_chg_bounds::UNCHANGED) || cg -> isFirst ();
00126
00127 CouNumber
00128 w = (*aux) (),
00129 x = (*xe) ();
00130
00131 enum auxSign aSign = cg -> Problem () -> Var (w_ind) -> sign ();
00132
00133
00134
00135
00136
00137 if (fabs (u-l) < COUENNE_EPS) {
00138
00139 CouNumber avg = 0.5 * (l+u),
00140 avg_k_1 = safe_pow (avg, k-1),
00141 lk = safe_pow (l, k),
00142 uk = safe_pow (u, k);
00143
00144 if (cL || cR)
00145 cg -> createCut (cs, u*lk - l*uk + avg * avg_k_1 * (1-k), aSign,
00146 w_ind, u - l + 1, x_ind, lk-uk - k * avg_k_1);
00147 return;
00148 }
00149
00150
00151
00152 int intk = 0;
00153
00154 if (k < - COUENNE_INFINITY) {
00155 if (cL || cR) cg -> createCut (cs, 0., aSign, w_ind, 1.);
00156 return;
00157 }
00158
00159 if (k > COUENNE_INFINITY)
00160 return;
00161
00162 if (fabs (k) < COUENNE_EPS) {
00163 if (cL || cR) cg -> createCut (cs, 1., aSign, w_ind, 1.);
00164 return;
00165 }
00166
00167 bool isInt = fabs (k - (double) (intk = COUENNE_round (k))) < COUENNE_EPS,
00168 isInvInt = !isInt && (fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS);
00169
00170
00171
00172 if ( (isInt || isInvInt)
00173 && (intk % 2)
00174 && (k > COUENNE_EPS)
00175
00176
00177 ) {
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187 CouNumber q = 0.;
00188
00189 if ((l<0.) && (u>0.)) {
00190
00191 Qroot qmap;
00192 q = qmap (intk);
00193 }
00194
00195 int sign;
00196
00197 if (isInvInt) {
00198 if (cg -> isFirst ()) {
00199 w = (l>0.) ? 1. : (u<0.) ? -1. : 0.;
00200 x = 0;
00201 }
00202 q = safe_pow (q, k);
00203 sign = -1;
00204 }
00205 else {
00206 if (cg -> isFirst ()) {
00207 x = (l>0) ? l : (u<0) ? u : 0;
00208 w = 0;
00209 }
00210 sign = 1;
00211 }
00212
00213
00214 CouNumber powThres = (k<=1) ? COUENNE_INFINITY : pow (COU_MAX_COEFF, 1./k);
00215
00216
00217 if ((aSign != expression::AUX_LEQ) && (l > -powThres)) {
00218 if (l>=0.) addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, u, sign);
00219 else if (u > q * l) {
00220 addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, q*l, u, sign);
00221 cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l,k), q*l, safe_pow (q*l,k), sign);
00222 } else {
00223 cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l,k), u, safe_pow (u, k), sign);
00224 }
00225 }
00226
00227
00228 if ((aSign != expression::AUX_GEQ) && (u < powThres)) {
00229 if (u<=0.) addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, u, -sign);
00230 else if (l < q * u) {
00231 addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, q*u, -sign);
00232 cg -> addSegment (cs, w_ind, x_ind, q*u, safe_pow (q*u,k), u, safe_pow (u,k), -sign);
00233 } else {
00234 cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l,k), u, safe_pow (u,k), -sign);
00235 }
00236 }
00237 }
00238 else {
00239
00240
00241
00242
00243
00244
00245
00246 if (!isInt
00247 && (!isInvInt || !(intk % 2))
00248 && (l < - COUENNE_EPS)
00249 && (u < (l=0)))
00250 return;
00251
00252
00253
00254
00255
00256 if ((k < 0) &&
00257 (l < - COUENNE_EPS) &&
00258 (u > COUENNE_EPS) &&
00259 aSign != expression::AUX_LEQ) {
00260
00261 if (!(intk % 2))
00262 cg -> addSegment (cs, w_ind, arglist_ [0] -> Index (),
00263 l, safe_pow (l,k), u, safe_pow (u,k), +1);
00264
00265
00266
00267
00268 return;
00269 }
00270
00271
00272
00273
00274
00275
00276 int sign = 1;
00277
00278
00279 if ( ((l < - COUENNE_EPS) && (intk % 2) && (k < -COUENNE_EPS))
00280 || ((u <= 0.) && (intk % 2) && (k > COUENNE_EPS))
00281 || (fabs (k-0.5) < 0.5 - COUENNE_EPS))
00282 sign = -1;
00283
00284 CouNumber powThres = CoinMin (COUENNE_INFINITY,
00285 pow (COU_MAX_COEFF, 1./k)),
00286 powStep = 1;
00287
00288
00289 if ((k < COUENNE_EPS) &&
00290 isInt && !(intk % 2) &&
00291 (l < -COUENNE_EPS) &&
00292 (u > COUENNE_EPS) &&
00293 (l > - powThres) &&
00294 (u < powThres) &&
00295 aSign != expression::AUX_LEQ)
00296
00297 cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l, k), u, safe_pow (u, k), 1);
00298
00299
00300 if (( (k > COUENNE_EPS)
00301 || (l > COUENNE_EPS)
00302 || (u < - COUENNE_EPS)) &&
00303 (l > - powThres) &&
00304 (u < powThres) &&
00305 (fabs (l+u) > COUENNE_EPS) &&
00306 (aSign != expression::AUX_GEQ))
00307
00308 cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l, k), u, safe_pow (u, k), -sign);
00309
00310
00311
00312 if (cg -> isFirst()) {
00313 x = (k<0) ? ((u<0) ? u : (l>0) ? l : 0) : 0;
00314 w = 0;
00315 }
00316
00317 if (k > COUENNE_EPS) {
00318
00319 if (u > powThres) u = CoinMax (x,l) + powStep;
00320 if (l < - powThres) l = CoinMin (x,u) - powStep;
00321 }
00322 else {
00323
00324 if (fabs (l) < COUENNE_EPS) l = 1. / powThres;
00325 if (fabs (u) < COUENNE_EPS) u = -1. / powThres;
00326 }
00327
00328 if ((sign == 1 && aSign != expression::AUX_LEQ) ||
00329 (sign == -1 && aSign != expression::AUX_GEQ))
00330 addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, u, sign);
00331 }
00332 }