From: Ian Jackson <ijackson@chiark.greenend.org.uk>
Date: Sun, 8 Apr 2018 16:17:44 +0000 (+0100)
Subject: Revert "simplex wip: compiles"
X-Git-Url: https://www.chiark.greenend.org.uk/ucgi/~ian/git?a=commitdiff_plain;h=85aa911142e91ee9340473846da91d84509164c3;p=moebius3.git

Revert "simplex wip: compiles"

This reverts commit 657118ffcf7627d9c97d66065fba358248bbf860.
---

diff --git a/curveopt.py b/curveopt.py
index ab27e44..3367021 100644
--- a/curveopt.py
+++ b/curveopt.py
@@ -48,7 +48,7 @@ class OptimisedCurve():
 
     oc._dbg(repr(fc_input))
 
-    findcurve_epsilon = 0.0001
+    findcurve_epsilon = 0.01
 
     cl = ['./findcurve', '%d' % (nt+1), '%.18g' % findcurve_epsilon]
     oc._dbg('STARTING FINDCURVE %s' % cl)
diff --git a/findcurve.c b/findcurve.c
index 53abb47..055018f 100644
--- a/findcurve.c
+++ b/findcurve.c
@@ -3,13 +3,11 @@
 
 #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"
 
@@ -28,13 +26,6 @@ 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;
@@ -53,11 +44,9 @@ static void prepare(double X[] /* startpoint */) {
 
 //#define DEBUG
 
-static double cb_f(const gsl_vector *xg, void *params) {
+static double cb_Efunc(void *xp) {
+  const double *X = xp;
   int P;
-
-  GET_X(xg);
-
   DECLARE_F_G;
   CALCULATE_F_G;
 
@@ -76,7 +65,7 @@ static double cb_f(const gsl_vector *xg, void *params) {
   return cost;
 }
 
-static void __attribute__((unused))  cb_step(const gsl_rng *rng, void *xp, double step_size) {
+static void cb_step(const gsl_rng *rng, void *xp, double step_size) {
   double *x = xp;
   int i;
   double step[NX];
@@ -88,7 +77,7 @@ static void __attribute__((unused))  cb_step(const gsl_rng *rng, void *xp, doubl
   //printf("]\n");
 }
 
-static double __attribute__((unused))  cb_metric(void *xp, void *yp) {
+static double cb_metric(void *xp, void *yp) {
   const double *x=xp, *y=yp;
   int i;
   double s;
@@ -100,9 +89,9 @@ static double __attribute__((unused))  cb_metric(void *xp, void *yp) {
 }
 
 static void __attribute__((unused)) cb_print(void *xp) {
-  GET_X(xp);
+  const double *x = xp;
   printf("\n");
-  printcore(X);
+  printcore(x);
 }
 
 static double scan1double(void) {
@@ -118,48 +107,42 @@ static double scan1double(void) {
 
 int main(int argc, const char *const *argv) {
   double epsilon;
-  int i, r;
+  int i;
 
   NP = atoi(argv[1]);
   epsilon = atof(argv[2]);
 
-  double startpoint_X[NX];
-  double input[NINPUT]; INPUT = input;
-
-  gsl_multimin_function func = {
-    .f = cb_f,
-    .n = NX,
-    .params = 0,
-  };
+  gsl_rng *rng = gsl_rng_alloc(gsl_rng_ranlxd2);
 
-  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);
+  double input[NINPUT]; INPUT = input;
+  double startpoint[NX];
 
   for (;;) {
     /* NINPUT + 1 doubles: startpoint, epsilon for residual */
     for (i=0; i<NINPUT; i++)
       INPUT[i] = scan1double();
 
-    prepare(startpoint_X);
-    for (i=0; i<NX; i++)
-      gsl_vector_set(current_gx, i, startpoint_X[i]);
+    gsl_rng_set(rng,0);
 
-    r = gsl_multimin_fminimizer_set(minimiser, &func, current_gx, step_size);
-    assert(!r);
+    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,
+    };
 
-    for (;;) {
-      r = gsl_multimin_fminimizer_iterate(minimiser);
-      assert(!r);
+    prepare(startpoint);
 
-      double size = gsl_multimin_fminimizer_size(minimiser);
-      if (size < epsilon) break;
-    }      
+    gsl_siman_solve(rng,
+		    startpoint,
+		    cb_Efunc, cb_step, cb_metric,
+		    0, // cb_print,
+		    0,0,0, sizeof(startpoint), siman_params);
 
-    cb_print(startpoint_X);
+    printcore(startpoint);
 
     printf("[]\n");
     fflush(stdout);