1 /* Copyright (c) 2007-2010 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_FORCE_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 return NLOPT_INVALID_ARGS;
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);
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
197 d->s = lds ? nlopt_sobol_create((unsigned) n) : NULL;
198 nlopt_sobol_skip(d->s, (unsigned) d->N, d->ps + 1);
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);
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);
211 nlopt_sobol_next(d->s, k + 1, lb, ub);
214 for (j = 0; j < n; ++j)
215 k[1 + j] = nlopt_urand(lb[j], ub[j]);
217 k[0] = f(n, k + 1, NULL, f_data);
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;
225 return NLOPT_SUCCESS;;
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 */
232 nlopt_stopping *stop,
233 int population, /* initial population (0=default) */
234 int lds) /* random or low-discrepancy seq. (lds) */
240 ret = crs_init(&d, n, x, lb, ub, stop, f, f_data, population, lds);
241 if (ret < 0) return ret;
243 best = rb_tree_min(&d.t);
245 memcpy(x, best->k + 1, sizeof(double) * n);
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;
258 memcpy(x, best->k + 1, sizeof(double) * n);
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;