/home/coin/SVN-release/OS-2.0.0/Couenne/src/branch/CouenneObject.cpp

Go to the documentation of this file.
00001 /* $Id: CouenneObject.cpp 141 2009-06-03 04:19:19Z pbelotti $ */
00002 /*
00003  * Name:    CouenneObject.cpp
00004  * Authors: Pierre Bonami, IBM Corp.
00005  *          Pietro Belotti, Carnegie Mellon University
00006  * Purpose: Base object for variables (to be used in branching)
00007  *
00008  * (C) Carnegie-Mellon University, 2006-09.
00009  * This file is licensed under the Common Public License (CPL)
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   // read options
00061   setParameters (base);
00062 
00063   // debug output ////////////////////////////////////////////
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   // read options
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   // a nonlinear constraint w = f(x) is violated. The infeasibility
00135   // should be given by something more elaborate than |w-f(x)|, that
00136   // is, it is the minimum, among the two branching nodes, of the
00137   // distance from the current optimum (w,x) and the optimum obtained
00138   // after convexifying the two subproblems. We call selectBranch for
00139   // the purpose, and save the output parameter into the branching
00140   // point that should be used later in createBranch.
00141 
00142   if (jnlst_ -> ProduceOutput (J_ITERSUMMARY, J_BRANCHING)) {
00143     printf ("CouObj::createBranch on ");
00144     reference_ -> print (); printf ("\n");
00145   }
00146 
00147   // copy current point into Couenne
00148   problem_ -> domain () -> push 
00149     (problem_ -> nVars (),
00150      info -> solution_,
00151      info -> lower_,
00152      info -> upper_);
00153 
00154   CouNumber  
00155     *brPts  = NULL,         // branching point(s)
00156     *brDist = NULL;         // distances from current LP point to each
00157                             // new convexification (usually two)
00158   expression *brVar = NULL; // branching variable
00159   int whichWay = 0;
00160 
00161   CouNumber improv;
00162 
00163   if (reference_ -> Image ())  // not used out of debug
00164     improv = reference_ -> Image () -> 
00165       selectBranch (this, info,                      // input parameters
00166                     brVar, brPts, brDist, whichWay); // result: who, where, distance, and direction
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); // MUST have a branching variable
00185 
00186   /*  if        (pseudoMultType_ == INTERVAL) {
00187     int index = brVar -> Index ();
00188     downEstimate_ = *brPts - info -> lower_ [index];
00189     upEstimate_   =          info -> upper_ [index] - *brPts;
00190   } else if (pseudoMultType_ == REV_INTERVAL) {
00191     int index = brVar -> Index ();
00192     downEstimate_ = *brPts - info -> upper_ [index];          // notice upper&lower inverted
00193     upEstimate_   =          info -> lower_ [index] - *brPts;
00194   } else 
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   // create branching object /////////////////////////////////////// 
00231   OsiBranchingObject *brObj = new CouenneBranchingObject 
00232     (si, this, jnlst_, brVar, way, *brPts, doFBBT_, doConvCuts_);
00233 
00234   problem_ -> domain () -> pop (); // Couenne discards current point
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; // 0.                                    // ]-inf,+inf[
00254     else                            return ((x < -COUENNE_EPS) ? (AGGR_MUL * (-1+x)) : // ]-inf,u]
00255                                             (x >  COUENNE_EPS) ? 0. : -AGGR_MUL);
00256   else
00257     if (u >  COUENNE_INFINITY / 10) return ((x >  COUENNE_EPS) ? (AGGR_MUL *  (1+x)) : // [l,+inf[
00258                                             (x < -COUENNE_EPS) ? 0. :  AGGR_MUL);
00259     else {                                                                             // [l,u]
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   // check if fval is a number (happens with e.g. w13 = w12/w5 and w5=0, see test/harker.nl)
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   //printf ("checkinf --> v=%e f=%e den=%e ret=%e ratio=%e\n", vval, fval, denom, retval, ratio);
00359 
00360   if ((ratio < 2)  && 
00361       (ratio > .5) &&
00362       ((retval /= denom) < CoinMin (COUENNE_EPS, feas_tolerance_)))
00363     retval = 0.;
00364 
00365   if (//(retval > 0.) &&
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   // Need to allow infeasibility for a variable without making the
00375   // whole problem infeasible with an infeasibility = 1e+50. Check
00376   // BonChooseVariable.cpp:382
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   // accept options for branching rules specific to each operator
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       // begin with default value in case specific exponent are not given
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       // read option
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   // now move it away from the bounds
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: // taken care of in selectBranch procedure
00541     break;
00542 
00543   default: 
00544     printf ("Couenne: invalid estimate setting procedure\n");
00545     exit (-1);
00546   }
00547 
00548 //   if (downEstimate_ <= 0.0 || upEstimate_ <= 0.0)
00549 //     printf ("%g [%g,%g] ---> [%g,%g]\n", 
00550 //          point,
00551 //          info -> lower_ [index],
00552 //          info -> upper_ [index],
00553 //          downEstimate_, upEstimate_);
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 }

Generated on Mon Aug 3 03:02:19 2009 by  doxygen 1.4.7