00001 /* $Id: minMaxDelta.cpp 708 2011-06-23 14:04:59Z pbelotti $ 00002 * 00003 * Name: minMaxDelta.cpp 00004 * Author: Pietro Belotti 00005 * Purpose: general function for computing best branching point based 00006 * on min max height of resulting convexifications (dychotomic 00007 * search) 00008 * 00009 * (C) Carnegie-Mellon University, 2007-10. 00010 * This file is licensed under the Eclipse Public License (EPL) 00011 */ 00012 00013 #include "CouenneObject.hpp" 00014 #include "CouenneFunTriplets.hpp" 00015 00016 namespace Couenne { 00017 00018 const int maxIter = 20; 00019 00021 CouNumber curvDistance (funtriplet *ft, CouNumber lb, CouNumber ub) { 00022 00023 // Consider the function f(x) between lb and ub. The slope of the 00024 // convexification on the concave side, y = alpha x + alpha0, is: 00025 00026 CouNumber alpha = (ft -> F (ub) - ft -> F (lb)) / (ub - lb); 00027 00028 // and the constant term, alpha0, is 00029 00030 CouNumber alpha0 = (ub * ft -> F (lb) - lb * ft -> F (ub)) / (ub - lb); 00031 00032 // The point at which f(.) has derivative equal to the slope is the 00033 // point of maximum height w.r.t the slope. The point z where 00034 // maximum of f(z) - (ax+b), where (ax+b) is the convexification 00035 // line (a=slope), is such that 00036 // 00037 // f'(z) - alpha = 0 ==> z = (f')^{-1} (alpha) 00038 00039 CouNumber z = ft -> FpInv (alpha); 00040 00041 // The real height is computed as [f(z) - (alpha * z + alpha0)] 00042 // divided by the norm sqrt (alpha^2 + 1) 00043 00044 return ((ft -> F (z) - (alpha * z + alpha0)) / sqrt (alpha * alpha + 1)); 00045 } 00046 00047 00049 CouNumber minMaxDelta (funtriplet *ft, CouNumber lb, CouNumber ub) { 00050 00051 CouNumber 00052 lbm = lb, // extremes of the interval where to look 00053 ubm = ub, 00054 b = 0.5 * (lbm + ubm); // starting point 00055 00056 for (int iter = 0; iter < maxIter; iter++) { 00057 00058 CouNumber distL = curvDistance (ft, lb, b), // max height at left 00059 distR = curvDistance (ft, b, ub), // max height at right 00060 delta = fabs (distL) - fabs (distR); 00061 00062 // fprintf (stderr, "%4d %10g %10g %10g %10g %10g %10g\n", 00063 // iter, lbm, ubm, b, distL, distR, delta); 00064 00065 if (fabs (delta) < COUENNE_EPS) 00066 break; 00067 00068 CouNumber oldb = b; 00069 00070 // choose a smarter b based on an estimate of the derivative of 00071 // the distance function at the current point, knowing it's null 00072 // at the extremes 00073 00074 b = 0.5 * (lbm + ubm); 00075 00076 if (delta > 0) ubm = oldb; // right max height is smaller, move left 00077 else lbm = oldb; // and viceversa 00078 } 00079 00080 return b; 00081 //return obj -> midInterval (b, lb, ub, info); 00082 } 00083 00084 00086 CouNumber maxHeight (funtriplet *ft, CouNumber lb, CouNumber ub) { 00087 /* fprintf (stderr,"slope is (%g - %g) / (%g - %g) = %g / %g = %g ----> inverse is %g\n", 00088 ft -> F (ub), 00089 ft -> F (lb), 00090 ub, lb, 00091 ft -> F (ub) - ft -> F (lb), 00092 (ub - lb), 00093 (ft -> F (ub) - ft -> F (lb)) / (ub - lb), 00094 ft -> FpInv ((ft -> F (ub) - ft -> F (lb)) / (ub - lb)));*/ 00095 return (ft -> FpInv ((ft -> F (ub) - ft -> F (lb)) / (ub - lb))); 00096 } 00097 00098 }