update path integral to use a fixed number of points

I noticed when fitting electrons that the cquad integration routine was not
very stable, i.e. it would return different results for *very* small changes in
the fit parameters which would cause the fit to stall.

Since it's very important for the minimizer that the likelihood function not
jump around, I am switching to integrating over the path by just using a fixed
number of points and using the trapezoidal rule. This seems to be a lot more
stable, and as a bonus I was able to combine the three integrals (direct
charge, indirect charge, and time) so that we only have to do a single loop.
This should hopefully make the speed comparable since the cquad routine was
fairly effective at only using as many function evaluations as needed.

Another benefit to this approach is that if needed, it will be easier to port
to a GPU.

6 files changed, 98 insertions, 74 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 1aa1653..ec58d11 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -22,11 +22,6 @@
 #include "proton.h"
 #include "id_particles.h"
 
-typedef struct intParams {
-    path *p;
-    int i;
-} intParams;
-
 particle *particle_init(int id, double T0, size_t n)
 {
     /* Returns a struct with arrays of the particle position and kinetic
@@ -160,9 +155,9 @@ void particle_free(particle *p)
     free(p);
 }
 
-double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, int reflected)
+static double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o)
 {
-    double pmt_dir[3], cos_theta, n, wavelength0, omega, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic;
+    double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge;
 
     z = 1.0;
 
@@ -182,31 +177,29 @@ double get_expected_charge(double x, double beta, double theta0, double *pos, do
 
     R = NORM(pos);
 
-    /* FIXME: I just calculate delta assuming 400 nm light. */
-    wavelength0 = 400.0;
-
     if (R <= AV_RADIUS) {
-        n = get_index_snoman_d2o(wavelength0);
+        n = n_d2o;
     } else {
-        n = get_index_snoman_h2o(wavelength0);
+        n = n_h2o;
     }
 
     if (beta < 1/n) return 0;
 
-    if (reflected)
-        f = get_weighted_pmt_reflectivity(acos(-cos_theta_pmt));
-    else
-        f = get_weighted_pmt_response(acos(-cos_theta_pmt));
+    theta_pmt = acos(-cos_theta_pmt);
+    f_reflec = get_weighted_pmt_reflectivity(theta_pmt);
+    f = get_weighted_pmt_response(theta_pmt);
 
     absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
     absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();
     absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic();
 
-    get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
-
     l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;
 
-    return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
+    charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
+
+    *reflected = f_reflec*charge;
+
+    return f*charge;
 }
 
 double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double sigma)
@@ -271,43 +264,41 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
     return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma));
 }
 
-static double gsl_time(double x, void *params)
+static void integrate_path(path *p, int pmt, double a, double b, size_t n, double *mu_direct, double *mu_indirect, double *time)
 {
-    intParams *pars = (intParams *) params;
-    double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o;
-
-    path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
+    size_t i;
+    double *qs, *q_indirects, *ts;
+    double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o, x;
 
-    get_path_length(pos,pmts[pars->i].pos,AV_RADIUS,&l_d2o,&l_h2o);
+    qs = malloc(sizeof(double)*n);
+    q_indirects = malloc(sizeof(double)*n);
+    ts = malloc(sizeof(double)*n);
 
     /* FIXME: I just calculate delta assuming 400 nm light. */
     wavelength0 = 400.0;
     n_d2o = get_index_snoman_d2o(wavelength0);
     n_h2o = get_index_snoman_h2o(wavelength0);
 
-    t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
+    for (i = 0; i < n; i++) {
+        x = a + i*(b-a)/(n-1);
 
-    return t*get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 0);
-}
+        path_eval(p, x, pos, dir, &beta, &t, &theta0);
 
-static double gsl_reflected_charge(double x, void *params)
-{
-    intParams *pars = (intParams *) params;
-    double dir[3], pos[3], t, theta0, beta;
+        get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
 
-    path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
+        t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
 
-    return get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 1);
-}
-
-static double gsl_charge(double x, void *params)
-{
-    intParams *pars = (intParams *) params;
-    double dir[3], pos[3], t, theta0, beta;
+        qs[i] = get_expected_charge(x, beta, theta0, pos, dir, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o);
+        ts[i] = t*qs[i];
+    }
 
-    path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
+    *mu_direct = trapz(qs,(b-a)/(n-1),n);
+    *mu_indirect = trapz(q_indirects,(b-a)/(n-1),n);
+    *time = trapz(ts,(b-a)/(n-1),n);
 
-    return get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 0);
+    free(qs);
+    free(q_indirects);
+    free(ts);
 }
 
 double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected)
@@ -532,7 +523,6 @@ static double getKineticEnergy(double x, void *p)
 double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast)
 {
     size_t i, j, nhit;
-    intParams params;
     double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
     particle *p;
     double pmt_dir[3], R, cos_theta, theta, wavelength0, n_d2o, n_h2o, theta_cerenkov, s, l_h2o, l_d2o;
@@ -540,19 +530,14 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
     double a, b;
     double tmp[3];
     double mu_reflected;
+    path *path;
 
     double mu_direct[MAX_PMTS];
     double ts[MAX_PMTS];
     double mu[MAX_PMTS];
     double mu_noise, mu_indirect;
 
-    gsl_integration_cquad_workspace *w = gsl_integration_cquad_workspace_alloc(100);
-    double result, error;
-
-    size_t nevals;
-
-    gsl_function F;
-    F.params = &params;
+    double result;
 
     p = particle_init(id, T0, 10000);
     range = p->range;
@@ -565,7 +550,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
     /* FIXME: is this formula valid for muons? */
     theta0 = get_scattering_rms(range/2,p0,beta0,1.0)/sqrt(range/2);
 
-    params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, p, z1, z2, n, p->mass);
+    path = path_init(pos, dir, T0, range, theta0, getKineticEnergy, p, z1, z2, n, p->mass);
 
     if (beta0 > BETA_MIN)
         get_smax(p, BETA_MIN, range, &smax);
