source: trunk/Couenne/src/convex/operators/exprMul-upperHull.cpp @ 312

Last change on this file since 312 was 312, checked in by pbelotti, 10 years ago

gave Couenne prefix to all header files

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 11.0 KB
Line 
1/* $Id: exprMul-upperHull.cpp 312 2010-03-10 23:43:21Z pbelotti $
2 *
3 * Name:    exprMul-upperHull.cpp
4 * Author:  Pietro Belotti
5 * Purpose: generates upper envelope of a product
6 *
7 * (C) Carnegie-Mellon University, 2010.
8 * This file is licensed under the Common Public License (CPL)
9 */
10
11#include "CouenneExprMul.hpp"
12#include "CouennePrecisions.hpp"
13#include "CouenneCutGenerator.hpp"
14
15//#define DEBUG
16
17int findIntersection (CouNumber  x0, CouNumber  y0,
18                      CouNumber  x1, CouNumber  y1,
19                      CouNumber *wl, CouNumber *wu,
20                      CouNumber *xA, CouNumber *yA,
21                      CouNumber *xB, CouNumber *yB);
22
23int genMulCoeff (CouNumber x1, CouNumber y1, 
24                  CouNumber x2, CouNumber y2, 
25                  char whichUse,
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}
38
39void upperEnvHull (const CouenneCutGenerator *cg, OsiCuts &cs, 
40                   int xi, CouNumber x0, CouNumber xl, CouNumber xu,
41                   int yi, CouNumber y0, CouNumber yl, CouNumber yu,
42                   int wi, CouNumber w0, CouNumber wl, CouNumber wu) {
43
44  //
45  // See forthcoming paper for explanation ;-)
46  //
47
48#ifdef DEBUG
49  printf ("entering points: ===================================================\n");
50  printf ("x [%d]: %9g\t[%9g\t%9g]\n", xi, x0, xl, xu);
51  printf ("y [%d]: %9g\t[%9g\t%9g]\n", yi, y0, yl, yu);
52  printf ("w [%d]: %9g\t[%9g\t%9g]\n", wi, w0, wl, wu);
53#endif
54
55  // Preprocess to reduce everything to a first-orthant problem
56
57  if ((wl <  0 && wu >  0)) // nothing to tighten
58    return;
59
60  // project back into bounding box
61  if (x0 < xl) x0 = xl;  if (x0 > xu) x0 = xu;
62  if (y0 < yl) y0 = yl;  if (y0 > yu) y0 = yu;
63
64
65  // preliminary bound tightening
66  if (wl >= 0) {
67    if        ((xl >= 0) || (yl >= 0) || (xl * yl < wl - COUENNE_EPS)) {
68      if (xl < 0) xl = 0;
69      if (yl < 0) yl = 0;
70    } else if ((xu <= 0) || (yu <= 0) || (xu * yu < wl - COUENNE_EPS)) {
71      if (xu > 0) xu = 0;
72      if (yu > 0) yu = 0;
73    }
74  } else {
75    if        ((xl >= 0) || (yu <= 0) || (xl * yu > wu + COUENNE_EPS)) {
76      if (xl < 0) xl = 0;
77      if (yu > 0) yu = 0;
78    } else if ((xu <= 0) || (yl >= 0) || (xu * yl > wu + COUENNE_EPS)) {
79      if (xu > 0) xu = 0;
80      if (yl < 0) yl = 0;
81    }
82  }
83
84  // Reduce
85
86  bool 
87    flipX = xl < 0,
88    flipY = yl < 0,
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;
109  }
110
111#ifdef DEBUG
112  printf ("reduced points:\n");
113  printf ("x: %9g\t[%9g\t%9g]\n", x0, xl, xu);
114  printf ("y: %9g\t[%9g\t%9g]\n", y0, yl, yu);
115  printf ("w: %9g\t[%9g\t%9g]\n", w0, wl, wu);
116#endif
117
118  // Check whether lower and upper curve both contain bounding box of
119  // x,y. If so, there is nothing to separate
120
121  if (((xl*yl >= wl) &&  // b.box totally contained between two curves
122       (xu*yu <= wu)) || //
123      (x0*y0 >= w0))     // or current point below curve
124    return;
125
126  // Find intersections of halfline from origin
127
128  CouNumber xLow, yLow, xUpp, yUpp;
129  if (findIntersection (0., 0., x0, y0, &wl, &wu, &xLow, &yLow, &xUpp, &yUpp))
130    return; // couldn't find proper point
131
132#ifdef DEBUG
133  printf ("intersections:\n");
134  printf ("lower: %9g\t%9g\tprod %9g\n", xLow, yLow, xLow*yLow);
135  printf ("upper: %9g\t%9g\tprod %9g\n", xUpp, yUpp, xUpp*yUpp);
136#endif
137
138  // Case 1: both are outside of bounding box, either NW or SE. McCormick's suffice.
139
140  if ((xLow <= xl && yUpp >= yu) ||
141      (yLow <= yl && xUpp >= xu))
142    return;
143
144  // OK, if you are here there will be at least one cut. Define
145  // coefficients and rhs --- will have to change them back if (flipX
146  // || flipY)
147
148  CouNumber
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?
155
156  if (xLow >= xl && xUpp <= xu &&
157      yLow >= yl && yUpp <= yu) {
158
159#ifdef DEBUG
160    printf ("easy lifting:\n");
161#endif
162
163    // Case 2: both are inside. Easy lifting...
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;
170
171  } else if (xLow >= xl && yLow >= yl && 
172             (xUpp > xu || yUpp > yu)) {
173
174#ifdef DEBUG
175    printf ("through lower, not upper:\n");
176#endif
177
178    // Case 3a and 3b: through lower curve, but not upper.
179
180    if (yUpp > yu) { // upper intersect is North; place it within box
181      yUpp = yu;
182      xUpp = wu / yu;
183    } else {         //                    East
184      xUpp = xu;
185      yUpp = wu / xu;
186    }
187
188    // find intersection on low curve on half line through new point and (x0,y0)
189    if (findIntersection (xUpp, yUpp, x0, y0, &wl, NULL, &xLow, &yLow, NULL, NULL))
190      return;
191
192    if (xLow < xl || yLow < yl) // McCormick's suffice
193      return;
194
195    // Otherwise, lift inequality on lower point
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;
202
203  } else if (xUpp <= xu && yUpp <= yu && 
204             (xLow < xl || yLow < yl)) {
205
206#ifdef DEBUG
207    printf ("through upper, not lower:\n");
208#endif
209
210    // Case 4a and 4b: viceversa (lift for validity)
211
212    if (yLow < yl) { // upper intersect is South; place it within box
213      yLow = yl;
214      xLow = wl / yl;
215    } else {         //                    West
216      xLow = xl;
217      yLow = wl / xl;
218    }
219
220    // find intersection on low curve on half line through new point and (x0,y0)
221    if (findIntersection (xLow, yLow, x0, y0, NULL, &wu, NULL, NULL, &xUpp, &yUpp))
222      return;
223
224    if (xUpp > xu || yUpp > yu) // McCormick's suffice
225      return;
226
227    upperCurve = true;
228
229    // Otherwise, lift inequality on UPPER point
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;
236
237  } else if ((xLow < xl && xUpp > xu) ||
238             (yLow < yl && yUpp > yu)) {
239
240#ifdef DEBUG
241    printf ("between lower and upper:\n");
242#endif
243
244    // Case 5: both outside of bounding box, N and S or W and
245    //         E. Separate both from lower and from upper
246
247    if (yLow < yl) { // upper intersect is South; place it within box
248      yLow = yl;      xLow = wl / yl;
249      yUpp = yu;      xUpp = wu / yu;
250    } else {         //                    West
251      xLow = xl;      yLow = wl / xl;
252      xUpp = xu;      yUpp = wu / xu;
253    }
254
255#ifdef DEBUG
256    printf ("New intersections:\n");
257    printf ("lower: %9g\t%9g\tprod %9g\n", xLow, yLow, xLow*yLow);
258    printf ("upper: %9g\t%9g\tprod %9g\n", xUpp, yUpp, xUpp*yUpp);
259#endif
260
261    // Nothing to find. Just separate two inequalities at the same
262    // point, just using different support
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; 
275
276    twoIneqs = true;
277
278  } else {
279
280#ifdef DEBUG
281    printf ("points are in a weird position:\n");
282    printf ("lower: %9g\t%9g\tprod %9g\n", xLow, yLow, xLow*yLow);
283    printf ("upper: %9g\t%9g\tprod %9g\n", xUpp, yUpp, xUpp*yUpp);
284#endif
285
286    return;
287  }
288
289  // Re-transform back into original variables
290
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
300
301  //cg -> createCut (cs, alpha*wb + 2*wb/xt, sign, wi, alpha, yi, 1., xi, wb/(xt*xt));
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  }
308}
309
310
311// finds intersections of a parametric line (x,y) = (x0,y0) + t [(x1,y1) - (x0,y0)]
312// on curves xy = wl and xy = yu
313
314int findIntersection (CouNumber  x0, CouNumber  y0,
315                      CouNumber  x1, CouNumber  y1,
316                      CouNumber *wl, CouNumber *wu,
317                      CouNumber *xA, CouNumber *yA,
318                      CouNumber *xB, CouNumber *yB) {
319
320  // The parametric line is of the form
321  //
322  //  x = x0 + t (x1-x0)
323  //  y = y0 + t (y1-y0)
324  //
325  // and for that to satisfy xy = wl and xy = wu we must have
326  //
327  // x0 y0 + t [x0 (y1-y0) + y0 (x1-x0)] + t^2 (x1-x0) (y1-y0) = wl
328  //                                                           = wu
329  //
330  // or a t^2 + b t + c - wl = 0 for proper values of a,b,c.
331  //    a t^2 + b t + c - wu = 0
332  //
333  // Because of the way this procedure will be used, of the two
334  // solutions found we must always use the minimum nonnegative one
335
336  CouNumber
337    a = (x1-x0) * (y1-y0),
338    c = x0 * y0,
339    b = x0*y1 + y0*x1 - 2*c; // x0 * (y1-y0) + y0 * (x1-x0)
340
341  if (fabs (a) < COUENNE_EPS)
342    return 1;
343
344  // These are the solution to the above equation.
345
346  CouNumber tL1, tL2, tU1, tU2, tL, tU;
347
348  if (wl) {
349    tL1 = (- b - sqrt (b*b - 4*a*(c-*wl))) / (2*a);
350    tL2 = (- b + sqrt (b*b - 4*a*(c-*wl))) / (2*a);
351    //printf ("solutions L: %g %g (b2-4ac=%g)\n", tL1, tL2, b*b - 4*a*(c-*wl));
352    tL = (tL1 < 0) ? tL2 : tL1;
353  }
354
355  if (wu) {
356    tU1 = (- b - sqrt (b*b - 4*a*(c-*wu))) / (2*a);
357    tU2 = (- b + sqrt (b*b - 4*a*(c-*wu))) / (2*a);
358    //printf ("solutions U: %g %g (b2-4ac=%g)\n", tU1, tU2, b*b - 4*a*(c-*wu));
359    tU = (tU1 < 0) ? tU2 : tU1;
360  }
361
362  if (xA) *xA = x0 + tL * (x1-x0);   if (yA) *yA = y0 + tL * (y1-y0);
363  if (xB) *xB = x0 + tU * (x1-x0);   if (yB) *yB = y0 + tU * (y1-y0);
364
365  return 0;
366}
367
368
369// generate coefficients for a plane through points (x1, y1, x1 y1)
370// and (x2, y2, x2 y2) such that intersecting it with one of them (the
371// first if whichUse==0, the second otherwise) gives a tangent to the
372// curve xy = k.
373
374int genMulCoeff (CouNumber x1, CouNumber y1, 
375                  CouNumber x2, CouNumber y2, 
376                  char whichUse,
377                  CouNumber &cX, CouNumber &cY, CouNumber &cW) {
378
379  // the x-y slope of this constraint must be tangent to a curve xy=k
380  // at (xD,yD). Easy:
381
382  CouNumber xD, yD, xO, yO;
383
384  if (!whichUse) {
385    xD = x1; xO = x2;
386    yD = y1; yO = y2;
387  } else {
388    xD = x2; xO = x1;
389    yD = y2; yO = y1;
390  }
391
392  cX = yD;
393  cY = xD;
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
401
402  // Now the hard part: lift it so that it touches the other curve
403
404  if (fabs (xO*yO - xD*xD) < COUENNE_EPS) 
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;
417
418  return 0;
419}
Note: See TracBrowser for help on using the repository browser.