chiark / gitweb /
added PLIS variable-metric methods
authorstevenj <stevenj@alum.mit.edu>
Mon, 3 Sep 2007 20:42:10 +0000 (16:42 -0400)
committerstevenj <stevenj@alum.mit.edu>
Mon, 3 Sep 2007 20:42:10 +0000 (16:42 -0400)
darcs-hash:20070903204210-c8de0-428d0689b27e43df180260d820e68000db8bdc3a.gz

api/nlopt.c
api/nlopt.h
luksan/luksan.h
luksan/plip.c
luksan/plis.c

index f3aaca59c47cce78d6d961d5f8569836b0b05682..99d06ee50c9831b5557459b9bfe0c64b9f2a5195 100644 (file)
@@ -55,6 +55,8 @@ static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = {
      "StoGO with randomized search (global, derivative-based)",
      "Low-storage BFGS (LBFGS) (local, derivative-based)",
      "Principal-axis, praxis (local, no-derivative)",
+     "Limited-memory variable-metric, rank 1 (local, derivative-based)",
+     "Limited-memory variable-metric, rank 2 (local, derivative-based)",
 };
 
 const char *nlopt_algorithm_name(nlopt_algorithm a)
@@ -304,6 +306,11 @@ static nlopt_result nlopt_minimize_(
         case NLOPT_LD_LBFGS: 
              return luksan_plis(n, f, f_data, lb, ub, x, minf, &stop);
 
+        case NLOPT_LD_VAR1: 
+        case NLOPT_LD_VAR2: 
+             return luksan_plip(n, f, f_data, lb, ub, x, minf, &stop,
+                  algorithm == NLOPT_LD_VAR1 ? 1 : 2);
+
         default:
              return NLOPT_INVALID_ARGS;
      }
index 8045cb651293b60cd8dbdfaab344a2c81a3f0462..b0265b21ea8acd4fc8482b74f84af6a70147bb56 100644 (file)
@@ -42,6 +42,9 @@ typedef enum {
 
      NLOPT_LN_PRAXIS,
 
+     NLOPT_LD_VAR1,
+     NLOPT_LD_VAR2,
+
      NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */
 } nlopt_algorithm;
 
index 8a9d3abc115c30399e276475dccab547cd5975aa..fc218387bdbbb8d2f32f78cbbd58e2896d915c00 100644 (file)
@@ -15,6 +15,13 @@ nlopt_result luksan_plis(int n, nlopt_func f, void *f_data,
                   double *minf,
                   nlopt_stopping *stop);
 
+nlopt_result luksan_plip(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 method);
+
 /*****************************  internal routines *************************/
 
 /* mssubs.c: */
@@ -118,6 +125,11 @@ typedef struct {
      int nfg, nfh;
 } stat_common;
 
+/* number of double variables that can be stored in scratch memory
+   ... it's >= 2007, and this is in the context of scientific computation,
+   so assume that at least 10M are available, and that sizeof(double)==8 */
+#define MEMAVAIL 1310720
+
 #ifdef __cplusplus
 }  /* extern "C" */
 #endif /* __cplusplus */
index e73c7c2a8949a0ab40672d8bd6c1e440eeea5e7c..0551c18f619c06293e8e18a64c223a54665d98f8 100644 (file)
@@ -420,3 +420,101 @@ L11190:
     return;
 } /* plip_ */
 
