1 /* $Id$
2 *
3 * Name: minMaxDelta.cpp
4 * Author: Pietro Belotti
5 * Purpose: general function for computing best branching point based
6 * on min max height of resulting convexifications (dychotomic
7 * search)
8 *
9 * (C) Carnegie-Mellon University, 2007-10.
10 * This file is licensed under the Eclipse Public License (EPL)
11 */
12
13 #include "CouenneObject.hpp"
14 #include "CouenneFunTriplets.hpp"
15
16 namespace Couenne {
17
18 const int maxIter = 20;
19
20 ///
curvDistance(funtriplet * ft,CouNumber lb,CouNumber ub)21 CouNumber curvDistance (funtriplet *ft, CouNumber lb, CouNumber ub) {
22
23 // Consider the function f(x) between lb and ub. The slope of the
24 // convexification on the concave side, y = alpha x + alpha0, is:
25
26 CouNumber alpha = (ft -> F (ub) - ft -> F (lb)) / (ub - lb);
27
28 // and the constant term, alpha0, is
29
30 CouNumber alpha0 = (ub * ft -> F (lb) - lb * ft -> F (ub)) / (ub - lb);
31
32 // The point at which f(.) has derivative equal to the slope is the
33 // point of maximum height w.r.t the slope. The point z where
34 // maximum of f(z) - (ax+b), where (ax+b) is the convexification
35 // line (a=slope), is such that
36 //
37 // f'(z) - alpha = 0 ==> z = (f')^{-1} (alpha)
38
39 CouNumber z = ft -> FpInv (alpha);
40
41 // The real height is computed as [f(z) - (alpha * z + alpha0)]
42 // divided by the norm sqrt (alpha^2 + 1)
43
44 return ((ft -> F (z) - (alpha * z + alpha0)) / sqrt (alpha * alpha + 1));
45 }
46
47
48 ///
minMaxDelta(funtriplet * ft,CouNumber lb,CouNumber ub)49 CouNumber minMaxDelta (funtriplet *ft, CouNumber lb, CouNumber ub) {
50
51 CouNumber
52 lbm = lb, // extremes of the interval where to look
53 ubm = ub,
54 b = 0.5 * (lbm + ubm); // starting point
55
56 for (int iter = 0; iter < maxIter; iter++) {
57
58 CouNumber distL = curvDistance (ft, lb, b), // max height at left
59 distR = curvDistance (ft, b, ub), // max height at right
60 delta = fabs (distL) - fabs (distR);
61
62 // fprintf (stderr, "%4d %10g %10g %10g %10g %10g %10g\n",
63 // iter, lbm, ubm, b, distL, distR, delta);
64
65 if (fabs (delta) < COUENNE_EPS)
66 break;
67
68 CouNumber oldb = b;
69
70 // choose a smarter b based on an estimate of the derivative of
71 // the distance function at the current point, knowing it's null
72 // at the extremes
73
74 b = 0.5 * (lbm + ubm);
75
76 if (delta > 0) ubm = oldb; // right max height is smaller, move left
77 else lbm = oldb; // and viceversa
78 }
79
80 return b;
81 //return obj -> midInterval (b, lb, ub, info);
82 }
83
84
85 ///
maxHeight(funtriplet * ft,CouNumber lb,CouNumber ub)86 CouNumber maxHeight (funtriplet *ft, CouNumber lb, CouNumber ub) {
87 /* fprintf (stderr,"slope is (%g - %g) / (%g - %g) = %g / %g = %g ----> inverse is %g\n",
88 ft -> F (ub),
89 ft -> F (lb),
90 ub, lb,
91 ft -> F (ub) - ft -> F (lb),
92 (ub - lb),
93 (ft -> F (ub) - ft -> F (lb)) / (ub - lb),
94 ft -> FpInv ((ft -> F (ub) - ft -> F (lb)) / (ub - lb)));*/
95 return (ft -> FpInv ((ft -> F (ub) - ft -> F (lb)) / (ub - lb)));
96 }
97
98 }
99