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      hull[nhull++] = n->k;
 249      if (xmin == xmax) return nhull;
 250
 251      /* set nmax = min mode with x == xmax */
 252 #if 0
 253      while (nmax->k[0] == xmax)
 254           nmax = rb_tree_pred(nmax); /* non-NULL since xmin != xmax */
 255      nmax = rb_tree_succ(nmax);
 256 #else
 257      /* performance hack (see also below) */
 258      {
 259           double kshift[2];
 260           kshift[0] = xmax * (1 - 1e-13);
 261           kshift[1] = -HUGE_VAL;
 262           nmax = rb_tree_find_gt(t, kshift); /* non-NULL since xmin != xmax */
 263      }
 264 #endif
 265
 266      ymaxmin = nmax->k[1];
 267      minslope = (ymaxmin - yminmin) / (xmax - xmin);
 268
 269      /* set n = first node with x != xmin */
 270 #if 0
 271      while (n->k[0] == xmin)
 272           n = rb_tree_succ(n); /* non-NULL since xmin != xmax */
 273 #else
 274      /* performance hack (see also below) */
 275      {
 276           double kshift[2];
 277           kshift[0] = xmin * (1 + 1e-13);
 278           kshift[1] = -HUGE_VAL;
 279           n = rb_tree_find_gt(t, kshift); /* non-NULL since xmin != xmax */
 280      }
 281 #endif
 282
 283      for (; n != nmax; n = rb_tree_succ(n)) {
 284           double *k = n->k;
 285           if (k[1] > yminmin + (k[0] - xmin) * minslope)
 286                continue;
 287
 288           /* performance hack: most of the points in DIRECT lie along
 289              vertical lines at a few x values, and we can exploit this */
 290           if (nhull && k[0] == hull[nhull - 1][0]) { /* x == previous x */
 291                if (k[1] > hull[nhull - 1][1]) {
 292                     double kshift[2];
 293                     /* because of the round to float in rect_diameter, above,
 294                        it shouldn't be possible for two diameters (x values)
 295                        to have a fractional difference < 1e-13.  Note
 296                        that k[0] > 0 always in DIRECT */
 297                     kshift[0] = k[0] * (1 + 1e-13);
 298                     kshift[1] = -HUGE_VAL;
 299                     n = rb_tree_pred(rb_tree_find_gt(t, kshift));
 300                     continue;
 301                }
 302                else { /* equal y values, add to hull */
 303                     hull[nhull++] = k;
 304                     continue;
 305                }
 306           }
 307
 308           /* remove points until we are making a "left turn" to k */
 309           while (nhull > 1) {
 310                double *t1 = hull[nhull - 1], *t2;
 311
 312                /* because we allow equal points in our hull, we have
 313                   to modify the standard convex-hull algorithm slightly:
 314                   we need to look backwards in the hull list until we
 315                   find a point t2 != t1 */
 316                int it2 = nhull - 2;
 317                do {
 318                     t2 = hull[it2--];
 319                } while (it2 >= 0 && t2[0] == t1[0] && t2[1] == t1[1]);
 320                if (it2 < 0) break;
 321
 322                /* cross product (t1-t2) x (k-t2) > 0 for a left turn: */
 323                if ((t1[0]-t2[0]) * (k[1]-t2[1])
 324                    - (t1[1]-t2[1]) * (k[0]-t2[0]) >= 0)
 325                     break;
 326                --nhull;
 327           }
 328           hull[nhull++] = k;
 329      }
 330      hull[nhull++] = nmax->k;
 331      return nhull;
 332 }
 333
 334 /***************************************************************************/
 335
 336 static int small(double *w, params *p)
 337 {
 338      int i;
 339      for (i = 0; i < p->n; ++i)
 340           if (w[i] > p->stop->xtol_abs[i] &&
 341               w[i] > (p->ub[i] - p->lb[i]) * p->stop->xtol_rel)
 342                return 0;
 343      return 1;
 344 }
 345
 346 static nlopt_result divide_good_rects(params *p)
 347 {
 348      const int n = p->n;
 349      double **hull;
 350      int nhull, i, xtol_reached = 1, divided_some = 0;
 351      double magic_eps = p->magic_eps;
 352
 353      if (p->hull_len < p->rtree.N) {
 354           p->hull_len += p->rtree.N;
 355           p->hull = (double **) realloc(p->hull, sizeof(double*)*p->hull_len);
 356           if (!p->hull) return NLOPT_OUT_OF_MEMORY;
 357      }
 358      nhull = convex_hull(&p->rtree, hull = p->hull);
 359  divisions:
 360      for (i = 0; i < nhull; ++i) {
 361           double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
 362           if (i > 0)
 363                K1 = (hull[i][1] - hull[i-1][1]) / (hull[i][0] - hull[i-1][0]);
 364           if (i < nhull-1)
 365                K1 = (hull[i][1] - hull[i+1][1]) / (hull[i][0] - hull[i+1][0]);
 366           K = MAX(K1, K2);
 367           if (hull[i][1] - K * hull[i][0]
 368               <= p->fmin - magic_eps * fabs(p->fmin)) {
 369                /* "potentially optimal" rectangle, so subdivide */
 370                divided_some = 1;
 371                nlopt_result ret = divide_rect(hull[i], p);
 372                if (ret != NLOPT_SUCCESS) return ret;
 373                xtol_reached = xtol_reached && small(hull[i] + 2+n, p);
 374           }
 375      }
 376      if (!divided_some) {
 377           if (magic_eps != 0) {
 378                magic_eps = 0;
 379                goto divisions; /* try again */
 380           }
 381           else { /* WTF? divide largest rectangle with smallest f */
 382                /* (note that this code actually gets called from time
 383                   to time, and the heuristic here seems to work well,
 384                   but I don't recall this situation being discussed in
 385                   the references?) */
 386                rb_node *max = rb_tree_max(&p->rtree);
 387                rb_node *pred = max;
 388                double wmax = max->k[0];
 389                do { /* note: this loop is O(N) worst-case time */
 390                     max = pred;
 391                     pred = rb_tree_pred(max);
 392                } while (pred && pred->k[0] == wmax);
 393                return divide_rect(max->k, p);
 394           }
 395      }
 396      return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
 397 }
 398
 399 /***************************************************************************/
 400
 401 /* lexographic sort order (d,f) of hyper-rects, for red-black tree */
 402 static int hyperrect_compare(double *a, double *b)
 403 {
 404      if (a[0] < b[0]) return -1;
 405      if (a[0] > b[0]) return +1;
 406      if (a[1] < b[1]) return -1;
 407      if (a[1] > b[1]) return +1;
 408      return (int) (a - b); /* tie-breaker */
 409 }
 410
 411 /***************************************************************************/
 412
 413 nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
 414                               const double *lb, const double *ub,
 415                               double *x,
 416                               double *fmin,
 417                               nlopt_stopping *stop,
 418                               double magic_eps, int which_alg)
 419 {
 420      params p;
 421      int i, x_center = 1;
 422      double *rnew;
 423      nlopt_result ret = NLOPT_OUT_OF_MEMORY;
 424
 425      p.magic_eps = magic_eps;
 426      p.which_diam = which_alg % 2;
 427      p.which_div = 0;
 428      p.lb = lb; p.ub = ub;
 429      p.stop = stop;
 430      p.n = n;
 431      p.L = RECT_LEN(n);
 432      p.f = f;
 433      p.f_data = f_data;
 434      p.xmin = x;
 435      p.fmin = f(n, x, NULL, f_data); stop->nevals++;
 436      p.work = 0;
 437      p.iwork = 0;
 438      p.hull = 0;
 439
 440      rb_tree_init(&p.rtree, hyperrect_compare);
 441
 442      p.work = (double *) malloc(sizeof(double) * (2*n));
 443      if (!p.work) goto done;
 444      p.iwork = (int *) malloc(sizeof(int) * n);
 445      if (!p.iwork) goto done;
 446      p.hull_len = 128; /* start with a reasonable number */
 447      p.hull = (double **) malloc(sizeof(double *) * p.hull_len);
 448      if (!p.hull) goto done;
 449
 450      if (!(rnew = (double *) malloc(sizeof(double) * p.L))) goto done;
 451      for (i = 0; i < n; ++i) {
 452           rnew[2+i] = 0.5 * (lb[i] + ub[i]);
 453           x_center = x_center
 454                && (fabs(rnew[2+i]-x[i]) < 1e-13*(1+fabs(x[i])));
 455           rnew[2+n+i] = ub[i] - lb[i];
 456      }
 457      rnew[0] = rect_diameter(n, rnew+2+n, &p);
 458      if (x_center)
 459           rnew[1] = p.fmin; /* avoid computing f(center) twice */
 460      else
 461           rnew[1] = function_eval(rnew+2, &p);
 462      if (!rb_tree_insert(&p.rtree, rnew)) {
 463           free(rnew);
 464           goto done;
 465      }
 466
 467      ret = divide_rect(rnew, &p);
 468      if (ret != NLOPT_SUCCESS) goto done;
 469
 470      while (1) {
 471           double fmin0 = p.fmin;
 472           ret = divide_good_rects(&p);
 473           if (ret != NLOPT_SUCCESS) goto done;
 474           if (p.fmin < fmin0 && nlopt_stop_f(p.stop, p.fmin, fmin0)) {
 475                ret = NLOPT_FTOL_REACHED;
 476                goto done;
 477           }
 478      }
 479
 480  done:
 481      rb_tree_destroy_with_keys(&p.rtree);
 482      free(p.hull);
 483      free(p.iwork);
 484      free(p.work);
 485
 486      *fmin = p.fmin;
 487      return ret;
 488 }
 489
 490 /* in the conventional DIRECT-type algorithm, we first rescale our
 491    coordinates to a unit hypercube ... we do this simply by
 492    wrapping cdirect() around cdirect_unscaled(). */
 493
 494 typedef struct {
 495      nlopt_func f;
 496      void *f_data;
 497      double *x;
 498      const double *lb, *ub;
 499 } uf_data;
 500 static double uf(int n, const double *xu, double *grad, void *d_)
 501 {
 502      uf_data *d = (uf_data *) d_;
 503      double f;
 504      int i;
 505      for (i = 0; i < n; ++i)
 506           d->x[i] = d->lb[i] + xu[i] * (d->ub[i] - d->lb[i]);
 507      f = d->f(n, d->x, grad, d->f_data);
 508      if (grad)
 509           for (i = 0; i < n; ++i)
 510                grad[i] *= d->ub[i] - d->lb[i];
 511      return f;
 512 }
 513
 514 nlopt_result cdirect(int n, nlopt_func f, void *f_data,
 515                      const double *lb, const double *ub,
 516                      double *x,
 517                      double *fmin,
 518                      nlopt_stopping *stop,
 519                      double magic_eps, int which_alg)
 520 {
 521      uf_data d;
 522      nlopt_result ret;
 523      const double *xtol_abs_save;
 524      int i;
 525
 526      d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
 527      d.x = (double *) malloc(sizeof(double) * n*4);
 528      if (!d.x) return NLOPT_OUT_OF_MEMORY;
 529
 530      for (i = 0; i < n; ++i) {
 531           x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
 532           d.x[n+i] = 0;
 533           d.x[2*n+i] = 1;
 534           d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
 535      }
 536      xtol_abs_save = stop->xtol_abs;
 537      stop->xtol_abs = d.x + 3*n;
 538      ret = cdirect_unscaled(n, uf, &d, d.x+n, d.x+2*n, x, fmin, stop,
 539                             magic_eps, which_alg);
 540      stop->xtol_abs = xtol_abs_save;
 541      for (i = 0; i < n; ++i)
 542           x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);
 543      free(d.x);
 544      return ret;
 545 }