{-# LANGUAGE MagicHash, UnboxedTuples #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}

{- |
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 (InterpolationMethod -> InterpolationMethod -> Bool
(InterpolationMethod -> InterpolationMethod -> Bool)
-> (InterpolationMethod -> InterpolationMethod -> Bool)
-> Eq InterpolationMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: InterpolationMethod -> InterpolationMethod -> Bool
$c/= :: InterpolationMethod -> InterpolationMethod -> Bool
== :: InterpolationMethod -> InterpolationMethod -> Bool
$c== :: InterpolationMethod -> InterpolationMethod -> Bool
Eq, Int -> InterpolationMethod -> ShowS
[InterpolationMethod] -> ShowS
InterpolationMethod -> String
(Int -> InterpolationMethod -> ShowS)
-> (InterpolationMethod -> String)
-> ([InterpolationMethod] -> ShowS)
-> Show InterpolationMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [InterpolationMethod] -> ShowS
$cshowList :: [InterpolationMethod] -> ShowS
show :: InterpolationMethod -> String
$cshow :: InterpolationMethod -> String
showsPrec :: Int -> InterpolationMethod -> ShowS
$cshowsPrec :: Int -> InterpolationMethod -> ShowS
Show, ReadPrec [InterpolationMethod]
ReadPrec InterpolationMethod
Int -> ReadS InterpolationMethod
ReadS [InterpolationMethod]
(Int -> ReadS InterpolationMethod)
-> ReadS [InterpolationMethod]
-> ReadPrec InterpolationMethod
-> ReadPrec [InterpolationMethod]
-> Read InterpolationMethod
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [InterpolationMethod]
$creadListPrec :: ReadPrec [InterpolationMethod]
readPrec :: ReadPrec InterpolationMethod
$creadPrec :: ReadPrec InterpolationMethod
readList :: ReadS [InterpolationMethod]
$creadList :: ReadS [InterpolationMethod]
readsPrec :: Int -> ReadS InterpolationMethod
$creadsPrec :: Int -> ReadS InterpolationMethod
Read)

methodToInt :: Integral a => InterpolationMethod -> a
methodToInt :: 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 :: Vector t -> Int
dim = Vector t -> Int
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size

-- FIXME
appVector :: (Ptr a -> a) -> Vector a -> a
appVector f :: Ptr a -> a
f x :: Vector a
x = IO a -> a
forall a. IO a -> a
unsafeInlinePerformIO (Vector a -> (Ptr a -> IO a) -> IO a
forall a b. Storable a => Vector a -> (Ptr a -> IO b) -> IO b
S.unsafeWith Vector a
x (a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (a -> IO a) -> (Ptr a -> a) -> Ptr a -> IO a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Ptr a -> a
f))

unsafeInlinePerformIO :: IO a -> a
unsafeInlinePerformIO (IO f :: State# RealWorld -> (# State# RealWorld, a #)
f) = case State# RealWorld -> (# State# RealWorld, a #)
f State# RealWorld
realWorld# of
    (# _, x :: a
x #) -> a
x

applyCFun :: String
-> String
-> (Ptr t -> Ptr a -> t -> t -> t -> Ptr p -> IO CInt)
-> InterpolationMethod
-> Vector t
-> Vector a
-> t
-> p
applyCFun hsname :: String
hsname cname :: String
cname fun :: Ptr t -> Ptr a -> t -> t -> t -> Ptr p -> IO CInt
fun mth :: InterpolationMethod
mth xs :: Vector t
xs ys :: Vector a
ys x :: t
x
  | Vector t -> Int
forall t. Numeric t => Vector t -> Int
dim Vector t
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Vector a -> Int
forall t. Numeric t => Vector t -> Int
dim Vector a
ys = String -> p
forall a. HasCallStack => String -> a
error (String -> p) -> String -> p
forall a b. (a -> b) -> a -> b
$
                         "Error: Vectors of unequal sizes " String -> ShowS
forall a. [a] -> [a] -> [a]
++
                         Int -> String
forall a. Show a => a -> String
show (Vector t -> Int
forall t. Numeric t => Vector t -> Int
dim Vector t
xs) String -> ShowS
forall a. [a] -> [a] -> [a]
++ " and " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show (Vector a -> Int
forall t. Numeric t => Vector t -> Int
dim Vector a
ys) String -> ShowS
forall a. [a] -> [a] -> [a]
++
                         " supplied to " String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
hsname
  | Bool
otherwise = IO p -> p
forall a. IO a -> a
unsafePerformIO (IO p -> p) -> IO p -> p
forall a b. (a -> b) -> a -> b
$
      ((Ptr t -> IO p) -> Vector t -> IO p)
-> Vector t -> (Ptr t -> IO p) -> IO p
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Ptr t -> IO p) -> Vector t -> IO p
forall a a. Storable a => (Ptr a -> a) -> Vector a -> a
appVector Vector t
xs ((Ptr t -> IO p) -> IO p) -> (Ptr t -> IO p) -> IO p
forall a b. (a -> b) -> a -> b
$ \xs' :: Ptr t
xs' ->
       ((Ptr a -> IO p) -> Vector a -> IO p)
-> Vector a -> (Ptr a -> IO p) -> IO p
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Ptr a -> IO p) -> Vector a -> IO p
forall a a. Storable a => (Ptr a -> a) -> Vector a -> a
appVector Vector a
ys ((Ptr a -> IO p) -> IO p) -> (Ptr a -> IO p) -> IO p
forall a b. (a -> b) -> a -> b
$ \ys' :: Ptr a
ys' ->
        (Ptr p -> IO p) -> IO p
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr p -> IO p) -> IO p) -> (Ptr p -> IO p) -> IO p
forall a b. (a -> b) -> a -> b
$ \y' :: Ptr p
y' -> do
          Ptr t -> Ptr a -> t -> t -> t -> Ptr p -> IO CInt
