summaryrefslogtreecommitdiff
path: root/packages/base/src/Data/Packed/Numeric.hs
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Data/Packed/Numeric.hs')
-rw-r--r--packages/base/src/Data/Packed/Numeric.hs288
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
27module 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
70import Data.Packed
71import Data.Packed.Internal.Numeric
72import Data.Complex
73import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
74import Data.Monoid(Monoid(mconcat))
75import Data.Packed.IO
76import Numeric.LinearAlgebra.Random
77
78------------------------------------------------------------------
79
80{- | Creates a real vector containing a range of values:
81
82>>> linspace 5 (-3,7::Double)
83fromList [-3.0,-0.5,2.0,4.5,7.0]@
84
85>>> linspace 5 (8,2+i) :: Vector (Complex Double)
86fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0]
87
88Logarithmic spacing can be defined as follows:
89
90@logspace n (a,b) = 10 ** linspace n (a,b)@
91-}
92linspace :: (Container Vector e) => Int -> (e, e) -> Vector e
93linspace 0 (a,b) = fromList[(a+b)/2]
94linspace 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
103Examples:
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
115matrix × matrix:
116
117>>> disp 2 (a <.> trans a)
1183x3
119 30 70 110
120 70 174 278
121110 278 446
122
123matrix × vector:
124
125>>> a <.> v
126fromList [3.0,11.0,19.0]
127
128dot product:
129
130>>> u <.> fromList[3,2,1::Double]
13110
132
133For complex vectors the first argument is conjugated:
134
135>>> fromList [1,i] <.> fromList[2*i+1,3]
1361.0 :+ (-1.0)
137
138>>> fromList [1,i,1-i] <.> complex a
139fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0]
140-}
141infixl 7 <.>
142(<.>) :: Contraction a b c => a -> b -> c
143(<.>) = contraction
144
145
146class Contraction a b c | a b -> c
147 where
148 -- | Matrix product, matrix - vector product, and dot product
149 contraction :: a -> b -> c
150
151instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where
152 u `contraction` v = conj u `udot` v
153
154instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where
155 contraction = mXv
156
157instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where
158 contraction v m = (conj v) `vXm` m
159
160instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where
161 contraction = mXm
162
163
164--------------------------------------------------------------------------------
165
166class 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
171instance Mul Matrix Matrix Matrix where
172 (<>) = mXm
173
174instance Mul Matrix Vector Vector where
175 (<>) m v = flatten $ m <> asColumn v
176
177instance 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)
183infixl 7 <\>
184(<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t
185(<\>) = linSolve
186
187class LSDiv c
188 where
189 linSolve :: Field t => Matrix t -> c t -> c t
190
191instance LSDiv Vector
192 where
193 linSolve m v = flatten (linearSolveSVD m (reshape 1 v))
194
195instance LSDiv Matrix
196 where
197 linSolve = linearSolveSVD
198
199--------------------------------------------------------------------------------
200
201class 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
215instance Container Vector e => Konst e Int Vector
216 where
217 konst = konst'
218
219instance Container Vector e => Konst e (Int,Int) Matrix
220 where
221 konst = konst'
222
223--------------------------------------------------------------------------------
224
225class 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
242instance Container Vector e => Build Int (e -> e) Vector e
243 where
244 build = build'
245
246instance 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
254infixl 7 ◇
255(◇) = contraction
256
257-- | dot product: @cdot u v = 'udot' ('conj' u) v@
258dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t
259dot u v = udot (conj u) v
260
261--------------------------------------------------------------------------------
262
263optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t
264optimiseMult = 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-}
278meanCov :: Matrix Double -> (Vector Double, Matrix Double)
279meanCov 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