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