5 #include "nlopt-util.h"
11 #define MIN(a,b) ((a) < (b) ? (a) : (b))
12 #define MAX(a,b) ((a) > (b) ? (a) : (b))
14 /***************************************************************************/
15 /* basic data structure:
17 * a hyper-rectangle is stored as an array of length L = 2n+3, where [1]
18 * is the value (f) of the function at the center, [0] is the "size"
19 * measure (d) of the rectangle, [3..n+2] are the coordinates of the
20 * center (c), [n+3..2n+2] are the widths of the sides (w), and [2]
21 * is an "age" measure for tie-breaking purposes.
23 * we store the hyper-rectangles in a red-black tree, sorted by (d,f)
24 * in lexographic order, to allow us to perform quick convex-hull
25 * calculations (in the future, we might make this data structure
26 * more sophisticated based on the dynamic convex-hull literature).
28 * n > 0 always, of course.
31 /* parameters of the search algorithm and various information that
32 needs to be passed around */
34 int n; /* dimension */
35 int L; /* size of each rectangle (2n+3) */
36 double magic_eps; /* Jones' epsilon parameter (1e-4 is recommended) */
37 int which_diam; /* which measure of hyper-rectangle diam to use:
38 0 = Jones, 1 = Gablonsky */
39 int which_div; /* which way to divide rects:
40 0: orig. Jones (divide all longest sides)
41 1: Gablonsky (cubes divide all, rects longest)
42 2: Jones Encyc. Opt.: pick random longest side */
43 int which_opt; /* which rects are considered "potentially optimal"
44 0: Jones (all pts on cvx hull, even equal pts)
45 1: Gablonsky DIRECT-L (pick one pt, if equal pts)
46 2: ~ 1, but pick points randomly if equal pts
47 ... 2 seems to suck compared to just picking oldest pt */
49 const double *lb, *ub;
50 nlopt_stopping *stop; /* stopping criteria */
51 nlopt_func f; void *f_data;
52 double *work; /* workspace, of length >= 2*n */
53 int *iwork; /* workspace, length >= n */
54 double fmin, *xmin; /* minimum so far */
56 /* red-black tree of hyperrects, sorted by (d,f,age) in
57 lexographical order */
59 int age; /* age for next new rect */
60 double **hull; /* array to store convex hull */
61 int hull_len; /* allocated length of hull array */
64 /***************************************************************************/
66 /* Evaluate the "diameter" (d) of a rectangle of widths w[n]
68 We round the result to single precision, which should be plenty for
69 the use we put the diameter to (rect sorting), to allow our
70 performance hack in convex_hull to work (in the Jones and Gablonsky
71 DIRECT algorithms, all of the rects fall into a few diameter
72 values, and we don't want rounding error to spoil this) */
73 static double rect_diameter(int n, const double *w, const params *p)
76 if (p->which_diam == 0) { /* Jones measure */
78 for (i = 0; i < n; ++i)
80 /* distance from center to a vertex */
81 return ((float) (sqrt(sum) * 0.5));
83 else { /* Gablonsky measure */
85 for (i = 0; i < n; ++i)
88 /* half-width of longest side */
89 return ((float) (maxw * 0.5));
93 #define ALLOC_RECT(rect, L) if (!(rect = (double*) malloc(sizeof(double)*(L)))) return NLOPT_OUT_OF_MEMORY
95 static double *fv_qsort = 0;
96 static int sort_fv_compare(const void *a_, const void *b_)
98 int a = *((const int *) a_), b = *((const int *) b_);
99 double fa = MIN(fv_qsort[2*a], fv_qsort[2*a+1]);
100 double fb = MIN(fv_qsort[2*b], fv_qsort[2*b+1]);
108 static void sort_fv(int n, double *fv, int *isort)
111 for (i = 0; i < n; ++i) isort[i] = i;
112 fv_qsort = fv; /* not re-entrant, sigh... */
113 qsort(isort, (unsigned) n, sizeof(int), sort_fv_compare);
117 static double function_eval(const double *x, params *p) {
118 double f = p->f(p->n, x, NULL, p->f_data);
121 memcpy(p->xmin, x, sizeof(double) * p->n);
126 #define FUNCTION_EVAL(fv,x,p,freeonerr) fv = function_eval(x, p); if (p->fmin < p->stop->fmin_max) { free(freeonerr); return NLOPT_FMIN_MAX_REACHED; } else if (nlopt_stop_evals((p)->stop)) { free(freeonerr); return NLOPT_MAXEVAL_REACHED; } else if (nlopt_stop_time((p)->stop)) { free(freeonerr); return NLOPT_MAXTIME_REACHED; }
128 #define THIRD (0.3333333333333333333333)
130 #define EQUAL_SIDE_TOL 5e-2 /* tolerance to equate side sizes */
132 /* divide rectangle idiv in the list p->rects */
133 static nlopt_result divide_rect(double *rdiv, params *p)
136 const const int n = p->n;
138 double *c = rdiv + 3; /* center of rect to divide */
139 double *w = c + n; /* widths of rect to divide */
141 int imax = 0, nlongest = 0;
144 for (i = 1; i < n; ++i)
147 for (i = 0; i < n; ++i)
148 if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL)
150 if (p->which_div == 1 || (p->which_div == 0 && nlongest == n)) {
151 /* trisect all longest sides, in increasing order of the average
152 function value along that direction */
153 double *fv = p->work;
154 int *isort = p->iwork;
155 for (i = 0; i < n; ++i) {
156 if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL) {
158 c[i] = csave - w[i] * THIRD;
159 FUNCTION_EVAL(fv[2*i], c, p, 0);
160 c[i] = csave + w[i] * THIRD;
161 FUNCTION_EVAL(fv[2*i+1], c, p, 0);
165 fv[2*i] = fv[2*i+1] = HUGE_VAL;
168 sort_fv(n, fv, isort);
169 if (!(node = rb_tree_find(&p->rtree, rdiv)))
170 return NLOPT_FAILURE;
171 for (i = 0; i < nlongest; ++i) {
173 w[isort[i]] *= THIRD;
174 rdiv[0] = rect_diameter(n, w, p);
176 node = rb_tree_resort(&p->rtree, node);
177 for (k = 0; k <= 1; ++k) {
180 memcpy(rnew, rdiv, sizeof(double) * L);
181 rnew[3 + isort[i]] += w[isort[i]] * (2*k-1);
182 rnew[1] = fv[2*isort[i]+k];
184 if (!rb_tree_insert(&p->rtree, rnew)) {
186 return NLOPT_OUT_OF_MEMORY;
193 if (nlongest > 1 && p->which_div == 2) {
194 /* randomly choose longest side */
195 i = nlopt_iurand(nlongest);
196 for (k = 0; k < n; ++k)
197 if (wmax - w[k] <= wmax * EQUAL_SIDE_TOL) {
198 if (!i) { i = k; break; }
203 i = imax; /* trisect longest side */
204 if (!(node = rb_tree_find(&p->rtree, rdiv)))
205 return NLOPT_FAILURE;
207 rdiv[0] = rect_diameter(n, w, p);
209 node = rb_tree_resort(&p->rtree, node);
210 for (k = 0; k <= 1; ++k) {
213 memcpy(rnew, rdiv, sizeof(double) * L);
214 rnew[3 + i] += w[i] * (2*k-1);
215 FUNCTION_EVAL(rnew[1], rnew + 3, p, rnew);
217 if (!rb_tree_insert(&p->rtree, rnew)) {
219 return NLOPT_OUT_OF_MEMORY;
223 return NLOPT_SUCCESS;
226 /***************************************************************************/
227 /* Convex hull algorithm, used later to find the potentially optimal
228 points. What we really have in DIRECT is a "dynamic convex hull"
229 problem, since we are dynamically adding/removing points and
230 updating the hull, but I haven't implemented any of the fancy
231 algorithms for this problem yet. */
233 /* Find the lower convex hull of a set of points (x,y) stored in a rb-tree
234 of pointers to {x,y} arrays sorted in lexographic order by (x,y).
236 Unlike standard convex hulls, we allow redundant points on the hull,
237 and even allow duplicate points if allow_dups is nonzero.
239 The return value is the number of points in the hull, with pointers
240 stored in hull[i] (should be an array of length >= t->N).
242 static int convex_hull(rb_tree *t, double **hull, int allow_dups)
246 double xmin, xmax, yminmin, ymaxmin;
249 /* Monotone chain algorithm [Andrew, 1979]. */
253 nmax = rb_tree_max(t);
260 do { /* include any duplicate points at (xmin,yminmin) */
261 hull[nhull++] = n->k;
263 } while (n && n->k[0] == xmin && n->k[1] == yminmin);
265 hull[nhull++] = n->k;
267 if (xmin == xmax) return nhull;
269 /* set nmax = min mode with x == xmax */
271 while (nmax->k[0] == xmax)
272 nmax = rb_tree_pred(nmax); /* non-NULL since xmin != xmax */
273 nmax = rb_tree_succ(nmax);
275 /* performance hack (see also below) */
278 kshift[0] = xmax * (1 - 1e-13);
279 kshift[1] = -HUGE_VAL;
280 nmax = rb_tree_find_gt(t, kshift); /* non-NULL since xmin != xmax */
284 ymaxmin = nmax->k[1];
285 minslope = (ymaxmin - yminmin) / (xmax - xmin);
287 /* set n = first node with x != xmin */
289 while (n->k[0] == xmin)
290 n = rb_tree_succ(n); /* non-NULL since xmin != xmax */
292 /* performance hack (see also below) */
295 kshift[0] = xmin * (1 + 1e-13);
296 kshift[1] = -HUGE_VAL;
297 n = rb_tree_find_gt(t, kshift); /* non-NULL since xmin != xmax */
301 for (; n != nmax; n = rb_tree_succ(n)) {
303 if (k[1] > yminmin + (k[0] - xmin) * minslope)
306 /* performance hack: most of the points in DIRECT lie along
307 vertical lines at a few x values, and we can exploit this */
308 if (nhull && k[0] == hull[nhull - 1][0]) { /* x == previous x */
309 if (k[1] > hull[nhull - 1][1]) {
311 /* because of the round to float in rect_diameter, above,
312 it shouldn't be possible for two diameters (x values)
313 to have a fractional difference < 1e-13. Note
314 that k[0] > 0 always in DIRECT */
315 kshift[0] = k[0] * (1 + 1e-13);
316 kshift[1] = -HUGE_VAL;
317 n = rb_tree_pred(rb_tree_find_gt(t, kshift));
320 else { /* equal y values, add to hull */
327 /* remove points until we are making a "left turn" to k */
329 double *t1 = hull[nhull - 1], *t2;
331 /* because we allow equal points in our hull, we have
332 to modify the standard convex-hull algorithm slightly:
333 we need to look backwards in the hull list until we
334 find a point t2 != t1 */
338 } while (it2 >= 0 && t2[0] == t1[0] && t2[1] == t1[1]);
341 /* cross product (t1-t2) x (k-t2) > 0 for a left turn: */
342 if ((t1[0]-t2[0]) * (k[1]-t2[1])
343 - (t1[1]-t2[1]) * (k[0]-t2[0]) >= 0)
351 do { /* include any duplicate points at (xmax,ymaxmin) */
352 hull[nhull++] = nmax->k;
353 nmax = rb_tree_succ(nmax);
354 } while (nmax && nmax->k[0] == xmax && n->k[1] == ymaxmin);
356 hull[nhull++] = nmax->k;
361 /***************************************************************************/
363 static int small(double *w, params *p)
366 for (i = 0; i < p->n; ++i)
367 if (w[i] > p->stop->xtol_abs[i] &&
368 w[i] > (p->ub[i] - p->lb[i]) * p->stop->xtol_rel)
373 static nlopt_result divide_good_rects(params *p)
377 int nhull, i, xtol_reached = 1, divided_some = 0;
378 double magic_eps = p->magic_eps;
380 if (p->hull_len < p->rtree.N) {
381 p->hull_len += p->rtree.N;
382 p->hull = (double **) realloc(p->hull, sizeof(double*)*p->hull_len);
383 if (!p->hull) return NLOPT_OUT_OF_MEMORY;
385 nhull = convex_hull(&p->rtree, hull = p->hull, p->which_opt != 1);
387 for (i = 0; i < nhull; ++i) {
388 double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
391 /* find unequal points before (im) and after (ip) to get slope */
392 for (im = i-1; im >= 0 && hull[im][0] == hull[i][0]; --im);
393 for (ip = i+1; ip < nhull && hull[ip][0] == hull[i][0]; ++ip);
396 K1 = (hull[i][1] - hull[im][1]) / (hull[i][0] - hull[im][0]);
398 K1 = (hull[i][1] - hull[ip][1]) / (hull[i][0] - hull[ip][0]);
400 if (hull[i][1] - K * hull[i][0]
401 <= p->fmin - magic_eps * fabs(p->fmin) || ip == nhull) {
402 /* "potentially optimal" rectangle, so subdivide */
403 nlopt_result ret = divide_rect(hull[i], p);
405 if (ret != NLOPT_SUCCESS) return ret;
406 xtol_reached = xtol_reached && small(hull[i] + 3+n, p);
409 /* for the DIRECT-L variant, we only divide one rectangle out
410 of all points with equal diameter and function values
411 ... note that for p->which_opt == 1, i == ip-1 should be a no-op
412 anyway, since we set allow_dups=0 in convex_hull above */
413 if (p->which_opt == 1)
414 i = ip - 1; /* skip to next unequal point for next iteration */
415 else if (p->which_opt == 2) /* like DIRECT-L but randomized */
416 i += nlopt_iurand(ip - i); /* possibly do another equal pt */
419 if (magic_eps != 0) {
421 goto divisions; /* try again */
423 else { /* WTF? divide largest rectangle with smallest f */
424 /* (note that this code actually gets called from time
425 to time, and the heuristic here seems to work well,
426 but I don't recall this situation being discussed in
428 rb_node *max = rb_tree_max(&p->rtree);
430 double wmax = max->k[0];
431 do { /* note: this loop is O(N) worst-case time */
433 pred = rb_tree_pred(max);
434 } while (pred && pred->k[0] == wmax);
435 return divide_rect(max->k, p);
438 return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
441 /***************************************************************************/
443 /* lexographic sort order (d,f,age) of hyper-rects, for red-black tree */
444 int cdirect_hyperrect_compare(double *a, double *b)
446 if (a[0] < b[0]) return -1;
447 if (a[0] > b[0]) return +1;
448 if (a[1] < b[1]) return -1;
449 if (a[1] > b[1]) return +1;
450 if (a[2] < b[2]) return -1;
451 if (a[2] > b[2]) return +1;
452 return (int) (a - b); /* tie-breaker, shouldn't be needed */
455 /***************************************************************************/
457 nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
458 const double *lb, const double *ub,
461 nlopt_stopping *stop,
462 double magic_eps, int which_alg)
467 nlopt_result ret = NLOPT_OUT_OF_MEMORY;
469 p.magic_eps = magic_eps;
470 p.which_diam = which_alg % 3;
471 p.which_div = (which_alg / 3) % 3;
472 p.which_opt = (which_alg / (3*3)) % 3;
473 p.lb = lb; p.ub = ub;
486 rb_tree_init(&p.rtree, cdirect_hyperrect_compare);
488 p.work = (double *) malloc(sizeof(double) * (2*n));
489 if (!p.work) goto done;
490 p.iwork = (int *) malloc(sizeof(int) * n);
491 if (!p.iwork) goto done;
492 p.hull_len = 128; /* start with a reasonable number */
493 p.hull = (double **) malloc(sizeof(double *) * p.hull_len);
494 if (!p.hull) goto done;
496 if (!(rnew = (double *) malloc(sizeof(double) * p.L))) goto done;
497 for (i = 0; i < n; ++i) {
498 rnew[3+i] = 0.5 * (lb[i] + ub[i]);
499 rnew[3+n+i] = ub[i] - lb[i];
501 rnew[0] = rect_diameter(n, rnew+3+n, &p);
502 rnew[1] = function_eval(rnew+3, &p);
504 if (!rb_tree_insert(&p.rtree, rnew)) {
509 ret = divide_rect(rnew, &p);
510 if (ret != NLOPT_SUCCESS) goto done;
513 double fmin0 = p.fmin;
514 ret = divide_good_rects(&p);
515 if (ret != NLOPT_SUCCESS) goto done;
516 if (p.fmin < fmin0 && nlopt_stop_f(p.stop, p.fmin, fmin0)) {
517 ret = NLOPT_FTOL_REACHED;
523 rb_tree_destroy_with_keys(&p.rtree);
532 /* in the conventional DIRECT-type algorithm, we first rescale our
533 coordinates to a unit hypercube ... we do this simply by
534 wrapping cdirect() around cdirect_unscaled(). */
536 double cdirect_uf(int n, const double *xu, double *grad, void *d_)
538 cdirect_uf_data *d = (cdirect_uf_data *) d_;
541 for (i = 0; i < n; ++i)
542 d->x[i] = d->lb[i] + xu[i] * (d->ub[i] - d->lb[i]);
543 f = d->f(n, d->x, grad, d->f_data);
545 for (i = 0; i < n; ++i)
546 grad[i] *= d->ub[i] - d->lb[i];
550 nlopt_result cdirect(int n, nlopt_func f, void *f_data,
551 const double *lb, const double *ub,
554 nlopt_stopping *stop,
555 double magic_eps, int which_alg)
559 const double *xtol_abs_save;
562 d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
563 d.x = (double *) malloc(sizeof(double) * n*4);
564 if (!d.x) return NLOPT_OUT_OF_MEMORY;
566 for (i = 0; i < n; ++i) {
567 x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
570 d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
572 xtol_abs_save = stop->xtol_abs;
573 stop->xtol_abs = d.x + 3*n;
574 ret = cdirect_unscaled(n, cdirect_uf, &d, d.x+n, d.x+2*n, x, fmin, stop,
575 magic_eps, which_alg);
576 stop->xtol_abs = xtol_abs_save;
577 for (i = 0; i < n; ++i)
578 x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);