From: Ian Jackson <ijackson@chiark.greenend.org.uk>
Date: Sun, 8 Apr 2018 15:10:33 +0000 (+0100)
Subject: simplex wip: compiles
X-Git-Url: https://www.chiark.greenend.org.uk/ucgi/~ian/git?a=commitdiff_plain;h=657118ffcf7627d9c97d66065fba358248bbf860;p=moebius3.git

simplex wip: compiles

Signed-off-by: Ian Jackson <ijackson@chiark.greenend.org.uk>
---

diff --git a/curveopt.py b/curveopt.py
index 8825da5..ee14c5d 100644
--- a/curveopt.py
+++ b/curveopt.py
@@ -48,7 +48,7 @@ class OptimisedCurve():
 
     oc._dbg(repr(fc_input))
 
-    findcurve_epsilon = 0.01
+    findcurve_epsilon = 0.0001
 
     cl = ['./findcurve', '%d' % (nt+1), '%.18g' % findcurve_epsilon]
     oc._dbg('STARTING FINDCURVE %s' % cl)
diff --git a/findcurve.c b/findcurve.c
index 055018f..53abb47 100644
--- a/findcurve.c
+++ b/findcurve.c
@@ -3,11 +3,13 @@
 
 #include <math.h>
 #include <string.h>
+#include <assert.h>
 
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_sf.h>
 #include <gsl/gsl_siman.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_multimin.h>
 
 #include "symbolic.c"
 
@@ -26,6 +28,13 @@ static int NP;
 static double *INPUT; /* dyanmic array, on main's stack */
 static double PREP[NPREP];
 
+#define GET_X(xg)					\
+  double X[NX];						\
+  ({ int get_x_i;					\
+     for (get_x_i=0; get_x_i<NX; get_x_i++)		\
+       X[get_x_i] = gsl_vector_get(xg, get_x_i);	\
+  })
+
 static void printcore(const double *X) {
   int i, j;
   DECLARE_F_G;
@@ -44,9 +53,11 @@ static void prepare(double X[] /* startpoint */) {
 
 //#define DEBUG
 
-static double cb_Efunc(void *xp) {
-  const double *X = xp;
+static double cb_f(const gsl_vector *xg, void *params) {
   int P;
+
+  GET_X(xg);
+
   DECLARE_F_G;
   CALCULATE_F_G;
 
@@ -65,7 +76,7 @@ static double cb_Efunc(void *xp) {
   return cost;
 }
 
-static void cb_step(const gsl_rng *rng, void *xp, double step_size) {
+static void __attribute__((unused))  cb_step(const gsl_rng *rng, void *xp, double step_size) {
   double *x = xp;
   int i;
   double step[NX];
@@ -77,7 +88,7 @@ static void cb_step(const gsl_rng *rng, void *xp, double step_size) {
   //printf("]\n");
 }
 
-static double cb_metric(void *xp, void *yp) {
+static double __attribute__((unused))  cb_metric(void *xp, void *yp) {
   const double *x=xp, *y=yp;
   int i;
   double s;
@@ -89,9 +100,9 @@ static double cb_metric(void *xp, void *yp) {
 }
 
 static void __attribute__((unused)) cb_print(void *xp) {
-  const double *x = xp;
+  GET_X(xp);
   printf("\n");
-  printcore(x);
+  printcore(X);
 }
 
 static double scan1double(void) {
@@ -107,42 +118,48 @@ static double scan1double(void) {
 
 int main(int argc, const char *const *argv) {
   double epsilon;
-  int i;
+  int i, r;
 
   NP = atoi(argv[1]);
   epsilon = atof(argv[2]);
 
-  gsl_rng *rng = gsl_rng_alloc(gsl_rng_ranlxd2);
-
+  double startpoint_X[NX];
   double input[NINPUT]; INPUT = input;
-  double startpoint[NX];
+
+  gsl_multimin_function func = {
+    .f = cb_f,
+    .n = NX,
+    .params = 0,
+  };
+
+  gsl_vector *current_gx = gsl_vector_alloc(NX);
+  gsl_vector *step_size = gsl_vector_alloc(NX);
+  gsl_vector_set_all(step_size, 0.01);
+
+  gsl_multimin_fminimizer *minimiser =
+    gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex2, NX);
 
   for (;;) {
     /* NINPUT + 1 doubles: startpoint, epsilon for residual */
     for (i=0; i<NINPUT; i++)
       INPUT[i] = scan1double();
 
-    gsl_rng_set(rng,0);
+    prepare(startpoint_X);
+    for (i=0; i<NX; i++)
+      gsl_vector_set(current_gx, i, startpoint_X[i]);
 
-    gsl_siman_params_t siman_params = {
-      .k = 1.0,
-      .t_initial = 0.5,
-      .mu_t = 1.001,
-      .t_min = epsilon * 1E-3,
-      .iters_fixed_T = 100,
-      .n_tries = 10,
-      .step_size = 0.05,
-    };
+    r = gsl_multimin_fminimizer_set(minimiser, &func, current_gx, step_size);
+    assert(!r);
 
-    prepare(startpoint);
+    for (;;) {
+      r = gsl_multimin_fminimizer_iterate(minimiser);
+      assert(!r);
 
-    gsl_siman_solve(rng,
-		    startpoint,
-		    cb_Efunc, cb_step, cb_metric,
-		    0, // cb_print,
-		    0,0,0, sizeof(startpoint), siman_params);
+      double size = gsl_multimin_fminimizer_size(minimiser);
+      if (size < epsilon) break;
+    }      
 
-    printcore(startpoint);
+    cb_print(startpoint_X);
 
     printf("[]\n");
     fflush(stdout);