From 83119fe9bb0c9c2a70ea5397aca9968bde7ffa07 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Fri, 31 Aug 2018 10:33:24 -0500 Subject: add interp1d function to do fast interpolation when the x values are evenly spaced --- src/misc.c | 16 ++++++++++++++++ src/misc.h | 1 + src/test.c | 47 ++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 63 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/misc.c b/src/misc.c index 6a67d58..2f02595 100644 --- a/src/misc.c +++ b/src/misc.c @@ -2,6 +2,22 @@ #include #include /* for size_t */ +double interp1d(double x, double *xp, double *yp, size_t n) +{ + /* A fast interpolation routine which assumes that the values in `xp` are + * evenly spaced. + * + * If x < xp[0] returns yp[0] and if x > xp[n-1] returns yp[n-1]. */ + size_t i; + + if (x < xp[0]) return yp[0]; + if (x > xp[n-1]) return yp[n-1]; + + i = (x-xp[0])/(xp[1]-xp[0]); + + return yp[i] + (yp[i+1]-yp[i])*(x-xp[i])/(xp[i+1]-xp[i]); +} + int isclose(double a, double b, double rel_tol, double abs_tol) { /* Returns 1 if a and b are "close". This algorithm is taken from Python's diff --git a/src/misc.h b/src/misc.h index 295e59c..c1897d9 100644 --- a/src/misc.h +++ b/src/misc.h @@ -3,6 +3,7 @@ #include /* for size_t */ +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); double logsumexp(double *a, size_t n); diff --git a/src/test.c b/src/test.c index cc054c0..2ea73d1 100644 --- a/src/test.c +++ b/src/test.c @@ -11,6 +11,7 @@ #include "path.h" #include "random.h" #include "pdg.h" +#include typedef int testFunction(char *err); @@ -396,6 +397,49 @@ int test_path(char *err) return 0; } +int test_interp1d(char *err) +{ + /* Tests the interp1d() function. */ + size_t i; + double xp[100], yp[100], range, x, y, expected; + gsl_interp_accel *acc; + gsl_spline *spline; + + range = 1.0; + + init_genrand(0); + + for (i = 0; i < sizeof(xp)/sizeof(xp[0]); i++) { + xp[i] = range*i/(100-1); + yp[i] = genrand_real2(); + } + + acc = gsl_interp_accel_alloc(); + spline = gsl_spline_alloc(gsl_interp_linear,100); + gsl_spline_init(spline,xp,yp,100); + + for (i = 0; i < 1000; i++) { + x = genrand_real2()*range; + expected = gsl_spline_eval(spline,x,acc); + y = interp1d(x,xp,yp,100); + if (!isclose(y, expected, 1e-9, 1e-9)) { + sprintf(err, "interp1d(%.2g) returned %.5g, but expected %.5g", x, y, expected); + goto err; + } + } + + gsl_interp_accel_free(acc); + gsl_spline_free(spline); + + return 0; + +err: + gsl_interp_accel_free(acc); + gsl_spline_free(spline); + + return 1; +} + struct tests { testFunction *test; char *name; @@ -408,7 +452,8 @@ struct tests { {test_muon_get_T, "test_muon_get_T"}, {test_logsumexp, "test_logsumexp"}, {test_sno_charge_integral, "test_sno_charge_integral"}, - {test_path, "test_path"} + {test_path, "test_path"}, + {test_interp1d, "test_interp1d"} }; int main(int argc, char **argv) -- cgit