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 "config.h"
   9
  10 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  11 #define MAX(a,b) ((a) > (b) ? (a) : (b))
  12
  13 /***************************************************************************/
  14 /* basic data structure:
  15  *
  16  * a hyper-rectangle is stored as an array of length 2n+2, where [0]
  17  * is the value of the function at the center, [1] is the "size"
  18  * measure of the rectangle, [2..n+1] are the coordinates of the
  19  * center, and [n+2..2n+1] are the widths of the sides.
  20  *
  21  * a list of rectangles is just an array of N hyper-rectangles
  22  * stored as an N x (2n+1) in row-major order.  Generally,
  23  * we allocate more than we need, allocating Na hyper-rectangles.
  24  *
  25  * n > 0 always
  26  */
  27
  28 #define RECT_LEN(n) (2*(n)+2) /* number of double values in a hyperrect */
  29
  30 /* parameters of the search algorithm and various information that
  31    needs to be passed around */
  32 typedef struct {
  33      double magic_eps; /* Jones' epsilon parameter (1e-4 is recommended) */
  34      int which_diam; /* which measure of hyper-rectangle diam to use:
  35                         0 = Jones, 1 = Gablonsky */
  36      int which_div; /* which way to divide rects:
  37                        0: Gablonsky (cubes divide all, rects longest)
  38                        1: orig. Jones (divide all longest sides)
  39                        2: Jones Encyc. Opt.: pick random longest side */
  40      const double *lb, *ub;
  41      nlopt_stopping *stop; /* stopping criteria */
  42      int n;
  43      nlopt_func f; void *f_data;
  44      double *work; /* workspace, of length >= 2*n */
  45      int *iwork, iwork_len; /* workspace, of length iwork_len >= n */
  46      double fmin, *xmin; /* minimum so far */
  47 } params;
  48
  49 /***************************************************************************/
  50
  51 /* evaluate the "diameter" (d) of a rectangle of widths w[n] */
  52 static double rect_diameter(int n, const double *w, const params *p)
  53 {
  54      int i;
  55      if (p->which_diam == 0) { /* Jones measure */
  56           double sum = 0;
  57           for (i = 0; i < n; ++i)
  58                sum += w[i] * w[i];
  59           return sqrt(sum) * 0.5; /* distance from center to a vertex */
  60      }
  61      else { /* Gablonsky measure */
  62           double maxw = 0;
  63           for (i = 0; i < n; ++i)
  64                if (w[i] > maxw)
  65                     maxw = w[i];
  66           return maxw * 0.5; /* half-width of longest side */
  67      }
  68 }
  69
  70 static double *alloc_rects(int n, int *Na, double *rects, int newN)
  71 {
  72      if (newN <= *Na)
  73           return rects;
  74      else {
  75           (*Na) += newN;
  76           return realloc(rects, sizeof(double) * RECT_LEN(n) * (*Na));
  77      }
  78 }
  79 #define ALLOC_RECTS(n, Nap, rects, newN) if (!(rects = alloc_rects(n, Nap, rects, newN))) return NLOPT_OUT_OF_MEMORY
  80
  81 static double *fv_qsort = 0;
  82 static int sort_fv_compare(const void *a_, const void *b_)
  83 {
  84      int a = *((const int *) a_), b = *((const int *) b_);
  85      double fa = MIN(fv_qsort[2*a], fv_qsort[2*a+1]);
  86      double fb = MIN(fv_qsort[2*b], fv_qsort[2*b+1]);
  87      if (fa < fb)
  88           return -1;
  89      else if (fa > fb)
  90           return +1;
  91      else
  92           return 0;
  93 }
  94 static void sort_fv(int n, double *fv, int *isort)
  95 {
  96      int i;
  97      for (i = 0; i < n; ++i) isort[i] = i;
  98      fv_qsort = fv; /* not re-entrant, sigh... */
  99      qsort(isort, (unsigned) n, sizeof(int), sort_fv_compare);
 100      fv_qsort = 0;
 101 }
 102
 103 static double function_eval(const double *x, params *p) {
 104      double f = p->f(p->n, x, NULL, p->f_data);
 105      if (f < p->fmin) {
 106           p->fmin = f;
 107           memcpy(p->xmin, x, sizeof(double) * p->n);
 108      }
 109      p->stop->nevals++;
 110      return f;
 111 }
 112 #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
 113
 114 #define THIRD (0.3333333333333333333333)
 115
 116 #define EQUAL_SIDE_TOL 5e-2 /* tolerance to equate side sizes */
 117
 118 /* divide rectangle idiv in the list rects */
 119 static nlopt_result divide_rect(int *N, int *Na, double **rects, int idiv,
 120                                 params *p)
 121 {
 122      int i;
 123      const const int n = p->n;
 124      const int L = RECT_LEN(n);
 125      double *r = *rects;
 126      double *c = r + L*idiv + 2; /* center of rect to divide */
 127      double *w = c + n; /* widths of rect to divide */
 128      double wmax = w[0];
 129      int imax = 0, nlongest = 0;
 130
 131      for (i = 1; i < n; ++i)
 132           if (w[i] > wmax)
 133                wmax = w[imax = i];
 134      for (i = 0; i < n; ++i)
 135           if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL)
 136                ++nlongest;
 137      if (p->which_div == 1 || (p->which_div == 0 && nlongest == n)) {
 138           /* trisect all longest sides, in increasing order of the average
 139              function value along that direction */
 140           double *fv = p->work;
 141           int *isort = p->iwork;
 142           for (i = 0; i < n; ++i) {
 143                if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL) {
 144                     double csave = c[i];
 145                     c[i] = csave - w[i] * THIRD;
 146                     FUNCTION_EVAL(fv[2*i], c, p);
 147                     c[i] = csave + w[i] * THIRD;
 148                     FUNCTION_EVAL(fv[2*i+1], c, p);
 149                     c[i] = csave;
 150                }
 151                else {
 152                     fv[2*i] = fv[2*i+1] = HUGE_VAL;
 153                }
 154           }
 155           sort_fv(n, fv, isort);
 156           ALLOC_RECTS(n, Na, r, (*N)+2*nlongest);
 157           *rects = r; c = r + L*idiv + 2; w = c + n;
 158           for (i = 0; i < nlongest; ++i) {
 159                int k;
 160                w[isort[i]] *= THIRD;
 161                r[L*idiv + 1] = rect_diameter(n, w, p);
 162                for (k = 0; k <= 1; ++k) {
 163                     memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
 164                     r[L*(*N) + 2 + isort[i]] += w[isort[i]] * (2*k-1);
 165                     r[L*(*N)] = fv[2*isort[i]+k];
 166                     ++(*N);
 167                }
 168           }
 169      }
 170      else {
 171           int k;
 172           if (nlongest > 1 && p->which_div == 2) {
 173                /* randomly choose longest side */
 174                i = nlopt_iurand(nlongest);
 175                for (k = 0; k < n; ++k)
 176                     if (wmax - w[k] <= wmax * EQUAL_SIDE_TOL) {
 177                          if (!i) { i = k; break; }
 178                          --i;
 179                     }
 180           }
 181           else
 182                i = imax; /* trisect longest side */
 183           ALLOC_RECTS(n, Na, r, (*N)+2);
 184           *rects = r; c = r + L*idiv + 2; w = c + n;
 185           w[i] *= THIRD;
 186           r[L*idiv + 1] = rect_diameter(n, w, p);
 187           for (k = 0; k <= 1; ++k) {
 188                memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
 189                r[L*(*N) + 2 + i] += w[i] * (2*k-1);
 190                FUNCTION_EVAL(r[L*(*N)], r + L*(*N) + 2, p);
 191                ++(*N);
 192           }
 193      }
 194      return NLOPT_SUCCESS;
 195 }
 196
 197 /***************************************************************************/
 198 /* O(N log N) convex hull algorithm, used later to find the potentially
 199    optimal points */
 200
 201 /* sort ihull by xy in lexographic order by x,y */
 202 static int s_qsort = 1; static double *xy_qsort = 0;
 203 static int sort_xy_compare(const void *a_, const void *b_)
 204 {
 205      int a = *((const int *) a_), b = *((const int *) b_);
 206      double xa = xy_qsort[a*s_qsort+1], xb = xy_qsort[b*s_qsort+1];
 207      if (xa < xb) return -1;
 208      else if (xb < xa) return +1;
 209      else {
 210           double ya = xy_qsort[a*s_qsort], yb = xy_qsort[b*s_qsort];
 211           if (ya < yb) return -1;
 212           else if (ya > yb) return +1;
 213           else return 0;
 214      }
 215 }
 216 static void sort_xy(int N, double *xy, int s, int *isort)
 217 {
 218      int i;
 219
 220      for (i = 0; i < N; ++i) isort[i] = i;
 221      s_qsort = s; xy_qsort = xy;
 222      qsort(isort, (unsigned) N, sizeof(int), sort_xy_compare);
 223      xy_qsort = 0;
 224 }
 225
 226 /* Find the lower convex hull of a set of points (xy[s*i+1], xy[s*i]), where
 227    0 <= i < N and s >= 2.
 228
 229    The return value is the number of points in the hull, with indices
 230    stored in ihull.  ihull and is should point to arrays of length >= N.
 231
 232    Note that we don't allow redundant points along the same line in the
 233    hull, similar to Gablonsky's version of DIRECT and differing from
 234    Jones'. */
 235 static int convex_hull(int N, double *xy, int s, int *ihull, int *is)
 236 {
 237      int minmin; /* first index (smallest y) with min x */
 238      int minmax; /* last index (largest y) with min x */
 239      int maxmin; /* first index (smallest y) with max x */
 240      int maxmax; /* last index (largest y) with max x */
 241      int i, nhull = 0;
 242      double minslope;
 243      double xmin, xmax;
 244
 245      /* Monotone chain algorithm [Andrew, 1979]. */
 246
 247      sort_xy(N, xy, s, is);
 248
 249      xmin = xy[s*is[minmin=0]+1]; xmax = xy[s*is[maxmax=N-1]+1];
 250
 251      if (xmin == xmax) { /* degenerate case */
 252           ihull[nhull++] = is[minmin];
 253           return nhull;
 254      }
 255
 256      for (minmax = minmin; minmax+1 < N && xy[s*is[minmax+1]+1]==xmin;
 257           ++minmax);
 258      for (maxmin = maxmax; maxmin-1>=0 && xy[s*is[maxmin-1]+1]==xmax;
 259           --maxmin);
 260
 261      minslope = (xy[s*is[maxmin]] - xy[s*is[minmin]]) / (xmax - xmin);
 262
 263      ihull[nhull++] = is[minmin];
 264      for (i = minmax + 1; i < maxmin; ++i) {
 265           int k = is[i];
 266           if (xy[s*k] > xy[s*is[minmin]] + (xy[s*k+1] - xmin) * minslope)
 267                continue;
 268           /* remove points until we are making a "left turn" to k */
 269           while (nhull > 1) {
 270                int t1 = ihull[nhull - 1], t2 = ihull[nhull - 2];
 271                /* cross product (t1-t2) x (k-t2) > 0 for a left turn: */
 272                if ((xy[s*t1+1]-xy[s*t2+1]) * (xy[s*k]-xy[s*t2])
 273                    - (xy[s*t1]-xy[s*t2]) * (xy[s*k+1]-xy[s*t2+1]) > 0)
 274                     break;
 275                --nhull;
 276           }
 277           ihull[nhull++] = k;
 278      }
 279      ihull[nhull++] = is[maxmin];
 280      return nhull;
 281 }
 282
 283 /***************************************************************************/
 284
 285 static int small(double *w, params *p)
 286 {
 287      int i;
 288      for (i = 0; i < p->n; ++i)
 289           if (w[i] > p->stop->xtol_abs[i] &&
 290               w[i] > (p->ub[i] - p->lb[i]) * p->stop->xtol_rel)
 291                return 0;
 292      return 1;
 293 }
 294
 295 static nlopt_result divide_good_rects(int *N, int *Na, double **rects,
 296                                       params *p)
 297 {
 298      const int n = p->n;
 299      const int L = RECT_LEN(n);
 300      int *ihull, nhull, i, xtol_reached = 1, divided_some = 0;
 301      double *r = *rects;
 302      double magic_eps = p->magic_eps;
 303
 304      if (p->iwork_len < n + 2*(*N)) {
 305           p->iwork_len = p->iwork_len + n + 2*(*N);
 306           p->iwork = (int *) realloc(p->iwork, sizeof(int) * p->iwork_len);
 307           if (!p->iwork)
 308                return NLOPT_OUT_OF_MEMORY;
 309      }
 310      ihull = p->iwork;
 311      nhull = convex_hull(*N, r, L, ihull, ihull + *N);
 312  divisions:
 313      for (i = 0; i < nhull; ++i) {
 314           double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
 315           if (i > 0)
 316                K1 = (r[L*ihull[i]] - r[L*ihull[i-1]]) /
 317                     (r[L*ihull[i]+1] - r[L*ihull[i-1]+1]);
 318           if (i < nhull-1)
 319                K1 = (r[L*ihull[i]] - r[L*ihull[i+1]]) /
 320                     (r[L*ihull[i]+1] - r[L*ihull[i+1]+1]);
 321           K = MAX(K1, K2);
 322           if (r[L*ihull[i]] - K * r[L*ihull[i]+1]
 323               <= p->fmin - magic_eps * fabs(p->fmin)) {
 324                /* "potentially optimal" rectangle, so subdivide */
 325                divided_some = 1;
 326                nlopt_result ret;
 327                ret = divide_rect(N, Na, rects, ihull[i], p);
 328                r = *rects; /* may have grown */
 329                if (ret != NLOPT_SUCCESS) return ret;
 330                xtol_reached = xtol_reached && small(r + L*ihull[i] + 2+n, p);
 331           }
 332      }
 333      if (!divided_some) {
 334           if (magic_eps != 0) {
 335                magic_eps = 0;
 336                goto divisions; /* try again */
 337           }
 338           else { /* WTF? divide largest rectangle */
 339                double wmax = r[1];
 340                int imax = 0;
 341                for (i = 1; i < *N; ++i)
 342                     if (r[L*i+1] > wmax)
 343                          wmax = r[L*(imax=i)+1];
 344                return divide_rect(N, Na, rects, imax, p);
 345           }
 346      }
 347      return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
 348 }
 349
 350 /***************************************************************************/
 351
 352 nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
 353                               const double *lb, const double *ub,
 354                               double *x,
 355                               double *fmin,
 356                               nlopt_stopping *stop,
 357                               double magic_eps, int which_alg)
 358 {
 359      params p;
 360      double *rects;
 361      int Na = 100, N = 1, i, x_center = 1;
 362      nlopt_result ret = NLOPT_OUT_OF_MEMORY;
 363
 364      p.magic_eps = magic_eps;
 365      p.which_diam = which_alg % 2;
 366      p.which_div = 0;
 367      p.lb = lb; p.ub = ub;
 368      p.stop = stop;
 369      p.n = n;
 370      p.f = f;
 371      p.f_data = f_data;
 372      p.xmin = x;
 373      p.fmin = f(n, x, NULL, f_data); stop->nevals++;
 374      p.work = 0;
 375      p.iwork = 0;
 376      rects = 0;
 377      p.work = (double *) malloc(sizeof(double) * 2*n);
 378      if (!p.work) goto done;
 379      rects = (double *) malloc(sizeof(double) * Na * RECT_LEN(n));
 380      if (!rects) goto done;
 381      p.iwork = (int *) malloc(sizeof(int) * (p.iwork_len = 2*Na + n));
 382      if (!p.iwork) goto done;
 383
 384      for (i = 0; i < n; ++i) {
 385           rects[2+i] = 0.5 * (lb[i] + ub[i]);
 386           x_center = x_center
 387                && (fabs(rects[2+i]-x[i]) < 1e-13*(1+fabs(x[i])));
 388           rects[2+n+i] = ub[i] - lb[i];
 389      }
 390      rects[1] = rect_diameter(n, rects+2+n, &p);
 391      if (x_center)
 392           rects[0] = p.fmin; /* avoid computing f(center) twice */
 393      else
 394           rects[0] = function_eval(rects+2, &p);
 395
 396      ret = divide_rect(&N, &Na, &rects, 0, &p);
 397      if (ret != NLOPT_SUCCESS) goto done;
 398
 399      while (1) {
 400           double fmin0 = p.fmin;
 401           ret = divide_good_rects(&N, &Na, &rects, &p);
 402           if (ret != NLOPT_SUCCESS) goto done;
 403           if (nlopt_stop_f(p.stop, p.fmin, fmin0)) {
 404                ret = NLOPT_FTOL_REACHED;
 405                goto done;
 406           }
 407      }
 408
 409  done:
 410      free(p.iwork);
 411      free(rects);
 412      free(p.work);
 413
 414      *fmin = p.fmin;
 415      return ret;
 416 }
 417
 418 /* in the conventional DIRECT-type algorithm, we first rescale our
 419    coordinates to a unit hypercube ... we do this simply by
 420    wrapping cdirect() around cdirect_unscaled(). */
 421
 422 typedef struct {
 423      nlopt_func f;
 424      void *f_data;
 425      double *x;
 426      const double *lb, *ub;
 427 } uf_data;
 428 static double uf(int n, const double *xu, double *grad, void *d_)
 429 {
 430      uf_data *d = (uf_data *) d_;
 431      double f;
 432      int i;
 433      for (i = 0; i < n; ++i)
 434           d->x[i] = d->lb[i] + xu[i] * (d->ub[i] - d->lb[i]);
 435      f = d->f(n, d->x, grad, d->f_data);
 436      if (grad)
 437           for (i = 0; i < n; ++i)
 438                grad[i] *= d->ub[i] - d->lb[i];
 439      return f;
 440 }
 441
 442 nlopt_result cdirect(int n, nlopt_func f, void *f_data,
 443                      const double *lb, const double *ub,
 444                      double *x,
 445                      double *fmin,
 446                      nlopt_stopping *stop,
 447                      double magic_eps, int which_alg)
 448 {
 449      uf_data d;
 450      nlopt_result ret;
 451      const double *xtol_abs_save;
 452      int i;
 453
 454      d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
 455      d.x = (double *) malloc(sizeof(double) * n*4);
 456      if (!d.x) return NLOPT_OUT_OF_MEMORY;
 457
 458      for (i = 0; i < n; ++i) {
 459           x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
 460           d.x[n+i] = 0;
 461           d.x[2*n+i] = 1;
 462           d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
 463      }
 464      xtol_abs_save = stop->xtol_abs;
 465      stop->xtol_abs = d.x + 3*n;
 466      ret = cdirect_unscaled(n, uf, &d, d.x+n, d.x+2*n, x, fmin, stop,
 467                             magic_eps, which_alg);
 468      stop->xtol_abs = xtol_abs_save;
 469      for (i = 0; i < n; ++i)
 470           x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);
 471      free(d.x);
 472      return ret;
 473 }