From 80163c517ff5defa5ac1da07554417799fd3d529 Mon Sep 17 00:00:00 2001 From: stevenj Date: Thu, 6 Sep 2007 21:52:48 -0400 Subject: [PATCH] added my own implementation of controlled random search (CRS) algorithm darcs-hash:20070907015248-c8de0-5f583ad3dac1a4d688b992a894c8ce75745bbbc7.gz --- Makefile.am | 4 +- api/Makefile.am | 2 +- api/nlopt.c | 10 +- api/nlopt.h | 2 + configure.ac | 1 + crs/Makefile.am | 6 ++ crs/README | 12 +++ crs/crs.c | 240 ++++++++++++++++++++++++++++++++++++++++++++++++ util/sobolseq.c | 10 +- 9 files changed, 279 insertions(+), 8 deletions(-) create mode 100644 crs/Makefile.am create mode 100644 crs/README create mode 100644 crs/crs.c diff --git a/Makefile.am b/Makefile.am index 86cbf1d..d0949db 100644 --- a/Makefile.am +++ b/Makefile.am @@ -3,13 +3,13 @@ lib_LTLIBRARIES = libnlopt.la ACLOCAL_AMFLAGS=-I ./m4 -SUBDIRS= util subplex direct cdirect stogo praxis luksan api octave . test +SUBDIRS= util subplex direct cdirect stogo praxis luksan crs api octave . test EXTRA_DIST=COPYRIGHT autogen.sh nlopt.pc.in m4 libnlopt_la_SOURCES = libnlopt_la_LIBADD = subplex/libsubplex.la direct/libdirect.la \ cdirect/libcdirect.la stogo/libstogo.la praxis/libpraxis.la \ -luksan/libluksan.la api/libapi.la util/libutil.la +luksan/libluksan.la crs/libcrs.la api/libapi.la util/libutil.la libnlopt_la_LDFLAGS = -no-undefined -version-info @SHARED_VERSION_INFO@ diff --git a/api/Makefile.am b/api/Makefile.am index 30ba1bf..63634e0 100644 --- a/api/Makefile.am +++ b/api/Makefile.am @@ -1,4 +1,4 @@ -AM_CPPFLAGS = -I$(top_srcdir)/cdirect -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/subplex -I$(top_srcdir)/praxis -I$(top_srcdir)/luksan -I$(top_srcdir)/util +AM_CPPFLAGS = -I$(top_srcdir)/cdirect -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/subplex -I$(top_srcdir)/praxis -I$(top_srcdir)/luksan -I$(top_srcdir)/crs -I$(top_srcdir)/util include_HEADERS = nlopt.h noinst_LTLIBRARIES = libapi.la diff --git a/api/nlopt.c b/api/nlopt.c index f0a2ab1..732a5d8 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -41,7 +41,7 @@ void nlopt_version(int *major, int *minor, int *bugfix) /*************************************************************************/ -static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = { +static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = { "DIRECT (global, no-derivative)", "DIRECT-L (global, no-derivative)", "Randomized DIRECT-L (global, no-derivative)", @@ -61,6 +61,9 @@ static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = { "Truncated Newton with restarting (local, derivative-based)", "Preconditioned truncated Newton (local, derivative-based)", "Preconditioned truncated Newton with restarting (local, derivative-based)", + "Controlled random search (CRS2) with local mutation (global, no-derivative)", + "Controlled quasi-random search (CRS2) with local mutation (global, no-derivative), using Sobol' LDS", + }; const char *nlopt_algorithm_name(nlopt_algorithm a) @@ -125,6 +128,8 @@ static double f_direct(int n, const double *x, int *undefined, void *data_) #include "luksan.h" +#include "crs.h" + /*************************************************************************/ /* for "hybrid" algorithms that combine global search with some @@ -323,6 +328,9 @@ static nlopt_result nlopt_minimize_( 1 + (algorithm - NLOPT_LD_TNEWTON) % 2, 1 + (algorithm - NLOPT_LD_TNEWTON) / 2); + case NLOPT_GN_CRS2_LM: + return crs_minimize(n, f, f_data, lb, ub, x, minf, &stop, 0); + default: return NLOPT_INVALID_ARGS; } diff --git a/api/nlopt.h b/api/nlopt.h index ab08264..750d308 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -50,6 +50,8 @@ typedef enum { NLOPT_LD_TNEWTON_PRECOND, NLOPT_LD_TNEWTON_PRECOND_RESTART, + NLOPT_GN_CRS2_LM, + NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */ } nlopt_algorithm; diff --git a/configure.ac b/configure.ac index 5d18814..e4ac3d7 100644 --- a/configure.ac +++ b/configure.ac @@ -181,6 +181,7 @@ AC_CONFIG_FILES([ subplex/Makefile praxis/Makefile luksan/Makefile + crs/Makefile test/Makefile ]) diff --git a/crs/Makefile.am b/crs/Makefile.am new file mode 100644 index 0000000..14ff7cb --- /dev/null +++ b/crs/Makefile.am @@ -0,0 +1,6 @@ +AM_CPPFLAGS = -I$(top_srcdir)/util -I$(top_srcdir)/api + +noinst_LTLIBRARIES = libcrs.la +libcrs_la_SOURCES = crs.c crs.h + +EXTRA_DIST = README diff --git a/crs/README b/crs/README new file mode 100644 index 0000000..3d44ed0 --- /dev/null +++ b/crs/README @@ -0,0 +1,12 @@ +This is my implementation of the "controlled random search" (CRS2) algorithm +with the "local mutation" modification, as defined by: + + P. Kaelo and M. M. Ali, "Some variants of the controlled random + search algorithm for global optimization," J. Optim. Theory Appl. + 130 (2), 253-264 (2006). + +It is under the same MIT license as the rest of my code in NLopt (see +../COPYRIGHT). + +Steven G. Johnson +September 2007 diff --git a/crs/crs.c b/crs/crs.c new file mode 100644 index 0000000..e5c0ba8 --- /dev/null +++ b/crs/crs.c @@ -0,0 +1,240 @@ +#include +#include + +#include "crs.h" +#include "redblack.h" + +/* Controlled Random Search 2 (CRS2) with "local mutation", as defined + by: + P. Kaelo and M. M. Ali, "Some variants of the controlled random + search algorithm for global optimization," J. Optim. Theory Appl. + 130 (2), 253-264 (2006). +*/ + +typedef struct { + int n; /* # dimensions */ + const double *lb, *ub; + nlopt_stopping *stop; /* stopping criteria */ + nlopt_func f; void *f_data; + + int N; /* # points in population */ + double *ps; /* population array N x (n+1) of tuples [f(x), x] */ + double *p; /* single point array (length n+1), for temp use */ + rb_tree t; /* red-black tree of population, sorted by f(x) */ + nlopt_sobol s; /* sobol data for LDS point generation, or NULL + to use pseudo-random numbers */ +} crs_data; + +/* sort order in red-black tree: keys [f(x), x] are sorted by f(x) */ +static int crs_compare(double *k1, double *k2) +{ + if (*k1 < *k2) return -1; + if (*k1 > *k2) return +1; + return k1 - k2; /* tie-breaker */ +} + +/* set x to a random trial value, as defined by CRS: + x = 2G - x_n, + where x_0 ... x_n are distinct points in the population + with x_0 the current best point and the other points are random, + and G is the centroid of x_0...x_{n-1} */ +static void random_trial(crs_data *d, double *x, rb_node *best) +{ + int n = d->n, n1 = n+1, i, k, i0, jn; + double *ps = d->ps, *xi; + + /* initialize x to x_0 = best point */ + memcpy(x, best->k + 1, sizeof(double) * n); + i0 = (best->k - ps) / n1; + + jn = nlopt_iurand(n); /* which of remaining n points is "x_n", + i.e. which to reflect through ... + this is necessary since we generate + the remaining points in order, so + just picking the last point would not + be very random */ + + /* use "method A" from + + Jeffrey Scott Vitter, "An efficient algorithm for + sequential random sampling," ACM Trans. Math. Soft. 13 (1), + 58--67 (1987). + + to randomly pick n distinct points out of the remaining N-1 (not + including i0!). (The same as "method S" in Knuth vol. 2.) + This method requires O(N) time, which is fine in our case + (there are better methods if n << N). */ + { + int Nleft = d->N - 1, nleft = n; + int Nfree = Nleft - nleft; + i = 0; i += i == i0; + while (nleft > 1) { + double q = ((double) Nfree) / Nleft; + double v = nlopt_urand(0., 1.); + while (q > v) { + ++i; i += i == i0; + --Nfree; --Nleft; + q = (q * Nfree) / Nleft; + } + xi = ps + n1 * i + 1; + if (jn-- == 0) /* point to reflect through */ + for (k = 0; k < n; ++k) x[k] -= xi[k] * (0.5*n); + else /* point to include in centroid */ + for (k = 0; k < n; ++k) x[k] += xi[k]; + ++i; i += i == i0; + --Nleft; --nleft; + } + i += nlopt_iurand(Nleft); i += i == i0; + xi = ps + n1 * i + 1; + if (jn-- == 0) /* point to reflect through */ + for (k = 0; k < n; ++k) x[k] -= xi[k] * (0.5*n); + else /* point to include in centroid */ + for (k = 0; k < n; ++k) x[k] += xi[k]; + } + for (k = 0; k < n; ++k) { + x[k] *= 2.0 / n; /* renormalize */ + if (x[k] > d->ub[k]) x[k] = d->ub[k]; + else if (x[k] < d->lb[k]) x[k] = d->lb[k]; + } +} + +#define NUM_MUTATION 1 /* # "local mutation" steps to try if trial fails */ + +static nlopt_result crs_trial(crs_data *d) +{ + rb_node *best = rb_tree_min(&d->t); + rb_node *worst = rb_tree_max(&d->t); + int mutation = NUM_MUTATION; + int i, n = d->n; + random_trial(d, d->p + 1, best); + do { + d->p[0] = d->f(n, d->p + 1, NULL, d->f_data); + d->stop->nevals++; + if (d->p[0] < worst->k[0]) break; + if (nlopt_stop_evals(d->stop)) return NLOPT_MAXEVAL_REACHED; + if (nlopt_stop_time(d->stop)) return NLOPT_MAXTIME_REACHED; + if (mutation) { + for (i = 0; i < n; ++i) { + double w = nlopt_urand(0.,1.); + d->p[1+i] = best->k[1+i] * (1 + w) - w * d->p[1+i]; + if (d->p[1+i] > d->ub[i]) d->p[1+i] = d->ub[i]; + else if (d->p[1+i] < d->lb[i]) d->p[1+i] = d->lb[i]; + } + mutation--; + } + else { + random_trial(d, d->p + 1, best); + mutation = NUM_MUTATION; + } + } while (1); + memcpy(worst->k, d->p, sizeof(double) * (n+1)); + rb_tree_resort(&d->t, worst); + return NLOPT_SUCCESS; +} + +static void crs_destroy(crs_data *d) +{ + nlopt_sobol_destroy(d->s); + rb_tree_destroy(&d->t); + free(d->ps); +} + +static nlopt_result crs_init(crs_data *d, int n, const double *x, + const double *lb, const double *ub, + nlopt_stopping *stop, nlopt_func f, void *f_data, + int lds) +{ + int i; + + /* TODO: how should we set the initial population size? + the Kaelo and Ali paper suggests 10*(n+1), but should + we add more random points if maxeval is large, or... ? */ + d->N = 10 * (n + 1); /* heuristic initial population size */ + + d->n = n; + d->stop = stop; + d->f = f; d->f_data = f_data; + d->ub = ub; d->lb = lb; + d->ps = (double *) malloc(sizeof(double) * (n + 1) * (d->N + 1)); + if (!d->ps) return NLOPT_OUT_OF_MEMORY; + d->p = d->ps + d->N * (n+1); + rb_tree_init(&d->t, crs_compare); + + /* we can either use pseudorandom points, as in the original CRS + algorithm, or use a low-discrepancy Sobol' sequence ... I tried + the latter, however, and it doesn't seem to help, probably + because we are only generating a small number of random points + to start with */ + d->s = lds ? nlopt_sobol_create((unsigned) n) : NULL; + nlopt_sobol_skip(d->s, (unsigned) d->N, d->ps + 1); + + /* generate initial points randomly, plus starting guess x */ + memcpy(d->ps + 1, x, sizeof(double) * n); + d->ps[0] = f(n, x, NULL, f_data); + stop->nevals++; + rb_tree_insert(&d->t, d->ps); + if (d->ps[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED; + if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED; + if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED; + for (i = 1; i < d->N; ++i) { + double *k = d->ps + i*(n+1); + if (d->s) + nlopt_sobol_next(d->s, k + 1, lb, ub); + else { + int j; + for (j = 0; j < n; ++j) + k[1 + j] = nlopt_urand(lb[j], ub[j]); + } + k[0] = f(n, k + 1, NULL, f_data); + stop->nevals++; + rb_tree_insert(&d->t, k); + if (k[0] < stop->minf_max) return NLOPT_MINF_MAX_REACHED; + if (nlopt_stop_evals(stop)) return NLOPT_MAXEVAL_REACHED; + if (nlopt_stop_time(stop)) return NLOPT_MAXTIME_REACHED; + } + + return NLOPT_SUCCESS;; +} + +nlopt_result crs_minimize(int n, nlopt_func f, void *f_data, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *minf, + nlopt_stopping *stop, + int lds) /* random or low-discrepancy seq. (lds) */ +{ + nlopt_result ret; + crs_data d; + rb_node *best; + + ret = crs_init(&d, n, x, lb, ub, stop, f, f_data, lds); + if (ret < 0) return ret; + + best = rb_tree_min(&d.t); + *minf = best->k[0]; + memcpy(x, best->k + 1, sizeof(double) * n); + + while (ret == NLOPT_SUCCESS) { + if (NLOPT_SUCCESS == (ret = crs_trial(&d))) { + best = rb_tree_min(&d.t); + if (best->k[0] < *minf) { + if (best->k[0] < stop->minf_max) + ret = NLOPT_MINF_MAX_REACHED; + else if (nlopt_stop_f(stop, best->k[0], *minf)) + ret = NLOPT_FTOL_REACHED; + else if (nlopt_stop_x(stop, best->k + 1, x)) + ret = NLOPT_XTOL_REACHED; + *minf = best->k[0]; + memcpy(x, best->k + 1, sizeof(double) * n); + } + if (ret != NLOPT_SUCCESS) { + if (nlopt_stop_evals(stop)) + ret = NLOPT_MAXEVAL_REACHED; + else if (nlopt_stop_time(stop)) + ret = NLOPT_MAXTIME_REACHED; + } + } + } + crs_destroy(&d); + return ret; +} diff --git a/util/sobolseq.c b/util/sobolseq.c index e390ac5..a662fc1 100644 --- a/util/sobolseq.c +++ b/util/sobolseq.c @@ -207,7 +207,7 @@ extern void nlopt_sobol_destroy(nlopt_sobol s) /* next vector x[sdim] in Sobol sequence, with each x[i] in (0,1) */ void nlopt_sobol_next01(nlopt_sobol s, double *x) { - if (!s || !sobol_gen(s, x)) { + if (!sobol_gen(s, x)) { /* fall back on pseudo random numbers in the unlikely event that we exceed 2^32-1 points */ unsigned i; @@ -232,9 +232,11 @@ void nlopt_sobol_next(nlopt_sobol s, double *x, points equal to the largest power of 2 smaller than n */ void nlopt_sobol_skip(nlopt_sobol s, unsigned n, double *x) { - unsigned k = 1; - while (k*2 < n) k *= 2; - while (k-- > 0) sobol_gen(s, x); + if (s) { + unsigned k = 1; + while (k*2 < n) k *= 2; + while (k-- > 0) sobol_gen(s, x); + } } /************************************************************************/ -- 2.30.2