crs/crs.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 <stdlib.h>
  24 #include <string.h>
  25
  26 #include "crs.h"
  27 #include "redblack.h"
  28
  29 /* Controlled Random Search 2 (CRS2) with "local mutation", as defined
  30    by:
  31        P. Kaelo and M. M. Ali, "Some variants of the controlled random
  32        search algorithm for global optimization," J. Optim. Theory Appl.
  33        130 (2), 253-264 (2006).
  34 */
  35
  36 typedef struct {
  37      int n; /* # dimensions */
  38      const double *lb, *ub;
  39      nlopt_stopping *stop; /* stopping criteria */
  40      nlopt_func f; void *f_data;
  41
  42      int N; /* # points in population */
  43      double *ps; /* population array N x (n+1) of tuples [f(x), x] */
  44      double *p; /* single point array (length n+1), for temp use */
  45      rb_tree t; /* red-black tree of population, sorted by f(x) */
  46      nlopt_sobol s; /* sobol data for LDS point generation, or NULL
  47                        to use pseudo-random numbers */
  48 } crs_data;
  49
  50 /* sort order in red-black tree: keys [f(x), x] are sorted by f(x) */
  51 static int crs_compare(double *k1, double *k2)
  52 {
  53      if (*k1 < *k2) return -1;
  54      if (*k1 > *k2) return +1;
  55      return k1 - k2; /* tie-breaker */
  56 }
  57
  58 /* set x to a random trial value, as defined by CRS:
  59      x = 2G - x_n,
  60    where x_0 ... x_n are distinct points in the population
  61    with x_0 the current best point and the other points are random,
  62    and G is the centroid of x_0...x_{n-1} */
  63 static void random_trial(crs_data *d, double *x, rb_node *best)
  64 {
  65      int n = d->n, n1 = n+1, i, k, i0, jn;
  66      double *ps = d->ps, *xi;
  67
  68      /* initialize x to x_0 = best point */
  69      memcpy(x, best->k + 1, sizeof(double) * n);
  70      i0 = (best->k - ps) / n1;
  71
  72      jn = nlopt_iurand(n); /* which of remaining n points is "x_n",
  73                               i.e. which to reflect through ...
  74                               this is necessary since we generate
  75                               the remaining points in order, so
  76                               just picking the last point would not
  77                               be very random */
  78
  79      /* use "method A" from
  80
  81            Jeffrey Scott Vitter, "An efficient algorithm for
  82            sequential random sampling," ACM Trans. Math. Soft. 13 (1),
  83            58--67 (1987).
  84
  85         to randomly pick n distinct points out of the remaining N-1 (not
  86         including i0!).  (The same as "method S" in Knuth vol. 2.)
  87         This method requires O(N) time, which is fine in our case
  88         (there are better methods if n << N). */
  89      {
  90           int Nleft = d->N - 1, nleft = n;
  91           int Nfree = Nleft - nleft;
  92           i = 0; i += i == i0;
  93           while (nleft > 1) {
  94                double q = ((double) Nfree) / Nleft;
  95                double v = nlopt_urand(0., 1.);
  96                while (q > v) {
  97                     ++i; i += i == i0;
  98                     --Nfree; --Nleft;
  99                     q = (q * Nfree) / Nleft;
 100                }
 101                xi = ps + n1 * i + 1;
 102                if (jn-- == 0) /* point to reflect through */
 103                     for (k = 0; k < n; ++k) x[k] -= xi[k] * (0.5*n);
 104                else /* point to include in centroid */
 105                     for (k = 0; k < n; ++k) x[k] += xi[k];
 106                ++i; i += i == i0;
 107                --Nleft; --nleft;
 108           }
 109           i += nlopt_iurand(Nleft); i += i == i0;
 110           xi = ps + n1 * i + 1;
 111           if (jn-- == 0) /* point to reflect through */
 112                for (k = 0; k < n; ++k) x[k] -= xi[k] * (0.5*n);
 113           else /* point to include in centroid */
 114                for (k = 0; k < n; ++k) x[k] += xi[k];
 115      }
 116      for (k = 0; k < n; ++k) {
 117           x[k] *= 2.0 / n; /* renormalize */
 118           if (x[k] > d->ub[k]) x[k] = d->ub[k];
 119           else if (x[k] < d->lb[k]) x[k] = d->lb[k];
 120      }
 121 }
 122
 123 #define NUM_MUTATION 1 /* # "local mutation" steps to try if trial fails */
 124
 125 static nlopt_result crs_trial(crs_data *d)
 126 {
 127      rb_node *best = rb_tree_min(&d->t);
 128      rb_node *worst = rb_tree_max(&d->t);
 129      int mutation = NUM_MUTATION;
 130      int i, n = d->n;
 131      random_trial(d, d->p + 1, best);
 132      do {
 133           d->p[0] = d->f(n, d->p + 1, NULL, d->f_data);
 134           d->stop->nevals++;
 135           if (nlopt_stop_forced(d->stop)) return NLOPT_FORCED_STOP;
 136           if (d->p[0] < worst->k[0]) break;
 137           if (nlopt_stop_evals(d->stop)) return NLOPT_MAXEVAL_REACHED;
 138           if (nlopt_stop_time(d->stop)) return NLOPT_MAXTIME_REACHED;
 139           if (mutation) {
 140                for (i = 0; i < n; ++i) {
 141                     double w = nlopt_urand(0.,1.);
 142                     d->p[1+i] = best->k[1+i] * (1 + w) - w * d->p[1+i];
 143                     if (d->p[1+i] > d->ub[i]) d->p[1+i] = d->ub[i];
 144                     else if (d->p[1+i] < d->lb[i]) d->p[1+i] = d->lb[i];
 145                }
 146                mutation--;
 147           }
 148           else {
 149                random_trial(d, d->p + 1, best);
 150                mutation = NUM_MUTATION;
 151           }
 152      } while (1);
 153      memcpy(worst->k, d->p, sizeof(double) * (n+1));
 154      rb_tree_resort(&d->t, worst);
 155      return NLOPT_SUCCESS;
 156 }
 157
 158 static void crs_destroy(crs_data *d)
 159 {
 160      nlopt_sobol_destroy(d->s);
 161      rb_tree_destroy(&d->t);
 162      free(d->ps);
 163 }
 164
 165 static nlopt_result crs_init(crs_data *d, int n, const double *x,
 166                              const double *lb, const double *ub,
 167                              nlopt_stopping *stop, nlopt_func f, void *f_data,
 168                              int population, int lds)
 169 {
 170      int i;
 171
 172      if (!population) {
 173           /* TODO: how should we set the default population size?
 174              the Kaelo and Ali paper suggests 10*(n+1), but should
 175              we add more random points if maxeval is large, or... ? */
 176           d->N = 10 * (n + 1); /* heuristic initial population size */
 177      }
 178      else
 179           d->N = population;
 180      if (d->N < n + 1) /* population must be big enough for a simplex */
 181           return NLOPT_INVALID_ARGS;
 182
 183      d->n = n;
 184      d->stop = stop;
 185      d->f = f; d->f_data = f_data;
 186      d->ub = ub; d->lb = lb;
 187      d->ps = (double *) malloc(sizeof(double) * (n + 1) * (d->N + 1));
 188      if (!d->ps) return NLOPT_OUT_OF_MEMORY;
 189      d->p = d->ps + d->N * (n+1);
 190      rb_tree_init(&d->t, crs_compare);
 191
 192      /* we can either use pseudorandom points, as in the original CRS
 193         algorithm, or use a low-discrepancy Sobol' sequence ... I tried
 194         the latter, however, and it doesn't seem to help, probably
 195         because we are only generating a small number of random points
 196         to start with */
 197      d->s = lds ? nlopt_sobol_create((unsigned) n) : NULL;
 198      nlopt_sobol_skip(d->s, (unsigned) d->N, d->ps + 1);
 199
 200      /* generate initial points randomly, plus starting guess x */
 201      memcpy(d->ps + 1, x, sizeof(double) * n);
 202      d->ps[0] = f(n, x, NULL, f_data);
 203      stop->nevals++;
 204      if (!rb_tree_insert(&d->t, d->ps)) return NLOPT_OUT_OF_MEMORY;
 205      if (d->ps[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
 206      if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
 207      if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
 208      for (i = 1; i < d->N; ++i) {
 209           double *k = d->ps + i*(n+1);
 210           if (d->s)
 211                nlopt_sobol_next(d->s, k + 1, lb, ub);
 212           else {
 213                int j;
 214                for (j = 0; j < n; ++j)
 215                     k[1 + j] = nlopt_urand(lb[j], ub[j]);
 216           }
 217           k[0] = f(n, k + 1, NULL, f_data);
 218           stop->nevals++;
 219           if (!rb_tree_insert(&d->t, k)) return NLOPT_OUT_OF_MEMORY;
 220           if (k[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED;
 221           if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED;
 222           if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED;
 223      }
 224
 225      return NLOPT_SUCCESS;;
 226 }
 227
 228 nlopt_result crs_minimize(int n, nlopt_func f, void *f_data,
 229                           const double *lb, const double *ub, /* bounds */
 230                           double *x, /* in: initial guess, out: minimizer */
 231                           double *minf,
 232                           nlopt_stopping *stop,
 233                           int population, /* initial population (0=default) */
 234                           int lds) /* random or low-discrepancy seq. (lds) */
 235 {
 236      nlopt_result ret;
 237      crs_data d;
 238      rb_node *best;
 239
 240      ret = crs_init(&d, n, x, lb, ub, stop, f, f_data, population, lds);
 241      if (ret < 0) return ret;
 242
 243      best = rb_tree_min(&d.t);
 244      *minf = best->k[0];
 245      memcpy(x, best->k + 1, sizeof(double) * n);
 246
 247      while (ret == NLOPT_SUCCESS) {
 248           if (NLOPT_SUCCESS == (ret = crs_trial(&d))) {
 249                best = rb_tree_min(&d.t);
 250                if (best->k[0] < *minf) {
 251                     if (best->k[0] < stop->minf_max)
 252                          ret = NLOPT_MINF_MAX_REACHED;
 253                     else if (nlopt_stop_f(stop, best->k[0], *minf))
 254                          ret = NLOPT_FTOL_REACHED;
 255                     else if (nlopt_stop_x(stop, best->k + 1, x))
 256                          ret = NLOPT_XTOL_REACHED;
 257                     *minf = best->k[0];
 258                     memcpy(x, best->k + 1, sizeof(double) * n);
 259                }
 260                if (ret != NLOPT_SUCCESS) {
 261                     if (nlopt_stop_evals(stop))
 262                          ret = NLOPT_MAXEVAL_REACHED;
 263                     else if (nlopt_stop_time(stop))
 264                          ret = NLOPT_MAXTIME_REACHED;
 265                }
 266           }
 267      }
 268      crs_destroy(&d);
 269      return ret;
 270 }