htswanalysis/MACS/lib/gsl/gsl-1.11/multimin/linear_minimize.c

   1 #include <gsl/gsl_math.h>
   2 #include <gsl/gsl_errno.h>
   3 #include <gsl/gsl_poly.h>
   4
   5 /* Find a minimum in x=[0,1] of the interpolating quadratic through
   6  * (0,f0) (1,f1) with derivative fp0 at x=0.  The interpolating
   7  * polynomial is q(x) = f0 + fp0 * z + (f1-f0-fp0) * z^2
   8  */
   9
  10 static double
  11 interp_quad (double f0, double fp0, double f1, double zl, double zh)
  12 {
  13   double fl = f0 + zl*(fp0 + zl*(f1 - f0 -fp0));
  14   double fh = f0 + zh*(fp0 + zh*(f1 - f0 -fp0));
  15   double c = 2 * (f1 - f0 - fp0);       /* curvature */
  16
  17   double zmin = zl, fmin = fl;
  18
  19   if (fh < fmin) { zmin = zh; fmin = fh; }
  20
  21   if (c > 0)  /* positive curvature required for a minimum */
  22     {
  23       double z = -fp0 / c;      /* location of minimum */
  24       if (z > zl && z < zh) {
  25         double f = f0 + z*(fp0 + z*(f1 - f0 -fp0));
  26         if (f < fmin) { zmin = z; fmin = f; };
  27       }
  28     }
  29
  30   return zmin;
  31 }
  32
  33 /* Find a minimum in x=[0,1] of the interpolating cubic through
  34  * (0,f0) (1,f1) with derivatives fp0 at x=0 and fp1 at x=1.
  35  *
  36  * The interpolating polynomial is:
  37  *
  38  * c(x) = f0 + fp0 * z + eta * z^2 + xi * z^3
  39  *
  40  * where eta=3*(f1-f0)-2*fp0-fp1, xi=fp0+fp1-2*(f1-f0).
  41  */
  42
  43 static double
  44 cubic (double c0, double c1, double c2, double c3, double z)
  45 {
  46   return c0 + z * (c1 + z * (c2 + z * c3));
  47 }
  48
  49 static void
  50 check_extremum (double c0, double c1, double c2, double c3, double z,
  51                 double *zmin, double *fmin)
  52 {
  53   /* could make an early return by testing curvature >0 for minimum */
  54
  55   double y = cubic (c0, c1, c2, c3, z);
  56
  57   if (y < *fmin)
  58     {
  59       *zmin = z;  /* accepted new point*/
  60       *fmin = y;
  61     }
  62 }
  63
  64 static double
  65 interp_cubic (double f0, double fp0, double f1, double fp1, double zl, double zh)
  66 {
  67   double eta = 3 * (f1 - f0) - 2 * fp0 - fp1;
  68   double xi = fp0 + fp1 - 2 * (f1 - f0);
  69   double c0 = f0, c1 = fp0, c2 = eta, c3 = xi;
  70   double zmin, fmin;
  71   double z0, z1;
  72
  73   zmin = zl; fmin = cubic(c0, c1, c2, c3, zl);
  74   check_extremum (c0, c1, c2, c3, zh, &zmin, &fmin);
  75
  76   {
  77     int n = gsl_poly_solve_quadratic (3 * c3, 2 * c2, c1, &z0, &z1);
  78
  79     if (n == 2)  /* found 2 roots */
  80       {
  81         if (z0 > zl && z0 < zh)
  82           check_extremum (c0, c1, c2, c3, z0, &zmin, &fmin);
  83         if (z1 > zl && z1 < zh)
  84           check_extremum (c0, c1, c2, c3, z1, &zmin, &fmin);
  85       }
  86     else if (n == 1)  /* found 1 root */
  87       {
  88         if (z0 > zl && z0 < zh)
  89           check_extremum (c0, c1, c2, c3, z0, &zmin, &fmin);
  90       }
  91   }
  92
  93   return zmin;
  94 }
  95
  96
  97 static double
  98 interpolate (double a, double fa, double fpa,
  99              double b, double fb, double fpb, double xmin, double xmax,
 100              int order)
 101 {
 102   /* Map [a,b] to [0,1] */
 103   double z, alpha, zmin, zmax;
 104
 105   zmin = (xmin - a) / (b - a);
 106   zmax = (xmax - a) / (b - a);
 107
 108   if (zmin > zmax)
 109     {
 110       double tmp = zmin;
 111       zmin = zmax;
 112       zmax = tmp;
 113     };
 114
 115   if (order > 2 && GSL_IS_REAL(fpb)) {
 116     z = interp_cubic (fa, fpa * (b - a), fb, fpb * (b - a), zmin, zmax);
 117   } else {
 118     z = interp_quad (fa, fpa * (b - a), fb, zmin, zmax);
 119   }
 120
 121   alpha = a + z * (b - a);
 122
 123   return alpha;
 124 }
 125
 126 /* recommended values from Fletcher are
 127    rho = 0.01, sigma = 0.1, tau1 = 9, tau2 = 0.05, tau3 = 0.5 */
 128
 129 static int
 130 minimize (gsl_function_fdf * fn, double rho, double sigma,
 131           double tau1, double tau2, double tau3,
 132           int order, double alpha1, double *alpha_new)
 133 {
 134   double f0, fp0, falpha, falpha_prev, fpalpha, fpalpha_prev, delta,
 135     alpha_next;
 136   double alpha = alpha1, alpha_prev = 0.0;
 137   double a, b, fa, fb, fpa, fpb;
 138   const size_t bracket_iters = 100, section_iters = 100;
 139   size_t i = 0;
 140
 141   GSL_FN_FDF_EVAL_F_DF (fn, 0.0, &f0, &fp0);
 142   falpha_prev = f0;
 143   fpalpha_prev = fp0;
 144
 145   /* Avoid uninitialized variables morning */
 146   a = 0.0; b = alpha;
 147   fa = f0; fb = 0.0;
 148   fpa = fp0; fpb = 0.0;
 149
 150   /* Begin bracketing */
 151
 152   while (i++ < bracket_iters)
 153     {
 154       falpha = GSL_FN_FDF_EVAL_F (fn, alpha);
 155
 156       /* Fletcher's rho test */
 157
 158       if (falpha > f0 + alpha * rho * fp0 || falpha >= falpha_prev)
 159         {
 160           a = alpha_prev; fa = falpha_prev; fpa = fpalpha_prev;
 161           b = alpha; fb = falpha; fpb = GSL_NAN;
 162           break;                /* goto sectioning */
 163         }
 164
 165       fpalpha = GSL_FN_FDF_EVAL_DF (fn, alpha);
 166
 167       /* Fletcher's sigma test */
 168
 169       if (fabs (fpalpha) <= -sigma * fp0)
 170         {
 171           *alpha_new = alpha;
 172           return GSL_SUCCESS;
 173         }
 174
 175       if (fpalpha >= 0)
 176         {
 177           a = alpha; fa = falpha; fpa = fpalpha;
 178           b = alpha_prev; fb = falpha_prev; fpb = fpalpha_prev;
 179           break;                /* goto sectioning */
 180         }
 181
 182       delta = alpha - alpha_prev;
 183
 184       {
 185         double lower = alpha + delta;
 186         double upper = alpha + tau1 * delta;
 187
 188         alpha_next = interpolate (alpha_prev, falpha_prev, fpalpha_prev,
 189                              alpha, falpha, fpalpha, lower, upper, order);
 190
 191       }
 192
 193       alpha_prev = alpha;
 194       falpha_prev = falpha;
 195       fpalpha_prev = fpalpha;
 196       alpha = alpha_next;
 197     }
 198
 199   /*  Sectioning of bracket [a,b] */
 200
 201   while (i++ < section_iters)
 202     {
 203       delta = b - a;
 204
 205       {
 206         double lower = a + tau2 * delta;
 207         double upper = b - tau3 * delta;
 208
 209         alpha = interpolate (a, fa, fpa, b, fb, fpb, lower, upper, order);
 210       }
 211
 212       falpha = GSL_FN_FDF_EVAL_F (fn, alpha);
 213
 214       if ((a-alpha)*fpa <= GSL_DBL_EPSILON) {
 215         /* roundoff prevents progress */
 216         return GSL_ENOPROG;
 217       };
 218
 219       if (falpha > f0 + rho * alpha * fp0 || falpha >= fa)
 220         {
 221           /*  a_next = a; */
 222           b = alpha; fb = falpha; fpb = GSL_NAN;
 223         }
 224       else
 225         {
 226           fpalpha = GSL_FN_FDF_EVAL_DF (fn, alpha);
 227
 228           if (fabs(fpalpha) <= -sigma * fp0)
 229             {
 230               *alpha_new = alpha;
 231               return GSL_SUCCESS;  /* terminate */
 232             }
 233
 234           if ( ((b-a) >= 0 && fpalpha >= 0) || ((b-a) <=0 && fpalpha <= 0))
 235             {
 236               b = a; fb = fa; fpb = fpa;
 237               a = alpha; fa = falpha; fpa = fpalpha;
 238             }
 239           else
 240             {
 241               a = alpha; fa = falpha; fpa = fpalpha;
 242             }
 243         }
 244     }
 245
 246   return GSL_SUCCESS;
 247 }