00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "CoinHelperFunctions.hpp"
00013 #include "CoinFinite.hpp"
00014
00015 #include "CouenneProblem.hpp"
00016 #include "CouenneObject.hpp"
00017 #include "CouenneBranchingObject.hpp"
00018
00019 const CouNumber default_alpha = 0.2;
00020 const CouNumber default_clamp = 0.2;
00021 const CouNumber max_pseudocost = 1000.;
00022
00024 CouenneObject::CouenneObject ():
00025
00026 OsiObject (),
00027 cutGen_ (NULL),
00028 problem_ (NULL),
00029 reference_ (NULL),
00030 strategy_ (MID_INTERVAL),
00031 jnlst_ (NULL),
00032 alpha_ (default_alpha),
00033 lp_clamp_ (default_clamp),
00034 feas_tolerance_ (feas_tolerance_default),
00035 doFBBT_ (true),
00036 doConvCuts_ (true),
00037 downEstimate_ (0.),
00038 upEstimate_ (0.),
00039 pseudoMultType_ (INFEASIBILITY) {}
00040
00041
00043 CouenneObject::CouenneObject (CouenneCutGenerator *cutgen,
00044 CouenneProblem *p,
00045 exprVar *ref,
00046 Bonmin::BabSetupBase *base,
00047 JnlstPtr jnlst):
00048 OsiObject (),
00049 cutGen_ (cutgen),
00050 problem_ (p),
00051 reference_ (ref),
00052 strategy_ (MID_INTERVAL),
00053 jnlst_ (jnlst),
00054 alpha_ (default_alpha),
00055 lp_clamp_ (default_clamp),
00056 feas_tolerance_ (feas_tolerance_default),
00057 doFBBT_ (true),
00058 doConvCuts_ (true),
00059 downEstimate_ (0.),
00060 upEstimate_ (0.),
00061 pseudoMultType_ (INFEASIBILITY) {
00062
00063
00064 setParameters (base);
00065
00066
00067
00068 if (reference_ &&
00069 (reference_ -> Type () == AUX) &&
00070 jnlst_ -> ProduceOutput (J_SUMMARY, J_BRANCHING)) {
00071
00072 printf ("created Expression Object: "); reference_ -> print ();
00073 if (reference_ -> Image ()) {
00074 printf (" := ");
00075 reference_ -> Image () -> print ();
00076 }
00077
00078 printf (" with %s strategy [clamp=%g, alpha=%g]\n",
00079 (strategy_ == LP_CLAMPED) ? "lp-clamped" :
00080 (strategy_ == LP_CENTRAL) ? "lp-central" :
00081 (strategy_ == BALANCED) ? "balanced" :
00082 (strategy_ == MIN_AREA) ? "min-area" :
00083 (strategy_ == MID_INTERVAL) ? "mid-point" :
00084 (strategy_ == NO_BRANCH) ? "no-branching (null infeasibility)" :
00085 "no strategy",
00086 lp_clamp_, alpha_);
00087 }
00088 }
00089
00090
00092 CouenneObject::CouenneObject (exprVar *ref,
00093 Bonmin::BabSetupBase *base,
00094 JnlstPtr jnlst):
00095
00096 OsiObject (),
00097 cutGen_ (NULL),
00098 problem_ (NULL),
00099 reference_ (ref),
00100 strategy_ (MID_INTERVAL),
00101 jnlst_ (jnlst),
00102 alpha_ (default_alpha),
00103 lp_clamp_ (default_clamp),
00104 feas_tolerance_ (feas_tolerance_default),
00105 doFBBT_ (true),
00106 doConvCuts_ (true),
00107 downEstimate_ (0.),
00108 upEstimate_ (0.),
00109 pseudoMultType_ (INFEASIBILITY) {
00110
00111
00112 setParameters (base);
00113 }
00114
00115
00117 CouenneObject::CouenneObject (const CouenneObject &src):
00118 OsiObject (src),
00119 cutGen_ (src.cutGen_),
00120 problem_ (src.problem_),
00121 reference_ (src.reference_),
00122 strategy_ (src.strategy_),
00123 jnlst_ (src.jnlst_),
00124 alpha_ (src.alpha_),
00125 lp_clamp_ (src.lp_clamp_),
00126 feas_tolerance_ (src.feas_tolerance_),
00127 doFBBT_ (src.doFBBT_),
00128 doConvCuts_ (src.doConvCuts_),
00129 downEstimate_ (src.downEstimate_),
00130 upEstimate_ (src.upEstimate_),
00131 pseudoMultType_ (src.pseudoMultType_) {}
00132
00133
00135 OsiBranchingObject *CouenneObject::createBranch (OsiSolverInterface *si,
00136 const OsiBranchingInformation *info,
00137 int way) const {
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147 if (jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_BRANCHING)) {
00148 printf ("CouObj::createBranch on ");
00149 reference_ -> print (); printf ("\n");
00150 }
00151
00152
00153 problem_ -> domain () -> push
00154 (problem_ -> nVars (),
00155 info -> solution_,
00156 info -> lower_,
00157 info -> upper_);
00158
00159 CouNumber
00160 *brPts = NULL,
00161 *brDist = NULL;
00162
00163 expression *brVar = NULL;
00164 int whichWay = 0;
00165
00166 CouNumber improv;
00167
00168 if (reference_ -> Image ())
00169 improv = reference_ -> Image () ->
00170 selectBranch (this, info,
00171 brVar, brPts, brDist, whichWay);
00172 else {
00173 brVar = reference_;
00174 brPts = (double *) realloc (brPts, sizeof (double));
00175 brDist = (double *) realloc (brDist, 2 * sizeof (double));
00176
00177 double point = info -> solution_ [reference_ -> Index ()];
00178
00179 *brPts = point;
00180 improv = 0.;
00181
00182 if (point > floor (point)) {improv = brDist [0] = point - floor (point);}
00183 if (point < ceil (point)) {improv = CoinMin (improv, brDist [1] = ceil (point) - point);}
00184
00185 point -= floor (point);
00186 whichWay = (point < 0.45) ? TWO_LEFT : (point > 0.55) ? TWO_RIGHT : TWO_RAND;
00187 }
00188
00189 assert (brVar);
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202 if (pseudoMultType_ == PROJECTDIST) {
00203 downEstimate_ = brDist [0];
00204 upEstimate_ = brDist [1];
00205 } else setEstimates (info, NULL, brPts);
00206
00208 if (jnlst_ -> ProduceOutput (J_MOREMATRIX, J_BRANCHING)) {
00209 printf ("brpts for "); reference_ -> print ();
00210 if (reference_ -> Image ()) {printf (" := "); reference_ -> Image () -> print ();}
00211 printf (" is on "); brVar -> print ();
00212 printf (" @ %.12g [%.12g,%.12g]\n", *brPts,
00213 problem_ -> Lb (brVar -> Index ()),
00214 problem_ -> Ub (brVar -> Index ()));
00215
00216 if (brVar) {
00217
00218 if (improv <= COUENNE_EPS) {
00219 printf ("### warning, infeas = %g for ", improv);
00220 reference_ -> print ();
00221 if (reference_ -> Image ()) {printf (":="); reference_ -> Image () -> print ();}
00222 printf ("\n");
00223 }
00224
00225 int index = brVar -> Index ();
00226 if (info -> lower_ [index] >=
00227 info -> upper_ [index] - COUENNE_EPS) {
00228 printf ("### warning, tiny bounding box [%g,%g] for x_%d\n",
00229 info -> lower_ [index],
00230 info -> upper_ [index], index);
00231 }
00232 }
00233 }
00234
00235
00236 OsiBranchingObject *brObj = new CouenneBranchingObject
00237 (si, this, jnlst_, cutGen_, problem_, brVar, way, *brPts, doFBBT_, doConvCuts_);
00238
00239 problem_ -> domain () -> pop ();
00240
00241 if (brPts) free (brPts);
00242 if (brDist) free (brDist);
00243
00244 return brObj;
00245 }
00246
00247
00249 CouNumber CouenneObject::midInterval (CouNumber x, CouNumber l, CouNumber u) const {
00250
00251 if (u < l + COUENNE_EPS)
00252 return (0.5 * (l + u));
00253
00254 if (x<l) x = l;
00255 else if (x>u) x = u;
00256
00257 if (l < -COUENNE_INFINITY / 10)
00258 if (u > COUENNE_INFINITY / 10) return x;
00259 else return ((x < -COUENNE_EPS) ? (AGGR_MUL * (-1+x)) :
00260 (x > COUENNE_EPS) ? 0. : -AGGR_MUL);
00261 else
00262 if (u > COUENNE_INFINITY / 10) return ((x > COUENNE_EPS) ? (AGGR_MUL * (1+x)) :
00263 (x < -COUENNE_EPS) ? 0. : AGGR_MUL);
00264 else {
00265 CouNumber point = alpha_ * x + (1. - alpha_) * (l + u) / 2.;
00266 if ((point-l) / (u-l) < closeToBounds) point = l + (u-l) * closeToBounds;
00267 else if ((u-point) / (u-l) < closeToBounds) point = u + (l-u) * closeToBounds;
00268 return point;
00269 }
00270 }
00271
00272
00274 CouNumber CouenneObject::getBrPoint (funtriplet *ft, CouNumber x0, CouNumber l, CouNumber u) const {
00275
00276 if ((l < -COUENNE_EPS) &&
00277 (u > COUENNE_EPS) &&
00278 (-l/u >= THRES_ZERO_SYMM) &&
00279 (-u/l >= THRES_ZERO_SYMM))
00280 return 0.;
00281
00282 CouNumber width = lp_clamp_ * (u-l);
00283
00284 switch (strategy_) {
00285
00286 case CouenneObject::MIN_AREA: return maxHeight (ft, l, u);
00287 case CouenneObject::BALANCED: return minMaxDelta (ft, l, u);
00288 case CouenneObject::LP_CLAMPED: return CoinMax (l + width, CoinMin (x0, u - width));
00289 case CouenneObject::LP_CENTRAL: return ((x0 < l + width) || (x0 > u - width)) ? (l+u)/2 : x0;
00290 case CouenneObject::MID_INTERVAL: return midInterval (x0, l, u);
00291 default:
00292 printf ("Couenne: unknown branching point selection strategy\n");
00293 exit (-1);
00294 }
00295 }
00296
00297
00299 double CouenneObject::infeasibility (const OsiBranchingInformation *info, int &way) const {
00300
00301 if (strategy_ == NO_BRANCH)
00302 return (upEstimate_ = downEstimate_ = 0.);
00303
00304 problem_ -> domain () -> push
00305 (problem_ -> nVars (),
00306 info -> solution_,
00307 info -> lower_,
00308 info -> upper_);
00309
00310 double retval = checkInfeasibility (info);
00311
00312 problem_ -> domain () -> pop ();
00313
00314 if (pseudoMultType_ == INFEASIBILITY) {
00315
00316 double point = info -> solution_ [reference_ -> Index ()];
00317
00318 if (reference_ -> isInteger ()) {
00319 if (retval < intInfeasibility (point)) {
00320 if (downEstimate_ < point - floor (point)) downEstimate_ = point - floor (point);
00321 if (upEstimate_ < ceil (point) - point) upEstimate_ = ceil (point) - point;
00322 retval = intInfeasibility (point);
00323 }
00324 }
00325 else upEstimate_ = downEstimate_ = retval;
00326 }
00327 else setEstimates (info, &retval, NULL);
00328
00329 return (reference_ -> isInteger ()) ?
00330 CoinMax (retval, intInfeasibility (info -> solution_ [reference_ -> Index ()])) :
00331 retval;
00332 }
00333
00334
00338 double CouenneObject::checkInfeasibility (const OsiBranchingInformation *info) const {
00339
00340 if (reference_ -> Type () == VAR)
00341 return (reference_ -> isInteger ()) ?
00342 intInfeasibility (info -> solution_ [reference_ -> Index ()]) : 0.;
00343
00344 double
00345 vval = info -> solution_ [reference_ -> Index ()],
00346 fval = (*(reference_ -> Image ())) (),
00347 denom = CoinMax (1., reference_ -> Image () -> gradientNorm (info -> solution_));
00348
00349
00350 if (CoinIsnan (fval)) {
00351 fval = vval + 1.;
00352 denom = 1.;
00353 }
00354
00355 if (fabs (fval) > COUENNE_INFINITY)
00356 fval = COUENNE_INFINITY;
00357
00358 double
00359 retval = fabs (vval - fval),
00360 ratio = (CoinMax (1., fabs (vval)) /
00361 CoinMax (1., fabs (fval)));
00362
00363
00364
00365 if ((ratio < 2) &&
00366 (ratio > .5) &&
00367 ((retval /= denom) < CoinMin (COUENNE_EPS, feas_tolerance_)))
00368 retval = 0.;
00369
00370 if (
00371 (jnlst_ -> ProduceOutput (J_DETAILED, J_BRANCHING))) {
00372
00373 printf (" infeas %g: ", retval);
00374 reference_ -> print ();
00375 if (reference_ -> Image ()) {printf (" := "); reference_ -> Image () -> print ();}
00376 printf ("\n");
00377 }
00378
00379
00380
00381
00382
00383 if (retval > 1.e40)
00384 retval = 1.e20;
00385
00386 return (reference_ -> isInteger ()) ?
00387 CoinMax (retval, intInfeasibility (info -> solution_ [reference_ -> Index ()])) :
00388 retval;
00389 }
00390
00391
00393 void CouenneObject::setParameters (Bonmin::BabSetupBase *base) {
00394
00395 if (!base) return;
00396
00397 std::string s;
00398
00399 base -> options () -> GetStringValue ("pseudocost_mult", s, "couenne.");
00400
00401 if (s == "interval_lp") pseudoMultType_ = INTERVAL_LP;
00402 else if (s == "interval_lp_rev") pseudoMultType_ = INTERVAL_LP_REV;
00403 else if (s == "interval_br") pseudoMultType_ = INTERVAL_BR;
00404 else if (s == "interval_br_rev") pseudoMultType_ = INTERVAL_BR_REV;
00405 else if (s == "infeasibility") pseudoMultType_ = INFEASIBILITY;
00406 else if (s == "projectDist") pseudoMultType_ = PROJECTDIST;
00407
00408 base -> options() -> GetStringValue ("branch_fbbt", s, "couenne."); doFBBT_ = (s=="yes");
00409 base -> options() -> GetStringValue ("branch_conv_cuts", s, "couenne."); doConvCuts_ = (s=="yes");
00410
00411 base -> options() -> GetNumericValue ("feas_tolerance", feas_tolerance_, "couenne.");
00412
00413 std::string brtype;
00414 base -> options () -> GetStringValue ("branch_pt_select", brtype, "couenne.");
00415
00416 if (brtype == "balanced") strategy_ = BALANCED;
00417 else if (brtype == "lp-clamped") strategy_ = LP_CLAMPED;
00418 else if (brtype == "lp-central") strategy_ = LP_CENTRAL;
00419 else if (brtype == "min-area") strategy_ = MIN_AREA;
00420 else if (brtype == "no-branch") strategy_ = NO_BRANCH;
00421 else if (brtype == "mid-point") {
00422 strategy_ = MID_INTERVAL;
00423 base -> options () -> GetNumericValue ("branch_midpoint_alpha", alpha_, "couenne.");
00424 }
00425
00426
00427 if (reference_ && reference_ -> Type () == AUX) {
00428
00429 std::string br_operator = "";
00430
00431 switch (reference_ -> Image () -> code ()) {
00432
00433 case COU_EXPRPOW: {
00434
00435
00436 base -> options () -> GetStringValue ("branch_pt_select_pow", brtype, "couenne.");
00437
00438 CouNumber expon = reference_ -> Image () -> ArgList () [1] -> Value ();
00439
00440 if (fabs (expon - 2.) < COUENNE_EPS) br_operator = "sqr";
00441 else if (fabs (expon - 3.) < COUENNE_EPS) br_operator = "cube";
00442 else if (expon < 0.) br_operator = "negpow";
00443 else br_operator = "pow";
00444 } break;
00445
00446 case COU_EXPRMUL:
00447 br_operator = (reference_ -> Image () -> ArgList () [0] -> Index () !=
00448 reference_ -> Image () -> ArgList () [1] -> Index ()) ?
00449 "prod" : "sqr";
00450 break;
00451 case COU_EXPRINV: br_operator = "negpow"; break;
00452 case COU_EXPRDIV: br_operator = "div"; break;
00453 case COU_EXPRLOG: br_operator = "log"; break;
00454 case COU_EXPREXP: br_operator = "exp"; break;
00455 case COU_EXPRSIN:
00456 case COU_EXPRCOS: br_operator = "trig"; break;
00457 default: break;
00458 }
00459
00460 if (br_operator != "") {
00461
00462 char select [40], sel_clamp [40];
00463 sprintf (select, "branch_pt_select_%s", br_operator.c_str ());
00464 sprintf (sel_clamp, "branch_lp_clamp_%s", br_operator.c_str ());
00465 base -> options () -> GetStringValue (select, brtype, "couenne.");
00466 base -> options () -> GetNumericValue (sel_clamp, lp_clamp_, "couenne.");
00467
00468 if (brtype == "balanced") strategy_ = BALANCED;
00469 else if (brtype == "lp-clamped") strategy_ = LP_CLAMPED;
00470 else if (brtype == "lp-central") strategy_ = LP_CENTRAL;
00471 else if (brtype == "min-area") strategy_ = MIN_AREA;
00472 else if (brtype == "no-branch") strategy_ = NO_BRANCH;
00473 else if (brtype == "mid-point") {
00474 strategy_ = MID_INTERVAL;
00475 base -> options () -> GetNumericValue ("branch_midpoint_alpha", alpha_, "couenne.");
00476 }
00477 }
00478 }
00479 }
00480
00481
00483 void CouenneObject::setEstimates (const OsiBranchingInformation *info,
00484 CouNumber *infeasibility,
00485 CouNumber *brpoint) const {
00486
00487 int index = reference_ -> Index ();
00488
00489 CouNumber
00490 *up = &upEstimate_,
00491 *down = &downEstimate_,
00492 point = 0.,
00493 lower = info -> lower_ [index],
00494 upper = info -> upper_ [index];
00495
00497
00498
00499
00500 if ((pseudoMultType_ == INTERVAL_LP_REV) ||
00501 (pseudoMultType_ == INTERVAL_BR_REV)) {
00502
00503 up = &downEstimate_;
00504 down = &upEstimate_;
00505 }
00506
00508 if (info &&
00509 ((pseudoMultType_ == INTERVAL_LP) ||
00510 (pseudoMultType_ == INTERVAL_LP_REV)))
00511
00512 point = info -> solution_ [index];
00513
00514 else if (brpoint &&
00515 ((pseudoMultType_ == INTERVAL_BR) ||
00516 (pseudoMultType_ == INTERVAL_BR_REV)))
00517
00518 point = *brpoint;
00519
00520
00521
00522 point = midInterval (point, lower, upper);
00523
00524 if ((lower > -COUENNE_INFINITY) &&
00525 (upper < COUENNE_INFINITY)) {
00526
00527 CouNumber delta = closeToBounds * (upper - lower);
00528
00529 if (point < lower + delta)
00530 point = lower + delta;
00531 else if (point > upper - delta)
00532 point = upper - delta;
00533 }
00534
00536 switch (pseudoMultType_) {
00537
00538 case INFEASIBILITY:
00539 if (infeasibility)
00540 upEstimate_ = downEstimate_ = *infeasibility;
00541 break;
00542
00543 case INTERVAL_LP:
00544 case INTERVAL_LP_REV:
00545 case INTERVAL_BR:
00546 case INTERVAL_BR_REV:
00547 assert (info);
00548 *up = CoinMin (max_pseudocost, upper - point);
00549 *down = CoinMin (max_pseudocost, point - lower);
00550 break;
00551
00552 case PROJECTDIST:
00553 break;
00554
00555 default:
00556 printf ("Couenne: invalid estimate setting procedure\n");
00557 exit (-1);
00558 }
00559
00560
00561
00562
00563
00564
00565
00566
00567 if (reference_ -> isInteger ()) {
00568
00569 CouNumber
00570 fracDn = point - floor (point),
00571 fracUp = ceil (point) - point;
00572
00573 if (downEstimate_ < fracDn) downEstimate_ = fracDn;
00574 if (upEstimate_ < fracUp) upEstimate_ = fracUp;
00575 }
00576
00577 assert (downEstimate_ >= 0. &&
00578 upEstimate_ >= 0.);
00579 }