summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src
diff options
context:
space:
mode:
Diffstat (limited to 'packages/hmatrix/src')
-rw-r--r--packages/hmatrix/src/Numeric/HMatrix.hs4
-rw-r--r--packages/hmatrix/src/Numeric/IO.hs120
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra.hs30
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs290
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs114
-rw-r--r--packages/hmatrix/src/Numeric/Random.hs3
6 files changed, 4 insertions, 557 deletions
diff --git a/packages/hmatrix/src/Numeric/HMatrix.hs b/packages/hmatrix/src/Numeric/HMatrix.hs
index 36fcf70..fcd3e02 100644
--- a/packages/hmatrix/src/Numeric/HMatrix.hs
+++ b/packages/hmatrix/src/Numeric/HMatrix.hs
@@ -129,8 +129,8 @@ module Numeric.HMatrix (
129 129
130import Numeric.HMatrix.Data 130import Numeric.HMatrix.Data
131 131
132import Numeric.Matrix() 132--import Numeric.Matrix()
133import Numeric.Vector() 133--import Numeric.Vector()
134import Numeric.Container 134import Numeric.Container
135import Numeric.LinearAlgebra.Algorithms 135import Numeric.LinearAlgebra.Algorithms
136import Numeric.LinearAlgebra.Util 136import Numeric.LinearAlgebra.Util
diff --git a/packages/hmatrix/src/Numeric/IO.hs b/packages/hmatrix/src/Numeric/IO.hs
index 58fa2b4..08d812b 100644
--- a/packages/hmatrix/src/Numeric/IO.hs
+++ b/packages/hmatrix/src/Numeric/IO.hs
@@ -20,128 +20,10 @@ module Numeric.IO (
20) where 20) where
21 21
22import Data.Packed 22import Data.Packed
23import Data.Packed.Development 23import Data.Packed.IO
24import System.Process(readProcess) 24import System.Process(readProcess)
25import Text.Printf(printf)
26import Data.List(intersperse)
27import Data.Complex
28import Numeric.GSL.Vector 25import Numeric.GSL.Vector
29 26
30{- | Creates a string from a matrix given a separator and a function to show each entry. Using
31this function the user can easily define any desired display function:
32
33@import Text.Printf(printf)@
34
35@disp = putStr . format \" \" (printf \"%.2f\")@
36
37-}
38format :: (Element t) => String -> (t -> String) -> Matrix t -> String
39format sep f m = table sep . map (map f) . toLists $ m
40
41{- | Show a matrix with \"autoscaling\" and a given number of decimal places.
42
43>>> putStr . disps 2 $ 120 * (3><4) [1..]
443x4 E3
45 0.12 0.24 0.36 0.48
46 0.60 0.72 0.84 0.96
47 1.08 1.20 1.32 1.44
48
49-}
50disps :: Int -> Matrix Double -> String
51disps d x = sdims x ++ " " ++ formatScaled d x
52
53{- | Show a matrix with a given number of decimal places.
54
55>>> dispf 2 (1/3 + ident 3)
56"3x3\n1.33 0.33 0.33\n0.33 1.33 0.33\n0.33 0.33 1.33\n"
57
58>>> putStr . dispf 2 $ (3><4)[1,1.5..]
593x4
601.00 1.50 2.00 2.50
613.00 3.50 4.00 4.50
625.00 5.50 6.00 6.50
63
64>>> putStr . unlines . tail . lines . dispf 2 . asRow $ linspace 10 (0,1)
650.00 0.11 0.22 0.33 0.44 0.56 0.67 0.78 0.89 1.00
66
67-}
68dispf :: Int -> Matrix Double -> String
69dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x
70
71sdims x = show (rows x) ++ "x" ++ show (cols x)
72
73formatFixed d x = format " " (printf ("%."++show d++"f")) $ x
74
75isInt = all lookslikeInt . toList . flatten
76
77formatScaled dec t = "E"++show o++"\n" ++ ss
78 where ss = format " " (printf fmt. g) t
79 g x | o >= 0 = x/10^(o::Int)
80 | otherwise = x*10^(-o)
81 o | rows t == 0 || cols t == 0 = 0
82 | otherwise = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t
83 fmt = '%':show (dec+3) ++ '.':show dec ++"f"
84
85{- | Show a vector using a function for showing matrices.
86
87>>> putStr . vecdisp (dispf 2) $ linspace 10 (0,1)
8810 |> 0.00 0.11 0.22 0.33 0.44 0.56 0.67 0.78 0.89 1.00
89
90-}
91vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String
92vecdisp f v
93 = ((show (dim v) ++ " |> ") ++) . (++"\n")
94 . unwords . lines . tail . dropWhile (not . (`elem` " \n"))
95 . f . trans . reshape 1
96 $ v
97
98{- | Tool to display matrices with latex syntax.
99
100>>> latexFormat "bmatrix" (dispf 2 $ ident 2)
101"\\begin{bmatrix}\n1 & 0\n\\\\\n0 & 1\n\\end{bmatrix}"
102
103-}
104latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc.
105 -> String -- ^ Formatted matrix, with elements separated by spaces and newlines
106 -> String
107latexFormat del tab = "\\begin{"++del++"}\n" ++ f tab ++ "\\end{"++del++"}"
108 where f = unlines . intersperse "\\\\" . map unwords . map (intersperse " & " . words) . tail . lines
109
110-- | Pretty print a complex number with at most n decimal digits.
111showComplex :: Int -> Complex Double -> String
112showComplex d (a:+b)
113 | isZero a && isZero b = "0"
114 | isZero b = sa
115 | isZero a && isOne b = s2++"i"
116 | isZero a = sb++"i"
117 | isOne b = sa++s3++"i"
118 | otherwise = sa++s1++sb++"i"
119 where
120 sa = shcr d a
121 sb = shcr d b
122 s1 = if b<0 then "" else "+"
123 s2 = if b<0 then "-" else ""
124 s3 = if b<0 then "-" else "+"
125
126shcr d a | lookslikeInt a = printf "%.0f" a
127 | otherwise = printf ("%."++show d++"f") a
128
129
130lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
131 where shx = show x
132
133isZero x = show x `elem` ["0.0","-0.0"]
134isOne x = show x `elem` ["1.0","-1.0"]
135
136-- | Pretty print a complex matrix with at most n decimal digits.
137dispcf :: Int -> Matrix (Complex Double) -> String
138dispcf d m = sdims m ++ "\n" ++ format " " (showComplex d) m
139
140--------------------------------------------------------------------
141
142-- | reads a matrix from a string containing a table of numbers.
143readMatrix :: String -> Matrix Double
144readMatrix = fromLists . map (map read). map words . filter (not.null) . lines
145 27
146{- | obtains the number of rows and columns in an ASCII data file 28{- | obtains the number of rows and columns in an ASCII data file
147 (provisionally using unix's wc). 29 (provisionally using unix's wc).
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra.hs b/packages/hmatrix/src/Numeric/LinearAlgebra.hs
deleted file mode 100644
index 1db860c..0000000
--- a/packages/hmatrix/src/Numeric/LinearAlgebra.hs
+++ /dev/null
@@ -1,30 +0,0 @@
1-----------------------------------------------------------------------------
2{- |
3Module : Numeric.LinearAlgebra
4Copyright : (c) Alberto Ruiz 2006-10
5License : GPL-style
6
7Maintainer : Alberto Ruiz (aruiz at um dot es)
8Stability : provisional
9Portability : uses ffi
10
11This module reexports all normally required functions for Linear Algebra applications.
12
13It also provides instances of standard classes 'Show', 'Read', 'Eq',
14'Num', 'Fractional', and 'Floating' for 'Vector' and 'Matrix'.
15In arithmetic operations one-component vectors and matrices automatically
16expand to match the dimensions of the other operand.
17
18-}
19-----------------------------------------------------------------------------
20{-# OPTIONS_HADDOCK hide #-}
21
22module Numeric.LinearAlgebra (
23 module Numeric.Container,
24 module Numeric.LinearAlgebra.Algorithms
25) where
26
27import Numeric.Container
28import Numeric.LinearAlgebra.Algorithms
29import Numeric.Matrix()
30import Numeric.Vector()
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
diff --git a/packages/hmatrix/src/Numeric/Random.hs b/packages/hmatrix/src/Numeric/Random.hs
index 7bf9e8b..0e722fd 100644
--- a/packages/hmatrix/src/Numeric/Random.hs
+++ b/packages/hmatrix/src/Numeric/Random.hs
@@ -21,8 +21,7 @@ module Numeric.Random (
21) where 21) where
22 22
23import Numeric.GSL.Vector 23import Numeric.GSL.Vector
24import Data.Packed 24import Numeric.Container
25import Data.Packed.Numeric
26import Numeric.LinearAlgebra.Algorithms 25import Numeric.LinearAlgebra.Algorithms
27import System.Random(randomIO) 26import System.Random(randomIO)
28 27
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-05 23:05:31 +0000