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.
29 /* Controlled Random Search 2 (CRS2) with "local mutation", as defined
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).
37 int n; /* # dimensions */
38 const double *lb, *ub;
39 nlopt_stopping *stop; /* stopping criteria */
40 nlopt_func f; void *f_data;
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 */
50 /* sort order in red-black tree: keys [f(x), x] are sorted by f(x) */
51 static int crs_compare(double *k1, double *k2)
53 if (*k1 < *k2) return -1;
54 if (*k1 > *k2) return +1;
55 return k1 - k2; /* tie-breaker */
58 /* set x to a random trial value, as defined by CRS:
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} */
63 static void random_trial(crs_data *d, double *x, rb_node *best)
65 int n = d->n, n1 = n+1, i, k, i0, jn;
66 double *ps = d->ps, *xi;
68 /* initialize x to x_0 = best point */
69 memcpy(x, best->k + 1, sizeof(double) * n);
70 i0 = (best->k - ps) / n1;
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
79 /* use "method A" from
81 Jeffrey Scott Vitter, "An efficient algorithm for
82 sequential random sampling," ACM Trans. Math. Soft. 13 (1),
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). */
90 int Nleft = d->N - 1, nleft = n;
91 int Nfree = Nleft - nleft;
94 double q = ((double) Nfree) / Nleft;
95 double v = nlopt_urand(0., 1.);
99 q = (q * Nfree) / Nleft;
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];
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];
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];
123 #define NUM_MUTATION 1 /* # "local mutation" steps to try if trial fails */
125 static nlopt_result crs_trial(crs_data *d)
127 rb_node *best = rb_tree_min(&d->t);
128 rb_node *worst = rb_tree_max(&d->t);
129 int mutation = NUM_MUTATION;
131 random_trial(d, d->p + 1, best);
133 d->p[0] = d->f(n, d->p + 1, NULL, d->f_data);
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;
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];
149 random_trial(d, d->p + 1, best);
150 mutation = NUM_MUTATION;
153 memcpy(worst->k, d->p, sizeof(double) * (n+1));
154 rb_tree_resort(&d->t, worst);
155 return NLOPT_SUCCESS;
158 static void crs_destroy(crs_data *d)
160 nlopt_sobol_destroy(d->s);
161 rb_tree_destroy(&d->t);
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)
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 */
180 if (d->N < n + 1) { /* population must be big enough for a simplex */
181 nlopt_stop_msg(stop, "population %d should be >= dimension + 1 = %d",
183 return NLOPT_INVALID_ARGS;
188 d->f = f; d->f_data = f_data;
189 d->ub = ub; d->lb = lb;
190 d->ps = (double *) malloc(sizeof(double) * (n + 1) * (d->N + 1));
191 if (!d->ps) return NLOPT_OUT_OF_MEMORY;
192 d->p = d->ps + d->N * (n+1);
193 rb_tree_init(&d->t, crs_compare);
195 /* we can either use pseudorandom points, as in the original CRS
196 algorithm, or use a low-discrepancy Sobol' sequence ... I tried
197 the latter, however, and it doesn't seem to help, probably
198 because we are only generating a small number of random points
200 d->s = lds ? nlopt_sobol_create((unsigned) n) : NULL;
201 nlopt_sobol_skip(d->s, (unsigned) d->N, d->ps + 1);
203 /* generate initial points randomly, plus starting guess x */
204 memcpy(d->ps + 1, x, sizeof(double) * n);
205 d->ps[0] = f(n, x, NULL, f_data);
207 if (!rb_tree_insert(&d->t, d->ps)) return NLOPT_OUT_OF_MEMORY;
208 if (d->ps[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
209 if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
210 if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
211 for (i = 1; i < d->N; ++i) {
212 double *k = d->ps + i*(n+1);
214 nlopt_sobol_next(d->s, k + 1, lb, ub);
217 for (j = 0; j < n; ++j)
218 k[1 + j] = nlopt_urand(lb[j], ub[j]);
220 k[0] = f(n, k + 1, NULL, f_data);
222 if (!rb_tree_insert(&d->t, k)) return NLOPT_OUT_OF_MEMORY;
223 if (k[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
224 if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
225 if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
228 return NLOPT_SUCCESS;;
231 nlopt_result crs_minimize(int n, nlopt_func f, void *f_data,
232 const double *lb, const double *ub, /* bounds */
233 double *x, /* in: initial guess, out: minimizer */
235 nlopt_stopping *stop,
236 int population, /* initial population (0=default) */
237 int lds) /* random or low-discrepancy seq. (lds) */
243 ret = crs_init(&d, n, x, lb, ub, stop, f, f_data, population, lds);
244 if (ret < 0) return ret;
246 best = rb_tree_min(&d.t);
248 memcpy(x, best->k + 1, sizeof(double) * n);
250 while (ret == NLOPT_SUCCESS) {
251 if (NLOPT_SUCCESS == (ret = crs_trial(&d))) {
252 best = rb_tree_min(&d.t);
253 if (best->k[0] < *minf) {
254 if (best->k[0] < stop->minf_max)
255 ret = NLOPT_MINF_MAX_REACHED;
256 else if (nlopt_stop_f(stop, best->k[0], *minf))
257 ret = NLOPT_FTOL_REACHED;
258 else if (nlopt_stop_x(stop, best->k + 1, x))
259 ret = NLOPT_XTOL_REACHED;
261 memcpy(x, best->k + 1, sizeof(double) * n);
263 if (ret != NLOPT_SUCCESS) {
264 if (nlopt_stop_evals(stop))
265 ret = NLOPT_MAXEVAL_REACHED;
266 else if (nlopt_stop_time(stop))
267 ret = NLOPT_MAXTIME_REACHED;