add a function to return the kahan sum of an array

For some reason the fit seems to have trouble with the kinetic energy.
Basically, it seems to "converge" even though when you run the minimization
again it finds a better minimum with a lower energy. I think this is likely due
to the fact that for muons the kinetic energy only really affects the range of
the muon and this is subject to error in the numerical integration.

I also thought that maybe it could be due to roundoff error in the likelihood
calculation, so I implemented the Kahan summation to try and reduce that. No
idea if it's actually improving things, but I should benchmark it later to see.

4 files changed, 56 insertions, 7 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 66e72e4..012243c 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -124,10 +124,10 @@ double getKineticEnergy(double x, double T0)
 
 double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel)
 {
-    size_t i, j;
+    size_t i, j, nhit;
     intParams params;
     double total_charge;
-    double logp[MAX_PE], nll, range, pmt_dir[3], R, x, cos_theta, theta, theta_cerenkov, theta0, E0, p0, beta0;
+    double logp[MAX_PE], nll[MAX_PMTS], range, pmt_dir[3], R, x, cos_theta, theta, theta_cerenkov, theta0, E0, p0, beta0;
     double tmean = 0.0;
     int npmt = 0;
 
@@ -236,7 +236,7 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
         mu[i] = mu_direct[i] + mu_indirect + mu_noise;
     }
 
-    nll = 0;
+    nhit = 0;
     for (i = 0; i < MAX_PMTS; i++) {
         if (ev->pmt_hits[i].flags || (pmts[i].pmt_type != PMT_NORMAL && pmts[i].pmt_type != PMT_OWL)) continue;
 
@@ -245,15 +245,15 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
             for (j = 1; j < MAX_PE; j++) {
                 logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log(mu[i]) - gsl_sf_lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], 1, &ts[i], tmean, 1.5);
             }
-            nll -= logsumexp(logp, sizeof(logp)/sizeof(double));
+            nll[nhit++] = -logsumexp(logp, sizeof(logp)/sizeof(double));
         } else {
             logp[0] = -mu[i];
             for (j = 1; j < MAX_PE; j++) {
                 logp[j] = log(get_pmiss(j)) - mu[i] + j*log(mu[i]) - gsl_sf_lnfact(j);
             }
-            nll -= logsumexp(logp, sizeof(logp)/sizeof(double));
+            nll[nhit++] = -logsumexp(logp, sizeof(logp)/sizeof(double));
         }
     }
 
-    return nll;
+    return kahan_sum(nll,nhit);
 }
diff --git a/src/misc.c b/src/misc.c
index 2f02595..bff0ba9 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -2,6 +2,26 @@
 #include <math.h>
 #include <stdlib.h> /* for size_t */
 
+double kahan_sum(double *x, size_t n)
+{
+    /* Returns the sum of the elements of `x` using the Kahan summation algorithm.
+     *
+     * See https://en.wikipedia.org/wiki/Kahan_summation_algorithm. */
+    size_t i;
+    double sum, c, y, t;
+
+    sum = 0.0;
+    c = 0.0;
+    for (i = 0; i < n; i++) {
+        y = x[i] - c;
+        t = sum + y;
+        c = (t - sum) - y;
+        sum = t;
+    }
+
+    return sum;
+}
+
 double interp1d(double x, double *xp, double *yp, size_t n)
 {
     /* A fast interpolation routine which assumes that the values in `xp` are
diff --git a/src/misc.h b/src/misc.h
index c1897d9..c93afba 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -3,6 +3,7 @@
 
 #include <stdlib.h> /* for size_t */
 
+double kahan_sum(double *x, size_t n);
 double interp1d(double x, double *xp, double *yp, size_t n);
 int isclose(double a, double b, double rel_tol, double abs_tol);
 int allclose(double *a, double *b, size_t n, double rel_tol, double abs_tol);
diff --git a/src/test.c b/src/test.c
index 2ea73d1..fb788ba 100644
--- a/src/test.c
+++ b/src/test.c
@@ -440,6 +440,33 @@ err:
     return 1;
 }
 
+int test_kahan_sum(char *err)
+{
+    /* Tests the kahan_sum() function. */
+    size_t i;
+    double x[100], sum, expected;
+
+    init_genrand(0);
+
+    expected = 0.0;
+    for (i = 0; i < sizeof(x)/sizeof(x[0]); i++) {
+        x[i] = genrand_real2();
+        expected += x[i];
+    }
+
+    sum = kahan_sum(x,sizeof(x)/sizeof(x[0]));
+
+    if (!isclose(sum, expected, 1e-9, 1e-9)) {
+        sprintf(err, "kahan_sum returned %.5g, but expected %.5g", sum, expected);
+        goto err;
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -453,7 +480,8 @@ struct tests {
     {test_logsumexp, "test_logsumexp"},
     {test_sno_charge_integral, "test_sno_charge_integral"},
     {test_path, "test_path"},
-    {test_interp1d, "test_interp1d"}
+    {test_interp1d, "test_interp1d"},
+    {test_kahan_sum, "test_kahan_sum"},
 };
 
 int main(int argc, char **argv)