path: root/packages/gsl/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