cdirect/cdirect.c

   1 #include <math.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4
   5 #include "nlopt-util.h"
   6 #include "nlopt.h"
   7 #include "cdirect.h"
   8 #include "redblack.h"
   9 #include "config.h"
  10
  11 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  12 #define MAX(a,b) ((a) > (b) ? (a) : (b))
  13
  14 /***************************************************************************/
  15 /* basic data structure:
  16  *
  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).
  21  *
  22  * a list of rectangles is just an array of N hyper-rectangles
  23  * stored as an N x L in row-major order.  Generally,
  24  * we allocate more than we need, allocating Na hyper-rectangles.
  25  *
  26  * n > 0 always
  27  */
  28
  29 #define RECT_LEN(n) (2*(n)+2) /* number of double values in a hyperrect */
  30
  31 /* parameters of the search algorithm and various information that
  32    needs to be passed around */
  33 typedef struct {
  34      int n; /* dimension */
  35      int L; /* RECT_LEN(n) */
  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: Gablonsky (cubes divide all, rects longest)
  41                        1: orig. Jones (divide all longest sides)
  42                        2: Jones Encyc. Opt.: pick random longest side */
  43      const double *lb, *ub;
  44      nlopt_stopping *stop; /* stopping criteria */
  45      nlopt_func f; void *f_data;
  46      double *work; /* workspace, of length >= 2*n */
  47      int *iwork; /* workspace, length >= n */
  48      double fmin, *xmin; /* minimum so far */
  49
  50      /* red-black tree of hyperrects, sorted by (d,f) in
  51         lexographical order */
  52      rb_tree rtree;
  53      double **hull; /* array to store convex hull */
  54      int hull_len; /* allocated length of hull array */
  55 } params;
  56
  57 /***************************************************************************/
  58
  59 /* evaluate the "diameter" (d) of a rectangle of widths w[n] */
  60 static double rect_diameter(int n, const double *w, const params *p)
  61 {
  62      int i;
  63      if (p->which_diam == 0) { /* Jones measure */
  64           double sum = 0;
  65           for (i = 0; i < n; ++i)
  66                sum += w[i] * w[i];
  67           return sqrt(sum) * 0.5; /* distance from center to a vertex */
  68      }
  69      else { /* Gablonsky measure */
  70           double maxw = 0;
  71           for (i = 0; i < n; ++i)
  72                if (w[i] > maxw)
  73                     maxw = w[i];
  74           return maxw * 0.5; /* half-width of longest side */
  75      }
  76 }
  77
  78 #define ALLOC_RECT(rect, L) if (!(rect = (double*) malloc(sizeof(double)*(L)))) return NLOPT_OUT_OF_MEMORY
  79
  80 static double *fv_qsort = 0;
  81 static int sort_fv_compare(const void *a_, const void *b_)
  82 {
  83      int a = *((const int *) a_), b = *((const int *) b_);
  84      double fa = MIN(fv_qsort[2*a], fv_qsort[2*a+1]);
  85      double fb = MIN(fv_qsort[2*b], fv_qsort[2*b+1]);
  86      if (fa < fb)
  87           return -1;
  88      else if (fa > fb)
  89           return +1;
  90      else
  91           return 0;
  92 }
  93 static void sort_fv(int n, double *fv, int *isort)
  94 {
  95      int i;
  96      for (i = 0; i < n; ++i) isort[i] = i;
  97      fv_qsort = fv; /* not re-entrant, sigh... */
  98      qsort(isort, (unsigned) n, sizeof(int), sort_fv_compare);
  99      fv_qsort = 0;
 100 }
 101
 102 static double function_eval(const double *x, params *p) {
 103      double f = p->f(p->n, x, NULL, p->f_data);
 104      if (f < p->fmin) {
 105           p->fmin = f;
 106           memcpy(p->xmin, x, sizeof(double) * p->n);
 107      }
 108      p->stop->nevals++;
 109      return f;
 110 }
 111 #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; }
 112
 113 #define THIRD (0.3333333333333333333333)
 114
 115 #define EQUAL_SIDE_TOL 5e-2 /* tolerance to equate side sizes */
 116
 117 /* divide rectangle idiv in the list p->rects */
 118 static nlopt_result divide_rect(double *rdiv, params *p)
 119 {
 120      int i;
 121      const const int n = p->n;
 122      const int L = p->L;
 123      double *c = rdiv + 2; /* center of rect to divide */
 124      double *w = c + n; /* widths of rect to divide */
 125      double wmax = w[0];
 126      int imax = 0, nlongest = 0;
 127      rb_node *node;
 128
 129      for (i = 1; i < n; ++i)
 130           if (w[i] > wmax)
 131                wmax = w[imax = i];
 132      for (i = 0; i < n; ++i)
 133           if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL)
 134                ++nlongest;
 135      if (p->which_div == 1 || (p->which_div == 0 && nlongest == n)) {
 136           /* trisect all longest sides, in increasing order of the average
 137              function value along that direction */
 138           double *fv = p->work;
 139           int *isort = p->iwork;
 140           for (i = 0; i < n; ++i) {
 141                if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL) {
 142                     double csave = c[i];
 143                     c[i] = csave - w[i] * THIRD;
 144                     FUNCTION_EVAL(fv[2*i], c, p, 0);
 145                     c[i] = csave + w[i] * THIRD;
 146                     FUNCTION_EVAL(fv[2*i+1], c, p, 0);
 147                     c[i] = csave;
 148                }
 149                else {
 150                     fv[2*i] = fv[2*i+1] = HUGE_VAL;
 151                }
 152           }
 153           sort_fv(n, fv, isort);
 154           if (!(node = rb_tree_find(&p->rtree, rdiv)))
 155                return NLOPT_FAILURE;
 156           for (i = 0; i < nlongest; ++i) {
 157                int k;
 158                w[isort[i]] *= THIRD;
 159                rdiv[0] = rect_diameter(n, w, p);
 160                node = rb_tree_resort(&p->rtree, node);
 161                for (k = 0; k <= 1; ++k) {
 162                     double *rnew;
 163                     ALLOC_RECT(rnew, L);
 164                     memcpy(rnew, rdiv, sizeof(double) * L);
 165                     rnew[2 + isort[i]] += w[isort[i]] * (2*k-1);
 166                     rnew[1] = fv[2*isort[i]+k];
 167                     if (!rb_tree_insert(&p->rtree, rnew)) {
 168                          free(rnew);
 169                          return NLOPT_FAILURE;
 170                     }
 171                }
 172           }
 173      }
 174      else {
 175           int k;
 176           if (nlongest > 1 && p->which_div == 2) {
 177                /* randomly choose longest side */
 178                i = nlopt_iurand(nlongest);
 179                for (k = 0; k < n; ++k)
 180                     if (wmax - w[k] <= wmax * EQUAL_SIDE_TOL) {
 181                          if (!i) { i = k; break; }
 182                          --i;
 183                     }
 184           }
 185           else
 186                i = imax; /* trisect longest side */
 187           if (!(node = rb_tree_find(&p->rtree, rdiv)))
 188                return NLOPT_FAILURE;
 189           w[i] *= THIRD;
 190           rdiv[0] = rect_diameter(n, w, p);
 191           node = rb_tree_resort(&p->rtree, node);
 192           for (k = 0; k <= 1; ++k) {
 193                double *rnew;
 194                ALLOC_RECT(rnew, L);
 195                memcpy(rnew, rdiv, sizeof(double) * L);
 196                rnew[2 + i] += w[i] * (2*k-1);
 197                FUNCTION_EVAL(rnew[1], rnew + 2, p, rnew);
 198                if (!rb_tree_insert(&p->rtree, rnew)) {
 199                     free(rnew);
 200                     return NLOPT_FAILURE;
 201                }
 202           }
 203      }
 204      return NLOPT_SUCCESS;
 205 }
 206
 207 /***************************************************************************/
 208 /* O(N log N) convex hull algorithm, used later to find the potentially
 209    optimal points */
 210
 211 /* Find the lower convex hull of a set of points (x,y) stored in a rb-tree
 212    of pointers to {x,y} arrays sorted in lexographic order by (x,y).
 213
 214    Unlike standard convex hulls, we allow redundant points on the hull.
 215
 216    The return value is the number of points in the hull, with pointers
 217    stored in hull[i] (should be an array of length >= t->N).
 218 */
 219 static int convex_hull(rb_tree *t, double **hull)
 220 {
 221      int nhull = 0;
 222      double minslope;
 223      double xmin, xmax, yminmin, ymaxmin;
 224      rb_node *n, *nmax;
 225
 226      /* Monotone chain algorithm [Andrew, 1979]. */
 227
 228      n = rb_tree_min(t);
 229      if (!n) return 0;
 230      nmax = rb_tree_max(t);
 231
 232      xmin = n->k[0];
 233      yminmin = n->k[1];
 234      xmax = nmax->k[0];
 235
 236      hull[nhull++] = n->k;
 237      if (xmin == xmax) return nhull;
 238
 239      /* set nmax = min mode with x == xmax */
 240      while (nmax->k[0] == xmax)
 241           nmax = rb_tree_pred(nmax); /* non-NULL since xmin != xmax */
 242      nmax = rb_tree_succ(nmax);
 243
 244      ymaxmin = nmax->k[1];
 245      minslope = (ymaxmin - yminmin) / (xmax - xmin);
 246
 247      /* set n = first node with x != xmin */
 248      while (n->k[0] == xmin)
 249           n = rb_tree_succ(n); /* non-NULL since xmin != xmax */
 250
 251      for (; n != nmax; n = rb_tree_succ(n)) {
 252           double *k = n->k;
 253           if (k[1] > yminmin + (k[0] - xmin) * minslope)
 254                continue;
 255           /* remove points until we are making a "left turn" to k */
 256           while (nhull > 1) {
 257                double *t1 = hull[nhull - 1], *t2 = hull[nhull - 2];
 258                /* cross product (t1-t2) x (k-t2) > 0 for a left turn: */
 259                if ((t1[0]-t2[0]) * (k[1]-t2[1])
 260                    - (t1[1]-t2[1]) * (k[0]-t2[0]) >= 0)
 261                     break;
 262                --nhull;
 263           }
 264           hull[nhull++] = k;
 265      }
 266      hull[nhull++] = nmax->k;
 267      return nhull;
 268 }
 269
 270 /***************************************************************************/
 271
 272 static int small(double *w, params *p)
 273 {
 274      int i;
 275      for (i = 0; i < p->n; ++i)
 276           if (w[i] > p->stop->xtol_abs[i] &&
 277               w[i] > (p->ub[i] - p->lb[i]) * p->stop->xtol_rel)
 278                return 0;
 279      return 1;
 280 }
 281
 282 static nlopt_result divide_good_rects(params *p)
 283 {
 284      const int n = p->n;
 285      double **hull;
 286      int nhull, i, xtol_reached = 1, divided_some = 0;
 287      double magic_eps = p->magic_eps;
 288
 289      if (p->hull_len < p->rtree.N) {
 290           p->hull_len += p->rtree.N;
 291           p->hull = (double **) realloc(p->hull, sizeof(double*)*p->hull_len);
 292           if (!p->hull) return NLOPT_OUT_OF_MEMORY;
 293      }
 294      nhull = convex_hull(&p->rtree, hull = p->hull);
 295  divisions:
 296      for (i = 0; i < nhull; ++i) {
 297           double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
 298           if (i > 0)
 299                K1 = (hull[i][1] - hull[i-1][1]) / (hull[i][0] - hull[i-1][0]);
 300           if (i < nhull-1)
 301                K1 = (hull[i][1] - hull[i+1][1]) / (hull[i][0] - hull[i+1][0]);
 302           K = MAX(K1, K2);
 303           if (hull[i][1] - K * hull[i][0]
 304               <= p->fmin - magic_eps * fabs(p->fmin)) {
 305                /* "potentially optimal" rectangle, so subdivide */
 306                divided_some = 1;
 307                nlopt_result ret = divide_rect(hull[i], p);
 308                if (ret != NLOPT_SUCCESS) return ret;
 309                xtol_reached = xtol_reached && small(hull[i] + 2+n, p);
 310           }
 311      }
 312      if (!divided_some) {
 313           if (magic_eps != 0) {
 314                magic_eps = 0;
 315                goto divisions; /* try again */
 316           }
 317           else { /* WTF? divide largest rectangle with smallest f */
 318                /* (note that this code actually gets called from time
 319                   to time, and the heuristic here seems to work well,
 320                   but I don't recall this situation being discussed in
 321                   the references?) */
 322                rb_node *max = rb_tree_max(&p->rtree);
 323                rb_node *pred = max;
 324                double wmax = max->k[0];
 325                do { /* note: this loop is O(N) worst-case time */
 326                     max = pred;
 327                     pred = rb_tree_pred(max);
 328                } while (pred && pred->k[0] == wmax);
 329                return divide_rect(max->k, p);
 330           }
 331      }
 332      return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
 333 }
 334
 335 /***************************************************************************/
 336
 337 /* lexographic sort order (d,f) of hyper-rects, for red-black tree */
 338 static int hyperrect_compare(double *a, double *b)
 339 {
 340      if (a[0] < b[0]) return -1;
 341      if (a[0] > b[0]) return +1;
 342      if (a[1] < b[1]) return -1;
 343      if (a[1] > b[1]) return +1;
 344      return (int) (a - b); /* tie-breaker */
 345 }
 346
 347 /***************************************************************************/
 348
 349 nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
 350                               const double *lb, const double *ub,
 351                               double *x,
 352                               double *fmin,
 353                               nlopt_stopping *stop,
 354                               double magic_eps, int which_alg)
 355 {
 356      params p;
 357      int i, x_center = 1;
 358      double *rnew;
 359      nlopt_result ret = NLOPT_OUT_OF_MEMORY;
 360
 361      p.magic_eps = magic_eps;
 362      p.which_diam = which_alg % 2;
 363      p.which_div = 0;
 364      p.lb = lb; p.ub = ub;
 365      p.stop = stop;
 366      p.n = n;
 367      p.L = RECT_LEN(n);
 368      p.f = f;
 369      p.f_data = f_data;
 370      p.xmin = x;
 371      p.fmin = f(n, x, NULL, f_data); stop->nevals++;
 372      p.work = 0;
 373      p.iwork = 0;
 374      p.hull = 0;
 375
 376      rb_tree_init(&p.rtree, hyperrect_compare);
 377
 378      p.work = (double *) malloc(sizeof(double) * (2*n));
 379      if (!p.work) goto done;
 380      p.iwork = (int *) malloc(sizeof(int) * n);
 381      if (!p.iwork) goto done;
 382      p.hull_len = 128; /* start with a reasonable number */
 383      p.hull = (double **) malloc(sizeof(double *) * p.hull_len);
 384      if (!p.hull) goto done;
 385
 386      if (!(rnew = (double *) malloc(sizeof(double) * p.L))) goto done;
 387      for (i = 0; i < n; ++i) {
 388           rnew[2+i] = 0.5 * (lb[i] + ub[i]);
 389           x_center = x_center
 390                && (fabs(rnew[2+i]-x[i]) < 1e-13*(1+fabs(x[i])));
 391           rnew[2+n+i] = ub[i] - lb[i];
 392      }
 393      rnew[0] = rect_diameter(n, rnew+2+n, &p);
 394      if (x_center)
 395           rnew[1] = p.fmin; /* avoid computing f(center) twice */
 396      else
 397           rnew[1] = function_eval(rnew+2, &p);
 398      if (!rb_tree_insert(&p.rtree, rnew)) {
 399           free(rnew);
 400           goto done;
 401      }
 402
 403      ret = divide_rect(rnew, &p);
 404      if (ret != NLOPT_SUCCESS) goto done;
 405
 406      while (1) {
 407           double fmin0 = p.fmin;
 408           ret = divide_good_rects(&p);
 409           if (ret != NLOPT_SUCCESS) goto done;
 410           if (nlopt_stop_f(p.stop, p.fmin, fmin0)) {
 411                ret = NLOPT_FTOL_REACHED;
 412                goto done;
 413           }
 414      }
 415
 416  done:
 417      rb_tree_destroy_with_keys(&p.rtree);
 418      free(p.hull);
 419      free(p.iwork);
 420      free(p.work);
 421
 422      *fmin = p.fmin;
 423      return ret;
 424 }
 425
 426 /* in the conventional DIRECT-type algorithm, we first rescale our
 427    coordinates to a unit hypercube ... we do this simply by
 428    wrapping cdirect() around cdirect_unscaled(). */
 429
 430 typedef struct {
 431      nlopt_func f;
 432      void *f_data;
 433      double *x;
 434      const double *lb, *ub;
 435 } uf_data;
 436 static double uf(int n, const double *xu, double *grad, void *d_)
 437 {
 438      uf_data *d = (uf_data *) d_;
 439      double f;
 440      int i;
 441      for (i = 0; i < n; ++i)
 442           d->x[i] = d->lb[i] + xu[i] * (d->ub[i] - d->lb[i]);
 443      f = d->f(n, d->x, grad, d->f_data);
 444      if (grad)
 445           for (i = 0; i < n; ++i)
 446                grad[i] *= d->ub[i] - d->lb[i];
 447      return f;
 448 }
 449
 450 nlopt_result cdirect(int n, nlopt_func f, void *f_data,
 451                      const double *lb, const double *ub,
 452                      double *x,
 453                      double *fmin,
 454                      nlopt_stopping *stop,
 455                      double magic_eps, int which_alg)
 456 {
 457      uf_data d;
 458      nlopt_result ret;
 459      const double *xtol_abs_save;
 460      int i;
 461
 462      d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
 463      d.x = (double *) malloc(sizeof(double) * n*4);
 464      if (!d.x) return NLOPT_OUT_OF_MEMORY;
 465
 466      for (i = 0; i < n; ++i) {
 467           x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
 468           d.x[n+i] = 0;
 469           d.x[2*n+i] = 1;
 470           d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
 471      }
 472      xtol_abs_save = stop->xtol_abs;
 473      stop->xtol_abs = d.x + 3*n;
 474      ret = cdirect_unscaled(n, uf, &d, d.x+n, d.x+2*n, x, fmin, stop,
 475                             magic_eps, which_alg);
 476      stop->xtol_abs = xtol_abs_save;
 477      for (i = 0; i < n; ++i)
 478           x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);
 479      free(d.x);
 480      return ret;
 481 }