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