1 /* Copyright (c) 2007-2014 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 <stdlib.h>
24 #include <string.h>
25
26 #include "crs.h"
27 #include "redblack.h"
28
29 /* Controlled Random Search 2 (CRS2) with "local mutation", as defined
30 by:
31 P. Kaelo and M. M. Ali, "Some variants of the controlled random
32 search algorithm for global optimization," J. Optim. Theory Appl.
33 130 (2), 253-264 (2006).
34 */
35
36 typedef struct {
37 int n; /* # dimensions */
38 const double *lb, *ub;
39 nlopt_stopping *stop; /* stopping criteria */
40 nlopt_func f; void *f_data;
41
42 int N; /* # points in population */
43 double *ps; /* population array N x (n+1) of tuples [f(x), x] */
44 double *p; /* single point array (length n+1), for temp use */
45 rb_tree t; /* red-black tree of population, sorted by f(x) */
46 nlopt_sobol s; /* sobol data for LDS point generation, or NULL
47 to use pseudo-random numbers */
48 } crs_data;
49
50 /* sort order in red-black tree: keys [f(x), x] are sorted by f(x) */
crs_compare(double * k1,double * k2)51 static int crs_compare(double *k1, double *k2)
52 {
53 if (*k1 < *k2) return -1;
54 if (*k1 > *k2) return +1;
55 return k1 - k2; /* tie-breaker */
56 }
57
58 /* set x to a random trial value, as defined by CRS:
59 x = 2G - x_n,
60 where x_0 ... x_n are distinct points in the population
61 with x_0 the current best point and the other points are random,
62 and G is the centroid of x_0...x_{n-1} */
random_trial(crs_data * d,double * x,rb_node * best)63 static void random_trial(crs_data *d, double *x, rb_node *best)
64 {
65 int n = d->n, n1 = n+1, i, k, i0, jn;
66 double *ps = d->ps, *xi;
67
68 /* initialize x to x_0 = best point */
69 memcpy(x, best->k + 1, sizeof(double) * n);
70 i0 = (best->k - ps) / n1;
71
72 jn = nlopt_iurand(n); /* which of remaining n points is "x_n",
73 i.e. which to reflect through ...
74 this is necessary since we generate
75 the remaining points in order, so
76 just picking the last point would not
77 be very random */
78
79 /* use "method A" from
80
81 Jeffrey Scott Vitter, "An efficient algorithm for
82 sequential random sampling," ACM Trans. Math. Soft. 13 (1),
83 58--67 (1987).
84
85 to randomly pick n distinct points out of the remaining N-1 (not
86 including i0!). (The same as "method S" in Knuth vol. 2.)
87 This method requires O(N) time, which is fine in our case
88 (there are better methods if n << N). */
89 {
90 int Nleft = d->N - 1, nleft = n;
91 int Nfree = Nleft - nleft;
92 i = 0; i += i == i0;
93 while (nleft > 1) {
94 double q = ((double) Nfree) / Nleft;
95 double v = nlopt_urand(0., 1.);
96 while (q > v) {
97 ++i; i += i == i0;
98 --Nfree; --Nleft;
99 q = (q * Nfree) / Nleft;
100 }
101 xi = ps + n1 * i + 1;
102 if (jn-- == 0) /* point to reflect through */
103 for (k = 0; k < n; ++k) x[k] -= xi[k] * (0.5*n);
104 else /* point to include in centroid */
105 for (k = 0; k < n; ++k) x[k] += xi[k];
106 ++i; i += i == i0;
107 --Nleft; --nleft;
108 }
109 i += nlopt_iurand(Nleft); i += i == i0;
110 xi = ps + n1 * i + 1;
111 if (jn-- == 0) /* point to reflect through */
112 for (k = 0; k < n; ++k) x[k] -= xi[k] * (0.5*n);
113 else /* point to include in centroid */
114 for (k = 0; k < n; ++k) x[k] += xi[k];
115 }
116 for (k = 0; k < n; ++k) {
117 x[k] *= 2.0 / n; /* renormalize */
118 if (x[k] > d->ub[k]) x[k] = d->ub[k];
119 else if (x[k] < d->lb[k]) x[k] = d->lb[k];
120 }
121 }
122
123 #define NUM_MUTATION 1 /* # "local mutation" steps to try if trial fails */
124
crs_trial(crs_data * d)125 static nlopt_result crs_trial(crs_data *d)
126 {
127 rb_node *best = rb_tree_min(&d->t);
128 rb_node *worst = rb_tree_max(&d->t);
129 int mutation = NUM_MUTATION;
130 int i, n = d->n;
131 random_trial(d, d->p + 1, best);
132 do {
133 d->p[0] = d->f(n, d->p + 1, NULL, d->f_data);
134 d->stop->nevals++;
135 if (nlopt_stop_forced(d->stop)) return NLOPT_FORCED_STOP;
136 if (d->p[0] < worst->k[0]) break;
137 if (nlopt_stop_evals(d->stop)) return NLOPT_MAXEVAL_REACHED;
138 if (nlopt_stop_time(d->stop)) return NLOPT_MAXTIME_REACHED;
139 if (mutation) {
140 for (i = 0; i < n; ++i) {
141 double w = nlopt_urand(0.,1.);
142 d->p[1+i] = best->k[1+i] * (1 + w) - w * d->p[1+i];
143 if (d->p[1+i] > d->ub[i]) d->p[1+i] = d->ub[i];
144 else if (d->p[1+i] < d->lb[i]) d->p[1+i] = d->lb[i];
145 }
146 mutation--;
147 }
148 else {
149 random_trial(d, d->p + 1, best);
150 mutation = NUM_MUTATION;
151 }
152 } while (1);
153 memcpy(worst->k, d->p, sizeof(double) * (n+1));
154 rb_tree_resort(&d->t, worst);
155 return NLOPT_SUCCESS;
156 }
157
crs_destroy(crs_data * d)158 static void crs_destroy(crs_data *d)
159 {
160 nlopt_sobol_destroy(d->s);
161 rb_tree_destroy(&d->t);
162 free(d->ps);
163 }
164
crs_init(crs_data * d,int n,const double * x,const double * lb,const double * ub,nlopt_stopping * stop,nlopt_func f,void * f_data,int population,int lds)165 static nlopt_result crs_init(crs_data *d, int n, const double *x,
166 const double *lb, const double *ub,
167 nlopt_stopping *stop, nlopt_func f, void *f_data,
168 int population, int lds)
169 {
170 int i;
171
172 if (!population) {
173 /* TODO: how should we set the default population size?
174 the Kaelo and Ali paper suggests 10*(n+1), but should
175 we add more random points if maxeval is large, or... ? */
176 d->N = 10 * (n + 1); /* heuristic initial population size */
177 }
178 else
179 d->N = population;
180 if (d->N < n + 1) /* population must be big enough for a simplex */
181 return NLOPT_INVALID_ARGS;
182
183 d->n = n;
184 d->stop = stop;
185 d->f = f; d->f_data = f_data;
186 d->ub = ub; d->lb = lb;
187 d->ps = (double *) malloc(sizeof(double) * (n + 1) * (d->N + 1));
188 if (!d->ps) return NLOPT_OUT_OF_MEMORY;
189 d->p = d->ps + d->N * (n+1);
190 rb_tree_init(&d->t, crs_compare);
191
192 /* we can either use pseudorandom points, as in the original CRS
193 algorithm, or use a low-discrepancy Sobol' sequence ... I tried
194 the latter, however, and it doesn't seem to help, probably
195 because we are only generating a small number of random points
196 to start with */
197 d->s = lds ? nlopt_sobol_create((unsigned) n) : NULL;
198 nlopt_sobol_skip(d->s, (unsigned) d->N, d->ps + 1);
199
200 /* generate initial points randomly, plus starting guess x */
201 memcpy(d->ps + 1, x, sizeof(double) * n);
202 d->ps[0] = f(n, x, NULL, f_data);
203 stop->nevals++;
204 if (!rb_tree_insert(&d->t, d->ps)) return NLOPT_OUT_OF_MEMORY;
205 if (d->ps[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
206 if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
207 if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
208 for (i = 1; i < d->N; ++i) {
209 double *k = d->ps + i*(n+1);
210 if (d->s)
211 nlopt_sobol_next(d->s, k + 1, lb, ub);
212 else {
213 int j;
214 for (j = 0; j < n; ++j)
215 k[1 + j] = nlopt_urand(lb[j], ub[j]);
216 }
217 k[0] = f(n, k + 1, NULL, f_data);
218 stop->nevals++;
219 if (!rb_tree_insert(&d->t, k)) return NLOPT_OUT_OF_MEMORY;
220 if (k[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
221 if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
222 if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
223 }
224
225 return NLOPT_SUCCESS;;
226 }
227
crs_minimize(int n,nlopt_func f,void * f_data,const double * lb,const double * ub,double * x,double * minf,nlopt_stopping * stop,int population,int lds)228 nlopt_result crs_minimize(int n, nlopt_func f, void *f_data,
229 const double *lb, const double *ub, /* bounds */
230 double *x, /* in: initial guess, out: minimizer */
231 double *minf,
232 nlopt_stopping *stop,
233 int population, /* initial population (0=default) */
234 int lds) /* random or low-discrepancy seq. (lds) */
235 {
236 nlopt_result ret;
237 crs_data d;
238 rb_node *best;
239
240 ret = crs_init(&d, n, x, lb, ub, stop, f, f_data, population, lds);
241 if (ret < 0) return ret;
242
243 best = rb_tree_min(&d.t);
244 *minf = best->k[0];
245 memcpy(x, best->k + 1, sizeof(double) * n);
246
247 while (ret == NLOPT_SUCCESS) {
248 if (NLOPT_SUCCESS == (ret = crs_trial(&d))) {
249 best = rb_tree_min(&d.t);
250 if (best->k[0] < *minf) {
251 if (best->k[0] < stop->minf_max)
252 ret = NLOPT_MINF_MAX_REACHED;
253 else if (nlopt_stop_f(stop, best->k[0], *minf))
254 ret = NLOPT_FTOL_REACHED;
255 else if (nlopt_stop_x(stop, best->k + 1, x))
256 ret = NLOPT_XTOL_REACHED;
257 *minf = best->k[0];
258 memcpy(x, best->k + 1, sizeof(double) * n);
259 }
260 if (ret != NLOPT_SUCCESS) {
261 if (nlopt_stop_evals(stop))
262 ret = NLOPT_MAXEVAL_REACHED;
263 else if (nlopt_stop_time(stop))
264 ret = NLOPT_MAXTIME_REACHED;
265 }
266 }
267 }
268 crs_destroy(&d);
269 return ret;
270 }
271