crs/crs.c

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