Merge pull request #122 from peddie/master

Add a new module containing the GSL interpolation interface

3 files changed, 396 insertions, 0 deletions

diff --git a/packages/gsl/src/Numeric/GSL.hs b/packages/gsl/src/Numeric/GSL.hs
index 61b8d7b..3f546bf 100644
--- a/packages/gsl/src/Numeric/GSL.hs
+++ b/packages/gsl/src/Numeric/GSL.hs
@@ -22,6 +22,7 @@ module Numeric.GSL (
 , module Numeric.GSL.Root
 , module Numeric.GSL.ODE
 , module Numeric.GSL.Fitting
+, module Numeric.GSL.Interpolation
 , module Data.Complex
 , setErrorHandlerOff
 ) where
@@ -34,6 +35,7 @@ import Numeric.GSL.Minimization
 import Numeric.GSL.Root
 import Numeric.GSL.ODE
 import Numeric.GSL.Fitting
+import Numeric.GSL.Interpolation
 import Data.Complex
 
 
diff --git a/packages/gsl/src/Numeric/GSL/Interpolation.hs b/packages/gsl/src/Numeric/GSL/Interpolation.hs
new file mode 100644
index 0000000..4d72ee2
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Interpolation.hs
@@ -0,0 +1,282 @@
+{- |
+Module      :  Numeric.GSL.Interpolation
+Copyright   :  (c) Matthew Peddie 2015
+License     :  GPL
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+Interpolation routines.
+
+<https://www.gnu.org/software/gsl/manual/html_node/Interpolation.html#Interpolation>
+
+The GSL routines @gsl_spline_eval@ and friends are used, but in spite
+of the names, they are not restricted to spline interpolation.  The
+functions in this module will work for any 'InterpolationMethod'.
+
+-}
+
+
+module Numeric.GSL.Interpolation (
+  -- * Interpolation methods
+  InterpolationMethod(..)
+  -- * Evaluation of interpolated functions
+  , evaluate
+  , evaluateV
+    -- * Evaluation of derivatives of interpolated functions
+  , evaluateDerivative
+  , evaluateDerivative2
+  , evaluateDerivativeV
+  , evaluateDerivative2V
+    -- * Evaluation of integrals of interpolated functions
+  , evaluateIntegral
+  , evaluateIntegralV
+) where
+
+import Data.Packed.Vector(Vector, fromList, dim)
+import Data.Packed.Foreign(appVector)
+import Foreign.C.Types
+import Foreign.Marshal.Alloc(alloca)
+import Foreign.Ptr(Ptr)
+import Foreign.Storable(peek)
+import Numeric.GSL.Internal
+import System.IO.Unsafe(unsafePerformIO)
+
+data InterpolationMethod = Linear
+                         | Polynomial
+                         | CSpline
+                         | CSplinePeriodic
+                         | Akima
+                         | AkimaPeriodic
+                         deriving (Eq, Show, Read)
+
+methodToInt :: Integral a => InterpolationMethod -> a
+methodToInt Linear = 0
+methodToInt Polynomial = 1
+methodToInt CSpline = 2
+methodToInt CSplinePeriodic = 3
+methodToInt Akima = 4
+methodToInt AkimaPeriodic = 5
+
+applyCFun hsname cname fun mth xs ys x
+  | dim xs /= dim ys = error $
+                         "Error: Vectors of unequal sizes " ++
+                         show (dim xs) ++ " and " ++ show (dim ys) ++
+                         " supplied to " ++ hsname
+  | otherwise = unsafePerformIO $
+      flip appVector xs $ \xs' ->
+       flip appVector ys $ \ys' ->
+        alloca $ \y' -> do
+          fun xs' ys'
+            (fromIntegral $ dim xs) x
+            (methodToInt mth) y' // check cname
+          peek y'
+
+foreign import ccall safe "spline_eval" c_spline_eval
+  :: Ptr Double -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
+
+--------------------------------------------------------------------
+{- | Evaluate a function by interpolating within the given dataset.  For
+example:
+
+>>> let xs = vector [1..10]
+>>> let ys = vector $ map (**2) [1..10]
+>>> evaluateV CSpline xs ys 2.2
+4.818867924528303
+
+To successfully @evaluateV xs ys x@, the vectors of corresponding
+domain-range values @xs@ and @ys@ must have identical lengths, and
+@xs@ must be monotonically increasing.  The evaluation point @x@ must
+lie between the smallest and largest values in @xs@.
+
+-}
+evaluateV :: InterpolationMethod  -- ^ What method to use to interpolate
+             -> Vector Double     -- ^ Data points sampling the domain of the function
+             -> Vector Double     -- ^ Data points sampling the range of the function
+             -> Double            -- ^ Point at which to evaluate the function
+             -> Double            -- ^ Interpolated result
+evaluateV = applyCFun "evaluateV" "spline_eval" c_spline_eval
+
+{- | Evaluate a function by interpolating within the given dataset.  For
+example:
+
+>>> let xs = [1..10]
+>>> let ys map (**2) [1..10]
+>>> evaluate Akima (zip xs ys) 2.2
+4.840000000000001
+
+To successfully @evaluate points x@, the domain (@x@) values in
+@points@ must be monotonically increasing.  The evaluation point @x@
+must lie between the smallest and largest values in the sampled
+domain.
+
+-}
+evaluate :: InterpolationMethod    -- ^ What method to use to interpolate
+            -> [(Double, Double)]  -- ^ (domain, range) values sampling the function
+            -> Double              -- ^ Point at which to evaluate the function
+            -> Double              -- ^ Interpolated result
+evaluate mth pts =
+  applyCFun "evaluate" "spline_eval" c_spline_eval_deriv
+  mth (fromList xs) (fromList ys)
+  where
+    (xs, ys) = unzip pts
+
+foreign import ccall safe "spline_eval_deriv" c_spline_eval_deriv
+  :: Ptr Double -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
+
+{- | Evaluate the derivative of a function by interpolating within the
+given dataset.  For example:
+
+>>> let xs = vector [1..10]
+>>> let ys = vector $ map (**2) [1..10]
+>>> evaluateDerivativeV CSpline xs ys 2.2
+4.338867924528302
+
+To successfully @evaluateDerivativeV xs ys x@, the vectors of
+corresponding domain-range values @xs@ and @ys@ must have identical
+lengths, and @xs@ must be monotonically increasing.  The interpolation
+point @x@ must lie between the smallest and largest values in @xs@.
+
+-}
+evaluateDerivativeV :: InterpolationMethod  -- ^ What method to use to interpolate
+                       -> Vector Double     -- ^ Data points @xs@ sampling the domain of the function
+                       -> Vector Double     -- ^ Data points @ys@ sampling the range of the function
+                       -> Double            -- ^ Point @x@ at which to evaluate the derivative
+                       -> Double            -- ^ Interpolated result
+evaluateDerivativeV =
+  applyCFun "evaluateDerivativeV" "spline_eval_deriv" c_spline_eval_deriv
+
+{- | Evaluate the derivative of a function by interpolating within the
+given dataset.  For example:
+
+>>> let xs = [1..10]
+>>> let ys map (**2) [1..10]
+>>> evaluateDerivative Akima (zip xs ys) 2.2
+4.4
+
+To successfully @evaluateDerivative points x@, the domain (@x@) values
+in @points@ must be monotonically increasing.  The evaluation point
+@x@ must lie between the smallest and largest values in the sampled
+domain.
+
+-}
+evaluateDerivative :: InterpolationMethod    -- ^ What method to use to interpolate
+                      -> [(Double, Double)]  -- ^ (domain, range) points sampling the function
+                      -> Double              -- ^ Point @x@ at which to evaluate the derivative
+                      -> Double              -- ^ Interpolated result
+evaluateDerivative mth pts =
+  applyCFun "evaluateDerivative" "spline_eval_deriv" c_spline_eval_deriv
+  mth (fromList xs) (fromList ys)
+  where
+    (xs, ys) = unzip pts
+
+foreign import ccall safe "spline_eval_deriv2" c_spline_eval_deriv2
+  :: Ptr Double -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
+
+{- | Evaluate the second derivative of a function by interpolating within the
+given dataset.  For example:
+
+>>> let xs = vector [1..10]
+>>> let ys = vector $ map (**2) [1..10]
+>>> evaluateDerivative2V CSpline xs ys 2.2
+2.4
+
+To successfully @evaluateDerivative2V xs ys x@, the vectors @xs@ and
+@ys@ must have identical lengths, and @xs@ must be monotonically
+increasing.  The evaluation point @x@ must lie between the smallest
+and largest values in @xs@.
+
+-}
+evaluateDerivative2V :: InterpolationMethod  -- ^ What method to use to interpolate
+                        -> Vector Double     -- ^ Data points @xs@ sampling the domain of the function
+                        -> Vector Double     -- ^ Data points @ys@ sampling the range of the function
+                        -> Double            -- ^ Point @x@ at which to evaluate the second derivative
+                        -> Double            -- ^ Interpolated result
+evaluateDerivative2V =
+  applyCFun "evaluateDerivative2V" "spline_eval_deriv2" c_spline_eval_deriv2
+
+{- | Evaluate the second derivative of a function by interpolating
+within the given dataset.  For example:
+
+>>> let xs = [1..10]
+>>> let ys map (**2) [1..10]
+>>> evaluateDerivative2 Akima (zip xs ys) 2.2
+2.0
+
+To successfully @evaluateDerivative2 points x@, the domain (@x@)
+values in @points@ must be monotonically increasing.  The evaluation
+point @x@ must lie between the smallest and largest values in the
+sampled domain.
+
+-}
+evaluateDerivative2 :: InterpolationMethod    -- ^ What method to use to interpolate
+                       -> [(Double, Double)]  -- ^ (domain, range) points sampling the function
+                       -> Double              -- ^ Point @x@ at which to evaluate the second derivative
+                       -> Double              -- ^ Interpolated result
+evaluateDerivative2 mth pts =
+  applyCFun "evaluateDerivative2" "spline_eval_deriv2" c_spline_eval_deriv2
+  mth (fromList xs) (fromList ys)
+  where
+    (xs, ys) = unzip pts
+
+foreign import ccall safe "spline_eval_integ" c_spline_eval_integ
+  :: Ptr Double -> Ptr Double -> CUInt -> Double -> Double -> CInt -> Ptr Double -> IO CInt
+
+applyCIntFun hsname cname fun mth xs ys a b
+  | dim xs /= dim ys = error $
+                         "Error: Vectors of unequal sizes " ++
+                         show (dim xs) ++ " and " ++ show (dim ys) ++
+                         " supplied to " ++ hsname
+  | otherwise = unsafePerformIO $
+      flip appVector xs $ \xs' ->
+       flip appVector ys $ \ys' ->
+        alloca $ \y' -> do
+          fun xs' ys'
+            (fromIntegral $ dim xs) a b
+            (methodToInt mth) y' // check cname
+          peek y'
+
+{- | Evaluate the definite integral of a function by interpolating
+within the given dataset.  For example:
+
+>>> let xs = vector [1..10]
+>>> let ys = vector $ map (**2) [1..10]
+>>> evaluateIntegralV CSpline xs ys 2.2 5.5
+51.89853207547169
+
+To successfully @evaluateIntegralV xs ys a b@, the vectors @xs@ and
+@ys@ must have identical lengths, and @xs@ must be monotonically
+increasing.  The integration bounds @a@ and @b@ must lie between the
+smallest and largest values in @xs@.
+
+-}
+evaluateIntegralV :: InterpolationMethod  -- ^ What method to use to interpolate
+                     -> Vector Double     -- ^ Data points @xs@ sampling the domain of the function
+                     -> Vector Double     -- ^ Data points @ys@ sampling the range of the function
+                     -> Double            -- ^ Lower integration bound @a@
+                     -> Double            -- ^ Upper integration bound @b@
+                     -> Double            -- ^ Resulting area
+evaluateIntegralV =
+  applyCIntFun "evaluateIntegralV" "spline_eval_integ" c_spline_eval_integ
+
+{- | Evaluate the definite integral of a function by interpolating
+within the given dataset.  For example:
+
+>>> let xs = [1..10]
+>>> let ys = map (**2) [1..10]
+>>> evaluateIntegralV CSpline (zip xs ys) (2.2, 5.5)
+51.909
+
+To successfully @evaluateIntegral points (a, b)@, the domain (@x@)
+values of @points@ must be monotonically increasing.  The integration
+bounds @a@ and @b@ must lie between the smallest and largest values in
+the sampled domain..
+-}
+evaluateIntegral :: InterpolationMethod    -- ^ What method to use to interpolate
+                    -> [(Double, Double)]  -- ^ (domain, range) points sampling the function
+                    -> (Double, Double)    -- ^ Integration bounds (@a@, @b@)
+                    -> Double              -- ^ Resulting area
+evaluateIntegral mth pts (a, b) =
+  applyCIntFun "evaluateIntegral" "spline_eval_integ" c_spline_eval_integ
+  mth (fromList xs) (fromList ys) a b
+  where
+    (xs, ys) = unzip pts
diff --git a/packages/gsl/src/Numeric/GSL/gsl-aux.c b/packages/gsl/src/Numeric/GSL/gsl-aux.c
index ffc5c20..e1b189c 100644
--- a/packages/gsl/src/Numeric/GSL/gsl-aux.c
+++ b/packages/gsl/src/Numeric/GSL/gsl-aux.c
@@ -34,6 +34,7 @@
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
 #include <gsl/gsl_roots.h>
