diff options
Diffstat (limited to 'packages/hmatrix/src/Numeric/LinearAlgebra')
-rw-r--r-- | packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs | 290 | ||||
-rw-r--r-- | packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs | 114 |
2 files changed, 0 insertions, 404 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 | |||
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs deleted file mode 100644 index 82de476..0000000 --- a/packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs +++ /dev/null | |||
@@ -1,114 +0,0 @@ | |||
1 | {-# LANGUAGE FlexibleContexts #-} | ||
2 | ----------------------------------------------------------------------------- | ||
3 | {- | | ||
4 | Module : Numeric.LinearAlgebra.Util.Convolution | ||
5 | Copyright : (c) Alberto Ruiz 2012 | ||
6 | License : GPL | ||
7 | |||
8 | Maintainer : Alberto Ruiz (aruiz at um dot es) | ||
9 | Stability : provisional | ||
10 | |||
11 | -} | ||
12 | ----------------------------------------------------------------------------- | ||
13 | |||
14 | module Numeric.LinearAlgebra.Util.Convolution( | ||
15 | corr, conv, corrMin, | ||
16 | corr2, conv2, separable | ||
17 | ) where | ||
18 | |||
19 | import Numeric.LinearAlgebra | ||
20 | |||
21 | |||
22 | vectSS :: Element t => Int -> Vector t -> Matrix t | ||
23 | vectSS n v = fromRows [ subVector k n v | k <- [0 .. dim v - n] ] | ||
24 | |||
25 | |||
26 | corr :: Product t => Vector t -- ^ kernel | ||
27 | -> Vector t -- ^ source | ||
28 | -> Vector t | ||
29 | {- ^ correlation | ||
30 | |||
31 | >>> corr (fromList[1,2,3]) (fromList [1..10]) | ||
32 | fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0] | ||
33 | |||
34 | -} | ||
35 | corr ker v | dim ker <= dim v = vectSS (dim ker) v <> ker | ||
36 | | otherwise = error $ "corr: dim kernel ("++show (dim ker)++") > dim vector ("++show (dim v)++")" | ||
37 | |||
38 | |||
39 | conv :: (Product t, Num t) => Vector t -> Vector t -> Vector t | ||
40 | {- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product) | ||
41 | |||
42 | >>> conv (fromList[1,1]) (fromList [-1,1]) | ||
43 | fromList [-1.0,0.0,1.0] | ||
44 | |||
45 | -} | ||
46 | conv ker v = corr ker' v' | ||
47 | where | ||
48 | ker' = (flatten.fliprl.asRow) ker | ||
49 | v' | dim ker > 1 = vjoin [z,v,z] | ||
50 | | otherwise = v | ||
51 | z = constant 0 (dim ker -1) | ||
52 | |||
53 | corrMin :: (Container Vector t, RealElement t, Product t) | ||
54 | => Vector t | ||
55 | -> Vector t | ||
56 | -> Vector t | ||
57 | -- ^ similar to 'corr', using 'min' instead of (*) | ||
58 | corrMin ker v = minEvery ss (asRow ker) <> ones | ||
59 | where | ||
60 | minEvery a b = cond a b a a b | ||
61 | ss = vectSS (dim ker) v | ||
62 | ones = konst 1 (dim ker) | ||
63 | |||
64 | |||
65 | |||
66 | matSS :: Element t => Int -> Matrix t -> [Matrix t] | ||
67 | matSS dr m = map (reshape c) [ subVector (k*c) n v | k <- [0 .. r - dr] ] | ||
68 | where | ||
69 | v = flatten m | ||
70 | c = cols m | ||
71 | r = rows m | ||
72 | n = dr*c | ||
73 | |||
74 | |||
75 | corr2 :: Product a => Matrix a -> Matrix a -> Matrix a | ||
76 | -- ^ 2D correlation | ||
77 | corr2 ker mat = dims | ||
78 | . concatMap (map (udot ker' . flatten) . matSS c . trans) | ||
79 | . matSS r $ mat | ||
80 | where | ||
81 | r = rows ker | ||
82 | c = cols ker | ||
83 | ker' = flatten (trans ker) | ||
84 | rr = rows mat - r + 1 | ||
85 | rc = cols mat - c + 1 | ||
86 | dims | rr > 0 && rc > 0 = (rr >< rc) | ||
87 | | otherwise = error $ "corr2: dim kernel ("++sz ker++") > dim matrix ("++sz mat++")" | ||
88 | sz m = show (rows m)++"x"++show (cols m) | ||
89 | |||
90 | conv2 :: (Num a, Product a, Container Vector a) => Matrix a -> Matrix a -> Matrix a | ||
91 | -- ^ 2D convolution | ||
92 | conv2 k m = corr2 (fliprl . flipud $ k) pm | ||
93 | where | ||
94 | pm | r == 0 && c == 0 = m | ||
95 | | r == 0 = fromBlocks [[z3,m,z3]] | ||
96 | | c == 0 = fromBlocks [[z2],[m],[z2]] | ||
97 | | otherwise = fromBlocks [[z1,z2,z1] | ||
98 | ,[z3, m,z3] | ||
99 | ,[z1,z2,z1]] | ||
100 | r = rows k - 1 | ||
101 | c = cols k - 1 | ||
102 | h = rows m | ||
103 | w = cols m | ||
104 | z1 = konst 0 (r,c) | ||
105 | z2 = konst 0 (r,w) | ||
106 | z3 = konst 0 (h,c) | ||
107 | |||
108 | -- TODO: could be simplified using future empty arrays | ||
109 | |||
110 | |||
111 | separable :: Element t => (Vector t -> Vector t) -> Matrix t -> Matrix t | ||
112 | -- ^ matrix computation implemented as separated vector operations by rows and columns. | ||
113 | separable f = fromColumns . map f . toColumns . fromRows . map f . toRows | ||
114 | |||