htswanalysis/MACS/lib/gsl/gsl-1.11/min/brent.c

   1 /* min/brent.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   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 /* brent.c -- brent minimum finding algorithm */
  21
  22 #include <config.h>
  23
  24 #include <stddef.h>
  25 #include <stdlib.h>
  26 #include <stdio.h>
  27 #include <math.h>
  28 #include <float.h>
  29
  30 #include <gsl/gsl_math.h>
  31 #include <gsl/gsl_errno.h>
  32 #include <gsl/gsl_min.h>
  33
  34 #include "min.h"
  35
  36 typedef struct
  37   {
  38     double d, e, v, w;
  39     double f_v, f_w;
  40   }
  41 brent_state_t;
  42
  43 static int brent_init (void *vstate, gsl_function * f, double x_minimum, double f_minimum, double x_lower, double f_lower, double x_upper, double f_upper);
  44 static int brent_iterate (void *vstate, gsl_function * f, double *x_minimum, double * f_minimum, double * x_lower, double * f_lower, double * x_upper, double * f_upper);
  45
  46 static int
  47 brent_init (void *vstate, gsl_function * f, double x_minimum, double f_minimum, double x_lower, double f_lower, double x_upper, double f_upper)
  48 {
  49   brent_state_t *state = (brent_state_t *) vstate;
  50
  51   const double golden = 0.3819660;      /* golden = (3 - sqrt(5))/2 */
  52
  53   double v = x_lower + golden * (x_upper - x_lower);
  54   double w = v;
  55
  56   double f_vw;
  57
  58   x_minimum = 0 ;  /* avoid warnings about unused varibles */
  59   f_minimum = 0 ;
  60   f_lower = 0 ;
  61   f_upper = 0 ;
  62
  63   state->v = v;
  64   state->w = w;
  65
  66   state->d = 0;
  67   state->e = 0;
  68
  69   SAFE_FUNC_CALL (f, v, &f_vw);
  70
  71   state->f_v = f_vw;
  72   state->f_w = f_vw;
  73
  74   return GSL_SUCCESS;
  75 }
  76
  77 static int
  78 brent_iterate (void *vstate, gsl_function * f, double *x_minimum, double * f_minimum, double * x_lower, double * f_lower, double * x_upper, double * f_upper)
  79 {
  80   brent_state_t *state = (brent_state_t *) vstate;
  81
  82   const double x_left = *x_lower;
  83   const double x_right = *x_upper;
  84
  85   const double z = *x_minimum;
  86   double d = state->e;
  87   double e = state->d;
  88   double u, f_u;
  89   const double v = state->v;
  90   const double w = state->w;
  91   const double f_v = state->f_v;
  92   const double f_w = state->f_w;
  93   const double f_z = *f_minimum;
  94
  95   const double golden = 0.3819660;      /* golden = (3 - sqrt(5))/2 */
  96
  97   const double w_lower = (z - x_left);
  98   const double w_upper = (x_right - z);
  99
 100   const double tolerance =  GSL_SQRT_DBL_EPSILON * fabs (z);
 101
 102   double p = 0, q = 0, r = 0;
 103
 104   const double midpoint = 0.5 * (x_left + x_right);
 105
 106   if (fabs (e) > tolerance)
 107     {
 108       /* fit parabola */
 109
 110       r = (z - w) * (f_z - f_v);
 111       q = (z - v) * (f_z - f_w);
 112       p = (z - v) * q - (z - w) * r;
 113       q = 2 * (q - r);
 114
 115       if (q > 0)
 116         {
 117           p = -p;
 118         }
 119       else
 120         {
 121           q = -q;
 122         }
 123
 124       r = e;
 125       e = d;
 126     }
 127
 128   if (fabs (p) < fabs (0.5 * q * r) && p < q * w_lower && p < q * w_upper)
 129     {
 130       double t2 = 2 * tolerance ;
 131
 132       d = p / q;
 133       u = z + d;
 134
 135       if ((u - x_left) < t2 || (x_right - u) < t2)
 136         {
 137           d = (z < midpoint) ? tolerance : -tolerance ;
 138         }
 139     }
 140   else
 141     {
 142       e = (z < midpoint) ? x_right - z : -(z - x_left) ;
 143       d = golden * e;
 144     }
 145
 146
 147   if (fabs (d) >= tolerance)
 148     {
 149       u = z + d;
 150     }
 151   else
 152     {
 153       u = z + ((d > 0) ? tolerance : -tolerance) ;
 154     }
 155
 156   state->e = e;
 157   state->d = d;
 158
 159   SAFE_FUNC_CALL(f, u, &f_u);
 160
 161   if (f_u <= f_z)
 162     {
 163       if (u < z)
 164         {
 165           *x_upper = z;
 166           *f_upper = f_z;
 167         }
 168       else
 169         {
 170           *x_lower = z;
 171           *f_lower = f_z;
 172         }
 173
 174       state->v = w;
 175       state->f_v = f_w;
 176       state->w = z;
 177       state->f_w = f_z;
 178       *x_minimum = u;
 179       *f_minimum = f_u;
 180       return GSL_SUCCESS;
 181     }
 182   else
 183     {
 184       if (u < z)
 185         {
 186           *x_lower = u;
 187           *f_lower = f_u;
 188         }
 189       else
 190         {
 191           *x_upper = u;
 192           *f_upper = f_u;
 193         }
 194
 195       if (f_u <= f_w || w == z)
 196         {
 197           state->v = w;
 198           state->f_v = f_w;
 199           state->w = u;
 200           state->f_w = f_u;
 201           return GSL_SUCCESS;
 202         }
 203       else if (f_u <= f_v || v == z || v == w)
 204         {
 205           state->v = u;
 206           state->f_v = f_u;
 207           return GSL_SUCCESS;
 208         }
 209     }
 210
 211   return GSL_SUCCESS;
 212 }
 213
 214
 215 static const gsl_min_fminimizer_type brent_type =
 216 {"brent",                       /* name */
 217  sizeof (brent_state_t),
 218  &brent_init,
 219  &brent_iterate};
 220
 221 const gsl_min_fminimizer_type *gsl_min_fminimizer_brent = &brent_type;