diff options
author | Alberto Ruiz <aruiz@um.es> | 2007-09-21 18:19:13 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2007-09-21 18:19:13 +0000 |
commit | d4cb2692f9dae748da23371057a983deca4b2f80 (patch) | |
tree | 3d3b3ef373bb0b6ed2a86c30accf85139397b819 /examples/tests.hs | |
parent | edfaf9e0d1dcfccc9015476510a23e8cf64288be (diff) |
improved tests
Diffstat (limited to 'examples/tests.hs')
-rw-r--r-- | examples/tests.hs | 345 |
1 files changed, 345 insertions, 0 deletions
diff --git a/examples/tests.hs b/examples/tests.hs new file mode 100644 index 0000000..ec838ca --- /dev/null +++ b/examples/tests.hs | |||
@@ -0,0 +1,345 @@ | |||
1 | {-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-} | ||
2 | |||
3 | module Main where | ||
4 | |||
5 | import Data.Packed.Internal((>|<), fdat, cdat, multiply', multiplyG, MatrixOrder(..),debug) | ||
6 | import GSL hiding (sin,cos,exp,choose) | ||
7 | import LinearAlgebra | ||
8 | import Test.QuickCheck hiding (test) | ||
9 | import Test.HUnit hiding ((~:),test) | ||
10 | import System.Random(randomRs,mkStdGen) | ||
11 | |||
12 | type RM = Matrix Double | ||
13 | type CM = Matrix (Complex Double) | ||
14 | |||
15 | dist :: (Normed t, Num t) => t -> t -> Double | ||
16 | dist a b = pnorm Infinity (a-b) | ||
17 | |||
18 | infixl 4 |~| | ||
19 | a |~| b = a :~8~: b | ||
20 | |||
21 | data Aprox a = (:~) a Int | ||
22 | |||
23 | (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool | ||
24 | a :~n~: b = dist a b < 10^^(-n) | ||
25 | |||
26 | |||
27 | maxdim = 10 | ||
28 | |||
29 | instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where | ||
30 | arbitrary = do | ||
31 | r <- arbitrary | ||
32 | i <- arbitrary | ||
33 | return (r:+i) | ||
34 | coarbitrary = undefined | ||
35 | |||
36 | instance (Field a, Arbitrary a) => Arbitrary (Matrix a) where | ||
37 | arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) | ||
38 | m <- choose (1,maxdim) | ||
39 | n <- choose (1,maxdim) | ||
40 | l <- vector (m*n) | ||
41 | ctype <- arbitrary | ||
42 | let h = if ctype then (m><n) else (m>|<n) | ||
43 | trMode <- arbitrary | ||
44 | let tr = if trMode then trans else id | ||
45 | return $ tr (h l) | ||
46 | coarbitrary = undefined | ||
47 | |||
48 | data PairM a = PairM (Matrix a) (Matrix a) deriving Show | ||
49 | instance (Num a, Field a, Arbitrary a) => Arbitrary (PairM a) where | ||
50 | arbitrary = do | ||
51 | a <- choose (1,maxdim) | ||
52 | b <- choose (1,maxdim) | ||
53 | c <- choose (1,maxdim) | ||
54 | l1 <- vector (a*b) | ||
55 | l2 <- vector (b*c) | ||
56 | return $ PairM ((a><b) (map fromIntegral (l1::[Int]))) ((b><c) (map fromIntegral (l2::[Int]))) | ||
57 | --return $ PairM ((a><b) l1) ((b><c) l2) | ||
58 | coarbitrary = undefined | ||
59 | |||
60 | data SqM a = SqM (Matrix a) deriving Show | ||
61 | sqm (SqM a) = a | ||
62 | instance (Field a, Arbitrary a) => Arbitrary (SqM a) where | ||
63 | arbitrary = do | ||
64 | n <- choose (1,maxdim) | ||
65 | l <- vector (n*n) | ||
66 | return $ SqM $ (n><n) l | ||
67 | coarbitrary = undefined | ||
68 | |||
69 | data Sym a = Sym (Matrix a) deriving Show | ||
70 | sym (Sym a) = a | ||
71 | instance (Linear Vector a, Arbitrary a) => Arbitrary (Sym a) where | ||
72 | arbitrary = do | ||
73 | SqM m <- arbitrary | ||
74 | return $ Sym (m + trans m) | ||
75 | coarbitrary = undefined | ||
76 | |||
77 | data Her = Her (Matrix (Complex Double)) deriving Show | ||
78 | her (Her a) = a | ||
79 | instance {-(Field a, Arbitrary a, Num a) =>-} Arbitrary Her where | ||
80 | arbitrary = do | ||
81 | SqM m <- arbitrary | ||
82 | return $ Her (m + conjTrans m) | ||
83 | coarbitrary = undefined | ||
84 | |||
85 | data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show | ||
86 | instance (Num a, Field a, Arbitrary a) => Arbitrary (PairSM a) where | ||
87 | arbitrary = do | ||
88 | a <- choose (1,maxdim) | ||
89 | c <- choose (1,maxdim) | ||
90 | l1 <- vector (a*a) | ||
91 | l2 <- vector (a*c) | ||
92 | return $ PairSM ((a><a) (map fromIntegral (l1::[Int]))) ((a><c) (map fromIntegral (l2::[Int]))) | ||
93 | --return $ PairSM ((a><a) l1) ((a><c) l2) | ||
94 | coarbitrary = undefined | ||
95 | |||
96 | instance (Field a, Arbitrary a) => Arbitrary (Vector a) where | ||
97 | arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) | ||
98 | m <- choose (1,maxdim^2) | ||
99 | l <- vector m | ||
100 | return $ fromList l | ||
101 | coarbitrary = undefined | ||
102 | |||
103 | data PairV a = PairV (Vector a) (Vector a) | ||
104 | instance (Field a, Arbitrary a) => Arbitrary (PairV a) where | ||
105 | arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) | ||
106 | m <- choose (1,maxdim^2) | ||
107 | l1 <- vector m | ||
108 | l2 <- vector m | ||
109 | return $ PairV (fromList l1) (fromList l2) | ||
110 | coarbitrary = undefined | ||
111 | |||
112 | ---------------------------------------------------------------------- | ||
113 | |||
114 | test str b = TestCase $ assertBool str b | ||
115 | |||
116 | ---------------------------------------------------------------------- | ||
117 | |||
118 | pseudorandomR seed (n,m) = reshape m $ fromList $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed | ||
119 | |||
120 | pseudorandomC seed (n,m) = toComplex (pseudorandomR seed (n,m), pseudorandomR (seed+1) (n,m)) | ||
121 | |||
122 | bigmat = m + trans m :: RM | ||
123 | where m = pseudorandomR 18 (1000,1000) | ||
124 | bigmatc = mc + conjTrans mc ::CM | ||
125 | where mc = pseudorandomC 19 (1000,1000) | ||
126 | |||
127 | ---------------------------------------------------------------------- | ||
128 | |||
129 | |||
130 | m = (3><3) | ||
131 | [ 1, 2, 3 | ||
132 | , 4, 5, 7 | ||
133 | , 2, 8, 4 :: Double | ||
134 | ] | ||
135 | |||
136 | mc = (3><3) | ||
137 | [ 1, 2, 3 | ||
138 | , 4, 5, 7 | ||
139 | , 2, 8, i | ||
140 | ] | ||
141 | |||
142 | |||
143 | mr = (3><4) | ||
144 | [ 1, 2, 3, 4, | ||
145 | 2, 4, 6, 8, | ||
146 | 1, 1, 1, 2:: Double | ||
147 | ] | ||
148 | |||
149 | mrc = (3><4) | ||
150 | [ 1, 2, 3, 4, | ||
151 | 2, 4, 6, 8, | ||
152 | i, i, i, 2 | ||
153 | ] | ||
154 | |||
155 | a = (3><4) | ||
156 | [ 1, 0, 0, 0 | ||
157 | , 0, 2, 0, 0 | ||
158 | , 0, 0, 0, 0 :: Double | ||
159 | ] | ||
160 | |||
161 | b = (3><4) | ||
162 | [ 1, 0, 0, 0 | ||
163 | , 0, 2, 3, 0 | ||
164 | , 0, 0, 4, 0 :: Double | ||
165 | ] | ||
166 | |||
167 | ac = (2><3) [1 .. 6::Double] | ||
168 | bc = (3><4) [7 .. 18::Double] | ||
169 | |||
170 | af = (2>|<3) [1,4,2,5,3,6::Double] | ||
171 | bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double] | ||
172 | |||
173 | ------------------------------------------------------- | ||
174 | |||
175 | detTest = det m == 26 && det mc == 38 :+ (-3) | ||
176 | |||
177 | invTest m = degenerate m || m <> inv m |~| ident (rows m) | ||
178 | |||
179 | pinvTest m = m <> p <> m |~| m | ||
180 | && p <> m <> p |~| p | ||
181 | && hermitian (m<>p) | ||
182 | && hermitian (p<>m) | ||
183 | where p = pinv m | ||
184 | |||
185 | square m = rows m == cols m | ||
186 | |||
187 | orthonormal m = square m && m <> ctrans m |~| ident (rows m) | ||
188 | |||
189 | hermitian m = m |~| ctrans m | ||
190 | |||
191 | svdTest svd m = u <> real d <> trans v |~| m | ||
192 | && orthonormal u && orthonormal v | ||
193 | where (u,d,v) = full svd m | ||
194 | |||
195 | svdTest' svd m = m |~| 0 || u <> real (diag s) <> trans v |~| m | ||
196 | where (u,s,v) = economy svd m | ||
197 | |||
198 | eigTest m = complex m <> v |~| v <> diag s | ||
199 | where (s, v) = eig m | ||
200 | |||
201 | eigTestSH m = m <> v |~| v <> real (diag s) | ||
202 | && orthonormal v | ||
203 | where (s, v) = eigSH m | ||
204 | |||
205 | rank m | m |~| 0 = 0 | ||
206 | | otherwise = dim s where (_,s,_) = economy svd m | ||
207 | |||
208 | zeros (r,c) = reshape c (constant 0 (r*c)) | ||
209 | |||
210 | ones (r,c) = zeros (r,c) + 1 | ||
211 | |||
212 | degenerate m = rank m < min (rows m) (cols m) | ||
213 | |||
214 | prec = 1E-15 | ||
215 | |||
216 | singular m = s1 < prec || s2/s1 < prec | ||
217 | where (_,ss,_) = svd m | ||
218 | s = toList ss | ||
219 | s1 = maximum s | ||
220 | s2 = minimum s | ||
221 | |||
222 | nullspaceTest m = null nl || m <> n |~| zeros (r,c) -- 0 | ||
223 | where nl = nullspacePrec 1 m | ||
224 | n = fromColumns nl | ||
225 | r = rows m | ||
226 | c = cols m - rank m | ||
227 | |||
228 | -------------------------------------------------------------------- | ||
229 | |||
230 | polyEval cs x = foldr (\c ac->ac*x+c) 0 cs | ||
231 | |||
232 | polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p)) | ||
233 | |||
234 | |||
235 | polySolveTest = test "polySolve" (polySolveTest' [1,2,3,4]) | ||
236 | |||
237 | --------------------------------------------------------------------- | ||
238 | |||
239 | quad f a b = fst $ integrateQAGS 1E-9 100 f a b | ||
240 | |||
241 | -- A multiple integral can be easily defined using partial application | ||
242 | quad2 f a b g1 g2 = quad h a b | ||
243 | where h x = quad (f x) (g1 x) (g2 x) | ||
244 | |||
245 | volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) | ||
246 | 0 r (const 0) (\x->sqrt (r*r-x*x)) | ||
247 | |||
248 | epsTol = 1E-8::Double | ||
249 | |||
250 | integrateTest = test "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < epsTol) | ||
251 | |||
252 | --------------------------------------------------------------------- | ||
253 | |||
254 | besselTest = test "bessel_J0_e" ( abs (r-expected) < e ) | ||
255 | where (r,e) = bessel_J0_e 5.0 | ||
256 | expected = -0.17759677131433830434739701 | ||
257 | |||
258 | exponentialTest = test "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 ) | ||
259 | where (v,e,err) = exp_e10_e 30.0 | ||
260 | expected = exp 30.0 | ||
261 | |||
262 | gammaTest = test "gamma" (gamma 5 == 24.0) | ||
263 | |||
264 | --------------------------------------------------------------------- | ||
265 | |||
266 | asFortran m = (rows m >|< cols m) $ toList (fdat m) | ||
267 | asC m = (rows m >< cols m) $ toList (cdat m) | ||
268 | |||
269 | mulC a b = multiply' RowMajor a b | ||
270 | mulF a b = multiply' ColumnMajor a b | ||
271 | |||
272 | --------------------------------------------------------------------- | ||
273 | |||
274 | tests = do | ||
275 | putStrLn "--------- internal -----" | ||
276 | quickCheck ((\m -> m == trans m).sym :: Sym Double -> Bool) | ||
277 | quickCheck ((\m -> m == trans m).sym :: Sym (Complex Double) -> Bool) | ||
278 | quickCheck $ \l -> null l || (toList . fromList) l == (l :: [Double]) | ||
279 | quickCheck $ \l -> null l || (toList . fromList) l == (l :: [Complex Double]) | ||
280 | quickCheck $ \m -> m == asC (m :: RM) | ||
281 | quickCheck $ \m -> m == asC (m :: CM) | ||
282 | quickCheck $ \m -> m == asFortran (m :: RM) | ||
283 | quickCheck $ \m -> m == asFortran (m :: CM) | ||
284 | quickCheck $ \m -> m == (asC . asFortran) (m :: RM) | ||
285 | quickCheck $ \m -> m == (asC . asFortran) (m :: CM) | ||
286 | putStrLn "--------- multiply ----" | ||
287 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: RM) | ||
288 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: CM) | ||
289 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: RM)) | ||
290 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: CM)) | ||
291 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: RM) | ||
292 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: CM) | ||
293 | putStrLn "--------- svd ---------" | ||
294 | quickCheck (svdTest svdR) | ||
295 | quickCheck (svdTest svdRdd) | ||
296 | quickCheck (svdTest svdC) | ||
297 | quickCheck (svdTest' svdR) | ||
298 | quickCheck (svdTest' svdRdd) | ||
299 | quickCheck (svdTest' svdC) | ||
300 | quickCheck (svdTest' svdg) | ||
301 | putStrLn "--------- eig ---------" | ||
302 | quickCheck (eigTest . sqm :: SqM Double -> Bool) | ||
303 | quickCheck (eigTest . sqm :: SqM (Complex Double) -> Bool) | ||
304 | quickCheck (eigTestSH . sym :: Sym Double -> Bool) | ||
305 | quickCheck (eigTestSH . her :: Her -> Bool) | ||
306 | putStrLn "--------- inv ------" | ||
307 | quickCheck (invTest . sqm :: SqM Double -> Bool) | ||
308 | quickCheck (invTest . sqm :: SqM (Complex Double) -> Bool) | ||
309 | putStrLn "--------- pinv ------" | ||
310 | quickCheck (pinvTest . sqm :: SqM Double -> Bool) | ||
311 | quickCheck (pinvTest . sqm :: SqM (Complex Double) -> Bool) | ||
312 | putStrLn "--------- nullspace ------" | ||
313 | quickCheck (nullspaceTest :: RM -> Bool) | ||
314 | quickCheck (nullspaceTest :: CM -> Bool) | ||
315 | putStrLn "--------- vector operations ------" | ||
316 | quickCheck $ (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM)) | ||
317 | quickCheck $ (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM)) | ||
318 | quickCheck $ (\u -> cos u * tan u |~| sin (u::RM)) | ||
319 | quickCheck $ (\u -> (cos u * tan u) :~6~: sin (u::CM)) | ||
320 | runTestTT $ TestList | ||
321 | [ test "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM) | ||
322 | , test "arith2" $ (((1+i) .* ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( (140*i-51).*1 :: CM) | ||
323 | , test "arith3" $ exp (i.*ones(10,10)*pi) + 1 |~| 0 | ||
324 | ] | ||
325 | putStrLn "--------- GSL ------" | ||
326 | quickCheck $ \v -> ifft (fft v) |~| v | ||
327 | runTestTT $ TestList | ||
328 | [ gammaTest | ||
329 | , besselTest | ||
330 | , exponentialTest | ||
331 | , integrateTest | ||
332 | , polySolveTest | ||
333 | , test "det" detTest | ||
334 | ] | ||
335 | |||
336 | bigtests = do | ||
337 | putStrLn "--------- big matrices -----" | ||
338 | runTestTT $ TestList | ||
339 | [ test "eigS" $ eigTestSH bigmat | ||
340 | , test "eigH" $ eigTestSH bigmatc | ||
341 | , test "eigR" $ eigTest bigmat | ||
342 | , test "eigC" $ eigTest bigmatc | ||
343 | ] | ||
344 | |||
345 | main = tests | ||