7 /* Controlled Random Search 2 (CRS2) with "local mutation", as defined
9 P. Kaelo and M. M. Ali, "Some variants of the controlled random
10 search algorithm for global optimization," J. Optim. Theory Appl.
11 130 (2), 253-264 (2006).
15 int n; /* # dimensions */
16 const double *lb, *ub;
17 nlopt_stopping *stop; /* stopping criteria */
18 nlopt_func f; void *f_data;
20 int N; /* # points in population */
21 double *ps; /* population array N x (n+1) of tuples [f(x), x] */
22 double *p; /* single point array (length n+1), for temp use */
23 rb_tree t; /* red-black tree of population, sorted by f(x) */
24 nlopt_sobol s; /* sobol data for LDS point generation, or NULL
25 to use pseudo-random numbers */
28 /* sort order in red-black tree: keys [f(x), x] are sorted by f(x) */
29 static int crs_compare(double *k1, double *k2)
31 if (*k1 < *k2) return -1;
32 if (*k1 > *k2) return +1;
33 return k1 - k2; /* tie-breaker */
36 /* set x to a random trial value, as defined by CRS:
38 where x_0 ... x_n are distinct points in the population
39 with x_0 the current best point and the other points are random,
40 and G is the centroid of x_0...x_{n-1} */
41 static void random_trial(crs_data *d, double *x, rb_node *best)
43 int n = d->n, n1 = n+1, i, k, i0, jn;
44 double *ps = d->ps, *xi;
46 /* initialize x to x_0 = best point */
47 memcpy(x, best->k + 1, sizeof(double) * n);
48 i0 = (best->k - ps) / n1;
50 jn = nlopt_iurand(n); /* which of remaining n points is "x_n",
51 i.e. which to reflect through ...
52 this is necessary since we generate
53 the remaining points in order, so
54 just picking the last point would not
57 /* use "method A" from
59 Jeffrey Scott Vitter, "An efficient algorithm for
60 sequential random sampling," ACM Trans. Math. Soft. 13 (1),
63 to randomly pick n distinct points out of the remaining N-1 (not
64 including i0!). (The same as "method S" in Knuth vol. 2.)
65 This method requires O(N) time, which is fine in our case
66 (there are better methods if n << N). */
68 int Nleft = d->N - 1, nleft = n;
69 int Nfree = Nleft - nleft;
72 double q = ((double) Nfree) / Nleft;
73 double v = nlopt_urand(0., 1.);
77 q = (q * Nfree) / Nleft;
80 if (jn-- == 0) /* point to reflect through */
81 for (k = 0; k < n; ++k) x[k] -= xi[k] * (0.5*n);
82 else /* point to include in centroid */
83 for (k = 0; k < n; ++k) x[k] += xi[k];
87 i += nlopt_iurand(Nleft); i += i == i0;
89 if (jn-- == 0) /* point to reflect through */
90 for (k = 0; k < n; ++k) x[k] -= xi[k] * (0.5*n);
91 else /* point to include in centroid */
92 for (k = 0; k < n; ++k) x[k] += xi[k];
94 for (k = 0; k < n; ++k) {
95 x[k] *= 2.0 / n; /* renormalize */
96 if (x[k] > d->ub[k]) x[k] = d->ub[k];
97 else if (x[k] < d->lb[k]) x[k] = d->lb[k];
101 #define NUM_MUTATION 1 /* # "local mutation" steps to try if trial fails */
103 static nlopt_result crs_trial(crs_data *d)
105 rb_node *best = rb_tree_min(&d->t);
106 rb_node *worst = rb_tree_max(&d->t);
107 int mutation = NUM_MUTATION;
109 random_trial(d, d->p + 1, best);
111 d->p[0] = d->f(n, d->p + 1, NULL, d->f_data);
113 if (d->p[0] < worst->k[0]) break;
114 if (nlopt_stop_evals(d->stop)) return NLOPT_MAXEVAL_REACHED;
115 if (nlopt_stop_time(d->stop)) return NLOPT_MAXTIME_REACHED;
117 for (i = 0; i < n; ++i) {
118 double w = nlopt_urand(0.,1.);
119 d->p[1+i] = best->k[1+i] * (1 + w) - w * d->p[1+i];
120 if (d->p[1+i] > d->ub[i]) d->p[1+i] = d->ub[i];
121 else if (d->p[1+i] < d->lb[i]) d->p[1+i] = d->lb[i];
126 random_trial(d, d->p + 1, best);
127 mutation = NUM_MUTATION;
130 memcpy(worst->k, d->p, sizeof(double) * (n+1));
131 rb_tree_resort(&d->t, worst);
132 return NLOPT_SUCCESS;
135 static void crs_destroy(crs_data *d)
137 nlopt_sobol_destroy(d->s);
138 rb_tree_destroy(&d->t);
142 static nlopt_result crs_init(crs_data *d, int n, const double *x,
143 const double *lb, const double *ub,
144 nlopt_stopping *stop, nlopt_func f, void *f_data,
149 /* TODO: how should we set the initial population size?
150 the Kaelo and Ali paper suggests 10*(n+1), but should
151 we add more random points if maxeval is large, or... ? */
152 d->N = 10 * (n + 1); /* heuristic initial population size */
156 d->f = f; d->f_data = f_data;
157 d->ub = ub; d->lb = lb;
158 d->ps = (double *) malloc(sizeof(double) * (n + 1) * (d->N + 1));
159 if (!d->ps) return NLOPT_OUT_OF_MEMORY;
160 d->p = d->ps + d->N * (n+1);
161 rb_tree_init(&d->t, crs_compare);
163 /* we can either use pseudorandom points, as in the original CRS
164 algorithm, or use a low-discrepancy Sobol' sequence ... I tried
165 the latter, however, and it doesn't seem to help, probably
166 because we are only generating a small number of random points
168 d->s = lds ? nlopt_sobol_create((unsigned) n) : NULL;
169 nlopt_sobol_skip(d->s, (unsigned) d->N, d->ps + 1);
171 /* generate initial points randomly, plus starting guess x */
172 memcpy(d->ps + 1, x, sizeof(double) * n);
173 d->ps[0] = f(n, x, NULL, f_data);
175 if (!rb_tree_insert(&d->t, d->ps)) return NLOPT_OUT_OF_MEMORY;
176 if (d->ps[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
177 if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
178 if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
179 for (i = 1; i < d->N; ++i) {
180 double *k = d->ps + i*(n+1);
182 nlopt_sobol_next(d->s, k + 1, lb, ub);
185 for (j = 0; j < n; ++j)
186 k[1 + j] = nlopt_urand(lb[j], ub[j]);
188 k[0] = f(n, k + 1, NULL, f_data);
190 if (!rb_tree_insert(&d->t, k)) return NLOPT_OUT_OF_MEMORY;
191 if (k[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
192 if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
193 if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
196 return NLOPT_SUCCESS;;
199 nlopt_result crs_minimize(int n, nlopt_func f, void *f_data,
200 const double *lb, const double *ub, /* bounds */
201 double *x, /* in: initial guess, out: minimizer */
203 nlopt_stopping *stop,
204 int lds) /* random or low-discrepancy seq. (lds) */
210 ret = crs_init(&d, n, x, lb, ub, stop, f, f_data, lds);
211 if (ret < 0) return ret;
213 best = rb_tree_min(&d.t);
215 memcpy(x, best->k + 1, sizeof(double) * n);
217 while (ret == NLOPT_SUCCESS) {
218 if (NLOPT_SUCCESS == (ret = crs_trial(&d))) {
219 best = rb_tree_min(&d.t);
220 if (best->k[0] < *minf) {
221 if (best->k[0] < stop->minf_max)
222 ret = NLOPT_MINF_MAX_REACHED;
223 else if (nlopt_stop_f(stop, best->k[0], *minf))
224 ret = NLOPT_FTOL_REACHED;
225 else if (nlopt_stop_x(stop, best->k + 1, x))
226 ret = NLOPT_XTOL_REACHED;
228 memcpy(x, best->k + 1, sizeof(double) * n);
230 if (ret != NLOPT_SUCCESS) {
231 if (nlopt_stop_evals(stop))
232 ret = NLOPT_MAXEVAL_REACHED;
233 else if (nlopt_stop_time(stop))
234 ret = NLOPT_MAXTIME_REACHED;