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