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  * 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).
  26  *
  27  * n > 0 always
  28  */
  29
  30 #define RECT_LEN(n) (2*(n)+2) /* number of double values in a hyperrect */
  31
  32 /* parameters of the search algorithm and various information that
  33    needs to be passed around */
  34 typedef struct {
  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 */
  50
  51      /* red-black tree of hyperrects, sorted by (d,f) in
  52         lexographical order */
  53      rb_tree rtree;
  54      double **hull; /* array to store convex hull */
  55      int hull_len; /* allocated length of hull array */
  56 } params;
  57
  58 /***************************************************************************/
  59
  60 /* Evaluate the "diameter" (d) of a rectangle of widths w[n]
  61
  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)
  68 {
  69      int i;
  70      if (p->which_diam == 0) { /* Jones measure */
  71           double sum = 0;
  72           for (i = 0; i < n; ++i)
  73                sum += w[i] * w[i];
  74           /* distance from center to a vertex */
  75           return ((float) (sqrt(sum) * 0.5));
  76      }
  77      else { /* Gablonsky measure */
  78           double maxw = 0;
  79           for (i = 0; i < n; ++i)
  80                if (w[i] > maxw)
  81                     maxw = w[i];
  82           /* half-width of longest side */
  83           return ((float) (maxw * 0.5));
  84      }
  85 }
  86
  87 #define ALLOC_RECT(rect, L) if (!(rect = (double*) malloc(sizeof(double)*(L)))) return NLOPT_OUT_OF_MEMORY
  88
  89 static double *fv_qsort = 0;
  90 static int sort_fv_compare(const void *a_, const void *b_)
  91 {
  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]);
  95      if (fa < fb)
  96           return -1;
  97      else if (fa > fb)
  98           return +1;
  99      else
 100           return 0;
 101 }
 102 static void sort_fv(int n, double *fv, int *isort)
 103 {
 104      int i;
 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);
 108      fv_qsort = 0;
 109 }
 110
 111 static double function_eval(const double *x, params *p) {
 112      double f = p->f(p->n, x, NULL, p->f_data);
 113      if (f < p->fmin) {
 114           p->fmin = f;
 115           memcpy(p->xmin, x, sizeof(double) * p->n);
 116      }
 117      p->stop->nevals++;
 118      return f;
 119 }
 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; }
 121
 122 #define THIRD (0.3333333333333333333333)
 123
 124 #define EQUAL_SIDE_TOL 5e-2 /* tolerance to equate side sizes */
 125
 126 /* divide rectangle idiv in the list p->rects */
 127 static nlopt_result divide_rect(double *rdiv, params *p)
 128 {
 129      int i;
 130      const const int n = p->n;
 131      const int L = p->L;
 132      double *c = rdiv + 2; /* center of rect to divide */
 133      double *w = c + n; /* widths of rect to divide */
 134      double wmax = w[0];
 135      int imax = 0, nlongest = 0;
 136      rb_node *node;
 137
 138      for (i = 1; i < n; ++i)
 139           if (w[i] > wmax)
 140                wmax = w[imax = i];
 141      for (i = 0; i < n; ++i)
 142           if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL)
 143                ++nlongest;
 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) {
 151                     double csave = c[i];
 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);
 156                     c[i] = csave;
 157                }
 158                else {
 159                     fv[2*i] = fv[2*i+1] = HUGE_VAL;
 160                }
 161           }
 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) {
 166                int k;
 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) {
 171                     double *rnew;
 172                     ALLOC_RECT(rnew, L);
 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)) {
 177                          free(rnew);
 178                          return NLOPT_OUT_OF_MEMORY;
 179                     }
 180                }
 181           }
 182      }
 183      else {
 184           int k;
 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; }
 191                          --i;
 192                     }
 193           }
 194           else
 195                i = imax; /* trisect longest side */
 196           if (!(node = rb_tree_find(&p->rtree, rdiv)))
 197                return NLOPT_FAILURE;
 198           w[i] *= THIRD;
 199           rdiv[0] = rect_diameter(n, w, p);
 200           node = rb_tree_resort(&p->rtree, node);
 201           for (k = 0; k <= 1; ++k) {
 202                double *rnew;
 203                ALLOC_RECT(rnew, L);
 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)) {
 208                     free(rnew);
 209                     return NLOPT_OUT_OF_MEMORY;
 210                }
 211           }
 212      }
 213      return NLOPT_SUCCESS;
 214 }
 215
 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. */
 222
 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).
 225
 226    Unlike standard convex hulls, we allow redundant points on the hull.
 227
 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).
 230 */
 231 static int convex_hull(rb_tree *t, double **hull)
 232 {
 233      int nhull = 0;
 234      double minslope;
 235      double xmin, xmax, yminmin, ymaxmin;
 236      rb_node *n, *nmax;
 237
 238      /* Monotone chain algorithm [Andrew, 1979]. */
 239
 240      n = rb_tree_min(t);
 241      if (!n) return 0;
 242      nmax = rb_tree_max(t);
 243
 244      xmin = n->k[0];
 245      yminmin = n->k[1];
 246      xmax = nmax->k[0];
 247
 248      do { /* include any duplicate points at (xmin,yminmin) */
 249           hull[nhull++] = n->k;
 250           n = rb_tree_succ(n);
 251      } while (n && n->k[0] == xmin && n->k[1] == yminmin);
 252      if (xmin == xmax) return nhull;
 253
 254      /* set nmax = min mode with x == xmax */
 255 #if 0
 256      while (nmax->k[0] == xmax)
 257           nmax = rb_tree_pred(nmax); /* non-NULL since xmin != xmax */
 258      nmax = rb_tree_succ(nmax);
 259 #else
 260      /* performance hack (see also below) */
 261      {
 262           double kshift[2];
 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 */
 266      }
 267 #endif
 268
 269      ymaxmin = nmax->k[1];
 270      minslope = (ymaxmin - yminmin) / (xmax - xmin);
 271
 272      /* set n = first node with x != xmin */
 273 #if 0
 274      while (n->k[0] == xmin)
 275           n = rb_tree_succ(n); /* non-NULL since xmin != xmax */
 276 #else
 277      /* performance hack (see also below) */
 278      {
 279           double kshift[2];
 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 */
 283      }
 284 #endif
 285
 286      for (; n != nmax; n = rb_tree_succ(n)) {
 287           double *k = n->k;
 288           if (k[1] > yminmin + (k[0] - xmin) * minslope)
 289                continue;
 290
 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]) {
 295                     double kshift[2];
 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));
 303                     continue;
 304                }
 305                else { /* equal y values, add to hull */
 306                     hull[nhull++] = k;
 307                     continue;
 308                }
 309           }
 310
 311           /* remove points until we are making a "left turn" to k */
 312           while (nhull > 1) {
 313                double *t1 = hull[nhull - 1], *t2;
 314
 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 */
 319                int it2 = nhull - 2;
 320                do {
 321                     t2 = hull[it2--];
 322                } while (it2 >= 0 && t2[0] == t1[0] && t2[1] == t1[1]);
 323                if (it2 < 0) break;
 324
 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)
 328                     break;
 329                --nhull;
 330           }
 331           hull[nhull++] = k;
 332      }
 333
 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);
 338
 339      return nhull;
 340 }
 341
 342 /***************************************************************************/
 343
 344 static int small(double *w, params *p)
 345 {
 346      int i;
 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)
 350                return 0;
 351      return 1;
 352 }
 353
 354 static nlopt_result divide_good_rects(params *p)
 355 {
 356      const int n = p->n;
 357      double **hull;
 358      int nhull, i, xtol_reached = 1, divided_some = 0;
 359      double magic_eps = p->magic_eps;
 360
 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;
 365      }
 366      nhull = convex_hull(&p->rtree, hull = p->hull);
 367  divisions:
 368      for (i = 0; i < nhull; ++i) {
 369           double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
 370           int im, ip;
 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);
 373           if (im >= 0)
 374                K1 = (hull[i][1] - hull[im][1]) / (hull[i][0] - hull[im][0]);
 375           if (ip < nhull)
 376                K1 = (hull[i][1] - hull[ip][1]) / (hull[i][0] - hull[ip][0]);
 377           K = MAX(K1, K2);
 378           if (hull[i][1] - K * hull[i][0]
 379               <= p->fmin - magic_eps * fabs(p->fmin)) {
 380                /* "potentially optimal" rectangle, so subdivide */
 381                nlopt_result ret = divide_rect(hull[i], p);
 382                divided_some = 1;
 383                if (ret != NLOPT_SUCCESS) return ret;
 384                xtol_reached = xtol_reached && small(hull[i] + 2+n, p);
 385           }
 386      }
 387      if (!divided_some) {
 388           if (magic_eps != 0) {
 389                magic_eps = 0;
 390                goto divisions; /* try again */
 391           }
 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
 396                   the references?) */
 397                rb_node *max = rb_tree_max(&p->rtree);
 398                rb_node *pred = max;
 399                double wmax = max->k[0];
 400                do { /* note: this loop is O(N) worst-case time */
 401                     max = pred;
 402                     pred = rb_tree_pred(max);
 403                } while (pred && pred->k[0] == wmax);
 404                return divide_rect(max->k, p);
 405           }
 406      }
 407      return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
 408 }
 409
 410 /***************************************************************************/
 411
 412 /* lexographic sort order (d,f) of hyper-rects, for red-black tree */
 413 static int hyperrect_compare(double *a, double *b)
 414 {
 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 */
 420 }
 421
 422 /***************************************************************************/
 423
 424 nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
 425                               const double *lb, const double *ub,
 426                               double *x,
 427                               double *fmin,
 428                               nlopt_stopping *stop,
 429                               double magic_eps, int which_alg)
 430 {
 431      params p;
 432      int i, x_center = 1;
 433      double *rnew;
 434      nlopt_result ret = NLOPT_OUT_OF_MEMORY;
 435
 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;
 440      p.stop = stop;
 441      p.n = n;
 442      p.L = RECT_LEN(n);
 443      p.f = f;
 444      p.f_data = f_data;
 445      p.xmin = x;
 446      p.fmin = f(n, x, NULL, f_data); stop->nevals++;
 447      p.work = 0;
 448      p.iwork = 0;
 449      p.hull = 0;
 450
 451      rb_tree_init(&p.rtree, hyperrect_compare);
 452
 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;
 460
 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           x_center = x_center
 465                && (fabs(rnew[2+i]-x[i]) < 1e-13*(1+fabs(x[i])));
 466           rnew[2+n+i] = ub[i] - lb[i];
 467      }
 468      rnew[0] = rect_diameter(n, rnew+2+n, &p);
 469      if (x_center)
 470           rnew[1] = p.fmin; /* avoid computing f(center) twice */
 471      else
 472           rnew[1] = function_eval(rnew+2, &p);
 473      if (!rb_tree_insert(&p.rtree, rnew)) {
 474           free(rnew);
 475           goto done;
 476      }
 477
 478      ret = divide_rect(rnew, &p);
 479      if (ret != NLOPT_SUCCESS) goto done;
 480
 481      while (1) {
 482           double fmin0 = p.fmin;
 483           ret = divide_good_rects(&p);
 484           if (ret != NLOPT_SUCCESS) goto done;
 485           if (p.fmin < fmin0 && nlopt_stop_f(p.stop, p.fmin, fmin0)) {
 486                ret = NLOPT_FTOL_REACHED;
 487                goto done;
 488           }
 489      }
 490
 491  done:
 492      rb_tree_destroy_with_keys(&p.rtree);
 493      free(p.hull);
 494      free(p.iwork);
 495      free(p.work);
 496
 497      *fmin = p.fmin;
 498      return ret;
 499 }
 500
 501 /* in the conventional DIRECT-type algorithm, we first rescale our
 502    coordinates to a unit hypercube ... we do this simply by
 503    wrapping cdirect() around cdirect_unscaled(). */
 504
 505 typedef struct {
 506      nlopt_func f;
 507      void *f_data;
 508      double *x;
 509      const double *lb, *ub;
 510 } uf_data;
 511 static double uf(int n, const double *xu, double *grad, void *d_)
 512 {
 513      uf_data *d = (uf_data *) d_;
 514      double f;
 515      int i;
 516      for (i = 0; i < n; ++i)
 517           d->x[i] = d->lb[i] + xu[i] * (d->ub[i] - d->lb[i]);
 518      f = d->f(n, d->x, grad, d->f_data);
 519      if (grad)
 520           for (i = 0; i < n; ++i)
 521                grad[i] *= d->ub[i] - d->lb[i];
 522      return f;
 523 }
 524
 525 nlopt_result cdirect(int n, nlopt_func f, void *f_data,
 526                      const double *lb, const double *ub,
 527                      double *x,
 528                      double *fmin,
 529                      nlopt_stopping *stop,
 530                      double magic_eps, int which_alg)
 531 {
 532      uf_data d;
 533      nlopt_result ret;
 534      const double *xtol_abs_save;
 535      int i;
 536
 537      d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
 538      d.x = (double *) malloc(sizeof(double) * n*4);
 539      if (!d.x) return NLOPT_OUT_OF_MEMORY;
 540
 541      for (i = 0; i < n; ++i) {
 542           x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
 543           d.x[n+i] = 0;
 544           d.x[2*n+i] = 1;
 545           d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
 546      }
 547      xtol_abs_save = stop->xtol_abs;
 548      stop->xtol_abs = d.x + 3*n;
 549      ret = cdirect_unscaled(n, uf, &d, d.x+n, d.x+2*n, x, fmin, stop,
 550                             magic_eps, which_alg);
 551      stop->xtol_abs = xtol_abs_save;
 552      for (i = 0; i < n; ++i)
 553           x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);
 554      free(d.x);
 555      return ret;
 556 }