cdirect/cdirect.c

   1 /* Copyright (c) 2007-2014 Massachusetts Institute of Technology
   2  *
   3  * Permission is hereby granted, free of charge, to any person obtaining
   4  * a copy of this software and associated documentation files (the
   5  * "Software"), to deal in the Software without restriction, including
   6  * without limitation the rights to use, copy, modify, merge, publish,
   7  * distribute, sublicense, and/or sell copies of the Software, and to
   8  * permit persons to whom the Software is furnished to do so, subject to
   9  * the following conditions:
  10  *
  11  * The above copyright notice and this permission notice shall be
  12  * included in all copies or substantial portions of the Software.
  13  *
  14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  16  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  17  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  18  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  19  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  20  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  21  */
  22
  23 #include <math.h>
  24 #include <stdlib.h>
  25 #include <string.h>
  26
  27 #include "nlopt-util.h"
  28 #include "nlopt.h"
  29 #include "cdirect.h"
  30 #include "redblack.h"
  31
  32 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  33 #define MAX(a,b) ((a) > (b) ? (a) : (b))
  34
  35 /***************************************************************************/
  36 /* basic data structure:
  37  *
  38  * a hyper-rectangle is stored as an array of length L = 2n+3, where [1]
  39  * is the value (f) of the function at the center, [0] is the "size"
  40  * measure (d) of the rectangle, [3..n+2] are the coordinates of the
  41  * center (c), [n+3..2n+2] are the widths of the sides (w), and [2]
  42  * is an "age" measure for tie-breaking purposes.
  43  *
  44  * we store the hyper-rectangles in a red-black tree, sorted by (d,f)
  45  * in lexographic order, to allow us to perform quick convex-hull
  46  * calculations (in the future, we might make this data structure
  47  * more sophisticated based on the dynamic convex-hull literature).
  48  *
  49  * n > 0 always, of course.
  50  */
  51
  52 /* parameters of the search algorithm and various information that
  53    needs to be passed around */
  54 typedef struct {
  55      int n; /* dimension */
  56      int L; /* size of each rectangle (2n+3) */
  57      double magic_eps; /* Jones' epsilon parameter (1e-4 is recommended) */
  58      int which_diam; /* which measure of hyper-rectangle diam to use:
  59                         0 = Jones, 1 = Gablonsky */
  60      int which_div; /* which way to divide rects:
  61                        0: orig. Jones (divide all longest sides)
  62                        1: Gablonsky (cubes divide all, rects longest)
  63                        2: Jones Encyc. Opt.: pick random longest side */
  64      int which_opt; /* which rects are considered "potentially optimal"
  65                        0: Jones (all pts on cvx hull, even equal pts)
  66                        1: Gablonsky DIRECT-L (pick one pt, if equal pts)
  67                        2: ~ 1, but pick points randomly if equal pts
  68                     ... 2 seems to suck compared to just picking oldest pt */
  69
  70      const double *lb, *ub;
  71      nlopt_stopping *stop; /* stopping criteria */
  72      nlopt_func f; void *f_data;
  73      double *work; /* workspace, of length >= 2*n */
  74      int *iwork; /* workspace, length >= n */
  75      double minf, *xmin; /* minimum so far */
  76
  77      /* red-black tree of hyperrects, sorted by (d,f,age) in
  78         lexographical order */
  79      rb_tree rtree;
  80      int age; /* age for next new rect */
  81      double **hull; /* array to store convex hull */
  82      int hull_len; /* allocated length of hull array */
  83 } params;
  84
  85 /***************************************************************************/
  86
  87 /* Evaluate the "diameter" (d) of a rectangle of widths w[n]
  88
  89    We round the result to single precision, which should be plenty for
  90    the use we put the diameter to (rect sorting), to allow our
  91    performance hack in convex_hull to work (in the Jones and Gablonsky
  92    DIRECT algorithms, all of the rects fall into a few diameter
  93    values, and we don't want rounding error to spoil this) */
  94 static double rect_diameter(int n, const double *w, const params *p)
  95 {
  96      int i;
  97      if (p->which_diam == 0) { /* Jones measure */
  98           double sum = 0;
  99           for (i = 0; i < n; ++i)
 100                sum += w[i] * w[i];
 101           /* distance from center to a vertex */
 102           return ((float) (sqrt(sum) * 0.5));
 103      }
 104      else { /* Gablonsky measure */
 105           double maxw = 0;
 106           for (i = 0; i < n; ++i)
 107                if (w[i] > maxw)
 108                     maxw = w[i];
 109           /* half-width of longest side */
 110           return ((float) (maxw * 0.5));
 111      }
 112 }
 113
 114 #define ALLOC_RECT(rect, L) if (!(rect = (double*) malloc(sizeof(double)*(L)))) return NLOPT_OUT_OF_MEMORY
 115
 116 static int sort_fv_compare(void *fv_, const void *a_, const void *b_)
 117 {
 118      const double *fv = (const double *) fv_;
 119      int a = *((const int *) a_), b = *((const int *) b_);
 120      double fa = MIN(fv[2*a], fv[2*a+1]);
 121      double fb = MIN(fv[2*b], fv[2*b+1]);
 122      if (fa < fb)
 123           return -1;
 124      else if (fa > fb)
 125           return +1;
 126      else
 127           return 0;
 128 }
 129 static void sort_fv(int n, double *fv, int *isort)
 130 {
 131      int i;
 132      for (i = 0; i < n; ++i) isort[i] = i;
 133      nlopt_qsort_r(isort, (unsigned) n, sizeof(int), fv, sort_fv_compare);
 134 }
 135
 136 static double function_eval(const double *x, params *p) {
 137      double f = p->f(p->n, x, NULL, p->f_data);
 138      if (f < p->minf) {
 139           p->minf = f;
 140           memcpy(p->xmin, x, sizeof(double) * p->n);
 141      }
 142      p->stop->nevals++;
 143      return f;
 144 }
 145 #define FUNCTION_EVAL(fv,x,p,freeonerr) fv = function_eval(x, p); if (nlopt_stop_forced((p)->stop)) { free(freeonerr); return NLOPT_FORCED_STOP; } else if (p->minf < p->stop->minf_max) { free(freeonerr); return NLOPT_MINF_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; }
 146
 147 #define THIRD (0.3333333333333333333333)
 148
 149 #define EQUAL_SIDE_TOL 5e-2 /* tolerance to equate side sizes */
 150
 151 /* divide rectangle idiv in the list p->rects */
 152 static nlopt_result divide_rect(double *rdiv, params *p)
 153 {
 154      int i;
 155      const int n = p->n;
 156      const int L = p->L;
 157      double *c = rdiv + 3; /* center of rect to divide */
 158      double *w = c + n; /* widths of rect to divide */
 159      double wmax = w[0];
 160      int imax = 0, nlongest = 0;
 161      rb_node *node;
 162
 163      for (i = 1; i < n; ++i)
 164           if (w[i] > wmax)
 165                wmax = w[imax = i];
 166      for (i = 0; i < n; ++i)
 167           if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL)
 168                ++nlongest;
 169      if (p->which_div == 1 || (p->which_div == 0 && nlongest == n)) {
 170           /* trisect all longest sides, in increasing order of the average
 171              function value along that direction */
 172           double *fv = p->work;
 173           int *isort = p->iwork;
 174           for (i = 0; i < n; ++i) {
 175                if (wmax - w[i] <= wmax * EQUAL_SIDE_TOL) {
 176                     double csave = c[i];
 177                     c[i] = csave - w[i] * THIRD;
 178                     FUNCTION_EVAL(fv[2*i], c, p, 0);
 179                     c[i] = csave + w[i] * THIRD;
 180                     FUNCTION_EVAL(fv[2*i+1], c, p, 0);
 181                     c[i] = csave;
 182                }
 183                else {
 184                     fv[2*i] = fv[2*i+1] = HUGE_VAL;
 185                }
 186           }
 187           sort_fv(n, fv, isort);
 188           if (!(node = rb_tree_find(&p->rtree, rdiv)))
 189                return NLOPT_FAILURE;
 190           for (i = 0; i < nlongest; ++i) {
 191                int k;
 192                w[isort[i]] *= THIRD;
 193                rdiv[0] = rect_diameter(n, w, p);
 194                rdiv[2] = p->age++;
 195                node = rb_tree_resort(&p->rtree, node);
 196                for (k = 0; k <= 1; ++k) {
 197                     double *rnew;
 198                     ALLOC_RECT(rnew, L);
 199                     memcpy(rnew, rdiv, sizeof(double) * L);
 200                     rnew[3 + isort[i]] += w[isort[i]] * (2*k-1);
 201                     rnew[1] = fv[2*isort[i]+k];
 202                     rnew[2] = p->age++;
 203                     if (!rb_tree_insert(&p->rtree, rnew)) {
 204                          free(rnew);
 205                          return NLOPT_OUT_OF_MEMORY;
 206                     }
 207                }
 208           }
 209      }
 210      else {
 211           int k;
 212           if (nlongest > 1 && p->which_div == 2) {
 213                /* randomly choose longest side */
 214                i = nlopt_iurand(nlongest);
 215                for (k = 0; k < n; ++k)
 216                     if (wmax - w[k] <= wmax * EQUAL_SIDE_TOL) {
 217                          if (!i) { i = k; break; }
 218                          --i;
 219                     }
 220           }
 221           else
 222                i = imax; /* trisect longest side */
 223           if (!(node = rb_tree_find(&p->rtree, rdiv)))
 224                return NLOPT_FAILURE;
 225           w[i] *= THIRD;
 226           rdiv[0] = rect_diameter(n, w, p);
 227           rdiv[2] = p->age++;
 228           node = rb_tree_resort(&p->rtree, node);
 229           for (k = 0; k <= 1; ++k) {
 230                double *rnew;
 231                ALLOC_RECT(rnew, L);
 232                memcpy(rnew, rdiv, sizeof(double) * L);
 233                rnew[3 + i] += w[i] * (2*k-1);
 234                FUNCTION_EVAL(rnew[1], rnew + 3, p, rnew);
 235                rnew[2] = p->age++;
 236                if (!rb_tree_insert(&p->rtree, rnew)) {
 237                     free(rnew);
 238                     return NLOPT_OUT_OF_MEMORY;
 239                }
 240           }
 241      }
 242      return NLOPT_SUCCESS;
 243 }
 244
 245 /***************************************************************************/
 246 /* Convex hull algorithm, used later to find the potentially optimal
 247    points.  What we really have in DIRECT is a "dynamic convex hull"
 248    problem, since we are dynamically adding/removing points and
 249    updating the hull, but I haven't implemented any of the fancy
 250    algorithms for this problem yet. */
 251
 252 /* Find the lower convex hull of a set of points (x,y) stored in a rb-tree
 253    of pointers to {x,y} arrays sorted in lexographic order by (x,y).
 254
 255    Unlike standard convex hulls, we allow redundant points on the hull,
 256    and even allow duplicate points if allow_dups is nonzero.
 257
 258    The return value is the number of points in the hull, with pointers
 259    stored in hull[i] (should be an array of length >= t->N).
 260 */
 261 static int convex_hull(rb_tree *t, double **hull, int allow_dups)
 262 {
 263      int nhull = 0;
 264      double minslope;
 265      double xmin, xmax, yminmin, ymaxmin;
 266      rb_node *n, *nmax;
 267
 268      /* Monotone chain algorithm [Andrew, 1979]. */
 269
 270      n = rb_tree_min(t);
 271      if (!n) return 0;
 272      nmax = rb_tree_max(t);
 273
 274      xmin = n->k[0];
 275      yminmin = n->k[1];
 276      xmax = nmax->k[0];
 277
 278      if (allow_dups)
 279           do { /* include any duplicate points at (xmin,yminmin) */
 280                hull[nhull++] = n->k;
 281                n = rb_tree_succ(n);
 282           } while (n && n->k[0] == xmin && n->k[1] == yminmin);
 283      else
 284           hull[nhull++] = n->k;
 285
 286      if (xmin == xmax) return nhull;
 287
 288      /* set nmax = min mode with x == xmax */
 289 #if 0
 290      while (nmax->k[0] == xmax)
 291           nmax = rb_tree_pred(nmax); /* non-NULL since xmin != xmax */
 292      nmax = rb_tree_succ(nmax);
 293 #else
 294      /* performance hack (see also below) */
 295      {
 296           double kshift[2];
 297           kshift[0] = xmax * (1 - 1e-13);
 298           kshift[1] = -HUGE_VAL;
 299           nmax = rb_tree_find_gt(t, kshift); /* non-NULL since xmin != xmax */
 300      }
 301 #endif
 302
 303      ymaxmin = nmax->k[1];
 304      minslope = (ymaxmin - yminmin) / (xmax - xmin);
 305
 306      /* set n = first node with x != xmin */
 307 #if 0
 308      while (n->k[0] == xmin)
 309           n = rb_tree_succ(n); /* non-NULL since xmin != xmax */
 310 #else
 311      /* performance hack (see also below) */
 312      {
 313           double kshift[2];
 314           kshift[0] = xmin * (1 + 1e-13);
 315           kshift[1] = -HUGE_VAL;
 316           n = rb_tree_find_gt(t, kshift); /* non-NULL since xmin != xmax */
 317      }
 318 #endif
 319
 320      for (; n != nmax; n = rb_tree_succ(n)) {
 321           double *k = n->k;
 322           if (k[1] > yminmin + (k[0] - xmin) * minslope)
 323                continue;
 324
 325           /* performance hack: most of the points in DIRECT lie along
 326              vertical lines at a few x values, and we can exploit this */
 327           if (nhull && k[0] == hull[nhull - 1][0]) { /* x == previous x */
 328                if (k[1] > hull[nhull - 1][1]) {
 329                     double kshift[2];
 330                     /* because of the round to float in rect_diameter, above,
 331                        it shouldn't be possible for two diameters (x values)
 332                        to have a fractional difference < 1e-13.  Note
 333                        that k[0] > 0 always in DIRECT */
 334                     kshift[0] = k[0] * (1 + 1e-13);
 335                     kshift[1] = -HUGE_VAL;
 336                     n = rb_tree_pred(rb_tree_find_gt(t, kshift));
 337                     continue;
 338                }
 339                else { /* equal y values, add to hull */
 340                     if (allow_dups)
 341                          hull[nhull++] = k;
 342                     continue;
 343                }
 344           }
 345
 346           /* remove points until we are making a "left turn" to k */
 347           while (nhull > 1) {
 348                double *t1 = hull[nhull - 1], *t2;
 349
 350                /* because we allow equal points in our hull, we have
 351                   to modify the standard convex-hull algorithm slightly:
 352                   we need to look backwards in the hull list until we
 353                   find a point t2 != t1 */
 354                int it2 = nhull - 2;
 355                do {
 356                     t2 = hull[it2--];
 357                } while (it2 >= 0 && t2[0] == t1[0] && t2[1] == t1[1]);
 358                if (it2 < 0) break;
 359
 360                /* cross product (t1-t2) x (k-t2) > 0 for a left turn: */
 361                if ((t1[0]-t2[0]) * (k[1]-t2[1])
 362                    - (t1[1]-t2[1]) * (k[0]-t2[0]) >= 0)
 363                     break;
 364                --nhull;
 365           }
 366           hull[nhull++] = k;
 367      }
 368
 369      if (allow_dups)
 370           do { /* include any duplicate points at (xmax,ymaxmin) */
 371                hull[nhull++] = nmax->k;
 372                nmax = rb_tree_succ(nmax);
 373           } while (nmax && nmax->k[0] == xmax && nmax->k[1] == ymaxmin);
 374      else
 375           hull[nhull++] = nmax->k;
 376
 377      return nhull;
 378 }
 379
 380 /***************************************************************************/
 381
 382 static int small(double *w, params *p)
 383 {
 384      int i;
 385      for (i = 0; i < p->n; ++i)
 386           if (w[i] > p->stop->xtol_abs[i] &&
 387               w[i] > (p->ub[i] - p->lb[i]) * p->stop->xtol_rel)
 388                return 0;
 389      return 1;
 390 }
 391
 392 static nlopt_result divide_good_rects(params *p)
 393 {
 394      const int n = p->n;
 395      double **hull;
 396      int nhull, i, xtol_reached = 1, divided_some = 0;
 397      double magic_eps = p->magic_eps;
 398
 399      if (p->hull_len < p->rtree.N) {
 400           p->hull_len += p->rtree.N;
 401           p->hull = (double **) realloc(p->hull, sizeof(double*)*p->hull_len);
 402           if (!p->hull) return NLOPT_OUT_OF_MEMORY;
 403      }
 404      nhull = convex_hull(&p->rtree, hull = p->hull, p->which_opt != 1);
 405  divisions:
 406      for (i = 0; i < nhull; ++i) {
 407           double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
 408           int im, ip;
 409
 410           /* find unequal points before (im) and after (ip) to get slope */
 411           for (im = i-1; im >= 0 && hull[im][0] == hull[i][0]; --im) ;
 412           for (ip = i+1; ip < nhull && hull[ip][0] == hull[i][0]; ++ip) ;
 413
 414           if (im >= 0)
 415                K1 = (hull[i][1] - hull[im][1]) / (hull[i][0] - hull[im][0]);
 416           if (ip < nhull)
 417                K2 = (hull[i][1] - hull[ip][1]) / (hull[i][0] - hull[ip][0]);
 418           K = MAX(K1, K2);
 419           if (hull[i][1] - K * hull[i][0]
 420               <= p->minf - magic_eps * fabs(p->minf) || ip == nhull) {
 421                /* "potentially optimal" rectangle, so subdivide */
 422                nlopt_result ret = divide_rect(hull[i], p);
 423                divided_some = 1;
 424                if (ret != NLOPT_SUCCESS) return ret;
 425                xtol_reached = xtol_reached && small(hull[i] + 3+n, p);
 426           }
 427
 428           /* for the DIRECT-L variant, we only divide one rectangle out
 429              of all points with equal diameter and function values
 430              ... note that for p->which_opt == 1, i == ip-1 should be a no-op
 431                  anyway, since we set allow_dups=0 in convex_hull above */
 432           if (p->which_opt == 1)
 433                i = ip - 1; /* skip to next unequal point for next iteration */
 434           else if (p->which_opt == 2) /* like DIRECT-L but randomized */
 435                i += nlopt_iurand(ip - i); /* possibly do another equal pt */
 436      }
 437      if (!divided_some) {
 438           if (magic_eps != 0) {
 439                magic_eps = 0;
 440                goto divisions; /* try again */
 441           }
 442           else { /* WTF? divide largest rectangle with smallest f */
 443                /* (note that this code actually gets called from time
 444                   to time, and the heuristic here seems to work well,
 445                   but I don't recall this situation being discussed in
 446                   the references?) */
 447                rb_node *max = rb_tree_max(&p->rtree);
 448                rb_node *pred = max;
 449                double wmax = max->k[0];
 450                do { /* note: this loop is O(N) worst-case time */
 451                     max = pred;
 452                     pred = rb_tree_pred(max);
 453                } while (pred && pred->k[0] == wmax);
 454                return divide_rect(max->k, p);
 455           }
 456      }
 457      return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
 458 }
 459
 460 /***************************************************************************/
 461
 462 /* lexographic sort order (d,f,age) of hyper-rects, for red-black tree */
 463 int cdirect_hyperrect_compare(double *a, double *b)
 464 {
 465      if (a[0] < b[0]) return -1;
 466      if (a[0] > b[0]) return +1;
 467      if (a[1] < b[1]) return -1;
 468      if (a[1] > b[1]) return +1;
 469      if (a[2] < b[2]) return -1;
 470      if (a[2] > b[2]) return +1;
 471      return (int) (a - b); /* tie-breaker, shouldn't be needed */
 472 }
 473
 474 /***************************************************************************/
 475
 476 nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
 477                               const double *lb, const double *ub,
 478                               double *x,
 479                               double *minf,
 480                               nlopt_stopping *stop,
 481                               double magic_eps, int which_alg)
 482 {
 483      params p;
 484      int i;
 485      double *rnew;
 486      nlopt_result ret = NLOPT_OUT_OF_MEMORY;
 487
 488      p.magic_eps = magic_eps;
 489      p.which_diam = which_alg % 3;
 490      p.which_div = (which_alg / 3) % 3;
 491      p.which_opt = (which_alg / (3*3)) % 3;
 492      p.lb = lb; p.ub = ub;
 493      p.stop = stop;
 494      p.n = n;
 495      p.L = 2*n+3;
 496      p.f = f;
 497      p.f_data = f_data;
 498      p.xmin = x;
 499      p.minf = HUGE_VAL;
 500      p.work = 0;
 501      p.iwork = 0;
 502      p.hull = 0;
 503      p.age = 0;
 504
 505      rb_tree_init(&p.rtree, cdirect_hyperrect_compare);
 506
 507      p.work = (double *) malloc(sizeof(double) * (2*n));
 508      if (!p.work) goto done;
 509      p.iwork = (int *) malloc(sizeof(int) * n);
 510      if (!p.iwork) goto done;
 511      p.hull_len = 128; /* start with a reasonable number */
 512      p.hull = (double **) malloc(sizeof(double *) * p.hull_len);
 513      if (!p.hull) goto done;
 514
 515      if (!(rnew = (double *) malloc(sizeof(double) * p.L))) goto done;
 516      for (i = 0; i < n; ++i) {
 517           rnew[3+i] = 0.5 * (lb[i] + ub[i]);
 518           rnew[3+n+i] = ub[i] - lb[i];
 519      }
 520      rnew[0] = rect_diameter(n, rnew+3+n, &p);
 521      rnew[1] = function_eval(rnew+3, &p);
 522      rnew[2] = p.age++;
 523      if (!rb_tree_insert(&p.rtree, rnew)) {
 524           free(rnew);
 525           goto done;
 526      }
 527
 528      ret = divide_rect(rnew, &p);
 529      if (ret != NLOPT_SUCCESS) goto done;
 530
 531      while (1) {
 532           double minf0 = p.minf;
 533           ret = divide_good_rects(&p);
 534           if (ret != NLOPT_SUCCESS) goto done;
 535           if (p.minf < minf0 && nlopt_stop_f(p.stop, p.minf, minf0)) {
 536                ret = NLOPT_FTOL_REACHED;
 537                goto done;
 538           }
 539      }
 540
 541  done:
 542      rb_tree_destroy_with_keys(&p.rtree);
 543      free(p.hull);
 544      free(p.iwork);
 545      free(p.work);
 546
 547      *minf = p.minf;
 548      return ret;
 549 }
 550
 551 /* in the conventional DIRECT-type algorithm, we first rescale our
 552    coordinates to a unit hypercube ... we do this simply by
 553    wrapping cdirect() around cdirect_unscaled(). */
 554
 555 double cdirect_uf(unsigned n, const double *xu, double *grad, void *d_)
 556 {
 557      cdirect_uf_data *d = (cdirect_uf_data *) d_;
 558      double f;
 559      unsigned i;
 560      for (i = 0; i < n; ++i)
 561           d->x[i] = d->lb[i] + xu[i] * (d->ub[i] - d->lb[i]);
 562      f = d->f(n, d->x, grad, d->f_data);
 563      if (grad)
 564           for (i = 0; i < n; ++i)
 565                grad[i] *= d->ub[i] - d->lb[i];
 566      return f;
 567 }
 568
 569 nlopt_result cdirect(int n, nlopt_func f, void *f_data,
 570                      const double *lb, const double *ub,
 571                      double *x,
 572                      double *minf,
 573                      nlopt_stopping *stop,
 574                      double magic_eps, int which_alg)
 575 {
 576      cdirect_uf_data d;
 577      nlopt_result ret;
 578      const double *xtol_abs_save;
 579      int i;
 580
 581      d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
 582      d.x = (double *) malloc(sizeof(double) * n*4);
 583      if (!d.x) return NLOPT_OUT_OF_MEMORY;
 584
 585      for (i = 0; i < n; ++i) {
 586           x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
 587           d.x[n+i] = 0;
 588           d.x[2*n+i] = 1;
 589           d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
 590      }
 591      xtol_abs_save = stop->xtol_abs;
 592      stop->xtol_abs = d.x + 3*n;
 593      ret = cdirect_unscaled(n, cdirect_uf, &d, d.x+n, d.x+2*n, x, minf, stop,
 594                             magic_eps, which_alg);
 595      stop->xtol_abs = xtol_abs_save;
 596      for (i = 0; i < n; ++i)
 597           x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);
 598      free(d.x);
 599      return ret;
 600 }