Changeset 525


Ignore:
Timestamp:
Apr 30, 2007 1:48:15 PM (12 years ago)
Author:
pbelotti
Message:

implied bounds: fix some of the infinity handling

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

Legend:

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

    r524 r525  
    170170    // propagate bounds to auxiliary variables
    171171
    172     //    printf ("===== PROPAGATE\n");
    173 
    174172    ntightened = problem_ -> tightenBounds (chg_bds);
    175173
     
    177175    // expression may have non-propagated bounds
    178176
    179     //    printf ("===== IMPLIED\n");
    180 
    181177    if (ntightened >= 0) // if last call didn't signal infeasibility
    182178      nbwtightened = problem_ -> impliedBounds (chg_bds);
    183179
    184     //    printf ("============\n");
    185180    if ((ntightened < 0) || (nbwtightened < 0)) {
    186181
  • branches/Couenne/Couenne/src/convex/operators/conv-exprExp.cpp

    r520 r525  
    3636
    3737  if ((   u < log (COUENNE_INFINITY) - 1)
    38       && (l > -    COUENNE_INFINITY  + 1)) {
     38      && (l > -    COUENNE_INFINITY/1e4  + 1)) { // tame lower bound
    3939
    4040    CouNumber expl     = exp (l),
  • branches/Couenne/Couenne/src/convex/operators/conv-exprLog.cpp

    r520 r525  
    1818
    1919#define LOG_STEP 10
    20 #define LOG_MININF 1e-300
     20#define LOG_MININF 1e-50
    2121
    2222// generate convexification cut for constraint w = this
     
    5656  else if (x > u) x = u;
    5757
    58   if (u > COUENNE_INFINITY - 1)
     58  if (u > 1e5 * log (COUENNE_INFINITY) - 1)
    5959    u = x + (LOG_STEP << cg -> nSamples ());
    6060
  • branches/Couenne/Couenne/src/expression/exprAux.cpp

    r524 r525  
    5050
    5151  //  if (!(cg -> isFirst ()))
    52   //if (j < cs.sizeRowCuts ())
    53   if (0)
     52  if (j < cs.sizeRowCuts ())
     53    if (0)
    5454    {
    5555      printf ("----------------Generated cut for ");
  • branches/Couenne/Couenne/src/expression/operators/exprLog.cpp

    r450 r525  
    2525
    2626  // [low|upp]er bound of w=log(x) is log (max (0, [low|upp]er (x)))
    27   lb = new exprLog (new exprMax (new exprConst (1e-100), lba));
    28   ub = new exprLog (new exprMax (new exprConst (1e-100), uba));
     27  //  lb = new exprLog (new exprMax (new exprConst (1e-100), lba));
     28  //  ub = new exprLog (new exprMax (new exprConst (1e-100), uba));
     29
     30  expression **all  = new expression * [4];
     31
     32  all [0] = new exprClone (lba); all [1] = new exprLog (lba);
     33  all [2] = new exprConst (0);   all [3] = new exprConst (- COUENNE_INFINITY);
     34  lb = new exprMax (all, 4);
     35
     36  expression **alu  = new expression * [4],
     37             **alum = new expression * [4];
     38
     39  alum [0] = new exprConst (COUENNE_INFINITY);
     40  alum [1] = new exprConst (COUENNE_INFINITY);
     41  alum [2] = new exprClone (uba);
     42  alum [3] = new exprLog (uba);
     43
     44  alu [0] = new exprClone (uba); alu [1] = new exprMin (alum, 4);
     45  alu [2] = new exprConst (0);   alu [3] = new exprConst (- COUENNE_INFINITY);
     46  ub = new exprMax (alu, 4);
    2947}
    3048
     
    5977  res      = updateBound ( 1, u + ind, exp (u [wind])) || res;
    6078
    61   if (res)
     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]);*/
    6284    chg [ind] = 1;
     85  }
    6386
    6487  return res;
  • branches/Couenne/Couenne/src/expression/operators/exprPow.cpp

    r520 r525  
    203203    if (k > 0.) { // simple, just follow bounds
    204204
    205       res = updateBound (-1, l + index, pow (wl, 1./k));
    206       res = updateBound (+1, u + index, pow (wu, 1./k)) || res;
     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);
     212
    207213    } else // slightly more complicated, resort to same method as in exprInv
    208214      res = invPowImplBounds (wind, index, l, u, 1./k);
     
    218224      if (bound > COUENNE_EPS) {
    219225
    220         res = updateBound (-1, l + index, - pow (bound, 1./k));
    221         res = updateBound (+1, u + index,   pow (bound, 1./k)) || res;
     226        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;
     229        } else {
     230          res = updateBound (-1, l + index, - COUENNE_INFINITY);
     231          res = updateBound (+1, u + index,   COUENNE_INFINITY) || res;
     232        }
    222233      }
    223234    } else { // x^k, k=(1/h), h integer and even, or x^k, neither k nor 1/k integer
    224235
    225         CouNumber lb = wl, ub = wu;
    226 
    227         if (k < 0) { // swap bounds as they swap on the curve x^k when
    228           lb = wu;
    229           ub = wl;
    230         }
     236      CouNumber lb = wl, ub = wu;
     237
     238      if (k < 0) { // swap bounds as they swap on the curve x^k when
     239        lb = wu;
     240        ub = wl;
     241      }
    231242       
    232         if (lb > COUENNE_EPS) res = updateBound (-1, l + index, pow (lb, 1./k));
     243      if (lb > COUENNE_EPS) res = updateBound (-1, l + index, pow (lb, 1./k));
     244
     245      if (fabs (ub) < COUENNE_INFINITY) {
    233246        if (ub > COUENNE_EPS) res = updateBound (+1, u + index, pow (ub, 1./k)) || res;
    234       }
     247      } else                  res = updateBound (+1, u + index, COUENNE_INFINITY) || res;       
     248    }
    235249
    236250  return res;
  • branches/Couenne/Couenne/src/expression/operators/exprSub.cpp

    r400 r525  
    112112  bool res = false;
    113113
     114  return false;
     115
    114116  // w >= l
    115117
  • branches/Couenne/Couenne/src/expression/operators/impliedBounds-exprMul.cpp

    r450 r525  
    2222      (arglist_ [ind=1] -> Type () <= CONST)) {
    2323
     24    // at least one constant in product w=cx:
     25    //
     26    // wl/c <= x <= wu/c, if c is positive
     27    // wu/c <= x <= wl/c, if c is negative
     28
    2429    CouNumber c = arglist_ [ind] -> Value ();
    2530
     
    3237    if (c > COUENNE_EPS) {
    3338
    34       res = updateBound (-1, l + ind, l [wind] / c);
    35       res = updateBound ( 1, u + ind, u [wind] / c) || res;
     39     
     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;
    3642    }
    3743    else if (c < - COUENNE_EPS) {
    3844
    39       res = updateBound (-1, l + ind, u [wind] / c);
    40       res = updateBound ( 1, u + ind, l [wind] / c) || res;
     45      //      return false;
     46
     47      //      printf ("w_%d [%g,%g] = %g x_%d [%g,%g]\n",
     48      //              wind, l [wind], u [wind], c, ind, l [ind], u [ind]);
     49
     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;
    4152    }
    4253
    43     if (res)
     54    if (res) {
    4455      chg [ind] = 1;
     56
     57      /*printf ("w_%d [%g,%g] -------> x_%d in [%g,%g] ",
     58              wind, l [wind], u [wind],
     59              ind,  l [ind],  u [ind]);*/
     60    }
    4561
    4662  } else {
  • branches/Couenne/Couenne/src/expression/operators/impliedBounds-exprSum.cpp

    r524 r525  
    110110  // now we have correct  I1, I2, C1, C2, ipos, ineg, and a0
    111111
    112   if (0) {
     112  /*if (0) {
    113113
    114114    printf ("w_%d = ", wind); print (std::cout);
     
    126126
    127127    printf ("\n");
    128   }
     128    }*/
    129129
    130130  // indices of the variable in I1 or I2 with infinite lower or upper
     
    252252      }
    253253
    254   if (0) {
     254  /*if (0) {
    255255
    256256    printf ("I1 (%d): ", ipos);
     
    263263
    264264    printf ("\n---------------------------------------------------------------\n");
    265   }
     265    }*/
    266266
    267267  // ...phew!
  • branches/Couenne/Couenne/src/include/operators/exprBDiv.h

    r524 r525  
    1313
    1414
    15 /// division that avoids NaN's
    16 inline CouNumber safeDiv (register CouNumber a, register CouNumber b) {
     15/// division that avoids NaN's and considers a sign when returning infinity
     16static inline CouNumber safeDiv (register CouNumber a, register CouNumber b, int sign) {
    1717
    18   if ((fabs (a) < 1e-10) && (fabs (b) < 1e-10)) return 0;
     18  if (fabs (a) < COUENNE_EPS) return 0;
     19  //    if (fabs (b) < COUENNE_EPS)) return 0;
     20  //    else return 0
    1921
    20   if (a >  1e20) return ((b < 0) ? -COUENNE_INFINITY :  COUENNE_INFINITY);
    21   if (a < -1e20) return ((b < 0) ?  COUENNE_INFINITY : -COUENNE_INFINITY);
     22  if (fabs (b) < COUENNE_EPS) return ((sign < 0) ? -COUENNE_INFINITY :  COUENNE_INFINITY);
     23
     24  if (a >  COUENNE_INFINITY) return ((sign < 0) ? -COUENNE_INFINITY :  COUENNE_INFINITY);
     25  if (a < -COUENNE_INFINITY) return ((sign < 0) ? -COUENNE_INFINITY :  COUENNE_INFINITY);
    2226
    2327  return a/b;
     
    6165  register CouNumber d = (*(arglist_ [2])) ();
    6266  register CouNumber D = (*(arglist_ [3])) ();
    63                                                      // (n,N,d,D)     lb
    64   if (d > COUENNE_EPS)                               // (?,?,+,+)
    65     if   (n > 0) return safeDiv (n,D);               // (+,+,+,+) --> n/D
    66     else         return safeDiv (n,d);               // (-,?,+,+) --> n/d
     67                                               // (n,N,d,D)     lb
     68  if (d > 0)                                   // (?,?,+,+)
     69    if   (n > 0)    return safeDiv (n,D,-1);      // (+,+,+,+) --> n/D
     70    else            return safeDiv (n,d,-1);      // (-,?,+,+) --> n/d
    6771  else { // d <= 0
    68     if      (D > COUENNE_EPS) return - COUENNE_INFINITY; // (?,?,-,+) --> unbounded
    69     else if (N > COUENNE_EPS) return safeDiv (N,D);      // (?,+,-,-) --> N/D
    70     else                      return safeDiv (N,d);      // (-,-,-,-) --> N/d
     72    if      (D > 0) return - COUENNE_INFINITY; // (?,?,-,+) --> unbounded
     73    else if (N > 0) return safeDiv (N,D,-1);      // (?,+,-,-) --> N/D
     74    else            return safeDiv (N,d,-1);      // (-,-,-,-) --> N/d
    7175  }
    7276}
     
    99103inline CouNumber exprUBDiv::operator () () {
    100104
    101   //  exprOp:: operator () ();
    102 
    103105  register CouNumber n = (*(arglist_ [0])) ();
    104106  register CouNumber N = (*(arglist_ [1])) ();
    105107  register CouNumber d = (*(arglist_ [2])) ();
    106108  register CouNumber D = (*(arglist_ [3])) ();
    107                                                        // (n,N,d,D)     lb
    108   if (d > COUENNE_EPS)
    109     if   (N < 0) return safeDiv (N,D);                 // (-,-,+,+) --> N/D
    110     else         return safeDiv (N,d);                 // (?,+,+,+) --> N/d
     109
     110  if (d > 0)                                   // (n,N,d,D)     lb
     111    if   (N < 0) return safeDiv (N,D,1);         // (-,-,+,+) --> N/D
     112    else         return safeDiv (N,d,1);         // (?,+,+,+) --> N/d
    111113  else { // d <= 0
    112     if      (D >   COUENNE_EPS) return + COUENNE_INFINITY; // (?,?,-,+) --> unbounded
    113     else if (n < - COUENNE_EPS) return safeDiv (n,D);  // (-,?,-,-) --> n/D
    114     else                        return safeDiv (n,d);  // (+,+,-,-) --> n/d
     114    if      (D > 0) return + COUENNE_INFINITY; // (?,?,-,+) --> unbounded
     115    else if (n < 0) return safeDiv (n,D,1);      // (-,?,-,-) --> n/D
     116    else            return safeDiv (n,d,1);      // (+,+,-,-) --> n/d
    115117  }
    116118}
  • branches/Couenne/Couenne/src/include/operators/exprPow.h

    r520 r525  
    8080                           register CouNumber exponent) {
    8181
    82   register int rndexp;
     82  if (base < 0) {
    8383
    84   if ((base < 0) &&
    85       ((fabs (exponent - (rndexp = COUENNE_round (exponent))) < COUENNE_EPS) ||
    86        ((fabs (exponent) > COUENNE_EPS) &&
    87         (fabs (1. / exponent - (rndexp = COUENNE_round (1. / exponent))) < COUENNE_EPS)))
    88       && (rndexp % 2))
    89     return (- pow (- base, exponent));
     84    register int rndexp;
     85
     86    if (((fabs (exponent - (rndexp = COUENNE_round (exponent))) < COUENNE_EPS) ||
     87         ((fabs (exponent) > COUENNE_EPS) &&
     88          (fabs (1. / exponent - (rndexp = COUENNE_round (1. / exponent))) < COUENNE_EPS)))
     89        && (rndexp % 2))
     90      return (- pow (- base, exponent));
     91  }
     92
     93  if (fabs (base) > COUENNE_INFINITY - 1) {
     94
     95    if (base < -COUENNE_INFINITY+1) {
     96
     97      register int intk = COUENNE_round (exponent);
     98
     99      if ((fabs (exponent - intk) < COUENNE_EPS) && (intk % 2))
     100        return -COUENNE_INFINITY;
     101    }
     102    else return COUENNE_INFINITY;
     103  }
    90104
    91105  return (pow (base, exponent));
  • branches/Couenne/Couenne/src/problem/impliedBounds.cpp

    r524 r525  
    1717      nvar = nVars ();
    1818
    19   /*for (int i=0; i < nVars () + nAuxs (); i++)
    20     printf ("x%3d: [%12.4f,%12.4f]\n", i, lb_ [i], ub_ [i]);*/
     19  //printf ("implied========================\n");
     20
     21  /*  for (int i=0; i < nVars () + nAuxs (); i++)
     22      printf ("x%d: [%g,%g]\n", i, lb_ [i], ub_ [i]);*/
    2123
    2224  for (int i=nAuxs (); i--;) {
    2325
    24     if (lb_ [nvar+i] > ub_ [nvar+i] + COUENNE_EPS)
     26    if (lb_ [nvar+i] > ub_ [nvar+i] + COUENNE_EPS) {
     27      /*printf ("w_%d has infeasible bounds [%g,%g]\n",
     28        i+nvar, lb_ [nvar+i], ub_ [nvar+i]);*/
    2529      return -1;
     30    }
    2631
    2732    //    if ((auxiliaries_ [i] -> Image () -> code () == COU_EXPRSUM) ||
    2833    //  (auxiliaries_ [i] -> Image () -> code () == COU_EXPRGROUP))
    29     if (auxiliaries_ [i] -> Image () -> impliedBound (nvar+i, lb_, ub_, chg_bds) > COUENNE_EPS)
     34
     35    CouNumber l0 = lb_ [nvar+i],
     36              u0 = ub_ [nvar+i];
     37
     38    if (auxiliaries_ [i] -> Image () -> impliedBound (nvar+i, lb_, ub_, chg_bds) > COUENNE_EPS) {
     39      /*printf ("implied: [%g,%g] -> [%g,%g] ", l0, u0, lb_ [nvar+i], ub_ [nvar+i]);
     40      auxiliaries_ [i] -> print (std::cout); printf (" := ");
     41      auxiliaries_ [i] -> Image () -> print (std::cout); printf ("\n");*/
    3042      nchg++;
     43    }
    3144  }
    3245
  • branches/Couenne/Couenne/src/problem/tightenBounds.cpp

    r524 r525  
    2828  // they depend on
    2929
    30   /*for (int i=0; i<nVars(); i++)
    31     printf (" [%g, %g]\n",
     30  /*printf ("tighten========================\n");
     31
     32  for (int i=0; i<nVars(); i++)
     33    printf ("x_%d [%g, %g]\n", i,
    3234            expression::Lbound (i),
    3335            expression::Ubound (i));*/
     
    3638       j < naux; j++) {
    3739
    38     /*printf (" [%g, %g] ",
     40    /*    printf ("w_%d [%g, %g] ", i+j,
    3941            expression::Lbound (i+j),
    4042            expression::Ubound (i+j));*/
     
    4345              uu = (*(Aux (j) -> Ub ())) ();
    4446
    45     /*printf (" ---> [%g, %g] ",
    46             expression::Lbound (i+j),
    47             expression::Ubound (i+j));*/
     47    //    printf (" ---> [%g, %g]\n ", ll, uu);
    4848
    4949    /*auxiliaries_ [j] -> print (std::cout);
     
    5252    printf ("\n");*/
    5353
    54     if (ll > uu + COUENNE_EPS)
     54    if (ll > uu + COUENNE_EPS) {
     55      /*printf ("w_%d has infeasible bounds [%g,%g]: ", i+j, ll, uu);
     56      Aux (j) -> Lb () -> print (std::cout); printf (" - ");
     57      Aux (j) -> Ub () -> print (std::cout); printf ("\n");*/
     58
    5559      return -1; // declare this node infeasible
     60    }
    5661
    5762    bool chg = false;
     
    6671        i+j, ll, lb_ [i+j]);*/
    6772
     73      /*printf (" propaga: [%g,(%g)] -> [%g,(%g)] ", lb_ [i+j], ub_ [i+j], ll, uu);
     74      Aux (j)             -> print (std::cout); printf (" := ");
     75      Aux (j) -> Image () -> print (std::cout); printf ("\n");*/
     76
    6877      lb_ [i+j] = ll;
    6978      chg = true;
     
    7786      /*printf ("update ubound %d: %g >= %g\n",
    7887        i+j, uu, ub_ [i+j]);*/
     88
     89      /*printf (" propaga: [(%g),%g] -> [(%g),%g] ", lb_ [i+j], ub_ [i+j], ll, uu);
     90      Aux (j)             -> print (std::cout); printf (" := ");
     91      Aux (j) -> Image () -> print (std::cout); printf ("\n");*/
    7992
    8093      ub_ [i+j] = uu;
Note: See TracChangeset for help on using the changeset viewer.