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 "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
00033
00034 exprAux *exprPow::standardize (CouenneProblem *p, bool addAux) {
00035
00036 expression *ret;
00037
00038 if (arglist_ [0] -> Type () == CONST) {
00039
00040 exprOp::standardize (p);
00041
00042 CouNumber base = arglist_ [0] -> Value ();
00043
00044 if (fabs (base - M_E) < COUENNE_EPS)
00045 ret = new exprExp (new exprClone (arglist_ [1]));
00046 else
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) {
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 {
00062
00063
00064
00065 exprOp::standardize (p);
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083 {
00084
00085
00086
00087
00088
00089
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
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
00107
00108
00109 CouNumber k = arglist_ [1] -> Value ();
00110
00111
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
00129
00130
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
00146
00147 int intk = 0;
00148
00149 if (k < - COUENNE_INFINITY) {
00150 if (cL || cR) cg -> createCut (cs, 0., 0, w_ind, 1.);
00151 return;
00152 }
00153
00154 if (k > COUENNE_INFINITY)
00155 return;
00156
00157 if (fabs (k) < COUENNE_EPS) {
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
00166
00167 if ( (isInt || isInvInt)
00168 && (intk % 2)
00169 && (k > COUENNE_EPS)
00170 && (l < - COUENNE_EPS)
00171 && (u > COUENNE_EPS)) {
00172
00173
00174
00175
00176
00177
00178
00179
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
00208 CouNumber powThres = (k<=1) ? COUENNE_INFINITY : pow (COU_MAX_COEFF, 1./k);
00209
00210
00211 if (l > -powThres) {
00212 if (l>0) addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, u, sign);
00213 else if (u > q * l) {
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
00220 if (u < powThres) {
00221 if (u<0) addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, u, -sign);
00222 else if (l < q * u) {
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
00231
00232
00233
00234
00235
00236 if (!isInt
00237 && (!isInvInt || !(intk % 2))
00238 && (l < - COUENNE_EPS)
00239 && (u < (l=0)))
00240 return;
00241
00242
00243
00244
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
00257
00258
00259
00260
00261 int sign = 1;
00262
00263
00264 if ( ((l < - COUENNE_EPS) && (intk % 2) && (k < -COUENNE_EPS))
00265 || ((u < - COUENNE_EPS) && (intk % 2) && (k > COUENNE_EPS))
00266 || (fabs (k-0.5) < 0.5 - COUENNE_EPS))
00267 sign = -1;
00268
00269 CouNumber powThres = CoinMin (COUENNE_INFINITY,
00270 pow (COU_MAX_COEFF, 1./k)),
00271 powStep = 1;
00272
00273
00274 if ((k < COUENNE_EPS) &&
00275 isInt && !(intk % 2) &&
00276 (l < -COUENNE_EPS) &&
00277 (u > COUENNE_EPS) &&
00278 (l > - powThres) &&
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
00285 if (( (k > COUENNE_EPS)
00286 || (l > COUENNE_EPS)
00287 || (u < - COUENNE_EPS)) &&
00288 (l > - powThres) &&
00289 (u < powThres) &&
00290 (fabs (l+u) > COUENNE_EPS))
00291
00292 cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l, k), u, safe_pow (u, k), -sign);
00293
00294
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;
00309 if (fabs (u) < COUENNE_EPS) u = -1. / powThres;
00310 }
00311
00312 addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, u, sign);
00313 }
00314 }