fun Ptr t
xs' Ptr a
ys'
            (Int -> t
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> t) -> Int -> t
forall a b. (a -> b) -> a -> b
$ Vector t -> Int
forall t. Numeric t => Vector t -> Int
dim Vector t
xs) t
x
            (InterpolationMethod -> t
forall a. Integral a => InterpolationMethod -> a
methodToInt InterpolationMethod
mth) Ptr p
y' IO CInt -> (IO CInt -> IO ()) -> IO ()
forall x y. x -> (x -> y) -> y
// String -> IO CInt -> IO ()
check String
cname
          Ptr p -> IO p
forall a. Storable a => Ptr a -> IO a
peek Ptr p
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 :: InterpolationMethod
-> Vector Double -> Vector Double -> Double -> Double
evaluateV = String
-> String
-> (Ptr Double
    -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt)
-> InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
forall t a p t t t.
(Numeric t, Numeric a, Storable p, Integral t, Num t) =>
String
-> String
-> (Ptr t -> Ptr a -> t -> t -> t -> Ptr p -> IO CInt)
-> InterpolationMethod
-> Vector t
-> Vector a
-> t
-> p
applyCFun "evaluateV" "spline_eval" Ptr Double
-> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
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 :: InterpolationMethod -> [(Double, Double)] -> Double -> Double
evaluate mth :: InterpolationMethod
mth pts :: [(Double, Double)]
pts =
  String
-> String
-> (Ptr Double
    -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt)
-> InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
forall t a p t t t.
(Numeric t, Numeric a, Storable p, Integral t, Num t) =>
String
-> String
-> (Ptr t -> Ptr a -> t -> t -> t -> Ptr p -> IO CInt)
-> InterpolationMethod
-> Vector t
-> Vector a
-> t
-> p
applyCFun "evaluate" "spline_eval" Ptr Double
-> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
c_spline_eval
  InterpolationMethod
mth ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
xs) ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
ys)
  where
    (xs :: [Double]
xs, ys :: [Double]
ys) = [(Double, Double)] -> ([Double], [Double])
forall a b. [(a, b)] -> ([a], [b])
unzip [(Double, Double)]
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 :: InterpolationMethod
-> Vector Double -> Vector Double -> Double -> Double
evaluateDerivativeV =
  String
-> String
-> (Ptr Double
    -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt)
-> InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
forall t a p t t t.
(Numeric t, Numeric a, Storable p, Integral t, Num t) =>
String
-> String
-> (Ptr t -> Ptr a -> t -> t -> t -> Ptr p -> IO CInt)
-> InterpolationMethod
-> Vector t
-> Vector a
-> t
-> p
applyCFun "evaluateDerivativeV" "spline_eval_deriv" Ptr Double
-> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
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 :: InterpolationMethod -> [(Double, Double)] -> Double -> Double
evaluateDerivative mth :: InterpolationMethod
mth pts :: [(Double, Double)]
pts =
  String
-> String
-> (Ptr Double
    -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt)
-> InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
forall t a p t t t.
(Numeric t, Numeric a, Storable p, Integral t, Num t) =>
String
-> String
-> (Ptr t -> Ptr a -> t -> t -> t -> Ptr p -> IO CInt)
-> InterpolationMethod
-> Vector t
-> Vector a
-> t
-> p
applyCFun "evaluateDerivative" "spline_eval_deriv" Ptr Double
-> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
c_spline_eval_deriv
  InterpolationMethod
mth ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
xs) ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
ys)
  where
    (xs :: [Double]
xs, ys :: [Double]
ys) = [(Double, Double)] -> ([Double], [Double])
forall a b. [(a, b)] -> ([a], [b])
unzip [(Double, Double)]
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 :: InterpolationMethod
-> Vector Double -> Vector Double -> Double -> Double
evaluateDerivative2V =
  String
-> String
-> (Ptr Double
    -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt)
-> InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
forall t a p t t t.
(Numeric t, Numeric a, Storable p, Integral t, Num t) =>
String
-> String
-> (Ptr t -> Ptr a -> t -> t -> t -> Ptr p -> IO CInt)
-> InterpolationMethod
-> Vector t
-> Vector a
-> t
-> p
applyCFun "evaluateDerivative2V" "spline_eval_deriv2" Ptr Double
-> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
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 :: InterpolationMethod -> [(Double, Double)] -> Double -> Double
evaluateDerivative2 mth :: InterpolationMethod
mth pts :: [(Double, Double)]
pts =
  String
-> String
-> (Ptr Double
    -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt)
-> InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
forall t a p t t t.
(Numeric t, Numeric a, Storable p, Integral t, Num t) =>
String
-> String
-> (Ptr t -> Ptr a -> t -> t -> t -> Ptr p -> IO CInt)
-> InterpolationMethod
-> Vector t
-> Vector a
-> t
-> p
applyCFun "evaluateDerivative2" "spline_eval_deriv2" Ptr Double
-> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
c_spline_eval_deriv2
  InterpolationMethod
mth ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
xs) ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
ys)
  where
    (xs :: [Double]
xs, ys :: [Double]
ys) = [(Double, Double)] -> ([Double], [Double])
forall a b. [(a, b)] -> ([a], [b])
unzip [(Double, Double)]
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 :: String
-> String
-> (Ptr t -> Ptr a -> t -> t -> t -> t -> Ptr p -> IO CInt)
-> InterpolationMethod
-> Vector t
-> Vector a
-> t
-> t
-> p
applyCIntFun hsname :: String
hsname cname :: String
cname fun :: Ptr t -> Ptr a -> t -> t -> t -> t -> Ptr p -> IO CInt
fun mth :: InterpolationMethod
mth xs :: Vector t
xs ys :: Vector a
ys a :: t
a b :: t
b
  | Vector t -> Int
forall t. Numeric t => Vector t -> Int
dim Vector t
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Vector a -> Int
forall t. Numeric t => Vector t -> Int
dim Vector a
ys = String -> p
forall a. HasCallStack => String -> a
error (String -> p) -> String -> p
forall a b. (a -> b) -> a -> b
$
                         "Error: Vectors of unequal sizes " String -> ShowS
forall a. [a] -> [a] -> [a]
++
                         Int -> String
forall a. Show a => a -> String
show (Vector t -> Int
forall t. Numeric t => Vector t -> Int
dim Vector t
xs) String -> ShowS
forall a. [a] -> [a] -> [a]
++ " and " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show (Vector a -> Int
forall t. Numeric t => Vector t -> Int
dim Vector a
ys) String -> ShowS
forall a. [a] -> [a] -> [a]
++
                         " supplied to " String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
hsname
  | Bool
otherwise = IO p -> p
forall a. IO a -> a
unsafePerformIO (IO p -> p) -> IO p -> p
forall a b. (a -> b) -> a -> b
$
      ((Ptr t -> IO p) -> Vector t -> IO p)
