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+2, where [1]
18 * is the value (f) of the function at the center, [0] is the "size"
19 * measure (d) of the rectangle, [2..n+1] are the coordinates of the
20 * center (c), and [n+2..2n+1] are the widths of the sides (w).
22 * we store the hyper-rectangles in a red-black tree, sorted by (d,f)
23 * in lexographic order, to allow us to perform quick convex-hull
24 * calculations (in the future, we might make this data structure
25 * more sophisticated based on the dynamic convex-hull literature).
30 #define RECT_LEN(n) (2*(n)+2) /* number of double values in a hyperrect */
32 /* parameters of the search algorithm and various information that
33 needs to be passed around */
35 int n; /* dimension */
36 int L; /* RECT_LEN(n) */
37 double magic_eps; /* Jones' epsilon parameter (1e-4 is recommended) */
38 int which_diam; /* which measure of hyper-rectangle diam to use:
39 0 = Jones, 1 = Gablonsky */
40 int which_div; /* which way to divide rects:
41 0: Gablonsky (cubes divide all, rects longest)
42 1: orig. Jones (divide all longest sides)
43 2: Jones Encyc. Opt.: pick random longest side */
44 const double *lb, *ub;
45 nlopt_stopping *stop; /* stopping criteria */
46 nlopt_func f; void *f_data;
47 double *work; /* workspace, of length >= 2*n */
48 int *iwork; /* workspace, length >= n */
49 double fmin, *xmin; /* minimum so far */
51 /* red-black tree of hyperrects, sorted by (d,f) in
52 lexographical order */
54 double **hull; /* array to store convex hull */
55 int hull_len; /* allocated length of hull array */
58 /***************************************************************************/
60 /* Evaluate the "diameter" (d) of a rectangle of widths w[n]
62 We round the result to single precision, which should be plenty for
63 the use we put the diameter to (rect sorting), to allow our
64 performance hack in convex_hull to work (in the Jones and Gablonsky
65 DIRECT algorithms, all of the rects fall into a few diameter
66 values, and we don't want rounding error to spoil this) */
67 static double rect_diameter(int n, const double *w, const params *p)
70 if (p->which_diam == 0) { /* Jones measure */
72 for (i = 0; i < n; ++i)
74 /* distance from center to a vertex */
75 return ((float) (sqrt(sum) * 0.5));
77 else { /* Gablonsky measure */
79 for (i = 0; i < n; ++i)
82 /* half-width of longest side */
83 return ((float) (maxw * 0.5));
87 #define ALLOC_RECT(rect, L) if (!(rect = (double*) malloc(sizeof(double)*(L)))) return NLOPT_OUT_OF_MEMORY
89 static double *fv_qsort = 0;
90 static int sort_fv_compare(const void *a_, const void *b_)
92 int a = *((const int *) a_), b = *((const int *) b_);
93 double fa = MIN(fv_qsort[2*a], fv_qsort[2*a+1]);
94 double fb = MIN(fv_qsort[2*b], fv_qsort[2*b+1]);
102 static void sort_fv(int n, double *fv, int *isort)
105 for (i = 0; i < n; ++i) isort[i] = i;
106 fv_qsort = fv; /* not re-entrant, sigh... */
107 qsort(isort, (unsigned) n, sizeof(int), sort_fv_compare);
111 static double function_eval(const double *x, params *p) {
112 double f = p->f(p->n, x, NULL, p->f_data);
115 memcpy(p->xmin, x, sizeof(double) * p->n);
120 #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; }
122 #define THIRD (0.3333333333333333333333)
124 #define EQUAL_SIDE_TOL 5e-2 /* tolerance to equate side sizes */
126 /* divide rectangle idiv in the list p->rects */
127 static nlopt_result divide_rect(double *rdiv, params *p)
130 const const int n = p->n;
132 double *c = rdiv + 2; /* center of rect to divide */
133 double *w = c + n; /* widths of rect to divide */
135 int imax = 0, nlongest = 0;
138 for (i = 1; i < n; ++i)
141 for (i = 0; i < n; ++i)
142 if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL)
144 if (p->which_div == 1 || (p->which_div == 0 && nlongest == n)) {
145 /* trisect all longest sides, in increasing order of the average
146 function value along that direction */
147 double *fv = p->work;
148 int *isort = p->iwork;
149 for (i = 0; i < n; ++i) {
150 if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL) {
152 c[i] = csave - w[i] * THIRD;
153 FUNCTION_EVAL(fv[2*i], c, p, 0);
154 c[i] = csave + w[i] * THIRD;
155 FUNCTION_EVAL(fv[2*i+1], c, p, 0);
159 fv[2*i] = fv[2*i+1] = HUGE_VAL;
162 sort_fv(n, fv, isort);
163 if (!(node = rb_tree_find(&p->rtree, rdiv)))
164 return NLOPT_FAILURE;
165 for (i = 0; i < nlongest; ++i) {
167 w[isort[i]] *= THIRD;
168 rdiv[0] = rect_diameter(n, w, p);
169 node = rb_tree_resort(&p->rtree, node);
170 for (k = 0; k <= 1; ++k) {
173 memcpy(rnew, rdiv, sizeof(double) * L);
174 rnew[2 + isort[i]] += w[isort[i]] * (2*k-1);
175 rnew[1] = fv[2*isort[i]+k];
176 if (!rb_tree_insert(&p->rtree, rnew)) {
178 return NLOPT_OUT_OF_MEMORY;
185 if (nlongest > 1 && p->which_div == 2) {
186 /* randomly choose longest side */
187 i = nlopt_iurand(nlongest);
188 for (k = 0; k < n; ++k)
189 if (wmax - w[k] <= wmax * EQUAL_SIDE_TOL) {
190 if (!i) { i = k; break; }
195 i = imax; /* trisect longest side */
196 if (!(node = rb_tree_find(&p->rtree, rdiv)))
197 return NLOPT_FAILURE;
199 rdiv[0] = rect_diameter(n, w, p);
200 node = rb_tree_resort(&p->rtree, node);
201 for (k = 0; k <= 1; ++k) {
204 memcpy(rnew, rdiv, sizeof(double) * L);
205 rnew[2 + i] += w[i] * (2*k-1);
206 FUNCTION_EVAL(rnew[1], rnew + 2, p, rnew);
207 if (!rb_tree_insert(&p->rtree, rnew)) {
209 return NLOPT_OUT_OF_MEMORY;
213 return NLOPT_SUCCESS;
216 /***************************************************************************/
217 /* Convex hull algorithm, used later to find the potentially optimal
218 points. What we really have in DIRECT is a "dynamic convex hull"
219 problem, since we are dynamically adding/removing points and
220 updating the hull, but I haven't implemented any of the fancy
221 algorithms for this problem yet. */
223 /* Find the lower convex hull of a set of points (x,y) stored in a rb-tree
224 of pointers to {x,y} arrays sorted in lexographic order by (x,y).
226 Unlike standard convex hulls, we allow redundant points on the hull.
228 The return value is the number of points in the hull, with pointers
229 stored in hull[i] (should be an array of length >= t->N).
231 static int convex_hull(rb_tree *t, double **hull)
235 double xmin, xmax, yminmin, ymaxmin;
238 /* Monotone chain algorithm [Andrew, 1979]. */
242 nmax = rb_tree_max(t);
248 do { /* include any duplicate points at (xmin,yminmin) */
249 hull[nhull++] = n->k;
251 } while (n && n->k[0] == xmin && n->k[1] == yminmin);
252 if (xmin == xmax) return nhull;
254 /* set nmax = min mode with x == xmax */
256 while (nmax->k[0] == xmax)
257 nmax = rb_tree_pred(nmax); /* non-NULL since xmin != xmax */
258 nmax = rb_tree_succ(nmax);
260 /* performance hack (see also below) */
263 kshift[0] = xmax * (1 - 1e-13);
264 kshift[1] = -HUGE_VAL;
265 nmax = rb_tree_find_gt(t, kshift); /* non-NULL since xmin != xmax */
269 ymaxmin = nmax->k[1];
270 minslope = (ymaxmin - yminmin) / (xmax - xmin);
272 /* set n = first node with x != xmin */
274 while (n->k[0] == xmin)
275 n = rb_tree_succ(n); /* non-NULL since xmin != xmax */
277 /* performance hack (see also below) */
280 kshift[0] = xmin * (1 + 1e-13);
281 kshift[1] = -HUGE_VAL;
282 n = rb_tree_find_gt(t, kshift); /* non-NULL since xmin != xmax */
286 for (; n != nmax; n = rb_tree_succ(n)) {
288 if (k[1] > yminmin + (k[0] - xmin) * minslope)
291 /* performance hack: most of the points in DIRECT lie along
292 vertical lines at a few x values, and we can exploit this */
293 if (nhull && k[0] == hull[nhull - 1][0]) { /* x == previous x */
294 if (k[1] > hull[nhull - 1][1]) {
296 /* because of the round to float in rect_diameter, above,
297 it shouldn't be possible for two diameters (x values)
298 to have a fractional difference < 1e-13. Note
299 that k[0] > 0 always in DIRECT */
300 kshift[0] = k[0] * (1 + 1e-13);
301 kshift[1] = -HUGE_VAL;
302 n = rb_tree_pred(rb_tree_find_gt(t, kshift));
305 else { /* equal y values, add to hull */
311 /* remove points until we are making a "left turn" to k */
313 double *t1 = hull[nhull - 1], *t2;
315 /* because we allow equal points in our hull, we have
316 to modify the standard convex-hull algorithm slightly:
317 we need to look backwards in the hull list until we
318 find a point t2 != t1 */
322 } while (it2 >= 0 && t2[0] == t1[0] && t2[1] == t1[1]);
325 /* cross product (t1-t2) x (k-t2) > 0 for a left turn: */
326 if ((t1[0]-t2[0]) * (k[1]-t2[1])
327 - (t1[1]-t2[1]) * (k[0]-t2[0]) >= 0)
334 do { /* include any duplicate points at (xmax,ymaxmin) */
335 hull[nhull++] = nmax->k;
336 nmax = rb_tree_succ(nmax);
337 } while (nmax && nmax->k[0] == xmax && n->k[1] == ymaxmin);
342 /***************************************************************************/
344 static int small(double *w, params *p)
347 for (i = 0; i < p->n; ++i)
348 if (w[i] > p->stop->xtol_abs[i] &&
349 w[i] > (p->ub[i] - p->lb[i]) * p->stop->xtol_rel)
354 static nlopt_result divide_good_rects(params *p)
358 int nhull, i, xtol_reached = 1, divided_some = 0;
359 double magic_eps = p->magic_eps;
361 if (p->hull_len < p->rtree.N) {
362 p->hull_len += p->rtree.N;
363 p->hull = (double **) realloc(p->hull, sizeof(double*)*p->hull_len);
364 if (!p->hull) return NLOPT_OUT_OF_MEMORY;
366 nhull = convex_hull(&p->rtree, hull = p->hull);
368 for (i = 0; i < nhull; ++i) {
369 double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
371 for (im = i-1; im >= 0 && hull[im][0] == hull[i][0]; --im);
372 for (ip = i+1; ip < nhull && hull[ip][0] == hull[i][0]; ++ip);
374 K1 = (hull[i][1] - hull[im][1]) / (hull[i][0] - hull[im][0]);
376 K1 = (hull[i][1] - hull[ip][1]) / (hull[i][0] - hull[ip][0]);
378 if (hull[i][1] - K * hull[i][0]
379 <= p->fmin - magic_eps * fabs(p->fmin) || ip == nhull) {
380 /* "potentially optimal" rectangle, so subdivide */
381 nlopt_result ret = divide_rect(hull[i], p);
383 if (ret != NLOPT_SUCCESS) return ret;
384 xtol_reached = xtol_reached && small(hull[i] + 2+n, p);
388 if (magic_eps != 0) {
390 goto divisions; /* try again */
392 else { /* WTF? divide largest rectangle with smallest f */
393 /* (note that this code actually gets called from time
394 to time, and the heuristic here seems to work well,
395 but I don't recall this situation being discussed in
397 rb_node *max = rb_tree_max(&p->rtree);
399 double wmax = max->k[0];
400 do { /* note: this loop is O(N) worst-case time */
402 pred = rb_tree_pred(max);
403 } while (pred && pred->k[0] == wmax);
404 return divide_rect(max->k, p);
407 return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
410 /***************************************************************************/
412 /* lexographic sort order (d,f) of hyper-rects, for red-black tree */
413 static int hyperrect_compare(double *a, double *b)
415 if (a[0] < b[0]) return -1;
416 if (a[0] > b[0]) return +1;
417 if (a[1] < b[1]) return -1;
418 if (a[1] > b[1]) return +1;
419 return (int) (a - b); /* tie-breaker */
422 /***************************************************************************/
424 nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
425 const double *lb, const double *ub,
428 nlopt_stopping *stop,
429 double magic_eps, int which_alg)
434 nlopt_result ret = NLOPT_OUT_OF_MEMORY;
436 p.magic_eps = magic_eps;
437 p.which_diam = which_alg % 10;
438 p.which_div = (which_alg / 10) % 10;
439 p.lb = lb; p.ub = ub;
451 rb_tree_init(&p.rtree, hyperrect_compare);
453 p.work = (double *) malloc(sizeof(double) * (2*n));
454 if (!p.work) goto done;
455 p.iwork = (int *) malloc(sizeof(int) * n);
456 if (!p.iwork) goto done;
457 p.hull_len = 128; /* start with a reasonable number */
458 p.hull = (double **) malloc(sizeof(double *) * p.hull_len);
459 if (!p.hull) goto done;
461 if (!(rnew = (double *) malloc(sizeof(double) * p.L))) goto done;
462 for (i = 0; i < n; ++i) {
463 rnew[2+i] = 0.5 * (lb[i] + ub[i]);
464 rnew[2+n+i] = ub[i] - lb[i];
466 rnew[0] = rect_diameter(n, rnew+2+n, &p);
467 rnew[1] = function_eval(rnew+2, &p);
468 if (!rb_tree_insert(&p.rtree, rnew)) {
473 ret = divide_rect(rnew, &p);
474 if (ret != NLOPT_SUCCESS) goto done;
477 double fmin0 = p.fmin;
478 ret = divide_good_rects(&p);
479 if (ret != NLOPT_SUCCESS) goto done;
480 if (p.fmin < fmin0 && nlopt_stop_f(p.stop, p.fmin, fmin0)) {
481 ret = NLOPT_FTOL_REACHED;
487 rb_tree_destroy_with_keys(&p.rtree);
496 /* in the conventional DIRECT-type algorithm, we first rescale our
497 coordinates to a unit hypercube ... we do this simply by
498 wrapping cdirect() around cdirect_unscaled(). */
504 const double *lb, *ub;
506 static double uf(int n, const double *xu, double *grad, void *d_)
508 uf_data *d = (uf_data *) d_;
511 for (i = 0; i < n; ++i)
512 d->x[i] = d->lb[i] + xu[i] * (d->ub[i] - d->lb[i]);
513 f = d->f(n, d->x, grad, d->f_data);
515 for (i = 0; i < n; ++i)
516 grad[i] *= d->ub[i] - d->lb[i];
520 nlopt_result cdirect(int n, nlopt_func f, void *f_data,
521 const double *lb, const double *ub,
524 nlopt_stopping *stop,
525 double magic_eps, int which_alg)
529 const double *xtol_abs_save;
532 d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
533 d.x = (double *) malloc(sizeof(double) * n*4);
534 if (!d.x) return NLOPT_OUT_OF_MEMORY;
536 for (i = 0; i < n; ++i) {
537 x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
540 d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
542 xtol_abs_save = stop->xtol_abs;
543 stop->xtol_abs = d.x + 3*n;
544 ret = cdirect_unscaled(n, uf, &d, d.x+n, d.x+2*n, x, fmin, stop,
545 magic_eps, which_alg);
546 stop->xtol_abs = xtol_abs_save;
547 for (i = 0; i < n; ++i)
548 x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);