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

Go to the documentation of this file.
00001 /* $Id: minMaxDelta.cpp 141 2009-06-03 04:19:19Z 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-08.
00010  * This file is licensed under the Common Public License (CPL)
00011  */
00012 
00013 #include "CouenneObject.hpp"
00014 #include "funtriplets.hpp"
00015 
00016 const int maxIter = 20;
00017 
00019 CouNumber curvDistance (funtriplet *ft, CouNumber lb, CouNumber ub) {
00020 
00021   // Consider the function f(x) between lb and ub. The slope of the
00022   // convexification on the concave side, y = alpha x + alpha0, is:
00023 
00024   CouNumber alpha = (ft -> F (ub) - ft -> F (lb)) / (ub - lb);
00025 
00026   // and the constant term, alpha0, is
00027 
00028   CouNumber alpha0 = (ub * ft -> F (lb) - lb * ft -> F (ub)) / (ub - lb);
00029 
00030   // The point at which f(.) has derivative equal to the slope is the
00031   // point of maximum height w.r.t the slope. The point z where
00032   // maximum of f(z) - (ax+b), where (ax+b) is the convexification
00033   // line (a=slope), is such that 
00034   //
00035   // f'(z) - alpha = 0   ==> z = (f')^{-1} (alpha)
00036 
00037   CouNumber z = ft -> FpInv (alpha);
00038 
00039   // The real height is computed as [f(z) - (alpha * z + alpha0)]
00040   // divided by the norm sqrt (alpha^2 + 1)
00041 
00042   return ((ft -> F (z) - (alpha * z + alpha0)) / sqrt (alpha * alpha + 1));
00043 }
00044 
00045 
00047 CouNumber minMaxDelta (funtriplet *ft, CouNumber lb, CouNumber ub) {
00048 
00049   CouNumber 
00050     lbm = lb,                // extremes of the interval where to look 
00051     ubm = ub,     
00052     b   = 0.5 * (lbm + ubm); // starting point
00053 
00054   for (int iter = 0; iter < maxIter; iter++) {
00055 
00056     CouNumber distL = curvDistance (ft, lb,  b),  // max height at left
00057               distR = curvDistance (ft,  b, ub),  // max height at right
00058               delta = fabs (distL) - fabs (distR);
00059 
00060     //    fprintf (stderr, "%4d %10g %10g %10g %10g %10g %10g\n", 
00061     //       iter, lbm, ubm, b, distL, distR, delta);
00062 
00063     if (fabs (delta) < COUENNE_EPS) 
00064       break;
00065 
00066     CouNumber oldb = b;
00067 
00068     // choose a smarter b based on an estimate of the derivative of
00069     // the distance function at the current point, knowing it's null
00070     // at the extremes
00071 
00072     b = 0.5 * (lbm + ubm);
00073 
00074     if (delta > 0) ubm = oldb; // right max height is smaller, move left
00075     else           lbm = oldb; // and viceversa
00076   }
00077 
00078   return b;
00079   //return obj -> midInterval (b, lb, ub);
00080 }
00081 
00082 
00084 CouNumber maxHeight (funtriplet *ft, CouNumber lb, CouNumber ub) {
00085   /* fprintf (stderr,"slope is (%g - %g) / (%g - %g) = %g / %g = %g ----> inverse is %g\n", 
00086           ft -> F (ub), 
00087           ft -> F (lb), 
00088           ub, lb,
00089           ft -> F (ub) - ft -> F (lb),
00090           (ub - lb),
00091           (ft -> F (ub) - ft -> F (lb)) / (ub - lb),
00092           ft -> FpInv ((ft -> F (ub) - ft -> F (lb)) / (ub - lb)));*/
00093   return (ft -> FpInv ((ft -> F (ub) - ft -> F (lb)) / (ub - lb)));
00094 }

Generated on Thu Oct 8 03:02:56 2009 by  doxygen 1.4.7