-> Vector t -> (Ptr t -> IO p) -> IO p
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Ptr t -> IO p) -> Vector t -> IO p
forall a a. Storable a => (Ptr a -> a) -> Vector a -> a
appVector Vector t
xs ((Ptr t -> IO p) -> IO p) -> (Ptr t -> IO p) -> IO p
forall a b. (a -> b) -> a -> b
$ \xs' :: Ptr t
xs' ->
       ((Ptr a -> IO p) -> Vector a -> IO p)
-> Vector a -> (Ptr a -> IO p) -> IO p
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Ptr a -> IO p) -> Vector a -> IO p
forall a a. Storable a => (Ptr a -> a) -> Vector a -> a
appVector Vector a
ys ((Ptr a -> IO p) -> IO p) -> (Ptr a -> IO p) -> IO p
forall a b. (a -> b) -> a -> b
$ \ys' :: Ptr a
ys' ->
        (Ptr p -> IO p) -> IO p
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr p -> IO p) -> IO p) -> (Ptr p -> IO p) -> IO p
forall a b. (a -> b) -> a -> b
$ \y' :: Ptr p
y' -> do
          Ptr t -> Ptr a -> t -> t -> t -> t -> Ptr p -> IO CInt
fun Ptr t
xs' Ptr a
ys'
            (Int -> t
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> t) -> Int -> t
forall a b. (a -> b) -> a -> b
$ Vector t -> Int
forall t. Numeric t => Vector t -> Int
dim Vector t
xs) t
a t
b
            (InterpolationMethod -> t
forall a. Integral a => InterpolationMethod -> a
methodToInt InterpolationMethod
mth) Ptr p
y' IO CInt -> (IO CInt -> IO ()) -> IO ()
forall x y. x -> (x -> y) -> y
// String -> IO CInt -> IO ()
check String
cname
          Ptr p -> IO p
forall a. Storable a => Ptr a -> IO a
peek Ptr p
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 :: InterpolationMethod
-> Vector Double -> Vector Double -> Double -> Double -> Double
evaluateIntegralV =
  String
-> String
-> (Ptr Double
    -> Ptr Double
    -> CUInt
    -> Double
    -> Double
    -> CInt
    -> Ptr Double
    -> IO CInt)
-> InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
-> Double
forall t a p t t t t.
(Numeric t, Numeric a, Storable p, Integral t, Num t) =>
String
-> String
-> (Ptr t -> Ptr a -> t -> t -> t -> t -> Ptr p -> IO CInt)
-> InterpolationMethod
-> Vector t
-> Vector a
-> t
-> t
-> p
applyCIntFun "evaluateIntegralV" "spline_eval_integ" Ptr Double
-> Ptr Double
-> CUInt
-> Double
-> Double
-> CInt
-> Ptr Double
-> IO CInt
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 :: InterpolationMethod
-> [(Double, Double)] -> (Double, Double) -> Double
evaluateIntegral mth :: InterpolationMethod
mth pts :: [(Double, Double)]
pts (a :: Double
a, b :: Double
b) =
  String
-> String
-> (Ptr Double
    -> Ptr Double
    -> CUInt
    -> Double
    -> Double
    -> CInt
    -> Ptr Double
    -> IO CInt)
-> InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
-> Double
forall t a p t t t t.
(Numeric t, Numeric a, Storable p, Integral t, Num t) =>
String
-> String
-> (Ptr t -> Ptr a -> t -> t -> t -> t -> Ptr p -> IO CInt)
-> InterpolationMethod
-> Vector t
-> Vector a
-> t
-> t
-> p
applyCIntFun "evaluateIntegral" "spline_eval_integ" Ptr Double
-> Ptr Double
-> CUInt
-> Double
-> Double
-> CInt
-> Ptr Double
-> IO CInt
c_spline_eval_integ
  InterpolationMethod
mth ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
xs) ([Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
ys) Double
a Double
b
  where
    (xs :: [Double]
xs, ys :: [Double]
ys) = [(Double, Double)] -> ([Double], [Double])
forall a b. [(a, b)] -> ([a], [b])
unzip [(Double, Double)]
pts