chiark / gitweb /
add initial re-implementation of DIRECT (still too slow because convex hull is re...
authorstevenj <stevenj@alum.mit.edu>
Sun, 26 Aug 2007 14:47:51 +0000 (10:47 -0400)
committerstevenj <stevenj@alum.mit.edu>
Sun, 26 Aug 2007 14:47:51 +0000 (10:47 -0400)
darcs-hash:20070826144751-c8de0-88f748e6b9430e210eac28a0258f5a03389b8865.gz

Makefile.am
api/Makefile.am
api/nlopt.c
cdirect/Makefile.am [new file with mode: 0644]
cdirect/cdirect.c [new file with mode: 0644]
cdirect/cdirect.h [new file with mode: 0644]
configure.ac

index a54a7fec8acdc26ab2f36340386e053db1386d4a..cb31dc30cb017995d6f4263c8f899174dc652ae4 100644 (file)
@@ -3,11 +3,13 @@ lib_LTLIBRARIES = libnlopt.la
 
 ACLOCAL_AMFLAGS=-I ./m4
 
-SUBDIRS= util lbfgs subplex direct stogo api . test
+SUBDIRS= util lbfgs subplex direct cdirect stogo api . test
 EXTRA_DIST=COPYRIGHT autogen.sh nlopt.pc.in m4
 
 libnlopt_la_SOURCES = 
-libnlopt_la_LIBADD = lbfgs/liblbfgs.la subplex/libsubplex.la direct/libdirect.la stogo/libstogo.la api/libapi.la util/libutil.la
+libnlopt_la_LIBADD = lbfgs/liblbfgs.la subplex/libsubplex.la   \
+direct/libdirect.la cdirect/libcdirect.la stogo/libstogo.la    \
+api/libapi.la util/libutil.la
 
 libnlopt_la_LDFLAGS = -no-undefined -version-info @SHARED_VERSION_INFO@
 
index 71a0fc2799782849b19a6e493b64a54f9d0341a5..9d3dbf21ec048cc119875377b24a132adf821efe 100644 (file)
@@ -1,4 +1,4 @@
-AM_CPPFLAGS = -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/lbfgs -I$(top_srcdir)/subplex -I$(top_srcdir)/util
+AM_CPPFLAGS = -I$(top_srcdir)/cdirect -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/lbfgs -I$(top_srcdir)/subplex -I$(top_srcdir)/util
 
 include_HEADERS = nlopt.h
 noinst_LTLIBRARIES = libapi.la
index 2a3258065901b88a745dc2234c8703095dc0d86b..f29e4ca87574965d4f80a40719dba923bcb97de2 100644 (file)
@@ -175,6 +175,8 @@ static double f_direct(int n, const double *x, int *undefined, void *data_)
 
 #include "l-bfgs-b.h"
 
+#include "cdirect.h"
+
 /*************************************************************************/
 
 /* same as nlopt_minimize, but xtol_abs is required to be non-NULL */
@@ -221,7 +223,12 @@ static nlopt_result nlopt_minimize_(
 
      switch (algorithm) {
         case NLOPT_GLOBAL_DIRECT:
-        case NLOPT_GLOBAL_DIRECT_L: {
+        case NLOPT_GLOBAL_DIRECT_L: 
+#if 0
+             return cdirect(n, f, f_data, lb, ub, x, fmin, &stop, 0, 
+                            algorithm == NLOPT_GLOBAL_DIRECT ? 0 : 1);
+#endif
+        {
              int iret;
              d.xtmp = (double *) malloc(sizeof(double) * n*2);
              if (!d.xtmp) return NLOPT_OUT_OF_MEMORY;
diff --git a/cdirect/Makefile.am b/cdirect/Makefile.am
new file mode 100644 (file)
index 0000000..968a0aa
--- /dev/null
@@ -0,0 +1,6 @@
+AM_CPPFLAGS = -I$(top_srcdir)/util -I$(top_srcdir)/api
+
+noinst_LTLIBRARIES = libcdirect.la
+libcdirect_la_SOURCES = cdirect.c cdirect.h
+
+EXTRA_DIST = 
diff --git a/cdirect/cdirect.c b/cdirect/cdirect.c
new file mode 100644 (file)
index 0000000..5b697ab
--- /dev/null
@@ -0,0 +1,461 @@
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "nlopt-util.h"
+#include "nlopt.h"
+#include "cdirect.h"
+#include "config.h"
+
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
+/***************************************************************************/
+/* basic data structure:
+ *
+ * a hyper-rectangle is stored as an array of length 2n+2, where [0]
+ * is the value of the function at the center, [1] is the "size"
+ * measure of the rectangle, [2..n+1] are the coordinates of the
+ * center, and [n+2..2n+1] are the widths of the sides.
+ *
+ * a list of rectangles is just an array of N hyper-rectangles
+ * stored as an N x (2n+1) in row-major order.  Generally,
+ * we allocate more than we need, allocating Na hyper-rectangles.
+ *
+ * n > 0 always
+ */
+
+#define RECT_LEN(n) (2*(n)+2) /* number of double values in a hyperrect */
+
+/* parameters of the search algorithm and various information that
+   needs to be passed around */
+typedef struct {
+     double magic_eps; /* Jones' epsilon parameter (1e-4 is recommended) */
+     int which_diam; /* which measure of hyper-rectangle diam to use:
+                       0 = Jones, 1 = Gablonsky */
+     const double *lb, *ub;
+     nlopt_stopping *stop; /* stopping criteria */
+     int n;
+     nlopt_func f; void *f_data;
+     double *work; /* workspace, of length >= 2*n */
+     int *iwork, iwork_len; /* workspace, of length iwork_len >= n */
+     double fmin, *xmin; /* minimum so far */
+} params;
+
+/***************************************************************************/
+
+/* evaluate the "diameter" (d) of a rectangle of widths w[n] */
+static double rect_diameter(int n, const double *w, const params *p)
+{
+     int i;
+     if (p->which_diam == 0) { /* Jones measure */
+         double sum = 0;
+         for (i = 0; i < n; ++i)
+              sum += w[i] * w[i];
+         return sqrt(sum) * 0.5; /* distance from center to a vertex */
+     }
+     else { /* Gablonsky measure */
+         double maxw = 0;
+         for (i = 0; i < n; ++i)
+              if (w[i] > maxw)
+                   maxw = w[i];
+         return w[i] * 0.5; /* half-width of longest side */
+     }
+}
+
+#define CUBE_TOL 5e-2 /* fractional tolerance to call something a "cube" */
+
+/* return true if the elements of w[n] (> 0) are all equal to within a
+   fractional tolerance of tol (i.e. they are the widths of a hypercube) */
+static int all_equal(int n, const double *w, double tol)
+{
+     double wmin, wmax;
+     int i;
+     wmin = wmax = w[0];
+     for (i = 1; i < n; ++i) {
+         if (w[i] < wmin) wmin = w[i];
+         if (w[i] > wmax) wmax = w[i];
+     }
+     return (wmax - wmin) < tol * wmax;
+}
+
+static double *alloc_rects(int n, int *Na, double *rects, int newN)
+{
+     if (newN <= *Na)
+         return rects;
+     else {
+         (*Na) += newN;
+         return realloc(rects, sizeof(double) * RECT_LEN(n) * (*Na));
+     }
+}
+#define ALLOC_RECTS(n, Nap, rects, newN) if (!(rects = alloc_rects(n, Nap, rects, newN))) return NLOPT_OUT_OF_MEMORY
+
+static double *fv_qsort = 0;
+static int sort_fv_compare(const void *a_, const void *b_)
+{
+     int a = *((const int *) a_), b = *((const int *) b_);
+     double fa = MIN(fv_qsort[2*a], fv_qsort[2*a+1]);
+     double fb = MIN(fv_qsort[2*b], fv_qsort[2*b+1]);
+     if (fa < fb)
+         return -1;
+     else if (fa > fb)
+         return +1;
+     else
+         return 0;
+}
+static void sort_fv(int n, double *fv, int *isort)
+{
+     int i;
+     for (i = 0; i < n; ++i) isort[i] = i;
+     fv_qsort = fv; /* not re-entrant, sigh... */
+     qsort(isort, (unsigned) n, sizeof(int), sort_fv_compare);
+     fv_qsort = 0;
+}
+
+static double function_eval(const double *x, params *p) {
+     double f = p->f(p->n, x, NULL, p->f_data);
+     if (f < p->fmin) {
+         p->fmin = f;
+         memcpy(p->xmin, x, sizeof(double) * p->n);
+     }
+     p->stop->nevals++;
+     return f;
+}
+#define FUNCTION_EVAL(fv,x,p) fv = function_eval(x, p); if (p->fmin < p->stop->fmin_max) return NLOPT_FMIN_MAX_REACHED; else if (nlopt_stop_evals((p)->stop)) return NLOPT_MAXEVAL_REACHED; else if (nlopt_stop_time((p)->stop)) return NLOPT_MAXTIME_REACHED
+
+#define THIRD (0.3333333333333333333333)
+
+
+/* divide rectangle idiv in the list rects */
+static nlopt_result divide_rect(int *N, int *Na, double **rects, int idiv,
+                               params *p)
+{
+     int i;
+     const const int n = p->n;
+     const int L = RECT_LEN(n);
+     double *r = *rects;
+     double *c = r + L*idiv + 2; /* center of rect to divide */
+     double *w = c + n; /* widths of rect to divide */
+
+     if (all_equal(n, w, CUBE_TOL)) { /* divide all dimensions */
+         double *fv = p->work;
+         int *isort = p->iwork;
+         for (i = 0; i < n; ++i) {
+              double csave = c[i];
+              c[i] = csave - w[i] * THIRD;
+              FUNCTION_EVAL(fv[2*i], c, p);
+              c[i] = csave + w[i] * THIRD;
+              FUNCTION_EVAL(fv[2*i+1], c, p);
+              c[i] = csave;
+         }
+         sort_fv(n, fv, isort);
+         ALLOC_RECTS(n, Na, r, (*N)+2*n); 
+         *rects = r; c = r + L*idiv + 2; w = c + n;
+         for (i = 0; i < n; ++i) {
+              int k;
+              w[isort[i]] *= THIRD;
+              r[L*idiv + 1] = rect_diameter(n, w, p);
+              for (k = 0; k <= 1; ++k) {
+                   memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
+                   r[L*(*N) + 2 + isort[i]] += w[isort[i]] * (2*k-1);
+                   r[L*(*N)] = fv[2*isort[i]+k];
+                   ++(*N);
+              }
+         }
+     }
+     else { /* divide longest side by 3 and split off 2 new rectangles */
+         int imax = 0;
+         double wmax = w[0];
+         for (i = 1; i < n; ++i)
+              if (w[i] > wmax)
+                   wmax = w[imax = i];
+         ALLOC_RECTS(n, Na, r, (*N)+2);
+         *rects = r; c = r + L*idiv + 2; w = c + n;
+         w[imax] *= THIRD;
+         r[L*idiv + 1] = rect_diameter(n, w, p);
+         for (i = 0; i <= 1; ++i) {
+              memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
+              r[L*(*N) + 2 + imax] += w[imax] * (2*i-1); /* move center */
+              ++(*N);
+              FUNCTION_EVAL(r[L*((*N)-1)], r + L*((*N)-1) + 2, p);
+         }
+     }
+     return NLOPT_SUCCESS;
+}
+
+/***************************************************************************/
+/* O(N log N) convex hull algorithm, used later to find the potentially
+   optimal points */
+
+/* sort ihull by xy in lexographic order by x,y */
+static int s_qsort = 1; static double *xy_qsort = 0;
+static int sort_xy_compare(const void *a_, const void *b_)
+{
+     int a = *((const int *) a_), b = *((const int *) b_);
+     double xa = xy_qsort[a*s_qsort+1], xb = xy_qsort[b*s_qsort+1];
+     if (xa < xb) return -1;
+     else if (xb < xa) return +1;
+     else {
+         double ya = xy_qsort[a*s_qsort], yb = xy_qsort[b*s_qsort];
+         if (ya < yb) return -1;
+         else if (ya > yb) return +1;
+         else return 0;
+     }
+}
+static void sort_xy(int N, double *xy, int s, int *isort)
+{
+     int i;
+
+     for (i = 0; i < N; ++i) isort[i] = i;
+     s_qsort = s; xy_qsort = xy;
+     qsort(isort, (unsigned) N, sizeof(int), sort_xy_compare);
+     xy_qsort = 0;
+}
+
+/* Find the lower convex hull of a set of points (xy[s*i+1], xy[s*i]), where
+   0 <= i < N and s >= 2.
+
+   The return value is the number of points in the hull, with indices
+   stored in ihull.  ihull and is should point to arrays of length >= N.
+
+   Note that we don't allow redundant points along the same line in the
+   hull, similar to Gablonsky's version of DIRECT and differing from
+   Jones'. */
+static int convex_hull(int N, double *xy, int s, int *ihull, int *is)
+{
+     int minmin; /* first index (smallest y) with min x */
+     int minmax; /* last index (largest y) with min x */
+     int maxmin; /* first index (smallest y) with max x */
+     int maxmax; /* last index (largest y) with max x */
+     int i, nhull = 0;
+     double minslope;
+     double xmin, xmax;
+
+     /* Monotone chain algorithm [Andrew, 1979]. */
+
+     sort_xy(N, xy, s, is);
+
+     xmin = xy[s*is[minmin=0]+1]; xmax = xy[s*is[maxmax=N-1]+1];
+
+     if (xmin == xmax) { /* degenerate case */
+         ihull[nhull++] = is[minmin];
+         return nhull;
+     }
+
+     for (minmax = minmin; minmax+1 < N && xy[s*is[minmax+1]+1]==xmin; 
+         ++minmax);
+     for (maxmin = maxmax; maxmin-1>=0 && xy[s*is[maxmin-1]+1]==xmax;
+         --maxmin);
+
+     minslope = (xy[s*is[maxmin]] - xy[s*is[minmin]]) / (xmax - xmin);
+
+     ihull[nhull++] = is[minmin];
+     for (i = minmax + 1; i < maxmin; ++i) {
+         int k = is[i];
+         if (xy[s*k] > xy[s*is[minmin]] + (xy[s*k+1] - xmin) * minslope)
+              continue;
+         /* remove points until we are making a "left turn" to k */
+         while (nhull > 1) {
+              int t1 = ihull[nhull - 1], t2 = ihull[nhull - 2];
+              /* cross product (t1-t2) x (k-t2) > 0 for a left turn: */
+              if ((xy[s*t1+1]-xy[s*t2+1]) * (xy[s*k]-xy[s*t2])
+                  - (xy[s*t1]-xy[s*t2]) * (xy[s*k+1]-xy[s*t2+1]) > 0)
+                   break;
+              --nhull;
+         }
+         ihull[nhull++] = k;
+     }
+     ihull[nhull++] = is[maxmin];
+     return nhull;
+}
+
+/***************************************************************************/
+
+static int small(double *w, params *p)
+{
+     int i;
+     for (i = 0; i < p->n; ++i)
+         if (w[i] > p->stop->xtol_abs[i] &&
+             w[i] > (p->ub[i] - p->lb[i]) * p->stop->xtol_rel)
+              return 0;
+     return 1;
+}
+
+static nlopt_result divide_good_rects(int *N, int *Na, double **rects, 
+                                     params *p)
+{
+     const int n = p->n;
+     const int L = RECT_LEN(n);
+     int *ihull, nhull, i, xtol_reached = 1, divided_some = 0;
+     double *r = *rects;
+     double magic_eps = p->magic_eps;
+
+     if (p->iwork_len < n + 2*(*N)) {
+         p->iwork_len = p->iwork_len + n + 2*(*N);
+         p->iwork = (int *) realloc(p->iwork, sizeof(int) * p->iwork_len);
+         if (!p->iwork)
+              return NLOPT_OUT_OF_MEMORY;
+     }
+     ihull = p->iwork;
+     nhull = convex_hull(*N, r, L, ihull, ihull + *N);
+ divisions:
+     for (i = 0; i < nhull; ++i) {
+         double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
+         if (i > 0)
+              K1 = (r[L*ihull[i]] - r[L*ihull[i-1]]) /
+                   (r[L*ihull[i]+1] - r[L*ihull[i-1]+1]);
+         if (i < nhull-1)
+              K1 = (r[L*ihull[i]] - r[L*ihull[i+1]]) /
+                   (r[L*ihull[i]+1] - r[L*ihull[i+1]+1]);
+         K = MAX(K1, K2);
+         if (r[L*ihull[i]] - K * r[L*ihull[i]+1]
+             <= p->fmin - magic_eps * fabs(p->fmin)) {
+              /* "potentially optimal" rectangle, so subdivide */
+              divided_some = 1;
+              nlopt_result ret;
+              ret = divide_rect(N, Na, rects, ihull[i], p);
+              r = *rects; /* may have grown */
+              if (ret != NLOPT_SUCCESS) return ret;
+              xtol_reached = xtol_reached && small(r + L*ihull[i] + 2+n, p);
+         }
+     }
+     if (!divided_some) {
+         if (magic_eps != 0) {
+              magic_eps = 0;
+              goto divisions; /* try again */
+         }
+         else { /* WTF? divide largest rectangle */
+              double wmax = r[1];
+              int imax = 0;
+              for (i = 1; i < *N; ++i)
+                   if (r[L*i+1] > wmax)
+                        wmax = r[L*(imax=i)+1];
+              return divide_rect(N, Na, rects, imax, p);
+         }
+     }
+     return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
+}
+
+/***************************************************************************/
+
+nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
+                             const double *lb, const double *ub,
+                             double *x,
+                             double *fmin,
+                             nlopt_stopping *stop,
+                             double magic_eps, int which_alg)
+{
+     params p;
+     double *rects;
+     int Na = 100, N = 1, i, x_center = 1;
+     nlopt_result ret = NLOPT_OUT_OF_MEMORY;
+
+     p.magic_eps = magic_eps;
+     p.which_diam = which_alg & 1;
+     p.lb = lb; p.ub = ub;
+     p.stop = stop;
+     p.n = n;
+     p.f = f;
+     p.f_data = f_data;
+     p.xmin = x;
+     p.fmin = f(n, x, NULL, f_data); stop->nevals++;
+     p.work = 0;
+     p.iwork = 0;
+     rects = 0;
+     p.work = (double *) malloc(sizeof(double) * 2*n);
+     if (!p.work) goto done;
+     rects = (double *) malloc(sizeof(double) * Na * RECT_LEN(n));
+     if (!rects) goto done;
+     p.iwork = (int *) malloc(sizeof(int) * (p.iwork_len = 2*Na + n));
+     if (!p.iwork) goto done;
+
+     for (i = 0; i < n; ++i) {
+         rects[2+i] = 0.5 * (lb[i] + ub[i]);
+         x_center = x_center
+              && (fabs(rects[2+i]-x[i]) < 1e-13*(1+fabs(x[i])));
+         rects[2+n+i] = ub[i] - lb[i];
+     }
+     rects[1] = rect_diameter(n, rects+2+n, &p);
+     if (x_center)
+         rects[0] = p.fmin; /* avoid computing f(center) twice */
+     else
+         rects[0] = function_eval(rects+2, &p);
+
+     ret = divide_rect(&N, &Na, &rects, 0, &p);
+     if (ret != NLOPT_SUCCESS) goto done;
+
+     while (1) {
+         double fmin0 = p.fmin;
+         ret = divide_good_rects(&N, &Na, &rects, &p);
+         if (ret != NLOPT_SUCCESS) goto done;
+         if (nlopt_stop_f(p.stop, p.fmin, fmin0)) {
+              ret = NLOPT_FTOL_REACHED;
+              goto done;
+         }
+     }
+
+ done:
+     free(p.iwork);
+     free(rects);
+     free(p.work);
+             
+     *fmin = p.fmin;
+     return ret;
+}
+
+/* in the conventional DIRECT-type algorithm, we first rescale our
+   coordinates to a unit hypercube ... we do this simply by
+   wrapping cdirect() around cdirect_unscaled(). */
+
+typedef struct {
+     nlopt_func f;
+     void *f_data;
+     double *x;
+     const double *lb, *ub;
+} uf_data;
+static double uf(int n, const double *xu, double *grad, void *d_)
+{
+     uf_data *d = (uf_data *) d_;
+     double f;
+     int i;
+     for (i = 0; i < n; ++i)
+         d->x[i] = d->lb[i] + xu[i] * (d->ub[i] - d->lb[i]);
+     f = d->f(n, d->x, grad, d->f_data);
+     if (grad)
+         for (i = 0; i < n; ++i)
+              grad[i] *= d->ub[i] - d->lb[i];
+     return f;
+}
+
+nlopt_result cdirect(int n, nlopt_func f, void *f_data,
+                     const double *lb, const double *ub,
+                     double *x,
+                     double *fmin,
+                     nlopt_stopping *stop,
+                     double magic_eps, int which_alg)
+{
+     uf_data d;
+     nlopt_result ret;
+     const double *xtol_abs_save;
+     int i;
+
+     d.f = f; d.f_data = f_data; d.lb = lb; d.ub = ub;
+     d.x = (double *) malloc(sizeof(double) * n*4);
+     if (!d.x) return NLOPT_OUT_OF_MEMORY;
+     
+     for (i = 0; i < n; ++i) {
+         x[i] = (x[i] - lb[i]) / (ub[i] - lb[i]);
+         d.x[n+i] = 0;
+         d.x[2*n+i] = 1;
+         d.x[3*n+i] = stop->xtol_abs[i] / (ub[i] - lb[i]);
+     }
+     xtol_abs_save = stop->xtol_abs;
+     stop->xtol_abs = d.x + 3*n;
+     ret = cdirect_unscaled(n, uf, &d, d.x+n, d.x+2*n, x, fmin, stop,
+                           magic_eps, which_alg);
+     stop->xtol_abs = xtol_abs_save;
+     for (i = 0; i < n; ++i)
+         x[i] = lb[i]+ x[i] * (ub[i] - lb[i]);
+     free(d.x);
+     return ret;
+}
diff --git a/cdirect/cdirect.h b/cdirect/cdirect.h
new file mode 100644 (file)
index 0000000..f101cda
--- /dev/null
@@ -0,0 +1,30 @@
+#ifndef CDIRECT_H
+#define CDIRECT_H
+
+#include "nlopt-util.h"
+#include "nlopt.h"
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif /* __cplusplus */
+
+extern nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
+                                    const double *lb, const double *ub,
+                                    double *x,
+                                    double *fmin,
+                                    nlopt_stopping *stop,
+                                    double magic_eps, int which_alg);
+
+extern nlopt_result cdirect(int n, nlopt_func f, void *f_data,
+                           const double *lb, const double *ub,
+                           double *x,
+                           double *fmin,
+                           nlopt_stopping *stop,
+                           double magic_eps, int which_alg);
+
+#ifdef __cplusplus
+}  /* extern "C" */
+#endif /* __cplusplus */
+
+#endif /* DIRECT_H */
index 377d96a425dac2ddd06284019b28408f4f297511..46ca6b6fc162cc8cab8638e6ef0cc7179ae3af5e 100644 (file)
@@ -96,6 +96,7 @@ AC_CONFIG_FILES([
    api/Makefile
    util/Makefile
    direct/Makefile
+   cdirect/Makefile
    stogo/Makefile
    lbfgs/Makefile
    subplex/Makefile