summaryrefslogtreecommitdiff
path: root/lib/Numeric/LinearAlgebra/Util.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/LinearAlgebra/Util.hs')
-rw-r--r--lib/Numeric/LinearAlgebra/Util.hs295
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{- |
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 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
58import Numeric.Container
59import Numeric.LinearAlgebra.Algorithms hiding (i)
60import Numeric.Matrix()
61import Numeric.Vector()
62
63import System.Random(randomIO)
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-- | pseudorandom matrix with uniform elements between 0 and 1
81randm :: RandDist
82 -> Int -- ^ rows
83 -> Int -- ^ columns
84 -> IO (Matrix Double)
85randm 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
90rand :: Int -> Int -> IO (Matrix Double)
91rand = randm Uniform
92
93{- | pseudorandom matrix with normal elements
94
95>>> x <- randn 3 5
96>>> disp 3 x
973x5
980.386 -1.141 0.491 -0.510 1.512
990.069 -0.919 1.022 -0.181 0.745
1000.313 -0.670 -0.097 -1.575 -0.583
101
102-}
103randn :: Int -> Int -> IO (Matrix Double)
104randn = 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-}
115diagl :: [Double] -> Matrix Double
116diagl = diag . fromList
117
118-- | a real matrix of zeros
119zeros :: Int -- ^ rows
120 -> Int -- ^ columns
121 -> Matrix Double
122zeros r c = konst 0 (r,c)
123
124-- | a real matrix of ones
125ones :: Int -- ^ rows
126 -> Int -- ^ columns
127 -> Matrix Double
128ones r c = konst 1 (r,c)
129
130-- | concatenation of real vectors
131infixl 3 &
132(&) :: Vector Double -> Vector Double -> Vector Double
133a & 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-}
146infixl 3 ¦
147(¦) :: Matrix Double -> Matrix Double -> Matrix Double
148a ¦ b = fromBlocks [[a,b]]
149
150-- | vertical concatenation of real matrices
151--
152-- (unicode 0x2014, em dash)
153(——) :: Matrix Double -> Matrix Double -> Matrix Double
154infixl 2 ——
155a —— b = fromBlocks [[a],[b]]
156
157(#) :: Matrix Double -> Matrix Double -> Matrix Double
158infixl 2 #
159a # b = fromBlocks [[a],[b]]
160
161-- | create a single row real matrix from a list
162row :: [Double] -> Matrix Double
163row = asRow . fromList
164
165-- | create a single column real matrix from a list
166col :: [Double] -> Matrix Double
167col = 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-}
178infixl 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-}
193infixl 9 ¿
194(¿) :: Element t => Matrix t -> [Int] -> Matrix t
195(¿)= flip extractColumns
196
197
198cross :: Vector Double -> Vector Double -> Vector Double
199-- ^ cross product (for three-element real vectors)
200cross 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
209norm :: Vector Double -> Double
210-- ^ 2-norm of real vector
211norm = pnorm PNorm2
212
213
214-- | Obtains a vector in the same direction with 2-norm=1
215unitary :: Vector Double -> Vector Double
216unitary v = v / scalar (norm v)
217
218-- | ('rows' &&& 'cols')
219size :: Matrix t -> (Int, Int)
220size m = (rows m, cols m)
221
222-- | trans . inv
223mt :: Matrix Double -> Matrix Double
224mt = trans . inv
225
226----------------------------------------------------------------------
227
228-- | Matrix of pairwise squared distances of row vectors
229-- (using the matrix product trick in blog.smola.org)
230pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double
231pairwiseD2 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
246rowOuters :: Matrix Double -> Matrix Double -> Matrix Double
247rowOuters 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
255null1 :: Matrix Double -> Vector Double
256null1 = last . toColumns . snd . rightSV
257
258-- | solution of overconstrained homogeneous symmetric linear system
259null1sym :: Matrix Double -> Vector Double
260null1sym = last . toColumns . snd . eigSH'
261
262--------------------------------------------------------------------------------
263
264vec :: Element t => Matrix t -> Vector t
265-- ^ stacking of columns
266vec = flatten . trans
267
268
269vech :: Element t => Matrix t -> Vector t
270-- ^ half-vectorization (of the lower triangular part)
271vech m = vjoin . zipWith f [0..] . toColumns $ m
272 where
273 f k v = subVector k (dim v - k) v
274
275
276dup :: (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)
278dup 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
289vtrans :: Element t => Int -> Matrix t -> Matrix t
290-- ^ generalized \"vector\" transposition: @'vtrans' 1 == 'trans'@, and @'vtrans' ('rows' m) m == 'asColumn' ('vec' m)@
291vtrans 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