+/* NLopt wrapper around plip_, handling dynamic allocation etc. */
+nlopt_result luksan_plip(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 method) /* 1 or 2, see below */
+{
+     int i, *ix, nb = 1;
+     double *work, *xl, *xu, *gf, *s, *xo, *go, *so, *xm, *xr, *gr;
+     double gmax, minf_est;
+     double xmax = 0; /* no maximum */
+     double tolg = 0; /* default gradient tolerance */
+     int iest = 0; /* we have no estimate of min function value */
+     int mit = 0; /* default no limit on #iterations */
+     int mfv = stop->maxeval;
+     stat_common stat;
+     int iterm;
+     int mf;
+
+     ix = (int*) malloc(sizeof(int) * n);
+     if (!ix) return NLOPT_OUT_OF_MEMORY;
+
+     /* FIXME: what should we set mf to?  The example program tlis.for
+        sets it to zero as far as I can tell, but it seems to greatly
+       improve convergence to make it > 0.  The computation time
+       per iteration, and of course the memory, seem to go as O(n * mf),
+       and we'll assume that the main limiting factor is the memory.
+       We'll assume that at least MEMAVAIL memory, or 4*n memory, whichever
+       is bigger, is available. */
+     mf = max(MEMAVAIL/n, 4);
+     if (stop->maxeval && stop->maxeval <= mf)
+         mf = max(stop->maxeval - 5, 1); /* mf > maxeval seems not good */
+
+ retry_alloc:
+     work = (double*) malloc(sizeof(double) * (n * 7 + max(n,n*mf) + 
+                                              max(n,mf)*2));
+     if (!work) { 
+         if (mf > 0) {
+              mf = 0; /* allocate minimal memory */
+              goto retry_alloc;
+         }
+         free(ix);
+         return NLOPT_OUT_OF_MEMORY;
+     }
+
+     xl = work; xu = xl + n;
+     gf = xu + n; s = gf + n; xo = s + n; go = xo + n; so = go + n;
+     xm = so + n;
+     xr = xm + max(n*mf,n); gr = xr + max(n,mf);
+
+     for (i = 0; i < n; ++i) {
+         int lbu = lb[i] <= -0.99 * HUGE_VAL; /* lb unbounded */
+         int ubu = ub[i] >= 0.99 * HUGE_VAL;  /* ub unbounded */
+         ix[i] = lbu ? (ubu ? 0 : 2) : (ubu ? 1 : (lb[i] == ub[i] ? 5 : 3));
+         xl[i] = lb[i];
+         xu[i] = ub[i];
+     }
+
+     /* ?  xo does not seem to be initialized in the
+       original Fortran code, but it is used upon
+       input to plip if mf > 0 ... perhaps ALLOCATE initializes
+       arrays to zero by default? */
+     memset(xo, 0, sizeof(double) * max(n,n*mf));
+
+     plip_(&n, &nb, x, ix, xl, xu, 
+          gf, s, xo, go, so, xm, xr, gr,
+          &xmax,
+
+          /* fixme: pass tol_rel and tol_abs and use NLopt check */
+          &stop->xtol_rel,
+          &stop->ftol_rel,
+          &stop->minf_max,
+          &tolg,
+          stop,
+
+          &minf_est, &gmax,
+          minf,
+          &mit, &mfv,
+          &iest,
+          &method, /* 1 == rank-one method VAR1, 2 == rank-two method VAR2 */
+          &mf,
+          &iterm, &stat,
+          f, f_data);
+
+     free(work);
+     free(ix);
+
+     switch (iterm) {
+        case 1: return NLOPT_XTOL_REACHED;
+        case 2: return NLOPT_FTOL_REACHED;
+        case 3: return NLOPT_MINF_MAX_REACHED;
+        case 4: return NLOPT_SUCCESS; /* gradient tolerance reached */
+        case 6: return NLOPT_SUCCESS;
+        case 12: case 13: return NLOPT_MAXEVAL_REACHED;
+        default: return NLOPT_FAILURE;
+     }
+}
index c31e959fa39dcb63f54ce23b5964f27abc65e9e3..3c6b71ad1b67c7a47834cb2ae1d0e0b375462f47 100644 (file)
@@ -415,11 +415,6 @@ L11190:
     return;
 } /* plis_ */
 
-/* number of double variables that can be stored in scratch memory
-   ... it's >= 2007, and this is in the context of scientific computation,
-   so assume that at least 10M are available, and that sizeof(double)==8 */
-#define MEMAVAIL 1310720
-
 /* NLopt wrapper around plis_, handling dynamic allocation etc. */
 nlopt_result luksan_plis(int n, nlopt_func f, void *f_data,
                  const double *lb, const double *ub, /* bounds */
@@ -427,7 +422,7 @@ nlopt_result luksan_plis(int n, nlopt_func f, void *f_data,
                  double *minf,
                  nlopt_stopping *stop)
 {
-     int i, *ix;
+     int i, *ix, nb = 1;
      double *work, *xl, *xu, *xo, *gf, *s, *go, *uo, *vo;
      double gmax, minf_est;
      double xmax = 0; /* no maximum */
@@ -458,7 +453,7 @@ nlopt_result luksan_plis(int n, nlopt_func f, void *f_data,
                                               max(n,mf)*2));
      if (!work) { 
          if (mf > 0) {
-              mf /= 4;
+              mf = 0; /* allocate minimal memory */
               goto retry_alloc;
          }
          free(ix);
@@ -483,7 +478,7 @@ nlopt_result luksan_plis(int n, nlopt_func f, void *f_data,
        arrays to zero by default? */
      memset(xo, 0, sizeof(double) * max(n,n*mf));
 
-     plis_(&n, &n, x, ix, xl, xu, 
+     plis_(&n, &nb, x, ix, xl, xu, 
           gf, s, xo, go, uo, vo,
           &xmax,