start by doing a series of "quick" minimizations

This commit updates the fit_event function to first do a series of "quick"
minimizations to try and find a good set of starting parameters for the "real"
fit.

I also updated the bounds on theta and phi since having hard bounds on these
angle coordinates could prevent the fitter from finding the minimum.

1 files changed, 63 insertions, 9 deletions

diff --git a/src/fit.c b/src/fit.c
index 87bf232..045141d 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -27,6 +27,26 @@ typedef struct fitParams {
     double epsrel;
 } fitParams;
 
+/* In order to start the fitter close to the minimum, we first do a series of
+ * "quick" minimizations starting at the following points. We keep track of the
+ * parameters with the best likelihood value and then start the "real"
+ * minimization from those parameters. */
+
+static struct startingParameters {
+    double x;
+    double y;
+    double z;
+    double theta;
+    double phi;
+} startingParameters[] = {
+    {0.0, 0.0, 0.0,    0.0, 0.0},
+    {0.0, 0.0, 0.0, M_PI/2, 0.0},
+    {0.0, 0.0, 0.0, M_PI/2, M_PI/2},
+    {0.0, 0.0, 0.0, M_PI/2, M_PI},
+    {0.0, 0.0, 0.0, M_PI/2, 3*M_PI/2},
+    {0.0, 0.0, 0.0,   M_PI, 0.0},
+};
+
 double nll(unsigned int n, const double *x, double *grad, void *params)
 {
     fitParams *fpars = (fitParams *) params;
@@ -59,7 +79,7 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
 
     long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
 
-    printf("%5zu %10.2f %7.2f %7.2f %7.2f %4.2f %4.2f %5.2f %5.2f %5.2f f() = %7.3e took %lld ms\n",
+    printf("%5zu %10.2f %7.2f %7.2f %7.2f %5.2f %5.2f %5.2f %5.2f %5.2f f() = %7.3e took %lld ms\n",
            iter++,
            x[0],
            x[1],
@@ -78,13 +98,13 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
 
 int fit_event(event *ev, double *xopt, double *fmin)
 {
+    size_t i;
     fitParams fpars;
     double x[9], ss[9], lb[9], ub[9], fval, n;
     int rv;
 
     nlopt_opt opt = nlopt_create(NLOPT_LN_BOBYQA, 9);
     nlopt_set_min_objective(opt,nll,&fpars);
-    nlopt_set_ftol_abs(opt, 1e-5);
 
     x[0] = 1200.0;
     x[1] = 0.0;
@@ -112,8 +132,8 @@ int fit_event(event *ev, double *xopt, double *fmin)
     lb[1] = -1000.0;
     lb[2] = -1000.0;
     lb[3] = -1000.0;
-    lb[4] = 0.0;
-    lb[5] = 0.0;
+    lb[4] = -INFINITY;
+    lb[5] = -INFINITY;
     lb[6] = 0.0;
     lb[7] = -10.0;
     lb[8] = -10.0;
@@ -122,8 +142,8 @@ int fit_event(event *ev, double *xopt, double *fmin)
     ub[1] = 1000.0;
     ub[2] = 1000.0;
     ub[3] = 1000.0;
-    ub[4] = M_PI;
-    ub[5] = 2*M_PI;
+    ub[4] = INFINITY;
+    ub[5] = INFINITY;
     ub[6] = 400.0;
     ub[7] = 10.0;
     ub[8] = 10.0;
@@ -134,10 +154,44 @@ int fit_event(event *ev, double *xopt, double *fmin)
     nlopt_set_initial_step(opt, ss);
 
     fpars.ev = ev;
-    fpars.epsrel = 1e-1;
 
     iter = 0;
-    rv = nlopt_optimize(opt,x,&fval);
+
+    /* First we do a set of "quick" minimizations to try and start the
+     * minimizer close to the minimum. To make these function evaluations
+     * faster, we set the absolute tolerance on the likelihood to 1.0, the
+     * maximum number of function evaluations to 100, and the relative
+     * tolerance on the numerical integration to 10%. */
+    fpars.epsrel = 1e-1;
+    nlopt_set_ftol_abs(opt, 1.0);
+    nlopt_set_maxeval(opt, 100);
+
+    for (i = 0; i < sizeof(startingParameters)/sizeof(startingParameters[0]); i++) {
+        x[0] = 1200.0;
+        x[1] = startingParameters[i].x;
+        x[2] = startingParameters[i].y;
+        x[3] = startingParameters[i].z;
+        x[4] = startingParameters[i].theta;
+        x[5] = startingParameters[i].phi;
+        x[6] = 130.0;
+        x[7] = 0.0;
+        x[8] = 0.0;
+
+        rv = nlopt_optimize(opt,x,&fval);
+
+        if (fval < *fmin || i == 0) {
+            *fmin = fval;
+            memcpy(xopt,x,sizeof(x));
+        }
+    }
+
+    memcpy(x,xopt,sizeof(x));
+
+    /* Now, we do the "real" minimization. */
+    fpars.epsrel = 1e-4;
+    nlopt_set_ftol_abs(opt, 1e-5);
+    nlopt_set_maxeval(opt, 1000);
+
 
     do {
         *fmin = fval;
@@ -270,7 +324,7 @@ int main(int argc, char **argv)
                     fit_event(&ev,xopt,&fmin);
 
                     if (fout) {
-                        fprintf(fout, "%" PRIu32 " %10.2f %7.2f %7.2f %7.2f %4.2f %4.2f %5.2f %5.2f %5.2f f() = %7.3e\n",
+                        fprintf(fout, "%" PRIu32 " %10.2f %7.2f %7.2f %7.2f %5.2f %5.2f %5.2f %5.2f %5.2f f() = %7.3e\n",
                                 ev.gtid,
                                 xopt[0],
                                 xopt[1],