@@ -580,8 +565,6 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
          * calculation. */
         if (fast && !ev->pmt_hits[i].hit) continue;
 
-        params.i = i;
-
         if (fast) {
             mu_direct[i] = get_total_charge_approx(T0, pos, dir, p, i, smax, theta0, &ts[i], &mu_reflected);
             ts[i] += t0;
@@ -640,20 +623,11 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
             if (a < 0.0) a = 0.0;
             if (b > range) b = range;
 
-            F.function = &gsl_charge;
-
-            gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
-            mu_direct[i] = result;
-
-            F.function = &gsl_reflected_charge;
-
-            gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
-            mu_indirect += result;
-
-            F.function = &gsl_time;
+            double q_indirect;
+            integrate_path(path, i, a, b, 100, mu_direct+i, &q_indirect, &result);
+            mu_indirect += q_indirect;
 
             if (mu_direct[i] > 1e-9) {
-                gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
                 ts[i] = t0 + result/mu_direct[i];
             } else {
                 /* If the expected number of PE is very small, our estimate of
@@ -678,12 +652,10 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
         }
     }
 
-    path_free(params.p);
+    path_free(path);
 
     particle_free(p);
 
-    gsl_integration_cquad_workspace_free(w);
-
     mu_noise = DARK_RATE*GTVALID*1e-9;
     mu_indirect *= CHARGE_FRACTION/10000.0;
 
diff --git a/src/likelihood.h b/src/likelihood.h
index 71f77f0..8328322 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -44,7 +44,6 @@ typedef struct particle {
 particle *particle_init(int id, double T0, size_t n);
 double particle_get_energy(double x, particle *p);
 void particle_free(particle *p);
-double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, int reflected);
 double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast);
 
 #endif
diff --git a/src/misc.c b/src/misc.c
index ba0a090..d20bc26 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -218,6 +218,26 @@ static struct {
     {100,363.73937555556347},
 };
 
+double trapz(const double *y, double dx, size_t n)
+{
+    /* Returns the integral of `y` using the trapezoidal rule assuming a
+     * constant grid spacing `dx`.
+     *
+     * See https://en.wikipedia.org/wiki/Trapezoidal_rule. */
+    size_t i;
+    double sum = 0.0;
+
+    if (n < 2) return 0.0;
+
+    sum = y[0];
+    for (i = 1; i < n-1; i++) {
+        sum += 2*y[i];
+    }
+    sum += y[n-1];
+
+    return sum*dx/2.0;
+}
+
 void get_path_length(double *pos1, double *pos2, double R, double *l1, double *l2)
 {
     /* Returns the path length inside and outside a circle of radius `R` for a
diff --git a/src/misc.h b/src/misc.h
index dadbbce..80a580f 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -3,9 +3,15 @@
 
 #include <stdlib.h> /* for size_t */
 
+/* Macro to compute the size of a static C array.
+ *
+ * See https://stackoverflow.com/questions/1598773. */
+#define LEN(x) ((sizeof(x)/sizeof(0[x]))/((size_t)(!(sizeof(x) % sizeof(0[x])))))
+
 #define LN_MAX 100
 #define LNFACT_MAX 100
 
+double trapz(const double *y, double dx, size_t n);
 void get_path_length(double *pos1, double *pos2, double R, double *l1, double *l2);
 double ln(unsigned int n);
 double lnfact(unsigned int n);
diff --git a/src/optics.c b/src/optics.c
index 1286490..2d28ba4 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -6,11 +6,7 @@
 #include <gsl/gsl_integration.h>
 #include <gsl/gsl_errno.h> /* for gsl_strerror() */
 #include <gsl/gsl_spline.h>
-
-/* Macro to compute the size of a static C array.
- *
- * See https://stackoverflow.com/questions/1598773. */
-#define LEN(x) ((sizeof(x)/sizeof(0[x]))/((size_t)(!(sizeof(x) % sizeof(0[x])))))
+#include "misc.h"
 
 /* Absorption coefficients for H2O and D2O as a function of wavelength from
  * SNOMAN.
diff --git a/src/test.c b/src/test.c
index c1da272..242e728 100644
--- a/src/test.c
+++ b/src/test.c
@@ -879,6 +879,36 @@ err:
     return 1;
 }
 
+int test_trapz(char *err)
+{
+    /* Tests the trapz() function. */
+    size_t i;
+    double y[100], integral, expected;
+
+    init_genrand(0);
+
+    for (i = 0; i < LEN(y); i++) {
+        y[i] = genrand_real2();
+    }
+
+    expected = 0.0;
+    for (i = 1; i < LEN(y); i++) {
+        expected += (y[i-1] + y[i])/2;
+    }
+
+    integral = trapz(y, 1.0, LEN(y));
+
+    if (!isclose(integral, expected, 0, 1e-9)) {
+        sprintf(err, "trapz() returned %.5g, but expected %.5g", integral, expected);
+        goto err;
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -906,6 +936,7 @@ struct tests {
     {test_proton_get_range, "test_proton_get_range"},
     {test_proton_get_energy, "test_proton_get_energy"},
     {test_log_norm, "test_log_norm"},
+    {test_trapz, "test_trapz"},
 };
 
 int main(int argc, char **argv)