1 /* Copyright (c) 2007-2012 Massachusetts Institute of Technology
3 * Permission is hereby granted, free of charge, to any person obtaining
4 * a copy of this software and associated documentation files (the
5 * "Software"), to deal in the Software without restriction, including
6 * without limitation the rights to use, copy, modify, merge, publish,
7 * distribute, sublicense, and/or sell copies of the Software, and to
8 * permit persons to whom the Software is furnished to do so, subject to
9 * the following conditions:
11 * The above copyright notice and this permission notice shall be
12 * included in all copies or substantial portions of the Software.
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
18 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
19 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
20 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
27 #include "neldermead.h"
30 /* Nelder-Mead simplex algorithm, used as a subroutine for the Rowan's
31 subplex algorithm. Modified to handle bound constraints ala
32 Richardson and Kuester (1973), as mentioned below. */
34 /* heuristic "strategy" constants: */
35 static const double alpha = 1, beta = 0.5, gamm = 2, delta = 0.5;
37 /* sort order in red-black tree: keys [f(x), x] are sorted by f(x) */
38 static int simplex_compare(double *k1, double *k2)
40 if (*k1 < *k2) return -1;
41 if (*k1 > *k2) return +1;
42 return k1 - k2; /* tie-breaker */
45 /* return 1 if a and b are approximately equal relative to floating-point
46 precision, 0 otherwise */
47 static int close(double a, double b)
49 return (fabs(a - b) <= 1e-13 * (fabs(a) + fabs(b)));
52 /* Perform the reflection xnew = c + scale * (c - xold),
53 returning 0 if xnew == c or xnew == xold (coincident points), 1 otherwise.
55 The reflected point xnew is "pinned" to the lower and upper bounds
56 (lb and ub), as suggested by J. A. Richardson and J. L. Kuester,
57 "The complex method for constrained optimization," Commun. ACM
58 16(8), 487-489 (1973). This is probably a suboptimal way to handle
59 bound constraints, but I don't know a better way. The main danger
60 with this is that the simplex might collapse into a
61 lower-dimensional hyperplane; this danger can be ameliorated by
62 restarting (as in subplex), however. */
63 static int reflectpt(int n, double *xnew,
64 const double *c, double scale, const double *xold,
65 const double *lb, const double *ub)
67 int equalc = 1, equalold = 1, i;
68 for (i = 0; i < n; ++i) {
69 double newx = c[i] + scale * (c[i] - xold[i]);
70 if (newx < lb[i]) newx = lb[i];
71 if (newx > ub[i]) newx = ub[i];
72 equalc = equalc && close(newx, c[i]);
73 equalold = equalold && close(newx, xold[i]);
76 return !(equalc || equalold);
79 #define CHECK_EVAL(xc,fc) \
81 if (nlopt_stop_forced(stop)) { ret=NLOPT_FORCED_STOP; goto done; } \
82 if ((fc) <= *minf) { \
83 *minf = (fc); memcpy(x, (xc), n * sizeof(double)); \
84 if (*minf < stop->minf_max) { ret=NLOPT_MINF_MAX_REACHED; goto done; } \
86 if (nlopt_stop_evals(stop)) { ret=NLOPT_MAXEVAL_REACHED; goto done; } \
87 if (nlopt_stop_time(stop)) { ret=NLOPT_MAXTIME_REACHED; goto done; }
89 /* Internal version of nldrmd_minimize, intended to be used as
90 a subroutine for the subplex method. Three differences compared
93 *minf should contain the value of f(x) (so that we don't have to
94 re-evaluate f at the starting x).
96 if psi > 0, then it *replaces* xtol and ftol in stop with the condition
97 that the simplex diameter |xl - xh| must be reduced by a factor of psi
98 ... this is for when nldrmd is used within the subplex method; for
99 ordinary termination tests, set psi = 0.
101 scratch should contain an array of length >= (n+1)*(n+1) + 2*n,
102 used as scratch workspace.
104 On output, *fdiff will contain the difference between the high
105 and low function values of the last simplex. */
106 nlopt_result nldrmd_minimize_(int n, nlopt_func f, void *f_data,
107 const double *lb, const double *ub, /* bounds */
108 double *x, /* in: initial guess, out: minimizer */
110 const double *xstep, /* initial step sizes */
111 nlopt_stopping *stop,
112 double psi, double *scratch,
115 double *pts; /* (n+1) x (n+1) array of n+1 points plus function val [0] */
116 double *c; /* centroid * n */
117 double *xcur; /* current point */
118 rb_tree t; /* red-black tree of simplex, sorted by f(x) */
120 double ninv = 1.0 / n;
121 nlopt_result ret = NLOPT_SUCCESS;
122 double init_diam = 0;
125 c = scratch + (n+1)*(n+1);
128 rb_tree_init(&t, simplex_compare);
132 /* initialize the simplex based on the starting xstep */
133 memcpy(pts+1, x, sizeof(double)*n);
135 if (*minf < stop->minf_max) { ret=NLOPT_MINF_MAX_REACHED; goto done; }
136 for (i = 0; i < n; ++i) {
137 double *pt = pts + (i+1)*(n+1);
138 memcpy(pt+1, x, sizeof(double)*n);
140 if (pt[1+i] > ub[i]) {
141 if (ub[i] - x[i] > fabs(xstep[i]) * 0.1)
143 else /* ub is too close to pt, go in other direction */
144 pt[1+i] = x[i] - fabs(xstep[i]);
146 if (pt[1+i] < lb[i]) {
147 if (x[i] - lb[i] > fabs(xstep[i]) * 0.1)
149 else {/* lb is too close to pt, go in other direction */
150 pt[1+i] = x[i] + fabs(xstep[i]);
151 if (pt[1+i] > ub[i]) /* go towards further of lb, ub */
152 pt[1+i] = 0.5 * ((ub[i] - x[i] > x[i] - lb[i] ?
153 ub[i] : lb[i]) + x[i]);
156 if (close(pt[1+i], x[i])) { ret=NLOPT_FAILURE; goto done; }
157 pt[0] = f(n, pt+1, NULL, f_data);
158 CHECK_EVAL(pt+1, pt[0]);
162 for (i = 0; i < n + 1; ++i)
163 if (!rb_tree_insert(&t, pts + i*(n+1))) {
164 ret = NLOPT_OUT_OF_MEMORY;
169 rb_node *low = rb_tree_min(&t);
170 rb_node *high = rb_tree_max(&t);
171 double fl = low->k[0], *xl = low->k + 1;
172 double fh = high->k[0], *xh = high->k + 1;
177 if (init_diam == 0) /* initialize diam. for psi convergence test */
178 for (i = 0; i < n; ++i) init_diam += fabs(xl[i] - xh[i]);
180 if (psi <= 0 && nlopt_stop_ftol(stop, fl, fh)) {
181 ret = NLOPT_FTOL_REACHED;
185 /* compute centroid ... if we cared about the perfomance of this,
186 we could do it iteratively by updating the centroid on
187 each step, but then we would have to be more careful about
188 accumulation of rounding errors... anyway n is unlikely to
189 be very large for Nelder-Mead in practical cases */
190 memset(c, 0, sizeof(double)*n);
191 for (i = 0; i < n + 1; ++i) {
192 double *xi = pts + i*(n+1) + 1;
194 for (j = 0; j < n; ++j)
197 for (i = 0; i < n; ++i) c[i] *= ninv;
199 /* x convergence check: find xcur = max radius from centroid */
200 memset(xcur, 0, sizeof(double)*n);
201 for (i = 0; i < n + 1; ++i) {
202 double *xi = pts + i*(n+1) + 1;
203 for (j = 0; j < n; ++j) {
204 double dx = fabs(xi[j] - c[j]);
205 if (dx > xcur[j]) xcur[j] = dx;
208 for (i = 0; i < n; ++i) xcur[i] += c[i];
211 for (i = 0; i < n; ++i) diam += fabs(xl[i] - xh[i]);
212 if (diam < psi * init_diam) {
213 ret = NLOPT_XTOL_REACHED;
217 else if (nlopt_stop_x(stop, c, xcur)) {
218 ret = NLOPT_XTOL_REACHED;
223 if (!reflectpt(n, xcur, c, alpha, xh, lb, ub)) {
224 ret=NLOPT_XTOL_REACHED; goto done;
226 fr = f(n, xcur, NULL, f_data);
227 CHECK_EVAL(xcur, fr);
229 if (fr < fl) { /* new best point, expand simplex */
230 if (!reflectpt(n, xh, c, gamm, xh, lb, ub)) {
231 ret=NLOPT_XTOL_REACHED; goto done;
233 fh = f(n, xh, NULL, f_data);
235 if (fh >= fr) { /* expanding didn't improve */
237 memcpy(xh, xcur, sizeof(double)*n);
240 else if (fr < rb_tree_pred(high)->k[0]) { /* accept new point */
241 memcpy(xh, xcur, sizeof(double)*n);
244 else { /* new worst point, contract */
246 if (!reflectpt(n,xcur,c, fh <= fr ? -beta : beta, xh, lb,ub)) {
247 ret=NLOPT_XTOL_REACHED; goto done;
249 fc = f(n, xcur, NULL, f_data);
250 CHECK_EVAL(xcur, fc);
251 if (fc < fr && fc < fh) { /* successful contraction */
252 memcpy(xh, xcur, sizeof(double)*n);
255 else { /* failed contraction, shrink simplex */
257 rb_tree_init(&t, simplex_compare);
258 for (i = 0; i < n+1; ++i) {
259 double *pt = pts + i * (n+1);
261 if (!reflectpt(n,pt+1, xl,-delta,pt+1, lb,ub)) {
262 ret = NLOPT_XTOL_REACHED;
265 pt[0] = f(n, pt+1, NULL, f_data);
266 CHECK_EVAL(pt+1, pt[0]);
274 rb_tree_resort(&t, high);
282 nlopt_result nldrmd_minimize(int n, nlopt_func f, void *f_data,
283 const double *lb, const double *ub, /* bounds */
284 double *x, /* in: initial guess, out: minimizer */
286 const double *xstep, /* initial step sizes */
287 nlopt_stopping *stop)
290 double *scratch, fdiff;
292 *minf = f(n, x, NULL, f_data);
294 if (nlopt_stop_forced(stop)) return NLOPT_FORCED_STOP;
295 if (*minf < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
296 if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
297 if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
299 scratch = (double*) malloc(sizeof(double) * ((n+1)*(n+1) + 2*n));
300 if (!scratch) return NLOPT_OUT_OF_MEMORY;
302 ret = nldrmd_minimize_(n, f, f_data, lb, ub, x, minf, xstep, stop,
303 0.0, scratch, &fdiff);