5 Copyright (c) 2008, 2010 by Attractive Chaos <attractor@live.co.uk>
7 Permission is hereby granted, free of charge, to any person obtaining
8 a copy of this software and associated documentation files (the
9 "Software"), to deal in the Software without restriction, including
10 without limitation the rights to use, copy, modify, merge, publish,
11 distribute, sublicense, and/or sell copies of the Software, and to
12 permit persons to whom the Software is furnished to do so, subject to
13 the following conditions:
15 The above copyright notice and this permission notice shall be
16 included in all copies or substantial portions of the Software.
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
19 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
20 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
21 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
22 BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
23 ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
24 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
28 /* Hooke-Jeeves algorithm for nonlinear minimization
30 Based on the pseudocodes by Bell and Pike (CACM 9(9):684-685), and
31 the revision by Tomlin and Smith (CACM 12(11):637-638). Both of the
32 papers are comments on Kaupe's Algorithm 178 "Direct Search" (ACM
33 6(6):313-314). The original algorithm was designed by Hooke and
34 Jeeves (ACM 8:212-229). This program is further revised according to
35 Johnson's implementation at Netlib (opt/hooke.c).
37 Hooke-Jeeves algorithm is very simple and it works quite well on a
38 few examples. However, it might fail to converge due to its heuristic
39 nature. A possible improvement, as is suggested by Johnson, may be to
40 choose a small r at the beginning to quickly approach to the minimum
41 and a large r at later step to hit the minimum.
49 static double __kmin_hj_aux(kmin_f func, int n, double *x1, void *data, double fx1, double *dx, int *n_calls)
53 for (k = 0; k != n; ++k) {
55 ftmp = func(n, x1, data); ++j;
56 if (ftmp < fx1) fx1 = ftmp;
57 else { /* search the opposite direction */
59 x1[k] += dx[k] + dx[k];
60 ftmp = func(n, x1, data); ++j;
61 if (ftmp < fx1) fx1 = ftmp;
62 else x1[k] -= dx[k]; /* back to the original x[k] */
66 return fx1; /* here: fx1=f(n,x1) */
69 double kmin_hj(kmin_f func, int n, double *x, void *data, double r, double eps, int max_calls)
71 double fx, fx1, *x1, *dx, radius;
73 x1 = (double*)calloc(n, sizeof(double));
74 dx = (double*)calloc(n, sizeof(double));
75 for (k = 0; k != n; ++k) { /* initial directions, based on MGJ */
76 dx[k] = fabs(x[k]) * r;
77 if (dx[k] == 0) dx[k] = r;
80 fx1 = fx = func(n, x, data); ++n_calls;
82 memcpy(x1, x, n * sizeof(double)); /* x1 = x */
83 fx1 = __kmin_hj_aux(func, n, x1, data, fx, dx, &n_calls);
85 for (k = 0; k != n; ++k) {
87 dx[k] = x1[k] > x[k]? fabs(dx[k]) : 0.0 - fabs(dx[k]);
89 x1[k] = x1[k] + x1[k] - t;
92 if (n_calls >= max_calls) break;
93 fx1 = func(n, x1, data); ++n_calls;
94 fx1 = __kmin_hj_aux(func, n, x1, data, fx1, dx, &n_calls);
96 for (k = 0; k != n; ++k)
97 if (fabs(x1[k] - x[k]) > .5 * fabs(dx[k])) break;
101 if (n_calls >= max_calls) break;
103 for (k = 0; k != n; ++k) dx[k] *= r;
104 } else break; /* converge */
110 // I copied this function somewhere several years ago with some of my modifications, but I forgot the source.
111 double kmin_brent(kmin1_f func, double a, double b, void *data, double tol, double *xmin)
113 double bound, u, r, q, fu, tmp, fa, fb, fc, c;
114 const double gold1 = 1.6180339887;
115 const double gold2 = 0.3819660113;
116 const double tiny = 1e-20;
117 const int max_iter = 100;
119 double e, d, w, v, mid, tol1, tol2, p, eold, fv, fw;
122 fa = func(a, data); fb = func(b, data);
123 if (fb > fa) { // swap, such that f(a) > f(b)
124 tmp = a; a = b; b = tmp;
125 tmp = fa; fa = fb; fb = tmp;
127 c = b + gold1 * (b - a), fc = func(c, data); // golden section extrapolation
129 bound = b + 100.0 * (c - b); // the farthest point where we want to go
130 r = (b - a) * (fb - fc);
131 q = (b - c) * (fb - fa);
132 if (fabs(q - r) < tiny) { // avoid 0 denominator
133 tmp = q > r? tiny : 0.0 - tiny;
135 u = b - ((b - c) * q - (b - a) * r) / (2.0 * tmp); // u is the parabolic extrapolation point
136 if ((b > u && u > c) || (b < u && u < c)) { // u lies between b and c
138 if (fu < fc) { // (b,u,c) bracket the minimum
139 a = b; b = u; fa = fb; fb = fu;
141 } else if (fu > fb) { // (a,b,u) bracket the minimum
145 u = c + gold1 * (c - b); fu = func(u, data); // golden section extrapolation
146 } else if ((c > u && u > bound) || (c < u && u < bound)) { // u lies between c and bound
148 if (fu < fc) { // fb > fc > fu
149 b = c; c = u; u = c + gold1 * (c - b);
150 fb = fc; fc = fu; fu = func(u, data);
151 } else { // (b,c,u) bracket the minimum
153 fa = fb; fb = fc; fc = fu;
156 } else if ((u > bound && bound > c) || (u < bound && bound < c)) { // u goes beyond the bound
157 u = bound; fu = func(u, data);
158 } else { // u goes the other way around, use golden section extrapolation
159 u = c + gold1 * (c - b); fu = func(u, data);
162 fa = fb; fb = fc; fc = fu;
164 if (a > c) u = a, a = c, c = u; // swap
166 // now, a<b<c, fa>fb and fb<fc, move on to Brent's algorithm
168 w = v = b; fv = fw = fb;
169 for (iter = 0; iter != max_iter; ++iter) {
171 tol2 = 2.0 * (tol1 = tol * fabs(b) + tiny);
172 if (fabs(b - mid) <= (tol2 - 0.5 * (c - a))) {
173 *xmin = b; return fb; // found
175 if (fabs(e) > tol1) {
176 // related to parabolic interpolation
177 r = (b - w) * (fb - fv);
178 q = (b - v) * (fb - fw);
179 p = (b - v) * q - (b - w) * r;
181 if (q > 0.0) p = 0.0 - p;
184 if (fabs(p) >= fabs(0.5 * q * eold) || p <= q * (a - b) || p >= q * (c - b)) {
185 d = gold2 * (e = (b >= mid ? a - b : c - b));
187 d = p / q; u = b + d; // actual parabolic interpolation happens here
188 if (u - a < tol2 || c - u < tol2)
189 d = (mid > b)? tol1 : 0.0 - tol1;
191 } else d = gold2 * (e = (b >= mid ? a - b : c - b)); // golden section interpolation
192 u = fabs(d) >= tol1 ? b + d : b + (d > 0.0? tol1 : -tol1);
194 if (fu <= fb) { // u is the minimum point so far
197 v = w; w = b; b = u; fv = fw; fw = fb; fb = fu;
198 } else { // adjust (a,c) and (u,v,w)
201 if (fu <= fw || w == b) {
204 } else if (fu <= fv || v == b || v == w) {