Changeset 557


Ignore:
Timestamp:
May 31, 2007 12:31:39 PM (12 years ago)
Author:
pbelotti
Message:

use a more elaborate structure to check changed bounds. Check bounds of other variables after applying OBBT for one. Remove exit from exprSub::impliedBound. Return changed bounds in exprPow. Added explicit handling of whichWay in CouenneObject?. Added branching rule for expr{Inv,Log} and exprMul (infinite bounds).

Location:
branches/Couenne/Couenne
Files:
45 edited

Legend:

Unmodified
Added
Removed
  • branches/Couenne/Couenne/TODO

    r552 r557  
     1- add cuts at subsequent calls at changed bounds also
    12- avoid very large coefficients
    23- Bilinear cuts
  • branches/Couenne/Couenne/src/branch/CouenneBranchingObject.cpp

    r556 r557  
    4343    alpha = CouenneBranchingObject::Alpha ();
    4444
    45   //  printf ("///////////brpt = %g, x = %Lf, l,u = [%Lf,%Lf]\n", brpoint, x, l, u);
    46 
    4745  if (fabs (brpoint) < COUENNE_INFINITY)
    4846    x = brpoint;
     
    7876    else                             value_ = 0;
    7977  else value_ = x;
    80 
    81 
    82   /*  if ((x > l + COUENNE_NEAR_BOUND) &&
    83       (x < u - COUENNE_NEAR_BOUND))      // x is not at the boundary
    84     value_ = x;
    85   else // current point is at one of the bounds
    86     if ((l > - COUENNE_INFINITY) &&
    87         (u <   COUENNE_INFINITY))
    88       // finite bounds, apply midpoint rule
    89       value_ = alpha * x + (1. - alpha) * (l + u) / 2.;
    90 
    91     else
    92       // infinite (one direction) bound interval, x is at the boundary
    93       // push it inwards
    94       // TODO: look for a proper displacement
    95       if (fabs (x-l) < COUENNE_EPS) value_ = l + (1+fabs (l)) / 2.;
    96       else                          value_ = u - (1+fabs (u)) / 2.;   */
    9778
    9879  if (0) {
     
    131112  }
    132113
    133   //if (0) printf (" [%.6e,%.6e]", l, u);
    134   //  if (way) solver -> setColLower (index_, reference_ -> isInteger () ? ceil  (value_) : value_);
    135   //  else     solver -> setColUpper (index_, reference_ -> isInteger () ? floor (value_) : value_);
    136 
    137114  if (!way) solver -> setColUpper (index_, integer_ ? floor (value_) : value_); // down branch
    138115  else      solver -> setColLower (index_, integer_ ? ceil  (value_) : value_); // up   branch
  • branches/Couenne/Couenne/src/branch/CouenneChooseVariable.cpp

    r556 r557  
    2222CouenneChooseVariable::CouenneChooseVariable (const CouenneChooseVariable &source):
    2323  OsiChooseVariable (source),
    24   problem_ (source.Problem ()) {}
     24  problem_ (source.problem_) {}
    2525
    2626/// Assignment operator
    2727CouenneChooseVariable & CouenneChooseVariable::operator= (const CouenneChooseVariable& rhs)
    28 {problem_ = rhs.Problem(); return *this;}
     28{problem_ = rhs.problem_; return *this;}
    2929
    3030/// Clone
  • branches/Couenne/Couenne/src/branch/CouenneObject.cpp

    r556 r557  
    2525
    2626  // whichWay should be set to which branch first (for two-way branching?)
     27  // if selectBranch not called, choose one at random
     28  whichWay_ = whichWay = TWO_RAND;
    2729
    2830  // infeasibility is always null for linear expressions
     
    110112       brVarInd_, brPts_, whichWay); // result: who, where, and how to branch
    111113
     114    whichWay_ = whichWay;
     115
    112116    if (brVarInd_ >= 0)
    113117      delta = improv;
     
    128132  if (0) {
    129133
    130     printf ("Inf |%+.4e - %+.4e| = %+.4e  %+.4e ",  ////[%.2f,%.2f]
     134    printf ("Inf |%+.4e - %+.4e| = %+.4e  %+.4e way %d, ind %d",  ////[%.2f,%.2f]
    131135            var, expr,
    132136            //      expression::Lbound (reference_ -> Index ()),
    133137            //      expression::Ubound (reference_ -> Index ()),
    134             fabs (var - expr), delta);
     138            fabs (var - expr), delta, whichWay, brVarInd_);
    135139
    136140    reference_             -> print (std::cout); std::cout << " = ";
     
    211215
    212216  // way has suggestion from CouenneObject::infeasibility()
     217  // Well, not exactly... it seems
    213218
    214219  bool isint = reference_ -> isInteger ();
     220
     221  way = whichWay_;
    215222
    216223  if (brVarInd_ >= 0) // if applied latest selectBranching
     
    225232    case THREE_RIGHT:
    226233    case THREE_RAND:
    227       //printf ("3 WAY Branch x%d at %g--%g [%d] (%d)\n", brVarInd_, *brPts_, brPts_ [1], way, isint);
     234      //printf ("3Way Branch x%d @ %g ][ %g [%d] (%d)\n", brVarInd_, *brPts_, brPts_ [1], way, isint);
    228235      return new CouenneThreeWayBranchObj (brVarInd_, brPts_ [0], brPts_ [1], way, isint);
    229236    default:
    230       printf ("CouenneObject::createBranch(): way has no sense\n");
     237      printf ("CouenneObject::createBranch(): way=%d has no sense\n", way);
    231238      exit (-1);
    232239    }
  • branches/Couenne/Couenne/src/branch/CouenneObject.hpp

    r537 r557  
    2121
    2222enum {TWO_LEFT,                 TWO_RIGHT,   TWO_RAND,
    23       THREE_LEFT, THREE_CENTER, THREE_RIGHT, THREE_RAND};
     23      THREE_LEFT, THREE_CENTER, THREE_RIGHT, THREE_RAND, BRANCH_NONE};
    2424
    2525
     
    3636    reference_ (ref),
    3737    brVarInd_  (-1),
    38     brPts_     (NULL) {}
     38    brPts_     (NULL),
     39    whichWay_  (BRANCH_NONE) {}
    3940
    4041  /// Destructor
     
    9495  /// ThreeWayBranching. Ends with a -COIN_DBL_MAX (not considered...)
    9596  mutable CouNumber *brPts_;
     97
     98  /// how to branch. This should be done automatically from
     99  /// ::infeasibility() and ::createBranch(), but for some reason
     100  /// obj->whichWay() in the last 20 lines of CbcNode.cpp returns 0,
     101  /// therefore we set our own whichWay_ within this object and use it
     102  /// between the two methods.
     103  mutable int whichWay_;
    96104};
    97105
  • branches/Couenne/Couenne/src/branch/CouenneThreeWayBranchObj.cpp

    r556 r557  
    8989  //        1 if ">= b"  node
    9090
    91   int way;//, ind = reference_ -> Index ();
     91  int way;
    9292
    9393  switch (branchIndex_) {
     
    120120    }
    121121  }
    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]);*/
    128122
    129123  branchIndex_++;
  • branches/Couenne/Couenne/src/branch/operators/branchExprAbs.cpp

    r556 r557  
    4040  way = TWO_RAND; // don't really care which subtree to visit first
    4141
    42   // exact distance from current point to the two subsequent
     42  // exact distance between current point and the two subsequent
    4343  // convexifications
    4444  return mymin (project (1., -1., 0., x0, y0, 0., COUENNE_INFINITY,  0, NULL, NULL),
  • branches/Couenne/Couenne/src/branch/operators/branchExprDiv.cpp

    r556 r557  
    1414#include <projections.h>
    1515
    16 
    17 #define BR_NEXT_ZERO 1e-3
    18 #define BR_MULT      1e-3
    1916
    2017/// set up branching object by evaluating many branching points for
  • branches/Couenne/Couenne/src/branch/operators/branchExprExp.cpp

    r537 r557  
    2323                                 double * &brpts,
    2424                                 int &way) {
    25   ind = -1;
    26   return 0.;
    2725
    2826  // two cases: inside or outside the belly.
     
    3331  // branch.
    3432  //
    35   // Outside: it suffices to project the current point on the line to
    36   // get the maxi-min displacement.
     33  // Outside: it suffices to project the current point on the line
     34  // (i.e. to get the closest point on the line) to get the maxi-min
     35  // displacement.
    3736  //
    38   // As happens for all monotonous functions, after chosing *brpts it
    39   // is equivalent to choose w's or x's index as ind, as the implied
    40   // bounds will do the rest.
     37  // As for all monotonous functions, after choosing *brpts it is
     38  // equivalent to choose w's or x's index as ind, as the implied- and
     39  // propagated bounds will do the rest.
    4140
    42   ind = argument_ -> Index ();
    43   int wi = w -> Index ();
     41  ind    = argument_ -> Index ();
     42  int wi = w         -> Index ();
    4443
    45   if ((ind < 0) || (wi < 0)) {printf ("Couenne, w=f(x): negative index\n"); exit (-1);}
     44  if ((ind < 0) || (wi < 0)) {printf ("Couenne, w=exp(x): negative index\n"); exit (-1);}
    4645
    47   brpts = (double *) realloc (brpts, sizeof (double));
     46  CouNumber y0 = info -> solution_ [wi],
     47            x0 = info -> solution_ [ind],
     48            l  = info -> lower_    [ind],
     49            u  = info -> upper_    [ind];
    4850
    49   CouNumber x0 = info -> solution_ [ind],
    50             y0 = info -> solution_ [wi];
     51  if (y0 < exp (x0)) {
    5152
    52   if (y0 < exp (x0)) { // A) Outside
     53    // Outside
     54
     55    brpts = (double *) realloc (brpts, sizeof (double));
    5356
    5457    *brpts = powNewton (x0, y0, exp, exp, exp);
     58
     59    if      (*brpts < l) *brpts = (u <   COUENNE_INFINITY) ? ((l+u)/2) : l + 1;
     60    else if (*brpts > u) *brpts = (l > - COUENNE_INFINITY) ? ((l+u)/2) : u - 1;
     61
    5562    way = TWO_RAND;
    5663    CouNumber dy = y0 - exp (*brpts);
     
    5865    return sqrt (x0*x0 + dy*dy);
    5966
    60   } else { // B) Inside
     67  } else {
    6168
    62     CouNumber l = info -> lower_ [ind],
    63               u = info -> upper_ [ind];
    64 
    65     // two cases:
     69    // Inside. Two cases:
    6670 
    6771    if ((l < -COUENNE_INFINITY) &&
    6872        (u >  COUENNE_INFINITY)) {
    6973
    70       // B1) bounds are infinite
     74      // 1) both bounds are infinite --> three way branching
    7175
    72       brpts = (double *) realloc (brpts, 2*sizeof (double));
    73       way = THREE_CENTER;
    74       *brpts = x0;
    75       brpts [1] = log (y0);
     76      brpts = (double *) realloc (brpts, 2 * sizeof (double));
     77      way = THREE_CENTER; // focus on central convexification first
     78
     79      brpts [0] = x0;       // draw vertical   from (x0,y0) south to curve y=exp(x)
     80      brpts [1] = log (y0); //      horizontal              east
    7681
    7782      CouNumber a = y0 - exp (x0), // sides of a triangle with (x0,y0)
     
    8388    } else {
    8489
     90      // 2) at least one of them is finite --> two way branching
     91
    8592      brpts = (double *) realloc (brpts, sizeof (double));
    86 
    87       // B2) at least one of them is finite
    8893
    8994      if (l < -COUENNE_INFINITY) {
    9095
    9196        *brpts = x0;
     97        if (*brpts > u - COUENNE_NEAR_BOUND) *brpts = u-1;
    9298        way = TWO_RIGHT;
    9399        return y0 - exp (x0);
     
    96102
    97103        *brpts = log (y0);
     104        if (*brpts < l + COUENNE_NEAR_BOUND) *brpts = l+1;
    98105        way = TWO_LEFT;
    99106        return x0 - log (y0);
     
    102109
    103110        *brpts = powNewton (x0, y0, exp, exp, exp);
     111
     112        if ((*brpts > u - COUENNE_NEAR_BOUND) ||
     113            (*brpts < l + COUENNE_NEAR_BOUND))
     114          *brpts = (l+u) / 2;
     115
    104116        way = TWO_RAND;
    105117
  • branches/Couenne/Couenne/src/branch/operators/branchExprInv.cpp

    r537 r557  
    22 * Name:    branchExprInv.cpp
    33 * Author:  Pietro Belotti
    4  * Purpose: return branch gain and branch object for 1/x
     4 * Purpose: return branch selection for 1/x
    55 *
    66 * (C) Pietro Belotti. This file is licensed under the Common Public License (CPL)
     
    88
    99#include <exprInv.h>
     10#include <exprDiv.h>
     11#include <exprPow.h>
    1012#include <CouennePrecisions.h>
    1113#include <CouenneTypes.h>
     
    1719/// set up branching object by evaluating many branching points for
    1820/// each expression's arguments
     21
    1922CouNumber exprInv::selectBranch (expression *w,
    2023                                 const OsiBranchingInformation *info,
     
    2225                                 double * &brpts,
    2326                                 int &way) {
    24   ind = -1;
    25   return 0.;
    2627
    27   /*ind = argument_ -> Index ();
     28  // two cases: inside or outside the bellies.
     29  //
     30  // Inside: the distance depends on the projection of the current
     31  // point onto the would-be upper envelopes, which forces us to look
     32  // at it numerically. If both bounds are infinite, create a ThreeWay
     33  // branch.
     34  //
     35  // Outside: it suffices to project the current point on the line
     36  // (i.e. to get the closest point on the line) to get the maxi-min
     37  // displacement.
     38  //
     39  // As for all monotonous functions, after choosing *brpts it is
     40  // equivalent to choose w's or x's index as ind, as the implied- and
     41  // propagated bounds will do the rest.
    2842
    29   if (ind < 0) {
    30     printf ("argument of exprAbs has negative index\n");
    31     exit (-1);
     43  ind    = argument_ -> Index ();
     44  int wi = w         -> Index ();
     45
     46  if ((ind < 0) || (wi < 0)) {printf ("Couenne, w=exp(x): negative index\n"); exit (-1);}
     47
     48  CouNumber y0 = info -> solution_ [wi],
     49            x0 = info -> solution_ [ind],
     50            l  = info -> lower_    [ind],
     51            u  = info -> upper_    [ind];
     52
     53  // two cases:
     54  //
     55  // 1) bounds include 0: three way branch to exclude 0 (refer to branchExprDiv.cpp)
     56  // 2) otherwise         two   way branch
     57
     58  if ((l < 0) && (u > 0)) {
     59
     60    brpts = (double *) realloc (brpts, 2 * sizeof (CouNumber));
     61    way = THREE_RIGHT;
     62
     63    brpts [0] = (l >= - BR_NEXT_ZERO - COUENNE_EPS) ? (l * BR_MULT) : - BR_NEXT_ZERO;
     64    brpts [1] = (u <=   BR_NEXT_ZERO + COUENNE_EPS) ? (u * BR_MULT) :   BR_NEXT_ZERO;
     65
     66    return ((fabs (x0) < COUENNE_EPS) ? 1. :
     67            (1 / x0 - info -> solution_ [wi]));
    3268  }
    3369
    34   brpts = (double *) malloc (sizeof (double));
    35   *brpts = 0.;
    36   way = TWO_LEFT;
     70  // case 2: look if inside or outside of belly (refer to branchExprExp.cpp)
    3771
    38   return mymin (project (1., -1., 0., x0, y0, 0., COUENNE_INFINITY,  0, NULL, NULL),
    39   project (1., +1., 0., x0, y0, -COUENNE_INFINITY, 0., 0, NULL, NULL));*/
    40 }
    4172
    42 /*
    43 /// distance covered by current point if branching rule applied to this expression
    44 double exprInv::BranchGain (expression *w, const OsiBranchingInformation *info) {
     73  if (x0*y0 < 1) { // outside bellies
    4574
    46   int xi = argument_ -> Index (),
    47       wi = w         -> Index ();
     75    brpts = (double *) realloc (brpts, sizeof (double));
    4876
    49   if (xi < 0 || wi < 0) {
    50     printf ("negative indices inbranchExprAbs\n");
    51     exit (-1);
     77    *brpts = powNewton (x0, y0, inv, oppInvSqr, inv_dblprime);
     78
     79    if      (*brpts < l) *brpts = (u <   COUENNE_INFINITY) ? ((l+u)/2) : l + 1;
     80    else if (*brpts > u) *brpts = (l > - COUENNE_INFINITY) ? ((l+u)/2) : u - 1;
     81
     82    way = (u < 0) ? TWO_RIGHT : TWO_LEFT; // explore finite interval first
     83
     84    CouNumber dy = y0 - 1. / *brpts;
     85    x0 -= *brpts;
     86    return sqrt (x0*x0 + dy*dy);
     87
    5288  }
    5389
    54   CouNumber x0 = expression::Variable (xi),
    55             y0 = expression::Variable (wi);
     90  // Inside. Two cases:
     91 
     92  if ((l <   COUENNE_EPS) && (u >   COUENNE_INFINITY) ||
     93      (u > - COUENNE_EPS) && (l < - COUENNE_INFINITY)) {
    5694
    57   return 0;
     95    // 1) bounds are infinite both horizontally and vertically
     96    //    (i.e. [-inf,0] or [0,inf]) --> three way branching
    5897
    59   //  return mymin (project (1, -1, 0, x0, y0, 0, COUENNE_INFINITY,  0, NULL, NULL),
    60   //            project (1, +1, 0, x0, y0, -COUENNE_INFINITY, 0, 0, NULL, NULL));
     98    brpts = (double *) realloc (brpts, 2 * sizeof (double));
     99    way = THREE_CENTER; // focus on central convexification first
     100
     101    brpts [0] = x0;        // draw vertical   from (x0,y0) south (north) to curve y=1/x
     102    brpts [1] = 1. / (y0); //      horizontal              west  (east)
     103
     104    CouNumber a = y0 - 1 / x0, // sides of a triangle with (x0,y0)
     105      b = x0 - 1 / y0, // as one of the vertices
     106      c = a * cos (atan (a/b));
     107
     108    return mymin (a, mymin (b, c));
     109  }
     110
     111  // 2) at least one of them is finite --> two way branching
     112
     113  brpts = (double *) realloc (brpts, sizeof (double));
     114
     115  if (l < - COUENNE_INFINITY) { // and u away from -0
     116
     117    *brpts = x0;
     118    if (*brpts > u - COUENNE_NEAR_BOUND) *brpts = u-1;
     119    way = TWO_RIGHT;
     120    return y0 - 1. / x0;
     121  }
     122
     123  if (u > COUENNE_INFINITY) { // and l away from +0
     124
     125    *brpts = x0;
     126    if (*brpts < l + COUENNE_NEAR_BOUND) *brpts = l+1;
     127    way = TWO_LEFT;
     128    return y0 - 1. / x0;
     129  }
     130
     131  // last case: nice finite interval and limited curve
     132
     133  *brpts = powNewton (x0, y0, inv, oppInvSqr, inv_dblprime);
     134
     135  if ((*brpts > u - COUENNE_NEAR_BOUND) ||
     136      (*brpts < l + COUENNE_NEAR_BOUND))
     137    *brpts = (l+u) / 2;
     138
     139  way = TWO_RAND;
     140
     141  CouNumber dx = x0 - *brpts,
     142    dy = y0 - 1. / *brpts;
     143
     144  return sqrt (dx*dx + dy*dy);
    61145}
    62 
    63 /// branching object best suited for this expression
    64 OsiBranchingObject *exprInv::BranchObject (expression *w, const OsiBranchingInformation *) {
    65 
    66   //  return new CouenneBranchingObject (Argument (), 0);
    67   return NULL;
    68 }
    69 */
  • branches/Couenne/Couenne/src/branch/operators/branchExprLog.cpp

    r537 r557  
    88
    99#include <exprLog.h>
     10#include <exprPow.h>
     11
    1012#include <CouennePrecisions.h>
    1113#include <CouenneTypes.h>
     
    2224                                 double * &brpts,
    2325                                 int &way) {
    24   ind = -1;
    25   return 0.;
    2626
    27   /*  ind = argument_ -> Index ();
     27  // quite similar to exprExp::selectBranch() (see branchExprExp.cpp)
    2828
    29   if (ind < 0) {
    30     printf ("argument of exprAbs has negative index\n");
    31     exit (-1);
     29  // two cases: inside or outside the belly.
     30  //
     31  // Inside: the distance depends on the projection of the current
     32  // point onto the would-be upper envelopes, which forces us to look
     33  // at it numerically. If both bounds are infinite, create a ThreeWay
     34  // branch.
     35  //
     36  // Outside: it suffices to project the current point on the line
     37  // (i.e. to get the closest point on the line) to get the maxi-min
     38  // displacement.
     39  //
     40  // As for all monotonous functions, after choosing *brpts it is
     41  // equivalent to choose w's or x's index as ind, as the implied- and
     42  // propagated bounds will do the rest.
     43
     44  ind    = argument_ -> Index ();
     45  int wi = w         -> Index ();
     46
     47  if ((ind < 0) || (wi < 0)) {printf ("Couenne, w=log(x): negative index\n"); exit (-1);}
     48
     49  CouNumber y0 = info -> solution_ [wi],
     50    x0 = info -> solution_ [ind],
     51    l  = info -> lower_    [ind],
     52    u  = info -> upper_    [ind];
     53
     54  if (y0 > log (x0)) {
     55
     56    if (x0 == 0) x0 = 1e-30;
     57
     58    // Outside
     59
     60    brpts = (double *) realloc (brpts, sizeof (double));
     61
     62    *brpts = powNewton (x0, y0, log, inv, oppInvSqr);
     63
     64    if      (*brpts < l) *brpts = (u < COUENNE_INFINITY) ? ((l+u)/2) : (10*l + 1);
     65    else if (*brpts > u) *brpts = (l+u) / 2;
     66
     67    way = TWO_LEFT;
     68    CouNumber dy = y0 - log (*brpts);
     69    x0 -= *brpts;
     70    return sqrt (x0*x0 + dy*dy);
     71
     72  }
     73
     74  // Inside. Two cases:
     75 
     76  if ((l < COUENNE_EPS * COUENNE_EPS) &&
     77      (u > COUENNE_INFINITY)) {
     78
     79    // 1) curve is unlimited in both senses --> three way branching
     80
     81    brpts = (double *) realloc (brpts, 2 * sizeof (double));
     82    way = THREE_CENTER; // focus on central convexification first
     83
     84    brpts [0] = exp (y0); // draw horizontal from (x0,y0) south to curve y=log(x)
     85    brpts [1] = x0;       //      vertical                east
     86
     87    CouNumber a = x0 - exp (y0), // sides of a triangle with (x0,y0)
     88      b = y0 - log (x0), // as one of the vertices
     89      c = a * cos (atan (a/b));
     90
     91    return mymin (a, mymin (b, c));
     92
     93  }
     94
     95  // 2) at least one of them is finite --> two way branching
     96
     97  brpts = (double *) realloc (brpts, sizeof (double));
     98
     99  if (l < COUENNE_EPS * COUENNE_EPS) {
     100
     101    *brpts = exp (y0);
     102    if ((*brpts > u - COUENNE_NEAR_BOUND) ||
     103        (*brpts < l + COUENNE_NEAR_BOUND))
     104      *brpts = (l+u) / 2;
     105
     106    way = TWO_RIGHT;
     107    return mymin (x0 - exp (y0), y0 - log (x0));
     108
    32109  }
     110 
     111  if (u > COUENNE_INFINITY) {
    33112
    34   brpts = (double *) malloc (sizeof (double));
    35   *brpts = 0.;
    36   way = TWO_LEFT;
     113    *brpts = x0;
     114    if (*brpts < l + COUENNE_NEAR_BOUND) *brpts = l+1;
     115    way = TWO_LEFT;
     116    return y0 - log (x0);
    37117
    38   return mymin (project (1., -1., 0., x0, y0, 0., COUENNE_INFINITY,  0, NULL, NULL),
    39   project (1., +1., 0., x0, y0, -COUENNE_INFINITY, 0., 0, NULL, NULL));*/
     118  }
     119
     120  // both are finite
     121
     122  *brpts = powNewton (x0, y0, log, inv, oppInvSqr);
     123
     124  if ((*brpts > u - COUENNE_NEAR_BOUND) ||
     125      (*brpts < l + COUENNE_NEAR_BOUND))
     126    *brpts = (l+u) / 2;
     127
     128  way = TWO_RAND;
     129
     130  CouNumber
     131    dx = x0 - *brpts,
     132    dy = y0 - log (*brpts);
     133
     134  return sqrt (dx*dx + dy*dy);
    40135}
    41 
    42 /*
    43 /// distance covered by current point if branching rule applied to this expression
    44 double exprLog::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 0;
    58   //  return mymin (project (1, -1, 0, x0, y0, 0, COUENNE_INFINITY,  0, NULL, NULL),
    59   //            project (1, +1, 0, x0, y0, -COUENNE_INFINITY, 0, 0, NULL, NULL));
    60 }
    61 
    62 /// branching object best suited for this expression
    63 OsiBranchingObject *exprLog::BranchObject (expression *w, const OsiBranchingInformation *) {
    64 
    65   // branching once on 0 is sufficient for expressions w=|x|
    66   //  return new CouenneBranchingObject (Argument (), 0);
    67   return NULL;
    68 }
    69 */
  • branches/Couenne/Couenne/src/branch/operators/branchExprMul.cpp

    r537 r557  
    22 * Name:    branchExprMul.cpp
    33 * Author:  Pietro Belotti
    4  * Purpose: return branch gain and branch object for multiplications
     4 * Purpose: return branch data for multiplications
    55 *
    66 * (C) Pietro Belotti. This file is licensed under the Common Public License (CPL)
     
    1515
    1616
     17#define LARGE_BOUND 9.9e12
     18#define BOUND_WING 100
     19
    1720/// set up branching object by evaluating many branching points for
    1821/// each expression's arguments
     
    2225                                 double * &brpts,
    2326                                 int &way) {
     27  int xi = arglist_ [0] -> Index (),
     28      yi = arglist_ [1] -> Index ();
     29
     30  if ((xi < 0) || (yi < 0)) {
     31    printf ("arguments of exprMul have negative index\n");
     32    exit (-1);
     33  }
     34
     35  CouNumber x0 = expression::Variable (xi),   y0 = expression::Variable (yi),
     36            xl = expression::Lbound   (xi),   yl = expression::Lbound   (yi), 
     37            xu = expression::Ubound   (xi),   yu = expression::Ubound   (yi); 
     38
     39  // First, try to avoid infinite bounds for multiplications, which
     40  // make them pretty hard to deal with
     41
     42  if ((xl < - COUENNE_INFINITY) &&
     43      (xu >   COUENNE_INFINITY)) { // first
     44
     45    ind = xi;
     46
     47    // branch around current point. If it is also at a crazy value,
     48    // reset it close to zero.
     49
     50    brpts = (double *) malloc (2 * sizeof (double));
     51    CouNumber curr = x0;
     52
     53    if (fabs (curr) >= LARGE_BOUND) curr = 0;
     54
     55    brpts [0] = curr - BOUND_WING;
     56    brpts [1] = curr + BOUND_WING;
     57
     58    way = THREE_CENTER;
     59
     60    return COUENNE_INFINITY;
     61  }
     62
     63  if ((yl < - COUENNE_INFINITY) &&
     64      (yu >   COUENNE_INFINITY)) { // and second factor
     65
     66    ind = yi;
     67
     68    // branch around current point. If it is also at a crazy value,
     69    // reset it close to zero.
     70
     71    brpts = (double *) malloc (2 * sizeof (double));
     72    CouNumber curr = y0;
     73
     74    if (fabs (curr) >= LARGE_BOUND) curr = 0;
     75
     76    brpts [0] = curr - BOUND_WING;
     77    brpts [1] = curr + BOUND_WING;
     78
     79    way = THREE_CENTER;
     80
     81    return COUENNE_INFINITY;
     82  }
     83
     84  // there is at least one finite corner of the bounding box
    2485
    2586  ind = -1;
    2687  return 0.;
    27 
    28   //  ind = argument_ -> Index ();
    29 
    30   /*if (ind < 0) {
    31     printf ("argument of exprAbs has negative index\n");
    32     exit (-1);
    33   }
    34 
    35   brpts = (double *) malloc (sizeof (double));
    36   *brpts = 0.;
    37   way = TWO_LEFT;
    38 
    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));*/
    4188}
    42 
    43 /*
    44 /// distance covered by current point if branching rule applied to this expression
    45 double exprMul::BranchGain (expression *w, const OsiBranchingInformation *info) {
    46 
    47   int xi = arglist_ [0] -> Index (),
    48       yi = arglist_ [1] -> Index (),
    49       wi = w            -> Index ();
    50 
    51   CouNumber x0 = expression::Variable (xi),
    52             y0 = expression::Variable (yi),
    53             w0 = expression::Variable (wi);
    54 
    55   return 0;
    56   //  return mymin (project (1, -1, 0, x0, y0, 0, COUENNE_INFINITY,  0, NULL, NULL),
    57   //            project (1, +1, 0, x0, y0, -COUENNE_INFINITY, 0, 0, NULL, NULL));
    58 }
    59 
    60 /// branching object best suited for this expression
    61 OsiBranchingObject *exprMul::BranchObject (expression *w, const OsiBranchingInformation *) {
    62 
    63   // branching once on 0 is sufficient for expressions w=|x|
    64   //  return new CouenneBranchingObject (Argument (), 0);
    65   return NULL;
    66 
    67 }
    68 */
  • branches/Couenne/Couenne/src/convex/CouenneCutGenerator.h

    r556 r557  
    180180  /// tighten bounds using propagation, implied bounds and reduced costs
    181181  bool boundTightening (const OsiSolverInterface &, OsiCuts &,
    182                         char *, Bonmin::BabInfo * = NULL) const;
     182                        t_chg_bounds *, Bonmin::BabInfo * = NULL) const;
    183183
    184184  /// Optimality Based Bound Tightening
    185   int obbt (const OsiSolverInterface &, OsiCuts &, char *, Bonmin::BabInfo *) const;
     185  int obbt (OsiSolverInterface *, OsiCuts &, t_chg_bounds *, Bonmin::BabInfo *) const;
    186186};
    187187
  • branches/Couenne/Couenne/src/convex/boundTightening.cpp

    r553 r557  
    2222bool CouenneCutGenerator::boundTightening (const OsiSolverInterface &si,
    2323                                           OsiCuts &cs,
    24                                            char *chg_bds,
     24                                           t_chg_bounds *chg_bds,
    2525                                           Bonmin::BabInfo * babInfo) const {
    2626
     
    4545      //      printf ("updating upper %g <-- %g\n", primal0, primal);
    4646      problem_ -> Ub (objInd) = UB;
    47       chg_bds [objInd] = 1;
     47      chg_bds [objInd]. upper = CHANGED;
    4848    }
    4949   
     
    5151        (LB   > dual0 + COUENNE_EPS)) { // update dual bound
    5252      problem_ -> Lb (objInd) = LB;
    53       chg_bds [objInd] = 1;
     53      chg_bds [objInd]. lower = CHANGED;
    5454    }
    5555
     
    8080            // can improve bound on variable w_i
    8181            problem_ -> Ub (i) = x + (UB-LB) / rc;
    82             chg_bds [i] = 1;
     82            chg_bds [i].upper = CHANGED;
    8383          }
    8484        }
     
    133133  // expression structure is not a tree.
    134134
     135  // TODO: limit should depend on number of constraints, that is,
     136  // bound transmission should be documented and the cycle should stop
     137  // as soon as no new constraint subgraph has benefited from bound
     138  // transmission.
     139
    135140  return true;
    136141}
  • branches/Couenne/Couenne/src/convex/generateCuts.cpp

    r556 r557  
    5252
    5353// translate changed bound sparse array into a dense one
    54 void sparse2dense (int ncols, char *chg_bds, int *&changed, int &nchanged) {
     54void sparse2dense (int ncols, t_chg_bounds *chg_bds, int *&changed, int &nchanged) {
    5555
    5656  // convert sparse chg_bds in something handier
     
    5858  nchanged = 0;
    5959
    60   for (register int i=ncols, j=0; i--; j++)
    61     if (*chg_bds++) {
     60  for (register int i=ncols, j=0; i--; j++, chg_bds++)
     61    if (chg_bds -> lower || chg_bds -> upper) {
    6262      *changed++ = j;
    6363      nchanged++;
     
    9191  // used with malloc/realloc/free
    9292
    93   char *chg_bds = new char [ncols];
     93  t_chg_bounds *chg_bds = new t_chg_bounds [ncols];
    9494
    9595  // fill it with zeros
    96   for (register int i = ncols; i--;)
    97     *chg_bds++ = 0;
     96  for (register int i = ncols; i--; chg_bds++)
     97    chg_bds -> lower = chg_bds -> upper = UNCHANGED;
    9898  chg_bds -= ncols;
    9999
     
    166166          const double * nowUpper = si.getColUpper();
    167167
    168           for (register int i=0; i < ncols; i++)
    169 
    170             if (beforeLower && (nowLower [i] >= beforeLower [i] + COUENNE_EPS) ||
    171                 beforeUpper && (nowUpper [i] <= beforeUpper [i] - COUENNE_EPS))
    172               chg_bds [i] = 1;
     168          for (register int i=0; i < ncols; i++) {
     169
     170            if (beforeLower && (nowLower [i] >= beforeLower [i] + COUENNE_EPS))
     171              chg_bds [i].lower = CHANGED;
     172            if (beforeUpper && (nowUpper [i] <= beforeUpper [i] - COUENNE_EPS))
     173              chg_bds [i].upper = CHANGED;
     174          }
    173175
    174176        } else printf ("WARNING: could not access parent's bounds\n");
     
    227229  if ((!firstcall_ || (info.pass > 0)) &&
    228230      (CoinDrand48 () < (double) COU_OBBT_CUTOFF_LEVEL / (info.level + 1))) {
     231
    229232    // apply OBBT at all levels up to the COU_OBBT_CUTOFF_LEVEL-th,
    230233    // and then with probability inversely proportional to the level
    231234
     235    OsiSolverInterface *csi = si.clone (true);
     236
     237    dynamic_cast <OsiClpSolverInterface *> (csi) -> setupForRepeatedUse ();
     238
    232239    int nImprov, nIter = 0;
    233240    bool repeat = true;
     
    235242    while (repeat &&
    236243           (nIter++ < MAX_OBBT_ITER) &&
    237            ((nImprov = obbt (si, cs, chg_bds, babInfo)) > 0))
     244           ((nImprov = obbt (csi, cs, chg_bds, babInfo)) > 0))
    238245
    239246      if (nImprov >= THRES_NBD_CHANGED) {
     
    244251        sparse2dense (ncols, chg_bds, changed, nchanged);
    245252       
    246         genColCuts (si, cs, nchanged, changed);
     253        genColCuts (*csi, cs, nchanged, changed);
    247254
    248255        int nCurCuts = cs.sizeRowCuts ();
    249         genRowCuts (si, cs, nchanged, changed);// !!!
     256        genRowCuts (*csi, cs, nchanged, changed);// !!!
    250257        repeat = nCurCuts < cs.sizeRowCuts (); // reapply only if new cuts available
    251258      }
     259
     260    delete csi;
    252261
    253262    if (nImprov < 0)
  • branches/Couenne/Couenne/src/convex/obbt.cpp

    r556 r557  
    1313
    1414#define OBBT_EPS 1e-3
     15#define MAX_OBBT_LP_ITERATION 100
     16
    1517
    1618/// reoptimize and change bound of a variable if needed
     
    2022                  bool isint) {            /// is this variable integer
    2123
     24  csi -> setObjSense (sense);
    2225  csi -> setDblParam (OsiDualObjectiveLimit,
    2326                      (sense == 1) ? COUENNE_INFINITY : -COUENNE_INFINITY);
    2427
    25   csi -> setObjSense (sense);
     28  csi -> setDblParam (OsiPrimalObjectiveLimit, bound);
    2629  csi -> resolve ();
    2730
     
    3841
    3942
    40 /// Optimality based bound tightening
    41 int CouenneCutGenerator::obbt (const OsiSolverInterface &si,
    42                                OsiCuts &cs,
    43                                char *chg_bds,
    44                                Bonmin::BabInfo * babInfo) const {
    45   int
    46     nImprov = 0,
    47     ncols   = si.getNumCols (),
    48     objind  = problem_ -> Obj (0) -> Body () -> Index ();
     43/// visit all variables (in one sense) and apply OBBT to each
     44int obbt_stage (const CouenneCutGenerator *cg,
     45                OsiSolverInterface *csi,
     46                OsiCuts &cs,
     47                t_chg_bounds *chg_bds,
     48                Bonmin::BabInfo * babInfo,
     49                int sense) {
     50
     51  CouenneProblem *p = cg -> Problem ();
     52
     53  int ncols   = csi -> getNumCols (),
     54      objind  = p -> Obj (0) -> Body () -> Index (),
     55      nImprov = 0;
    4956
    5057  double *objcoe = (double *) malloc (ncols * sizeof (double));
     
    5562  objcoe -= ncols;
    5663
    57   // create clone of current solver interface
    58   // TODO: clone only at pass 0 and delete at the end
    59   OsiSolverInterface *csi = si.clone (true);
    60 
    61   // apply all (row+column) cuts to formulation
    62   csi -> applyCuts (cs);
    63 
    6464  // for all (original+auxiliary) variables x_i,
     65  ////////////////////////////////////////////////////
    6566  for (int i=0; i<ncols; i++)
    6667    //  for (int i=0; i<problem_ ->nVars(); i++)
    6768
    68     if (i != objind) { // do not improve objective's bounds
     69    if ((i != objind) &&
     70        ((sense > 0) && !(chg_bds [i].lower & EXACT) ||
     71         (sense < 0) && !(chg_bds [i].upper & EXACT))) { // do not improve objective's bounds
    6972
    70       int nOrig = problem_ -> nVars ();
    71       bool chg  = false,
    72         isInt  = (i < nOrig) ?
    73         (problem_ -> Var (i)       -> isInteger ()) :
    74         (problem_ -> Aux (i-nOrig) -> isInteger ());
     73      int nOrig  = p -> nVars ();
     74
     75      bool isInt = (i < nOrig) ?
     76        (p -> Var (i)       -> isInteger ()) :
     77        (p -> Aux (i-nOrig) -> isInteger ());
    7578
    7679      objcoe [i] = 1;
    7780      csi -> setObjective (objcoe);
    7881
    79       // minimize...
     82      CouNumber &bound = (sense == 1) ? (p -> Lb (i)) : (p -> Ub (i));
    8083
    81       if (updateBound (csi, 1, problem_ -> Lb (i), isInt)) {
    82         csi -> setColLower (i, problem_ -> Lb (i));
    83         chg = true;
    84       }
    85      
    86       // ...and then maximize x_i on csi.
     84      // m{in,ax}imize xi on csi
    8785
    88       if (updateBound (csi, -1, problem_ -> Ub (i), isInt)) {
    89         csi -> setColUpper  (i, problem_ -> Ub (i));
    90         chg = true;
    91       }
     86      if (updateBound (csi, sense, bound, isInt)) {
    9287
    93       if (chg) {
     88        const double *sol = csi -> getColSolution ();
    9489
    95         if (!(boundTightening (si, cs, chg_bds, babInfo))) {
    96           delete csi;
     90        if (sense==1) {csi -> setColLower (i, bound); chg_bds [i].lower |= CHANGED | EXACT;}
     91        else          {csi -> setColUpper (i, bound); chg_bds [i].upper |= CHANGED | EXACT;}
     92
     93        // check value and bounds of other variables
     94
     95        for (int j=0; j<ncols; j++)
     96          if ((j!=i) && (j!=objind)) {
     97
     98            if (sol [j] <= p -> Lb (j) + COUENNE_EPS) chg_bds [j].lower |= EXACT;
     99            if (sol [j] >= p -> Ub (j) - COUENNE_EPS) chg_bds [j].upper |= EXACT;
     100          }
     101
     102        // re-check considering reduced costs (more expensive)
     103
     104        CouNumber *redcost = NULL;
     105
     106        // first, compute reduced cost when c = c - e_i, where e_i is
     107        // a vector with all zero except a one in position i. This
     108        // serves as a base to compute modified reduced cost below.
     109
     110        if (0)
     111        for (int j=0; j<ncols; j++)
     112          if ((j!=i) && (j!=objind)) {
     113
     114            // fake a change in the objective function and compute
     115            // reduced cost. If resulting vector is all positive
     116            // (negative), this solution is also optimal for the
     117            // minimization (maximization) of x_j and the
     118            // corresponding chg_bds[j] can be set to EXACT.
     119
     120            if (sol [j] <= p -> Lb (j) + COUENNE_EPS) chg_bds [j].lower |= EXACT;
     121            if (sol [j] >= p -> Ub (j) - COUENNE_EPS) chg_bds [j].upper |= EXACT;
     122          }
     123
     124        // re-apply bound tightening
     125
     126        if (!(cg -> boundTightening (*csi, cs, chg_bds, babInfo)))
    97127          return -1; // tell caller this is infeasible
    98         }
    99128
    100         chg_bds [i] = 1;
    101129        nImprov++;
    102130      }
     
    106134    }
    107135
    108   delete csi;
    109136  return nImprov;
    110137}
     138
     139/// Optimality based bound tightening
     140int CouenneCutGenerator::obbt (OsiSolverInterface *csi,
     141                               OsiCuts &cs,
     142                               t_chg_bounds *chg_bds,
     143                               Bonmin::BabInfo * babInfo) const {
     144  int nImprov;
     145
     146  csi -> setIntParam (OsiMaxNumIteration, MAX_OBBT_LP_ITERATION);
     147
     148  // apply all (row+column) cuts to formulation
     149  csi -> applyCuts (cs);
     150
     151  nImprov =  obbt_stage (this, csi, cs, chg_bds, babInfo,  1); // minimization
     152  nImprov += obbt_stage (this, csi, cs, chg_bds, babInfo, -1); // maximization
     153
     154  return nImprov;
     155}
  • branches/Couenne/Couenne/src/expression/CouenneTypes.h

    r545 r557  
    3434enum convexity {NONCONVEX, CONVEX, CONCAVE, AFFINE};
    3535
     36/** (un)changed bounds */
     37
     38#define UNCHANGED 0
     39#define CHANGED   1
     40#define EXACT     2
     41
     42typedef struct chg_bounds {
     43  char lower;
     44  char upper;
     45} t_chg_bounds;
     46
    3647typedef double CouNumber;
    3748typedef CouNumber (*unary_function) (CouNumber);
  • branches/Couenne/Couenne/src/expression/exprCopy.h

    r543 r557  
    121121
    122122  /// implied bound processing
    123   bool impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg)
     123  bool impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg)
    124124    {return copy_ -> impliedBound (wind, l, u, chg);}
    125125};
  • branches/Couenne/Couenne/src/expression/exprVar.cpp

    r476 r557  
    5454/// implied bound processing. Expression w = x, upon change in lower
    5555/// or upper bound of w, whose index is wind
    56 bool exprVar::impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg) {
     56bool exprVar::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    5757
    58   bool res = updateBound (-1, l + varIndex_, l [wind]);
    59   res      = updateBound (+1, u + varIndex_, u [wind]) || res;
     58  bool res;
    6059
    61   if (res)
    62     chg [varIndex_] = 1;
     60  if (updateBound (-1, l + varIndex_, l [wind])) {res = true; chg [varIndex_].lower = CHANGED;}
     61  if (updateBound (+1, u + varIndex_, u [wind])) {res = true; chg [varIndex_].upper = CHANGED;}
    6362
    6463  return res;
  • branches/Couenne/Couenne/src/expression/exprVar.h

    r543 r557  
    9999
    100100  /// implied bound processing
    101   virtual bool impliedBound (int, CouNumber *, CouNumber *, char *);
     101  virtual bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *);
    102102
    103103  /// rank of an original variable is always one
  • branches/Couenne/Couenne/src/expression/expression.h

    r556 r557  
    223223  /// bound. The method returns the best bound improvement obtained on
    224224  /// all variables of the expression.
    225   virtual bool impliedBound (int, CouNumber *, CouNumber *, char *)
     225  virtual bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *)
    226226    {return false;}
    227227
  • branches/Couenne/Couenne/src/expression/operators/exprAbs.cpp

    r543 r557  
    5555
    5656
    57 /// printing
    58 
    59 //void exprAbs::print (std::ostream& out) const {
    60 //  exprUnary::print (out, "abs", PRE);
    61   /*
    62   out << "|";
    63   argument_ -> print (out);
    64   out << "|";
    65   */
    66 //}
    67 
    68 
    6957/// implied bound processing for expression w = |x|, upon change in
    7058/// lower- and/or upper bound of w, whose index is wind
    7159
    72 bool exprAbs::impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg) {
     60bool exprAbs::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    7361
    7462  int index = argument_ -> Index ();
     
    8472  if (wl > 0) {
    8573
    86     if (*xl > 0) tighter = updateBound (-1, xl,  wl);
    87     if (*xu < 0) tighter = updateBound (+1, xu, -wl) || tighter;
     74    if      (*xl > 0) {if (updateBound (-1, xl,  wl)) {tighter = true; chg [index].lower = CHANGED;}}
     75    else if (*xu < 0) {if (updateBound (+1, xu, -wl)) {tighter = true; chg [index].upper = CHANGED;}}
    8876  }
    8977
    9078  // w <= u (if u < 0 the problem is infeasible)
    9179
    92   tighter = updateBound (-1, xl, -wu) || tighter;
    93   tighter = updateBound (+1, xu,  wu) || tighter;
    94 
    95   if (tighter)
    96     chg [index] = 1;
     80  if (updateBound (-1, xl, -wu)) {tighter = true; chg [index].lower = CHANGED;}
     81  if (updateBound (+1, xu,  wu)) {tighter = true; chg [index].upper = CHANGED;}
    9782
    9883  return tighter;
  • branches/Couenne/Couenne/src/expression/operators/exprAbs.h

    r543 r557  
    5353
    5454  /// implied bound processing
    55   bool impliedBound (int, CouNumber *, CouNumber *, char *);
     55  bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *);
    5656
    5757  /// set up branching object by evaluating many branching points for
  • branches/Couenne/Couenne/src/expression/operators/exprCos.h

    r551 r557  
    3838    {return "cos";}
    3939
    40   /// print position
    41   //  virtual const std::string printPos ()
    42   //    {return PRE;}
    43 
    44 
    45   /// print "cos" and argument
    46   //  void print (std::ostream&) const;
    47 
    4840  /// obtain derivative of expression
    4941  expression *differentiate (int index);
     
    6759/// convex envelope for sine/cosine
    6860
    69 //void addHexagon (const CouenneCutGenerator *, // pointer to the caller cut generator
    70 //               OsiCuts &,      // cut set to be enriched
    71 //               unary_function, // sine or cosine
    72 //               exprAux *,      // auxiliary variable
    73 //               expression *);  // argument of cos/sin (should be a variable)
    74 
    7561#endif
  • branches/Couenne/Couenne/src/expression/operators/exprDiv.h

    r543 r557  
    1212#include <exprOp.h>
    1313#include <CouennePrecisions.h>
     14
     15#define BR_NEXT_ZERO 1e-3
     16#define BR_MULT      1e-3
    1417
    1518
     
    7881
    7982  /// implied bound processing
    80   bool impliedBound (int, CouNumber *, CouNumber *, char *);
     83  bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *);
    8184
    8285  /// set up branching object by evaluating many branching points for
  • branches/Couenne/Couenne/src/expression/operators/exprExp.cpp

    r543 r557  
    11/*
    2  * Name:    exprExp.C
     2 * Name:    exprExp.cpp
    33 * Author:  Pietro Belotti
    44 * Purpose: definition of the exponential
     
    4444/// implied bound processing for expression w = exp(x), upon change in
    4545/// lower- and/or upper bound of w, whose index is wind
    46 bool exprExp::impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg) {
     46bool exprExp::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    4747
    48   bool res = false;
     48  bool resU, resL = resU = false;
    4949  int ind = argument_ -> Index ();
    5050
     
    5252
    5353  if ((b = l [wind]) >= COUENNE_EPS) // lower bound
    54     res = updateBound (-1, l + ind, log (b));
     54    resL = updateBound (-1, l + ind, log (b));
    5555
    5656  if ((b = u [wind]) >= COUENNE_EPS) // upper bound
    57     res = updateBound ( 1, u + ind, log (b)) || res;
     57    resU = updateBound ( 1, u + ind, log (b));
    5858  else if (b < - COUENNE_EPS) {
    5959    // make it infeasible
    60     res = updateBound ( 1, u + ind, -1) || res;
    61     res = updateBound (-1, l + ind,  1) || res;
     60    resU = updateBound ( 1, u + ind, -1) || resU;
     61    resL = updateBound (-1, l + ind,  1) || resL;
    6262  }
    6363
    64   if (res)
    65     chg [ind] = 1;
     64  if (resL) chg [ind].lower = CHANGED;
     65  if (resU) chg [ind].upper = CHANGED;
    6666
    67   return res;
     67  return (resL || resU);
    6868}
  • branches/Couenne/Couenne/src/expression/operators/exprExp.h

    r543 r557  
    6464
    6565  /// implied bound processing
    66   bool impliedBound (int, CouNumber *, CouNumber *, char *);
     66  bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *);
    6767
    6868  /// set up branching object by evaluating many branching points for
  • branches/Couenne/Couenne/src/expression/operators/exprInv.cpp

    r545 r557  
    3939
    4040
    41 /// general function (see below)
    42 bool invPowImplBounds (int, int, CouNumber *, CouNumber *, CouNumber);
     41/// general function to tighten implied bounds of a function w = x^k,
     42/// k negative, integer or inverse integer, and odd
     43
     44void invPowImplBounds (int wind, int index,
     45                       CouNumber *l, CouNumber *u, CouNumber k,
     46                       bool &resL, bool &resU) {
     47
     48  CouNumber wl = l [wind],
     49            wu = u [wind];
     50
     51  // 0 <= l <= w <= u
     52
     53  if (wl >= 0.) {
     54    if (wu > COUENNE_EPS) {
     55      if (wu < COUENNE_INFINITY - 1) resL = updateBound (-1, l + index, pow (wu, k));
     56      else                           resL = updateBound (-1, l + index, 0.);
     57    }
     58    if (wl > COUENNE_EPS)            resU = updateBound (+1, u + index, pow (wl, k));
     59  }
     60
     61  // l <= w <= u <= 0
     62
     63  if (wu <= -0.) {
     64    if (wl < - COUENNE_EPS) {
     65      if (wl > - COUENNE_INFINITY + 1) resU = updateBound (+1, u + index, pow (wl, k)) || resU;
     66      else                             resU = updateBound (+1, u + index, 0.)          || resU;
     67    }
     68    if (wu < - COUENNE_EPS)            resL = updateBound (-1, l + index, pow (wu, k)) || resL;
     69  }
     70}
    4371
    4472
     
    4674/// lower- and/or upper bound of w, whose index is wind
    4775
    48 bool exprInv::impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg) {
     76bool exprInv::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    4977
    5078  // Expression w = 1/x: we can only improve the bounds if
     
    5785  int index = argument_ -> Index ();
    5886
    59   bool res = invPowImplBounds (wind, index, l, u, -1.);
     87  bool resL, resU = resL = false;
    6088
    61   if (res)
    62     chg [index] = 1;
     89  invPowImplBounds (wind, index, l, u, -1., resL, resU);
    6390
    64   return res;
     91  if (resL) chg [index].lower = CHANGED;
     92  if (resU) chg [index].upper = CHANGED;
     93
     94  return (resL || resU);
    6595}
    6696
    67 
    68 /// general function to tighten implied bounds of a function w = x^k,
    69 /// k negative, integer or inverse integer, and odd
    70 
    71 bool invPowImplBounds (int wind, int index, CouNumber *l, CouNumber *u, CouNumber k) {
    72 
    73   CouNumber wl = l [wind],
    74             wu = u [wind];
    75 
    76   bool res = false;
    77 
    78   // 0 <= l <= w <= u
    79 
    80   if (wl >= 0.) {
    81     if (wu > COUENNE_EPS) {
    82       if (wu < COUENNE_INFINITY - 1) res = updateBound (-1, l + index, pow (wu, k));
    83       else                           res = updateBound (-1, l + index, 0.);
    84     }
    85     if (wl > COUENNE_EPS)            res = updateBound (+1, u + index, pow (wl, k)) || res;
    86   }
    87 
    88   // l <= w <= u <= 0
    89 
    90   if (wu <= -0.) {
    91     if (wl < - COUENNE_EPS) {
    92       if (wl > - COUENNE_INFINITY + 1) res = updateBound (+1, u + index, pow (wl, k)) || res;
    93       else                             res = updateBound (+1, u + index, 0.)          || res;
    94     }
    95     if (wu < - COUENNE_EPS)            res = updateBound (-1, l + index, pow (wu, k)) || res;
    96   }
    97 
    98   return res;
    99 }
  • branches/Couenne/Couenne/src/expression/operators/exprInv.h

    r543 r557  
    6868
    6969  /// implied bound processing
    70   bool impliedBound (int, CouNumber *, CouNumber *, char *);
     70  bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *);
    7171
    7272  /// set up branching object by evaluating many branching points for
  • branches/Couenne/Couenne/src/expression/operators/exprLog.cpp

    r543 r557  
    6161
    6262
    63 /// printing
    64 
    65 //void exprLog::print (std::ostream& out) const
    66 //{exprUnary::print (out, "log", PRE);}
    67 
    68 
    6963/// implied bound processing for expression w = log(x), upon change in
    7064/// lower- and/or upper bound of w, whose index is wind
    7165
    72 bool exprLog::impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg) {
     66bool exprLog::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    7367
    7468  int ind = argument_ -> Index ();
     69  bool res;
    7570
    76   bool res = updateBound (-1, l + ind, exp (l [wind]));
    77   res      = updateBound ( 1, u + ind, exp (u [wind])) || res;
    78 
    79   if (res) {
    80 
    81     /*printf ("w_%d [%g,%g] -------> x_%d in [%g,%g] ",
    82             wind, l [wind], u [wind],
    83             ind,  l [ind],  u [ind]);*/
    84     chg [ind] = 1;
    85   }
     71  if (updateBound (-1, l + ind, exp (l [wind]))) {res = true; chg [ind].lower = CHANGED;}
     72  if (updateBound ( 1, u + ind, exp (u [wind]))) {res = true; chg [ind].upper = CHANGED;}
    8673
    8774  return res;
  • branches/Couenne/Couenne/src/expression/operators/exprLog.h

    r543 r557  
    6262                     OsiCuts &cs, const CouenneCutGenerator *cg);
    6363
    64   ///
     64  /// code for comparisons
    6565  virtual enum expr_type code () {return COU_EXPRINV;}
    6666
    6767  /// implied bound processing
    68   bool impliedBound (int, CouNumber *, CouNumber *, char *);
     68  bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *);
    6969
    7070  /// set up branching object by evaluating many branching points for
  • branches/Couenne/Couenne/src/expression/operators/exprMul.h

    r543 r557  
    6767
    6868  /// implied bound processing
    69   bool impliedBound (int, CouNumber *, CouNumber *, char *);
     69  bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *);
    7070
    7171  /// set up branching object by evaluating many branching points for
  • branches/Couenne/Couenne/src/expression/operators/exprOpp.cpp

    r543 r557  
    2828
    2929
    30 // printing
    31 
    32 //void exprOpp::print (std::ostream& out) const
    33 //{exprUnary::print (out, "-", PRE);}
    34 
    35 
    3630/// implied bound processing for expression w = -x, upon change in
    3731/// lower- and/or upper bound of w, whose index is wind
    3832
    39 bool exprOpp::impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg) {
     33bool exprOpp::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    4034
    4135  int ind = argument_ -> Index ();
     36  bool res;
    4237
    43   bool res = updateBound (-1, l + ind, - u [wind]);
    44   res      = updateBound ( 1, u + ind, - l [wind]) || res;
    45 
    46   if (res)
    47     chg [ind] = 1;
     38  if (updateBound (-1, l + ind, - u [wind])) {res = true; chg [ind].lower = CHANGED;}
     39  if (updateBound ( 1, u + ind, - l [wind])) {res = true; chg [ind].upper = CHANGED;}
    4840
    4941  return res;
  • branches/Couenne/Couenne/src/expression/operators/exprOpp.h

    r543 r557  
    5757                     OsiCuts &cs, const CouenneCutGenerator *cg);
    5858
    59   ///
     59  /// code for comparisons
    6060  virtual enum expr_type code () {return COU_EXPROPP;}
    6161
    6262  /// implied bound processing
    63   bool impliedBound (int, CouNumber *, CouNumber *, char *);
     63  bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *);
    6464};
    6565
  • branches/Couenne/Couenne/src/expression/operators/exprPow.cpp

    r543 r557  
    160160/// inverse integer, and odd
    161161
    162 bool invPowImplBounds (int, int, CouNumber *, CouNumber *, CouNumber);
     162void invPowImplBounds (int, int, CouNumber *, CouNumber *, CouNumber, bool &, bool &);
    163163
    164164
     
    166166/// lower- and/or upper bound of w, whose index is wind
    167167
    168 bool exprPow::impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg) {
    169 
    170   bool res = false;
     168bool exprPow::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
     169
     170  bool resL, resU = resL = false;
    171171
    172172  if (arglist_ [0] -> Type () <= CONST)   // base is constant or zero, nothing to do
     
    203203    if (k > 0.) { // simple, just follow bounds
    204204
    205       if (wl > - COUENNE_INFINITY)
    206         res    = updateBound (-1, l + index, pow (wl, 1./k));
    207       else res = updateBound (-1, l + index, - COUENNE_INFINITY);
    208 
    209       if (wu < COUENNE_INFINITY)
    210         res    = updateBound (+1, u + index, pow (wu, 1./k)) || res;
    211       else res = updateBound (+1, u + index, COUENNE_INFINITY);
     205      resL = (wl > - COUENNE_INFINITY) ?
     206        updateBound (-1, l + index, pow (wl, 1./k)) :
     207        updateBound (-1, l + index, - COUENNE_INFINITY);
     208
     209      resU = (wu < COUENNE_INFINITY) ?
     210        updateBound (+1, u + index, pow (wu, 1./k)) :
     211        updateBound (+1, u + index, COUENNE_INFINITY);
    212212
    213213    } else // slightly more complicated, resort to same method as in exprInv
    214       res = invPowImplBounds (wind, index, l, u, 1./k);
     214      invPowImplBounds (wind, index, l, u, 1./k, resL, resU);
    215215  }
    216216  else
     
    225225
    226226        if (fabs (bound) < COUENNE_INFINITY) {
    227           res = updateBound (-1, l + index, - pow (bound, 1./k));
    228           res = updateBound (+1, u + index,   pow (bound, 1./k)) || res;
     227          resL = updateBound (-1, l + index, - pow (bound, 1./k));
     228          resU = updateBound (+1, u + index,   pow (bound, 1./k));
    229229        } else {
    230           res = updateBound (-1, l + index, - COUENNE_INFINITY);
    231           res = updateBound (+1, u + index,   COUENNE_INFINITY) || res;
     230          resL = updateBound (-1, l + index, - COUENNE_INFINITY);
     231          resU = updateBound (+1, u + index,   COUENNE_INFINITY);
    232232        }
    233233      }
     
    240240        ub = wl;
    241241      }
    242        
    243       if (lb > COUENNE_EPS) res = updateBound (-1, l + index, pow (lb, 1./k));
     242
     243      if (lb > COUENNE_EPS) resL = updateBound (-1, l + index, pow (lb, 1./k));
    244244
    245245      if (fabs (ub) < COUENNE_INFINITY) {
    246         if (ub > COUENNE_EPS) res = updateBound (+1, u + index, pow (ub, 1./k)) || res;
    247       } else                  res = updateBound (+1, u + index, COUENNE_INFINITY) || res;       
    248     }
    249 
    250   return res;
    251 }
     246        if (ub > COUENNE_EPS) resU = updateBound (+1, u + index, pow (ub, 1./k));
     247      } else                  resU = updateBound (+1, u + index, COUENNE_INFINITY);
     248    }
     249
     250  if (resL) chg [index].lower = CHANGED;
     251  if (resU) chg [index].upper = CHANGED;
     252
     253  return (resL || resU);
     254}
  • branches/Couenne/Couenne/src/expression/operators/exprPow.h

    r543 r557  
    7575
    7676  /// implied bound processing
    77   bool impliedBound (int, CouNumber *, CouNumber *, char *);
     77  bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *);
    7878
    7979  /// set up branching object by evaluating many branching points for
  • branches/Couenne/Couenne/src/expression/operators/exprSub.cpp

    r543 r557  
    9292/// lower- and/or upper bound of w, whose index is wind
    9393
    94 bool exprSub::impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg) {
     94bool exprSub::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    9595
    9696  // caution, xi or yi might be -1
     
    112112  bool res = false;
    113113
    114   return false;
     114  /// !!! REMOVED A return false;
    115115
    116116  // w >= l
    117117
    118   if ((xi >= 0) && (res = updateBound (-1, l + xi, yl + wl)))        chg [xi] = 1;
    119   if ((yi >= 0) && (res = updateBound (+1, u + yi, xu - wl) || res)) chg [yi] = 1;
     118  if ((xi >= 0) && (updateBound (-1, l + xi, yl + wl))) {res = true; chg [xi].lower = CHANGED;}
     119  if ((yi >= 0) && (updateBound (+1, u + yi, xu - wl))) {res = true; chg [yi].upper = CHANGED;}
    120120
    121121  // w <= u
    122122
    123   if ((xi >= 0) && (res = updateBound (+1, u + xi, yu + wu) || res)) chg [xi] = 1;
    124   if ((yi >= 0) && (res = updateBound (-1, l + yi, xl - wu) || res)) chg [yi] = 1;
     123  if ((xi >= 0) && (updateBound (+1, u + xi, yu + wu))) {res = true; chg [xi].upper = CHANGED;}
     124  if ((yi >= 0) && (updateBound (-1, l + yi, xl - wu))) {res = true; chg [yi].lower = CHANGED;}
    125125
    126126  return res;
  • branches/Couenne/Couenne/src/expression/operators/exprSub.h

    r543 r557  
    7373
    7474  /// implied bound processing
    75   bool impliedBound (int, CouNumber *, CouNumber *, char *);
     75  bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *);
    7676};
    7777
  • branches/Couenne/Couenne/src/expression/operators/exprSum.h

    r543 r557  
    6464
    6565  /// implied bound processing
    66   virtual bool impliedBound (int, CouNumber *, CouNumber *, char *);
     66  virtual bool impliedBound (int, CouNumber *, CouNumber *, t_chg_bounds *);
    6767
    6868  /// compute best variable to branch on (nonsense here, as there is
  • branches/Couenne/Couenne/src/expression/operators/impliedBounds-exprDiv.cpp

    r524 r557  
    1515/// lower- and/or upper bound of w, whose index is wind
    1616
    17 bool exprDiv::impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg) {
     17bool exprDiv::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    1818
    1919  bool resx, resy = resx = false;
     
    4040    if (c > COUENNE_EPS) {
    4141
    42       resx = updateBound (-1, l + ind, l [wind] * c);
    43       resx = updateBound ( 1, u + ind, u [wind] * c) || resx;
     42      if (updateBound (-1, l + ind, l [wind] * c)) {resx = true; chg [ind].lower = CHANGED;}
     43      if (updateBound ( 1, u + ind, u [wind] * c)) {resx = true; chg [ind].upper = CHANGED;}
    4444    }
    4545    else if (c < - COUENNE_EPS) {
    4646
    47       resx = updateBound (-1, l + ind, u [wind] * c);
    48       resx = updateBound ( 1, u + ind, l [wind] * c) || resx;
     47      if (updateBound (-1, l + ind, u [wind] * c)) {resx = true; chg [ind].lower = CHANGED;}
     48      if (updateBound ( 1, u + ind, l [wind] * c)) {resx = true; chg [ind].upper = CHANGED;}
    4949    }
    50 
    51     if (resx)
    52       chg [ind] = 1;
    53 
    5450  } else {
    55 
    56     int xi = arglist_ [0] -> Index (),
    57         yi = arglist_ [1] -> Index ();
    58 
    59     CouNumber x0;
    6051
    6152    // deal with all other cases
     
    9384    // may improve).
    9485 
     86    int xi = arglist_ [0] -> Index (),
     87        yi = arglist_ [1] -> Index ();
     88
     89    CouNumber x0;
     90
    9591    CouNumber *xl = l + xi, *yl = l + yi, wl = l [wind],
    9692              *xu = u + xi, *yu = u + yi, wu = u [wind];
     
    115111    }
    116112
     113    bool resxL, resxU, resyL,
     114      resyU = resxL = resxU = resyL = false;
     115
    117116    // general case
    118117
     
    121120      // point C: (xl,yl)
    122121
    123       resy = (*yl<0) && (*yl > *xl/wl) && updateBound (-1, yl, 0) || resy;//
     122      resyL = (*yl<0) && (*yl > *xl/wl) && updateBound (-1, yl, 0) || resyL; //
    124123
    125124      if ((*yl>0) && (*yl < *xl/wl)) { // point C violates x/y >= wl, down
    126         resx = updateBound (-1, xl, *yu*wl) || resx;//
    127         resy = updateBound (-1, yl, *xu/wl) || resy;//
     125        resxL = updateBound (-1, xl, *yu*wl) || resxL; //
     126        resyL = updateBound (-1, yl, *xu/wl) || resyL; //
    128127      }
    129128
     
    131130
    132131      if ((*yu<0) && (*yu > *xu/wl)) { // point B violates x/y >= wl, down
    133         resx = updateBound (+1, xu, *yl*wl) || resx;
    134         resy = updateBound (+1, yu, *xl/wl) || resy;
    135       }
    136 
    137       resy = (*yu>0) && (*yu < *xu/wl) && updateBound (+1, yu, 0) || resy;
     132        resxU = updateBound (+1, xu, *yl*wl) || resxU;
     133        resyU = updateBound (+1, yu, *xl/wl) || resyU;
     134      }
     135
     136      resyU = (*yu>0) && (*yu < *xu/wl) && updateBound (+1, yu, 0) || resyU;
    138137
    139138    } else if (wl >   COUENNE_EPS) { // w >= wl, wl positive
     
    141140      //
    142141
    143       resy = (*yl<0) && (*yl < *xl/wl) && updateBound (-1, yl, mymin (*xl/wl, 0)) || resy;
    144       resx = (*yl>0) && (*yl > *xl/wl) && updateBound (-1, xl, *yl*wl)            || resx;
    145 
    146       //
    147 
    148       resx = (*yu<0) && (*yu < *xu/wl) && updateBound (+1, xu, *yu*wl)            || resx;
    149       resy = (*yu>0) && (*yu > *xu/wl) && updateBound (+1, yu, mymax (*xu/wl, 0)) || resy;
     142      resyL = (*yl<0) && (*yl < *xl/wl) && updateBound (-1, yl, mymin (*xl/wl, 0)) || resyL;
     143      resxL = (*yl>0) && (*yl > *xl/wl) && updateBound (-1, xl, *yl*wl)            || resxL;
     144
     145      //
     146
     147      resxU = (*yu<0) && (*yu < *xu/wl) && updateBound (+1, xu, *yu*wl)            || resxU;
     148      resyU = (*yu>0) && (*yu > *xu/wl) && updateBound (+1, yu, mymax (*xu/wl, 0)) || resyU;
    150149    }
    151150
     
    156155      //
    157156
    158       resy = (*yl<0) && (*yl > *xu/wu) && updateBound (-1, yl, 0)      || resy;
     157      resyL = (*yl<0) && (*yl > *xu/wu) && updateBound (-1, yl, 0)      || resyL;
    159158
    160159      if ((*yl>0) && (*yl < *xu/wu)) {
    161         resx = updateBound (+1, xu, *yu*wu) || resx;
    162         resy = updateBound (-1, yl, *xl/wu) || resy;
     160        resxU = updateBound (+1, xu, *yu*wu) || resxU;
     161        resyL = updateBound (-1, yl, *xl/wu) || resyL;
    163162      }
    164163
     
    166165
    167166      if ((*yu<0) && (*yu > *xl/wu)) {
    168         resx = updateBound (-1, xl, *yl*wu) || resx;
    169         resy = updateBound (+1, yu, *xu/wu) || resy;
    170       }
    171 
    172       resy = (*yu>0) && (*yu < *xl/wu) && updateBound (+1, yu, 0)      || resy;
     167        resxL = updateBound (-1, xl, *yl*wu) || resxL;
     168        resyU = updateBound (+1, yu, *xu/wu) || resyU;
     169      }
     170
     171      resyU = (*yu>0) && (*yu < *xl/wu) && updateBound (+1, yu, 0)      || resyU;
    173172
    174173    } else if (wu < - COUENNE_EPS) { // w <= wu, wu negative
     
    176175      //
    177176
    178       resy = (*yl<0) && (*yl < *xu/wu) && updateBound (-1, yl, mymin (*xu/wu,0))  || resy;//
    179       resx = (*yu<0) && (*yu < *xl/wu) && updateBound (-1, xl, *yu*wu)            || resx;
    180 
    181       //
    182 
    183       resx = (*yl>0) && (*yl > *xu/wu) && updateBound (+1, xu, *yl*wu)            || resx;
    184       resy = (*yu>0) && (*yu > *xl/wu) && updateBound (+1, yu, mymax (*xl/wu,0))  || resy;
    185     }
    186 
    187     if (resx) chg [xi] = 1;
    188     if (resy) chg [yi] = 1;
     177      resyL = (*yl<0) && (*yl < *xu/wu) && updateBound (-1, yl, mymin (*xu/wu,0))  || resyL;//
     178      resxL = (*yu<0) && (*yu < *xl/wu) && updateBound (-1, xl, *yu*wu)            || resxL;
     179
     180      //
     181
     182      resxU = (*yl>0) && (*yl > *xu/wu) && updateBound (+1, xu, *yl*wu)            || resxU;
     183      resyU = (*yu>0) && (*yu > *xl/wu) && updateBound (+1, yu, mymax (*xl/wu,0))  || resyU;
     184    }
     185
     186    if (resxL) chg [xi].lower = CHANGED;
     187    if (resxU) chg [xi].upper = CHANGED;
     188    if (resyL) chg [yi].lower = CHANGED;
     189    if (resyU) chg [yi].upper = CHANGED;
    189190
    190191    /*if (resx || resy)
     
    192193              wind, wl, wu, xi, *xl, *xu, yi, *yl, *yu);
    193194              else printf ("                                                 \r");*/
     195
     196    resx = resxL || resxU;
     197    resy = resyL || resyU;
    194198  }
    195199
  • branches/Couenne/Couenne/src/expression/operators/impliedBounds-exprMul.cpp

    r525 r557  
    1414/// lower- and/or upper bound of w, whose index is wind
    1515
    16 bool exprMul::impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg) {
     16bool exprMul::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    1717
    18   bool res = false;
     18  bool resL, resU = resL = false;
    1919  int ind;
    2020
     
    3838
    3939     
    40       res = (l [wind] > - COUENNE_INFINITY) && updateBound (-1, l + ind, l [wind] / c);
    41       res = (u [wind] <   COUENNE_INFINITY) && updateBound ( 1, u + ind, u [wind] / c) || res;
     40      resL = (l [wind] > - COUENNE_INFINITY) && updateBound (-1, l + ind, l [wind] / c);
     41      resU = (u [wind] <   COUENNE_INFINITY) && updateBound ( 1, u + ind, u [wind] / c);
    4242    }
    4343    else if (c < - COUENNE_EPS) {
    44 
    45       //      return false;
    4644
    4745      //      printf ("w_%d [%g,%g] = %g x_%d [%g,%g]\n",
    4846      //              wind, l [wind], u [wind], c, ind, l [ind], u [ind]);
    4947
    50       res = (u [wind] <   COUENNE_INFINITY) && updateBound (-1, l + ind, u [wind] / c);
    51       res = (l [wind] > - COUENNE_INFINITY) && updateBound ( 1, u + ind, l [wind] / c) || res;
     48      resL = (u [wind] <   COUENNE_INFINITY) && updateBound (-1, l + ind, u [wind] / c);
     49      resU = (l [wind] > - COUENNE_INFINITY) && updateBound ( 1, u + ind, l [wind] / c);
    5250    }
    5351
    54     if (res) {
    55       chg [ind] = 1;
     52    if (resL) chg [ind].lower = CHANGED;
     53    if (resU) chg [ind].upper = CHANGED;
    5654
    5755      /*printf ("w_%d [%g,%g] -------> x_%d in [%g,%g] ",
    5856              wind, l [wind], u [wind],
    5957              ind,  l [ind],  u [ind]);*/
    60     }
    61 
    6258  } else {
    6359
     
    7773    // w's lower bound
    7874
    79     bool resx = false, resy = false;
     75    bool resxL, resxU, resyL, resyU =
     76      resxL = resxU = resyL = false;
    8077
    8178    if (wl >= 0.) {
     
    8481
    8582      if (*xu * *yu < wl) {
    86         resx = (*xu * *yl < wl) && updateBound (+1, xu, wl / *yl) || resx;
    87         resy = (*xl * *yu < wl) && updateBound (+1, yu, wl / *xl) || resy;
     83        resxU = (*xu * *yl < wl) && updateBound (+1, xu, wl / *yl);
     84        resyU = (*xl * *yu < wl) && updateBound (+1, yu, wl / *xl);
    8885      }
    8986
     
    9188
    9289      if (*xl * *yl < wl) {
    93         resx = (*xl * *yu < wl) && updateBound (-1, xl, wl / *yu) || resx;
    94         resy = (*xu * *yl < wl) && updateBound (-1, yl, wl / *xu) || resy;
     90        resxL = (*xl * *yu < wl) && updateBound (-1, xl, wl / *yu);
     91        resyL = (*xu * *yl < wl) && updateBound (-1, yl, wl / *xu);
    9592      }
    9693    } else {
     
    9996
    10097      // upper left
    101       resx = (*xl * *yl < wl) && (*yl > 0.) && updateBound (-1, xl, wl / *yl) || resx; // point C
    102       resy = (*xu * *yu < wl) && (*yu > 0.) && updateBound (+1, yu, wl / *xu) || resy; // point B
     98      resxL = (*xl * *yl < wl) && (*yl > 0.) && updateBound (-1, xl, wl / *yl); // point C
     99      resyU = (*xu * *yu < wl) && (*yu > 0.) && updateBound (+1, yu, wl / *xu); // point B
    103100
    104101      // lower right
    105       resy = (*xl * *yl < wl) && (*yl < 0.) && updateBound (-1, yl, wl / *xl) || resy; // point C
    106       resx = (*xu * *yu < wl) && (*yu < 0.) && updateBound (+1, xu, wl / *yu) || resx; // point B
     102      resyL = (*xl * *yl < wl) && (*yl < 0.) && updateBound (-1, yl, wl / *xl); // point C
     103      resxU = (*xu * *yu < wl) && (*yu < 0.) && updateBound (+1, xu, wl / *yu); // point B
    107104    }
    108105
     
    114111
    115112      // upper right
    116       resx = (*xu * *yl > wu) && (*yl > 0.) && updateBound (+1, xu, wu / *yl) || resx; // point D
    117       resy = (*xl * *yu > wu) && (*yu > 0.) && updateBound (+1, yu, wu / *xl) || resy; // point A
     113      resxU = (*xu * *yl > wu) && (*yl > 0.) && updateBound (+1, xu, wu / *yl) || resxU; // point D
     114      resyU = (*xl * *yu > wu) && (*yu > 0.) && updateBound (+1, yu, wu / *xl) || resyU; // point A
    118115
    119116      // lower left
    120       resx = (*xl * *yu > wu) && (*yu < 0.) && updateBound (-1, xl, wu / *yu) || resx; // point A
    121       resy = (*xu * *yl > wu) && (*yl < 0.) && updateBound (-1, yl, wu / *xu) || resy; // point D
     117      resxL = (*xl * *yu > wu) && (*yu < 0.) && updateBound (-1, xl, wu / *yu) || resxL; // point A
     118      resyL = (*xu * *yl > wu) && (*yl < 0.) && updateBound (-1, yl, wu / *xu) || resyL; // point D
    122119
    123120    } else {
     
    126123
    127124      if (*xu * *yl > wu) {
    128         resx = (*xu * *yu > wu) && updateBound (+1, xu, wu / *yu) || resx;
    129         resy = (*xl * *yl > wu) && updateBound (-1, yl, wu / *xl) || resy;
     125        resxU = (*xu * *yu > wu) && updateBound (+1, xu, wu / *yu) || resxU;
     126        resyL = (*xl * *yl > wu) && updateBound (-1, yl, wu / *xl) || resyL;
    130127      }
    131128
     
    133130
    134131      if (*xl * *yu > wu) {
    135         resx = (*xl * *yl > wu) && updateBound (-1, xl, wu / *yl) || resx;
    136         resy = (*xu * *yu > wu) && updateBound (+1, yu, wu / *xu) || resy;
     132        resxL = (*xl * *yl > wu) && updateBound (-1, xl, wu / *yl) || resxL;
     133        resyU = (*xu * *yu > wu) && updateBound (+1, yu, wu / *xu) || resyU;
    137134      }
    138135    }
     
    143140              else printf ("                                                 \r");*/
    144141
    145     if (resx) chg [xi] = 1;
    146     if (resy) chg [yi] = 1;
     142    if (resxL) chg [xi].lower = CHANGED;
     143    if (resxU) chg [xi].upper = CHANGED;
     144    if (resyL) chg [yi].lower = CHANGED;
     145    if (resyU) chg [yi].upper = CHANGED;
    147146
    148     res = resx || resy;
     147    resL = resxL || resyL;
     148    resU = resxU || resyU;
    149149  }
    150150
    151   return res;
     151  return (resL || resU);
    152152}
  • branches/Couenne/Couenne/src/expression/operators/impliedBounds-exprSum.cpp

    r525 r557  
    2121#define MALLOC_STEP 1000
    2222
    23 bool exprSum::impliedBound (int wind, CouNumber *l, CouNumber *u, char *chg) {
     23bool exprSum::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
    2424
    2525  /**
     
    200200      int ind = I1 [i];
    201201      if (tighter = updateBound (+1, u + ind, (wu - lower) / C1 [i] + lc [ind]) || tighter)
    202         chg [ind] = 1;
     202        chg [ind].upper = CHANGED;
    203203    }
    204204
     
    207207      int ind = I2 [i];
    208208      if (tighter = updateBound (-1, l + ind, (wu - lower) / C2 [i] + uc [ind]) || tighter)
    209         chg [ind] = 1;
     209        chg [ind].lower = CHANGED;
    210210    }
    211211  } else
     
    214214      int ind = I1 [infLo1];
    215215      if (tighter = updateBound (+1, u + ind, (wu - lower) / C1 [infLo1]) || tighter)
    216         chg [ind] = 1;
     216        chg [ind].upper = CHANGED;
    217217    }
    218218    else
     
    220220        int ind = I2 [infUp2];
    221221        if (tighter = updateBound (-1, l + ind, (wu - lower) / C2 [infUp2]) || tighter)
    222           chg [ind] = 1;
     222          chg [ind].lower = CHANGED;
    223223      }
    224224
     
    230230      int ind = I1 [i];
    231231      if (tighter = updateBound (-1, l + ind, (wl - upper) / C1 [i] + uc [ind]) || tighter)
    232         chg [ind] = 1;
     232        chg [ind].lower = CHANGED;
    233233    }
    234234
     
    236236      int ind = I2 [i];
    237237      if (tighter = updateBound (+1, u + ind, (wl - upper) / C2 [i] + lc [ind]) || tighter)
    238         chg [ind] = 1;
     238        chg [ind].upper = CHANGED;
    239239    }
    240240  } else
     
    243243      int ind = I1 [infUp1];
    244244      if (tighter = updateBound (-1, l + ind, (wl - upper) / C1 [infUp1]) || tighter)
    245         chg [ind] = 1;
     245        chg [ind].lower = CHANGED;
    246246    }
    247247    else
     
    249249        int ind = I2 [infLo2];
    250250        if (tighter = updateBound (+1, u + ind, (wl - upper) / C2 [infLo2]) || tighter)
    251           chg [ind] = 1;
     251          chg [ind].upper = CHANGED;
    252252      }
    253253
  • branches/Couenne/Couenne/src/problem/CouenneProblem.h

    r545 r557  
    145145
    146146  /// bound tightening
    147   int tightenBounds (char *) const;
     147  int tightenBounds (t_chg_bounds *) const;
    148148
    149149  /// search for implied bounds
    150   int impliedBounds (char *) const;
     150  int impliedBounds (t_chg_bounds *) const;
    151151};
    152152
  • branches/Couenne/Couenne/src/problem/impliedBounds.cpp

    r556 r557  
    1212/// Bound tightening for auxiliary variables
    1313
    14 int CouenneProblem::impliedBounds (char *chg_bds) const {
     14int CouenneProblem::impliedBounds (t_chg_bounds *chg_bds) const {
    1515
    1616  int nchg = 0, //< number of bounds changed for propagation
  • branches/Couenne/Couenne/src/problem/tightenBounds.cpp

    r553 r557  
    1515/// Bound propagation for auxiliary variables
    1616
    17 int CouenneProblem::tightenBounds (char *chg_bds) const {
     17int CouenneProblem::tightenBounds (t_chg_bounds *chg_bds) const {
    1818
    1919  int nchg = 0; //< number of bounds changed for propagation
     
    6767      }
    6868
    69       bool chg = false;
     69      //      bool chg = false;
    7070
    7171      // check if lower bound got higher   
     
    8484
    8585        lb_ [i+j] = ll;
    86         chg = true;
     86        chg_bds [i+j].lower = CHANGED;
     87        nchg++;
    8788      }
    8889
     
    101102
    102103        ub_ [i+j] = uu;
    103         chg = true;
    104       }
    105 
    106       if (chg && chg_bds && !(chg_bds [i+j])) {
     104        chg_bds [i+j].upper = CHANGED;
    107105        nchg++;
    108         chg_bds [i+j] = 1;
    109106      }
    110107
Note: See TracChangeset for help on using the changeset viewer.