support calling both original and new DIRECT code, make direct limits more dynamic...

darcs-hash:20070829234111-c8de0-2ce2867d6444a75b3434ec6b971c0aafefd3d6b2.gz

diff --git a/api/nlopt.c b/api/nlopt.c
index 1e89fada456712eb685684ab5731039280341188..00ac1748ab0e7798bf061ac47d1b17ba1293c69b 100644 (file)
--- a/api/nlopt.c
+++ b/api/nlopt.c
@@ -43,6 +43,7 @@ void nlopt_version(int *major, int *minor, int *bugfix)
 static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = {
      "DIRECT (global)",
      "DIRECT-L (global)",
+     "Randomized DIRECT-L (global)",
      "Original DIRECT version (global)",
      "Original DIRECT-L version (global)",
      "Subplex (local)",
@@ -230,8 +231,10 @@ static nlopt_result nlopt_minimize_(
      switch (algorithm) {
         case NLOPT_GLOBAL_DIRECT:
         case NLOPT_GLOBAL_DIRECT_L: 
-             return cdirect(n, f, f_data, lb, ub, x, fmin, &stop, 1e-4, 
-                            algorithm == NLOPT_GLOBAL_DIRECT ? 0 : 1);
+        case NLOPT_GLOBAL_DIRECT_L_RANDOMIZED: 
+             return cdirect(n, f, f_data, lb, ub, x, fmin, &stop, 0.0, 
+                            (algorithm == NLOPT_GLOBAL_DIRECT ? 0 : 1)
+                            + 10 * (algorithm == NLOPT_GLOBAL_DIRECT_L_RANDOMIZED ? 2 : 0));
 
         case NLOPT_GLOBAL_ORIG_DIRECT:
         case NLOPT_GLOBAL_ORIG_DIRECT_L: 
@@ -241,9 +244,9 @@ static nlopt_result nlopt_minimize_(
              if (!d.xtmp) return NLOPT_OUT_OF_MEMORY;
              memcpy(d.xtmp + n, x, sizeof(double) * n); d.x0 = d.xtmp + n;
              iret = direct_optimize(f_direct, &d, n, lb, ub, x, fmin,
-                                    maxeval, 500, 1e-4, 1e-4,
-                                    xtol_rel, xtol_rel,
-                                    DIRECT_UNKNOWN_FGLOBAL, -1.0,
+                                    maxeval, 999, 0.0, 0.0,
+                                    pow(xtol_rel, (double) n), -1.0,
+                                    stop.fmin_max, 0.0,
                                     NULL, 
                                     algorithm == NLOPT_GLOBAL_ORIG_DIRECT
                                     ? DIRECT_ORIGINAL
@@ -263,7 +266,7 @@ static nlopt_result nlopt_minimize_(
                  case DIRECT_MAXITER_EXCEEDED:
                       return NLOPT_MAXEVAL_REACHED;
                  case DIRECT_GLOBAL_FOUND:
-                      return NLOPT_SUCCESS;
+                      return NLOPT_FMIN_MAX_REACHED;
                  case DIRECT_VOLTOL:
                  case DIRECT_SIGMATOL:
                       return NLOPT_XTOL_REACHED;

diff --git a/api/nlopt.h b/api/nlopt.h
index 516931a051d9f5fcd9dc9f43a7efd49103e892ea..8f9e30d029fcdd9af2beb769ad7ee1ffe68eff24 100644 (file)
--- a/api/nlopt.h
+++ b/api/nlopt.h
@@ -15,6 +15,7 @@ typedef enum {
 
      NLOPT_GLOBAL_DIRECT = 0,
      NLOPT_GLOBAL_DIRECT_L,
+     NLOPT_GLOBAL_DIRECT_L_RANDOMIZED,
      NLOPT_GLOBAL_ORIG_DIRECT,
      NLOPT_GLOBAL_ORIG_DIRECT_L,
 

diff --git a/cdirect/cdirect.c b/cdirect/cdirect.c
index bd799976c71d4182ee5154c8d174ea18da798c6c..aa21b96c3788efe9d6127955393449707b908a5a 100644 (file)
--- a/cdirect/cdirect.c
+++ b/cdirect/cdirect.c
@@ -245,7 +245,10 @@ static int convex_hull(rb_tree *t, double **hull)
      yminmin = n->k[1];
      xmax = nmax->k[0];
 
-     hull[nhull++] = n->k;
+     do { /* include any duplicate points at (xmin,yminmin) */
+         hull[nhull++] = n->k;
+         n = rb_tree_succ(n);
+     } while (n && n->k[0] == xmin && n->k[1] == yminmin);
      if (xmin == xmax) return nhull;
 
      /* set nmax = min mode with x == xmax */
@@ -327,7 +330,12 @@ static int convex_hull(rb_tree *t, double **hull)
          }
          hull[nhull++] = k;
      }
-     hull[nhull++] = nmax->k;
+
+     do { /* include any duplicate points at (xmax,ymaxmin) */
+         hull[nhull++] = nmax->k;
+         nmax = rb_tree_succ(nmax);
+     } while (nmax && nmax->k[0] == xmax && n->k[1] == ymaxmin);
+
      return nhull;
 }
 
@@ -359,16 +367,19 @@ static nlopt_result divide_good_rects(params *p)
  divisions:
      for (i = 0; i < nhull; ++i) {
          double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
-         if (i > 0)
-              K1 = (hull[i][1] - hull[i-1][1]) / (hull[i][0] - hull[i-1][0]);
-         if (i < nhull-1)
-              K1 = (hull[i][1] - hull[i+1][1]) / (hull[i][0] - hull[i+1][0]);
+         int im, ip;
+         for (im = i-1; im >= 0 && hull[im][0] == hull[i][0]; --im);
+         for (ip = i+1; ip < nhull && hull[ip][0] == hull[i][0]; ++ip);
+         if (im >= 0)
+              K1 = (hull[i][1] - hull[im][1]) / (hull[i][0] - hull[im][0]);
+         if (ip < nhull)
+              K1 = (hull[i][1] - hull[ip][1]) / (hull[i][0] - hull[ip][0]);
          K = MAX(K1, K2);
          if (hull[i][1] - K * hull[i][0]
              <= p->fmin - magic_eps * fabs(p->fmin)) {
               /* "potentially optimal" rectangle, so subdivide */
-              divided_some = 1;
               nlopt_result ret = divide_rect(hull[i], p);
+              divided_some = 1;
               if (ret != NLOPT_SUCCESS) return ret;
               xtol_reached = xtol_reached && small(hull[i] + 2+n, p);
          }
@@ -423,8 +434,8 @@ nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
      nlopt_result ret = NLOPT_OUT_OF_MEMORY;
 
      p.magic_eps = magic_eps;
-     p.which_diam = which_alg % 2;
-     p.which_div = 0;
+     p.which_diam = which_alg % 10;
+     p.which_div = (which_alg / 10) % 10;
      p.lb = lb; p.ub = ub;
      p.stop = stop;
      p.n = n;

diff --git a/direct/DIRect.c b/direct/DIRect.c
index 60d43aa1f2f4a32ae1a8db675e8739976ad23851..994ee050aa7f631260d47b3aff43f8236664252b 100644 (file)
--- a/direct/DIRect.c
+++ b/direct/DIRect.c
@@ -54,7 +54,7 @@
     doublereal d__1;
 
     /* constants (FIXME: change to variable?) */
-    integer MAXFUNC = *maxf <= 0 ? 100050 : *maxf + 50;
+    integer MAXFUNC = *maxf <= 0 ? 101000 : (*maxf + 1000 + *maxf / 2);
     const integer MAXDEEP = 1000;
     const integer MAXDIV = 5000;
 
@@ -99,6 +99,7 @@
     MY_ALLOC(c__, doublereal, MAXFUNC * (*n));
     MY_ALLOC(length, integer, MAXFUNC * (*n));
     MY_ALLOC(f, doublereal, MAXFUNC * 2);
+    if (*maxf <= 0) *maxf = MAXFUNC - 1000;
 
     MY_ALLOC(s, integer, MAXDIV * 2);
     MY_ALLOC(w, doublereal, (*n));

diff --git a/direct/direct.h b/direct/direct.h
index 4ceda349e978b31992d3950b41c72dd2f86615cd..8b6f916d8e7fee6cdd9814b6110da54b153dd570 100644 (file)
--- a/direct/direct.h
+++ b/direct/direct.h
@@ -44,7 +44,7 @@ extern direct_return_code direct_optimize(
      double *x, double *fmin, 
 
      int max_feval, int max_iter,
-     double reltol, double abstol,
+     double magic_eps, double magic_eps_abs,
      double volume_reltol, double sigma_reltol,
 
      double fglobal,

diff --git a/direct/direct_wrap.c b/direct/direct_wrap.c
index 406a44da221098edde824cb0d9760bec5ba97a61..3a26838be4a3e21f04ef6ebae568e8e2debdc97b 100644 (file)
--- a/direct/direct_wrap.c
+++ b/direct/direct_wrap.c
@@ -20,8 +20,12 @@
    x: an array of length dimension, set to optimum variables upon return
    fmin: on return, set to minimum f value
 
+   magic_eps, magic_eps_abs: Jones' "magic" epsilon parameter, and
+                             also an absolute version of the same
+                            (not multipled by fmin).  Jones suggests
+                            setting this to 1e-4, but 0 also works...
+
    max_feval, max_iter: maximum number of function evaluations & DIRECT iters
-   reltol, abstol: relative and absolute tolerances (0 if none)
    volume_reltol: relative tolerance on hypercube volume (0 if none)
    sigma_reltol: relative tolerance on hypercube "measure" (??) (0 if none)
 
@@ -44,7 +48,7 @@ direct_return_code direct_optimize(
      double *x, double *fmin, 
 
      int max_feval, int max_iter,
-     double reltol, double abstol,
+     double magic_eps, double magic_eps_abs,
      double volume_reltol, double sigma_reltol,
 
      double fglobal,
@@ -80,7 +84,7 @@ direct_return_code direct_optimize(
          u[i] = upper_bounds[i];
      }
      
-     direct_direct_(f, x, &dimension, &reltol, abstol,
+     direct_direct_(f, x, &dimension, &magic_eps, magic_eps_abs,
                    &max_feval, &max_iter,
                    fmin,
                    l, u,