diff options
author | Alberto Ruiz <aruiz@um.es> | 2007-06-09 12:10:58 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2007-06-09 12:10:58 +0000 |
commit | 7931a9b18ea84ed5f49e2803ba596f190567d9d8 (patch) | |
tree | 64a08a62b2bffcf48becbab03933f3c7b4527a73 | |
parent | e21f42f742959ec9452add9b6c6e08d30d9584ed (diff) |
more tests
-rw-r--r-- | examples/tests.hs | 149 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Matrix.hs | 17 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Vector.hs | 1 |
3 files changed, 163 insertions, 4 deletions
diff --git a/examples/tests.hs b/examples/tests.hs new file mode 100644 index 0000000..f167b92 --- /dev/null +++ b/examples/tests.hs | |||
@@ -0,0 +1,149 @@ | |||
1 | -- | ||
2 | -- QuickCheck tests | ||
3 | -- | ||
4 | |||
5 | ----------------------------------------------------------------------------- | ||
6 | |||
7 | import Data.Packed.Internal.Vector | ||
8 | import Data.Packed.Internal.Matrix | ||
9 | import LAPACK | ||
10 | import Test.QuickCheck | ||
11 | import Complex | ||
12 | |||
13 | {- | ||
14 | -- Bravo por quickCheck! | ||
15 | |||
16 | pinvProp1 tol m = (rank m == cols m) ==> pinv m <> m ~~ ident (cols m) | ||
17 | where infix 2 ~~ | ||
18 | (~~) = approxEqual tol | ||
19 | |||
20 | pinvProp2 tol m = 0 < r && r <= c ==> (r==c) `trivial` (m <> pinv m <> m ~~ m) | ||
21 | where r = rank m | ||
22 | c = cols m | ||
23 | infix 2 ~~ | ||
24 | (~~) = approxEqual tol | ||
25 | |||
26 | nullspaceProp tol m = cr > 0 ==> m <> nt ~~ zeros | ||
27 | where nt = trans (nullspace m) | ||
28 | cr = corank m | ||
29 | r = rows m | ||
30 | zeros = create [r,cr] $ replicate (r*cr) 0 | ||
31 | |||
32 | -} | ||
33 | |||
34 | r >< c = f where | ||
35 | f l | dim v == r*c = matrixFromVector RowMajor c v | ||
36 | | otherwise = error $ "inconsistent list size = " | ||
37 | ++show (dim v) ++"in ("++show r++"><"++show c++")" | ||
38 | where v = fromList l | ||
39 | |||
40 | r >|< c = f where | ||
41 | f l | dim v == r*c = matrixFromVector ColumnMajor c v | ||
42 | | otherwise = error $ "inconsistent list size = " | ||
43 | ++show (dim v) ++"in ("++show r++"><"++show c++")" | ||
44 | where v = fromList l | ||
45 | |||
46 | ac = (2><3) [1 .. 6::Double] | ||
47 | bc = (3><4) [7 .. 18::Double] | ||
48 | |||
49 | mz = (2 >< 3) [1,2,3,4,5,6:+(1::Double)] | ||
50 | |||
51 | af = (2>|<3) [1,4,2,5,3,6::Double] | ||
52 | bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double] | ||
53 | |||
54 | a |=| b = rows a == rows b && | ||
55 | cols a == cols b && | ||
56 | toList (cdat a) == toList (cdat b) && | ||
57 | toList (fdat a) == toList (fdat b) | ||
58 | |||
59 | aprox fun a b = rows a == rows b && | ||
60 | cols a == cols b && | ||
61 | eps > aproxL fun (toList (t a)) (toList (t b)) | ||
62 | where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat | ||
63 | |||
64 | aproxL fun v1 v2 = sum (zipWith (\a b-> fun (a-b)) v1 v2) / fromIntegral (length v1) | ||
65 | |||
66 | (|~|) = aprox abs | ||
67 | (|~~|) = aprox magnitude | ||
68 | |||
69 | eps = 1E-8::Double | ||
70 | |||
71 | asFortran m = (rows m >|< cols m) $ toList (fdat m) | ||
72 | asC m = (rows m >< cols m) $ toList (cdat m) | ||
73 | |||
74 | mulC a b = multiply RowMajor a b | ||
75 | mulF a b = multiply ColumnMajor a b | ||
76 | |||
77 | cc = mulC ac bf | ||
78 | cf = mulF af bc | ||
79 | |||
80 | r = mulC cc (trans cf) | ||
81 | |||
82 | ident n = diag (constant n 1) | ||
83 | |||
84 | rd = (2><2) | ||
85 | [ 43492.0, 50572.0 | ||
86 | , 102550.0, 119242.0 :: Double] | ||
87 | |||
88 | instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where | ||
89 | arbitrary = do | ||
90 | r <- arbitrary | ||
91 | i <- arbitrary | ||
92 | return (r:+i) | ||
93 | coarbitrary = undefined | ||
94 | |||
95 | instance (Field a, Arbitrary a) => Arbitrary (Matrix a) where | ||
96 | arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) | ||
97 | m <- choose (1,10) | ||
98 | n <- choose (1,10) | ||
99 | l <- vector (m*n) | ||
100 | ctype <- arbitrary | ||
101 | let h = if ctype then (m><n) else (m>|<n) | ||
102 | trMode <- arbitrary | ||
103 | let tr = if trMode then trans else id | ||
104 | return $ tr (h l) | ||
105 | coarbitrary = undefined | ||
106 | |||
107 | data PairM a = PairM (Matrix a) (Matrix a) deriving Show | ||
108 | instance (Num a, Field a, Arbitrary a) => Arbitrary (PairM a) where | ||
109 | arbitrary = do | ||
110 | a <- choose (1,10) | ||
111 | b <- choose (1,10) | ||
112 | c <- choose (1,10) | ||
113 | l1 <- vector (a*b) | ||
114 | l2 <- vector (b*c) | ||
115 | return $ PairM ((a><b) (map fromIntegral (l1::[Int]))) ((b><c) (map fromIntegral (l2::[Int]))) | ||
116 | --return $ PairM ((a><b) l1) ((b><c) l2) | ||
117 | coarbitrary = undefined | ||
118 | |||
119 | type BaseType = Double | ||
120 | |||
121 | |||
122 | svdTestR fun prod m = u <> s <> trans v |~| m | ||
123 | && u <> trans u |~| ident (rows m) | ||
124 | && v <> trans v |~| ident (cols m) | ||
125 | where (u,s,v) = fun m | ||
126 | (<>) = prod | ||
127 | |||
128 | |||
129 | svdTestC fun prod m = u <> s' <> (trans v) |~~| m | ||
130 | && u <> (liftMatrix conj) (trans u) |~~| ident (rows m) | ||
131 | && v <> (liftMatrix conj) (trans v) |~~| ident (cols m) | ||
132 | where (u,s,v) = fun m | ||
133 | (<>) = prod | ||
134 | s' = liftMatrix comp s | ||
135 | |||
136 | comp v = toComplex (v,constant (dim v) 0) | ||
137 | |||
138 | main = do | ||
139 | quickCheck $ \l -> null l || (toList . fromList) l == (l :: [BaseType]) | ||
140 | quickCheck $ \m -> m |=| asC (m :: Matrix BaseType) | ||
141 | quickCheck $ \m -> m |=| asFortran (m :: Matrix BaseType) | ||
142 | quickCheck $ \m -> m |=| (asC . asFortran) (m :: Matrix BaseType) | ||
143 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| mulF m1 (m2 :: Matrix BaseType) | ||
144 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| trans (mulF (trans m2) (trans m1 :: Matrix BaseType)) | ||
145 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| multiplyG m1 (m2 :: Matrix BaseType) | ||
146 | quickCheck (svdTestR svdR mulC) | ||
147 | quickCheck (svdTestR svdR mulF) | ||
148 | quickCheck (svdTestC svdC mulC) | ||
149 | quickCheck (svdTestC svdC mulF) | ||
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index 4383e79..0664e9b 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs | |||
@@ -179,6 +179,15 @@ joinVert ms = case common cols ms of | |||
179 | joinHoriz :: Field t => [Matrix t] -> Matrix t | 179 | joinHoriz :: Field t => [Matrix t] -> Matrix t |
180 | joinHoriz ms = trans. joinVert . map trans $ ms | 180 | joinHoriz ms = trans. joinVert . map trans $ ms |
181 | 181 | ||
182 | -- | creates a complex vector from vectors with real and imaginary parts | ||
183 | toComplex :: (Vector Double, Vector Double) -> Vector (Complex Double) | ||
184 | toComplex (r,i) = asComplex $ cdat $ fromColumns [r,i] | ||
185 | |||
186 | -- | obtains the complex conjugate of a complex vector | ||
187 | conj :: Vector (Complex Double) -> Vector (Complex Double) | ||
188 | conj v = asComplex $ cdat $ reshape 2 (asReal v) `mulC` diag (fromList [1,-1]) | ||
189 | where mulC = multiply RowMajor | ||
190 | |||
182 | ------------------------------------------------------------------------------ | 191 | ------------------------------------------------------------------------------ |
183 | 192 | ||
184 | -- | Reverse rows | 193 | -- | Reverse rows |
@@ -191,7 +200,7 @@ fliprl m = fromColumns . reverse . toColumns $ m | |||
191 | 200 | ||
192 | ----------------------------------------------------------------- | 201 | ----------------------------------------------------------------- |
193 | 202 | ||
194 | liftMatrix f m = m { dat = f dat, tdat = f tdat } -- check sizes | 203 | liftMatrix f m = m { dat = f (dat m), tdat = f (tdat m) } -- check sizes |
195 | 204 | ||
196 | ------------------------------------------------------------------ | 205 | ------------------------------------------------------------------ |
197 | 206 | ||
@@ -291,7 +300,7 @@ subMatrixG (r0,c0) (rt,ct) x = reshape ct $ fromList $ concat $ map (subList c0 | |||
291 | diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do | 300 | diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do |
292 | m <- createMatrix RowMajor n n | 301 | m <- createMatrix RowMajor n n |
293 | fun // vec v // mat dat m // check msg [dat m] | 302 | fun // vec v // mat dat m // check msg [dat m] |
294 | return m | 303 | return m {tdat = dat m} |
295 | 304 | ||
296 | -- | diagonal matrix from a real vector | 305 | -- | diagonal matrix from a real vector |
297 | diagR :: Vector Double -> Matrix Double | 306 | diagR :: Vector Double -> Matrix Double |
@@ -319,6 +328,6 @@ diagG v = reshape c $ fromList $ [ l!!(i-1) * delta k i | k <- [1..c], i <- [1.. | |||
319 | diagRect s r c | 328 | diagRect s r c |
320 | | dim s < min r c = error "diagRect" | 329 | | dim s < min r c = error "diagRect" |
321 | | r == c = diag s | 330 | | r == c = diag s |
322 | | r < c = joinHoriz [diag s , zeros (r,c-r)] | 331 | | r < c = trans $ diagRect s c r |
323 | | otherwise = joinVert [diag s , zeros (r-c,c)] | 332 | | r > c = joinVert [diag s , zeros (r-c,c)] |
324 | where zeros (r,c) = reshape c $ constant (r*c) 0 | 333 | where zeros (r,c) = reshape c $ constant (r*c) 0 |
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs index 36d5df7..8f4e6a4 100644 --- a/lib/Data/Packed/Internal/Vector.hs +++ b/lib/Data/Packed/Internal/Vector.hs | |||
@@ -145,6 +145,7 @@ asReal v = V { dim = 2*dim v, fptr = castForeignPtr (fptr v), ptr = castPtr (pt | |||
145 | asComplex :: Vector Double -> Vector (Complex Double) | 145 | asComplex :: Vector Double -> Vector (Complex Double) |
146 | asComplex v = V { dim = dim v `div` 2, fptr = castForeignPtr (fptr v), ptr = castPtr (ptr v) } | 146 | asComplex v = V { dim = dim v `div` 2, fptr = castForeignPtr (fptr v), ptr = castPtr (ptr v) } |
147 | 147 | ||
148 | |||
148 | constantG n x = fromList (replicate n x) | 149 | constantG n x = fromList (replicate n x) |
149 | 150 | ||
150 | constantR :: Int -> Double -> Vector Double | 151 | constantR :: Int -> Double -> Vector Double |