5 #include "nlopt-util.h"
10 #define MIN(a,b) ((a) < (b) ? (a) : (b))
11 #define MAX(a,b) ((a) > (b) ? (a) : (b))
13 /***************************************************************************/
14 /* basic data structure:
16 * a hyper-rectangle is stored as an array of length 2n+2, where [0]
17 * is the value of the function at the center, [1] is the "size"
18 * measure of the rectangle, [2..n+1] are the coordinates of the
19 * center, and [n+2..2n+1] are the widths of the sides.
21 * a list of rectangles is just an array of N hyper-rectangles
22 * stored as an N x (2n+1) in row-major order. Generally,
23 * we allocate more than we need, allocating Na hyper-rectangles.
28 #define RECT_LEN(n) (2*(n)+2) /* number of double values in a hyperrect */
30 /* parameters of the search algorithm and various information that
31 needs to be passed around */
33 double magic_eps; /* Jones' epsilon parameter (1e-4 is recommended) */
34 int which_diam; /* which measure of hyper-rectangle diam to use:
35 0 = Jones, 1 = Gablonsky */
36 const double *lb, *ub;
37 nlopt_stopping *stop; /* stopping criteria */
39 nlopt_func f; void *f_data;
40 double *work; /* workspace, of length >= 2*n */
41 int *iwork, iwork_len; /* workspace, of length iwork_len >= n */
42 double fmin, *xmin; /* minimum so far */
45 /***************************************************************************/
47 /* evaluate the "diameter" (d) of a rectangle of widths w[n] */
48 static double rect_diameter(int n, const double *w, const params *p)
51 if (p->which_diam == 0) { /* Jones measure */
53 for (i = 0; i < n; ++i)
55 return sqrt(sum) * 0.5; /* distance from center to a vertex */
57 else { /* Gablonsky measure */
59 for (i = 0; i < n; ++i)
62 return maxw * 0.5; /* half-width of longest side */
66 #define CUBE_TOL 5e-2 /* fractional tolerance to call something a "cube" */
68 /* return true if the elements of w[n] (> 0) are all equal to within a
69 fractional tolerance of tol (i.e. they are the widths of a hypercube) */
70 static int all_equal(int n, const double *w, double tol)
75 for (i = 1; i < n; ++i) {
76 if (w[i] < wmin) wmin = w[i];
77 if (w[i] > wmax) wmax = w[i];
79 return (wmax - wmin) < tol * wmax;
82 static double *alloc_rects(int n, int *Na, double *rects, int newN)
88 return realloc(rects, sizeof(double) * RECT_LEN(n) * (*Na));
91 #define ALLOC_RECTS(n, Nap, rects, newN) if (!(rects = alloc_rects(n, Nap, rects, newN))) return NLOPT_OUT_OF_MEMORY
93 static double *fv_qsort = 0;
94 static int sort_fv_compare(const void *a_, const void *b_)
96 int a = *((const int *) a_), b = *((const int *) b_);
97 double fa = MIN(fv_qsort[2*a], fv_qsort[2*a+1]);
98 double fb = MIN(fv_qsort[2*b], fv_qsort[2*b+1]);
106 static void sort_fv(int n, double *fv, int *isort)
109 for (i = 0; i < n; ++i) isort[i] = i;
110 fv_qsort = fv; /* not re-entrant, sigh... */
111 qsort(isort, (unsigned) n, sizeof(int), sort_fv_compare);
115 static double function_eval(const double *x, params *p) {
116 double f = p->f(p->n, x, NULL, p->f_data);
119 memcpy(p->xmin, x, sizeof(double) * p->n);
124 #define FUNCTION_EVAL(fv,x,p) fv = function_eval(x, p); if (p->fmin < p->stop->fmin_max) return NLOPT_FMIN_MAX_REACHED; else if (nlopt_stop_evals((p)->stop)) return NLOPT_MAXEVAL_REACHED; else if (nlopt_stop_time((p)->stop)) return NLOPT_MAXTIME_REACHED
126 #define THIRD (0.3333333333333333333333)
129 /* divide rectangle idiv in the list rects */
130 static nlopt_result divide_rect(int *N, int *Na, double **rects, int idiv,
134 const const int n = p->n;
135 const int L = RECT_LEN(n);
137 double *c = r + L*idiv + 2; /* center of rect to divide */
138 double *w = c + n; /* widths of rect to divide */
140 if (all_equal(n, w, CUBE_TOL)) { /* divide all dimensions */
141 double *fv = p->work;
142 int *isort = p->iwork;
143 for (i = 0; i < n; ++i) {
145 c[i] = csave - w[i] * THIRD;
146 FUNCTION_EVAL(fv[2*i], c, p);
147 c[i] = csave + w[i] * THIRD;
148 FUNCTION_EVAL(fv[2*i+1], c, p);
151 sort_fv(n, fv, isort);
152 ALLOC_RECTS(n, Na, r, (*N)+2*n);
153 *rects = r; c = r + L*idiv + 2; w = c + n;
154 for (i = 0; i < n; ++i) {
156 w[isort[i]] *= THIRD;
157 r[L*idiv + 1] = rect_diameter(n, w, p);
158 for (k = 0; k <= 1; ++k) {
159 memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
160 r[L*(*N) + 2 + isort[i]] += w[isort[i]] * (2*k-1);
161 r[L*(*N)] = fv[2*isort[i]+k];
166 else { /* divide longest side by 3 and split off 2 new rectangles */
169 for (i = 1; i < n; ++i)
172 ALLOC_RECTS(n, Na, r, (*N)+2);
173 *rects = r; c = r + L*idiv + 2; w = c + n;
175 r[L*idiv + 1] = rect_diameter(n, w, p);
176 for (i = 0; i <= 1; ++i) {
177 memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
178 r[L*(*N) + 2 + imax] += w[imax] * (2*i-1); /* move center */
180 FUNCTION_EVAL(r[L*((*N)-1)], r + L*((*N)-1) + 2, p);
183 return NLOPT_SUCCESS;
186 /***************************************************************************/
187 /* O(N log N) convex hull algorithm, used later to find the potentially
190 /* sort ihull by xy in lexographic order by x,y */
191 static int s_qsort = 1; static double *xy_qsort = 0;
192 static int sort_xy_compare(const void *a_, const void *b_)
194 int a = *((const int *) a_), b = *((const int *) b_);
195 double xa = xy_qsort[a*s_qsort+1], xb = xy_qsort[b*s_qsort+1];
196 if (xa < xb) return -1;
197 else if (xb < xa) return +1;
199 double ya = xy_qsort[a*s_qsort], yb = xy_qsort[b*s_qsort];
200 if (ya < yb) return -1;
201 else if (ya > yb) return +1;
205 static void sort_xy(int N, double *xy, int s, int *isort)
209 for (i = 0; i < N; ++i) isort[i] = i;
210 s_qsort = s; xy_qsort = xy;
211 qsort(isort, (unsigned) N, sizeof(int), sort_xy_compare);
215 /* Find the lower convex hull of a set of points (xy[s*i+1], xy[s*i]), where
216 0 <= i < N and s >= 2.
218 The return value is the number of points in the hull, with indices
219 stored in ihull. ihull and is should point to arrays of length >= N.
221 Note that we don't allow redundant points along the same line in the
222 hull, similar to Gablonsky's version of DIRECT and differing from
224 static int convex_hull(int N, double *xy, int s, int *ihull, int *is)
226 int minmin; /* first index (smallest y) with min x */
227 int minmax; /* last index (largest y) with min x */
228 int maxmin; /* first index (smallest y) with max x */
229 int maxmax; /* last index (largest y) with max x */
234 /* Monotone chain algorithm [Andrew, 1979]. */
236 sort_xy(N, xy, s, is);
238 xmin = xy[s*is[minmin=0]+1]; xmax = xy[s*is[maxmax=N-1]+1];
240 if (xmin == xmax) { /* degenerate case */
241 ihull[nhull++] = is[minmin];
245 for (minmax = minmin; minmax+1 < N && xy[s*is[minmax+1]+1]==xmin;
247 for (maxmin = maxmax; maxmin-1>=0 && xy[s*is[maxmin-1]+1]==xmax;
250 minslope = (xy[s*is[maxmin]] - xy[s*is[minmin]]) / (xmax - xmin);
252 ihull[nhull++] = is[minmin];
253 for (i = minmax + 1; i < maxmin; ++i) {
255 if (xy[s*k] > xy[s*is[minmin]] + (xy[s*k+1] - xmin) * minslope)
257 /* remove points until we are making a "left turn" to k */
259 int t1 = ihull[nhull - 1], t2 = ihull[nhull - 2];
260 /* cross product (t1-t2) x (k-t2) > 0 for a left turn: */
261 if ((xy[s*t1+1]-xy[s*t2+1]) * (xy[s*k]-xy[s*t2])
262 - (xy[s*t1]-xy[s*t2]) * (xy[s*k+1]-xy[s*t2+1]) > 0)
268 ihull[nhull++] = is[maxmin];
272 /***************************************************************************/
274 static int small(double *w, params *p)
277 for (i = 0; i < p->n; ++i)
278 if (w[i] > p->stop->xtol_abs[i] &&
279 w[i] > (p->ub[i] - p->lb[i]) * p->stop->xtol_rel)
284 static nlopt_result divide_good_rects(int *N, int *Na, double **rects,
288 const int L = RECT_LEN(n);
289 int *ihull, nhull, i, xtol_reached = 1, divided_some = 0;
291 double magic_eps = p->magic_eps;
293 if (p->iwork_len < n + 2*(*N)) {
294 p->iwork_len = p->iwork_len + n + 2*(*N);
295 p->iwork = (int *) realloc(p->iwork, sizeof(int) * p->iwork_len);
297 return NLOPT_OUT_OF_MEMORY;
300 nhull = convex_hull(*N, r, L, ihull, ihull + *N);
302 for (i = 0; i < nhull; ++i) {
303 double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
305 K1 = (r[L*ihull[i]] - r[L*ihull[i-1]]) /
306 (r[L*ihull[i]+1] - r[L*ihull[i-1]+1]);
308 K1 = (r[L*ihull[i]] - r[L*ihull[i+1]]) /
309 (r[L*ihull[i]+1] - r[L*ihull[i+1]+1]);
311 if (r[L*ihull[i]] - K * r[L*ihull[i]+1]
312 <= p->fmin - magic_eps * fabs(p->fmin)) {
313 /* "potentially optimal" rectangle, so subdivide */
316 ret = divide_rect(N, Na, rects, ihull[i], p);
317 r = *rects; /* may have grown */
318 if (ret != NLOPT_SUCCESS) return ret;
319 xtol_reached = xtol_reached && small(r + L*ihull[i] + 2+n, p);
323 if (magic_eps != 0) {
325 goto divisions; /* try again */
327 else { /* WTF? divide largest rectangle */
330 for (i = 1; i < *N; ++i)
332 wmax = r[L*(imax=i)+1];
333 return divide_rect(N, Na, rects, imax, p);
336 return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
339 /***************************************************************************/
341 nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
342 const double *lb, const double *ub,
345 nlopt_stopping *stop,
346 double magic_eps, int which_alg)
350 int Na = 100, N = 1, i, x_center = 1;
351 nlopt_result ret = NLOPT_OUT_OF_MEMORY;
353 p.magic_eps = magic_eps;
354 p.which_diam = which_alg & 1;
355 p.lb = lb; p.ub = ub;
361 p.fmin = f(n, x, NULL, f_data); stop->nevals++;
365 p.work = (double *) malloc(sizeof(double) * 2*n);
366 if (!p.work) goto done;
367 rects = (double *) malloc(sizeof(double) * Na * RECT_LEN(n));
368 if (!rects) goto done;
369 p.iwork = (int *) malloc(sizeof(int) * (p.iwork_len = 2*Na + n));
370 if (!p.iwork) goto done;
372 for (i = 0; i < n; ++i) {
373 rects[2+i] = 0.5 * (lb[i] + ub[i]);
375 && (fabs(rects[2+i]-x[i]) < 1e-13*(1+fabs(x[i])));
376 rects[2+n+i] = ub[i] - lb[i];
378 rects[1] = rect_diameter(n, rects+2+n, &p);
380 rects[0] = p.fmin; /* avoid computing f(center) twice */
382 rects[0] = function_eval(rects+2, &p);
384 ret = divide_rect(&N, &Na, &rects, 0, &p);
385 if (ret != NLOPT_SUCCESS) goto done;
388 double fmin0 = p.fmin;
389 ret = divide_good_rects(&N, &Na, &rects, &p);
390 if (ret != NLOPT_SUCCESS) goto done;
391 if (nlopt_stop_f(p.stop, p.fmin, fmin0)) {
392 ret = NLOPT_FTOL_REACHED;
406 /* in the conventional DIRECT-type algorithm, we first rescale our
407 coordinates to a unit hypercube ... we do this simply by
408 wrapping cdirect() around cdirect_unscaled(). */
414 const double *lb, *ub;
416 static double uf(int n, const double *xu, double *grad, void *d_)
418 uf_data *d = (uf_data *) d_;
421 for (i = 0; i < n; ++i)
422 d->x[i] = d->lb[i] + xu[i] * (d->ub[i] - d->lb[i]);
423 f = d->f(n, d->x, grad, d->f_data);
425 for (i = 0; i < n; ++i)
426 grad[i] *= d->ub[i] - d->lb[i];
430 nlopt_result cdirect(int n, nlopt_func f, void *f_data,
431 const double *lb, const double *ub,
434 nlopt_stopping *stop,
435 double magic_eps, int which_alg)
439 const double *xtol_abs_save;
442 d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
443 d.x = (double *) malloc(sizeof(double) * n*4);
444 if (!d.x) return NLOPT_OUT_OF_MEMORY;
446 for (i = 0; i < n; ++i) {
447 x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
450 d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
452 xtol_abs_save = stop->xtol_abs;
453 stop->xtol_abs = d.x + 3*n;
454 ret = cdirect_unscaled(n, uf, &d, d.x+n, d.x+2*n, x, fmin, stop,
455 magic_eps, which_alg);
456 stop->xtol_abs = xtol_abs_save;
457 for (i = 0; i < n; ++i)
458 x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);