diff options
| -rw-r--r-- | src/likelihood.c | 12 | ||||
| -rw-r--r-- | src/misc.c | 20 | ||||
| -rw-r--r-- | src/misc.h | 1 | ||||
| -rw-r--r-- | src/test.c | 30 |
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); } @@ -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 @@ -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); @@ -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) |