00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "CoinHelperFunctions.hpp"
00012
00013 #include "CouenneExprPow.hpp"
00014 #include "CouenneObject.hpp"
00015 #include "CouenneBranchingObject.hpp"
00016 #include "CouenneProjections.hpp"
00017 #include "CouenneFunTriplets.hpp"
00018
00019 using namespace Couenne;
00020
00022 CouNumber negPowSelectBranch (const CouenneObject *obj,
00023 const OsiBranchingInformation *info,
00024 double * &brpts,
00025 double * &brDist,
00026
00027 int &way,
00028 CouNumber k,
00029 CouNumber x0, CouNumber y0,
00030 CouNumber l, CouNumber u);
00031
00032
00035 CouNumber exprPow::selectBranch (const CouenneObject *obj,
00036 const OsiBranchingInformation *info,
00037 expression *&var,
00038 double * &brpts,
00039 double * &brDist,
00040
00041 int &way) {
00042
00043
00044
00045
00046 var = arglist_ [0];
00047
00048 int
00049 ind = arglist_ [0] -> Index (),
00050 wi = obj -> Reference () -> Index ();
00051
00052 assert ((ind >= 0) && (wi >= 0) && (arglist_ [1] -> Type () == CONST));
00053
00054 double k = arglist_ [1] -> Value ();
00055
00056 CouNumber y0 = info -> solution_ [wi],
00057 x0 = info -> solution_ [ind],
00058 l = info -> lower_ [ind],
00059 u = info -> upper_ [ind];
00060
00061
00062
00063
00064
00065
00066
00067 if (x0 < l) x0 = l;
00068 else if (x0 > u) x0 = u;
00069
00070
00071 if (fabs (u-l) < COUENNE_EPS) {
00072 brpts = (CouNumber *) realloc (brpts, sizeof (CouNumber));
00073 *brpts = 0.5*(l+u);
00074 way = TWO_RAND;
00075 brDist = (double *) realloc (brDist, 2 * sizeof (double));
00076 return (brDist [0] = brDist [1] = fabs (y0 - pow (x0, k)));
00077 }
00078
00079
00080 if (k<0) return negPowSelectBranch (obj, info, brpts, brDist, way, k, x0, y0, l, u);
00081
00082 brDist = (double *) realloc (brDist, 2 * sizeof (double));
00083
00084 int intk = 0;
00085
00086 bool isInt = fabs (k - (double) (intk = COUENNE_round (k))) < COUENNE_EPS,
00087 isInvInt = !isInt && (k != 0.) && fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS;
00088
00089
00090
00091 if (isInt && !(intk % 2)) {
00092
00093 if ((l < - COUENNE_INFINITY) &&
00094 (u > COUENNE_INFINITY)) {
00095
00096 if (y0 < pow (x0, k)) {
00097
00098 brpts = (double *) realloc (brpts, sizeof (double));
00099 powertriplet pt (k);
00100 *brpts = powNewton (x0, y0, &pt);
00101
00102 way = (x0 > 0) ? TWO_RIGHT : TWO_LEFT;
00103
00104 x0 -= *brpts;
00105 y0 -= pow (*brpts, k);
00106
00107 return (brDist [0] = brDist [1] = sqrt (x0*x0 + y0*y0));
00108 }
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 brpts = (double *) realloc (brpts, sizeof (double));
00132 *brpts = x0;
00133 way = TWO_RAND;
00134
00135 return (brDist [0] = brDist [1] = fabs (y0 - pow (x0,k)));
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152 }
00153
00154
00155
00156 brpts = (double *) realloc (brpts, sizeof (double));
00157
00158 if (l < -COUENNE_INFINITY) {
00159
00160
00161 *brpts = obj -> midInterval (-safe_pow (y0, 1. / k), l, u, info);
00162 way = TWO_RIGHT;
00163
00164
00165
00166 return CoinMin (brDist [0] = x0 - *brpts,
00167 brDist [1] = projectSeg
00168 (x0, y0,
00169 *brpts, safe_pow (*brpts, k),
00170 u, safe_pow (u, k), -1));
00171 }
00172
00173 if (u > COUENNE_INFINITY) {
00174
00175
00176
00177 *brpts = obj -> midInterval (safe_pow (y0, 1. / k), l, u, info);
00178 way = TWO_LEFT;
00179
00180
00181
00182 return CoinMin (brDist [1] = *brpts - x0,
00183 brDist [0] = projectSeg
00184 (x0, y0,
00185 l, safe_pow (l, k),
00186 *brpts, safe_pow (*brpts, k), -1));
00187 }
00188
00189
00190
00191 powertriplet ft (k);
00192
00193 *brpts = obj -> getBrPoint (&ft, x0, l, u, info);
00194
00195 way = (x0 < *brpts) ? TWO_LEFT : TWO_RIGHT;
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 CouNumber retval =
00206 CoinMin (brDist [0] = projectSeg (x0,y0, l, safe_pow(l,k), *brpts, safe_pow(*brpts,k), 0),
00207 brDist [1] = projectSeg (x0,y0, *brpts, safe_pow(*brpts,k), u, safe_pow(u,k), 0));
00208
00209
00210
00211 return retval;
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221 }
00222
00223
00224
00225 brpts = (double *) realloc (brpts, sizeof (double));
00226 *brpts = x0;
00227 CouNumber pow0 = pow (x0, k);
00228
00229
00230
00231 if (isInt && (intk % 2)) {
00232
00233 way = (x0 > 0.) ? TWO_RIGHT : TWO_LEFT;
00234
00235 if (((l < - COUENNE_INFINITY) && (u > COUENNE_INFINITY)) ||
00236 ((l < - COUENNE_INFINITY) && (y0 < pow0)) ||
00237 ((u > COUENNE_INFINITY) && (y0 > pow0))) {
00238
00239 if (((y0 > 0) && (y0 < pow0)) ||
00240 ((y0 < 0) && (y0 > pow0))) {
00241
00242 *brpts = 0;
00243 return (brDist [0] = brDist [1] = fabs (pow0 - y0));
00244
00245 } else {
00246
00247 *brpts = pow (y0, 1./k);
00248
00249 return (brDist [0] = brDist [1] = (y0 > 0) ?
00250 projectSeg (x0, y0, x0, CoinMax (pow0, 0.), *brpts, y0, 0) :
00251 projectSeg (x0, y0, x0, CoinMin (pow0, 0.), *brpts, y0, 0));
00252 }
00253 }
00254
00255
00256
00257
00258 powertriplet pt (k);
00259 *brpts = obj -> getBrPoint (&pt, x0, l, u, info);
00260
00261
00262
00263
00264 if (y0 < pow0)
00265 *brpts = -*brpts;
00266
00267 if (*brpts > x0) {
00268
00269 brDist [0] = y0 - safe_pow (*brpts, k);
00270 brDist [1] = sqrt (brDist [0] * brDist [0] + (x0 - *brpts) * (x0 - *brpts));
00271
00272 } else {
00273
00274 brDist [1] = y0 - safe_pow (*brpts, k);
00275 brDist [0] = sqrt (brDist [1] * brDist [1] + (x0 - *brpts) * (x0 - *brpts));
00276 }
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287 }
00288
00289
00290
00291 if (isInvInt && (intk % 2)) {
00292
00293 way = (x0 > 0.) ? TWO_RIGHT : TWO_LEFT;
00294
00295 if (((l < - COUENNE_INFINITY) && (u > COUENNE_INFINITY)) ||
00296 ((l < - COUENNE_INFINITY) && (y0 < pow0)) ||
00297 ((u > COUENNE_INFINITY) && (y0 > pow0))) {
00298
00299 if (((x0 > 0.) && (y0 > pow0)) ||
00300 ((x0 < 0.) && (y0 < pow0))) {
00301
00302
00303
00304 *brpts = 0.;
00305 return (brDist [0] = brDist [1] = fabs (pow0 - y0));
00306
00307 } else {
00308
00309 *brpts = x0;
00310
00311
00312
00313 return (brDist [0] = brDist [1] = (x0 > 0) ?
00314 projectSeg (x0,y0,x0, pow0, CoinMax (0., pow (y0, 1./k)), y0, 0) :
00315 projectSeg (x0,y0,x0, pow0, CoinMin (0., pow (y0, 1./k)), y0, 0));
00316 }
00317 }
00318
00319
00320
00321
00322 powertriplet pt (k);
00323 *brpts = obj -> getBrPoint (&pt, x0, l, u, info);
00324
00325
00326
00327
00328 if (y0 > pow0) *brpts = -*brpts;
00329 else if (y0 < pow0) *brpts = -*brpts;
00330
00331 if (*brpts > x0) {
00332
00333 brDist [0] = y0 - safe_pow (*brpts, k);
00334 brDist [1] = sqrt (brDist [0] * brDist [0] + (x0 - *brpts) * (x0 - *brpts));
00335
00336 } else {
00337
00338 brDist [1] = y0 - safe_pow (*brpts, k);
00339 brDist [0] = sqrt (brDist [1] * brDist [1] + (x0 - *brpts) * (x0 - *brpts));
00340 }
00341
00342 return CoinMax (brDist [0], brDist [1]);
00343
00344
00345
00346
00347
00348
00349 }
00350
00351 if (k>1) {
00352
00353 if (y0 < pow0) {
00354
00355
00356 powertriplet pt (k);
00357 *brpts = powNewton (x0, y0, &pt);
00358
00359 way = TWO_LEFT;
00360
00361 x0 -= *brpts;
00362 y0 -= pow (*brpts, k);
00363
00364 return (brDist [0] = brDist [1] = sqrt (x0*x0 + y0*y0));
00365
00366 } else {
00367
00368
00369
00370
00371 if (obj -> Strategy () == CouenneObject::MID_INTERVAL)
00372 *brpts = 0.5 * (x0 + pow (y0, 1. / k));
00373 else {
00374 powertriplet pt (k);
00375 *brpts = obj -> getBrPoint (&pt, x0, l, u, info);
00376 }
00377
00378 way = TWO_LEFT;
00379
00380 if (l < 0.) l = 0.;
00381
00382 CouNumber
00383 powbpt = pow (*brpts, k),
00384 projL = projectSeg (x0, y0, l, pow (l, k), *brpts, powbpt, -1);
00385
00386 return (brDist[0] = brDist[1] = (u > COUENNE_INFINITY) ?
00387 CoinMin (projL, *brpts - x0) :
00388 CoinMin (projL, projectSeg (x0,y0, *brpts,powbpt, u, pow(u,k), -1)));
00389 }
00390
00391 } else {
00392
00393 if (y0 < pow0) {
00394
00395
00396
00397
00398
00399 if (obj -> Strategy () == CouenneObject::MID_INTERVAL)
00400 *brpts = 0.5 * (x0 + pow (y0, 1. / k));
00401 else {
00402 powertriplet pt (k);
00403 *brpts = obj -> getBrPoint (&pt, x0, l, u, info);
00404 }
00405
00406
00407 way = TWO_LEFT;
00408
00409 if (l < 0) l = 0;
00410
00411 CouNumber
00412 powbpt = pow (*brpts, k),
00413 projL = projectSeg (x0, y0, l, pow (l, k), *brpts, powbpt, +1);
00414
00415 return (brDist[0] = brDist[1] = (u > COUENNE_INFINITY) ?
00416 CoinMin (projL, powbpt - y0) :
00417 CoinMin (projL, projectSeg (x0, y0, *brpts, powbpt, u, pow (u,k), +1)));
00418
00419 } else {
00420
00421 powertriplet pt (k);
00422 *brpts = powNewton (x0, y0, &pt);
00423
00424 way = TWO_LEFT;
00425
00426 x0 -= *brpts;
00427 y0 -= pow (*brpts, k);
00428
00429 return (brDist [0] = brDist [1] = sqrt (x0*x0 + y0*y0));
00430 }
00431 }
00432
00433
00434
00435
00436 assert (false);
00437
00438 var = NULL;
00439 return 0.;
00440 }