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

   1 /* multimin/directional_minimize.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 static void
  21 take_step (const gsl_vector * x, const gsl_vector * p,
  22            double step, double lambda, gsl_vector * x1, gsl_vector * dx)
  23 {
  24   gsl_vector_set_zero (dx);
  25   gsl_blas_daxpy (-step * lambda, p, dx);
  26
  27   gsl_vector_memcpy (x1, x);
  28   gsl_blas_daxpy (1.0, dx, x1);
  29 }
  30
  31 static void
  32 intermediate_point (gsl_multimin_function_fdf * fdf,
  33                     const gsl_vector * x, const gsl_vector * p,
  34                     double lambda,
  35                     double pg,
  36                     double stepa, double stepc,
  37                     double fa, double fc,
  38                     gsl_vector * x1, gsl_vector * dx, gsl_vector * gradient,
  39                     double * step, double * f)
  40 {
  41   double stepb, fb;
  42
  43 trial:
  44   {
  45     double u = fabs (pg * lambda * stepc);
  46     stepb = 0.5 * stepc * u / ((fc - fa) + u);
  47   }
  48
  49   take_step (x, p, stepb, lambda, x1, dx);
  50
  51   fb = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);
  52
  53 #ifdef DEBUG
  54   printf ("trying stepb = %g  fb = %.18e\n", stepb, fb);
  55 #endif
  56
  57   if (fb >= fa  && stepb > 0.0)
  58     {
  59       /* downhill step failed, reduce step-size and try again */
  60       fc = fb;
  61       stepc = stepb;
  62       goto trial;
  63     }
  64 #ifdef DEBUG
  65   printf ("ok!\n");
  66 #endif
  67
  68   *step = stepb;
  69   *f = fb;
  70   GSL_MULTIMIN_FN_EVAL_DF(fdf, x1, gradient);
  71 }
  72
  73 static void
  74 minimize (gsl_multimin_function_fdf * fdf,
  75           const gsl_vector * x, const gsl_vector * p,
  76           double lambda,
  77           double stepa, double stepb, double stepc,
  78           double fa, double fb, double fc, double tol,
  79           gsl_vector * x1, gsl_vector * dx1,
  80           gsl_vector * x2, gsl_vector * dx2, gsl_vector * gradient,
  81           double * step, double * f, double * gnorm)
  82 {
  83   /* Starting at (x0, f0) move along the direction p to find a minimum
  84      f(x0 - lambda * p), returning the new point x1 = x0-lambda*p,
  85      f1=f(x1) and g1 = grad(f) at x1.  */
  86
  87   double u = stepb;
  88   double v = stepa;
  89   double w = stepc;
  90   double fu = fb;
  91   double fv = fa;
  92   double fw = fc;
  93
  94   double old2 = fabs(w - v);
  95   double old1 = fabs(v - u);
  96
  97   double stepm, fm, pg, gnorm1;
  98
  99   int iter = 0;
 100
 101   gsl_vector_memcpy (x2, x1);
 102   gsl_vector_memcpy (dx2, dx1);
 103
 104   *f = fb;
 105   *step = stepb;
 106   *gnorm = gsl_blas_dnrm2 (gradient);
 107
 108 mid_trial:
 109
 110   iter++;
 111
 112   if (iter > 10)
 113     {
 114       return;  /* MAX ITERATIONS */
 115     }
 116
 117   {
 118     double dw = w - u;
 119     double dv = v - u;
 120     double du = 0.0;
 121
 122     double e1 = ((fv - fu) * dw * dw + (fu - fw) * dv * dv);
 123     double e2 = 2.0 * ((fv - fu) * dw + (fu - fw) * dv);
 124
 125     if (e2 != 0.0)
 126       {
 127         du = e1 / e2;
 128       }
 129
 130     if (du > 0.0 && du < (stepc - stepb) && fabs(du) < 0.5 * old2)
 131       {
 132         stepm = u + du;
 133       }
 134     else if (du < 0.0 && du > (stepa - stepb) && fabs(du) < 0.5 * old2)
 135       {
 136         stepm = u + du;
 137       }
 138     else if ((stepc - stepb) > (stepb - stepa))
 139       {
 140         stepm = 0.38 * (stepc - stepb) + stepb;
 141       }
 142     else
 143       {
 144         stepm = stepb - 0.38 * (stepb - stepa);
 145       }
 146   }
 147
 148   take_step (x, p, stepm, lambda, x1, dx1);
 149
 150   fm = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);
 151
 152 #ifdef DEBUG
 153   printf ("trying stepm = %g  fm = %.18e\n", stepm, fm);
 154 #endif
 155
 156   if (fm > fb)
 157     {
 158       if (fm < fv)
 159         {
 160           w = v;
 161           v = stepm;
 162           fw = fv;
 163           fv = fm;
 164         }
 165       else if (fm < fw)
 166         {
 167           w = stepm;
 168           fw = fm;
 169         }
 170
 171       if (stepm < stepb)
 172         {
 173           stepa = stepm;
 174           fa = fm;
 175         }
 176       else
 177         {
 178           stepc = stepm;
 179           fc = fm;
 180         }
 181       goto mid_trial;
 182     }
 183   else if (fm <= fb)
 184     {
 185       old2 = old1;
 186       old1 = fabs(u - stepm);
 187       w = v;
 188       v = u;
 189       u = stepm;
 190       fw = fv;
 191       fv = fu;
 192       fu = fm;
 193
 194       gsl_vector_memcpy (x2, x1);
 195       gsl_vector_memcpy (dx2, dx1);
 196
 197       GSL_MULTIMIN_FN_EVAL_DF (fdf, x1, gradient);
 198       gsl_blas_ddot (p, gradient, &pg);
 199       gnorm1 = gsl_blas_dnrm2 (gradient);
 200
 201 #ifdef DEBUG
 202       printf ("p: "); gsl_vector_fprintf(stdout, p, "%g");
 203       printf ("g: "); gsl_vector_fprintf(stdout, gradient, "%g");
 204       printf ("gnorm: %.18e\n", gnorm1);
 205       printf ("pg: %.18e\n", pg);
 206       printf ("orth: %g\n", fabs (pg * lambda/ gnorm1));
 207 #endif
 208       *f = fm;
 209       *step = stepm;
 210       *gnorm = gnorm1;
 211
 212       if (fabs (pg * lambda / gnorm1) < tol)
 213         {
 214 #ifdef DEBUG
 215           printf("ok!\n");
 216 #endif
 217           return; /* SUCCESS */
 218         }
 219
 220       if (stepm < stepb)
 221         {
 222           stepc = stepb;
 223           fc = fb;
 224           stepb = stepm;
 225           fb = fm;
 226         }
 227       else
 228         {
 229           stepa = stepb;
 230           fa = fb;
 231           stepb = stepm;
 232           fb = fm;
 233         }
 234       goto mid_trial;
 235     }
 236 }