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.
27 #include "nlopt-util.h"
32 /* Hybrid algorithm, inspired by DIRECT, that uses another local
33 * optimization algorithm within each rectangle, and then looks
34 * in the largest remaining rectangle (breaking ties by minimum
35 * function value and then by age.
37 * Each hyperrect is represented by an array of length 3*n+3 consisting
38 * of (d, -f, -a, x, c, w), where d=diameter, f=f(x), a=age, x=local optimum
43 int n; /* dimension */
45 const double *lb, *ub;
46 nlopt_stopping *stop; /* stopping criteria */
47 nlopt_func f; void *f_data;
48 double minf, *xmin; /* min so far */
49 rb_tree rtree; /* red-black tree of rects, sorted by (d,-f,-a) */
50 int age; /* age for next new rect */
51 double *work; /* workspace of length >= 2*n */
53 nlopt_algorithm local_alg; /* local search algorithm */
54 int local_maxeval; /* max # local iterations (0 if unlimited) */
56 int randomized_div; /* 1 to use randomized division algorithm */
59 #define MIN(a,b) ((a) < (b) ? (a) : (b))
61 #define THIRD (0.3333333333333333333333) /* 1/3 */
63 /************************************************************************/
65 static double fcount(int n, const double *x, double *grad, void *p_)
67 params *p = (params *) p_;
69 return p->f(n, x, grad, p->f_data);
72 static nlopt_result optimize_rect(double *r, params *p)
75 double *lb = p->work, *ub = lb + n;
76 double *x = r + 3, *c = x + n, *w = c + n;
77 double t = nlopt_seconds();
79 nlopt_stopping *stop = p->stop;
82 if (stop->maxeval > 0 &&
83 stop->nevals >= stop->maxeval) return NLOPT_MAXEVAL_REACHED;
84 if (stop->maxtime > 0 &&
85 t - stop->start >= stop->maxtime) return NLOPT_MAXTIME_REACHED;
87 for (i = 0; i < n; ++i) {
88 lb[i] = c[i] - 0.5 * w[i];
89 ub[i] = c[i] + 0.5 * w[i];
91 ret = nlopt_minimize(p->local_alg, n, fcount, p,
93 stop->minf_max, stop->ftol_rel, stop->ftol_abs,
94 stop->xtol_rel, stop->xtol_abs,
95 p->local_maxeval > 0 ?
97 stop->maxeval - stop->nevals)
98 : stop->maxeval - stop->nevals,
99 stop->maxtime - (t - stop->start));
102 if (minf < p->minf) {
104 memcpy(p->xmin, x, sizeof(double) * n);
105 if (ret == NLOPT_MINF_MAX_REACHED) return ret;
107 return NLOPT_SUCCESS;
112 /* given a hyperrect r, randomize the starting guess within the middle
113 third of the box (don't guess too close to edges) */
114 static void randomize_x(int n, double *r)
116 double *x = r + 3, *c = x + n, *w = c + n;
118 for (i = 0; i < n; ++i)
119 x[i] = nlopt_urand(c[i] - w[i]*(0.5*THIRD),
120 c[i] + w[i]*(0.5*THIRD));
123 /************************************************************************/
125 static double longest(int n, const double *w)
127 double wmax = w[n-1];
128 for (n = n-2; n >= 0; n--) if (w[n] > wmax) wmax = w[n];
132 #define EQUAL_SIDE_TOL 5e-2 /* tolerance to equate side sizes */
134 static nlopt_result divide_largest(params *p)
138 rb_node *node = rb_tree_max(&p->rtree); /* just using it as a heap */
139 double minf_start = p->minf;
140 double *r = node->k, *rnew = NULL;
141 double *x = r + 3, *c = x + n, *w = c + n;
142 const double *lb = p->lb, *ub = p->ub;
147 /* printf("rect:, %d, %g, %g, %g, %g\n", p->stop->nevals, c[0], c[1], w[0], w[1]); */
150 for (i = 0; i < n; ++i)
151 if (w[i] > p->stop->xtol_rel * (ub[i] - lb[i])
152 && w[i] > p->stop->xtol_abs[i])
154 if (i == n) return NLOPT_XTOL_REACHED;
156 if (p->randomized_div) { /* randomly pick among ~largest sides */
158 wmax = longest(n, w);
159 for (i = 0; i < n; ++i)
160 if (wmax - w[i] < EQUAL_SIDE_TOL * wmax) ++nlongest;
161 i = 1 + nlopt_iurand(nlongest);
162 for (idiv = 0; idiv < n; ++idiv) {
163 if (wmax - w[idiv] < EQUAL_SIDE_TOL * wmax) --i;
167 else { /* just pick first largest side */
169 for (i = 1; i < n; ++i) if (w[i] > wmax) wmax = w[idiv = i];
172 if (fabs(x[idiv] - c[idiv]) > (0.5 * THIRD) * w[idiv]) { /* bisect */
173 double deltac = (x[idiv] > c[idiv] ? 0.25 : -0.25) * w[idiv];
176 r[0] = longest(n, w); /* new diameter */
177 /* r[1] unchanged since still contains local optimum x */
179 node = rb_tree_resort(&p->rtree, node);
181 rnew = (double *) malloc(sizeof(double) * L);
182 if (!rnew) return NLOPT_OUT_OF_MEMORY;
183 memcpy(rnew, r, sizeof(double) * L);
185 rnew[3+n+idiv] -= deltac*2;
186 if (p->randomized_div)
187 randomize_x(n, rnew);
189 memcpy(rnew+3, rnew+3+n, sizeof(double) * n); /* x = c */
190 ret = optimize_rect(rnew, p);
191 if (ret != NLOPT_SUCCESS) { free(rnew); return ret; }
192 if (!rb_tree_insert(&p->rtree, rnew)) {
193 free(rnew); return NLOPT_OUT_OF_MEMORY;
198 r[0] = longest(n, w);
199 /* r[1] unchanged since still contains local optimum x */
201 node = rb_tree_resort(&p->rtree, node);
203 for (i = -1; i <= +1; i += 2) {
204 rnew = (double *) malloc(sizeof(double) * L);
205 if (!rnew) return NLOPT_OUT_OF_MEMORY;
206 memcpy(rnew, r, sizeof(double) * L);
208 rnew[3+n+idiv] += w[i] * i;
209 if (p->randomized_div)
210 randomize_x(n, rnew);
212 memcpy(rnew+3, rnew+3+n, sizeof(double) * n); /* x = c */
213 ret = optimize_rect(rnew, p);
214 if (ret != NLOPT_SUCCESS) { free(rnew); return ret; }
215 if (!rb_tree_insert(&p->rtree, rnew)) {
216 free(rnew); return NLOPT_OUT_OF_MEMORY;
220 if (p->minf < minf_start && nlopt_stop_f(p->stop, p->minf, minf_start))
221 return NLOPT_FTOL_REACHED;
222 return NLOPT_SUCCESS;
225 /************************************************************************/
227 nlopt_result cdirect_hybrid_unscaled(int n, nlopt_func f, void *f_data,
228 const double *lb, const double *ub,
231 nlopt_stopping *stop,
232 nlopt_algorithm local_alg,
239 nlopt_result ret = NLOPT_OUT_OF_MEMORY;
243 p.lb = lb; p.ub = ub;
251 p.local_alg = local_alg;
252 p.local_maxeval = local_maxeval;
253 p.randomized_div = randomized_div;
255 rb_tree_init(&p.rtree, cdirect_hyperrect_compare);
256 p.work = (double *) malloc(sizeof(double) * (2*n));
257 if (!p.work) goto done;
259 if (!(rnew = (double *) malloc(sizeof(double) * p.L))) goto done;
260 for (i = 0; i < n; ++i) {
261 rnew[3+i] = rnew[3+n+i] = 0.5 * (lb[i] + ub[i]);
262 rnew[3+2*n+i] = ub[i] - lb[i];
264 rnew[0] = longest(n, rnew+2*n);
266 ret = optimize_rect(rnew, &p);
267 if (ret != NLOPT_SUCCESS) { free(rnew); goto done; }
268 if (!rb_tree_insert(&p.rtree, rnew)) { free(rnew); goto done; }
271 ret = divide_largest(&p);
272 } while (ret == NLOPT_SUCCESS);
275 rb_tree_destroy_with_keys(&p.rtree);
282 /* rescaled to unit hypercube so that all x[i] are weighted equally */
283 nlopt_result cdirect_hybrid(int n, nlopt_func f, void *f_data,
284 const double *lb, const double *ub,
287 nlopt_stopping *stop,
288 nlopt_algorithm local_alg,
294 const double *xtol_abs_save;
297 d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
298 d.x = (double *) malloc(sizeof(double) * n*4);
299 if (!d.x) return NLOPT_OUT_OF_MEMORY;
301 for (i = 0; i < n; ++i) {
302 x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
305 d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
307 xtol_abs_save = stop->xtol_abs;
308 stop->xtol_abs = d.x + 3*n;
309 ret = cdirect_hybrid_unscaled(n, cdirect_uf, &d, d.x+n, d.x+2*n,
310 x, minf, stop, local_alg, local_maxeval,
312 stop->xtol_abs = xtol_abs_save;
313 for (i = 0; i < n; ++i)
314 x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);