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 [0]
  18  * is the value (f) of the function at the center, [1] 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 *rects; /* the hyper-rectangles */
  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, iwork_len; /* workspace, of length iwork_len >= n */
  49      double fmin, *xmin; /* minimum so far */
  50
  51      /* red-black tree of hyperrect indices, sorted by (d,f) in
  52         lexographical order */
  53      rb_tree rtree;
  54 } params;
  55
  56 /***************************************************************************/
  57
  58 /* evaluate the "diameter" (d) of a rectangle of widths w[n] */
  59 static double rect_diameter(int n, const double *w, const params *p)
  60 {
  61      int i;
  62      if (p->which_diam == 0) { /* Jones measure */
  63           double sum = 0;
  64           for (i = 0; i < n; ++i)
  65                sum += w[i] * w[i];
  66           return sqrt(sum) * 0.5; /* distance from center to a vertex */
  67      }
  68      else { /* Gablonsky measure */
  69           double maxw = 0;
  70           for (i = 0; i < n; ++i)
  71                if (w[i] > maxw)
  72                     maxw = w[i];
  73           return maxw * 0.5; /* half-width of longest side */
  74      }
  75 }
  76
  77 static double *alloc_rects(int n, int *Na, double *rects, int newN)
  78 {
  79      if (newN <= *Na)
  80           return rects;
  81      else {
  82           (*Na) += newN;
  83           return realloc(rects, sizeof(double) * RECT_LEN(n) * (*Na));
  84      }
  85 }
  86 #define ALLOC_RECTS(n, Nap, rects, newN) if (!(rects = alloc_rects(n, Nap, rects, newN))) 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) fv = function_eval(x, p); if (p->fmin < p->stop->fmin_max) return NLOPT_FMIN_MAX_REACHED; else if (nlopt_stop_evals((p)->stop)) return NLOPT_MAXEVAL_REACHED; else if (nlopt_stop_time((p)->stop)) 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(int *N, int *Na, int idiv, params *p)
 127 {
 128      int i;
 129      const const int n = p->n;
 130      const int L = p->L;
 131      double *r = p->rects;
 132      double *c = r + L*idiv + 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);
 154                     c[i] = csave + w[i] * THIRD;
 155                     FUNCTION_EVAL(fv[2*i+1], c, p);
 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           ALLOC_RECTS(n, Na, r, (*N)+2*nlongest);
 164           p->rects = r; c = r + L*idiv + 2; w = c + n;
 165           for (i = 0; i < nlongest; ++i) {
 166                int k;
 167                if (!(node = rb_tree_find(&p->rtree, idiv)))
 168                     return NLOPT_FAILURE;
 169                w[isort[i]] *= THIRD;
 170                r[L*idiv + 1] = rect_diameter(n, w, p);
 171                rb_tree_resort(&p->rtree, node);
 172                for (k = 0; k <= 1; ++k) {
 173                     memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
 174                     r[L*(*N) + 2 + isort[i]] += w[isort[i]] * (2*k-1);
 175                     r[L*(*N)] = fv[2*isort[i]+k];
 176                     if (!rb_tree_insert(&p->rtree, *N))
 177                          return NLOPT_FAILURE;
 178                     ++(*N);
 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           ALLOC_RECTS(n, Na, r, (*N)+2);
 196           p->rects = r; c = r + L*idiv + 2; w = c + n;
 197           if (!(node = rb_tree_find(&p->rtree, idiv)))
 198                return NLOPT_FAILURE;
 199           w[i] *= THIRD;
 200           r[L*idiv + 1] = rect_diameter(n, w, p);
 201           rb_tree_resort(&p->rtree, node);
 202           for (k = 0; k <= 1; ++k) {
 203                memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
 204                r[L*(*N) + 2 + i] += w[i] * (2*k-1);
 205                FUNCTION_EVAL(r[L*(*N)], r + L*(*N) + 2, p);
 206                if (!rb_tree_insert(&p->rtree, *N))
 207                     return NLOPT_FAILURE;
 208                ++(*N);
 209           }
 210      }
 211      return NLOPT_SUCCESS;
 212 }
 213
 214 /***************************************************************************/
 215 /* O(N log N) convex hull algorithm, used later to find the potentially
 216    optimal points */
 217
 218 /* Find the lower convex hull of a set of points (xy[s*i+1], xy[s*i]), where
 219    0 <= i < N and s >= 2.
 220
 221    The return value is the number of points in the hull, with indices
 222    stored in ihull.  ihull should point to arrays of length >= N.
 223    rb_tree should be a red-black tree of indices (keys == i) sorted
 224    in lexographic order by (xy[s*i+1], xy[s*i]).
 225 */
 226 static int convex_hull(int N, double *xy, int s, int *ihull, rb_tree *t)
 227 {
 228      int nhull;
 229      double minslope;
 230      double xmin, xmax, yminmin, ymaxmin;
 231      rb_node *n, *nmax;
 232
 233      if (N == 0) return 0;
 234
 235      /* Monotone chain algorithm [Andrew, 1979]. */
 236
 237      n = rb_tree_min(t);
 238      nmax = rb_tree_max(t);
 239
 240      xmin = xy[s*(n->k)+1];
 241      yminmin = xy[s*(n->k)];
 242      xmax = xy[s*(nmax->k)+1];
 243
 244      ihull[nhull = 1] = n->k;
 245      if (xmin == xmax) return nhull;
 246
 247      /* set nmax = min mode with x == xmax */
 248      while (xy[s*(nmax->k)+1] == xmax)
 249           nmax = rb_tree_pred(t, nmax); /* non-NULL since xmin != xmax */
 250      nmax = rb_tree_succ(t, nmax);
 251
 252      ymaxmin = xy[s*(nmax->k)];
 253      minslope = (ymaxmin - yminmin) / (xmax - xmin);
 254
 255      /* set n = first node with x != xmin */
 256      while (xy[s*(n->k)+1] == xmin)
 257           n = rb_tree_succ(t, n); /* non-NULL since xmin != xmax */
 258
 259      for (; n != nmax; n = rb_tree_succ(t, n)) {
 260           int k = n->k;
 261           if (xy[s*k] > yminmin + (xy[s*k+1] - xmin) * minslope)
 262                continue;
 263           /* remove points until we are making a "left turn" to k */
 264           while (nhull > 1) {
 265                int t1 = ihull[nhull - 1], t2 = ihull[nhull - 2];
 266                /* cross product (t1-t2) x (k-t2) > 0 for a left turn: */
 267                if ((xy[s*t1+1]-xy[s*t2+1]) * (xy[s*k]-xy[s*t2])
 268                    - (xy[s*t1]-xy[s*t2]) * (xy[s*k+1]-xy[s*t2+1]) > 0)
 269                     break;
 270                --nhull;
 271           }
 272           ihull[nhull++] = k;
 273      }
 274      ihull[nhull++] = nmax->k;
 275      return nhull;
 276 }
 277
 278 /***************************************************************************/
 279
 280 static int small(double *w, params *p)
 281 {
 282      int i;
 283      for (i = 0; i < p->n; ++i)
 284           if (w[i] > p->stop->xtol_abs[i] &&
 285               w[i] > (p->ub[i] - p->lb[i]) * p->stop->xtol_rel)
 286                return 0;
 287      return 1;
 288 }
 289
 290 static nlopt_result divide_good_rects(int *N, int *Na, params *p)
 291 {
 292      const int n = p->n;
 293      const int L = p->L;
 294      int *ihull, nhull, i, xtol_reached = 1, divided_some = 0;
 295      double *r = p->rects;
 296      double magic_eps = p->magic_eps;
 297
 298      if (p->iwork_len < *N) {
 299           p->iwork_len = p->iwork_len + *N;
 300           p->iwork = (int *) realloc(p->iwork, sizeof(int) * p->iwork_len);
 301           if (!p->iwork)
 302                return NLOPT_OUT_OF_MEMORY;
 303      }
 304      ihull = p->iwork;
 305      nhull = convex_hull(*N, r, L, ihull, &p->rtree);
 306  divisions:
 307      for (i = 0; i < nhull; ++i) {
 308           double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
 309           if (i > 0)
 310                K1 = (r[L*ihull[i]] - r[L*ihull[i-1]]) /
 311                     (r[L*ihull[i]+1] - r[L*ihull[i-1]+1]);
 312           if (i < nhull-1)
 313                K1 = (r[L*ihull[i]] - r[L*ihull[i+1]]) /
 314                     (r[L*ihull[i]+1] - r[L*ihull[i+1]+1]);
 315           K = MAX(K1, K2);
 316           if (r[L*ihull[i]] - K * r[L*ihull[i]+1]
 317               <= p->fmin - magic_eps * fabs(p->fmin)) {
 318                /* "potentially optimal" rectangle, so subdivide */
 319                divided_some = 1;
 320                nlopt_result ret;
 321                ret = divide_rect(N, Na, ihull[i], p);
 322                r = p->rects; /* may have grown */
 323                if (ret != NLOPT_SUCCESS) return ret;
 324                xtol_reached = xtol_reached && small(r + L*ihull[i] + 2+n, p);
 325           }
 326      }
 327      if (!divided_some) {
 328           if (magic_eps != 0) {
 329                magic_eps = 0;
 330                goto divisions; /* try again */
 331           }
 332           else { /* WTF? divide largest rectangle */
 333                double wmax = r[1];
 334                int imax = 0;
 335                for (i = 1; i < *N; ++i)
 336                     if (r[L*i+1] > wmax)
 337                          wmax = r[L*(imax=i)+1];
 338                return divide_rect(N, Na, imax, p);
 339           }
 340      }
 341      return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
 342 }
 343
 344 /***************************************************************************/
 345
 346 /* lexographic sort order (d,f) of hyper-rects, for red-black tree */
 347 static int hyperrect_compare(int i, int j, void *p_)
 348 {
 349      params *p = (params *) p_;
 350      int L = p->L;
 351      double *r = p->rects;
 352      double di = r[i*L+1], dj = r[j*L+1], fi, fj;
 353      if (di < dj) return -1;
 354      if (dj < di) return +1;
 355      fi = r[i*L]; fj = r[j*L];
 356      if (fi < fj) return -1;
 357      if (fj < fi) return +1;
 358      return 0;
 359 }
 360
 361 /***************************************************************************/
 362
 363 nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
 364                               const double *lb, const double *ub,
 365                               double *x,
 366                               double *fmin,
 367                               nlopt_stopping *stop,
 368                               double magic_eps, int which_alg)
 369 {
 370      params p;
 371      int Na = 100, N = 1, i, x_center = 1;
 372      nlopt_result ret = NLOPT_OUT_OF_MEMORY;
 373
 374      p.magic_eps = magic_eps;
 375      p.which_diam = which_alg % 2;
 376      p.which_div = 0;
 377      p.lb = lb; p.ub = ub;
 378      p.stop = stop;
 379      p.n = n;
 380      p.L = RECT_LEN(n);
 381      p.f = f;
 382      p.f_data = f_data;
 383      p.xmin = x;
 384      p.fmin = f(n, x, NULL, f_data); stop->nevals++;
 385      p.work = 0;
 386      p.iwork = 0;
 387      p.rects = 0;
 388
 389      if (!rb_tree_init(&p.rtree, hyperrect_compare, &p)) goto done;
 390      p.work = (double *) malloc(sizeof(double) * 2*n);
 391      if (!p.work) goto done;
 392      p.rects = (double *) malloc(sizeof(double) * Na * RECT_LEN(n));
 393      if (!p.rects) goto done;
 394      p.iwork = (int *) malloc(sizeof(int) * (p.iwork_len = Na + n));
 395      if (!p.iwork) goto done;
 396
 397      for (i = 0; i < n; ++i) {
 398           p.rects[2+i] = 0.5 * (lb[i] + ub[i]);
 399           x_center = x_center
 400                && (fabs(p.rects[2+i]-x[i]) < 1e-13*(1+fabs(x[i])));
 401           p.rects[2+n+i] = ub[i] - lb[i];
 402      }
 403      p.rects[1] = rect_diameter(n, p.rects+2+n, &p);
 404      if (x_center)
 405           p.rects[0] = p.fmin; /* avoid computing f(center) twice */
 406      else
 407           p.rects[0] = function_eval(p.rects+2, &p);
 408      if (!rb_tree_insert(&p.rtree, 0)) {
 409           ret = NLOPT_FAILURE;
 410           goto done;
 411      }
 412
 413      ret = divide_rect(&N, &Na, 0, &p);
 414      if (ret != NLOPT_SUCCESS) goto done;
 415
 416      while (1) {
 417           double fmin0 = p.fmin;
 418           ret = divide_good_rects(&N, &Na, &p);
 419           if (ret != NLOPT_SUCCESS) goto done;
 420           if (nlopt_stop_f(p.stop, p.fmin, fmin0)) {
 421                ret = NLOPT_FTOL_REACHED;
 422                goto done;
 423           }
 424      }
 425
 426  done:
 427      rb_tree_destroy(&p.rtree);
 428      free(p.iwork);
 429      free(p.rects);
 430      free(p.work);
 431
 432      *fmin = p.fmin;
 433      return ret;
 434 }
 435
 436 /* in the conventional DIRECT-type algorithm, we first rescale our
 437    coordinates to a unit hypercube ... we do this simply by
 438    wrapping cdirect() around cdirect_unscaled(). */
 439
 440 typedef struct {
 441      nlopt_func f;
 442      void *f_data;
 443      double *x;
 444      const double *lb, *ub;
 445 } uf_data;
 446 static double uf(int n, const double *xu, double *grad, void *d_)
 447 {
 448      uf_data *d = (uf_data *) d_;
 449      double f;
 450      int i;
 451      for (i = 0; i < n; ++i)
 452           d->x[i] = d->lb[i] + xu[i] * (d->ub[i] - d->lb[i]);
 453      f = d->f(n, d->x, grad, d->f_data);
 454      if (grad)
 455           for (i = 0; i < n; ++i)
 456                grad[i] *= d->ub[i] - d->lb[i];
 457      return f;
 458 }
 459
 460 nlopt_result cdirect(int n, nlopt_func f, void *f_data,
 461                      const double *lb, const double *ub,
 462                      double *x,
 463                      double *fmin,
 464                      nlopt_stopping *stop,
 465                      double magic_eps, int which_alg)
 466 {
 467      uf_data d;
 468      nlopt_result ret;
 469      const double *xtol_abs_save;
 470      int i;
 471
 472      d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
 473      d.x = (double *) malloc(sizeof(double) * n*4);
 474      if (!d.x) return NLOPT_OUT_OF_MEMORY;
 475
 476      for (i = 0; i < n; ++i) {
 477           x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
 478           d.x[n+i] = 0;
 479           d.x[2*n+i] = 1;
 480           d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
 481      }
 482      xtol_abs_save = stop->xtol_abs;
 483      stop->xtol_abs = d.x + 3*n;
 484      ret = cdirect_unscaled(n, uf, &d, d.x+n, d.x+2*n, x, fmin, stop,
 485                             magic_eps, which_alg);
 486      stop->xtol_abs = xtol_abs_save;
 487      for (i = 0; i < n; ++i)
 488           x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);
 489      free(d.x);
 490      return ret;
 491 }