diff options
Diffstat (limited to 'packages/base/src/Data/Packed/Numeric.hs')
-rw-r--r-- | packages/base/src/Data/Packed/Numeric.hs | 288 |
1 files changed, 288 insertions, 0 deletions
diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs new file mode 100644 index 0000000..01cf6c5 --- /dev/null +++ b/packages/base/src/Data/Packed/Numeric.hs | |||
@@ -0,0 +1,288 @@ | |||
1 | {-# LANGUAGE TypeFamilies #-} | ||
2 | {-# LANGUAGE FlexibleContexts #-} | ||
3 | {-# LANGUAGE FlexibleInstances #-} | ||
4 | {-# LANGUAGE MultiParamTypeClasses #-} | ||
5 | {-# LANGUAGE FunctionalDependencies #-} | ||
6 | {-# LANGUAGE UndecidableInstances #-} | ||
7 | |||
8 | ----------------------------------------------------------------------------- | ||
9 | -- | | ||
10 | -- Module : Data.Packed.Numeric | ||
11 | -- Copyright : (c) Alberto Ruiz 2010-14 | ||
12 | -- License : BSD3 | ||
13 | -- Maintainer : Alberto Ruiz | ||
14 | -- Stability : provisional | ||
15 | -- | ||
16 | -- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines. | ||
17 | -- | ||
18 | -- The 'Container' class is used to define optimized generic functions which work | ||
19 | -- on 'Vector' and 'Matrix' with real or complex elements. | ||
20 | -- | ||
21 | -- Some of these functions are also available in the instances of the standard | ||
22 | -- numeric Haskell classes provided by "Numeric.LinearAlgebra". | ||
23 | -- | ||
24 | ----------------------------------------------------------------------------- | ||
25 | {-# OPTIONS_HADDOCK hide #-} | ||
26 | |||
27 | module Data.Packed.Numeric ( | ||
28 | -- * Basic functions | ||
29 | module Data.Packed, | ||
30 | Konst(..), Build(..), | ||
31 | linspace, | ||
32 | diag, ident, | ||
33 | ctrans, | ||
34 | -- * Generic operations | ||
35 | Container(..), | ||
36 | -- add, mul, sub, divide, equal, scaleRecip, addConstant, | ||
37 | scalar, conj, scale, arctan2, cmap, | ||
38 | atIndex, minIndex, maxIndex, minElement, maxElement, | ||
39 | sumElements, prodElements, | ||
40 | step, cond, find, assoc, accum, | ||
41 | Transposable(..), Linear(..), | ||
42 | -- * Matrix product | ||
43 | Product(..), udot, dot, (◇), | ||
44 | Mul(..), | ||
45 | Contraction(..),(<.>), | ||
46 | optimiseMult, | ||
47 | mXm,mXv,vXm,LSDiv,(<\>), | ||
48 | outer, kronecker, | ||
49 | -- * Random numbers | ||
50 | RandDist(..), | ||
51 | randomVector, | ||
52 | gaussianSample, | ||
53 | uniformSample, | ||
54 | meanCov, | ||
55 | -- * Element conversion | ||
56 | Convert(..), | ||
57 | Complexable(), | ||
58 | RealElement(), | ||
59 | |||
60 | RealOf, ComplexOf, SingleOf, DoubleOf, | ||
61 | |||
62 | IndexOf, | ||
63 | module Data.Complex, | ||
64 | -- * IO | ||
65 | module Data.Packed.IO, | ||
66 | -- * Misc | ||
67 | Testable(..) | ||
68 | ) where | ||
69 | |||
70 | import Data.Packed | ||
71 | import Data.Packed.Internal.Numeric | ||
72 | import Data.Complex | ||
73 | import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD) | ||
74 | import Data.Monoid(Monoid(mconcat)) | ||
75 | import Data.Packed.IO | ||
76 | import Numeric.LinearAlgebra.Random | ||
77 | |||
78 | ------------------------------------------------------------------ | ||
79 | |||
80 | {- | Creates a real vector containing a range of values: | ||
81 | |||
82 | >>> linspace 5 (-3,7::Double) | ||
83 | fromList [-3.0,-0.5,2.0,4.5,7.0]@ | ||
84 | |||
85 | >>> linspace 5 (8,2+i) :: Vector (Complex Double) | ||
86 | fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0] | ||
87 | |||
88 | Logarithmic spacing can be defined as follows: | ||
89 | |||
90 | @logspace n (a,b) = 10 ** linspace n (a,b)@ | ||
91 | -} | ||
92 | linspace :: (Container Vector e) => Int -> (e, e) -> Vector e | ||
93 | linspace 0 (a,b) = fromList[(a+b)/2] | ||
94 | linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1] | ||
95 | where s = (b-a)/fromIntegral (n-1) | ||
96 | |||
97 | -------------------------------------------------------- | ||
98 | |||
99 | {- | Matrix product, matrix - vector product, and dot product (equivalent to 'contraction') | ||
100 | |||
101 | (This operator can also be written using the unicode symbol ◇ (25c7).) | ||
102 | |||
103 | Examples: | ||
104 | |||
105 | >>> let a = (3><4) [1..] :: Matrix Double | ||
106 | >>> let v = fromList [1,0,2,-1] :: Vector Double | ||
107 | >>> let u = fromList [1,2,3] :: Vector Double | ||
108 | |||
109 | >>> a | ||
110 | (3><4) | ||
111 | [ 1.0, 2.0, 3.0, 4.0 | ||
112 | , 5.0, 6.0, 7.0, 8.0 | ||
113 | , 9.0, 10.0, 11.0, 12.0 ] | ||
114 | |||
115 | matrix × matrix: | ||
116 | |||
117 | >>> disp 2 (a <.> trans a) | ||
118 | 3x3 | ||
119 | 30 70 110 | ||
120 | 70 174 278 | ||
121 | 110 278 446 | ||
122 | |||
123 | matrix × vector: | ||
124 | |||
125 | >>> a <.> v | ||
126 | fromList [3.0,11.0,19.0] | ||
127 | |||
128 | dot product: | ||
129 | |||
130 | >>> u <.> fromList[3,2,1::Double] | ||
131 | 10 | ||
132 | |||
133 | For complex vectors the first argument is conjugated: | ||
134 | |||
135 | >>> fromList [1,i] <.> fromList[2*i+1,3] | ||
136 | 1.0 :+ (-1.0) | ||
137 | |||
138 | >>> fromList [1,i,1-i] <.> complex a | ||
139 | fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0] | ||
140 | -} | ||
141 | infixl 7 <.> | ||
142 | (<.>) :: Contraction a b c => a -> b -> c | ||
143 | (<.>) = contraction | ||
144 | |||
145 | |||
146 | class Contraction a b c | a b -> c | ||
147 | where | ||
148 | -- | Matrix product, matrix - vector product, and dot product | ||
149 | contraction :: a -> b -> c | ||
150 | |||
151 | instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where | ||
152 | u `contraction` v = conj u `udot` v | ||
153 | |||
154 | instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where | ||
155 | contraction = mXv | ||
156 | |||
157 | instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where | ||
158 | contraction v m = (conj v) `vXm` m | ||
159 | |||
160 | instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where | ||
161 | contraction = mXm | ||
162 | |||
163 | |||
164 | -------------------------------------------------------------------------------- | ||
165 | |||
166 | class Mul a b c | a b -> c where | ||
167 | infixl 7 <> | ||
168 | -- | Matrix-matrix, matrix-vector, and vector-matrix products. | ||
169 | (<>) :: Product t => a t -> b t -> c t | ||
170 | |||
171 | instance Mul Matrix Matrix Matrix where | ||
172 | (<>) = mXm | ||
173 | |||
174 | instance Mul Matrix Vector Vector where | ||
175 | (<>) m v = flatten $ m <> asColumn v | ||
176 | |||
177 | instance Mul Vector Matrix Vector where | ||
178 | (<>) v m = flatten $ asRow v <> m | ||
179 | |||
180 | -------------------------------------------------------------------------------- | ||
181 | |||
182 | -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD) | ||
183 | infixl 7 <\> | ||
184 | (<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t | ||
185 | (<\>) = linSolve | ||
186 | |||
187 | class LSDiv c | ||
188 | where | ||
189 | linSolve :: Field t => Matrix t -> c t -> c t | ||
190 | |||
191 | instance LSDiv Vector | ||
192 | where | ||
193 | linSolve m v = flatten (linearSolveSVD m (reshape 1 v)) | ||
194 | |||
195 | instance LSDiv Matrix | ||
196 | where | ||
197 | linSolve = linearSolveSVD | ||
198 | |||
199 | -------------------------------------------------------------------------------- | ||
200 | |||
201 | class Konst e d c | d -> c, c -> d | ||
202 | where | ||
203 | -- | | ||
204 | -- >>> konst 7 3 :: Vector Float | ||
205 | -- fromList [7.0,7.0,7.0] | ||
206 | -- | ||
207 | -- >>> konst i (3::Int,4::Int) | ||
208 | -- (3><4) | ||
209 | -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 | ||
210 | -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 | ||
211 | -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ] | ||
212 | -- | ||
213 | konst :: e -> d -> c e | ||
214 | |||
215 | instance Container Vector e => Konst e Int Vector | ||
216 | where | ||
217 | konst = konst' | ||
218 | |||
219 | instance Container Vector e => Konst e (Int,Int) Matrix | ||
220 | where | ||
221 | konst = konst' | ||
222 | |||
223 | -------------------------------------------------------------------------------- | ||
224 | |||
225 | class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f | ||
226 | where | ||
227 | -- | | ||
228 | -- >>> build 5 (**2) :: Vector Double | ||
229 | -- fromList [0.0,1.0,4.0,9.0,16.0] | ||
230 | -- | ||
231 | -- Hilbert matrix of order N: | ||
232 | -- | ||
233 | -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double | ||
234 | -- >>> putStr . dispf 2 $ hilb 3 | ||
235 | -- 3x3 | ||
236 | -- 1.00 0.50 0.33 | ||
237 | -- 0.50 0.33 0.25 | ||
238 | -- 0.33 0.25 0.20 | ||
239 | -- | ||
240 | build :: d -> f -> c e | ||
241 | |||
242 | instance Container Vector e => Build Int (e -> e) Vector e | ||
243 | where | ||
244 | build = build' | ||
245 | |||
246 | instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e | ||
247 | where | ||
248 | build = build' | ||
249 | |||
250 | -------------------------------------------------------------------------------- | ||
251 | |||
252 | -- | alternative unicode symbol (25c7) for 'contraction' | ||
253 | (◇) :: Contraction a b c => a -> b -> c | ||
254 | infixl 7 ◇ | ||
255 | (◇) = contraction | ||
256 | |||
257 | -- | dot product: @cdot u v = 'udot' ('conj' u) v@ | ||
258 | dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t | ||
259 | dot u v = udot (conj u) v | ||
260 | |||
261 | -------------------------------------------------------------------------------- | ||
262 | |||
263 | optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t | ||
264 | optimiseMult = mconcat | ||
265 | |||
266 | -------------------------------------------------------------------------------- | ||
267 | |||
268 | |||
269 | {- | Compute mean vector and covariance matrix of the rows of a matrix. | ||
270 | |||
271 | >>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3]) | ||
272 | (fromList [4.010341078059521,5.0197204699640405], | ||
273 | (2><2) | ||
274 | [ 1.9862461923890056, -1.0127225830525157e-2 | ||
275 | , -1.0127225830525157e-2, 3.0373954915729318 ]) | ||
276 | |||
277 | -} | ||
278 | meanCov :: Matrix Double -> (Vector Double, Matrix Double) | ||
279 | meanCov x = (med,cov) where | ||
280 | r = rows x | ||
281 | k = 1 / fromIntegral r | ||
282 | med = konst k r `vXm` x | ||
283 | meds = konst 1 r `outer` med | ||
284 | xc = x `sub` meds | ||
285 | cov = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc) | ||
286 | |||
287 | -------------------------------------------------------------------------------- | ||
288 | |||