update fit to use nlopt for minimization

The GSL library only has the Nelder Mead Simplex algorithm for doing
multidimensional minimization without gradient information. The nlopt library
has lots of different minimization algorithms so it's easier to switch between
them to see which one works best.

2 files changed, 52 insertions, 82 deletions

diff --git a/src/Makefile b/src/Makefile
index 2359e66..ae8807b 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,5 +1,5 @@
 CFLAGS=-O0 -Wall -g -DSWAP_BYTES
-LDLIBS=-lm -lgsl -lgslcblas
+LDLIBS=-lm -lgsl -lgslcblas -lnlopt_cxx -lstdc++
 
 all: test test-vector test-likelihood fit test-charge
 
diff --git a/src/fit.c b/src/fit.c
index 1f90f12..bfb5f49 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -4,39 +4,51 @@
 #include "Record_Info.h"
 #include "event.h"
 #include "zdab_utils.h"
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_multimin.h>
 #include "scattering.h"
 #include "pmt.h"
 #include "sno_charge.h"
 #include "db.h"
 #include "dqxx.h"
+#include <nlopt.h>
+#include <math.h> /* for sin(), cos(), etc. */
 
 #define EV_RECORD    0x45562020 // 'EV  '  (as written to ZDAB file)
 
-double nll(const gsl_vector *v, void *params)
+double nll(unsigned int n, const double *x, double *grad, void *params)
 {
     double T, theta, phi, t0;
     double pos[3], dir[3];
+    static size_t iter;
+    double fval;
 
-    T = gsl_vector_get(v,0);
+    T = x[0];
 
-    /* FIXME: should use a minimization algorithm with bounds. */
-    if (T < 0) return 1e9;
+    pos[0] = x[1];
+    pos[1] = x[2];
+    pos[2] = x[3];
 
-    pos[0] = gsl_vector_get(v,1);
-    pos[1] = gsl_vector_get(v,2);
-    pos[2] = gsl_vector_get(v,3);
-
-    theta = gsl_vector_get(v,4);
-    phi = gsl_vector_get(v,5);
+    theta = x[4];
+    phi = x[5];
     dir[0] = sin(theta)*cos(phi);
     dir[1] = sin(theta)*sin(phi);
     dir[2] = cos(theta);
 
-    t0 = gsl_vector_get(v,6);
+    t0 = x[6];
+
+    fval = nll_muon((event *) params, T, pos, dir, t0);
 
-    return nll_muon((event *) params, T, pos, dir, t0);
+    printf("%5zu %.2g %.2g %.2g %.2g %.2g %.2g %.2g f() = %7.3e\n",
+           iter++,
+           x[0],
+           x[1],
+           x[2],
+           x[3],
+           x[4],
+           x[5],
+           x[6],
+           fval);
+
+    return fval;
 }
 
 int main(int argc, char **argv)
@@ -50,6 +62,9 @@ int main(int argc, char **argv)
     event ev = {0};
     int crate, card, channel;
     int id;
+    double x[7];
+    double ss[7];
+    double opt_f;
 
     if (argc < 2) {
         fprintf(stderr, "usage: fit [filename]\n");
@@ -72,23 +87,8 @@ int main(int argc, char **argv)
 
     ev.run = -1;
 
-    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2;
-    gsl_multimin_fminimizer *s = NULL;
-    gsl_vector *ss, *x;
-    gsl_multimin_function minex_func;
-
-    size_t iter = 0;
-    int status;
-    double size;
-
-    /* Starting point */
-    x = gsl_vector_alloc(7);
-
-    ss = gsl_vector_alloc(7);
-
-    minex_func.n = 7;
-    minex_func.f = nll;
-    minex_func.params = &ev;
+    nlopt_opt opt = nlopt_create(NLOPT_LN_SBPLX, 7);
+    nlopt_set_min_objective(opt,nll,&ev);
 
     init_interpolation();
     init_charge();
@@ -132,54 +132,25 @@ int main(int argc, char **argv)
             if (ev.run != -1) {
                 if (ev.run != bev.run || ev.gtid != bev.gtr_id) {
                     /* New event, so we need to fit the old event. */
-                    iter = 0;
-
-                    gsl_vector_set(x,0,500.0);
-                    gsl_vector_set(x,1,100.0);
-                    gsl_vector_set(x,2,100.0);
-                    gsl_vector_set(x,3,100.0);
-                    gsl_vector_set(x,4,1.57);
-                    gsl_vector_set(x,5,0.0);
-                    gsl_vector_set(x,6,120.0);
-
-                    gsl_vector_set(ss,0,10.0);
-                    gsl_vector_set(ss,1,10.0);
-                    gsl_vector_set(ss,2,10.0);
-                    gsl_vector_set(ss,3,10.0);
-                    gsl_vector_set(ss,4,0.1);
-                    gsl_vector_set(ss,5,0.1);
-                    gsl_vector_set(ss,6,1.0);
-
-                    s = gsl_multimin_fminimizer_alloc(T, 7);
-                    gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
-
-                    do {
-                        iter++;
-                        status = gsl_multimin_fminimizer_iterate(s);
-
-                        if (status) break;
-
-                        size = gsl_multimin_fminimizer_size(s);
-                        status = gsl_multimin_test_size(size,1e-2);
-
-                        if (status == GSL_SUCCESS) {
-                            printf("converged to minimum at\n");
-                        }
-
-                        printf("%5zu %.2g %.2g %.2g %.2g %.2g %.2g %.2g f() = %7.3e size = %.3f\n",
-                               iter,
-                               gsl_vector_get(s->x,0),
-                               gsl_vector_get(s->x,1),
-                               gsl_vector_get(s->x,2),
-                               gsl_vector_get(s->x,3),
-                               gsl_vector_get(s->x,4),
-                               gsl_vector_get(s->x,5),
-                               gsl_vector_get(s->x,6),
-                               s->fval,
-                               size);
-                    } while (status == GSL_CONTINUE && iter < 1000);
-
-                    gsl_multimin_fminimizer_free(s);
+                    x[0] = 500.0;
+                    x[1] = 100.0;
+                    x[2] = 100.0;
+                    x[3] = 100.0;
+                    x[4] = 1.57;
+                    x[5] = 0.0;
+                    x[6] = 120.0;
+
+                    ss[0] = 10.0;
+                    ss[1] = 10.0;
+                    ss[2] = 10.0;
+                    ss[3] = 10.0;
+                    ss[4] = 0.1;
+                    ss[5] = 0.1;
+                    ss[6] = 1.0;
+
+                    nlopt_set_initial_step(opt, ss);
+
+                    nlopt_optimize(opt,x,&opt_f);
                 }
             }
             ev.run = bev.run;
@@ -191,8 +162,7 @@ int main(int argc, char **argv)
 
     db_free(db);
 
-    gsl_vector_free(x);
-    gsl_vector_free(ss);
+    nlopt_destroy(opt);
 
     zebra_close(f);