hmatrix-gsl-0.18.0.1: src/Numeric/GSL/Interpolation.hs
{-# LANGUAGE MagicHash, UnboxedTuples #-}
{- |
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 Numeric.LinearAlgebra(Vector, fromList, size, Numeric)
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)
-- FIXME
import qualified Data.Vector.Storable as S
import GHC.Base (IO(..), realWorld#)
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
dim :: Numeric t => Vector t -> Int
dim = size
-- FIXME
appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f))
unsafeInlinePerformIO (IO f) = case f realWorld# of
(# _, x #) -> x
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
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