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