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.

1 files changed, 31 insertions, 0 deletions

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)