chiark / gitweb /
added nlopt_force_stop termination
[nlopt.git] / neldermead / sbplx.c
1 /* Copyright (c) 2007-2010 Massachusetts Institute of Technology
2  *
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:
10  * 
11  * The above copyright notice and this permission notice shall be
12  * included in all copies or substantial portions of the Software.
13  * 
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. 
21  */
22
23 #include <math.h>
24 #include <stdlib.h>
25 #include <stdio.h>
26 #include <string.h>
27
28 #include "neldermead.h"
29
30 /* subplex strategy constants: */
31 static const double psi = 0.25, omega = 0.1;
32 static const int nsmin = 2, nsmax = 5;
33
34 int sbplx_verbose = 0; /* for debugging */
35
36 /* qsort_r comparison function for sorting indices into delta-x array */
37 static int p_compare(void *dx_, const void *i_, const void *j_)
38 {
39      const double *dx = (const double *) dx_;
40      int i = *((const int *) i_), j = *((const int *) j_);
41      double dxi = fabs(dx[i]), dxj = fabs(dx[j]);
42      return (dxi > dxj ? -1 : (dxi < dxj ? +1 : 0));
43 }
44
45 typedef struct {
46      const int *p; /* subspace index permutation */
47      int is; /* starting index for this subspace */
48      int n; /* dimension of underlying space */
49      double *x; /* current x vector */
50      nlopt_func f; void *f_data; /* the "actual" underlying function */
51 } subspace_data;
52
53 /* wrapper around objective function for subspace optimization */
54 static double subspace_func(int ns, const double *xs, double *grad, void *data)
55 {
56      subspace_data *d = (subspace_data *) data;
57      int i, is = d->is;
58      const int *p = d->p;
59      double *x = d->x;
60
61      (void) grad; /* should always be NULL here */
62      for (i = is; i < is + ns; ++i) x[p[i]] = xs[i-is];
63      return d->f(d->n, x, NULL, d->f_data);
64 }
65
66 nlopt_result sbplx_minimize(int n, nlopt_func f, void *f_data,
67                             const double *lb, const double *ub, /* bounds */
68                             double *x, /* in: initial guess, out: minimizer */
69                             double *minf,
70                             const double *xstep0, /* initial step sizes */
71                             nlopt_stopping *stop)
72 {
73      nlopt_result ret = NLOPT_SUCCESS;
74      double *xstep, *xprev, *dx, *xs, *lbs, *ubs, *xsstep, *scratch;
75      int *p; /* permuted indices of x sorted by decreasing magnitude |dx| */
76      int i;
77      subspace_data sd;
78      double fprev;
79
80      *minf = f(n, x, NULL, f_data);
81      stop->nevals++;
82      if (nlopt_stop_forced(stop)) return NLOPT_FORCE_STOP;
83      if (*minf < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
84      if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
85      if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
86
87      xstep = (double*)malloc(sizeof(double) * (n*3 + nsmax*4
88                                                + (nsmax+1)*(nsmax+1)+2*nsmax));
89      if (!xstep) return NLOPT_OUT_OF_MEMORY;
90      xprev = xstep + n; dx = xprev + n;
91      xs = dx + n; xsstep = xs + nsmax; 
92      lbs = xsstep + nsmax; ubs = lbs + nsmax;
93      scratch = ubs + nsmax;
94      p = (int *) malloc(sizeof(int) * n);
95      if (!p) { free(xstep); return NLOPT_OUT_OF_MEMORY; }
96
97      memcpy(xstep, xstep0, n * sizeof(double));
98      memset(dx, 0, n * sizeof(double));
99
100      sd.p = p;
101      sd.n = n;
102      sd.x = x;
103      sd.f = f;
104      sd.f_data = f_data;
105
106      while (1) {
107           double normi = 0;
108           double normdx = 0;
109           int ns, nsubs = 0;
110           int nevals = stop->nevals;
111           double fdiff, fdiff_max = 0;
112
113           memcpy(xprev, x, n * sizeof(double));
114           fprev = *minf;
115
116           /* sort indices into the progress vector dx by decreasing
117              order of magnitude |dx| */
118           for (i = 0; i < n; ++i) p[i] = i;
119           nlopt_qsort_r(p, (size_t) n, sizeof(int), dx, p_compare);
120
121           /* find the subspaces, and perform nelder-mead on each one */
122           for (i = 0; i < n; ++i) normdx += fabs(dx[i]); /* L1 norm */
123           for (i = 0; i + nsmin < n; i += ns) {
124                /* find subspace starting at index i */
125                int k, nk;
126                double ns_goodness = -HUGE_VAL, norm = normi;
127                nk = i+nsmax > n ? n : i+nsmax; /* max k for this subspace */
128                for (k = i; k < i+nsmin-1; ++k) norm += fabs(dx[p[k]]);
129                ns = nsmin;
130                for (k = i+nsmin-1; k < nk; ++k) {
131                     double goodness;
132                     norm += fabs(dx[p[k]]);
133                     /* remaining subspaces must be big enough to partition */
134                     if (n-(k+1) < nsmin) continue;
135                     /* maximize figure of merit defined by Rowan thesis:
136                        look for sudden drops in average |dx| */
137                     if (k+1 < n)
138                          goodness = norm/(k+1) - (normdx-norm)/(n-(k+1));
139                     else
140                          goodness = normdx/n;
141                     if (goodness > ns_goodness) {
142                          ns_goodness = goodness;
143                          ns = (k+1)-i;
144                     }
145                }
146                for (k = i; k < i+ns; ++k) normi += fabs(dx[p[k]]);
147                /* do nelder-mead on subspace of dimension ns starting w/i */
148                sd.is = i;
149                for (k = i; k < i+ns; ++k) {
150                     xs[k-i] = x[p[k]];
151                     xsstep[k-i] = xstep[p[k]];
152                     lbs[k-i] = lb[p[k]];
153                     ubs[k-i] = ub[p[k]];
154                }
155                ++nsubs;
156                nevals = stop->nevals;
157                ret = nldrmd_minimize_(ns, subspace_func, &sd, lbs,ubs,xs, minf,
158                                       xsstep, stop, psi, scratch, &fdiff);
159                if (fdiff > fdiff_max) fdiff_max = fdiff;
160                if (sbplx_verbose)
161                     printf("%d NM iterations for (%d,%d) subspace\n",
162                            stop->nevals - nevals, sd.is, ns);
163                for (k = i; k < i+ns; ++k) x[p[k]] = xs[k-i];
164                if (ret == NLOPT_FAILURE) { ret=NLOPT_XTOL_REACHED; goto done; }
165                if (ret != NLOPT_XTOL_REACHED) goto done;
166           }
167           /* nelder-mead on last subspace */
168           ns = n - i;
169           sd.is = i;
170           for (; i < n; ++i) {
171                xs[i-sd.is] = x[p[i]];
172                xsstep[i-sd.is] = xstep[p[i]];
173                lbs[i-sd.is] = lb[p[i]];
174                ubs[i-sd.is] = ub[p[i]];
175           }
176           ++nsubs;
177           nevals = stop->nevals;
178           ret = nldrmd_minimize_(ns, subspace_func, &sd, lbs,ubs,xs, minf,
179                                  xsstep, stop, psi, scratch, &fdiff);
180           if (fdiff > fdiff_max) fdiff_max = fdiff;
181           if (sbplx_verbose)
182                printf("sbplx: %d NM iterations for (%d,%d) subspace\n",
183                       stop->nevals - nevals, sd.is, ns);
184           for (i = sd.is; i < n; ++i) x[p[i]] = xs[i-sd.is];
185           if (ret == NLOPT_FAILURE) { ret=NLOPT_XTOL_REACHED; goto done; }
186           if (ret != NLOPT_XTOL_REACHED) goto done;
187
188           /* termination tests: */
189           if (nlopt_stop_ftol(stop, *minf, *minf + fdiff_max)) {
190                ret = NLOPT_FTOL_REACHED;
191                goto done;
192           }
193           if (nlopt_stop_x(stop, x, xprev)) {
194                int j;
195                /* as explained in Rowan's thesis, it is important
196                   to check |xstep| as well as |x-xprev|, since if
197                   the step size is too large (in early iterations),
198                   the inner Nelder-Mead may not make much progress */
199                for (j = 0; j < n; ++j)
200                     if (fabs(xstep[j]) * psi > stop->xtol_abs[j]
201                         && fabs(xstep[j]) * psi > stop->xtol_rel * fabs(x[j]))
202                          break;
203                if (j == n) {
204                     ret = NLOPT_XTOL_REACHED;
205                     goto done;
206                }
207           }
208
209           /* compute change in optimal point */
210           for (i = 0; i < n; ++i) dx[i] = x[i] - xprev[i];
211
212           /* setting stepsizes */
213           {
214                double scale;
215                if (nsubs == 1)
216                     scale = psi;
217                else {
218                     double stepnorm = 0, dxnorm = 0;
219                     for (i = 0; i < n; ++i) {
220                          stepnorm += fabs(xstep[i]);
221                          dxnorm += fabs(dx[i]);
222                     }
223                     scale = dxnorm / stepnorm;
224                     if (scale < omega) scale = omega;
225                     if (scale > 1/omega) scale = 1/omega;
226                }
227                if (sbplx_verbose)
228                     printf("sbplx: stepsize scale factor = %g\n", scale);
229                for (i = 0; i < n; ++i) 
230                     xstep[i] = (dx[i] == 0) ? -(xstep[i] * scale)
231                          : copysign(xstep[i] * scale, dx[i]);
232           }
233      }
234
235  done:
236      free(p);
237      free(xstep);
238      return ret;
239 }