GRASS GIS 8 Programmer's Manual 8.2.0(2022)-exported
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
43typedef struct
44{
45 double d, e, v, w;
46 double f_v, f_w;
47}
48brent_state_t;
49
50
51
52static int
53brent(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
180double 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
double r
struct state state
Definition: parser.c:103