chiark / gitweb /
added qsort_r replacement
[nlopt.git] / cdirect / cdirect.c
1 /* Copyright (c) 2007-2008 Massachusetts Institute of Technology
2  *
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:
10  * 
11  * The above copyright notice and this permission notice shall be
12  * included in all copies or substantial portions of the Software.
13  * 
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. 
21  */
22
23 #include <math.h>
24 #include <stdlib.h>
25 #include <string.h>
26
27 #include "nlopt-util.h"
28 #include "nlopt.h"
29 #include "cdirect.h"
30 #include "redblack.h"
31 #include "config.h"
32
33 #define MIN(a,b) ((a) < (b) ? (a) : (b))
34 #define MAX(a,b) ((a) > (b) ? (a) : (b))
35
36 /***************************************************************************/
37 /* basic data structure:
38  *
39  * a hyper-rectangle is stored as an array of length L = 2n+3, where [1]
40  * is the value (f) of the function at the center, [0] is the "size"
41  * measure (d) of the rectangle, [3..n+2] are the coordinates of the
42  * center (c), [n+3..2n+2] are the widths of the sides (w), and [2]
43  * is an "age" measure for tie-breaking purposes.
44  *
45  * we store the hyper-rectangles in a red-black tree, sorted by (d,f)
46  * in lexographic order, to allow us to perform quick convex-hull
47  * calculations (in the future, we might make this data structure
48  * more sophisticated based on the dynamic convex-hull literature).
49  *
50  * n > 0 always, of course.
51  */
52
53 /* parameters of the search algorithm and various information that
54    needs to be passed around */
55 typedef struct {
56      int n; /* dimension */
57      int L; /* size of each rectangle (2n+3) */
58      double magic_eps; /* Jones' epsilon parameter (1e-4 is recommended) */
59      int which_diam; /* which measure of hyper-rectangle diam to use:
60                         0 = Jones, 1 = Gablonsky */
61      int which_div; /* which way to divide rects:
62                        0: orig. Jones (divide all longest sides)
63                        1: Gablonsky (cubes divide all, rects longest)
64                        2: Jones Encyc. Opt.: pick random longest side */
65      int which_opt; /* which rects are considered "potentially optimal"
66                        0: Jones (all pts on cvx hull, even equal pts)
67                        1: Gablonsky DIRECT-L (pick one pt, if equal pts)
68                        2: ~ 1, but pick points randomly if equal pts 
69                     ... 2 seems to suck compared to just picking oldest pt */
70   
71      const double *lb, *ub;
72      nlopt_stopping *stop; /* stopping criteria */
73      nlopt_func f; void *f_data;
74      double *work; /* workspace, of length >= 2*n */
75      int *iwork; /* workspace, length >= n */
76      double minf, *xmin; /* minimum so far */
77      
78      /* red-black tree of hyperrects, sorted by (d,f,age) in
79         lexographical order */
80      rb_tree rtree;
81      int age; /* age for next new rect */
82      double **hull; /* array to store convex hull */
83      int hull_len; /* allocated length of hull array */
84 } params;
85
86 /***************************************************************************/
87
88 /* Evaluate the "diameter" (d) of a rectangle of widths w[n] 
89
90    We round the result to single precision, which should be plenty for
91    the use we put the diameter to (rect sorting), to allow our
92    performance hack in convex_hull to work (in the Jones and Gablonsky
93    DIRECT algorithms, all of the rects fall into a few diameter
94    values, and we don't want rounding error to spoil this) */
95 static double rect_diameter(int n, const double *w, const params *p)
96 {
97      int i;
98      if (p->which_diam == 0) { /* Jones measure */
99           double sum = 0;
100           for (i = 0; i < n; ++i)
101                sum += w[i] * w[i];
102           /* distance from center to a vertex */
103           return ((float) (sqrt(sum) * 0.5)); 
104      }
105      else { /* Gablonsky measure */
106           double maxw = 0;
107           for (i = 0; i < n; ++i)
108                if (w[i] > maxw)
109                     maxw = w[i];
110           /* half-width of longest side */
111           return ((float) (maxw * 0.5));
112      }
113 }
114
115 #define ALLOC_RECT(rect, L) if (!(rect = (double*) malloc(sizeof(double)*(L)))) return NLOPT_OUT_OF_MEMORY
116
117 static int sort_fv_compare(void *fv_, const void *a_, const void *b_)
118 {
119      const double *fv = (const double *) fv_;
120      int a = *((const int *) a_), b = *((const int *) b_);
121      double fa = MIN(fv[2*a], fv[2*a+1]);
122      double fb = MIN(fv[2*b], fv[2*b+1]);
123      if (fa < fb)
124           return -1;
125      else if (fa > fb)
126           return +1;
127      else
128           return 0;
129 }
130 static void sort_fv(int n, double *fv, int *isort)
131 {
132      int i;
133      for (i = 0; i < n; ++i) isort[i] = i;
134      nlopt_qsort_r(isort, (unsigned) n, sizeof(int), fv, sort_fv_compare);
135 }
136
137 static double function_eval(const double *x, params *p) {
138      double f = p->f(p->n, x, NULL, p->f_data);
139      if (f < p->minf) {
140           p->minf = f;
141           memcpy(p->xmin, x, sizeof(double) * p->n);
142      }
143      p->stop->nevals++;
144      return f;
145 }
146 #define FUNCTION_EVAL(fv,x,p,freeonerr) fv = function_eval(x, p); if (p->minf < p->stop->minf_max) { free(freeonerr); return NLOPT_MINF_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; }
147
148 #define THIRD (0.3333333333333333333333)
149
150 #define EQUAL_SIDE_TOL 5e-2 /* tolerance to equate side sizes */
151
152 /* divide rectangle idiv in the list p->rects */
153 static nlopt_result divide_rect(double *rdiv, params *p)
154 {
155      int i;
156      const const int n = p->n;
157      const int L = p->L;
158      double *c = rdiv + 3; /* center of rect to divide */
159      double *w = c + n; /* widths of rect to divide */
160      double wmax = w[0];
161      int imax = 0, nlongest = 0;
162      rb_node *node;
163
164      for (i = 1; i < n; ++i)
165           if (w[i] > wmax)
166                wmax = w[imax = i];
167      for (i = 0; i < n; ++i)
168           if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL)
169                ++nlongest;
170      if (p->which_div == 1 || (p->which_div == 0 && nlongest == n)) {
171           /* trisect all longest sides, in increasing order of the average
172              function value along that direction */
173           double *fv = p->work;
174           int *isort = p->iwork;
175           for (i = 0; i < n; ++i) {
176                if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL) {
177                     double csave = c[i];
178                     c[i] = csave - w[i] * THIRD;
179                     FUNCTION_EVAL(fv[2*i], c, p, 0);
180                     c[i] = csave + w[i] * THIRD;
181                     FUNCTION_EVAL(fv[2*i+1], c, p, 0);
182                     c[i] = csave;
183                }
184                else {
185                     fv[2*i] = fv[2*i+1] = HUGE_VAL;
186                }
187           }
188           sort_fv(n, fv, isort);
189           if (!(node = rb_tree_find(&p->rtree, rdiv)))
190                return NLOPT_FAILURE;
191           for (i = 0; i < nlongest; ++i) {
192                int k;
193                w[isort[i]] *= THIRD;
194                rdiv[0] = rect_diameter(n, w, p);
195                rdiv[2] = p->age++;
196                node = rb_tree_resort(&p->rtree, node);
197                for (k = 0; k <= 1; ++k) {
198                     double *rnew;
199                     ALLOC_RECT(rnew, L);
200                     memcpy(rnew, rdiv, sizeof(double) * L);
201                     rnew[3 + isort[i]] += w[isort[i]] * (2*k-1);
202                     rnew[1] = fv[2*isort[i]+k];
203                     rnew[2] = p->age++;
204                     if (!rb_tree_insert(&p->rtree, rnew)) {
205                          free(rnew);
206                          return NLOPT_OUT_OF_MEMORY;
207                     }
208                }
209           }
210      }
211      else {
212           int k;
213           if (nlongest > 1 && p->which_div == 2) { 
214                /* randomly choose longest side */
215                i = nlopt_iurand(nlongest);
216                for (k = 0; k < n; ++k)
217                     if (wmax - w[k] <= wmax * EQUAL_SIDE_TOL) {
218                          if (!i) { i = k; break; }
219                          --i;
220                     }
221           }
222           else
223                i = imax; /* trisect longest side */
224           if (!(node = rb_tree_find(&p->rtree, rdiv)))
225                return NLOPT_FAILURE;
226           w[i] *= THIRD;
227           rdiv[0] = rect_diameter(n, w, p);
228           rdiv[2] = p->age++;
229           node = rb_tree_resort(&p->rtree, node);
230           for (k = 0; k <= 1; ++k) {
231                double *rnew;
232                ALLOC_RECT(rnew, L);
233                memcpy(rnew, rdiv, sizeof(double) * L);
234                rnew[3 + i] += w[i] * (2*k-1);
235                FUNCTION_EVAL(rnew[1], rnew + 3, p, rnew);
236                rnew[2] = p->age++;
237                if (!rb_tree_insert(&p->rtree, rnew)) {
238                     free(rnew);
239                     return NLOPT_OUT_OF_MEMORY;
240                }
241           }
242      }
243      return NLOPT_SUCCESS;
244 }
245
246 /***************************************************************************/
247 /* Convex hull algorithm, used later to find the potentially optimal
248    points.  What we really have in DIRECT is a "dynamic convex hull"
249    problem, since we are dynamically adding/removing points and
250    updating the hull, but I haven't implemented any of the fancy
251    algorithms for this problem yet. */
252
253 /* Find the lower convex hull of a set of points (x,y) stored in a rb-tree
254    of pointers to {x,y} arrays sorted in lexographic order by (x,y).
255
256    Unlike standard convex hulls, we allow redundant points on the hull,
257    and even allow duplicate points if allow_dups is nonzero.
258
259    The return value is the number of points in the hull, with pointers
260    stored in hull[i] (should be an array of length >= t->N).
261 */
262 static int convex_hull(rb_tree *t, double **hull, int allow_dups)
263 {
264      int nhull = 0;
265      double minslope;
266      double xmin, xmax, yminmin, ymaxmin;
267      rb_node *n, *nmax;
268
269      /* Monotone chain algorithm [Andrew, 1979]. */
270
271      n = rb_tree_min(t);
272      if (!n) return 0;
273      nmax = rb_tree_max(t);
274
275      xmin = n->k[0];
276      yminmin = n->k[1];
277      xmax = nmax->k[0];
278
279      if (allow_dups)
280           do { /* include any duplicate points at (xmin,yminmin) */
281                hull[nhull++] = n->k;
282                n = rb_tree_succ(n);
283           } while (n && n->k[0] == xmin && n->k[1] == yminmin);
284      else
285           hull[nhull++] = n->k;
286
287      if (xmin == xmax) return nhull;
288
289      /* set nmax = min mode with x == xmax */
290 #if 0
291      while (nmax->k[0] == xmax)
292           nmax = rb_tree_pred(nmax); /* non-NULL since xmin != xmax */
293      nmax = rb_tree_succ(nmax);
294 #else
295      /* performance hack (see also below) */
296      {
297           double kshift[2];
298           kshift[0] = xmax * (1 - 1e-13);
299           kshift[1] = -HUGE_VAL;
300           nmax = rb_tree_find_gt(t, kshift); /* non-NULL since xmin != xmax */
301      }
302 #endif
303
304      ymaxmin = nmax->k[1];
305      minslope = (ymaxmin - yminmin) / (xmax - xmin);
306
307      /* set n = first node with x != xmin */
308 #if 0
309      while (n->k[0] == xmin)
310           n = rb_tree_succ(n); /* non-NULL since xmin != xmax */
311 #else
312      /* performance hack (see also below) */
313      {
314           double kshift[2];
315           kshift[0] = xmin * (1 + 1e-13);
316           kshift[1] = -HUGE_VAL;
317           n = rb_tree_find_gt(t, kshift); /* non-NULL since xmin != xmax */
318      }
319 #endif
320
321      for (; n != nmax; n = rb_tree_succ(n)) { 
322           double *k = n->k;
323           if (k[1] > yminmin + (k[0] - xmin) * minslope)
324                continue;
325
326           /* performance hack: most of the points in DIRECT lie along
327              vertical lines at a few x values, and we can exploit this */
328           if (nhull && k[0] == hull[nhull - 1][0]) { /* x == previous x */
329                if (k[1] > hull[nhull - 1][1]) {
330                     double kshift[2];
331                     /* because of the round to float in rect_diameter, above,
332                        it shouldn't be possible for two diameters (x values)
333                        to have a fractional difference < 1e-13.  Note
334                        that k[0] > 0 always in DIRECT */
335                     kshift[0] = k[0] * (1 + 1e-13);
336                     kshift[1] = -HUGE_VAL;
337                     n = rb_tree_pred(rb_tree_find_gt(t, kshift));
338                     continue;
339                }
340                else { /* equal y values, add to hull */
341                     if (allow_dups)
342                          hull[nhull++] = k;
343                     continue;
344                }
345           }
346
347           /* remove points until we are making a "left turn" to k */
348           while (nhull > 1) {
349                double *t1 = hull[nhull - 1], *t2;
350
351                /* because we allow equal points in our hull, we have
352                   to modify the standard convex-hull algorithm slightly:
353                   we need to look backwards in the hull list until we
354                   find a point t2 != t1 */
355                int it2 = nhull - 2;
356                do {
357                     t2 = hull[it2--];
358                } while (it2 >= 0 && t2[0] == t1[0] && t2[1] == t1[1]);
359                if (it2 < 0) break;
360
361                /* cross product (t1-t2) x (k-t2) > 0 for a left turn: */
362                if ((t1[0]-t2[0]) * (k[1]-t2[1])
363                    - (t1[1]-t2[1]) * (k[0]-t2[0]) >= 0)
364                     break;
365                --nhull;
366           }
367           hull[nhull++] = k;
368      }
369
370      if (allow_dups)
371           do { /* include any duplicate points at (xmax,ymaxmin) */
372                hull[nhull++] = nmax->k;
373                nmax = rb_tree_succ(nmax);
374           } while (nmax && nmax->k[0] == xmax && n->k[1] == ymaxmin);
375      else
376           hull[nhull++] = nmax->k;
377
378      return nhull;
379 }
380
381 /***************************************************************************/
382
383 static int small(double *w, params *p)
384 {
385      int i;
386      for (i = 0; i < p->n; ++i)
387           if (w[i] > p->stop->xtol_abs[i] &&
388               w[i] > (p->ub[i] - p->lb[i]) * p->stop->xtol_rel)
389                return 0;
390      return 1;
391 }
392
393 static nlopt_result divide_good_rects(params *p)
394 {
395      const int n = p->n;
396      double **hull;
397      int nhull, i, xtol_reached = 1, divided_some = 0;
398      double magic_eps = p->magic_eps;
399
400      if (p->hull_len < p->rtree.N) {
401           p->hull_len += p->rtree.N;
402           p->hull = (double **) realloc(p->hull, sizeof(double*)*p->hull_len);
403           if (!p->hull) return NLOPT_OUT_OF_MEMORY;
404      }
405      nhull = convex_hull(&p->rtree, hull = p->hull, p->which_opt != 1);
406  divisions:
407      for (i = 0; i < nhull; ++i) {
408           double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
409           int im, ip;
410
411           /* find unequal points before (im) and after (ip) to get slope */
412           for (im = i-1; im >= 0 && hull[im][0] == hull[i][0]; --im);
413           for (ip = i+1; ip < nhull && hull[ip][0] == hull[i][0]; ++ip);
414
415           if (im >= 0)
416                K1 = (hull[i][1] - hull[im][1]) / (hull[i][0] - hull[im][0]);
417           if (ip < nhull)
418                K1 = (hull[i][1] - hull[ip][1]) / (hull[i][0] - hull[ip][0]);
419           K = MAX(K1, K2);
420           if (hull[i][1] - K * hull[i][0]
421               <= p->minf - magic_eps * fabs(p->minf) || ip == nhull) {
422                /* "potentially optimal" rectangle, so subdivide */
423                nlopt_result ret = divide_rect(hull[i], p);
424                divided_some = 1;
425                if (ret != NLOPT_SUCCESS) return ret;
426                xtol_reached = xtol_reached && small(hull[i] + 3+n, p);
427           }
428
429           /* for the DIRECT-L variant, we only divide one rectangle out
430              of all points with equal diameter and function values
431              ... note that for p->which_opt == 1, i == ip-1 should be a no-op
432                  anyway, since we set allow_dups=0 in convex_hull above */
433           if (p->which_opt == 1)
434                i = ip - 1; /* skip to next unequal point for next iteration */
435           else if (p->which_opt == 2) /* like DIRECT-L but randomized */
436                i += nlopt_iurand(ip - i); /* possibly do another equal pt */
437      }
438      if (!divided_some) {
439           if (magic_eps != 0) {
440                magic_eps = 0;
441                goto divisions; /* try again */
442           }
443           else { /* WTF? divide largest rectangle with smallest f */
444                /* (note that this code actually gets called from time
445                   to time, and the heuristic here seems to work well,
446                   but I don't recall this situation being discussed in
447                   the references?) */
448                rb_node *max = rb_tree_max(&p->rtree);
449                rb_node *pred = max;
450                double wmax = max->k[0];
451                do { /* note: this loop is O(N) worst-case time */
452                     max = pred;
453                     pred = rb_tree_pred(max);
454                } while (pred && pred->k[0] == wmax);
455                return divide_rect(max->k, p);
456           }
457      }
458      return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
459 }
460
461 /***************************************************************************/
462
463 /* lexographic sort order (d,f,age) of hyper-rects, for red-black tree */
464 int cdirect_hyperrect_compare(double *a, double *b)
465 {
466      if (a[0] < b[0]) return -1;
467      if (a[0] > b[0]) return +1;
468      if (a[1] < b[1]) return -1;
469      if (a[1] > b[1]) return +1;
470      if (a[2] < b[2]) return -1;
471      if (a[2] > b[2]) return +1;
472      return (int) (a - b); /* tie-breaker, shouldn't be needed */
473 }
474
475 /***************************************************************************/
476
477 nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
478                               const double *lb, const double *ub,
479                               double *x,
480                               double *minf,
481                               nlopt_stopping *stop,
482                               double magic_eps, int which_alg)
483 {
484      params p;
485      int i;
486      double *rnew;
487      nlopt_result ret = NLOPT_OUT_OF_MEMORY;
488
489      p.magic_eps = magic_eps;
490      p.which_diam = which_alg % 3;
491      p.which_div = (which_alg / 3) % 3;
492      p.which_opt = (which_alg / (3*3)) % 3;
493      p.lb = lb; p.ub = ub;
494      p.stop = stop;
495      p.n = n;
496      p.L = 2*n+3;
497      p.f = f;
498      p.f_data = f_data;
499      p.xmin = x;
500      p.minf = HUGE_VAL;
501      p.work = 0;
502      p.iwork = 0;
503      p.hull = 0;
504      p.age = 0;
505
506      rb_tree_init(&p.rtree, cdirect_hyperrect_compare);
507
508      p.work = (double *) malloc(sizeof(double) * (2*n));
509      if (!p.work) goto done;
510      p.iwork = (int *) malloc(sizeof(int) * n);
511      if (!p.iwork) goto done;
512      p.hull_len = 128; /* start with a reasonable number */
513      p.hull = (double **) malloc(sizeof(double *) * p.hull_len);
514      if (!p.hull) goto done;
515
516      if (!(rnew = (double *) malloc(sizeof(double) * p.L))) goto done;
517      for (i = 0; i < n; ++i) {
518           rnew[3+i] = 0.5 * (lb[i] + ub[i]);
519           rnew[3+n+i] = ub[i] - lb[i];
520      }
521      rnew[0] = rect_diameter(n, rnew+3+n, &p);
522      rnew[1] = function_eval(rnew+3, &p);
523      rnew[2] = p.age++;
524      if (!rb_tree_insert(&p.rtree, rnew)) {
525           free(rnew);
526           goto done;
527      }
528
529      ret = divide_rect(rnew, &p);
530      if (ret != NLOPT_SUCCESS) goto done;
531
532      while (1) {
533           double minf0 = p.minf;
534           ret = divide_good_rects(&p);
535           if (ret != NLOPT_SUCCESS) goto done;
536           if (p.minf < minf0 && nlopt_stop_f(p.stop, p.minf, minf0)) {
537                ret = NLOPT_FTOL_REACHED;
538                goto done;
539           }
540      }
541
542  done:
543      rb_tree_destroy_with_keys(&p.rtree);
544      free(p.hull);
545      free(p.iwork);
546      free(p.work);
547               
548      *minf = p.minf;
549      return ret;
550 }
551
552 /* in the conventional DIRECT-type algorithm, we first rescale our
553    coordinates to a unit hypercube ... we do this simply by
554    wrapping cdirect() around cdirect_unscaled(). */
555
556 double cdirect_uf(int n, const double *xu, double *grad, void *d_)
557 {
558      cdirect_uf_data *d = (cdirect_uf_data *) d_;
559      double f;
560      int i;
561      for (i = 0; i < n; ++i)
562           d->x[i] = d->lb[i] + xu[i] * (d->ub[i] - d->lb[i]);
563      f = d->f(n, d->x, grad, d->f_data);
564      if (grad)
565           for (i = 0; i < n; ++i)
566                grad[i] *= d->ub[i] - d->lb[i];
567      return f;
568 }
569
570 nlopt_result cdirect(int n, nlopt_func f, void *f_data,
571                      const double *lb, const double *ub,
572                      double *x,
573                      double *minf,
574                      nlopt_stopping *stop,
575                      double magic_eps, int which_alg)
576 {
577      cdirect_uf_data d;
578      nlopt_result ret;
579      const double *xtol_abs_save;
580      int i;
581
582      d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
583      d.x = (double *) malloc(sizeof(double) * n*4);
584      if (!d.x) return NLOPT_OUT_OF_MEMORY;
585      
586      for (i = 0; i < n; ++i) {
587           x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
588           d.x[n+i] = 0;
589           d.x[2*n+i] = 1;
590           d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
591      }
592      xtol_abs_save = stop->xtol_abs;
593      stop->xtol_abs = d.x + 3*n;
594      ret = cdirect_unscaled(n, cdirect_uf, &d, d.x+n, d.x+2*n, x, minf, stop,
595                             magic_eps, which_alg);
596      stop->xtol_abs = xtol_abs_save;
597      for (i = 0; i < n; ++i)
598           x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);
599      free(d.x);
600      return ret;
601 }