+#include <gsl/gsl_spline.h>
 #include <gsl/gsl_multifit_nlin.h>
 #include <string.h>
 #include <stdio.h>
@@ -140,6 +141,117 @@ int deriv(int code, double f(double, void*), double x, double h, double * result
     return 0;
 }
 
+int spline_eval(const double xa[], const double ya[], unsigned int size,
+                double x, int method, double *y) {
+  DEBUGMSG("spline_eval");
+  const gsl_interp_type *T;
+  switch (method) {
+    case 0: { T = gsl_interp_linear; break; }
+    case 1: { T = gsl_interp_polynomial; break; }
+    case 2: { T = gsl_interp_cspline; break; }
+    case 3: { T = gsl_interp_cspline_periodic; break; }
+    case 4: { T = gsl_interp_akima; break; }
+    case 5: { T = gsl_interp_akima_periodic; break; }
+    default: ERROR(BAD_CODE);
+  }
+
+  gsl_spline *spline = gsl_spline_alloc(T, size);
+  if (NULL == spline) ERROR(MEM);
+  const int initres = gsl_spline_init(spline, xa, ya, size);
+  CHECK(initres,initres);
+  gsl_interp_accel *acc = gsl_interp_accel_alloc();
+  if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); };
+
+  const int res = gsl_spline_eval_e(spline, x, acc, y);
+  CHECK(res,res);
+  gsl_interp_accel_free(acc);
+  gsl_spline_free(spline);
+  OK
+}
+
+int spline_eval_deriv(const double xa[], const double ya[], unsigned int size,
+                      double x, int method, double *y) {
+  DEBUGMSG("spline_eval_deriv");
+  const gsl_interp_type *T;
+  switch (method) {
+    case 0: { T = gsl_interp_linear; break; }
+    case 1: { T = gsl_interp_polynomial; break; }
+    case 2: { T = gsl_interp_cspline; break; }
+    case 3: { T = gsl_interp_cspline_periodic; break; }
+    case 4: { T = gsl_interp_akima; break; }
+    case 5: { T = gsl_interp_akima_periodic; break; }
+    default: ERROR(BAD_CODE);
+  }
+
+  gsl_spline *spline = gsl_spline_alloc(T, size);
+  if (NULL == spline) ERROR(MEM);
+  const int initres = gsl_spline_init(spline, xa, ya, size);
+  CHECK(initres,initres);
+  gsl_interp_accel *acc = gsl_interp_accel_alloc();
+  if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); };
+
+  const int res = gsl_spline_eval_deriv_e(spline, x, acc, y);
+  CHECK(res,res);
+  gsl_interp_accel_free(acc);
+  gsl_spline_free(spline);
+  OK
+}
+
+int spline_eval_deriv2(const double xa[], const double ya[], unsigned int size,
+                       double x, int method, double *y) {
+  DEBUGMSG("spline_eval_deriv2");
+  const gsl_interp_type *T;
+  switch (method) {
+    case 0: { T = gsl_interp_linear; break; }
+    case 1: { T = gsl_interp_polynomial; break; }
+    case 2: { T = gsl_interp_cspline; break; }
+    case 3: { T = gsl_interp_cspline_periodic; break; }
+    case 4: { T = gsl_interp_akima; break; }
+    case 5: { T = gsl_interp_akima_periodic; break; }
+    default: ERROR(BAD_CODE);
+  }
+
+  gsl_spline *spline = gsl_spline_alloc(T, size);
+  if (NULL == spline) ERROR(MEM);
+  const int initres = gsl_spline_init(spline, xa, ya, size);
+  CHECK(initres,initres);
+  gsl_interp_accel *acc = gsl_interp_accel_alloc();
+  if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); };
+
+  const int res = gsl_spline_eval_deriv2_e(spline, x, acc, y);
+  CHECK(res,res);
+  gsl_interp_accel_free(acc);
+  gsl_spline_free(spline);
+  OK
+}
+
+int spline_eval_integ(const double xa[], const double ya[], unsigned int size,
+                      double a, double b, int method, double *y) {
+  DEBUGMSG("spline_eval_integ");
+  const gsl_interp_type *T;
+  switch (method) {
+    case 0: { T = gsl_interp_linear; break; }
+    case 1: { T = gsl_interp_polynomial; break; }
+    case 2: { T = gsl_interp_cspline; break; }
+    case 3: { T = gsl_interp_cspline_periodic; break; }
+    case 4: { T = gsl_interp_akima; break; }
+    case 5: { T = gsl_interp_akima_periodic; break; }
+    default: ERROR(BAD_CODE);
+  }
+
+  gsl_spline *spline = gsl_spline_alloc(T, size);
+  if (NULL == spline) ERROR(MEM);
+  const int initres = gsl_spline_init(spline, xa, ya, size);
+  CHECK(initres,initres);
+  gsl_interp_accel *acc = gsl_interp_accel_alloc();
+  if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); };
+
+  const int res = gsl_spline_eval_integ_e(spline, a, b, acc, y);
+  CHECK(res,res);
+  gsl_interp_accel_free(acc);
+  gsl_spline_free(spline);
+  OK
+}
 
 int integrate_qng(double f(double, void*), double a, double b, double aprec, double prec,
                    double *result, double*error) {