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