From a25d6f45d7131bd86508f82870d2e24606a811e8 Mon Sep 17 00:00:00 2001 From: Matthew Peddie Date: Sun, 26 Apr 2015 12:17:00 +1000 Subject: Add a new module containing the GSL interpolation interface --- packages/gsl/hmatrix-gsl.cabal | 1 + packages/gsl/src/Numeric/GSL.hs | 2 + packages/gsl/src/Numeric/GSL/Interpolation.hs | 282 ++++++++++++++++++++++++++ packages/gsl/src/Numeric/GSL/gsl-aux.c | 112 ++++++++++ 4 files changed, 397 insertions(+) create mode 100644 packages/gsl/src/Numeric/GSL/Interpolation.hs (limited to 'packages') diff --git a/packages/gsl/hmatrix-gsl.cabal b/packages/gsl/hmatrix-gsl.cabal index 95fa3ed..2f6f51b 100644 --- a/packages/gsl/hmatrix-gsl.cabal +++ b/packages/gsl/hmatrix-gsl.cabal @@ -43,6 +43,7 @@ library Numeric.GSL.ODE, Numeric.GSL, Numeric.GSL.LinearAlgebra, + Numeric.GSL.Interpolation, Graphics.Plot other-modules: Numeric.GSL.Internal, Numeric.GSL.Vector, 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. + + + +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 #include #include +#include #include #include #include @@ -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) { -- cgit v1.2.3