source: trunk/Couenne/src/branch/operators/branchExprMul.cpp @ 591

Last change on this file since 591 was 591, checked in by pbelotti, 9 years ago

probing now applied even in absence of NLP solution

  • Property svn:keywords set to Author Date Id Revision
File size: 10.2 KB
Line 
1/* $Id: branchExprMul.cpp 591 2011-05-31 17:03:03Z pbelotti $
2 *
3 * Name:    branchExprMul.cpp
4 * Author:  Pietro Belotti
5 * Purpose: return branch data for multiplications
6 *
7 * (C) Carnegie-Mellon University, 2006-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11#include "CouennePrecisions.hpp"
12#include "CouenneTypes.hpp"
13#include "CouenneObject.hpp"
14
15#include "CouenneExprMul.hpp"
16#include "CouenneFunTriplets.hpp"
17#include "CouenneProjections.hpp"
18
19using namespace Couenne;
20
21/// set up branching object by evaluating many branching points for
22/// each expression's arguments
23CouNumber exprMul::selectBranch (const CouenneObject *obj,
24                                 const OsiBranchingInformation *info,
25                                 expression *&var,
26                                 double * &brpts, 
27                                 double * &brDist, // distance of current LP
28                                                   // point to new convexifications
29                                 int &way) {
30
31  if (brDist) {free (brDist); brDist = NULL;} // clear it, computeMulBrDist will fill it
32
33  int xi = arglist_ [0] -> Index (),
34      yi = arglist_ [1] -> Index (),
35      wi = obj -> Reference () -> Index ();
36
37  assert ((xi >= 0) && (yi >= 0) && (wi >= 0));
38
39  CouNumber
40    x0 = info -> solution_  [xi], y0 = info -> solution_  [yi],
41    xl = info -> lower_     [xi], yl = info -> lower_     [yi],
42    xu = info -> upper_     [xi], yu = info -> upper_     [yi];
43  //w0 = info -> solution_  [wi];
44
45#ifdef DEBUG
46  printf ("    branch MUL: %g [%g,%g] %g [%g,%g]\n", 
47          x0, xl, xu, y0, yl, yu);
48#endif
49
50  brpts = (double *) realloc (brpts, sizeof (double));
51
52  // Constant x and/or y //////////////////////////////////////////////////////////
53
54  if (fabs (xu-xl) < COUENNE_EPS) { // x almost constant
55
56    if (fabs (yu-yl) < COUENNE_EPS) { // both almost constant, return null result
57
58      var = NULL;
59      return 0.;
60
61    } else { // x constant, branch on y
62
63      var = arglist_ [1];
64      *brpts = 0.5 * (yl+yu);
65      brDist = (double *) realloc (brDist, 2 * sizeof (double));
66
67      brDist [0] = projectSeg (x0, y0, yl, xl*yl, *brpts, *brpts * xl, 0);
68      brDist [1] = projectSeg (x0, y0, *brpts, *brpts * xl, yu, xl*yu, 0);
69
70      //return fabs (w0 - x0*y0);
71      return CoinMin (brDist [0], brDist [1]);
72    }
73
74  } else if (fabs (yu-yl) < COUENNE_EPS) { // y constant, branch on x
75
76    var = arglist_ [0];
77    *brpts = 0.5 * (xl+xu);
78    brDist = (double *) realloc (brDist, 2 * sizeof (double));
79
80    brDist [0] = projectSeg (x0, y0, xl, xl*yl, *brpts, *brpts * yl, 0);
81    brDist [1] = projectSeg (x0, y0, *brpts, *brpts * yl, xu, xu*yl, 0);
82
83    //return fabs (w0 - x0*y0);
84    return CoinMin (brDist [0], brDist [1]);
85  }
86
87  // Unbounded x and/or y /////////////////////////////////////////////////////////
88
89  if ((((var = arglist_ [0]) -> Index() >= 0) && (xl < -COUENNE_INFINITY) && (xu > COUENNE_INFINITY)) ||
90      (((var = arglist_ [1]) -> Index() >= 0) && (yl < -COUENNE_INFINITY) && (yu > COUENNE_INFINITY))) {
91
92    *brpts = 0.;
93    brDist = computeMulBrDist (info, xi, yi, wi, var -> Index (), brpts);
94    way = (info -> solution_ [var -> Index ()] > *brpts) ? TWO_RIGHT : TWO_LEFT;
95
96    return CoinMin (brDist [0], brDist [1]);
97  }
98
99  // TODO: don't privilege xi over yi
100
101  // at most one bound is infinite ///////////////////////////////////////////////
102
103  int ind = -1;
104
105  if      (xl < -COUENNE_INFINITY)                              // x unbounded below
106    {ind = xi; *brpts = obj -> midInterval (((x0 < 0.) ? 2 : 0.5) * x0, xl, xu); way = TWO_RIGHT;}
107
108  else if (xu >  COUENNE_INFINITY)                              // x unbounded above
109    {ind = xi; *brpts = obj -> midInterval (((x0 > 0.) ? 2 : 0.5) * x0, xl, xu); way = TWO_LEFT;} 
110
111  else if (yl < -COUENNE_INFINITY)                              // y unbounded below
112    {ind = yi; *brpts = obj -> midInterval (((y0 < 0.) ? 2 : 0.5) * y0, yl, yu); way = TWO_RIGHT;}
113
114  else if (yu >  COUENNE_INFINITY)                              // y unbounded above
115    {ind = yi; *brpts = obj -> midInterval (((y0 > 0.) ? 2 : 0.5) * y0, yl, yu) ;way = TWO_LEFT;} 
116
117  else { // both are bounded
118
119    // Check if, though bounded, they are a bit too large
120
121    CouNumber delta = (yu-yl) - (xu-xl);
122
123    if      (delta > +COUENNE_EPS) ind = yi;
124    else if (delta < -COUENNE_EPS) ind = xi;
125    else ind = (CoinDrand48 () < 0.5) ? xi : yi;
126
127    CouNumber
128      pt = info -> solution_  [ind],
129      lb = info -> lower_     [ind],
130      ub = info -> upper_     [ind],
131      margin = obj -> lp_clamp () * (ub - lb);
132
133    if ((lb < -COUENNE_EPS) && 
134        (ub >  COUENNE_EPS) && 
135        (-lb/ub >= THRES_ZERO_SYMM) &&
136        (-ub/lb >= THRES_ZERO_SYMM))
137      // interval is fairly symmetric around 0, branch on it
138      *brpts = 0.;
139
140#define LARGE_VALUE 1e8
141
142    else if ((lb < - LARGE_VALUE) && 
143             (ub >   LARGE_VALUE) &&
144             (fabs (pt) > LARGE_VALUE / 10))
145      *brpts = 0.;
146
147    else switch (obj -> Strategy ()) {
148      case CouenneObject::LP_CENTRAL:   *brpts = pt; if ((pt < lb + margin) || 
149                                                         (pt > ub - margin)) 
150                                                       pt = .5 * (lb+ub);                     break;
151      case CouenneObject::LP_CLAMPED:   *brpts = CoinMax (lb + margin, 
152                                                 CoinMin (ub - margin, pt));                  break;
153      case CouenneObject::MID_INTERVAL: *brpts = obj -> midInterval (pt, lb, ub);             break;
154      case CouenneObject::BALANCED:     *brpts = balancedMul (info, (ind == xi) ? 0 : 1, wi); break;
155      case CouenneObject::MIN_AREA: // in products, the minimum volume
156                                    // subdivision is at the middle of
157                                    // the interval
158      default:                          *brpts = (0.5 * (lb+ub));                             break;
159    }
160
161    way = (pt > *brpts) ? TWO_RIGHT : TWO_LEFT;
162  }
163
164  assert (ind >= 0);
165
166  var = arglist_ [(ind == xi) ? 0 : 1];
167
168  brDist = computeMulBrDist (info, xi, yi, wi, ind, brpts);
169
170#ifdef DEBUG
171  printf ("    MUL: br on x_%d %g [%g,%g] [%g,%g] (%g,%g)\n", 
172          ind, *brpts, xl, xu, yl, yu, x0, y0);
173#endif
174
175  return CoinMin (brDist [0], brDist [1]);
176  //return fabs (w0 - x0*y0);
177}
178
179
180// branching point for multiplication according to the balanced strategy
181CouNumber exprMul::balancedMul (const OsiBranchingInformation *info, int index, int wind) {
182
183  // first of all, make sure both arguments are variables
184
185  int other;
186
187  if (index==0) {
188    index = arglist_ [0] -> Index ();
189    other = arglist_ [1] -> Index ();
190  } else {
191    index = arglist_ [1] -> Index ();
192    other = arglist_ [0] -> Index ();
193  }
194
195  assert ((index >= 0) && (other >= 0));
196
197  CouNumber
198    xl = info -> lower_    [index],  yl = info -> lower_    [other],
199    xu = info -> upper_    [index],  yu = info -> upper_    [other],
200    x0 = info -> solution_ [index],  y0 = info -> solution_ [other],
201    w0 = info -> solution_ [wind];
202   
203  // It is quite tricky to implement a balanced strategy for products,
204  // because it is more difficult to measure "balancedness" for binary
205  // operators than it is for univariate functions.
206  //
207  // As a rule of thumb, we therefore apply the usual balanced
208  // strategy for the univariate function resulting from constraining
209  // (x,y) to a segment crossing their bounding box [xl,xu] X [yl,yu].
210  //
211  // Said segment is the set of points between (xl,yl) and (xu,yu) if
212  // the current point is above the curve w:=xy, otherwise it is the
213  // other diagonal, i.e. the set of points between (xl,yu) and
214  // (xu,yl).
215  //
216  // In the two cases, we have the point
217  //
218  // (above) P(t) = (xp,yp) := (xl + t (xu-xl), yl + t (yu-yl))
219  // (below) P(t) = (xp,yp) := (xu + t (xl-xu), yl + t (yu-yl))
220  //
221  // with t in [0,1], which forms the following second degree
222  // polynomial when multiplying the coordinates:
223  //
224  // (above) f(t) = xp*yp = (yu-yl)*(xu-xl)*t^2 + [yl(xu-xl) + xl(yu-yl)]*t + xl*yl
225  // (below) f(t) = xp*yp = (yu-yl)*(xl-xu)*t^2 + [yl(xl-xu) + xu(yu-yl)]*t + xu*yl
226  //
227  // which is a quadratic function that can be expressed in the form
228  // of f'(z) = z^2 + c if we apply an affine transformation to t:
229  //
230  // t = mz + q
231  //
232  // such that the resulting coefficients of the quadratic- and the
233  // linear terms are one and zero, respectively. Thus:
234  //
235  // (above) f(z) = (yu-yl)*(xu-xl)*(mz+q)^2               + [yl(xu-xl)-xl(yu-yl)]*(mz+q) + xl*yl =
236  //              = (yu-yl)*(xu-xl)*(m^2 z^2 + 2qmz + q^2) + [yl(xu-xl)-xl(yu-yl)]*(mz+q) + xl*yl =
237  //              = z^2 + c
238  //
239  // (below) f(z) = (yu-yl)*(xl-xu)*(mz+q)^2               + [yl(xl-xu)-xu(yu-yl)]*(mz+q) + xu*yl =
240  //              = (yu-yl)*(xl-xu)*(m^2 z^2 + 2qmz + q^2) + [yl(xl-xu)-xu(yu-yl)]*(mz+q) + xu*yl =
241  //              = -z^2 + c
242  //
243  // if and only if
244  //
245  // (above) ((yu-yl)*(xu-xl)) * m^2 =  1   )
246  //                                        } <====>  m = 1 / sqrt ((yu-yl)*(xu-xl))
247  // (below) ((yu-yl)*(xl-xu)) * m^2 = -1   )
248  //
249  // (above)  2qm*(yu-yl)*(xu-xl) + m[yl(xu-xl)-xl(yu-yl)] = 0   <====>
250  //          q = -[yl(xu-xl)-xl(yu-yl)] / (2(yu-yl)*(xu-xl))
251  //
252  // (below)  2qm*(yu-yl)*(xl-xu) + m[yl(xl-xu)-xu(yu-yl)] = 0   <====>
253  //          q = -[yl(xl-xu)-xu(yu-yl)] / (2(yu-yl)*(xl-xu))
254  //
255  // If the point is below the curve, a very similar reasoning applies
256  // (simply swap xl with xu).
257  //
258  // Hence, we simply apply the balanced strategy to the function
259  // f(z)=z^2 with z bounded between -q/m and (1-q)/m. The returning
260  // value z_opt must be transformed to get
261  //
262  // t_opt = m z_opt + q
263  //
264  // and the branching point is xl + t_opt (xu-xl)
265  //
266
267  powertriplet ft (2);
268
269  // A: above
270  // B: below
271
272  bool above = (w0 > x0*y0);
273
274  CouNumber
275    dx     = xu-xl,
276    dy     = yu-yl,
277    area   = dx*dy,          // coefficient of t^2
278    bA     =  yl*dx - xl*dy, // coefficient of t
279    bB     = -yl*dx - xu*dy, // coefficient of t
280    m      = 1. / sqrt (area),
281    qA     = -bA / (2*area),
282    qB     =  bB / (2*area),
283    z_opt = above ? 
284      minMaxDelta (&ft, -qA/m, (1-qA)/m): 
285      minMaxDelta (&ft, -qB/m, (1-qB)/m),
286    t_optA = m*z_opt + qA,
287    t_optB = m*z_opt + qB;
288
289  /*printf ("------------------\n(%d,%d): [%g,%g], [%g,%g]\n", index, other, xl, xu, yl, yu);
290  printf ("dx = %g, dy = %g, area = %g, bA = %g, bB = %g, m = %g\n", dx, dy, area, bA, bB, m);
291  printf ("qA = %g, qB = %g, z_opt = %g, tA = %g, tB = %g\n",
292          qA, qB, z_opt, t_optA, t_optB);
293  printf ("w = %g %c xy = %g*%g = %g --> %g\n", w0,
294          (w0 > x0*y0) ? '>' : '<', x0, y0, x0*y0,
295          (w0 > x0*y0) ? (xl + dx * t_optA) : (xu - dx * t_optB));*/
296
297  return (w0 > x0*y0) ? 
298    (xl + dx * t_optA) : 
299    (xu - dx * t_optB);
300}
Note: See TracBrowser for help on using the repository browser.