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