summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src/Numeric/LinearAlgebra
diff options
context:
space:
mode:
Diffstat (limited to 'packages/hmatrix/src/Numeric/LinearAlgebra')
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs290
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs114
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{- |
4Module : Numeric.LinearAlgebra.Util
5Copyright : (c) Alberto Ruiz 2013
6License : GPL
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10
11-}
12-----------------------------------------------------------------------------
13{-# OPTIONS_HADDOCK hide #-}
14
15module 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
58import Numeric.Container
59import Numeric.IO
60import Numeric.LinearAlgebra.Algorithms hiding (i)
61import Numeric.Matrix()
62import Numeric.Vector()
63
64import Numeric.LinearAlgebra.Util.Convolution
65import 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
712x2
720.33333 0.00000
730.00000 0.33333
74
75-}
76disp :: Int -> Matrix Double -> IO ()
77
78disp 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-}
90diagl :: [Double] -> Matrix Double
91diagl = diag . fromList
92
93-- | a real matrix of zeros
94zeros :: Int -- ^ rows
95 -> Int -- ^ columns
96 -> Matrix Double
97zeros r c = konst 0 (r,c)
98
99-- | a real matrix of ones
100ones :: Int -- ^ rows
101 -> Int -- ^ columns
102 -> Matrix Double
103ones r c = konst 1 (r,c)
104
105-- | concatenation of real vectors
106infixl 3 &
107(&) :: Vector Double -> Vector Double -> Vector Double
108a & 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-}
121infixl 3 ¦
122(¦) :: Matrix Double -> Matrix Double -> Matrix Double
123a ¦ b = fromBlocks [[a,b]]
124
125-- | vertical concatenation of real matrices
126--
127-- (unicode 0x2014, em dash)
128(——) :: Matrix Double -> Matrix Double -> Matrix Double
129infixl 2 ——
130a —— b = fromBlocks [[a],[b]]
131
132(#) :: Matrix Double -> Matrix Double -> Matrix Double
133infixl 2 #
134a # b = fromBlocks [[a],[b]]
135
136-- | create a single row real matrix from a list
137row :: [Double] -> Matrix Double
138row = asRow . fromList
139
140-- | create a single column real matrix from a list
141col :: [Double] -> Matrix Double
142col = 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-}
153infixl 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-}
168infixl 9 ¿
169(¿) :: Element t => Matrix t -> [Int] -> Matrix t
170(¿)= flip extractColumns
171
172
173cross :: Vector Double -> Vector Double -> Vector Double
174-- ^ cross product (for three-element real vectors)
175cross 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
184norm :: Vector Double -> Double
185-- ^ 2-norm of real vector
186norm = pnorm PNorm2
187
188
189-- | Obtains a vector in the same direction with 2-norm=1
190unitary :: Vector Double -> Vector Double
191unitary v = v / scalar (norm v)
192
193-- | ('rows' &&& 'cols')
194size :: Matrix t -> (Int, Int)
195size m = (rows m, cols m)
196
197-- | trans . inv
198mt :: Matrix Double -> Matrix Double
199mt = 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-}
212meanCov :: Matrix Double -> (Vector Double, Matrix Double)
213meanCov 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)
225pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double
226pairwiseD2 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
241rowOuters :: Matrix Double -> Matrix Double -> Matrix Double
242rowOuters 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
250null1 :: Matrix Double -> Vector Double
251null1 = last . toColumns . snd . rightSV
252
253-- | solution of overconstrained homogeneous symmetric linear system
254null1sym :: Matrix Double -> Vector Double
255null1sym = last . toColumns . snd . eigSH'
256
257--------------------------------------------------------------------------------
258
259vec :: Element t => Matrix t -> Vector t
260-- ^ stacking of columns
261vec = flatten . trans
262
263
264vech :: Element t => Matrix t -> Vector t
265-- ^ half-vectorization (of the lower triangular part)
266vech m = vjoin . zipWith f [0..] . toColumns $ m
267 where
268 f k v = subVector k (dim v - k) v
269
270
271dup :: (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)
273dup 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
284vtrans :: Element t => Int -> Matrix t -> Matrix t
285-- ^ generalized \"vector\" transposition: @'vtrans' 1 == 'trans'@, and @'vtrans' ('rows' m) m == 'asColumn' ('vec' m)@
286vtrans 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{- |
4Module : Numeric.LinearAlgebra.Util.Convolution
5Copyright : (c) Alberto Ruiz 2012
6License : GPL
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10
11-}
12-----------------------------------------------------------------------------
13
14module Numeric.LinearAlgebra.Util.Convolution(
15 corr, conv, corrMin,
16 corr2, conv2, separable
17) where
18
19import Numeric.LinearAlgebra
20
21
22vectSS :: Element t => Int -> Vector t -> Matrix t
23vectSS n v = fromRows [ subVector k n v | k <- [0 .. dim v - n] ]
24
25
26corr :: Product t => Vector t -- ^ kernel
27 -> Vector t -- ^ source
28 -> Vector t
29{- ^ correlation
30
31>>> corr (fromList[1,2,3]) (fromList [1..10])
32fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0]
33
34-}
35corr 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
39conv :: (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])
43fromList [-1.0,0.0,1.0]
44
45-}
46conv 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
53corrMin :: (Container Vector t, RealElement t, Product t)
54 => Vector t
55 -> Vector t
56 -> Vector t
57-- ^ similar to 'corr', using 'min' instead of (*)
58corrMin 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
66matSS :: Element t => Int -> Matrix t -> [Matrix t]
67matSS 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
75corr2 :: Product a => Matrix a -> Matrix a -> Matrix a
76-- ^ 2D correlation
77corr2 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
90conv2 :: (Num a, Product a, Container Vector a) => Matrix a -> Matrix a -> Matrix a
91-- ^ 2D convolution
92conv2 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
111separable :: Element t => (Vector t -> Vector t) -> Matrix t -> Matrix t
112-- ^ matrix computation implemented as separated vector operations by rows and columns.
113separable f = fromColumns . map f . toColumns . fromRows . map f . toRows
114