diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-16 13:35:35 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-16 13:35:35 +0200 |
commit | 1838c4248679b7476bb8716a76171712dc3cd335 (patch) | |
tree | e3909ac3354eaf917bb1ebc5a7366412b6ab8f0f /packages/base/src/Numeric/LinearAlgebra/Util.hs | |
parent | a2d99e7d0e83fcedf3a856cdb927309e28a8eddd (diff) |
linear algebra moved
Diffstat (limited to 'packages/base/src/Numeric/LinearAlgebra/Util.hs')
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Util.hs | 291 |
1 files changed, 291 insertions, 0 deletions
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs new file mode 100644 index 0000000..533c54b --- /dev/null +++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs | |||
@@ -0,0 +1,291 @@ | |||
1 | {-# LANGUAGE FlexibleContexts #-} | ||
2 | ----------------------------------------------------------------------------- | ||
3 | {- | | ||
4 | Module : Numeric.LinearAlgebra.Util | ||
5 | Copyright : (c) Alberto Ruiz 2013 | ||
6 | License : GPL | ||
7 | |||
8 | Maintainer : Alberto Ruiz (aruiz at um dot es) | ||
9 | Stability : provisional | ||
10 | |||
11 | -} | ||
12 | ----------------------------------------------------------------------------- | ||
13 | {-# OPTIONS_HADDOCK hide #-} | ||
14 | |||
15 | module Numeric.LinearAlgebra.Util( | ||
16 | |||
17 | -- * Convenience functions | ||
18 | size, disp, | ||
19 | zeros, ones, | ||
20 | diagl, | ||
21 | row, | ||
22 | col, | ||
23 | (&), (¦), (——), (#), | ||
24 | (?), (¿), | ||
25 | cross, | ||
26 | norm, | ||
27 | unitary, | ||
28 | mt, | ||
29 | pairwiseD2, | ||
30 | meanCov, | ||
31 | rowOuters, | ||
32 | null1, | ||
33 | null1sym, | ||
34 | -- * Convolution | ||
35 | -- ** 1D | ||
36 | corr, conv, corrMin, | ||
37 | -- ** 2D | ||
38 | corr2, conv2, separable, | ||
39 | -- * Tools for the Kronecker product | ||
40 | -- | ||
41 | -- | (see A. Fusiello, A matter of notation: Several uses of the Kronecker product in | ||
42 | -- 3d computer vision, Pattern Recognition Letters 28 (15) (2007) 2127-2132) | ||
43 | |||
44 | -- | ||
45 | -- | @`vec` (a \<> x \<> b) == ('trans' b ` 'kronecker' ` a) \<> 'vec' x@ | ||
46 | vec, | ||
47 | vech, | ||
48 | dup, | ||
49 | vtrans, | ||
50 | {- -- * Plot | ||
51 | mplot, | ||
52 | plot, parametricPlot, | ||
53 | splot, mesh, meshdom, | ||
54 | matrixToPGM, imshow, | ||
55 | gnuplotX, gnuplotpdf, gnuplotWin | ||
56 | -} | ||
57 | ) where | ||
58 | |||
59 | import Numeric.Container | ||
60 | import Data.Packed.IO | ||
61 | import Numeric.LinearAlgebra.Algorithms hiding (i) | ||
62 | import Numeric.Matrix() | ||
63 | import Numeric.Vector() | ||
64 | |||
65 | import Numeric.LinearAlgebra.Util.Convolution | ||
66 | --import Graphics.Plot | ||
67 | |||
68 | |||
69 | {- | print a real matrix with given number of digits after the decimal point | ||
70 | |||
71 | >>> disp 5 $ ident 2 / 3 | ||
72 | 2x2 | ||
73 | 0.33333 0.00000 | ||
74 | 0.00000 0.33333 | ||
75 | |||
76 | -} | ||
77 | disp :: Int -> Matrix Double -> IO () | ||
78 | |||
79 | disp n = putStrLn . dispf n | ||
80 | |||
81 | |||
82 | {- | create a real diagonal matrix from a list | ||
83 | |||
84 | >>> diagl [1,2,3] | ||
85 | (3><3) | ||
86 | [ 1.0, 0.0, 0.0 | ||
87 | , 0.0, 2.0, 0.0 | ||
88 | , 0.0, 0.0, 3.0 ] | ||
89 | |||
90 | -} | ||
91 | diagl :: [Double] -> Matrix Double | ||
92 | diagl = diag . fromList | ||
93 | |||
94 | -- | a real matrix of zeros | ||
95 | zeros :: Int -- ^ rows | ||
96 | -> Int -- ^ columns | ||
97 | -> Matrix Double | ||
98 | zeros r c = konst 0 (r,c) | ||
99 | |||
100 | -- | a real matrix of ones | ||
101 | ones :: Int -- ^ rows | ||
102 | -> Int -- ^ columns | ||
103 | -> Matrix Double | ||
104 | ones r c = konst 1 (r,c) | ||
105 | |||
106 | -- | concatenation of real vectors | ||
107 | infixl 3 & | ||
108 | (&) :: Vector Double -> Vector Double -> Vector Double | ||
109 | a & b = vjoin [a,b] | ||
110 | |||
111 | {- | horizontal concatenation of real matrices | ||
112 | |||
113 | (unicode 0x00a6, broken bar) | ||
114 | |||
115 | >>> ident 3 ¦ konst 7 (3,4) | ||
116 | (3><7) | ||
117 | [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0 | ||
118 | , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0 | ||
119 | , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ] | ||
120 | |||
121 | -} | ||
122 | infixl 3 ¦ | ||
123 | (¦) :: Matrix Double -> Matrix Double -> Matrix Double | ||
124 | a ¦ b = fromBlocks [[a,b]] | ||
125 | |||
126 | -- | vertical concatenation of real matrices | ||
127 | -- | ||
128 | -- (unicode 0x2014, em dash) | ||
129 | (——) :: Matrix Double -> Matrix Double -> Matrix Double | ||
130 | infixl 2 —— | ||
131 | a —— b = fromBlocks [[a],[b]] | ||
132 | |||
133 | (#) :: Matrix Double -> Matrix Double -> Matrix Double | ||
134 | infixl 2 # | ||
135 | a # b = fromBlocks [[a],[b]] | ||
136 | |||
137 | -- | create a single row real matrix from a list | ||
138 | row :: [Double] -> Matrix Double | ||
139 | row = asRow . fromList | ||
140 | |||
141 | -- | create a single column real matrix from a list | ||
142 | col :: [Double] -> Matrix Double | ||
143 | col = asColumn . fromList | ||
144 | |||
145 | {- | extract rows | ||
146 | |||
147 | >>> (20><4) [1..] ? [2,1,1] | ||
148 | (3><4) | ||
149 | [ 9.0, 10.0, 11.0, 12.0 | ||
150 | , 5.0, 6.0, 7.0, 8.0 | ||
151 | , 5.0, 6.0, 7.0, 8.0 ] | ||
152 | |||
153 | -} | ||
154 | infixl 9 ? | ||
155 | (?) :: Element t => Matrix t -> [Int] -> Matrix t | ||
156 | (?) = flip extractRows | ||
157 | |||
158 | {- | extract columns | ||
159 | |||
160 | (unicode 0x00bf, inverted question mark, Alt-Gr ?) | ||
161 | |||
162 | >>> (3><4) [1..] ¿ [3,0] | ||
163 | (3><2) | ||
164 | [ 4.0, 1.0 | ||
165 | , 8.0, 5.0 | ||
166 | , 12.0, 9.0 ] | ||
167 | |||
168 | -} | ||
169 | infixl 9 ¿ | ||
170 | (¿) :: Element t => Matrix t -> [Int] -> Matrix t | ||
171 | (¿)= flip extractColumns | ||
172 | |||
173 | |||
174 | cross :: Vector Double -> Vector Double -> Vector Double | ||
175 | -- ^ cross product (for three-element real vectors) | ||
176 | cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3] | ||
177 | | otherwise = error $ "cross ("++show x++") ("++show y++")" | ||
178 | where | ||
179 | [x1,x2,x3] = toList x | ||
180 | [y1,y2,y3] = toList y | ||
181 | z1 = x2*y3-x3*y2 | ||
182 | z2 = x3*y1-x1*y3 | ||
183 | z3 = x1*y2-x2*y1 | ||
184 | |||
185 | norm :: Vector Double -> Double | ||
186 | -- ^ 2-norm of real vector | ||
187 | norm = pnorm PNorm2 | ||
188 | |||
189 | |||
190 | -- | Obtains a vector in the same direction with 2-norm=1 | ||
191 | unitary :: Vector Double -> Vector Double | ||
192 | unitary v = v / scalar (norm v) | ||
193 | |||
194 | -- | ('rows' &&& 'cols') | ||
195 | size :: Matrix t -> (Int, Int) | ||
196 | size m = (rows m, cols m) | ||
197 | |||
198 | -- | trans . inv | ||
199 | mt :: Matrix Double -> Matrix Double | ||
200 | mt = trans . inv | ||
201 | |||
202 | -------------------------------------------------------------------------------- | ||
203 | |||
204 | {- | Compute mean vector and covariance matrix of the rows of a matrix. | ||
205 | |||
206 | >>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3]) | ||
207 | (fromList [4.010341078059521,5.0197204699640405], | ||
208 | (2><2) | ||
209 | [ 1.9862461923890056, -1.0127225830525157e-2 | ||
210 | , -1.0127225830525157e-2, 3.0373954915729318 ]) | ||
211 | |||
212 | -} | ||
213 | meanCov :: Matrix Double -> (Vector Double, Matrix Double) | ||
214 | meanCov x = (med,cov) where | ||
215 | r = rows x | ||
216 | k = 1 / fromIntegral r | ||
217 | med = konst k r `vXm` x | ||
218 | meds = konst 1 r `outer` med | ||
219 | xc = x `sub` meds | ||
220 | cov = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc) | ||
221 | |||
222 | -------------------------------------------------------------------------------- | ||
223 | |||
224 | -- | Matrix of pairwise squared distances of row vectors | ||
225 | -- (using the matrix product trick in blog.smola.org) | ||
226 | pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double | ||
227 | pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y | ||
228 | | otherwise = error $ "pairwiseD2 with different number of columns: " | ||
229 | ++ show (size x) ++ ", " ++ show (size y) | ||
230 | where | ||
231 | ox = one (rows x) | ||
232 | oy = one (rows y) | ||
233 | oc = one (cols x) | ||
234 | one k = constant 1 k | ||
235 | x2 = x * x <> oc | ||
236 | y2 = y * y <> oc | ||
237 | ok = cols x == cols y | ||
238 | |||
239 | -------------------------------------------------------------------------------- | ||
240 | |||
241 | -- | outer products of rows | ||
242 | rowOuters :: Matrix Double -> Matrix Double -> Matrix Double | ||
243 | rowOuters a b = a' * b' | ||
244 | where | ||
245 | a' = kronecker a (ones 1 (cols b)) | ||
246 | b' = kronecker (ones 1 (cols a)) b | ||
247 | |||
248 | -------------------------------------------------------------------------------- | ||
249 | |||
250 | -- | solution of overconstrained homogeneous linear system | ||
251 | null1 :: Matrix Double -> Vector Double | ||
252 | null1 = last . toColumns . snd . rightSV | ||
253 | |||
254 | -- | solution of overconstrained homogeneous symmetric linear system | ||
255 | null1sym :: Matrix Double -> Vector Double | ||
256 | null1sym = last . toColumns . snd . eigSH' | ||
257 | |||
258 | -------------------------------------------------------------------------------- | ||
259 | |||
260 | vec :: Element t => Matrix t -> Vector t | ||
261 | -- ^ stacking of columns | ||
262 | vec = flatten . trans | ||
263 | |||
264 | |||
265 | vech :: Element t => Matrix t -> Vector t | ||
266 | -- ^ half-vectorization (of the lower triangular part) | ||
267 | vech m = vjoin . zipWith f [0..] . toColumns $ m | ||
268 | where | ||
269 | f k v = subVector k (dim v - k) v | ||
270 | |||
271 | |||
272 | dup :: (Num t, Num (Vector t), Element t) => Int -> Matrix t | ||
273 | -- ^ duplication matrix (@'dup' k \<> 'vech' m == 'vec' m@, for symmetric m of 'dim' k) | ||
274 | dup k = trans $ fromRows $ map f es | ||
275 | where | ||
276 | rs = zip [0..] (toRows (ident (k^(2::Int)))) | ||
277 | es = [(i,j) | j <- [0..k-1], i <- [0..k-1], i>=j ] | ||
278 | f (i,j) | i == j = g (k*j + i) | ||
279 | | otherwise = g (k*j + i) + g (k*i + j) | ||
280 | g j = v | ||
281 | where | ||
282 | Just v = lookup j rs | ||
283 | |||
284 | |||
285 | vtrans :: Element t => Int -> Matrix t -> Matrix t | ||
286 | -- ^ generalized \"vector\" transposition: @'vtrans' 1 == 'trans'@, and @'vtrans' ('rows' m) m == 'asColumn' ('vec' m)@ | ||
287 | vtrans p m | r == 0 = fromBlocks . map (map asColumn . takesV (replicate q p)) . toColumns $ m | ||
288 | | otherwise = error $ "vtrans " ++ show p ++ " of matrix with " ++ show (rows m) ++ " rows" | ||
289 | where | ||
290 | (q,r) = divMod (rows m) p | ||
291 | |||