chiark / gitweb /
detect a couple more out-of-memory conditions in CRS
[nlopt.git] / 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 }