# Changeset 556

Ignore:
Timestamp:
May 28, 2007 6:59:50 PM (12 years ago)
Message:

set OBBT (optimality based bound tightening) as default. Use NLP solution, when available, to add convexification cuts. Created specialized branching rules for divisions. Added extra cuts for products with positive lower- or negative upper bounds

Location:
branches/Couenne/Couenne/src
Files:
16 edited

Unmodified
Removed
• ## branches/Couenne/Couenne/src/branch/CouenneBranchingObject.cpp

 r539 integer_ (isint) { firstBranch_ = (way == TWO_LEFT)  ? 0 : (way == TWO_RIGHT) ? 1 : (CoinDrand48 () < 0.5) ? 0 : 1; firstBranch_ =  (way == TWO_LEFT)      ? 0 : ((way == TWO_RIGHT)    ? 1 : ((CoinDrand48 () < 0.5) ? 0 : 1)); if (index_ < 0) { alpha = CouenneBranchingObject::Alpha (); if (brpoint > - COIN_DBL_MAX) //  printf ("///////////brpt = %g, x = %Lf, l,u = [%Lf,%Lf]\n", brpoint, x, l, u); if (fabs (brpoint) < COUENNE_INFINITY) x = brpoint; if ((x-l < COUENNE_NEAR_BOUND) || (u-x < COUENNE_NEAR_BOUND)) value_ = alpha * x + (1. - alpha) * (l + u) / 2.; if (u < COUENNE_INFINITY) if (l > - COUENNE_INFINITY)    value_ = alpha * x + (1. - alpha) * (l + u) / 2.; else                           value_ = (u<0) ? -(1+u) : u/2; else if (l > - COUENNE_INFINITY) value_ = (l>0) ?  (1+l) : l/2; else                             value_ = 0; else value_ = x; /*  if ((x > l + COUENNE_NEAR_BOUND) && else                          value_ = u - (1+fabs (u)) / 2.;   */ /*  if (0) { printf ("Branch::constructor: "); reference_ -> print (std::cout); printf (" on %.6e (%.6e) [%.6e,%.6e]\n", if (0) { printf ("=== x%d branches on %g (at %g) [%g,%g]\n", index_, value_, expression::Variable (reference_ -> Index ()), expression::Lbound   (reference_ -> Index ()), expression::Ubound   (reference_ -> Index ())); }*/ expression::Variable (index_), expression::Lbound   (index_), expression::Ubound   (index_)); } } else      solver -> setColLower (index_, integer_ ? ceil  (value_) : value_); // up   branch // TODO: apply bound tightening // TODO: apply bound tightening to evaluate change in dual bound if (0) { //    printf (" --> [%.6e,%.6e]\n", l, u); printf ("### Branch: x%d %c= %.12f\n", printf ("### Branch: x%d %c= %g\n", index_, way ? '>' : '<', value_);
• ## branches/Couenne/Couenne/src/branch/CouenneChooseVariable.cpp

 r515 int CouenneChooseVariable::setupList (OsiBranchingInformation *info, bool initialize) { /// LEAVE THIS HERE. /// LEAVE THIS HERE. Need update of current point to evaluate infeasibility problem_ -> update ((CouNumber *) (info -> solution_), (CouNumber *) (info -> lower_),
• ## branches/Couenne/Couenne/src/branch/CouenneObject.cpp

 r539 //expression::Variable (reference_ -> Index ()); if (0) { /* if (1) { printf ("Inf: = |%.6e - %.6e| = %.6e  ",  ////[%.2f,%.2f] reference_ -> Image () -> print (std::cout); printf ("\n"); } */ CouNumber delta = fabs (var - expr); // that should be used later in createBranch. // TODO: how to avoid this part of the code when called from // CbcModel::setSolution() ? CouNumber improv = reference_ -> Image () -> selectBranch (reference_, info,             // input parameters delta = improv; if (fabs (delta) < COUENNE_EPS) delta = 0.; // This should make getFixVar() useless if used in exprMul or // exprDiv, i.e., the only non-unary operators. // TODO: how to avoid this part of the code when called from // CbcModel::setSolution() ? // make delta a function of the variable's rank and multiplicity + WEI_MULT * (1. - 1. / fixvar -> Multiplicity ());*/ } if (0) { printf ("Inf |%+.4e - %+.4e| = %+.4e  %+.4e ",  ////[%.2f,%.2f] var, expr, //      expression::Lbound (reference_ -> Index ()), //      expression::Ubound (reference_ -> Index ()), fabs (var - expr), delta); reference_             -> print (std::cout); std::cout << " = "; reference_ -> Image () -> print (std::cout); printf ("\n"); } return delta; // fix all variables upon which this auxiliary depends if (expr -> Argument ()){ // unary function if (expr -> Argument ()) { // unary function index = expr -> Argument () -> Index (); bool isint = reference_ -> isInteger (); if (0) if (brVarInd_ >= 0) if (brVarInd_ >= 0) // if applied latest selectBranching switch (way) { case TWO_LEFT: case TWO_RIGHT: case TWO_RAND: //printf ("2 Branch x%d at %g [%d] (%d)\n", brVarInd_, *brPts_, way, isint); return new CouenneBranchingObject (brVarInd_, way, *brPts_, isint); case THREE_LEFT: case THREE_RIGHT: case THREE_RAND: //printf ("3 WAY Branch x%d at %g--%g [%d] (%d)\n", brVarInd_, *brPts_, brPts_ [1], way, isint); return new CouenneThreeWayBranchObj (brVarInd_, brPts_ [0], brPts_ [1], way, isint); default: exit (-1); } // if selectBranch returned -1, apply default branching rule if (0) { // the current point to the two new convexifications, which is the // function that we want to maximize. //  CouNumber delta = reference_ -> Image () -> BranchGain (reference_, info); // operator-dependent branching object //  OsiBranchingObject *opobj = reference_ -> Image () -> BranchObject (reference_, info); //  if (opobj) //    return opobj; expression *depvar = reference_ -> Image () -> getFixVar ();
• ## branches/Couenne/Couenne/src/branch/CouenneThreeWayBranchObj.cpp

 r537 first_ = 0; }*/ if (0) { printf ("=3= x%d branches on %g and %g (at %g) [%g,%g]\n", index_, lcrop_, rcrop_, expression::Variable (index_), expression::Lbound   (index_), expression::Ubound   (index_)); } } // TODO: apply bound tightening /*if (0) { printf (" --> [%.6e,%.6e]\n", l, u); printf ("### Branch: x%d %c= %.12f\n", reference_ -> Index (), way ? '>' : '<', value_); for (int i=0; i < solver -> getNumCols(); i++) printf (" %3d [%.3e,%.3e]\n", i, solver -> getColLower () [i], solver -> getColUpper () [i]); }*/ if (0) { switch (way) { case -1: printf ("#3# Branch: x%d <= %g\n",               index_, lcrop_); break; // left case  0: printf ("#3# Branch: %g <= x%d <= %g\n", lcrop_, index_, rcrop_); break; // center case  1: printf ("#3# Branch: x%d >= %g\n",               index_, rcrop_); break; // right default: printf ("Couenne: branching on nonsense way %d\n", way); } } /*printf (" --> [%.6e,%.6e]\n", l, u); printf ("### Branch: x%d %c= %.12f\n", reference_ -> Index (), way ? '>' : '<', value_); for (int i=0; i < solver -> getNumCols(); i++) printf (" %3d [%.3e,%.3e]\n", i, solver -> getColLower () [i], solver -> getColUpper () [i]);*/ branchIndex_++;
• ## branches/Couenne/Couenne/src/branch/operators/branchExprAbs.cpp

 r537 * Name:    branchExprAbs.cpp * Author:  Pietro Belotti * Purpose: return branch gain and branch object for exprAbs * Purpose: return branch suggestion for exprAbs * * (C) Pietro Belotti. This file is licensed under the Common Public License (CPL) ind = argument_ -> Index (); CouNumber x0 = info -> solution_ [ind], y0 = info -> solution_ [w -> Index ()]; if (ind < 0) { printf ("argument of exprAbs has negative index\n"); } CouNumber x0 = info -> solution_ [ind], y0 = info -> solution_ [w -> Index ()]; brpts = (double *) realloc (brpts, sizeof (double)); // the smartest branching point for |x| is 0, as the two subproblems // will have exact convexifications *brpts = 0.; way = TWO_RAND; // don't really care which subtree to visit first // exact distance from current point to the two subsequent // convexifications return mymin (project (1., -1., 0., x0, y0, 0., COUENNE_INFINITY,  0, NULL, NULL), project (1., +1., 0., x0, y0, -COUENNE_INFINITY, 0., 0, NULL, NULL)); } /* /// distance covered by current point if branching rule applied to this expression double exprAbs::BranchGain (expression *w, const OsiBranchingInformation *info) { int xi = argument_ -> Index (), wi = w         -> Index (); if (xi < 0 || wi < 0) { printf ("negative indices inbranchExprAbs\n"); exit (-1); } CouNumber x0 = expression::Variable (xi), y0 = expression::Variable (wi); return } /// branching object best suited for this expression OsiBranchingObject *exprAbs::BranchObject (expression *w, const OsiBranchingInformation *) { // branching once on 0 is sufficient for expressions w=|x| return new CouenneBranchingObject (Argument (), 0); } */
• ## branches/Couenne/Couenne/src/branch/operators/branchExprDiv.cpp

 r537 * Name:    branchExprDiv.cpp * Author:  Pietro Belotti * Purpose: return branch (expected) gain and branch object for divisions * Purpose: select branch for divisions * * (C) Pietro Belotti. This file is licensed under the Common Public License (CPL) #define BR_NEXT_ZERO 1e-3 #define BR_MULT      1e-3 /// set up branching object by evaluating many branching points for /// each expression's arguments int &way) { ind = -1; return 0.; int xi = arglist_ [0] -> Index (), yi = arglist_ [1] -> Index (), wi = w            -> Index (); //  ind = argument_ -> Index (); /* if (ind < 0) { printf ("argument of exprAbs has negative index\n"); if ((xi < 0) || (yi < 0)) { printf ("arguments of exprDiv have negative index\n"); exit (-1); } brpts = (double *) malloc (sizeof (double)); *brpts = 0.; way = TWO_LEFT; // choosing branching variable and -point is difficult, use // proportion in bound intervals return mymin (project (1., -1., 0., x0, y0, 0., COUENNE_INFINITY,  0, NULL, NULL), project (1., +1., 0., x0, y0, -COUENNE_INFINITY, 0., 0, NULL, NULL));*/ CouNumber yl = info -> lower_    [yi], yu = info -> upper_    [yi], y0 = info -> solution_ [yi]; // if [yl,yu] contains 0, create three nodes, with a small central // one tight around zero -- it has a bad convexification and will be // explored last if ((yl < 0) && (yu > 0)) { ind = yi; brpts = (double *) realloc (brpts, 2 * sizeof (CouNumber)); way = THREE_RIGHT; brpts [0] = (yl >= -BR_NEXT_ZERO - COUENNE_EPS) ? (yl * BR_MULT) : -BR_NEXT_ZERO; brpts [1] = (yu <=  BR_NEXT_ZERO + COUENNE_EPS) ? (yu * BR_MULT) :  BR_NEXT_ZERO; return ((fabs (y0) < COUENNE_EPS) ? 1. : (info -> solution_ [xi] / y0 - info -> solution_ [w -> Index ()])); } // [yl,yu] can only be unlimited in one sense, and interval does not // contain 0. // // As convexification is still carried out by applying McCormick // rules to x=w*y (where original operator is w=x/y), try to get // closer to a point where both y and w are bounded, if necessary by // branching on w. // // First, branch on y if unbounded, and then on w. As a result of // bound tightening, if both y and w are bounded, x is, too. if ((yl < -COUENNE_INFINITY) || (yu >  COUENNE_INFINITY)) { ind = yi; brpts = (double *) realloc (brpts, sizeof (CouNumber)); // if y0 close to bounds, branch away from it if      (fabs (y0-yl) < COUENNE_NEAR_BOUND) *brpts = y0 + 1 + yl*10; else if (fabs (y0-yu) < COUENNE_NEAR_BOUND) *brpts = y0 - 1 + yu*10; else                                        *brpts = y0; way = (y0 > 0) ? TWO_LEFT : TWO_RIGHT; return ((fabs (y0) < COUENNE_EPS) ? 1. : (info -> solution_ [xi] / y0 - info -> solution_ [w -> Index ()])); } // y is bounded, and y0 should not be 0; if w is unbounded, it is // better to bound w first as x will be too. CouNumber wl = info -> lower_    [wi], wu = info -> upper_    [wi], w0 = info -> solution_ [wi], x0 = info -> solution_ [xi]; if ((wl < -COUENNE_INFINITY) || (wu >  COUENNE_INFINITY)) { ind = wi; if ((wl < -COUENNE_INFINITY) && (wu >  COUENNE_INFINITY)) { // unbounded in two directions CouNumber wreal = x0 / y0, wmin = w0, wmax = wreal; // (x0,y0,w0) is below surface if (wreal < w0) { // (x0,y0,w0) is above surface wmin = wreal; wmax = w0; } // unbounded in both directions: three way branching brpts = (double *) realloc (brpts, 2 * sizeof (CouNumber)); brpts [0] = wmin; brpts [1] = wmax; way = THREE_CENTER; } else { // unbounded in one direction only, use two way branching brpts = (double *) realloc (brpts, sizeof (CouNumber)); // if y0 close to bounds, branch away from it if      (fabs (w0-wl) < COUENNE_NEAR_BOUND) *brpts = w0 + 1 + wl*10; else if (fabs (w0-wu) < COUENNE_NEAR_BOUND) *brpts = w0 - 1 + wu*10; else                                        *brpts = w0; way = (wl < - COUENNE_INFINITY) ? TWO_RIGHT : TWO_LEFT; } return ((fabs (y0) < COUENNE_EPS) ? 1. : (x0/y0 - w0)); } // w and y are bounded (and so is x). Choose between x, y, z // depending on intervals first and then to vicinity to midpoint CouNumber xl = info -> lower_ [xi], xu = info -> upper_ [xi]; CouNumber dx = xu-xl, dy = yu-yl, dw = wu-wl; brpts = (double *) realloc (brpts, sizeof (CouNumber)); // Check largest interval and choose branch variable accordingly. // Branching point depends on where the current point is, but for // now just focus on the width of the intervals way = TWO_RAND; if (dx > dy) if (dx > dw) {ind = xi; *brpts = (xl + xu) / 2.; return fabs (x0 - y0*w0);} // dx maximum else         {ind = wi; *brpts = (wl + wu) / 2.; return fabs (w0 - x0/y0);} // dw else if (dy > dw) {ind = yi; *brpts = (yl + yu) / 2.; return fabs (y0 - x0/w0);} // dy else         {ind = wi; *brpts = (wl + wu) / 2.; return fabs (w0 - x0/y0);} // dw } /* /// distance covered by current point if branching rule applied to this expression double exprDiv::BranchGain (expression *w, const OsiBranchingInformation *info) { return 0; } /// branching object best suited for this expression OsiBranchingObject *exprDiv::BranchObject (expression *w, const OsiBranchingInformation *) { return NULL; } */
• ## branches/Couenne/Couenne/src/convex/CouenneCutGenerator.h

 r555 /// should we add the violated cuts only (true), or all of them (false)? bool addviolated_; mutable bool addviolated_; /// what kind of sampling should be performed? Bonmin::Bab *BabPtr_; /// signal infeasibility of current node (found through /// signal infeasibility of current node (found through bound tightening) mutable bool infeasNode_;
• ## branches/Couenne/Couenne/src/convex/genRowCuts.cpp

 r545 /*expression * image = problem_ -> Aux (i) -> Image (); if (   (image -> Linearity () > LINEAR)          // if linear, no need to cut twice && (image -> dependsOn (changed, nchanged))  // if expression does not depend on )*/                                            // changed variables, do not cut && (image -> dependsOn (changed, nchanged)) // if expression does not depend // on changed variables, do not cut || !changed)*/                              // or if no changed variable is passed if (problem_ -> Aux (i) -> Multiplicity () > 0) problem_ -> Aux (i) -> generateCuts (si, cs, this);
• ## branches/Couenne/Couenne/src/convex/generateCuts.cpp

 r555 // maximum number of obbt iterations #define MAX_OBBT_ITER 4 #define MAX_OBBT_ITER 1 #define LARGE_TOL (LARGE_BOUND/1e6) // // set and lift bound for auxiliary variable associated with objective // function void fictitiousBound (OsiCuts &cs, CouenneProblem *p, else if (sense>0) {if (fabs (p->Ub(ind_obj)-LARGE_BOUND)Ub(ind_obj) = COUENNE_INFINITY;} else          if (fabs (p->Lb(ind_obj)+LARGE_BOUND)Lb(ind_obj) = -COUENNE_INFINITY; else          if (fabs (p->Lb(ind_obj)+LARGE_BOUND)Lb(ind_obj) =-COUENNE_INFINITY; } const CglTreeInfo info) const { //  if (firstcall_) //    fictitiousBound (si, problem_, true); int nInitCuts = cs.sizeRowCuts (); //  printf (":::::::::::::::::::::::: level = %d, pass = %d, intree=%d\n Bounds:\n", //  printf (":::::::::: level = %d, pass = %d, intree=%d\n",// Bounds:\n", //  info.level, info.pass, info.inTree); Bonmin::BabInfo * babInfo = dynamic_cast (si.getAuxiliaryInfo ()); if (babInfo){ babInfo -> setFeasibleNode (); // PIERRE -> PIETRO //HERE IS THE CODE TO OBTAIN SOLUTION JUST FOUND BY NLP SOLVER (AUXILIARIES SHOULD BE //CORRECT. SOLUTION SHOULD BE THE ONE FOUND AT THE NODE EVEN IF IT IS NOT AS GOOD AS THE // BEST KNOWN. // YOU CAN DELETE THE COMMENT std::cout<<"Nlp solved: "<< (babInfo -> nlpSolution() != NULL) < nlpSolution(); babInfo -> setHasNlpSolution(false);//Have to reset it after use otherwise will stay true at next processed node. } if (babInfo) babInfo -> setFeasibleNode (); double now   = CoinCpuTime (); int    ncols = problem_ -> nVars () + problem_ -> nAuxs (); sparse2dense (ncols, chg_bds, changed, nchanged); ////////////////////// genRowCuts (si, cs, nchanged, changed); ////////////////////// /* if (firstcall_) { // set trivial dual bound to objective function, if there is none int ind_obj = problem_ -> Obj (0) -> Body () -> Index (); if (ind_obj >= 0) { if (problem_ -> Obj (0) -> Sense () == MINIMIZE) { if (problem_ -> Lb (ind_obj) < - LARGE_BOUND) problem_ -> Lb (ind_obj) = - LARGE_BOUND; } else if (problem_ -> Ub (ind_obj) > LARGE_BOUND) problem_ -> Ub (ind_obj) = LARGE_BOUND; } } */ double *nlpSol; //-------------------------------------------- if (babInfo && ((nlpSol = const_cast (babInfo -> nlpSolution ())))) { // obtain solution just found by nlp solver // Auxiliaries should be correct. solution should be the one found // at the node even if it is not as good as the best known. // save violation flag and disregard it while adding cut at NLP // point (which are not violated by the current, NLP, solution) bool save_av = addviolated_; addviolated_ = false; // update problem current point with NLP solution problem_ -> update (nlpSol, NULL, NULL); genRowCuts (si, cs, nchanged, changed);  // add cuts addviolated_ = save_av;     // restore previous value babInfo -> setHasNlpSolution (false); // reset it after use } else genRowCuts (si, cs, nchanged, changed); //--------------------------------------------- // change tightened bounds through OsiCuts int nCurCuts = cs.sizeRowCuts (); genRowCuts (si, cs, nchanged, changed); genRowCuts (si, cs, nchanged, changed);// !!! repeat = nCurCuts < cs.sizeRowCuts (); // reapply only if new cuts available }
• ## branches/Couenne/Couenne/src/convex/obbt.cpp

 r553 // create clone of current solver interface // TODO: clone only at pass 0 and delete at the end OsiSolverInterface *csi = si.clone (true); // for all (original+auxiliary) variables x_i, for (int i=0; inVars(); i++) if (i != objind) { // do not improve objective's bounds
• ## branches/Couenne/Couenne/src/convex/operators/conv-exprInv.cpp

 r551 aux -> getBounds (wle, wue); CouNumber x; //  CouNumber x; int w_ind = aux       -> Index (),