GRASS Programmer's Manual
6.4.2(2012)
|
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 }