summaryrefslogtreecommitdiff
path: root/packages/base/src/Internal
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Internal')
-rw-r--r--packages/base/src/Internal/Container.hs285
1 files changed, 285 insertions, 0 deletions
diff --git a/packages/base/src/Internal/Container.hs b/packages/base/src/Internal/Container.hs
new file mode 100644
index 0000000..216e31e
--- /dev/null
+++ b/packages/base/src/Internal/Container.hs
@@ -0,0 +1,285 @@
1{-# LANGUAGE FlexibleContexts #-}
2{-# LANGUAGE FlexibleInstances #-}
3{-# LANGUAGE MultiParamTypeClasses #-}
4{-# LANGUAGE FunctionalDependencies #-}
5{-# LANGUAGE UndecidableInstances #-}
6
7-----------------------------------------------------------------------------
8-- |
9-- Module : Internal.Container
10-- Copyright : (c) Alberto Ruiz 2010-14
11-- License : BSD3
12-- Maintainer : Alberto Ruiz
13-- Stability : provisional
14--
15-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines.
16--
17-- The 'Container' class is used to define optimized generic functions which work
18-- on 'Vector' and 'Matrix' with real or complex elements.
19--
20-- Some of these functions are also available in the instances of the standard
21-- numeric Haskell classes provided by "Numeric.LinearAlgebra".
22--
23-----------------------------------------------------------------------------
24
25module Internal.Container where
26
27import Internal.Tools
28import Internal.Vector
29import Internal.Matrix
30import Internal.Element
31import Internal.Numeric
32import Data.Complex
33import Internal.Algorithms(Field,linearSolveSVD)
34import Data.Vector.Storable(fromList)
35
36------------------------------------------------------------------
37
38{- | Creates a real vector containing a range of values:
39
40>>> linspace 5 (-3,7::Double)
41fromList [-3.0,-0.5,2.0,4.5,7.0]@
42
43>>> linspace 5 (8,2+i) :: Vector (Complex Double)
44fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0]
45
46Logarithmic spacing can be defined as follows:
47
48@logspace n (a,b) = 10 ** linspace n (a,b)@
49-}
50linspace :: (Fractional e, Container Vector e) => Int -> (e, e) -> Vector e
51linspace 0 _ = fromList[]
52linspace 1 (a,b) = fromList[(a+b)/2]
53linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1]
54 where s = (b-a)/fromIntegral (n-1)
55
56--------------------------------------------------------------------------------
57
58infixl 7 <.>
59-- | An infix synonym for 'dot'
60(<.>) :: Numeric t => Vector t -> Vector t -> t
61(<.>) = dot
62
63
64infixr 8 <·>, #>
65
66{- | infix synonym for 'dot'
67
68>>> vector [1,2,3,4] <·> vector [-2,0,1,1]
695.0
70
71>>> let 𝑖 = 0:+1 :: ℂ
72>>> fromList [1+𝑖,1] <·> fromList [1,1+𝑖]
732.0 :+ 0.0
74
75(the dot symbol "·" is obtained by Alt-Gr .)
76
77-}
78(<·>) :: Numeric t => Vector t -> Vector t -> t
79(<·>) = dot
80
81
82{- | infix synonym for 'app'
83
84>>> let m = (2><3) [1..]
85>>> m
86(2><3)
87 [ 1.0, 2.0, 3.0
88 , 4.0, 5.0, 6.0 ]
89
90>>> let v = vector [10,20,30]
91
92>>> m #> v
93fromList [140.0,320.0]
94
95-}
96(#>) :: Numeric t => Matrix t -> Vector t -> Vector t
97(#>) = mXv
98
99-- | dense matrix-vector product
100app :: Numeric t => Matrix t -> Vector t -> Vector t
101app = (#>)
102
103infixl 8 <#
104-- | dense vector-matrix product
105(<#) :: Numeric t => Vector t -> Matrix t -> Vector t
106(<#) = vXm
107
108--------------------------------------------------------------------------------
109
110class Mul a b c | a b -> c where
111 infixl 7 <>
112 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
113 (<>) :: Product t => a t -> b t -> c t
114
115instance Mul Matrix Matrix Matrix where
116 (<>) = mXm
117
118instance Mul Matrix Vector Vector where
119 (<>) m v = flatten $ m <> asColumn v
120
121instance Mul Vector Matrix Vector where
122 (<>) v m = flatten $ asRow v <> m
123
124--------------------------------------------------------------------------------
125
126{- | Least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
127
128@
129a = (3><2)
130 [ 1.0, 2.0
131 , 2.0, 4.0
132 , 2.0, -1.0 ]
133@
134
135@
136v = vector [13.0,27.0,1.0]
137@
138
139>>> let x = a <\> v
140>>> x
141fromList [3.0799999999999996,5.159999999999999]
142
143>>> a #> x
144fromList [13.399999999999999,26.799999999999997,1.0]
145
146It also admits multiple right-hand sides stored as columns in a matrix.
147
148-}
149infixl 7 <\>
150(<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t
151(<\>) = linSolve
152
153class LSDiv c
154 where
155 linSolve :: Field t => Matrix t -> c t -> c t
156
157instance LSDiv Vector
158 where
159 linSolve m v = flatten (linearSolveSVD m (reshape 1 v))
160
161instance LSDiv Matrix
162 where
163 linSolve = linearSolveSVD
164
165--------------------------------------------------------------------------------
166
167class Konst e d c | d -> c, c -> d
168 where
169 -- |
170 -- >>> konst 7 3 :: Vector Float
171 -- fromList [7.0,7.0,7.0]
172 --
173 -- >>> konst i (3::Int,4::Int)
174 -- (3><4)
175 -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
176 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
177 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ]
178 --
179 konst :: e -> d -> c e
180
181instance Container Vector e => Konst e Int Vector
182 where
183 konst = konst'
184
185instance (Num e, Container Vector e) => Konst e (Int,Int) Matrix
186 where
187 konst = konst'
188
189--------------------------------------------------------------------------------
190
191class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
192 where
193 -- |
194 -- >>> build 5 (**2) :: Vector Double
195 -- fromList [0.0,1.0,4.0,9.0,16.0]
196 --
197 -- Hilbert matrix of order N:
198 --
199 -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
200 -- >>> putStr . dispf 2 $ hilb 3
201 -- 3x3
202 -- 1.00 0.50 0.33
203 -- 0.50 0.33 0.25
204 -- 0.33 0.25 0.20
205 --
206 build :: d -> f -> c e
207
208instance Container Vector e => Build Int (e -> e) Vector e
209 where
210 build = build'
211
212instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e
213 where
214 build = build'
215
216--------------------------------------------------------------------------------
217
218-- @dot u v = 'udot' ('conj' u) v@
219dot :: (Numeric t) => Vector t -> Vector t -> t
220dot u v = udot (conj u) v
221
222--------------------------------------------------------------------------------
223
224optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t
225optimiseMult = mconcat
226
227--------------------------------------------------------------------------------
228
229
230{- | Compute mean vector and covariance matrix of the rows of a matrix.
231
232>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3])
233(fromList [4.010341078059521,5.0197204699640405],
234(2><2)
235 [ 1.9862461923890056, -1.0127225830525157e-2
236 , -1.0127225830525157e-2, 3.0373954915729318 ])
237
238-}
239meanCov :: Matrix Double -> (Vector Double, Matrix Double)
240meanCov x = (med,cov) where
241 r = rows x
242 k = 1 / fromIntegral r
243 med = konst k r `vXm` x
244 meds = konst 1 r `outer` med
245 xc = x `sub` meds
246 cov = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc)
247
248--------------------------------------------------------------------------------
249
250class ( Container Vector t
251 , Container Matrix t
252 , Konst t Int Vector
253 , Konst t (Int,Int) Matrix
254 , Product t
255 ) => Numeric t
256
257instance Numeric Double
258instance Numeric (Complex Double)
259instance Numeric Float
260instance Numeric (Complex Float)
261instance Numeric I
262
263--------------------------------------------------------------------------------
264
265sortVector :: (Ord t, Element t) => Vector t -> Vector t
266sortVector = sortV
267
268sortIndex :: (Ord t, Element t) => Vector t -> Vector I
269sortIndex = sortI
270
271ccompare :: (Ord t, Container c t) => c t -> c t -> c I
272ccompare = ccompare'
273
274cselect :: (Container c t) => c I -> c t -> c t -> c t -> c t
275cselect = cselect'
276
277remap :: Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t
278remap i j m
279 | minElement i >= 0 && maxElement i < fromIntegral (rows m) &&
280 minElement j >= 0 && maxElement j < fromIntegral (cols m) = remapM i' j' m
281 | otherwise = error $ "out of range index in rmap"
282 where
283 [i',j'] = conformMs [i,j]
284
285