1 /* Copyright (c) 2007-2014 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])) {
157 nlopt_stop_msg(stop, "starting step size led to simplex that was too small in dimension %d: %g is too close to x[%d]=%g",
158 i, pt[1+i], i, x[i]);
162 pt[0] = f(n, pt+1, NULL, f_data);
163 CHECK_EVAL(pt+1, pt[0]);
167 for (i = 0; i < n + 1; ++i)
168 if (!rb_tree_insert(&t, pts + i*(n+1))) {
169 ret = NLOPT_OUT_OF_MEMORY;
174 rb_node *low = rb_tree_min(&t);
175 rb_node *high = rb_tree_max(&t);
176 double fl = low->k[0], *xl = low->k + 1;
177 double fh = high->k[0], *xh = high->k + 1;
182 if (init_diam == 0) /* initialize diam. for psi convergence test */
183 for (i = 0; i < n; ++i) init_diam += fabs(xl[i] - xh[i]);
185 if (psi <= 0 && nlopt_stop_ftol(stop, fl, fh)) {
186 ret = NLOPT_FTOL_REACHED;
190 /* compute centroid ... if we cared about the perfomance of this,
191 we could do it iteratively by updating the centroid on
192 each step, but then we would have to be more careful about
193 accumulation of rounding errors... anyway n is unlikely to
194 be very large for Nelder-Mead in practical cases */
195 memset(c, 0, sizeof(double)*n);
196 for (i = 0; i < n + 1; ++i) {
197 double *xi = pts + i*(n+1) + 1;
199 for (j = 0; j < n; ++j)
202 for (i = 0; i < n; ++i) c[i] *= ninv;
204 /* x convergence check: find xcur = max radius from centroid */
205 memset(xcur, 0, sizeof(double)*n);
206 for (i = 0; i < n + 1; ++i) {
207 double *xi = pts + i*(n+1) + 1;
208 for (j = 0; j < n; ++j) {
209 double dx = fabs(xi[j] - c[j]);
210 if (dx > xcur[j]) xcur[j] = dx;
213 for (i = 0; i < n; ++i) xcur[i] += c[i];
216 for (i = 0; i < n; ++i) diam += fabs(xl[i] - xh[i]);
217 if (diam < psi * init_diam) {
218 ret = NLOPT_XTOL_REACHED;
222 else if (nlopt_stop_x(stop, c, xcur)) {
223 ret = NLOPT_XTOL_REACHED;
228 if (!reflectpt(n, xcur, c, alpha, xh, lb, ub)) {
229 ret=NLOPT_XTOL_REACHED; goto done;
231 fr = f(n, xcur, NULL, f_data);
232 CHECK_EVAL(xcur, fr);
234 if (fr < fl) { /* new best point, expand simplex */
235 if (!reflectpt(n, xh, c, gamm, xh, lb, ub)) {
236 ret=NLOPT_XTOL_REACHED; goto done;
238 fh = f(n, xh, NULL, f_data);
240 if (fh >= fr) { /* expanding didn't improve */
242 memcpy(xh, xcur, sizeof(double)*n);
245 else if (fr < rb_tree_pred(high)->k[0]) { /* accept new point */
246 memcpy(xh, xcur, sizeof(double)*n);
249 else { /* new worst point, contract */
251 if (!reflectpt(n,xcur,c, fh <= fr ? -beta : beta, xh, lb,ub)) {
252 ret=NLOPT_XTOL_REACHED; goto done;
254 fc = f(n, xcur, NULL, f_data);
255 CHECK_EVAL(xcur, fc);
256 if (fc < fr && fc < fh) { /* successful contraction */
257 memcpy(xh, xcur, sizeof(double)*n);
260 else { /* failed contraction, shrink simplex */
262 rb_tree_init(&t, simplex_compare);
263 for (i = 0; i < n+1; ++i) {
264 double *pt = pts + i * (n+1);
266 if (!reflectpt(n,pt+1, xl,-delta,pt+1, lb,ub)) {
267 ret = NLOPT_XTOL_REACHED;
270 pt[0] = f(n, pt+1, NULL, f_data);
271 CHECK_EVAL(pt+1, pt[0]);
279 rb_tree_resort(&t, high);
287 nlopt_result nldrmd_minimize(int n, nlopt_func f, void *f_data,
288 const double *lb, const double *ub, /* bounds */
289 double *x, /* in: initial guess, out: minimizer */
291 const double *xstep, /* initial step sizes */
292 nlopt_stopping *stop)
295 double *scratch, fdiff;
297 *minf = f(n, x, NULL, f_data);
299 if (nlopt_stop_forced(stop)) return NLOPT_FORCED_STOP;
300 if (*minf < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
301 if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
302 if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
304 scratch = (double*) malloc(sizeof(double) * ((n+1)*(n+1) + 2*n));
305 if (!scratch) return NLOPT_OUT_OF_MEMORY;
307 ret = nldrmd_minimize_(n, f, f_data, lb, ub, x, minf, xstep, stop,
308 0.0, scratch, &fdiff);