cdirect/hybrid.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 /* Hybrid algorithm, inspired by DIRECT, that uses another local
  12  * optimization algorithm within each rectangle, and then looks
  13  * in the largest remaining rectangle (breaking ties by minimum
  14  * function value and then by age.
  15  *
  16  * Each hyperrect is represented by an array of length 3*n+3 consisting
  17  * of (d, -f, -a, x, c, w), where d=diameter, f=f(x), a=age, x=local optimum
  18  * c=center, w=widths.
  19  */
  20
  21 typedef struct {
  22      int n; /* dimension */
  23      int L; /* 3*n+3 */
  24      const double *lb, *ub;
  25      nlopt_stopping *stop; /* stopping criteria */
  26      nlopt_func f; void *f_data;
  27      double minf, *xmin; /* min so far */
  28      rb_tree rtree; /* red-black tree of rects, sorted by (d,-f,-a) */
  29      int age; /* age for next new rect */
  30      double *work; /* workspace of length >= 2*n */
  31
  32      nlopt_algorithm local_alg; /* local search algorithm */
  33      int local_maxeval; /* max # local iterations (0 if unlimited) */
  34
  35      int randomized_div; /* 1 to use randomized division algorithm */
  36 } params;
  37
  38 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  39
  40 #define THIRD (0.3333333333333333333333) /* 1/3 */
  41
  42 /************************************************************************/
  43
  44 static double fcount(int n, const double *x, double *grad, void *p_)
  45 {
  46      params *p = (params *) p_;
  47      p->stop->nevals++;
  48      return p->f(n, x, grad, p->f_data);
  49 }
  50
  51 static nlopt_result optimize_rect(double *r, params *p)
  52 {
  53      int i, n = p->n;
  54      double *lb = p->work, *ub = lb + n;
  55      double *x = r + 3, *c = x + n, *w = c + n;
  56      double t = nlopt_seconds();
  57      double minf;
  58      nlopt_stopping *stop = p->stop;
  59      nlopt_result ret;
  60
  61      if (stop->maxeval > 0 &&
  62          stop->nevals >= stop->maxeval) return NLOPT_MAXEVAL_REACHED;
  63      if (stop->maxtime > 0 &&
  64          t - stop->start >= stop->maxtime) return NLOPT_MAXTIME_REACHED;
  65
  66      for (i = 0; i < n; ++i) {
  67           lb[i] = c[i] - 0.5 * w[i];
  68           ub[i] = c[i] + 0.5 * w[i];
  69      }
  70      ret = nlopt_minimize(p->local_alg, n, fcount, p,
  71                           lb, ub, x, &minf,
  72                           stop->minf_max, stop->ftol_rel, stop->ftol_abs,
  73                           stop->xtol_rel, stop->xtol_abs,
  74                           p->local_maxeval > 0 ?
  75                           MIN(p->local_maxeval,
  76                               stop->maxeval - stop->nevals)
  77                           : stop->maxeval - stop->nevals,
  78                           stop->maxtime - (t - stop->start));
  79      r[1] = -minf;
  80      if (ret > 0) {
  81           if (minf < p->minf) {
  82                p->minf = minf;
  83                memcpy(p->xmin, x, sizeof(double) * n);
  84                if (ret == NLOPT_MINF_MAX_REACHED) return ret;
  85           }
  86           return NLOPT_SUCCESS;
  87      }
  88      return ret;
  89 }
  90
  91 /* given a hyperrect r, randomize the starting guess within the middle
  92    third of the box (don't guess too close to edges) */
  93 static void randomize_x(int n, double *r)
  94 {
  95      double *x = r + 3, *c = x + n, *w = c + n;
  96      int i;
  97      for (i = 0; i < n; ++i)
  98           x[i] = nlopt_urand(c[i] - w[i]*(0.5*THIRD),
  99                              c[i] + w[i]*(0.5*THIRD));
 100 }
 101
 102 /************************************************************************/
 103
 104 static double longest(int n, const double *w)
 105 {
 106      double wmax = w[n-1];
 107      for (n = n-2; n >= 0; n--) if (w[n] > wmax) wmax = w[n];
 108      return wmax;
 109 }
 110
 111 #define EQUAL_SIDE_TOL 5e-2 /* tolerance to equate side sizes */
 112
 113 static nlopt_result divide_largest(params *p)
 114 {
 115      int L = p->L;
 116      int n = p->n;
 117      rb_node *node = rb_tree_max(&p->rtree); /* just using it as a heap */
 118      double minf_start = p->minf;
 119      double *r = node->k, *rnew = NULL;
 120      double *x = r + 3, *c = x + n, *w = c + n;
 121      const double *lb = p->lb, *ub = p->ub;
 122      int i, idiv;
 123      double wmax;
 124      nlopt_result ret;
 125
 126      /* printf("rect:, %d, %g, %g, %g, %g\n", p->stop->nevals, c[0], c[1], w[0], w[1]); */
 127
 128      /* check xtol */
 129      for (i = 0; i < n; ++i)
 130           if (w[i] > p->stop->xtol_rel * (ub[i] - lb[i])
 131               && w[i] > p->stop->xtol_abs[i])
 132                break;
 133      if (i == n) return NLOPT_XTOL_REACHED;
 134
 135      if (p->randomized_div) { /* randomly pick among ~largest sides */
 136           int nlongest = 0;
 137           wmax = longest(n, w);
 138           for (i = 0; i < n; ++i)
 139                if (wmax - w[i] < EQUAL_SIDE_TOL * wmax) ++nlongest;
 140           i = 1 + nlopt_iurand(nlongest);
 141           for (idiv = 0; idiv < n; ++idiv) {
 142                if (wmax - w[idiv] < EQUAL_SIDE_TOL * wmax) --i;
 143                if (!i) break;
 144           }
 145      }
 146      else { /* just pick first largest side */
 147           wmax = w[idiv = 0];
 148           for (i = 1; i < n; ++i) if (w[i] > wmax) wmax = w[idiv = i];
 149      }
 150
 151      if (fabs(x[idiv] - c[idiv]) > (0.5 * THIRD) * w[idiv]) { /* bisect */
 152           double deltac = (x[idiv] > c[idiv] ? 0.25 : -0.25) * w[idiv];
 153           w[idiv] *= 0.5;
 154           c[idiv] += deltac;
 155           r[0] = longest(n, w); /* new diameter */
 156           /* r[1] unchanged since still contains local optimum x */
 157           r[2] = p->age--;
 158           node = rb_tree_resort(&p->rtree, node);
 159
 160           rnew = (double *) malloc(sizeof(double) * L);
 161           if (!rnew) return NLOPT_OUT_OF_MEMORY;
 162           memcpy(rnew, r, sizeof(double) * L);
 163           rnew[2] = p->age--;
 164           rnew[3+n+idiv] -= deltac*2;
 165           if (p->randomized_div)
 166                randomize_x(n, rnew);
 167           else
 168                memcpy(rnew+3, rnew+3+n, sizeof(double) * n); /* x = c */
 169           ret = optimize_rect(rnew, p);
 170           if (ret != NLOPT_SUCCESS) { free(rnew); return ret; }
 171           if (!rb_tree_insert(&p->rtree, rnew)) {
 172                free(rnew); return NLOPT_OUT_OF_MEMORY;
 173           }
 174      }
 175      else { /* trisect */
 176           w[idiv] *= THIRD;
 177           r[0] = longest(n, w);
 178           /* r[1] unchanged since still contains local optimum x */
 179           r[2] = p->age--;
 180           node = rb_tree_resort(&p->rtree, node);
 181
 182           for (i = -1; i <= +1; i += 2) {
 183                rnew = (double *) malloc(sizeof(double) * L);
 184                if (!rnew) return NLOPT_OUT_OF_MEMORY;
 185                memcpy(rnew, r, sizeof(double) * L);
 186                rnew[2] = p->age--;
 187                rnew[3+n+idiv] += w[i] * i;
 188                if (p->randomized_div)
 189                     randomize_x(n, rnew);
 190                else
 191                     memcpy(rnew+3, rnew+3+n, sizeof(double) * n); /* x = c */
 192                ret = optimize_rect(rnew, p);
 193                if (ret != NLOPT_SUCCESS) { free(rnew); return ret; }
 194                if (!rb_tree_insert(&p->rtree, rnew)) {
 195                     free(rnew); return NLOPT_OUT_OF_MEMORY;
 196                }
 197           }
 198      }
 199      if (p->minf < minf_start && nlopt_stop_f(p->stop, p->minf, minf_start))
 200           return NLOPT_FTOL_REACHED;
 201      return NLOPT_SUCCESS;
 202 }
 203
 204 /************************************************************************/
 205
 206 nlopt_result cdirect_hybrid_unscaled(int n, nlopt_func f, void *f_data,
 207                                      const double *lb, const double *ub,
 208                                      double *x,
 209                                      double *minf,
 210                                      nlopt_stopping *stop,
 211                                      nlopt_algorithm local_alg,
 212                                      int local_maxeval,
 213                                      int randomized_div)
 214 {
 215      params p;
 216      int i;
 217      double *rnew;
 218      nlopt_result ret = NLOPT_OUT_OF_MEMORY;
 219
 220      p.n = n;
 221      p.L = 3*n+3;
 222      p.lb = lb; p.ub = ub;
 223      p.stop = stop;
 224      p.f = f;
 225      p.f_data = f_data;
 226      p.minf = HUGE_VAL;
 227      p.xmin = x;
 228      p.age = 0;
 229      p.work = 0;
 230      p.local_alg = local_alg;
 231      p.local_maxeval = local_maxeval;
 232      p.randomized_div = randomized_div;
 233
 234      rb_tree_init(&p.rtree, cdirect_hyperrect_compare);
 235      p.work = (double *) malloc(sizeof(double) * (2*n));
 236      if (!p.work) goto done;
 237
 238      if (!(rnew = (double *) malloc(sizeof(double) * p.L))) goto done;
 239      for (i = 0; i < n; ++i) {
 240           rnew[3+i] = rnew[3+n+i] = 0.5 * (lb[i] + ub[i]);
 241           rnew[3+2*n+i] = ub[i] - lb[i];
 242      }
 243      rnew[0] = longest(n, rnew+2*n);
 244      rnew[2] = p.age--;
 245      ret = optimize_rect(rnew, &p);
 246      if (ret != NLOPT_SUCCESS) { free(rnew); goto done; }
 247      if (!rb_tree_insert(&p.rtree, rnew)) { free(rnew); goto done; }
 248
 249      do {
 250           ret = divide_largest(&p);
 251      } while (ret == NLOPT_SUCCESS);
 252
 253  done:
 254      rb_tree_destroy_with_keys(&p.rtree);
 255      free(p.work);
 256
 257      *minf = p.minf;
 258      return ret;
 259 }
 260
 261 /* rescaled to unit hypercube so that all x[i] are weighted equally  */
 262 nlopt_result cdirect_hybrid(int n, nlopt_func f, void *f_data,
 263                             const double *lb, const double *ub,
 264                             double *x,
 265                             double *minf,
 266                             nlopt_stopping *stop,
 267                             nlopt_algorithm local_alg,
 268                             int local_maxeval,
 269                             int randomized_div)
 270 {
 271      cdirect_uf_data d;
 272      nlopt_result ret;
 273      const double *xtol_abs_save;
 274      int i;
 275
 276      d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
 277      d.x = (double *) malloc(sizeof(double) * n*4);
 278      if (!d.x) return NLOPT_OUT_OF_MEMORY;
 279
 280      for (i = 0; i < n; ++i) {
 281           x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
 282           d.x[n+i] = 0;
 283           d.x[2*n+i] = 1;
 284           d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
 285      }
 286      xtol_abs_save = stop->xtol_abs;
 287      stop->xtol_abs = d.x + 3*n;
 288      ret = cdirect_hybrid_unscaled(n, cdirect_uf, &d, d.x+n, d.x+2*n,
 289                                    x, minf, stop, local_alg, local_maxeval,
 290                                    randomized_div);
 291      stop->xtol_abs = xtol_abs_save;
 292      for (i = 0; i < n; ++i)
 293           x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);
 294      free(d.x);
 295      return ret;
 296 }