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