summaryrefslogtreecommitdiff
path: root/packages/gsl/src/Numeric/GSL
diff options
context:
space:
mode:
authorMatthew Peddie <mpeddie@gmail.com>2015-04-26 12:17:00 +1000
committerMatthew Peddie <mpeddie@gmail.com>2015-04-26 12:17:00 +1000
commita25d6f45d7131bd86508f82870d2e24606a811e8 (patch)
tree52df8c1e8dda157896d4884bbd0dbfb439132d21 /packages/gsl/src/Numeric/GSL
parente282d2b4c5f389941225d02df94e7baf38c4f38e (diff)
Add a new module containing the GSL interpolation interface
Diffstat (limited to 'packages/gsl/src/Numeric/GSL')
-rw-r--r--packages/gsl/src/Numeric/GSL/Interpolation.hs282
-rw-r--r--packages/gsl/src/Numeric/GSL/gsl-aux.c112
2 files changed, 394 insertions, 0 deletions
diff --git a/packages/gsl/src/Numeric/GSL/Interpolation.hs b/packages/gsl/src/Numeric/GSL/Interpolation.hs
new file mode 100644
index 0000000..4d72ee2
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Interpolation.hs
@@ -0,0 +1,282 @@
1{- |
2Module : Numeric.GSL.Interpolation
3Copyright : (c) Matthew Peddie 2015
4License : GPL
5Maintainer : Alberto Ruiz
6Stability : provisional
7
8Interpolation routines.
9
10<https://www.gnu.org/software/gsl/manual/html_node/Interpolation.html#Interpolation>
11
12The GSL routines @gsl_spline_eval@ and friends are used, but in spite
13of the names, they are not restricted to spline interpolation. The
14functions in this module will work for any 'InterpolationMethod'.
15
16-}
17
18
19module Numeric.GSL.Interpolation (
20 -- * Interpolation methods
21 InterpolationMethod(..)
22 -- * Evaluation of interpolated functions
23 , evaluate
24 , evaluateV
25 -- * Evaluation of derivatives of interpolated functions
26 , evaluateDerivative
27 , evaluateDerivative2
28 , evaluateDerivativeV
29 , evaluateDerivative2V
30 -- * Evaluation of integrals of interpolated functions
31 , evaluateIntegral
32 , evaluateIntegralV
33) where
34
35import Data.Packed.Vector(Vector, fromList, dim)
36import Data.Packed.Foreign(appVector)
37import Foreign.C.Types
38import Foreign.Marshal.Alloc(alloca)
39import Foreign.Ptr(Ptr)
40import Foreign.Storable(peek)
41import Numeric.GSL.Internal
42import System.IO.Unsafe(unsafePerformIO)
43
44data InterpolationMethod = Linear
45 | Polynomial
46 | CSpline
47 | CSplinePeriodic
48 | Akima
49 | AkimaPeriodic
50 deriving (Eq, Show, Read)
51
52methodToInt :: Integral a => InterpolationMethod -> a
53methodToInt Linear = 0
54methodToInt Polynomial = 1
55methodToInt CSpline = 2
56methodToInt CSplinePeriodic = 3
57methodToInt Akima = 4
58methodToInt AkimaPeriodic = 5
59
60applyCFun hsname cname fun mth xs ys x
61 | dim xs /= dim ys = error $
62 "Error: Vectors of unequal sizes " ++
63 show (dim xs) ++ " and " ++ show (dim ys) ++
64 " supplied to " ++ hsname
65 | otherwise = unsafePerformIO $
66 flip appVector xs $ \xs' ->
67 flip appVector ys $ \ys' ->
68 alloca $ \y' -> do
69 fun xs' ys'
70 (fromIntegral $ dim xs) x
71 (methodToInt mth) y' // check cname
72 peek y'
73
74foreign import ccall safe "spline_eval" c_spline_eval
75 :: Ptr Double -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
76
77--------------------------------------------------------------------
78{- | Evaluate a function by interpolating within the given dataset. For
79example:
80
81>>> let xs = vector [1..10]
82>>> let ys = vector $ map (**2) [1..10]
83>>> evaluateV CSpline xs ys 2.2
844.818867924528303
85
86To successfully @evaluateV xs ys x@, the vectors of corresponding
87domain-range values @xs@ and @ys@ must have identical lengths, and
88@xs@ must be monotonically increasing. The evaluation point @x@ must
89lie between the smallest and largest values in @xs@.
90
91-}
92evaluateV :: InterpolationMethod -- ^ What method to use to interpolate
93 -> Vector Double -- ^ Data points sampling the domain of the function
94 -> Vector Double -- ^ Data points sampling the range of the function
95 -> Double -- ^ Point at which to evaluate the function
96 -> Double -- ^ Interpolated result
97evaluateV = applyCFun "evaluateV" "spline_eval" c_spline_eval
98
99{- | Evaluate a function by interpolating within the given dataset. For
100example:
101
102>>> let xs = [1..10]
103>>> let ys map (**2) [1..10]
104>>> evaluate Akima (zip xs ys) 2.2
1054.840000000000001
106
107To successfully @evaluate points x@, the domain (@x@) values in
108@points@ must be monotonically increasing. The evaluation point @x@
109must lie between the smallest and largest values in the sampled
110domain.
111
112-}
113evaluate :: InterpolationMethod -- ^ What method to use to interpolate
114 -> [(Double, Double)] -- ^ (domain, range) values sampling the function
115 -> Double -- ^ Point at which to evaluate the function
116 -> Double -- ^ Interpolated result
117evaluate mth pts =
118 applyCFun "evaluate" "spline_eval" c_spline_eval_deriv
119 mth (fromList xs) (fromList ys)
120 where
121 (xs, ys) = unzip pts
122
123foreign import ccall safe "spline_eval_deriv" c_spline_eval_deriv
124 :: Ptr Double -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
125
126{- | Evaluate the derivative of a function by interpolating within the
127given dataset. For example:
128
129>>> let xs = vector [1..10]
130>>> let ys = vector $ map (**2) [1..10]
131>>> evaluateDerivativeV CSpline xs ys 2.2
1324.338867924528302
133
134To successfully @evaluateDerivativeV xs ys x@, the vectors of
135corresponding domain-range values @xs@ and @ys@ must have identical
136lengths, and @xs@ must be monotonically increasing. The interpolation
137point @x@ must lie between the smallest and largest values in @xs@.
138
139-}
140evaluateDerivativeV :: InterpolationMethod -- ^ What method to use to interpolate
141 -> Vector Double -- ^ Data points @xs@ sampling the domain of the function
142 -> Vector Double -- ^ Data points @ys@ sampling the range of the function
143 -> Double -- ^ Point @x@ at which to evaluate the derivative
144 -> Double -- ^ Interpolated result
145evaluateDerivativeV =
146 applyCFun "evaluateDerivativeV" "spline_eval_deriv" c_spline_eval_deriv
147
148{- | Evaluate the derivative of a function by interpolating within the
149given dataset. For example:
150
151>>> let xs = [1..10]
152>>> let ys map (**2) [1..10]
153>>> evaluateDerivative Akima (zip xs ys) 2.2
1544.4
155
156To successfully @evaluateDerivative points x@, the domain (@x@) values
157in @points@ must be monotonically increasing. The evaluation point
158@x@ must lie between the smallest and largest values in the sampled
159domain.
160
161-}
162evaluateDerivative :: InterpolationMethod -- ^ What method to use to interpolate
163 -> [(Double, Double)] -- ^ (domain, range) points sampling the function
164 -> Double -- ^ Point @x@ at which to evaluate the derivative
165 -> Double -- ^ Interpolated result
166evaluateDerivative mth pts =
167 applyCFun "evaluateDerivative" "spline_eval_deriv" c_spline_eval_deriv
168 mth (fromList xs) (fromList ys)
169 where
170 (xs, ys) = unzip pts
171
172foreign import ccall safe "spline_eval_deriv2" c_spline_eval_deriv2
173 :: Ptr Double -> Ptr Double -> CUInt -> Double -> CInt -> Ptr Double -> IO CInt
174
175{- | Evaluate the second derivative of a function by interpolating within the
176given dataset. For example:
177
178>>> let xs = vector [1..10]
179>>> let ys = vector $ map (**2) [1..10]
180>>> evaluateDerivative2V CSpline xs ys 2.2
1812.4
182
183To successfully @evaluateDerivative2V xs ys x@, the vectors @xs@ and
184@ys@ must have identical lengths, and @xs@ must be monotonically
185increasing. The evaluation point @x@ must lie between the smallest
186and largest values in @xs@.
187
188-}
189evaluateDerivative2V :: InterpolationMethod -- ^ What method to use to interpolate
190 -> Vector Double -- ^ Data points @xs@ sampling the domain of the function
191 -> Vector Double -- ^ Data points @ys@ sampling the range of the function
192 -> Double -- ^ Point @x@ at which to evaluate the second derivative
193 -> Double -- ^ Interpolated result
194evaluateDerivative2V =
195 applyCFun "evaluateDerivative2V" "spline_eval_deriv2" c_spline_eval_deriv2
196
197{- | Evaluate the second derivative of a function by interpolating
198within the given dataset. For example:
199
200>>> let xs = [1..10]
201>>> let ys map (**2) [1..10]
202>>> evaluateDerivative2 Akima (zip xs ys) 2.2
2032.0
204
205To successfully @evaluateDerivative2 points x@, the domain (@x@)
206values in @points@ must be monotonically increasing. The evaluation
207point @x@ must lie between the smallest and largest values in the
208sampled domain.
209
210-}
211evaluateDerivative2 :: InterpolationMethod -- ^ What method to use to interpolate
212 -> [(Double, Double)] -- ^ (domain, range) points sampling the function
213 -> Double -- ^ Point @x@ at which to evaluate the second derivative
214 -> Double -- ^ Interpolated result
215evaluateDerivative2 mth pts =
216 applyCFun "evaluateDerivative2" "spline_eval_deriv2" c_spline_eval_deriv2
217 mth (fromList xs) (fromList ys)
218 where
219 (xs, ys) = unzip pts
220
221foreign import ccall safe "spline_eval_integ" c_spline_eval_integ
222 :: Ptr Double -> Ptr Double -> CUInt -> Double -> Double -> CInt -> Ptr Double -> IO CInt
223
224applyCIntFun hsname cname fun mth xs ys a b
225 | dim xs /= dim ys = error $
226 "Error: Vectors of unequal sizes " ++
227 show (dim xs) ++ " and " ++ show (dim ys) ++
228 " supplied to " ++ hsname
229 | otherwise = unsafePerformIO $
230 flip appVector xs $ \xs' ->
231 flip appVector ys $ \ys' ->
232 alloca $ \y' -> do
233 fun xs' ys'
234 (fromIntegral $ dim xs) a b
235 (methodToInt mth) y' // check cname
236 peek y'
237
238{- | Evaluate the definite integral of a function by interpolating
239within the given dataset. For example:
240
241>>> let xs = vector [1..10]
242>>> let ys = vector $ map (**2) [1..10]
243>>> evaluateIntegralV CSpline xs ys 2.2 5.5
24451.89853207547169
245
246To successfully @evaluateIntegralV xs ys a b@, the vectors @xs@ and
247@ys@ must have identical lengths, and @xs@ must be monotonically
248increasing. The integration bounds @a@ and @b@ must lie between the
249smallest and largest values in @xs@.
250
251-}
252evaluateIntegralV :: InterpolationMethod -- ^ What method to use to interpolate
253 -> Vector Double -- ^ Data points @xs@ sampling the domain of the function
254 -> Vector Double -- ^ Data points @ys@ sampling the range of the function
255 -> Double -- ^ Lower integration bound @a@
256 -> Double -- ^ Upper integration bound @b@
257 -> Double -- ^ Resulting area
258evaluateIntegralV =
259 applyCIntFun "evaluateIntegralV" "spline_eval_integ" c_spline_eval_integ
260
261{- | Evaluate the definite integral of a function by interpolating
262within the given dataset. For example:
263
264>>> let xs = [1..10]
265>>> let ys = map (**2) [1..10]
266>>> evaluateIntegralV CSpline (zip xs ys) (2.2, 5.5)
26751.909
268
269To successfully @evaluateIntegral points (a, b)@, the domain (@x@)
270values of @points@ must be monotonically increasing. The integration
271bounds @a@ and @b@ must lie between the smallest and largest values in
272the sampled domain..
273-}
274evaluateIntegral :: InterpolationMethod -- ^ What method to use to interpolate
275 -> [(Double, Double)] -- ^ (domain, range) points sampling the function
276 -> (Double, Double) -- ^ Integration bounds (@a@, @b@)
277 -> Double -- ^ Resulting area
278evaluateIntegral mth pts (a, b) =
279 applyCIntFun "evaluateIntegral" "spline_eval_integ" c_spline_eval_integ
280 mth (fromList xs) (fromList ys) a b
281 where
282 (xs, ys) = unzip pts
diff --git a/packages/gsl/src/Numeric/GSL/gsl-aux.c b/packages/gsl/src/Numeric/GSL/gsl-aux.c
index ffc5c20..e1b189c 100644
--- a/packages/gsl/src/Numeric/GSL/gsl-aux.c
+++ b/packages/gsl/src/Numeric/GSL/gsl-aux.c
@@ -34,6 +34,7 @@
34#include <gsl/gsl_rng.h> 34#include <gsl/gsl_rng.h>
35#include <gsl/gsl_randist.h> 35#include <gsl/gsl_randist.h>
36#include <gsl/gsl_roots.h> 36#include <gsl/gsl_roots.h>
37#include <gsl/gsl_spline.h>
37#include <gsl/gsl_multifit_nlin.h> 38#include <gsl/gsl_multifit_nlin.h>
38#include <string.h> 39#include <string.h>
39#include <stdio.h> 40#include <stdio.h>
@@ -140,6 +141,117 @@ int deriv(int code, double f(double, void*), double x, double h, double * result
140 return 0; 141 return 0;
141} 142}
142 143
144int spline_eval(const double xa[], const double ya[], unsigned int size,
145 double x, int method, double *y) {
146 DEBUGMSG("spline_eval");
147 const gsl_interp_type *T;
148 switch (method) {
149 case 0: { T = gsl_interp_linear; break; }
150 case 1: { T = gsl_interp_polynomial; break; }
151 case 2: { T = gsl_interp_cspline; break; }
152 case 3: { T = gsl_interp_cspline_periodic; break; }
153 case 4: { T = gsl_interp_akima; break; }
154 case 5: { T = gsl_interp_akima_periodic; break; }
155 default: ERROR(BAD_CODE);
156 }
157
158 gsl_spline *spline = gsl_spline_alloc(T, size);
159 if (NULL == spline) ERROR(MEM);
160 const int initres = gsl_spline_init(spline, xa, ya, size);
161 CHECK(initres,initres);
162 gsl_interp_accel *acc = gsl_interp_accel_alloc();
163 if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); };
164
165 const int res = gsl_spline_eval_e(spline, x, acc, y);
166 CHECK(res,res);
167 gsl_interp_accel_free(acc);
168 gsl_spline_free(spline);
169 OK
170}
171
172int spline_eval_deriv(const double xa[], const double ya[], unsigned int size,
173 double x, int method, double *y) {
174 DEBUGMSG("spline_eval_deriv");
175 const gsl_interp_type *T;
176 switch (method) {
177 case 0: { T = gsl_interp_linear; break; }
178 case 1: { T = gsl_interp_polynomial; break; }
179 case 2: { T = gsl_interp_cspline; break; }
180 case 3: { T = gsl_interp_cspline_periodic; break; }
181 case 4: { T = gsl_interp_akima; break; }
182 case 5: { T = gsl_interp_akima_periodic; break; }
183 default: ERROR(BAD_CODE);
184 }
185
186 gsl_spline *spline = gsl_spline_alloc(T, size);
187 if (NULL == spline) ERROR(MEM);
188 const int initres = gsl_spline_init(spline, xa, ya, size);
189 CHECK(initres,initres);
190 gsl_interp_accel *acc = gsl_interp_accel_alloc();
191 if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); };
192
193 const int res = gsl_spline_eval_deriv_e(spline, x, acc, y);
194 CHECK(res,res);
195 gsl_interp_accel_free(acc);
196 gsl_spline_free(spline);
197 OK
198}
199
200int spline_eval_deriv2(const double xa[], const double ya[], unsigned int size,
201 double x, int method, double *y) {
202 DEBUGMSG("spline_eval_deriv2");
203 const gsl_interp_type *T;
204 switch (method) {
205 case 0: { T = gsl_interp_linear; break; }
206 case 1: { T = gsl_interp_polynomial; break; }
207 case 2: { T = gsl_interp_cspline; break; }
208 case 3: { T = gsl_interp_cspline_periodic; break; }
209 case 4: { T = gsl_interp_akima; break; }
210 case 5: { T = gsl_interp_akima_periodic; break; }
211 default: ERROR(BAD_CODE);
212 }
213
214 gsl_spline *spline = gsl_spline_alloc(T, size);
215 if (NULL == spline) ERROR(MEM);
216 const int initres = gsl_spline_init(spline, xa, ya, size);
217 CHECK(initres,initres);
218 gsl_interp_accel *acc = gsl_interp_accel_alloc();
219 if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); };
220
221 const int res = gsl_spline_eval_deriv2_e(spline, x, acc, y);
222 CHECK(res,res);
223 gsl_interp_accel_free(acc);
224 gsl_spline_free(spline);
225 OK
226}
227
228int spline_eval_integ(const double xa[], const double ya[], unsigned int size,
229 double a, double b, int method, double *y) {
230 DEBUGMSG("spline_eval_integ");
231 const gsl_interp_type *T;
232 switch (method) {
233 case 0: { T = gsl_interp_linear; break; }
234 case 1: { T = gsl_interp_polynomial; break; }
235 case 2: { T = gsl_interp_cspline; break; }
236 case 3: { T = gsl_interp_cspline_periodic; break; }
237 case 4: { T = gsl_interp_akima; break; }
238 case 5: { T = gsl_interp_akima_periodic; break; }
239 default: ERROR(BAD_CODE);
240 }
241
242 gsl_spline *spline = gsl_spline_alloc(T, size);
243 if (NULL == spline) ERROR(MEM);
244 const int initres = gsl_spline_init(spline, xa, ya, size);
245 CHECK(initres,initres);
246 gsl_interp_accel *acc = gsl_interp_accel_alloc();
247 if (NULL == acc) { gsl_spline_free(spline); ERROR(MEM); };
248
249 const int res = gsl_spline_eval_integ_e(spline, a, b, acc, y);
250 CHECK(res,res);
251 gsl_interp_accel_free(acc);
252 gsl_spline_free(spline);
253 OK
254}
143 255
144int integrate_qng(double f(double, void*), double a, double b, double aprec, double prec, 256int integrate_qng(double f(double, void*), double a, double b, double aprec, double prec,
145 double *result, double*error) { 257 double *result, double*error) {