Changeset 506


Ignore:
Timestamp:
Apr 25, 2007 6:09:37 PM (12 years ago)
Author:
pbelotti
Message:

avoid cuts with huge coefficients. Added powNewton to compute deepest cut for powers

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

Legend:

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

    r476 r506  
    2525  // TODO: remove check of !firstcall_ if point is available already
    2626
     27  // Add tangent in any case
     28
    2729  if (((!firstcall_) || ((x >= l) && (x <= u)))
    2830      && (fabs (opp_slope) < COUENNE_INFINITY))
     
    4648      opp_slope = - fprime (sample);
    4749
    48       if (fabs (opp_slope) < COUENNE_INFINITY)
     50      if ((fabs (opp_slope) < COUENNE_INFINITY) &&
     51          (fabs (sample-x) > COUENNE_EPS)) // do not add twice cut at current point
    4952          createCut (cs, f (sample) + opp_slope * sample, sign,
    5053                     w_ind, 1.,
    5154                     x_ind, opp_slope,
    52                      -1, 0.,
    53                      is_global);
     55                     -1, 0., is_global);
    5456
    5557        //      printf ("  Uniform %d: ", i); cut -> print ();
     
    5860    }
    5961  }
    60   /* else if (convtype_ == CURRENT_ONLY) {
    6162
    62     if (fabs (opp_slope) < COUENNE_INFINITY)
    63       createCut (cs, f (x) + opp_slope * x, sign, w_ind, 1.,
    64                  x_ind, opp_slope, -1, 0., is_global);
    65 
    66       //      addTangent (cs, w_ind, x_ind, x, f (x), fprime (x), sign);
    67   } */
    6863  else if (convtype_ != CURRENT_ONLY) {
    6964
  • branches/Couenne/Couenne/src/convex/generateCuts.cpp

    r501 r506  
    2121                                        OsiCuts &cs,
    2222                                        const CglTreeInfo info) const {
    23  
     23
    2424  Bonmin::BabInfo * babInfo = dynamic_cast <Bonmin::BabInfo *> (si.getAuxiliaryInfo ());
    2525
    2626  if (babInfo)
    2727    babInfo -> setFeasibleNode ();
    28 
    29   // lift bound on objective auxiliary to avoid overly strict implied
    30   // bounds
    3128
    3229  double now = CoinCpuTime ();
     
    134131  int objInd = problem_ -> Obj (0) -> Body () -> Index ();
    135132  if (objInd < 0) objInd = 0;
    136   /*
    137   if (babInfo) {
    138    
    139     CouNumber primal = babInfo -> mipBound (),
    140               dual   = babInfo -> bestObjectiveValue ();
     133
     134  if (babInfo && (babInfo -> babPtr ())) {
     135
     136    CouNumber primal  = babInfo  -> babPtr () -> model (). getObjValue(),
     137              dual    = babInfo  -> babPtr () -> model (). getBestPossibleObjValue (),
     138              primal0 = problem_ -> Ub (objInd),
     139              dual0   = problem_ -> Lb (objInd);
    141140
    142141    // Bonmin assumes minimization. Hence, primal (dual) is an UPPER
    143142    // (LOWER) bound.
    144143
    145     if (problem_ -> Ub (objInd) > primal) { // update primal bound (MIP)
     144    if ((primal <   COUENNE_INFINITY - 1) && (primal < primal0)) { // update primal bound (MIP)
    146145      problem_ -> Ub (objInd) = primal;
    147       if (primal < COUENNE_INFINITY - 1) chg_bds [objInd] = 1;
    148     }
    149 
    150     if (problem_ -> Lb (objInd) < dual) { // update dual bound
     146      chg_bds [objInd] = 1;
     147    }
     148
     149    if ((dual   > - COUENNE_INFINITY + 1) && (dual   > dual0)) { // update dual bound
    151150      problem_ -> Lb (objInd) = dual;
    152       if (dual < - COUENNE_INFINITY + 1) chg_bds [objInd] = 1;
    153     }
    154   }
    155   */
    156 
    157   if (BabPtr_) { // update primal bound with best feasible solution object
    158 
    159     int objInd = problem_ -> Obj (0) -> Body () -> Index ();
    160 
    161     if (objInd >= 0) {
    162 
    163       CouNumber bestObj = babInfo -> babPtr() -> model().getObjValue();
    164      
    165       // Bonmin assumes minimization. Bonmin::Bab::bestObj () should
    166       // be considered an UPPER bound.
    167 
    168       if (problem_ -> Ub (objInd) > bestObj) {
    169 
    170         problem_ -> Ub (objInd) = bestObj;
    171         if (bestObj < COUENNE_INFINITY - 1)
    172           chg_bds [objInd] = 1;
    173       }
    174     }
    175   }
    176 
     151      chg_bds [objInd] = 1;
     152    }
     153  }
     154 
    177155
    178156  //////////////////////// BOUND TIGHTENING ///////////////////////////////////
     
    188166    ntightened = problem_ -> tightenBounds (chg_bds);
    189167
    190     // implied bounds. Only do it after other bounds have changed,
    191     // i.e. after branching or if upper bound found
    192 
    193     if (!firstcall_ || chg_bds [objInd]) {
    194       //    if (!firstcall_) {
    195 
    196       nbwtightened = problem_ -> impliedBounds (chg_bds);
     168    // implied bounds. Call also at the beginning, as some common
     169    // expression may have non-propagated bounds
     170
     171    nbwtightened = problem_ -> impliedBounds (chg_bds);
    197172     
    198       for (register int i=0; i < ncols; i++)
    199         if (expression::Lbound (i) >= expression::Ubound (i) + COUENNE_EPS) {
    200 
    201           /*printf ("Couenne: infeasible bounds on w_%d [%.12e,%.12e]\n",
    202             i, expression::Lbound (i), expression::Ubound (i));*/
    203 
    204           /*OsiColCut *infeascut = new OsiColCut;
    205           if (infeascut) {
    206             double upper = -1., lower = +1.;
    207             infeascut -> setLbs (1, &i, &lower);
    208             infeascut -> setUbs (1, &i, &upper);
    209             cs.insert (infeascut);
    210             delete infeascut;
    211             }*/
    212 
    213           //      infeasNode () = true; // make this node infeasible
    214 
    215           if (babInfo)
    216             babInfo -> setInfeasibleNode ();
    217           else printf ("warning, no babinfo\n");
    218 
    219           goto end_genCuts;
    220         }
    221     }
     173    for (register int i=0; i < ncols; i++)
     174      if (expression::Lbound (i) >= expression::Ubound (i) + COUENNE_EPS) {
     175
     176        if (babInfo)
     177          babInfo -> setInfeasibleNode ();
     178
     179        goto end_genCuts;
     180      }
    222181  } while (ntightened && nbwtightened && (niter++ < 10));
    223182  // Not ((ntightened || nbwtightened) && (niter++ < 10)), as we need
  • branches/Couenne/Couenne/src/convex/operators/conv-exprPow-envelope.cpp

    r480 r506  
    3333
    3434
     35// function k*(k-1)*x^(k-2)
     36inline static CouNumber power_k_dblprime (CouNumber x)
     37{return exponent * (exponent-1) * pow (x, exponent-2);}
     38
     39
    3540// adds convex (upper/lower) envelope to a power function
    3641
    3742void addPowEnvelope (const CouenneCutGenerator *cg, OsiCuts &cs,
    3843                     int wi, int xi,
    39                      CouNumber x, CouNumber k,
     44                     CouNumber x, CouNumber y,
     45                     CouNumber k,
    4046                     CouNumber l, CouNumber u,
    4147                     int sign) {
    4248
    43   //  int ns = cg -> nSamples ();
     49  exponent = k;
    4450
    45   exponent = k;
    46   /*
    47   if (k <= 0) k = 1;
     51  // set x to get a deeper cut (so that we get a tangent which is
     52  // orthogonal with line through current- and tangent point)
    4853
    49   if (l < - COUENNE_INFINITY + 1) {
    50     if (u > COUENNE_INFINITY - 1) {
     54  //  x = powNewton (x, y, power_k, power_k_prime, power_k_dblprime);
    5155
    52       l = - pow (POW_FACTOR, 1. / (1. + k)) * ns/2;
    53       u =   pow (POW_FACTOR, 1. / (1. + k)) * ns/2;
    54     } else
    55       l = u - pow (POW_FACTOR, 1. / (1. + k)) * ns;
    56 
    57   } else if (u > COUENNE_INFINITY - 1)
    58     u = l + pow (POW_FACTOR,   1. / (1. + k)) * ns;
    59   */
     56  if      (x<l) x=l;
     57  else if (x>u) x=u;
    6058
    6159  // limit the bounds for the envelope
  • branches/Couenne/Couenne/src/convex/operators/conv-exprPow.cpp

    r497 r506  
    3535    CouNumber base = arglist_ [0] -> Value ();
    3636
    37     if (fabs (base - M_E) < COUENNE_EPS)
     37    if (fabs (base - M_E) < COUENNE_EPS) // is base e = 2.7182818...
    3838      return p -> addAuxiliary (new exprExp (new exprClone (arglist_ [1])));
    39     else
     39    else // no? convert k^x to e^(x log (k))
    4040      return p -> addAuxiliary
    4141        (new exprExp (new exprClone
     
    4444  }
    4545  else
    46 
    47     if (arglist_ [1] -> Type () != CONST) // expression is x^y, reduce
    48                                           // to exp (y*log(x));
     46    if (arglist_ [1] -> Type () != CONST) //  x^y, convert to exp (y*log(x));
    4947      return p -> addAuxiliary
    5048        (new exprExp (new exprClone (p -> addAuxiliary
     
    5351                        new exprClone (p -> addAuxiliary
    5452                                       (new exprLog (new exprClone (arglist_ [0])))))))));
    55     else                                  // expression is x^k, return
    56                                           // as it is
    57         return p -> addAuxiliary (this);
     53
     54    else return p -> addAuxiliary (this);  // expression is x^k, return as it is
    5855}
    5956
     
    8178  int x_ind = xe  -> Index ();
    8279
    83   CouNumber x = (*xe)  (),
     80  CouNumber w = (*aux) (),
     81            x = (*xe)  (),
    8482            l = (*xle) (),
    8583            u = (*xue) ();
     
    121119  bool isInt    =            fabs (k    - (double) (intk = FELINE_round (k)))    < COUENNE_EPS,
    122120       isInvInt = !isInt && (fabs (1./k - (double) (intk = FELINE_round (1./k))) < COUENNE_EPS);
    123    
     121
    124122  // two macro-cases:
    125123
     
    152150
    153151    if (l > -COUENNE_INFINITY) {
    154       if (u > q * l) { // upper x is after "turning point", add also lower envelope
    155         addPowEnvelope (cg, cs, w_ind, x_ind, x, k, q*l, u, sign);
    156         cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l,k), q*l, safe_pow (q*l,k), sign);
    157       }
    158     else
    159       cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l,k), u, safe_pow (u,k), sign);
    160     }
    161 
    162     // check if upper part needs a concave envelope
     152      if (u > q * l) { // upper x is after "turning point", add lower envelope
     153        addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, q*l, u, sign);
     154        cg      -> addSegment (cs, w_ind, x_ind, l, safe_pow (l,k), q*l, safe_pow (q*l,k), sign);
     155      } else cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l,k), u,   safe_pow (u,k),   sign);
     156    }
     157
     158    // upper envelope
    163159
    164160    if (u < COUENNE_INFINITY) {
    165       if (l < q * u) {
    166         addPowEnvelope (cg, cs, w_ind, x_ind, x, k, l, q*u, -sign);
    167         cg    -> addSegment (cs, w_ind, x_ind, q*u, safe_pow (q*u,k), u, safe_pow (u,k), -sign);
    168       }
    169       else cg -> addSegment (cs, w_ind, x_ind, l,   safe_pow (l,k),   u, safe_pow (u,k), -sign);
     161      if (l < q * u) { // lower x is before "turning point", add upper envelope
     162        addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, q*u, -sign);
     163        cg      -> addSegment (cs, w_ind, x_ind, q*u, safe_pow (q*u,k), u, safe_pow (u,k), -sign);
     164      } else cg -> addSegment (cs, w_ind, x_ind, l,   safe_pow (l,k),   u, safe_pow (u,k), -sign);
    170165    }
    171166  }
    172167  else {
    173168
     169    printf ("==== exprPow: x=%.5f  [%.5f %.5f]\n", x, l, u);
     170   
    174171    // 2) all other cases.
    175172
     
    188185    // tails of the asymptotes.
    189186
    190     if ((k < COUENNE_EPS) && (l < - COUENNE_EPS) && (u > COUENNE_EPS)) {
     187    if ((k < 0) &&
     188        (l < - COUENNE_EPS) &&
     189        (u >   COUENNE_EPS)) {
    191190
    192191      if (!(intk % 2))
    193192        cg -> addSegment (cs, w_ind, arglist_ [0] -> Index (),
    194                     l, safe_pow (l,k), u, safe_pow (u,k), +1);
     193                          l, safe_pow (l,k), u, safe_pow (u,k), +1);
    195194      return;
    196195    }
     
    199198    // be enveloped. Standard segment and tangent cuts can be applied.
    200199
    201     // create envelope. Choose sign based on k
    202 
    203     int sign = +1;
    204 
    205     // invert sign if (k is odd negative and l<0) or (k is in [0,1])
    206     if ((l < - COUENNE_EPS) && (k < 0) && (intk % 2)
    207         || (fabs (k-0.5) < 0.5 - COUENNE_EPS) // k in [0,1]
    208         || ((u < - COUENNE_EPS) && (intk % 2) && (k > COUENNE_EPS)))
     200    // create upper envelope (segment)
     201
     202    int sign = 1; // sign based on k
     203
     204    // invert sign if
     205    if (   ((l < - COUENNE_EPS) && (intk % 2) && (k < -COUENNE_EPS)) // k<0 odd, l<0
     206        || ((u < - COUENNE_EPS) && (intk % 2) && (k >  COUENNE_EPS)) // k>0 odd, u<0
     207        || (fabs (k-0.5) < 0.5 - COUENNE_EPS))                       // k in [0,1]
    209208      sign = -1;
    210209
    211     // upper envelope -- when k negative, add only if bounds are far from 0
     210    // upper envelope -- when k negative, add only if
     211
     212    CouNumber powThreshold = pow (1e20, 1./k);
    212213
    213214    if ((  (k > COUENNE_EPS)
    214         || (l > COUENNE_EPS)
     215        || (l > COUENNE_EPS)             // bounds do not contain 0
    215216        || (u < - COUENNE_EPS)) &&
    216         (l > - COUENNE_INFINITY + 1) &&
    217         (u <   COUENNE_INFINITY - 1))
    218 
    219       cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l,k), u, safe_pow (u,k), -sign);
     217        (l > - powThreshold) &&  // and are finite
     218        (u <   powThreshold))
     219
     220      cg -> addSegment (cs, w_ind, x_ind, l, safe_pow (l, k), u, safe_pow (u, k), -sign);
    220221
    221222    // similarly, pay attention not to add infinite slopes
     
    225226    if (k > COUENNE_EPS) {
    226227
    227       if (u >   COUENNE_INFINITY - 1) u = x + POWER_RANGE;
    228       if (l < - COUENNE_INFINITY + 1) l = x - POWER_RANGE;
     228      if (u >   powThreshold) u = x + POWER_RANGE;
     229      if (l < - powThreshold) l = x - POWER_RANGE;
    229230    }
    230231    else {
     
    234235    }
    235236
    236     addPowEnvelope (cg, cs, w_ind, x_ind, x, k, l, u, sign);
     237    addPowEnvelope (cg, cs, w_ind, x_ind, x, w, k, l, u, sign);
    237238  }
    238239
  • branches/Couenne/Couenne/src/convex/operators/powNewton.cpp

    r497 r506  
    7272  expon = atof (argv [1]);
    7373
    74   //  for (register int i=10000; i--;)
    75   r = powNewton (xc, yc, f, fp, fpp);
     74  for (register int i=10000; i--;)
     75    r = powNewton (xc, yc, f, fp, fpp);
    7676
    7777  printf ("xc = %.14f: xk = %.15f, slope %.15f -- %.15f ==> %.15f\n",
  • branches/Couenne/Couenne/src/expression/exprAux.cpp

    r498 r506  
    4747  if ((!(cg -> isFirst ())) &&
    4848      (fabs ((l = expression::Ubound (varIndex_)) -
    49              expression::Lbound (varIndex_)) < COUENNE_EPS))
    50 
     49                  expression::Lbound (varIndex_)) < COUENNE_EPS))
    5150    cg -> createCut (cs, l, 0, varIndex_, 1.);
    5251  else image_ -> generateCuts (this, si, cs, cg);
     
    5554  //  if (j < cs.sizeRowCuts ())
    5655  if (0)
    57   //  if (varIndex_ == 60)
    5856    {
    5957      printf ("----------------Generated cut for ");
     
    6765  //////////////////////////////////////////////////////////////
    6866
    69   if (0) {// [cool!] print graph-readable output for displaying
    70           // inequalities on a Cartesian plane
     67  if (0) { // [cool!] print graph-readable output for displaying
     68           // inequalities on a Cartesian plane
    7169
    72     if ((image_ -> code () == COU_EXPRSIN) ||
     70    if (1 || (image_ -> code () == COU_EXPRSIN) ||
    7371        (image_ -> code () == COU_EXPRPOW) ||
     72        (image_ -> code () == COU_EXPREXP) ||
    7473        (image_ -> code () == COU_EXPRLOG) ||
    7574        (image_ -> code () == COU_EXPRCOS)) {
     
    118117            if (y < minY) minY = y;
    119118
    120             printf ("#=# %.3f %.3f\n", x, y);
     119            printf ("#=# %.12e %.12e\n", x, y);
    121120          }
    122121        }
    123122       
    124         lb -= 1;
    125         ub += 1;
     123        //lb -= 1;
     124        //ub += 1;
    126125
    127126        // plot lines defining constraint (only for cuts involving at
     
    142141          printf ("#=# #m=2,S=%d\n", (cs.rowCutPtr (jj) -> sense () == 'L') ? 10:11);
    143142
    144           printf ("#=# %.3f %.3f\n", lb0, (rhs - el [1] * lb0) / el [0]);
    145           printf ("#=# %.3f %.3f\n", ub0, (rhs - el [1] * ub0) / el [0]);
     143          printf ("#=# %.12e %.12e\n", lb0, (rhs - el [1] * lb0) / el [0]);
     144          printf ("#=# %.12e %.12e\n", ub0, (rhs - el [1] * ub0) / el [0]);
    146145        }
    147146
    148         exit (0);
    149147        cg -> Problem () -> X () [xi] = curx;
    150148      }
  • branches/Couenne/Couenne/src/include/operators/exprPow.h

    r450 r506  
    9494
    9595// compute power
    96 
    9796inline CouNumber exprPow::operator () () {
    9897
     
    107106
    108107// add upper/lower envelope to power in convex/concave areas
     108void addPowEnvelope (const CouenneCutGenerator *, OsiCuts &, int, int,
     109                     CouNumber, CouNumber, CouNumber, CouNumber, CouNumber, int);
    109110
    110 void addPowEnvelope (const CouenneCutGenerator *, OsiCuts &, int, int,
    111                      CouNumber, CouNumber, CouNumber, CouNumber, int);
     111
     112// find proper tangent point to add deepest tangent cut
     113CouNumber powNewton (CouNumber, CouNumber, unary_function, unary_function, unary_function);
     114
    112115#endif
Note: See TracChangeset for help on using the changeset viewer.