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