Changeset 303


Ignore:
Timestamp:
Feb 13, 2010 5:57:58 PM (11 years ago)
Author:
pbelotti
Message:

separation of cuts for bilinear terms

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Couenne/src/convex/operators/exprMul-upperHull.cpp

    r296 r303  
    2424                  CouNumber x2, CouNumber y2,
    2525                  char whichUse,
    26                   CouNumber &cX, CouNumber &cY, CouNumber &cW, CouNumber &c0);
     26                  CouNumber &cX, CouNumber &cY, CouNumber &cW);
     27
     28
     29// invert interval bounds and current point
     30inline void invertInterval (register double &l, register double &u, register double x) {
     31
     32  register double tmp = l;
     33  l  = -u;
     34  u  = -tmp;
     35
     36  x  = -x;
     37}
    2738
    2839void upperEnvHull (const CouenneCutGenerator *cg, OsiCuts &cs,
     
    3243
    3344  //
    34   // See upcoming paper for explanation ;-)
     45  // See forthcoming paper for explanation ;-)
    3546  //
    3647
     
    5162  if (y0 < yl) y0 = yl;  if (y0 > yu) y0 = yu;
    5263
     64
     65  // preliminary bound tightening
    5366  if (wl >= 0) {
    5467    if        ((xl >= 0) || (yl >= 0) || (xl * yl < wl - COUENNE_EPS)) {
     
    7487    flipX = xl < 0,
    7588    flipY = yl < 0,
    76     flipW;
    77 
    78   {
    79     CouNumber tmp, tmp2;
    80 
    81     if (flipX && flipY) { // swap bounds on x and y
    82 
    83       tmp =  xl;     tmp2 =  yl;
    84       xl  = -xu;     yl   = -yu;
    85       xu  = -tmp;    yu   = -tmp;
    86 
    87       x0  = -x0;     y0   = -y0;
    88 
    89     } else if (flipX) { // swap bounds on x and w only
    90 
    91       tmp =  xl;     tmp2 =  wl;
    92       xl  = -xu;     wl   = -wu;
    93       xu  = -tmp;    wu   = -tmp;
    94 
    95       x0  = -x0;     w0   = -w0;
    96 
    97       flipW = true;
    98 
    99     } else if (flipY) { // swap bounds on y and w only
    100 
    101       tmp =  yl;     tmp2 =  wl;
    102       yl  = -yu;     wl   = -wu;
    103       yu  = -tmp;    wu   = -tmp;
    104 
    105       y0  = -y0;     w0   = -w0;
    106 
    107       flipW = true;
    108     }
     89    flipW = false;
     90
     91  if (flipX && flipY) { // swap bounds on x and y
     92
     93    invertInterval (xl,xu,x0);
     94    invertInterval (yl,yu,y0);
     95
     96  } else if (flipX) { // swap bounds on x and w only
     97
     98    invertInterval (xl,xu,x0);
     99    invertInterval (wl,wu,w0);
     100
     101    flipW = true;
     102
     103  } else if (flipY) { // swap bounds on y and w only
     104
     105    invertInterval (yl,yu,y0);
     106    invertInterval (wl,wu,w0);
     107
     108    flipW = true;
    109109  }
    110110
     
    147147
    148148  CouNumber
    149     cX, cY, cW, c0,
    150     cXp, cYp, cWp, c0p; // these only if two inequalities are added
    151 
    152   bool twoIneqs = false;
     149    cX,  cY,  cW,  c0,  c0X,  c0Y,  c0W,
     150    cXp, cYp, cWp, c0p, c0Xp, c0Yp, c0Wp; // these only if two inequalities are added
     151
     152  bool
     153    twoIneqs   = false, // generate two inequalities for weird case
     154    upperCurve = false; // has lifting taken place from upper curve?
    153155
    154156  if (xLow >= xl && xUpp <= xu &&
    155157      yLow >= yl && yUpp <= yu) {
    156158
    157     //printf ("easy lifting:\n");
     159#ifdef DEBUG
     160    printf ("easy lifting:\n");
     161#endif
    158162
    159163    // Case 2: both are inside. Easy lifting...
    160     if (genMulCoeff (xLow, yLow, xUpp, yUpp, 0, cX, cY, cW, c0))
    161       return;
     164    if (genMulCoeff (xLow, yLow, xUpp, yUpp, 0, cX, cY, cW))
     165      return;
     166
     167    c0X = cX * xLow;
     168    c0Y = cY * yLow;
     169    c0W = cW * wl;
    162170
    163171  } else if (xLow >= xl && yLow >= yl &&
    164172             (xUpp > xu || yUpp > yu)) {
    165173
    166     //printf ("through lower, not upper:\n");
     174#ifdef DEBUG
     175    printf ("through lower, not upper:\n");
     176#endif
    167177
    168178    // Case 3a and 3b: through lower curve, but not upper.
     
    184194
    185195    // Otherwise, lift inequality on lower point
    186     if (genMulCoeff (xLow, yLow, xUpp, yUpp, 0, cX, cY, cW, c0))
    187       return;
     196    if (genMulCoeff (xLow, yLow, xUpp, yUpp, 0, cX, cY, cW))
     197      return;
     198
     199    c0X = cX * xLow;
     200    c0Y = cY * yLow;
     201    c0W = cW * wl;
    188202
    189203  } else if (xUpp <= xu && yUpp <= yu &&
    190204             (xLow < xl || yLow < yl)) {
    191205
    192     //printf ("through upper, not lower:\n");
     206#ifdef DEBUG
     207    printf ("through upper, not lower:\n");
     208#endif
     209
    193210    // Case 4a and 4b: viceversa (lift for validity)
    194211
     
    208225      return;
    209226
     227    upperCurve = true;
     228
    210229    // Otherwise, lift inequality on UPPER point
    211     if (genMulCoeff (xLow, yLow, xUpp, yUpp, 0, cX, cY, cW, c0))
    212       return;
     230    if (genMulCoeff (xLow, yLow, xUpp, yUpp, 1, cX, cY, cW))
     231      return;
     232
     233    c0X = cX * xUpp;
     234    c0Y = cY * yUpp;
     235    c0W = cW * wu;
    213236
    214237  } else if ((xLow < xl && xUpp > xu) ||
    215238             (yLow < yl && yUpp > yu)) {
    216239
    217     //printf ("between lower and upper:\n");
     240#ifdef DEBUG
     241    printf ("between lower and upper:\n");
     242#endif
    218243
    219244    // Case 5: both outside of bounding box, N and S or W and
     
    236261    // Nothing to find. Just separate two inequalities at the same
    237262    // point, just using different support
    238     if (genMulCoeff (xLow, yLow, xUpp, yUpp, 0, cX,  cY,  cW,  c0) ||
    239         genMulCoeff (xLow, yLow, xUpp, yUpp, 1, cXp, cYp, cWp, c0p))
    240       return;
     263    if (genMulCoeff (xLow, yLow, xUpp, yUpp, 0, cX,  cY,  cW) ||
     264        genMulCoeff (xLow, yLow, xUpp, yUpp, 1, cXp, cYp, cWp))
     265      return;
     266
     267#ifdef DEBUG
     268    printf ("coeffs: (%g,%g,%g) [(%g,%g,%g)]\n",
     269            cX,cY,cW, cXp, cYp, cWp);
     270#endif
     271
     272    c0X = cX * xLow;    c0Xp = cXp * xUpp;
     273    c0Y = cY * yLow;    c0Yp = cYp * yUpp;
     274    c0W = cW * wl;      c0Wp = cWp * wu; 
    241275
    242276    twoIneqs = true;
     
    255289  // Re-transform back into original variables
    256290
    257   if (flipX) {cX = -cX; cXp = -cXp;}
    258   if (flipY) {cY = -cY; cYp = -cYp;}
    259   if (flipW) {cW = -cW; cWp = -cWp;}
     291  if (flipX) {cX = -cX; cXp = -cXp; c0X = -c0X; c0Xp = -c0Xp;}
     292  if (flipY) {cY = -cY; cYp = -cYp; c0Y = -c0Y; c0Yp = -c0Yp;}
     293  if (flipW) {cW = -cW; cWp = -cWp; c0W = -c0W; c0Wp = -c0Wp;}
     294
     295  c0  = c0X  + c0Y  + c0W;
     296
     297#ifdef DEBUG
     298  printf ("there are cuts\n");
     299#endif
    260300
    261301  //cg -> createCut (cs, alpha*wb + 2*wb/xt, sign, wi, alpha, yi, 1., xi, wb/(xt*xt));
    262   cg   -> createCut (cs, -c0,  flipW ? +1 : -1, wi, cW,  yi, -cY,  xi, -cX);
    263 
    264   if (twoIneqs)
    265     cg -> createCut (cs, -c0p, flipW ? +1 : -1, wi, cWp, yi, -cYp, xi, -cXp);
     302  cg   -> createCut (cs, c0,  1, wi, cW,  yi, cY,  xi, cX);
     303
     304  if (twoIneqs) {
     305    c0p = c0Xp + c0Yp + c0Wp;
     306    cg -> createCut (cs, c0p, 1, wi, cWp, yi, cYp, xi, cXp);
     307  }
    266308}
    267309
     
    276318                      CouNumber *xB, CouNumber *yB) {
    277319
    278   // The parametric curve is of the form
     320  // The parametric line is of the form
    279321  //
    280322  //  x = x0 + t (x1-x0)
     
    333375                  CouNumber x2, CouNumber y2,
    334376                  char whichUse,
    335                   CouNumber &cX, CouNumber &cY, CouNumber &cW, CouNumber &c0) {
     377                  CouNumber &cX, CouNumber &cY, CouNumber &cW) {
    336378
    337379  // the x-y slope of this constraint must be tangent to a curve xy=k
     
    350392  cX = yD;
    351393  cY = xD;
    352   c0 = 2*xD*yD;
     394  //c0 = 2*xD*yD;
     395
     396#ifdef DEBUG
     397    printf ("points: (%g,%g) (%g,%g), cW = (%g - %g) / %g = %g\n",
     398            xD,yD, xO,yO, 2*xD*yD, (cX*xO + cY*yO), (xO*yO - xD*yD),
     399            (2*xD*yD - (cX*xO + cY*yO)) / (xO*yO - xD*yD));
     400#endif
    353401
    354402  // Now the hard part: lift it so that it touches the other curve
    355403
    356404  if (fabs (xO*yO - xD*xD) < COUENNE_EPS)
    357     return 1;
    358 
    359   cW = (c0 - (cX*xO + cY*yO)) / (xO*yO - xD*xD);
    360 
    361   c0 += cW * xD*yD;
    362 
    363   if ((xO < xD) || (yO < yD))
    364     cW = -cW;
     405    return 1; // no cut if the two points are on the same curve
     406
     407  // should ALWAYS be negative
     408  cW = (2*xD*yD - (cX*xO + cY*yO)) / (xO*yO - xD*yD);
     409
     410  //c0 += cW * xD*yD;
     411
     412  //if ((xO < xD) || (yO < yD))
     413  //cW = -cW;
     414
     415  //if (xO*yO < xD*yD)
     416  //cW = -cW;
    365417
    366418  return 0;
Note: See TracChangeset for help on using the changeset viewer.