/home/coin/SVN-release/OS-2.4.0/Couenne/src/branch/operators/minMaxDelta.cpp

Go to the documentation of this file.
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 }

Generated on Thu Sep 22 03:05:56 2011 by  doxygen 1.4.7