Changeset 520


Ignore:
Timestamp:
Apr 27, 2007 3:39:58 PM (13 years ago)
Author:
pbelotti
Message:

limit power upper/lower bounds to get reasonable coefficients

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

Legend:

Unmodified
Added
Removed
  • branches/Couenne/Couenne/src/convex/addEnvelope.cpp

    r506 r520  
    2020                                       CouNumber x, CouNumber l, CouNumber u,
    2121                                       bool is_global) const {
    22 
    2322  CouNumber opp_slope = - fprime (x);
    2423
  • branches/Couenne/Couenne/src/convex/generateCuts.cpp

    r515 r520  
    231231    for (int i=0, j = problem_ -> nAuxs (); j--; i++) {
    232232
    233       expression * image = problem_ -> Aux (i) -> Image ();
     233      // TODO: check if list contains all and only aux's to cut
     234      /*expression * image = problem_ -> Aux (i) -> Image ();
    234235      if (   (image -> Linearity () > LINEAR)          // if linear, no need to cut twice
    235236          && (image -> dependsOn (changed, nchanged))  // if expression does not depend on
    236           )                                            // changed variables, do not cut
     237          )*/                                            // changed variables, do not cut
    237238        problem_ -> Aux (i) -> generateCuts (si, cs, this);
    238239    }
  • branches/Couenne/Couenne/src/convex/operators/conv-exprExp.cpp

    r480 r520  
    1111#include <exprConst.h>
    1212#include <exprAux.h>
     13#include <exprPow.h>
    1314
    1415#include <CouenneProblem.h>
     
    2425  argument_ -> getBounds (le, ue);
    2526
    26   CouNumber x = (*argument_) (),
     27  CouNumber x = (cg -> isFirst ()) ?
     28                 0 : powNewton ((*argument_) (), (*aux) (), exp, exp, exp),
    2729            l = (*le) (),
    2830            u = (*ue) ();
     
    4850  int ns = cg -> nSamples ();
    4951
    50   // change bounds to get finite coefficients
    51 
    52   //  CouNumber fact = 2 * ns;
    53 
    54   /*if (x > 0) {
    55     if (l < log (COUENNE_EPS))      l = x / fact - 1;
    56     if (u > log (COUENNE_INFINITY)) u = x * fact + 1;
    57   }
    58   else {
    59     if (l < log (COUENNE_EPS))      l = x * fact - 1;
    60     if (u > log (COUENNE_INFINITY)) u = x / fact + 1;
    61     }*/
    62 
    6352  // add tangents with finite coefficients
    6453  if (l < log (COUENNE_EPS))      l = x - ns;
  • branches/Couenne/Couenne/src/convex/operators/conv-exprGroup.cpp

    r460 r520  
    9393  OsiRowCut *cut   = new OsiRowCut;
    9494
    95   CouNumber rhs = c0_;
     95  CouNumber rhs = - c0_;
    9696
    9797  // first, make room for aux variable
     
    111111
    112112    if (curr -> Type () == CONST) // constant term in sum
    113       rhs += curr -> Value ();
     113      rhs -= curr -> Value ();
    114114    else {                        // variable
    115115      coeff [++nterms] = 1.;
     
    119119
    120120  cut -> setRow (nterms+1, index, coeff);
    121 
    122   rhs = - rhs;
    123121
    124122  cut -> setUb (rhs);
  • branches/Couenne/Couenne/src/convex/operators/conv-exprInv.cpp

    r476 r520  
    4646
    4747
    48 // derivative of inv (x)
    49 
    50 inline CouNumber oppInvSqr (register CouNumber x)
    51 {return (- inv (x*x));}
    52 
    53 
    5448#define MIN_DENOMINATOR 1e-10
    5549
     
    7973  aux -> getBounds (wle, wue);
    8074
    81   expression *xe = argument_;
    82 
    83   CouNumber x = (*xe)  ();
     75  CouNumber x;
    8476
    8577  int w_ind = aux       -> Index (),
     
    10799  }
    108100
    109   int sign = (l > 0) ? +1 : -1;
    110 
    111101  // bound
    112 
    113   cg -> addEnvelope (cs, sign, inv, oppInvSqr, w_ind, x_ind, x, l, u);
     102  cg -> addEnvelope
     103    (cs, (l > 0) ? +1 : -1,
     104     inv, oppInvSqr, w_ind, x_ind,
     105     (cg -> isFirst ()) ?
     106       // place it somewhere in the interval (we won't care)
     107       ((l > COUENNE_EPS) ? l : u) :
     108       // not first call, gotta replace it where it gives deepest cut
     109       powNewton ((*argument_) (), (*aux) (), inv, oppInvSqr, inv_dblprime),
     110     l, u);
    114111
    115112  delete xle;
  • branches/Couenne/Couenne/src/convex/operators/conv-exprLog.cpp

    r480 r520  
    99#include <CouenneTypes.h>
    1010#include <exprLog.h>
     11#include <exprInv.h>
     12#include <exprPow.h>
    1113#include <exprConst.h>
    1214
     
    2729  argument_ -> getBounds (le, ue);
    2830
    29   CouNumber x = (*argument_) (),
     31  CouNumber x = (cg -> isFirst ()) ?
     32                 1 : powNewton ((*argument_) (), (*aux) (), log, inv, oppInvSqr),
    3033            l = (*le) (),
    3134            u = (*ue) ();
  • branches/Couenne/Couenne/src/convex/operators/conv-exprPow-envelope.cpp

    r506 r520  
    5252  // orthogonal with line through current- and tangent point)
    5353
    54   //  x = powNewton (x, y, power_k, power_k_prime, power_k_dblprime);
     54  if (!(cg -> isFirst ()))
     55    x = powNewton (x, y, power_k, power_k_prime, power_k_dblprime);
    5556
    5657  if      (x<l) x=l;
     
    5960  // limit the bounds for the envelope
    6061
    61   CouNumber step = 1 + log (1. + (double) (cg -> nSamples ()));
     62  CouNumber step     = 1 + log (1. + (double) (cg -> nSamples ())),
     63            powThres = pow (COU_MAX_COEFF, 1./k);
    6264
    63   if (l < - COUENNE_INFINITY + 1) {
     65  if (l < - powThres + 1) {
    6466    l = x - step;
    65     if (u > COUENNE_INFINITY - 1)
     67    if (u > powThres - 1)
    6668      u = x + step;
    6769  } else
    68     if (u > COUENNE_INFINITY - 1)
     70    if (u > powThres - 1)
    6971      u = x + step;
    7072
  • branches/Couenne/Couenne/src/convex/operators/conv-exprPow-getBounds.cpp

    r476 r520  
    6262      int rndexp;
    6363
    64       bool isInt    =  fabs (expon - (rndexp = FELINE_round (expon))) < COUENNE_EPS;
     64      bool isInt    =  fabs (expon - (rndexp = COUENNE_round (expon))) < COUENNE_EPS;
    6565      bool isInvInt = !isInt && 
    6666                      ((fabs (expon) > COUENNE_EPS) &&
    67                        (fabs (1/expon - (rndexp = FELINE_round (1/expon))) < COUENNE_EPS));
     67                       (fabs (1/expon - (rndexp = COUENNE_round (1/expon))) < COUENNE_EPS));
    6868
    6969      if ((isInt || isInvInt) && (rndexp % 2) && (rndexp > 0)) {
  • branches/Couenne/Couenne/src/convex/operators/conv-exprPow.cpp

    r507 r520  
    7878  int x_ind = xe  -> Index ();
    7979
    80   CouNumber w = (*aux) (),
    81             x = (*xe)  (),
     80  CouNumber w, x,
    8281            l = (*xle) (),
    8382            u = (*xue) ();
     
    117116  }
    118117
    119   bool isInt    =            fabs (k    - (double) (intk = FELINE_round (k)))    < COUENNE_EPS,
    120        isInvInt = !isInt && (fabs (1./k - (double) (intk = FELINE_round (1./k))) < COUENNE_EPS);
     118  bool isInt    =            fabs (k    - (double) (intk = COUENNE_round (k)))    < COUENNE_EPS,
     119       isInvInt = !isInt && (fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS);
    121120
    122121  // two macro-cases:
     
    142141
    143142    if (isInvInt) {
     143      if (cg -> isFirst ()) {
     144        w = (l>0) ? 1 : (u<0) ? -1 : 0;
     145        x = 0;
     146      }
    144147      q = safe_pow (q, k);
    145148      sign = -1;
    146149    }
    147     else sign = 1;
     150    else {
     151      if (cg -> isFirst ()) {
     152        x = (l>0) ? l : (u<0) ? u : 0;
     153        w = 0;
     154      }
     155      sign = 1;
     156    }
     157
     158    CouNumber powThres = pow (COU_MAX_COEFF, 1./k); // don't want big coefficients
    148159
    149160    // lower envelope
    150161
    151     if (l > -COUENNE_INFINITY) {
     162    if (l > -powThres) {
    152163      if (u > q * l) { // upper x is after "turning point", add lower envelope
    153164        addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, q*l, u, sign);
     
    158169    // upper envelope
    159170
    160     if (u < COUENNE_INFINITY) {
     171    if (u < powThres) {
    161172      if (l < q * u) { // lower x is before "turning point", add upper envelope
    162173        addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, q*u, -sign);
     
    176187        && (!isInvInt || !(intk % 2))
    177188        && (l < - COUENNE_EPS)
    178         && (u < (l=0)))        // CAUTION! l is updated here, if negative
     189        && (u < (l=0)))        // CAUTION! l updated if negative (k real)
    179190      return;
    180191
     
    208219    // upper envelope -- when k negative, add only if
    209220
    210     CouNumber powThreshold = pow (1e20, 1./k);
     221    CouNumber powThres = pow (COU_MAX_COEFF, 1./k), // don't want big coefficients
     222              powStep  = 1;
    211223
    212224    if ((  (k > COUENNE_EPS)
    213         || (l > COUENNE_EPS)             // bounds do not contain 0
     225        || (l > COUENNE_EPS)       // bounds do not contain 0
    214226        || (u < - COUENNE_EPS)) &&
    215         (l > - powThreshold) &&  // and are finite
    216         (u <   powThreshold))
     227        (l > - powThres) &&    // and are finite
     228        (u <   powThres))
    217229
    218230      cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l, k), u, safe_pow (u, k), -sign);
     
    220232    // similarly, pay attention not to add infinite slopes
    221233
    222 #define POWER_RANGE 1e2;
     234    if (cg -> isFirst()) {
     235      x = (k<0) ? ((u<0) ? u : (l>0) ? l : 0) : 0;
     236      w = 0;
     237    }
    223238
    224239    if (k > COUENNE_EPS) {
    225240
    226       if (u >   powThreshold) u = x + POWER_RANGE;
    227       if (l < - powThreshold) l = x - POWER_RANGE;
     241      if (u >   powThres) u = x + powStep;
     242      if (l < - powThres) l = x - powStep;
    228243    }
    229244    else {
    230245
    231       if (fabs (l) < COUENNE_EPS) l =  1. / POWER_RANGE; // l --> 0+
    232       if (fabs (u) < COUENNE_EPS) u = -1. / POWER_RANGE; // u --> 0-
     246      if (fabs (l) < COUENNE_EPS) l =  1. / powThres; // l --> 0+
     247      if (fabs (u) < COUENNE_EPS) u = -1. / powThres; // u --> 0-
    233248    }
    234249
  • branches/Couenne/Couenne/src/convex/operators/conv-exprSub.cpp

    r476 r520  
    3737  expression *y = arglist_ [1];
    3838
    39   int wind = w -> Index ();
    40   int xind = x -> Index ();
    41   int yind = y -> Index ();
     39  int wi = w -> Index ();
     40  int xi = x -> Index ();
     41  int yi = y -> Index ();
    4242
    43   if (x->Type () == CONST) { // (c - y) or (c - d)
    44 
    45     CouNumber c = x -> Value ();
    46 
    47     if (y->Type() == CONST) cg -> createCut (cs, c - y->Value(), 0, wind, 1.,   -1, 0., -1, 0., true);
    48     else                    cg -> createCut (cs, c,              0, wind, 1., yind, 1., -1, 0., true);
    49   }
     43  if (x->Type () == CONST) // (c - y) or (c - d)
     44    if (y->Type() == CONST) cg -> createCut (cs, x->Value()-y->Value(), 0, wi, 1, -1, 0, -1, 0, true);
     45    else                    cg -> createCut (cs, x->Value(),            0, wi, 1, yi, 1, -1, 0, true);
    5046  else // (x - y) or (x - d)
    51     if (y->Type() == CONST) cg -> createCut (cs, y->Value(),     0, wind, -1., xind, 1., -1, 0.,true);
    52     else                    cg -> createCut (cs, 0.,          0, wind, -1., xind, 1., yind, -1.,true);
     47    if (y->Type() == CONST) cg -> createCut (cs, y->Value(),     0, wi, -1., xi, 1., -1, 0.,true);
     48    else                    cg -> createCut (cs, 0.,          0, wi, -1., xi, 1., yi, -1.,true);
    5349}
  • branches/Couenne/Couenne/src/convex/operators/conv-exprSum.cpp

    r450 r520  
    1515
    1616
    17 // generate convexification cut for constraint w = this
     17/// generate convexification cut for constraint w = this
    1818
    1919void exprSum::generateCuts (exprAux *w, const OsiSolverInterface &si,
     
    2929  CouNumber rhs = 0;
    3030
    31   // first, make room for aux variable
     31  /// first, make room for aux variable
    3232  coeff [0] = -1.; index [0] = w -> Index ();
    3333
    3434  int nv = 1;
    3535
    36   // scan arglist for (aux) variables and constants
     36  /// scan arglist for (aux) variables and constants
    3737  for (int i=0; i<nargs_; i++) {
    3838
    3939    if (arglist_ [i] -> Type () == CONST)
    40       rhs += arglist_ [i] -> Value ();
     40      rhs -= arglist_ [i] -> Value ();
    4141    else {
    4242      coeff [nv]   = 1.;
     
    4646
    4747  cut -> setRow (nv, index, coeff);
    48   cut -> setUb (-rhs);
    49   cut -> setLb (-rhs);
     48  cut -> setUb (rhs);
     49  cut -> setLb (rhs);
    5050
    5151  delete [] index;
    5252  delete [] coeff;
    5353
    54   // added only once, it is global
     54  /// added only once, it is global
    5555  cut -> setGloballyValid ();
    5656
  • branches/Couenne/Couenne/src/expression/exprAux.cpp

    r515 r520  
    5454      print (std::cout);  printf (" := ");
    5555      image_ -> print (std::cout);
     56
     57      printf (" [%.3e,%.3e] ---> ",
     58              expression::Lbound (varIndex_),
     59              expression::Ubound (varIndex_));
     60
     61      int index;
     62      if ((image_ -> Argument ()) &&
     63          ((index = image_ -> Argument () -> Index ()) >= 0))
     64        printf ("[%.3e,%.3e] ",
     65                expression::Lbound (index),
     66                expression::Ubound (index));
     67      else if (image_ -> ArgList ())
     68        for (int i=0; i<image_ -> nArgs (); i++)
     69          if ((index = image_ -> ArgList () [i] -> Index ()) >= 0)
     70        printf ("[%.3e,%.3e] ",
     71                expression::Lbound (index),
     72                expression::Ubound (index));
     73       
    5674      printf("\n");
    5775      for (int jj=j; jj < cs.sizeRowCuts ();jj++)
  • branches/Couenne/Couenne/src/expression/operators/exprPow.cpp

    r498 r520  
    123123    double exponent = arglist_ [1] -> Value ();
    124124
    125     if (fabs (exponent - FELINE_round (exponent)) > COUENNE_EPS)
     125    if (fabs (exponent - COUENNE_round (exponent)) > COUENNE_EPS)
    126126      return NONLINEAR;
    127127
    128128    if (arglist_ [1] -> Type () == CONST) {
    129129
    130       int expInt = (int) FELINE_round (exponent);
     130      int expInt = (int) COUENNE_round (exponent);
    131131
    132132      if (arglist_ [0] -> Linearity () == LINEAR) {
     
    192192  int intk; // integer (or integer inverse of) exponent
    193193
    194   bool isint    =           (k    - (intk = FELINE_round (k))    < COUENNE_EPS), // k   is integer
    195        isinvint = !isint && (1./k - (intk = FELINE_round (1./k)) < COUENNE_EPS); // 1/k is integer
     194  bool isint    =           (k    - (intk = COUENNE_round (k))    < COUENNE_EPS), // k   is integer
     195       isinvint = !isint && (1./k - (intk = COUENNE_round (1./k)) < COUENNE_EPS); // 1/k is integer
    196196
    197197  CouNumber wl = l [wind], // lower w
  • branches/Couenne/Couenne/src/include/CouennePrecisions.h

    r427 r520  
    1313#include <math.h>
    1414
    15 #define COUENNE_EPS       1e-7 /* keep it at least 1e-7 or strange things happen */
     15/* keep it at least 1e-7 or strange things happen */
     16#define COUENNE_EPS       1e-7
    1617#define COUENNE_EPS_SIMPL 1e-20
    1718#define COUENNE_INFINITY  1e+50
     19#define COU_MAX_COEFF     1e6
    1820
    19 #define FELINE_round(x) ((int) (floor ((x)+0.5)))
     21#define COUENNE_round(x) ((int) (floor ((x)+0.5)))
    2022
    2123#endif
  • branches/Couenne/Couenne/src/include/operators/exprInv.h

    r395 r520  
    1313
    1414
    15 // the operator itself
    16 
     15/// the operator itself
    1716inline CouNumber inv (register CouNumber arg)
    1817{return 1.0 / arg;}
    1918
     19/// derivative of inv (x)
     20inline CouNumber oppInvSqr (register CouNumber x)
     21{return (- inv (x*x));}
    2022
    21 // class inverse (1/f(x))
     23
     24/// inv_dblprime, second derivative of inv (x)
     25inline CouNumber inv_dblprime (register CouNumber x)
     26{return (2 * inv (x*x*x));}
     27
     28
     29/// class inverse (1/f(x))
    2230
    2331class exprInv: public exprUnary {
     
    2533 public:
    2634
    27   // Constructors, destructor
     35  /// Constructors, destructor
    2836  exprInv  (expression *al):
    2937    exprUnary (al) {} //< non-leaf expression, with argument list
    3038
    31   // cloning method
     39  /// cloning method
    3240  expression *clone () const
    3341    {return new exprInv (argument_ -> clone ());}
     
    3644  inline unary_function F () {return inv;}
    3745
    38   // output "1/argument"
     46  /// output "1/argument"
    3947  void print (std::ostream&) const;
    4048
    41   // differentiation
     49  /// differentiation
    4250  expression *differentiate (int index);
    4351
    44   // get a measure of "how linear" the expression is (see CouenneTypes.h)
     52  /// get a measure of "how linear" the expression is (see CouenneTypes.h)
    4553  virtual inline int Linearity () {
    4654    if (argument_ -> Type () == CONST) return CONSTANT;
     
    4856  }
    4957
    50   // Get lower and upper bound of an expression (if any)
     58  /// Get lower and upper bound of an expression (if any)
    5159  void getBounds (expression *&, expression *&);
    5260
    53   // generate equality between *this and *w
     61  /// generate equality between *this and *w
    5462  void generateCuts (exprAux *w, const OsiSolverInterface &si,
    5563                     OsiCuts &cs, const CouenneCutGenerator *cg);
    5664
    57   ///
     65  /// code for comparisons
    5866  virtual enum expr_type code () {return COU_EXPRINV;}
    5967
  • branches/Couenne/Couenne/src/include/operators/exprPow.h

    r506 r520  
    8383
    8484  if ((base < 0) &&
    85       ((fabs (exponent - (rndexp = FELINE_round (exponent))) < COUENNE_EPS) ||
     85      ((fabs (exponent - (rndexp = COUENNE_round (exponent))) < COUENNE_EPS) ||
    8686       ((fabs (exponent) > COUENNE_EPS) &&
    87         (fabs (1. / exponent - (rndexp = FELINE_round (1. / exponent))) < COUENNE_EPS)))
     87        (fabs (1. / exponent - (rndexp = COUENNE_round (1. / exponent))) < COUENNE_EPS)))
    8888      && (rndexp % 2))
    8989    return (- pow (- base, exponent));
  • branches/Couenne/Couenne/src/readnl/readnl.cpp

    r448 r520  
    116116
    117117    expression *subst = body -> simplify ();
    118     if (subst) body = subst;
     118    if (subst) {
     119      delete body; // VALGRIND
     120      body = subst;
     121    }
    119122
    120123    // ThirdParty/ASL/solvers/asl.h, line 336: 0 is minimization, 1 is maximization
     
    248251
    249252    expression *subst = body -> simplify ();
    250     if (subst) body = subst;
     253    if (subst) {
     254      delete body; // VALGRIND
     255      body = subst;
     256    }
    251257
    252258    // add them (and set lower-upper bound)
Note: See TracChangeset for help on using the changeset viewer.