1 /* Copyright (c) 2007-2008 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 (d->p[0] < worst->k[0]) break;
136 if (nlopt_stop_evals(d->stop)) return NLOPT_MAXEVAL_REACHED;
137 if (nlopt_stop_time(d->stop)) return NLOPT_MAXTIME_REACHED;
139 for (i = 0; i < n; ++i) {
140 double w = nlopt_urand(0.,1.);
141 d->p[1+i] = best->k[1+i] * (1 + w) - w * d->p[1+i];
142 if (d->p[1+i] > d->ub[i]) d->p[1+i] = d->ub[i];
143 else if (d->p[1+i] < d->lb[i]) d->p[1+i] = d->lb[i];
148 random_trial(d, d->p + 1, best);
149 mutation = NUM_MUTATION;
152 memcpy(worst->k, d->p, sizeof(double) * (n+1));
153 rb_tree_resort(&d->t, worst);
154 return NLOPT_SUCCESS;
157 static void crs_destroy(crs_data *d)
159 nlopt_sobol_destroy(d->s);
160 rb_tree_destroy(&d->t);
164 static nlopt_result crs_init(crs_data *d, int n, const double *x,
165 const double *lb, const double *ub,
166 nlopt_stopping *stop, nlopt_func f, void *f_data,
171 /* TODO: how should we set the initial population size?
172 the Kaelo and Ali paper suggests 10*(n+1), but should
173 we add more random points if maxeval is large, or... ? */
174 d->N = 10 * (n + 1); /* heuristic initial population size */
178 d->f = f; d->f_data = f_data;
179 d->ub = ub; d->lb = lb;
180 d->ps = (double *) malloc(sizeof(double) * (n + 1) * (d->N + 1));
181 if (!d->ps) return NLOPT_OUT_OF_MEMORY;
182 d->p = d->ps + d->N * (n+1);
183 rb_tree_init(&d->t, crs_compare);
185 /* we can either use pseudorandom points, as in the original CRS
186 algorithm, or use a low-discrepancy Sobol' sequence ... I tried
187 the latter, however, and it doesn't seem to help, probably
188 because we are only generating a small number of random points
190 d->s = lds ? nlopt_sobol_create((unsigned) n) : NULL;
191 nlopt_sobol_skip(d->s, (unsigned) d->N, d->ps + 1);
193 /* generate initial points randomly, plus starting guess x */
194 memcpy(d->ps + 1, x, sizeof(double) * n);
195 d->ps[0] = f(n, x, NULL, f_data);
197 if (!rb_tree_insert(&d->t, d->ps)) return NLOPT_OUT_OF_MEMORY;
198 if (d->ps[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
199 if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
200 if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
201 for (i = 1; i < d->N; ++i) {
202 double *k = d->ps + i*(n+1);
204 nlopt_sobol_next(d->s, k + 1, lb, ub);
207 for (j = 0; j < n; ++j)
208 k[1 + j] = nlopt_urand(lb[j], ub[j]);
210 k[0] = f(n, k + 1, NULL, f_data);
212 if (!rb_tree_insert(&d->t, k)) return NLOPT_OUT_OF_MEMORY;
213 if (k[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
214 if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
215 if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
218 return NLOPT_SUCCESS;;
221 nlopt_result crs_minimize(int n, nlopt_func f, void *f_data,
222 const double *lb, const double *ub, /* bounds */
223 double *x, /* in: initial guess, out: minimizer */
225 nlopt_stopping *stop,
226 int lds) /* random or low-discrepancy seq. (lds) */
232 ret = crs_init(&d, n, x, lb, ub, stop, f, f_data, lds);
233 if (ret < 0) return ret;
235 best = rb_tree_min(&d.t);
237 memcpy(x, best->k + 1, sizeof(double) * n);
239 while (ret == NLOPT_SUCCESS) {
240 if (NLOPT_SUCCESS == (ret = crs_trial(&d))) {
241 best = rb_tree_min(&d.t);
242 if (best->k[0] < *minf) {
243 if (best->k[0] < stop->minf_max)
244 ret = NLOPT_MINF_MAX_REACHED;
245 else if (nlopt_stop_f(stop, best->k[0], *minf))
246 ret = NLOPT_FTOL_REACHED;
247 else if (nlopt_stop_x(stop, best->k + 1, x))
248 ret = NLOPT_XTOL_REACHED;
250 memcpy(x, best->k + 1, sizeof(double) * n);
252 if (ret != NLOPT_SUCCESS) {
253 if (nlopt_stop_evals(stop))
254 ret = NLOPT_MAXEVAL_REACHED;
255 else if (nlopt_stop_time(stop))
256 ret = NLOPT_MAXTIME_REACHED;