From 83119fe9bb0c9c2a70ea5397aca9968bde7ffa07 Mon Sep 17 00:00:00 2001
From: tlatorre <tlatorre@uchicago.edu>
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 <math.h>
 #include <stdlib.h> /* 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 <stdlib.h> /* 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 <gsl/gsl_spline.h>
 
 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