GRASS GIS 7 Programmer's Manual  7.5.svn(2017)-r71915
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
brent.c
Go to the documentation of this file.
1 /* min/brent.c
2  *
3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
4  *
5  * Taken from 'GSL - The GNU Scientific Library':
6  * "One dimensional Minimization"
7  * http://sources.redhat.com/gsl/
8  * modified by Stefano Menegon 2004
9  *
10  * This program is free software; you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation; either version 2 of the License, or (at
13  * your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful, but
16  * WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  * General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with this program; if not, write to the Free Software
22  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23  */
24 
25 #include <stdlib.h>
26 #include <stdio.h>
27 #include <math.h>
28 #include <float.h>
29 
30 #define GSL_SQRT_DBL_EPSILON 1.e-4
31 #define GSL_DBL_EPSILON 1.e-8
32 
33 /*
34  #define SAFE_FUNC_CALL(f, x, yp) \
35  do { \
36  *yp = GSL_FN_EVAL(f,x); \
37  if (!finite(*yp)) \
38  fprintf(stderr,"function not continuous\n");\
39  } while (0)
40  */
41 
42 
43 typedef struct
44 {
45  double d, e, v, w;
46  double f_v, f_w;
47 }
48 brent_state_t;
49 
50 
51 
52 static int
53 brent(void *vstate, double (*f) (), double *x_minimum, double *f_minimum,
54  double *x_lower, double *f_lower, double *x_upper, double *f_upper)
55 {
56  brent_state_t *state = (brent_state_t *) vstate;
57 
58  const double x_left = *x_lower;
59  const double x_right = *x_upper;
60 
61  const double z = *x_minimum;
62  double d = state->e;
63  double e = state->d;
64  double u, f_u;
65  const double v = state->v;
66  const double w = state->w;
67  const double f_v = state->f_v;
68  const double f_w = state->f_w;
69  const double f_z = *f_minimum;
70 
71  const double golden = 0.3819660; /* golden = (3 - sqrt(5))/2 */
72 
73  const double w_lower = (z - x_left);
74  const double w_upper = (x_right - z);
75 
76  const double tolerance = GSL_SQRT_DBL_EPSILON * fabs(z);
77 
78  double p = 0, q = 0, r = 0;
79 
80  const double midpoint = 0.5 * (x_left + x_right);
81 
82  if (fabs(e) > tolerance) {
83  /* fit parabola */
84 
85  r = (z - w) * (f_z - f_v);
86  q = (z - v) * (f_z - f_w);
87  p = (z - v) * q - (z - w) * r;
88  q = 2 * (q - r);
89 
90  if (q > 0) {
91  p = -p;
92  }
93  else {
94  q = -q;
95  }
96 
97  r = e;
98  e = d;
99  }
100 
101  if (fabs(p) < fabs(0.5 * q * r) && p < q * w_lower && p < q * w_upper) {
102  double t2 = 2 * tolerance;
103 
104  d = p / q;
105  u = z + d;
106 
107  if ((u - x_left) < t2 || (x_right - z) < t2) {
108  d = (z < midpoint) ? tolerance : -tolerance;
109  }
110  }
111  else {
112  e = (z < midpoint) ? x_right - z : -(z - x_left);
113  d = golden * e;
114  }
115 
116 
117  if (fabs(d) >= tolerance) {
118  u = z + d;
119  }
120  else {
121  u = z + ((d > 0) ? tolerance : -tolerance);
122  }
123 
124  state->e = e;
125  state->d = d;
126 
127  /* SAFE_FUNC_CALL(f, u, &f_u); */
128  f_u = (*f) (u);
129 
130  if (f_u > f_z) {
131  if (u < z) {
132  *x_lower = u;
133  *f_lower = f_u;
134  return 0;
135  }
136  else {
137  *x_upper = u;
138  *f_upper = f_u;
139  return 0;
140  }
141  }
142  else if (f_u < f_z) {
143  if (u < z) {
144  *x_upper = z;
145  *f_upper = f_z;
146  }
147  else {
148  *x_lower = z;
149  *f_lower = f_z;
150  }
151 
152  state->v = w;
153  state->f_v = f_w;
154  state->w = z;
155  state->f_w = f_z;
156  *x_minimum = u;
157  *f_minimum = f_u;
158  return 0;
159  }
160  else if (f_u <= f_w || w == z) {
161  state->v = w;
162  state->f_v = f_w;
163  state->w = u;
164  state->f_w = f_u;
165  return 0;
166  }
167  else if (f_u <= f_v || v == z || v == w) {
168  state->v = u;
169  state->f_v = f_u;
170  return 0;
171  }
172  else {
173  return -1;
174  }
175 }
176 
177 
178 /* Code modified by Stefano Menegon 1st of February 2004 */
179 
180 double brent_iterate(double (*f) (), double x_lower, double x_upper,
181  int maxiter)
182 {
183  int i;
184  double x_minimum = (x_upper + x_lower) / 2.;
185  double f_minimum;
186  double f_lower;
187  double f_upper;
188 
189  /* stato iterazione */
190  brent_state_t itstate;
191  const double golden = 0.3819660; /* golden = (3 - sqrt(5))/2 */
192  double v = x_lower + golden * (x_upper - x_lower);
193  double w = v;
194  double f_vw;
195 
196  f_lower = (*f) (x_lower);
197  f_upper = (*f) (x_upper);
198  f_minimum = (*f) (x_minimum);
199 
200  itstate.v = v;
201  itstate.w = w;
202 
203  itstate.d = 0.;
204  itstate.e = 0.;
205 
206  /* SAFE_FUNC_CALL (f, v, &f_vw); */
207 
208  f_vw = (*f) (v);
209 
210  itstate.f_v = f_vw;
211  itstate.f_w = f_vw;
212 
213  for (i = 0; i < maxiter; i++) {
214  brent(&itstate, f, &x_minimum, &f_minimum, &x_lower, &f_lower,
215  &x_upper, &f_upper);
216  if (fabs(f_upper - f_lower) < GSL_DBL_EPSILON * fabs(f_minimum)) {
217  return x_minimum;
218  }
219  }
220 
221  return x_minimum;
222 }
#define GSL_SQRT_DBL_EPSILON
Definition: brent.c:30
#define GSL_DBL_EPSILON
Definition: brent.c:31
double brent_iterate(double(*f)(), double x_lower, double x_upper, int maxiter)
Definition: brent.c:180
struct state state
Definition: parser.c:102
double r
Definition: r_raster.c:39