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