GRASS Programmer's Manual  6.4.2(2012)
brent.c
Go to the documentation of this file.
00001 /* min/brent.c
00002  * 
00003  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
00004  *
00005  * Taken from 'GSL - The GNU Scientific Library':
00006  *             "One dimensional Minimization"
00007  *             http://sources.redhat.com/gsl/
00008  * modified by Stefano Menegon 2004
00009  *
00010  * This program is free software; you can redistribute it and/or modify
00011  * it under the terms of the GNU General Public License as published by
00012  * the Free Software Foundation; either version 2 of the License, or (at
00013  * your option) any later version.
00014  * 
00015  * This program is distributed in the hope that it will be useful, but
00016  * WITHOUT ANY WARRANTY; without even the implied warranty of
00017  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018  * General Public License for more details.
00019  * 
00020  * You should have received a copy of the GNU General Public License
00021  * along with this program; if not, write to the Free Software
00022  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
00023  */
00024 
00025 #include <stdlib.h>
00026 #include <stdio.h>
00027 #include <math.h>
00028 #include <float.h>
00029 
00030 #define GSL_SQRT_DBL_EPSILON   1.e-4
00031 #define GSL_DBL_EPSILON  1.e-8
00032 
00033 /*
00034    #define SAFE_FUNC_CALL(f, x, yp) \
00035    do { \
00036    *yp = GSL_FN_EVAL(f,x); \
00037    if (!finite(*yp)) \
00038    fprintf(stderr,"function not continuous\n");\
00039    } while (0)
00040  */
00041 
00042 
00043 typedef struct
00044 {
00045     double d, e, v, w;
00046     double f_v, f_w;
00047 }
00048 brent_state_t;
00049 
00050 
00051 
00052 static int
00053 brent(void *vstate, double (*f) (), double *x_minimum, double *f_minimum,
00054       double *x_lower, double *f_lower, double *x_upper, double *f_upper)
00055 {
00056     brent_state_t *state = (brent_state_t *) vstate;
00057 
00058     const double x_left = *x_lower;
00059     const double x_right = *x_upper;
00060 
00061     const double z = *x_minimum;
00062     double d = state->e;
00063     double e = state->d;
00064     double u, f_u;
00065     const double v = state->v;
00066     const double w = state->w;
00067     const double f_v = state->f_v;
00068     const double f_w = state->f_w;
00069     const double f_z = *f_minimum;
00070 
00071     const double golden = 0.3819660;    /* golden = (3 - sqrt(5))/2 */
00072 
00073     const double w_lower = (z - x_left);
00074     const double w_upper = (x_right - z);
00075 
00076     const double tolerance = GSL_SQRT_DBL_EPSILON * fabs(z);
00077 
00078     double p = 0, q = 0, r = 0;
00079 
00080     const double midpoint = 0.5 * (x_left + x_right);
00081 
00082     if (fabs(e) > tolerance) {
00083         /* fit parabola */
00084 
00085         r = (z - w) * (f_z - f_v);
00086         q = (z - v) * (f_z - f_w);
00087         p = (z - v) * q - (z - w) * r;
00088         q = 2 * (q - r);
00089 
00090         if (q > 0) {
00091             p = -p;
00092         }
00093         else {
00094             q = -q;
00095         }
00096 
00097         r = e;
00098         e = d;
00099     }
00100 
00101     if (fabs(p) < fabs(0.5 * q * r) && p < q * w_lower && p < q * w_upper) {
00102         double t2 = 2 * tolerance;
00103 
00104         d = p / q;
00105         u = z + d;
00106 
00107         if ((u - x_left) < t2 || (x_right - z) < t2) {
00108             d = (z < midpoint) ? tolerance : -tolerance;
00109         }
00110     }
00111     else {
00112         e = (z < midpoint) ? x_right - z : -(z - x_left);
00113         d = golden * e;
00114     }
00115 
00116 
00117     if (fabs(d) >= tolerance) {
00118         u = z + d;
00119     }
00120     else {
00121         u = z + ((d > 0) ? tolerance : -tolerance);
00122     }
00123 
00124     state->e = e;
00125     state->d = d;
00126 
00127     /*  SAFE_FUNC_CALL(f, u, &f_u); */
00128     f_u = (*f) (u);
00129 
00130     if (f_u > f_z) {
00131         if (u < z) {
00132             *x_lower = u;
00133             *f_lower = f_u;
00134             return 0;
00135         }
00136         else {
00137             *x_upper = u;
00138             *f_upper = f_u;
00139             return 0;
00140         }
00141     }
00142     else if (f_u < f_z) {
00143         if (u < z) {
00144             *x_upper = z;
00145             *f_upper = f_z;
00146         }
00147         else {
00148             *x_lower = z;
00149             *f_lower = f_z;
00150         }
00151 
00152         state->v = w;
00153         state->f_v = f_w;
00154         state->w = z;
00155         state->f_w = f_z;
00156         *x_minimum = u;
00157         *f_minimum = f_u;
00158         return 0;
00159     }
00160     else if (f_u <= f_w || w == z) {
00161         state->v = w;
00162         state->f_v = f_w;
00163         state->w = u;
00164         state->f_w = f_u;
00165         return 0;
00166     }
00167     else if (f_u <= f_v || v == z || v == w) {
00168         state->v = u;
00169         state->f_v = f_u;
00170         return 0;
00171     }
00172     else {
00173         return -1;
00174     }
00175 }
00176 
00177 
00178 /* Code modified by Stefano Menegon 1st of February 2004 */
00179 
00180 double brent_iterate(double (*f) (), double x_lower, double x_upper,
00181                      int maxiter)
00182 {
00183     int i;
00184     double x_minimum = (x_upper + x_lower) / 2.;
00185     double f_minimum;
00186     double f_lower;
00187     double f_upper;
00188 
00189     /* stato iterazione */
00190     brent_state_t itstate;
00191     const double golden = 0.3819660;    /* golden = (3 - sqrt(5))/2 */
00192     double v = x_lower + golden * (x_upper - x_lower);
00193     double w = v;
00194     double f_vw;
00195 
00196     f_lower = (*f) (x_lower);
00197     f_upper = (*f) (x_upper);
00198     f_minimum = (*f) (x_minimum);
00199 
00200     itstate.v = v;
00201     itstate.w = w;
00202 
00203     itstate.d = 0.;
00204     itstate.e = 0.;
00205 
00206     /*  SAFE_FUNC_CALL (f, v, &f_vw); */
00207 
00208     f_vw = (*f) (v);
00209 
00210     itstate.f_v = f_vw;
00211     itstate.f_w = f_vw;
00212 
00213     for (i = 0; i < maxiter; i++) {
00214         brent(&itstate, f, &x_minimum, &f_minimum, &x_lower, &f_lower,
00215               &x_upper, &f_upper);
00216         if (fabs(f_upper - f_lower) < GSL_DBL_EPSILON * fabs(f_minimum)) {
00217             return x_minimum;
00218         }
00219     }
00220 
00221     return x_minimum;
00222 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines