From: Ian Jackson <ijackson@chiark.greenend.org.uk>
Date: Sun, 19 Nov 2017 17:10:39 +0000 (+0000)
Subject: findcurve: wip
X-Git-Url: https://www.chiark.greenend.org.uk/ucgi/~ian/git?a=commitdiff_plain;h=ea203f95f6809ca97fa2dfbbd6b79b95634a9486;p=moebius3.git

findcurve: wip

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

diff --git a/.gitignore b/.gitignore
index f9f89d2..e1cabef 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,3 +5,4 @@ moebius-mesh.scad
 commitid.scad
 *.tmp
 *.auto.scad
+symbolic.c
diff --git a/Makefile b/Makefile
index 0cd2b89..1a5948f 100644
--- a/Makefile
+++ b/Makefile
@@ -16,6 +16,11 @@ moebius-mesh.scad: meshscad $(PYLIBS)
 
 AUTO_TOPLEVELS := $(foreach m,$(USING_AUTOS),$(shell $(PLAY)/toplevel-find $m))
 
+symbolic.c:	symbolic.py
+		./symbolic.py -q >$@.tmp && mv -f $@.tmp $@
+
+findcurve.o:	findcurve.c symbolic.c
+
 autoincs:	$(AUTO_INCS)
 scads:		$(addsuffix .auto.scad, $(AUTO_TOPLEVELS))
 stls:		$(addsuffix .auto.stl, $(AUTO_TOPLEVELS))
diff --git a/findcurve.c b/findcurve.c
new file mode 100644
index 0000000..6ba51ef
--- /dev/null
+++ b/findcurve.c
@@ -0,0 +1,85 @@
+/*
+ */
+
+#include <gsl/gsl_multiroots.h>
+
+#include "symbolic.c"
+
+#define X(i) gsl_vector_get(x,i)
+#define J_END_COL(i) \
+  for (j=0; j<N; j++) gsl_matrix_set(f,i,j,J_COL[j]);
+
+static int cb_f(const gsl_vector *x, void *params, gsl_vector *f) {
+  double *target = params;
+  double F[6];
+
+  int i;
+  X_EXTRACT;
+  F_POPULATE;
+  for (i=0; i<N; i++) gsl_vector_set(f,i, F(i) - target[i]);
+  return 0;
+}
+
+static int cb_df(const gsl_vector *x, void *params, gsl_matrix *J) {
+  double J_COL[N];
+  int j;
+  X_EXTRACT;
+  J_POPULATE;
+}
+
+static int cb_fdf(const gsl_vector *x, void *params,
+		  gsl_vector *f, gsl_matrix *J) {
+  cb_f(x,params,f);
+  cb_df(x,params,df);
+}
+
+static double *scan1double(void) {
+  double v;
+  int r;
+
+  r = scanf("%lf",&v);
+  if (ferror(stdin)) { perror("stdin"); exit(-1); }
+  if (feof(stdin)) exit(0);
+  if (r!=1) { fputs("bad input\n",stderr); exit(-1); }
+  return v;
+}
+
+static int main(void) {
+  double target[N], epsilon;
+
+  gsl_multiroot_function_fdf function;
+  function.f = cb_f;
+  function.df = cb_df;
+  function.fdf = cb_fdf;
+  function.n = N;
+  function.params = target;
+
+  gsl_vector *startpoint = gsl_vector_alloc(N);
+
+  gsl_multiroot_fdfsolver *solver
+    = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_hybrid, N);
+
+  for (;;) {
+    /* 2N+1 doubles: target, initial guess, epsilon for residual */
+    for (i=0; i<N; i++)
+      target[i] = scan1double();
+    for (i=0; i<N; i++)
+      gsl_vector_set(startpoint,i,scan1double());
+    epsilon = scan1double();
+
+    gsl_multiroot_fdfsolver_set(solver, function, startpoint);
+
+    for (;;) {
+      r = gsl_multiroot_fdfsolver_iterate(solver);
+      if (r) break;
+
+      r = gsl_multiroot_test_residual(solver->f, epsilon);
+      if (r != GSL_CONTINUE) break;
+    }
+
+    if (r) {
+      fprintf(stderr,"ERROR %s\n",gsl_strerror(gsl_errno));
+      exit(-1);
+    }
+  }
+}
diff --git a/symbolic.py b/symbolic.py
index a7dc79e..ecd94b5 100755
--- a/symbolic.py
+++ b/symbolic.py
@@ -204,7 +204,7 @@ def gen_diff(current, smalls):
     dbg('j')
     j = cse_prep_cprint(j, 'jtmp')
     for ix in range(0, j.cols):
-      cprint(ourccode(j.col(ix), 'J(%d)' % ix))
+      cprint(ourccode(j.col(ix), 'J_COL'))
       cprint('J_END_COL(%d)' % ix)
   else:
     small = smalls[0]
@@ -220,10 +220,14 @@ def gen_diff(current, smalls):
     )
     cprint('} /* %s small */' % small)
 
+def gen_misc():
+  cprintraw('// AUTOGENERATED - DO NOT EDIT\n')
+  cprintraw('#define N %d\n' % len(params))
+
 def gen_x_extract():
-  cprint('#define X_EXTRACT \\')
+  cprint('#define X_EXTRACT')
   for ix in range(0, len(params)):
-    cprint('double %s = X(%d)' % (params[ix], ix))
+    cprint('double %s = X(%d);' % (params[ix], ix))
   cprintraw()
 
 def gen_f_populate():
@@ -236,6 +240,10 @@ def gen_j_populate():
   gen_diff(result_dirnscaled, (sh*sh*la, th*th*la))
   cprintraw('')
 
-gen_x_extract()
-gen_f_populate()
-gen_j_populate()
+def gen_C():
+  gen_misc()
+  gen_x_extract()
+  gen_f_populate()
+  gen_j_populate()
+
+gen_C()