{-# LANGUAGE MagicHash, UnboxedTuples #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
module Numeric.GSL.Interpolation (
InterpolationMethod(..)
, evaluate
, evaluateV
, evaluateDerivative
, evaluateDerivative2
, evaluateDerivativeV
, evaluateDerivative2V
, 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)
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
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
evaluateV :: InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
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 :: InterpolationMethod
-> [(Double, Double)]
-> Double
-> Double
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
evaluateDerivativeV :: InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
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
evaluateDerivative :: InterpolationMethod
-> [(Double, Double)]
-> Double
-> Double
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
evaluateDerivative2V :: InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
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
evaluateDerivative2 :: InterpolationMethod
-> [(Double, Double)]
-> Double
-> Double
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'
evaluateIntegralV :: InterpolationMethod
-> Vector Double
-> Vector Double
-> Double
-> Double
-> Double
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
evaluateIntegral :: InterpolationMethod
-> [(Double, Double)]
-> (Double, Double)
-> Double
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