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, 6 insertions, 0 deletions

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);