summaryrefslogtreecommitdiff
path: root/packages/gsl/src/Numeric/GSL/Interpolation.hs
blob: d060468f0e90d9d4749c6865a9a30071bc2738d6 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
{- |
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)

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

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