Changeset 556


Ignore:
Timestamp:
May 28, 2007 6:59:50 PM (12 years ago)
Author:
pbelotti
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

Legend:

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

    r539 r556  
    2828  integer_ (isint) {
    2929
    30   firstBranch_ = (way == TWO_LEFT)  ? 0 :
    31                  (way == TWO_RIGHT) ? 1 :
    32                  (CoinDrand48 () < 0.5) ? 0 : 1;
     30  firstBranch_ =  (way == TWO_LEFT)      ? 0 :
     31                 ((way == TWO_RIGHT)    ? 1 :
     32                 ((CoinDrand48 () < 0.5) ? 0 : 1));
    3333
    3434  if (index_ < 0) {
     
    4343    alpha = CouenneBranchingObject::Alpha ();
    4444
    45   if (brpoint > - COIN_DBL_MAX)
     45  //  printf ("///////////brpt = %g, x = %Lf, l,u = [%Lf,%Lf]\n", brpoint, x, l, u);
     46
     47  if (fabs (brpoint) < COUENNE_INFINITY)
    4648    x = brpoint;
    4749
     
    7072  if ((x-l < COUENNE_NEAR_BOUND) ||
    7173      (u-x < COUENNE_NEAR_BOUND))
    72     value_ = alpha * x + (1. - alpha) * (l + u) / 2.;
     74    if (u < COUENNE_INFINITY)
     75      if (l > - COUENNE_INFINITY)    value_ = alpha * x + (1. - alpha) * (l + u) / 2.;
     76      else                           value_ = (u<0) ? -(1+u) : u/2;
     77    else if (l > - COUENNE_INFINITY) value_ = (l>0) ?  (1+l) : l/2;
     78    else                             value_ = 0;
    7379  else value_ = x;
     80
    7481
    7582  /*  if ((x > l + COUENNE_NEAR_BOUND) &&
     
    8996      else                          value_ = u - (1+fabs (u)) / 2.;   */
    9097
    91   /*  if (0) {
    92     printf ("Branch::constructor: ");
    93     reference_ -> print (std::cout);
    94     printf (" on %.6e (%.6e) [%.6e,%.6e]\n",
     98  if (0) {
     99    printf ("=== x%d branches on %g (at %g) [%g,%g]\n",
     100            index_,
    95101            value_,
    96             expression::Variable (reference_ -> Index ()),
    97             expression::Lbound   (reference_ -> Index ()),
    98             expression::Ubound   (reference_ -> Index ()));
    99             }*/
     102            expression::Variable (index_),
     103            expression::Lbound   (index_),
     104            expression::Ubound   (index_));
     105  }
    100106}
    101107
     
    132138  else      solver -> setColLower (index_, integer_ ? ceil  (value_) : value_); // up   branch
    133139
    134   // TODO: apply bound tightening
     140  // TODO: apply bound tightening to evaluate change in dual bound
    135141
    136142  if (0) {
    137143
    138144    //    printf (" --> [%.6e,%.6e]\n", l, u);
    139     printf ("### Branch: x%d %c= %.12f\n",
     145    printf ("### Branch: x%d %c= %g\n",
    140146            index_, way ? '>' : '<', value_);
    141147           
  • branches/Couenne/Couenne/src/branch/CouenneChooseVariable.cpp

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

    r539 r556  
    4747
    4848  //expression::Variable (reference_ -> Index ());
    49  if (0) {
     49  /* if (1) {
    5050
    5151    printf ("Inf: = |%.6e - %.6e| = %.6e  ",  ////[%.2f,%.2f]
     
    5757    reference_ -> Image () -> print (std::cout); printf ("\n");
    5858  }
     59  */
    5960
    6061  CouNumber delta = fabs (var - expr);
     
    102103    // that should be used later in createBranch.
    103104
     105    // TODO: how to avoid this part of the code when called from
     106    // CbcModel::setSolution() ?
     107
    104108    CouNumber improv = reference_ -> Image () -> selectBranch
    105109      (reference_, info,             // input parameters
     
    109113      delta = improv;
    110114
     115    if (fabs (delta) < COUENNE_EPS)
     116      delta = 0.;
     117
    111118    // This should make getFixVar() useless if used in exprMul or
    112119    // exprDiv, i.e., the only non-unary operators.
    113 
    114     // TODO: how to avoid this part of the code when called from
    115     // CbcModel::setSolution() ?
    116120
    117121    // make delta a function of the variable's rank and multiplicity
     
    121125            + WEI_MULT * (1. - 1. / fixvar -> Multiplicity ());*/
    122126  }
     127
     128  if (0) {
     129
     130    printf ("Inf |%+.4e - %+.4e| = %+.4e  %+.4e ",  ////[%.2f,%.2f]
     131            var, expr,
     132            //      expression::Lbound (reference_ -> Index ()),
     133            //      expression::Ubound (reference_ -> Index ()),
     134            fabs (var - expr), delta);
     135
     136    reference_             -> print (std::cout); std::cout << " = ";
     137    reference_ -> Image () -> print (std::cout); printf ("\n");
     138  }
    123139
    124140  return delta;
     
    146162  // fix all variables upon which this auxiliary depends
    147163
    148   if (expr -> Argument ()){ // unary function
     164  if (expr -> Argument ()) { // unary function
    149165
    150166    index = expr -> Argument () -> Index ();
     
    198214  bool isint = reference_ -> isInteger ();
    199215
    200   if (0)
    201   if (brVarInd_ >= 0)
     216  if (brVarInd_ >= 0) // if applied latest selectBranching
    202217    switch (way) {
    203218    case TWO_LEFT:
    204219    case TWO_RIGHT:
    205220    case TWO_RAND:
     221      //printf ("2 Branch x%d at %g [%d] (%d)\n", brVarInd_, *brPts_, way, isint);
    206222      return new CouenneBranchingObject (brVarInd_, way, *brPts_, isint);
    207223    case THREE_LEFT:
     
    209225    case THREE_RIGHT:
    210226    case THREE_RAND:
     227      //printf ("3 WAY Branch x%d at %g--%g [%d] (%d)\n", brVarInd_, *brPts_, brPts_ [1], way, isint);
    211228      return new CouenneThreeWayBranchObj (brVarInd_, brPts_ [0], brPts_ [1], way, isint);
    212229    default:
     
    214231      exit (-1);
    215232    }
     233
     234  // if selectBranch returned -1, apply default branching rule
    216235
    217236  if (0) {
     
    238257  // the current point to the two new convexifications, which is the
    239258  // function that we want to maximize.
    240 
    241   //  CouNumber delta = reference_ -> Image () -> BranchGain (reference_, info);
    242 
    243   // operator-dependent branching object
    244   //  OsiBranchingObject *opobj = reference_ -> Image () -> BranchObject (reference_, info);
    245 
    246   //  if (opobj)
    247   //    return opobj;
    248259
    249260  expression *depvar = reference_ -> Image () -> getFixVar ();
  • branches/Couenne/Couenne/src/branch/CouenneThreeWayBranchObj.cpp

    r537 r556  
    6464        first_ = 0;
    6565        }*/
     66
     67  if (0) {
     68    printf ("=3= x%d branches on %g and %g (at %g) [%g,%g]\n",
     69            index_,
     70            lcrop_,
     71            rcrop_,
     72            expression::Variable (index_),
     73            expression::Lbound   (index_),
     74            expression::Ubound   (index_));
     75  }
    6676}
    6777
     
    102112  // TODO: apply bound tightening
    103113
    104   /*if (0) {
    105   printf (" --> [%.6e,%.6e]\n", l, u);
    106   printf ("### Branch: x%d %c= %.12f\n",
    107   reference_ -> Index (), way ? '>' : '<', value_);
    108   for (int i=0; i < solver -> getNumCols(); i++)
    109   printf (" %3d [%.3e,%.3e]\n", i, solver -> getColLower () [i],
    110   solver -> getColUpper () [i]);
    111   }*/
     114  if (0) {
     115    switch (way) {
     116    case -1: printf ("#3# Branch: x%d <= %g\n",               index_, lcrop_); break; // left
     117    case  0: printf ("#3# Branch: %g <= x%d <= %g\n", lcrop_, index_, rcrop_); break; // center
     118    case  1: printf ("#3# Branch: x%d >= %g\n",               index_, rcrop_); break; // right
     119    default: printf ("Couenne: branching on nonsense way %d\n", way);
     120    }
     121  }
     122    /*printf (" --> [%.6e,%.6e]\n", l, u);
     123    printf ("### Branch: x%d %c= %.12f\n",
     124            reference_ -> Index (), way ? '>' : '<', value_);
     125    for (int i=0; i < solver -> getNumCols(); i++)
     126      printf (" %3d [%.3e,%.3e]\n", i, solver -> getColLower () [i],
     127      solver -> getColUpper () [i]);*/
    112128
    113129  branchIndex_++;
  • branches/Couenne/Couenne/src/branch/operators/branchExprAbs.cpp

    r537 r556  
    22 * Name:    branchExprAbs.cpp
    33 * Author:  Pietro Belotti
    4  * Purpose: return branch gain and branch object for exprAbs
     4 * Purpose: return branch suggestion for exprAbs
    55 *
    66 * (C) Pietro Belotti. This file is licensed under the Common Public License (CPL)
     
    2424  ind = argument_ -> Index ();
    2525
    26   CouNumber x0 = info -> solution_ [ind],
    27             y0 = info -> solution_ [w -> Index ()];
    28 
    2926  if (ind < 0) {
    3027    printf ("argument of exprAbs has negative index\n");
     
    3229  }
    3330
     31  CouNumber x0 = info -> solution_ [ind],
     32            y0 = info -> solution_ [w -> Index ()];
     33
    3434  brpts = (double *) realloc (brpts, sizeof (double));
     35
     36  // the smartest branching point for |x| is 0, as the two subproblems
     37  // will have exact convexifications
    3538  *brpts = 0.;
     39
    3640  way = TWO_RAND; // don't really care which subtree to visit first
    3741
     42  // exact distance from current point to the two subsequent
     43  // convexifications
    3844  return mymin (project (1., -1., 0., x0, y0, 0., COUENNE_INFINITY,  0, NULL, NULL),
    3945                project (1., +1., 0., x0, y0, -COUENNE_INFINITY, 0., 0, NULL, NULL));
    4046}
    41 
    42 /*
    43 /// distance covered by current point if branching rule applied to this expression
    44 double exprAbs::BranchGain (expression *w, const OsiBranchingInformation *info) {
    45 
    46   int xi = argument_ -> Index (),
    47       wi = w         -> Index ();
    48 
    49   if (xi < 0 || wi < 0) {
    50     printf ("negative indices inbranchExprAbs\n");
    51     exit (-1);
    52   }
    53 
    54   CouNumber x0 = expression::Variable (xi),
    55             y0 = expression::Variable (wi);
    56 
    57   return
    58 }
    59 
    60 /// branching object best suited for this expression
    61 OsiBranchingObject *exprAbs::BranchObject (expression *w, const OsiBranchingInformation *) {
    62 
    63   // branching once on 0 is sufficient for expressions w=|x|
    64   return new CouenneBranchingObject (Argument (), 0);
    65 }
    66 */
  • branches/Couenne/Couenne/src/branch/operators/branchExprDiv.cpp

    r537 r556  
    22 * Name:    branchExprDiv.cpp
    33 * Author:  Pietro Belotti
    4  * Purpose: return branch (expected) gain and branch object for divisions
     4 * Purpose: select branch for divisions
    55 *
    66 * (C) Pietro Belotti. This file is licensed under the Common Public License (CPL)
     
    1515
    1616
     17#define BR_NEXT_ZERO 1e-3
     18#define BR_MULT      1e-3
     19
    1720/// set up branching object by evaluating many branching points for
    1821/// each expression's arguments
     
    2326                                 int &way) {
    2427
    25   ind = -1;
    26   return 0.;
     28  int xi = arglist_ [0] -> Index (),
     29      yi = arglist_ [1] -> Index (),
     30      wi = w            -> Index ();
    2731
    28   //  ind = argument_ -> Index ();
    29   /*
    30   if (ind < 0) {
    31     printf ("argument of exprAbs has negative index\n");
     32  if ((xi < 0) || (yi < 0)) {
     33    printf ("arguments of exprDiv have negative index\n");
    3234    exit (-1);
    3335  }
    3436
    35   brpts = (double *) malloc (sizeof (double));
    36   *brpts = 0.;
    37   way = TWO_LEFT;
     37  // choosing branching variable and -point is difficult, use
     38  // proportion in bound intervals
    3839
    39   return mymin (project (1., -1., 0., x0, y0, 0., COUENNE_INFINITY,  0, NULL, NULL),
    40   project (1., +1., 0., x0, y0, -COUENNE_INFINITY, 0., 0, NULL, NULL));*/
     40  CouNumber yl = info -> lower_    [yi],
     41            yu = info -> upper_    [yi],
     42            y0 = info -> solution_ [yi];
    4143
     44  // if [yl,yu] contains 0, create three nodes, with a small central
     45  // one tight around zero -- it has a bad convexification and will be
     46  // explored last
     47
     48  if ((yl < 0) && (yu > 0)) {
     49
     50    ind = yi;
     51    brpts = (double *) realloc (brpts, 2 * sizeof (CouNumber));
     52    way = THREE_RIGHT;
     53
     54    brpts [0] = (yl >= -BR_NEXT_ZERO - COUENNE_EPS) ? (yl * BR_MULT) : -BR_NEXT_ZERO;
     55    brpts [1] = (yu <=  BR_NEXT_ZERO + COUENNE_EPS) ? (yu * BR_MULT) :  BR_NEXT_ZERO;
     56
     57    return ((fabs (y0) < COUENNE_EPS) ? 1. :
     58            (info -> solution_ [xi] / y0 - info -> solution_ [w -> Index ()]));
     59  }
     60
     61  // [yl,yu] can only be unlimited in one sense, and interval does not
     62  // contain 0.
     63  //
     64  // As convexification is still carried out by applying McCormick
     65  // rules to x=w*y (where original operator is w=x/y), try to get
     66  // closer to a point where both y and w are bounded, if necessary by
     67  // branching on w.
     68  //
     69  // First, branch on y if unbounded, and then on w. As a result of
     70  // bound tightening, if both y and w are bounded, x is, too.
     71
     72  if ((yl < -COUENNE_INFINITY) ||
     73      (yu >  COUENNE_INFINITY)) {
     74
     75    ind = yi;
     76    brpts = (double *) realloc (brpts, sizeof (CouNumber));
     77
     78    // if y0 close to bounds, branch away from it
     79    if      (fabs (y0-yl) < COUENNE_NEAR_BOUND) *brpts = y0 + 1 + yl*10;
     80    else if (fabs (y0-yu) < COUENNE_NEAR_BOUND) *brpts = y0 - 1 + yu*10;
     81    else                                        *brpts = y0;
     82
     83    way = (y0 > 0) ? TWO_LEFT : TWO_RIGHT;
     84
     85    return ((fabs (y0) < COUENNE_EPS) ? 1. :
     86            (info -> solution_ [xi] / y0 - info -> solution_ [w -> Index ()]));
     87  }
     88
     89  // y is bounded, and y0 should not be 0; if w is unbounded, it is
     90  // better to bound w first as x will be too.
     91
     92  CouNumber wl = info -> lower_    [wi],
     93            wu = info -> upper_    [wi],
     94            w0 = info -> solution_ [wi],
     95            x0 = info -> solution_ [xi];
     96
     97  if ((wl < -COUENNE_INFINITY) ||
     98      (wu >  COUENNE_INFINITY)) {
     99
     100    ind = wi;
     101
     102    if ((wl < -COUENNE_INFINITY) &&
     103        (wu >  COUENNE_INFINITY)) {
     104
     105      // unbounded in two directions
     106
     107      CouNumber wreal = x0 / y0,
     108        wmin = w0, wmax = wreal; // (x0,y0,w0) is below surface
     109
     110      if (wreal < w0) { // (x0,y0,w0) is above surface
     111        wmin = wreal;
     112        wmax = w0;
     113      }
     114
     115      // unbounded in both directions: three way branching
     116      brpts = (double *) realloc (brpts, 2 * sizeof (CouNumber));
     117      brpts [0] = wmin;
     118      brpts [1] = wmax;
     119      way = THREE_CENTER;
     120
     121    } else {
     122
     123      // unbounded in one direction only, use two way branching
     124
     125      brpts = (double *) realloc (brpts, sizeof (CouNumber));
     126
     127      // if y0 close to bounds, branch away from it
     128      if      (fabs (w0-wl) < COUENNE_NEAR_BOUND) *brpts = w0 + 1 + wl*10;
     129      else if (fabs (w0-wu) < COUENNE_NEAR_BOUND) *brpts = w0 - 1 + wu*10;
     130      else                                        *brpts = w0;
     131
     132      way = (wl < - COUENNE_INFINITY) ? TWO_RIGHT : TWO_LEFT;
     133    }
     134
     135    return ((fabs (y0) < COUENNE_EPS) ? 1. : (x0/y0 - w0));
     136  }
     137
     138  // w and y are bounded (and so is x). Choose between x, y, z
     139  // depending on intervals first and then to vicinity to midpoint
     140  CouNumber xl = info -> lower_ [xi],
     141            xu = info -> upper_ [xi];
     142
     143  CouNumber dx = xu-xl,
     144            dy = yu-yl,
     145            dw = wu-wl;
     146
     147  brpts = (double *) realloc (brpts, sizeof (CouNumber));
     148
     149  // Check largest interval and choose branch variable accordingly.
     150  // Branching point depends on where the current point is, but for
     151  // now just focus on the width of the intervals
     152
     153  way = TWO_RAND;
     154
     155  if (dx > dy)
     156    if (dx > dw) {ind = xi; *brpts = (xl + xu) / 2.; return fabs (x0 - y0*w0);} // dx maximum
     157    else         {ind = wi; *brpts = (wl + wu) / 2.; return fabs (w0 - x0/y0);} // dw
     158  else
     159    if (dy > dw) {ind = yi; *brpts = (yl + yu) / 2.; return fabs (y0 - x0/w0);} // dy
     160    else         {ind = wi; *brpts = (wl + wu) / 2.; return fabs (w0 - x0/y0);} // dw
    42161}
    43 
    44 /*
    45 /// distance covered by current point if branching rule applied to this expression
    46 double exprDiv::BranchGain (expression *w, const OsiBranchingInformation *info) {
    47 
    48   return 0;
    49 }
    50 
    51 /// branching object best suited for this expression
    52 OsiBranchingObject *exprDiv::BranchObject (expression *w, const OsiBranchingInformation *) {
    53 
    54   return NULL;
    55 }
    56 */
  • branches/Couenne/Couenne/src/convex/CouenneCutGenerator.h

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

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

    r555 r556  
    2424
    2525// maximum number of obbt iterations
    26 #define MAX_OBBT_ITER 4
     26#define MAX_OBBT_ITER 1
    2727
    2828#define LARGE_TOL (LARGE_BOUND/1e6)
    2929
    30 //
     30// set and lift bound for auxiliary variable associated with objective
     31// function
    3132void fictitiousBound (OsiCuts &cs,
    3233                      CouenneProblem *p,
     
    4647  else
    4748    if (sense>0) {if (fabs (p->Ub(ind_obj)-LARGE_BOUND)<LARGE_TOL) p->Ub(ind_obj) = COUENNE_INFINITY;}
    48     else          if (fabs (p->Lb(ind_obj)+LARGE_BOUND)<LARGE_TOL) p->Lb(ind_obj) = -COUENNE_INFINITY;
     49    else          if (fabs (p->Lb(ind_obj)+LARGE_BOUND)<LARGE_TOL) p->Lb(ind_obj) =-COUENNE_INFINITY;
    4950}
    5051
     
    7374                                        const CglTreeInfo info) const {
    7475
    75   //  if (firstcall_)
    76   //    fictitiousBound (si, problem_, true);
    77 
    7876  int nInitCuts = cs.sizeRowCuts ();
    7977
    80   //  printf (":::::::::::::::::::::::: level = %d, pass = %d, intree=%d\n Bounds:\n",
     78  //  printf (":::::::::: level = %d, pass = %d, intree=%d\n",// Bounds:\n",
    8179  //  info.level, info.pass, info.inTree);
    8280
    8381  Bonmin::BabInfo * babInfo = dynamic_cast <Bonmin::BabInfo *> (si.getAuxiliaryInfo ());
    8482
    85   if (babInfo){
    86     babInfo -> setFeasibleNode ();
    87    
    88     // PIERRE -> PIETRO
    89     //HERE IS THE CODE TO OBTAIN SOLUTION JUST FOUND BY NLP SOLVER (AUXILIARIES SHOULD BE
    90     //CORRECT. SOLUTION SHOULD BE THE ONE FOUND AT THE NODE EVEN IF IT IS NOT AS GOOD AS THE
    91     // BEST KNOWN.
    92     // YOU CAN DELETE THE COMMENT
    93     std::cout<<"Nlp solved: "<< (babInfo -> nlpSolution() != NULL) <<std::endl;
    94     const double * nlpSol = babInfo -> nlpSolution();
    95     babInfo -> setHasNlpSolution(false);//Have to reset it after use otherwise will stay true at next processed node.
    96   }
     83  if (babInfo)
     84    babInfo -> setFeasibleNode ();   
     85
    9786  double now   = CoinCpuTime ();
    9887  int    ncols = problem_ -> nVars () + problem_ -> nAuxs ();
     
    202191  sparse2dense (ncols, chg_bds, changed, nchanged);
    203192
    204   //////////////////////
    205   genRowCuts (si, cs, nchanged, changed);
    206   //////////////////////
    207 
    208   /*
    209   if (firstcall_) {
    210 
    211     // set trivial dual bound to objective function, if there is none
    212 
    213     int ind_obj = problem_ -> Obj (0) -> Body () -> Index ();
    214 
    215     if (ind_obj >= 0) {
    216       if (problem_ -> Obj (0) -> Sense () == MINIMIZE) {
    217        if (problem_ -> Lb (ind_obj) < - LARGE_BOUND)
    218          problem_ -> Lb (ind_obj) = - LARGE_BOUND;
    219       }
    220       else
    221        if (problem_ -> Ub (ind_obj) > LARGE_BOUND)
    222          problem_ -> Ub (ind_obj) = LARGE_BOUND;
    223     }
    224   }
    225   */
     193  double *nlpSol;
     194
     195  //--------------------------------------------
     196
     197  if (babInfo && ((nlpSol = const_cast <double *> (babInfo -> nlpSolution ())))) {
     198
     199    // obtain solution just found by nlp solver
     200
     201    // Auxiliaries should be correct. solution should be the one found
     202    // at the node even if it is not as good as the best known.
     203
     204    // save violation flag and disregard it while adding cut at NLP
     205    // point (which are not violated by the current, NLP, solution)
     206    bool save_av = addviolated_;
     207    addviolated_ = false;
     208
     209    // update problem current point with NLP solution
     210    problem_ -> update (nlpSol, NULL, NULL);
     211    genRowCuts (si, cs, nchanged, changed);  // add cuts
     212
     213    addviolated_ = save_av;     // restore previous value
     214
     215    babInfo -> setHasNlpSolution (false); // reset it after use
     216  }
     217  else genRowCuts (si, cs, nchanged, changed);
     218
     219  //---------------------------------------------
    226220
    227221  // change tightened bounds through OsiCuts
     
    253247
    254248        int nCurCuts = cs.sizeRowCuts ();
    255         genRowCuts (si, cs, nchanged, changed);
     249        genRowCuts (si, cs, nchanged, changed);// !!!
    256250        repeat = nCurCuts < cs.sizeRowCuts (); // reapply only if new cuts available
    257251      }
  • branches/Couenne/Couenne/src/convex/obbt.cpp

    r553 r556  
    5656
    5757  // create clone of current solver interface
     58  // TODO: clone only at pass 0 and delete at the end
    5859  OsiSolverInterface *csi = si.clone (true);
    5960
     
    6263
    6364  // for all (original+auxiliary) variables x_i,
    64   for (int i=0; i<ncols; i++)
     65  for (int i=0; i<ncols; i++)
     66    //  for (int i=0; i<problem_ ->nVars(); i++)
    6567
    6668    if (i != objind) { // do not improve objective's bounds
  • branches/Couenne/Couenne/src/convex/operators/conv-exprInv.cpp

    r551 r556  
    7373  aux -> getBounds (wle, wue);
    7474
    75   CouNumber x;
     75  //  CouNumber x;
    7676
    7777  int w_ind = aux       -> Index (),
  • branches/Couenne/Couenne/src/convex/operators/conv-exprMul-genCuts.cpp

    r537 r556  
    99#include <CouenneTypes.h>
    1010#include <exprMul.h>
     11#include <exprPow.h>
    1112#include <exprDiv.h>
    1213#include <CouenneProblem.h>
     
    1617void addImplTangent (const CouenneCutGenerator *, OsiCuts &,
    1718                     CouNumber, CouNumber, CouNumber, int, int, int, int);
     19
     20/// used below to specify k in curve y=k/x
     21static CouNumber multInv;
     22
     23inline CouNumber kinv   (register CouNumber x)
     24{return multInv / x;}
     25
     26inline CouNumber kinvp  (register CouNumber x)
     27{return -multInv / (x*x);}
     28
     29inline CouNumber kinvpp (register CouNumber x)
     30{return 2*multInv / (x*x*x);}
     31
     32
     33/// Add cut around curve x*y=k
     34
     35void contourCut (const CouenneCutGenerator *cg,
     36                 OsiCuts &cs,
     37                 CouNumber xp, CouNumber yp, // current point
     38                 CouNumber wb,               // bound on w
     39                 int sign,                   // is wb lower or upper?
     40                 CouNumber x0, CouNumber y0, // (allegedly) outside point
     41                 CouNumber x1, CouNumber y1, //             inside
     42                 int xi, int yi, int wi) {   // indices of the variables
     43
     44  // Upper right corner of the bounding box of (x,y) is feasible,
     45  // the opposite corner is not, hence there is a cut violated by
     46  // (x0,y0).
     47
     48  // If (x0,y0) is not in the same orthant as the contour in
     49  // question, move it in so that we can apply a Newton step to
     50  // find closest point on contour.
     51
     52  int xsign = (x1 >= 0) ? 1 : -1, // define orthant where the "inside
     53      ysign = (y1 >= 0) ? 1 : -1; // point" lies
     54
     55  if      (((xsign > 0) ? xp : -xp) <= COUENNE_EPS)
     56    if    (((ysign > 0) ? yp : -yp) <= COUENNE_EPS) {
     57
     58      // opposite orthant, put in the right one where constraint is violated
     59      xp = yp = sqrt (fabs (wb))/2;
     60      if (xsign<0) xp = -xp;
     61      if (ysign<0) yp = -yp;
     62    }                                                // otherwise, must cross one axis only:
     63    else                                            {xp = sqrt (fabs(wb/yp)); if (xsign<0) xp=-xp;}//y
     64  else if (((ysign > 0) ? yp : -yp) <= COUENNE_EPS) {yp = sqrt (fabs(wb/xp)); if (ysign<0) yp=-yp;}//x
     65
     66  multInv = wb;
     67
     68  CouNumber
     69    // tangent point closest to current point
     70    xt    = powNewton (xp, yp, kinv, kinvp, kinvpp),
     71    // coefficient of w in the lifted cut
     72    alpha = ((fabs (x1) < COUENNE_INFINITY) &&
     73             (fabs (y1) < COUENNE_INFINITY)) ?
     74       ((2*wb/xt - y1 - wb*x1 / (xt*xt)) / (x1*y1 - wb)) : 0;
     75
     76  //  printf ("+++++ %d %d %d. [%c] xp (%g,%g) wb %g out(%g,%g) in(%g,%g) --> [%g,%g] alpha %g\n",
     77  //             xi, yi, wi, (sign<0) ? '-' : '+', xp, yp, wb, x0, y0, x1, y1, xt, wb/xt, alpha);
     78
     79  if (alpha != 0)
     80    cg     -> createCut (cs, alpha*wb + 2*wb/xt, sign, wi, alpha, yi, 1., xi, wb/(xt*xt));
     81  else  cg -> createCut (cs,            2*wb/xt, sign,            yi, 1., xi, wb/(xt*xt));
     82}
     83
    1884
    1985/// generate convexification cut for constraint w = x*y
     
    140206  if (is_boundbox_regular (yu, xl)) cg -> createCut (cs, yu*xl, +1, wi, -1., xi, yu, yi, xl);
    141207
    142   // extra cuts for the two cases w = xy >= l > 0 and w = xy <= u < 0
    143   //
    144   // These cuts are added at the tangent of the curve xy=k, on the
    145   // intersection with the bounding box. These are in fact NOT implied
    146   // by the above cuts (as happens for division, for instance) and may
    147   // be of help.
    148 
    149   // TODO: are these really useful/correct?
    150 
    151   if (wu < - COUENNE_EPS) {
    152     // check points A and B: second orthant intersections
    153     if ((xu*yl > wu) && (xl*yu <= wu)) {
    154       addImplTangent (cg, cs, xl, yl, wu, xi, yi, +1, +1); // A
    155       addImplTangent (cg, cs, xu, yu, wu, xi, yi, -1, +1); // B
     208  // add different cuts, to cut out current point in bounding box but
     209  // out of the hyperbola's belly
     210
     211  CouNumber x0 = (*(arglist_ [0])) (),
     212            y0 = (*(arglist_ [1])) ();
     213
     214  // If w=xy and w >= l > 0 (resp. w <= u < 0) are "tight" bounds
     215  // (i.e. they are tighter than those obtained through propagation of
     216  // x and y's bounds), McCormick's convexification is not tight as
     217  // the surface has a curve contour at w=l (resp. w=u).
     218  //
     219  // McCormick rules induce a tangent to this contour at the bounds of
     220  // both variables, but it may be useful to add further cuts along
     221  // the contour to eliminate infeasible point (x0,y0), whose product
     222  // is in the convexification but out of the contour (on its "convex"
     223  // side, or "out of the belly").
     224  //
     225  // Suppose P (xt,l/xt) (resp. (xt,u/xt) is the point on the contour
     226  // closest to (x0,y0), found through a Newton method. The cut is
     227  // tangent to the contour in P and has the form
     228  //
     229  //        y - l/xt >= -l/(xt^2) (x-xt)   if xl*yl < l and xu*yu > l
     230  //        y - l/xt <= -l/(xt^2) (x-xt)   if xl*yl > l and xu*yu < l
     231  //
     232  // (resp. y - u/xt <= -u/(xt^2) (x-xt)   if xl*yu > u and xu*yl < u
     233  //        y - u/xt >= -u/(xt^2) (x-xt)   if xl*yu < u and xu*yl > u)
     234  //
     235  // These can be lifted to satisfy, at equality, the point
     236  // (xu,yu,wu=xu*yu) (resp. (xl,yl,wl=xl*yl)), where xl and xu are
     237  // lower and upper bound of x, etc.
     238  //
     239  //        alpha (w - l) + y - l/xt >= -l/(xt^2) (x-xt) ...
     240  //
     241  // where alpha is such that the relation holds at equality at the
     242  // point (xu,yu,xu*yu):
     243  //
     244  //    alpha = [-yu + l/xt - l/(xt^2)(xu-xt)] / (xu*yu - l)
     245
     246  if ((x0 > xl + COUENNE_EPS) && (y0 > yl + COUENNE_EPS) &&
     247      (x0 < xu + COUENNE_EPS) && (y0 < yu + COUENNE_EPS)) {
     248
     249    if ((wl > 0) && (x0*y0 < wl)) { // that is, if (x0,y0) is out of the contour
     250
     251      CouNumber xyl = xl * yl;
     252
     253      // first and third orthant
     254      if      ((xyl <  wl) && (xu*yu >=wl)) contourCut (cg,cs, x0,y0, wl, +1, xl,yl, xu,yu, xi,yi,wi);
     255      else if ((xyl >= wl) && (xu*yu < wl)) contourCut (cg,cs, x0,y0, wl, -1, xu,yu, xl,yl, xi,yi,wi);
    156256    }
    157     else
    158       // check points C and D: fourth orthant intersections
    159       if ((xl*yu > wu) && (xu*yl <= wu)) {
    160         addImplTangent (cg, cs, xl, yl, wu, xi, yi, -1, -1); // C
    161         addImplTangent (cg, cs, xu, yu, wu, xi, yi, +1, -1); // D
    162       }
    163   }
    164   else
    165     if (wl > COUENNE_EPS) {
    166       // check points A and B: third orthant intersections
    167       if ((xl*yl >= wl) && (xu*yu < wl)) {
    168         addImplTangent (cg, cs, xl, yu, wl, xi, yi, -1, -1); // A
    169         addImplTangent (cg, cs, xu, yl, wl, xi, yi, +1, -1); // B
    170       }
    171       else
    172         // check points C and D: first orthant intersections
    173         if ((xu*yu >= wl) && (xl*yl < wl)) {
    174           addImplTangent (cg, cs, xl, yu, wl, xi, yi, +1, +1); // C
    175           addImplTangent (cg, cs, xu, yl, wl, xi, yi, -1, +1); // D
    176         }
     257
     258  // Similarly for w <= u < 0
     259
     260    if ((wu < 0) && (x0*y0 > wu)) { // that is, if (x0,y0) is out of the contour
     261
     262      CouNumber xuyl = xl * yu;
     263
     264      // second and fourth orthant
     265      if      ((xuyl > wu) && (xl*yu <=wu)) contourCut (cg,cs, x0,y0, wu, +1, xu,yl, xl,yu, xi,yi,wi);
     266      else if ((xuyl <=wu) && (xl*yu > wu)) contourCut (cg,cs, x0,y0, wu, -1, xl,yu, xu,yl, xi,yi,wi);
    177267    }
     268  }
    178269}
    179 
    180 
    181 /// add tangent to feasibility region of w=x*y at intersection with bounding box
    182 
    183 void addImplTangent (const CouenneCutGenerator *cg, OsiCuts &cs,
    184                      CouNumber xb,    // x bound
    185                      CouNumber yb,    // y
    186                      CouNumber wb,    // w
    187                      int xi, int yi,  // indices of variables passed to createCut
    188                      int sign_check,  // get maximum or minimum of abscissa
    189                      int sign_ineq) { // lower or upper half-plane
    190 
    191   CouNumber xp, yp; // intersection point
    192 
    193   //  point is (xb, wb/xb) if xb*yb > wb, (wb/yb,yb) otherwise
    194 
    195   CouNumber check = xb*yb - wb; // violation (signed)
    196 
    197   if (sign_check < 0) // intersection point depends on violation check
    198     check = -check;
    199 
    200   if (check > 0) {xp = xb;    yp = wb/xb;}
    201   else           {xp = wb/yb; yp = yb;}
    202 
    203   // infinite or null bounds give vertical or horizontal cuts, useless
    204   if ((fabs (xp) < COUENNE_EPS) ||
    205       (fabs (yp) < COUENNE_EPS) ||
    206       (fabs (xp) > COUENNE_INFINITY) ||
    207       (fabs (yp) > COUENNE_INFINITY))
    208     return;
    209 
    210   CouNumber w_xp = wb / xp;
    211 
    212   cg -> createCut (cs, yp+w_xp, sign_ineq, yi, 1., xi, w_xp/xp);
    213 }
  • branches/Couenne/Couenne/src/convex/operators/powNewton.cpp

    r534 r556  
    1010#include <CouenneTypes.h>
    1111
    12 #define MAX_ITER 100
     12#define MAX_ITER 10
    1313#define COU_POW_TOLERANCE 1e-12
    1414
     
    3939
    4040  // Newton loop. Tolerance is set above
    41   for (register int k = MAX_ITER; (fabs (F) > COU_POW_TOLERANCE) && k--;) {
     41  for (int k = MAX_ITER; k--;) {
    4242
    4343    xk -= F / Fp;
     
    4646    fpk = fp (xk);
    4747    F   = xk - xc + fpk * fk;
     48    if (fabs (F) > COU_POW_TOLERANCE) break;
    4849    Fp  = 1 + fpp (xk) * fk + fpk * fpk;
    4950  }
  • branches/Couenne/Couenne/src/expression/expression.h

    r545 r556  
    231231
    232232  /// set up branching object by evaluating many branching points for
    233   /// each expression's arguments
    234   CouNumber selectBranch (expression *w,
    235                           const OsiBranchingInformation *info,
    236                           int &ind,
    237                           double * &brpts,
    238                           int &way)
     233  /// each expression's arguments. Return estimated improvement in
     234  /// objective function
     235  virtual CouNumber selectBranch (expression *w,
     236                                  const OsiBranchingInformation *info,
     237                                  int &ind,
     238                                  double * &brpts,
     239                                  int &way)
    239240    {ind = -1; return 0.;}
    240241
     
    249250
    250251
    251 /// updates maximum violation. Used with all impliedBound. Returns 1
    252 /// if a bound has been modified, 0 otherwise
     252/// updates maximum violation. Used with all impliedBound. Returns true
     253/// if a bound has been modified, false otherwise
    253254
    254255inline bool updateBound (int sign, CouNumber *dst, CouNumber src) {
  • branches/Couenne/Couenne/src/problem/impliedBounds.cpp

    r553 r556  
    3939    //  (auxiliaries_ [i] -> Image () -> code () == COU_EXPRGROUP))
    4040
    41     CouNumber l0 = lb_ [nvar+i],
    42               u0 = ub_ [nvar+i];
     41    //    CouNumber l0 = lb_ [nvar+i],
     42    //              u0 = ub_ [nvar+i];
    4343
    4444   if (auxiliaries_ [i] -> Image () -> impliedBound (nvar+i, lb_, ub_, chg_bds) > COUENNE_EPS) {
  • branches/Couenne/Couenne/src/problem/problem.cpp

    r497 r556  
    235235
    236236  // copy arrays (do not simply make x_ point to x)
    237   for (register int i = (n==-1) ? nvars : n; i--;) {
    238 
    239     x_  [i] = x [i];
    240     lb_ [i] = l [i];
    241     ub_ [i] = u [i];
    242   }
     237
     238  if (x) for (register int i = (n==-1) ? nvars : n; i--;)  x_  [i] = x [i];
     239  if (l) for (register int i = (n==-1) ? nvars : n; i--;)  lb_ [i] = l [i];
     240  if (u) for (register int i = (n==-1) ? nvars : n; i--;)  ub_ [i] = u [i];
    243241
    244242  expression::update (x_, lb_, ub_);
Note: See TracChangeset for help on using the changeset viewer.