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