diff options
author | Alberto Ruiz <aruiz@um.es> | 2008-02-27 10:57:29 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2008-02-27 10:57:29 +0000 |
commit | 2c99bcb1de3fae6e2e075316126cc70658e20ac9 (patch) | |
tree | c80dfb5933b88efb3cd3d913240a690fc59ae0c3 | |
parent | 500f5fca244dadab494655aa73d7183df1c87c50 (diff) |
tests reorganized
-rw-r--r-- | examples/tests.hs | 490 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests.hs | 144 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests/Properties.hs | 34 |
3 files changed, 177 insertions, 491 deletions
diff --git a/examples/tests.hs b/examples/tests.hs index b6c9a36..cd923cd 100644 --- a/examples/tests.hs +++ b/examples/tests.hs | |||
@@ -1,138 +1,13 @@ | |||
1 | {-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-} | ||
2 | |||
3 | module Main where | 1 | module Main where |
4 | 2 | ||
5 | import Numeric.GSL hiding (sin,cos,exp,choose) | ||
6 | import Numeric.LinearAlgebra | 3 | import Numeric.LinearAlgebra |
7 | import Numeric.LinearAlgebra.LAPACK | 4 | import Numeric.LinearAlgebra.Tests |
8 | import qualified Numeric.GSL.Matrix as GSL | ||
9 | import Test.QuickCheck hiding (test) | ||
10 | import Test.HUnit hiding ((~:),test) | ||
11 | import System.Random(randomRs,mkStdGen) | 5 | import System.Random(randomRs,mkStdGen) |
12 | import System.Info | 6 | import Test.HUnit hiding (test) |
13 | import Data.List(foldl1', transpose) | ||
14 | import System(getArgs) | 7 | import System(getArgs) |
15 | import Debug.Trace(trace) | ||
16 | |||
17 | debug x = trace (show x) x | ||
18 | |||
19 | type RM = Matrix Double | ||
20 | type CM = Matrix (Complex Double) | ||
21 | |||
22 | -- relative error | ||
23 | dist :: (Normed t, Num t) => t -> t -> Double | ||
24 | dist a b = r | ||
25 | where norm = pnorm Infinity | ||
26 | na = norm a | ||
27 | nb = norm b | ||
28 | nab = norm (a-b) | ||
29 | mx = max na nb | ||
30 | mn = min na nb | ||
31 | r = if mn < eps | ||
32 | then mx | ||
33 | else nab/mx | ||
34 | |||
35 | infixl 4 |~| | ||
36 | a |~| b = a :~10~: b | ||
37 | |||
38 | data Aprox a = (:~) a Int | ||
39 | |||
40 | (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool | ||
41 | a :~n~: b = dist a b < 10^^(-n) | ||
42 | |||
43 | |||
44 | maxdim = 10 | ||
45 | |||
46 | instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where | ||
47 | arbitrary = do | ||
48 | r <- arbitrary | ||
49 | i <- arbitrary | ||
50 | return (r:+i) | ||
51 | coarbitrary = undefined | ||
52 | |||
53 | instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where | ||
54 | arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) | ||
55 | m <- choose (1,maxdim) | ||
56 | n <- choose (1,maxdim) | ||
57 | l <- vector (m*n) | ||
58 | ctype <- arbitrary | ||
59 | let h = if ctype then (m><n) else (m>|<n) | ||
60 | trMode <- arbitrary | ||
61 | let tr = if trMode then trans else id | ||
62 | return $ tr (h l) | ||
63 | coarbitrary = undefined | ||
64 | |||
65 | data PairM a = PairM (Matrix a) (Matrix a) deriving Show | ||
66 | instance (Num a, Element a, Arbitrary a) => Arbitrary (PairM a) where | ||
67 | arbitrary = do | ||
68 | a <- choose (1,maxdim) | ||
69 | b <- choose (1,maxdim) | ||
70 | c <- choose (1,maxdim) | ||
71 | l1 <- vector (a*b) | ||
72 | l2 <- vector (b*c) | ||
73 | return $ PairM ((a><b) (map fromIntegral (l1::[Int]))) ((b><c) (map fromIntegral (l2::[Int]))) | ||
74 | --return $ PairM ((a><b) l1) ((b><c) l2) | ||
75 | coarbitrary = undefined | ||
76 | |||
77 | data SqM a = SqM (Matrix a) deriving Show | ||
78 | sqm (SqM a) = a | ||
79 | instance (Element a, Arbitrary a) => Arbitrary (SqM a) where | ||
80 | arbitrary = do | ||
81 | n <- choose (1,maxdim) | ||
82 | l <- vector (n*n) | ||
83 | return $ SqM $ (n><n) l | ||
84 | coarbitrary = undefined | ||
85 | |||
86 | data Sym a = Sym (Matrix a) deriving Show | ||
87 | sym (Sym a) = a | ||
88 | instance (Linear Vector a, Arbitrary a) => Arbitrary (Sym a) where | ||
89 | arbitrary = do | ||
90 | SqM m <- arbitrary | ||
91 | return $ Sym (m + trans m) | ||
92 | coarbitrary = undefined | ||
93 | 8 | ||
94 | data Her = Her (Matrix (Complex Double)) deriving Show | ||
95 | her (Her a) = a | ||
96 | instance {-(Field a, Arbitrary a, Num a) =>-} Arbitrary Her where | ||
97 | arbitrary = do | ||
98 | SqM m <- arbitrary | ||
99 | return $ Her (m + ctrans m) | ||
100 | coarbitrary = undefined | ||
101 | 9 | ||
102 | data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show | 10 | pseudorandomR seed (n,m) = reshape m $ fromList $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed |
103 | instance (Num a, Field a, Arbitrary a) => Arbitrary (PairSM a) where | ||
104 | arbitrary = do | ||
105 | a <- choose (1,maxdim) | ||
106 | c <- choose (1,maxdim) | ||
107 | l1 <- vector (a*a) | ||
108 | l2 <- vector (a*c) | ||
109 | return $ PairSM ((a><a) (map fromIntegral (l1::[Int]))) ((a><c) (map fromIntegral (l2::[Int]))) | ||
110 | --return $ PairSM ((a><a) l1) ((a><c) l2) | ||
111 | coarbitrary = undefined | ||
112 | |||
113 | instance (Field a, Arbitrary a) => Arbitrary (Vector a) where | ||
114 | arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) | ||
115 | m <- choose (1,maxdim^2) | ||
116 | l <- vector m | ||
117 | return $ fromList l | ||
118 | coarbitrary = undefined | ||
119 | |||
120 | data PairV a = PairV (Vector a) (Vector a) | ||
121 | instance (Field a, Arbitrary a) => Arbitrary (PairV a) where | ||
122 | arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) | ||
123 | m <- choose (1,maxdim^2) | ||
124 | l1 <- vector m | ||
125 | l2 <- vector m | ||
126 | return $ PairV (fromList l1) (fromList l2) | ||
127 | coarbitrary = undefined | ||
128 | |||
129 | ---------------------------------------------------------------------- | ||
130 | |||
131 | test str b = TestCase $ assertBool str b | ||
132 | |||
133 | ---------------------------------------------------------------------- | ||
134 | |||
135 | pseudorandomR seed (n,m) = reshape m $ fromList $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed | ||
136 | 11 | ||
137 | pseudorandomC seed (n,m) = toComplex (pseudorandomR seed (n,m), pseudorandomR (seed+1) (n,m)) | 12 | pseudorandomC seed (n,m) = toComplex (pseudorandomR seed (n,m), pseudorandomR (seed+1) (n,m)) |
138 | 13 | ||
@@ -141,366 +16,23 @@ bigmat = m + trans m :: RM | |||
141 | bigmatc = mc + ctrans mc ::CM | 16 | bigmatc = mc + ctrans mc ::CM |
142 | where mc = pseudorandomC 19 (1000,1000) | 17 | where mc = pseudorandomC 19 (1000,1000) |
143 | 18 | ||
144 | ---------------------------------------------------------------------- | 19 | utest str b = TestCase $ assertBool str b |
145 | |||
146 | |||
147 | m = (3><3) | ||
148 | [ 1, 2, 3 | ||
149 | , 4, 5, 7 | ||
150 | , 2, 8, 4 :: Double | ||
151 | ] | ||
152 | |||
153 | mc = (3><3) | ||
154 | [ 1, 2, 3 | ||
155 | , 4, 5, 7 | ||
156 | , 2, 8, i | ||
157 | ] | ||
158 | |||
159 | |||
160 | mr = (3><4) | ||
161 | [ 1, 2, 3, 4, | ||
162 | 2, 4, 6, 8, | ||
163 | 1, 1, 1, 2:: Double | ||
164 | ] | ||
165 | |||
166 | mrc = (3><4) | ||
167 | [ 1, 2, 3, 4, | ||
168 | 2, 4, 6, 8, | ||
169 | i, i, i, 2 | ||
170 | ] | ||
171 | |||
172 | a = (3><4) | ||
173 | [ 1, 0, 0, 0 | ||
174 | , 0, 2, 0, 0 | ||
175 | , 0, 0, 0, 0 :: Double | ||
176 | ] | ||
177 | |||
178 | b = (3><4) | ||
179 | [ 1, 0, 0, 0 | ||
180 | , 0, 2, 3, 0 | ||
181 | , 0, 0, 4, 0 :: Double | ||
182 | ] | ||
183 | |||
184 | ac = (2><3) [1 .. 6::Double] | ||
185 | bc = (3><4) [7 .. 18::Double] | ||
186 | |||
187 | af = (2>|<3) [1,4,2,5,3,6::Double] | ||
188 | bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double] | ||
189 | |||
190 | ------------------------------------------------------- | ||
191 | 20 | ||
192 | feye n = flipud (ident n) :: Matrix Double | 21 | feye n = flipud (ident n) :: Matrix Double |
193 | 22 | ||
194 | |||
195 | luTest1 m = m |~| p <> l <> u | ||
196 | where (l,u,p,_) = lu m | ||
197 | |||
198 | detTest1 = det m == 26 | ||
199 | && det mc == 38 :+ (-3) | ||
200 | && det (feye 2) == -1 | ||
201 | |||
202 | detTest2 m = s d1 |~| s d2 | ||
203 | where d1 = det m | ||
204 | d2 = det' m * det q | ||
205 | det' m = product $ toList $ takeDiag r | ||
206 | (q,r) = qr m | ||
207 | s x = fromList [x] | ||
208 | |||
209 | invTest m = degenerate m || m <> inv m |~| ident (rows m) | ||
210 | |||
211 | pinvTest m = m <> p <> m |~| m | ||
212 | && p <> m <> p |~| p | ||
213 | && hermitian (m<>p) | ||
214 | && hermitian (p<>m) | ||
215 | where p = pinv m | ||
216 | |||
217 | square m = rows m == cols m | ||
218 | |||
219 | unitary m = square m && m <> ctrans m |~| ident (rows m) | ||
220 | |||
221 | hermitian m = m |~| ctrans m | ||
222 | |||
223 | upperTriang m = rows m == 1 || down == z | ||
224 | where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m)) | ||
225 | z = constant 0 (dim down) | ||
226 | |||
227 | upperHessenberg m = rows m < 3 || down == z | ||
228 | where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m)) | ||
229 | z = constant 0 (dim down) | ||
230 | |||
231 | svdTest svd m = u <> real d <> trans v |~| m | ||
232 | && unitary u && unitary v | ||
233 | where (u,d,v) = full svd m | ||
234 | |||
235 | svdTest' svd m = m |~| 0 || u <> real (diag s) <> trans v |~| m | ||
236 | where (u,s,v) = economy svd m | ||
237 | |||
238 | eigTest m = complex m <> v |~| v <> diag s | ||
239 | where (s, v) = eig m | ||
240 | |||
241 | eigTestSH m = m <> v |~| v <> real (diag s) | ||
242 | && unitary v | ||
243 | && m |~| v <> real (diag s) <> ctrans v | ||
244 | where (s, v) = eigSH m | ||
245 | |||
246 | zeros (r,c) = reshape c (constant 0 (r*c)) | ||
247 | |||
248 | ones (r,c) = zeros (r,c) + 1 | ||
249 | |||
250 | degenerate m = rank m < min (rows m) (cols m) | ||
251 | |||
252 | prec = 1E-15 | ||
253 | |||
254 | singular m = s1 < prec || s2/s1 < prec | ||
255 | where (_,ss,_) = svd m | ||
256 | s = toList ss | ||
257 | s1 = maximum s | ||
258 | s2 = minimum s | ||
259 | |||
260 | nullspaceTest m = null nl || m <> n |~| zeros (r,c) -- 0 | ||
261 | where nl = nullspacePrec 1 m | ||
262 | n = fromColumns nl | ||
263 | r = rows m | ||
264 | c = cols m - rank m | ||
265 | |||
266 | -------------------------------------------------------------------- | ||
267 | |||
268 | polyEval cs x = foldr (\c ac->ac*x+c) 0 cs | ||
269 | |||
270 | polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p)) | ||
271 | |||
272 | |||
273 | polySolveTest = test "polySolve" (polySolveTest' [1,2,3,4]) | ||
274 | |||
275 | --------------------------------------------------------------------- | ||
276 | |||
277 | quad f a b = fst $ integrateQAGS 1E-9 100 f a b | ||
278 | |||
279 | -- A multiple integral can be easily defined using partial application | ||
280 | quad2 f a b g1 g2 = quad h a b | ||
281 | where h x = quad (f x) (g1 x) (g2 x) | ||
282 | |||
283 | volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) | ||
284 | 0 r (const 0) (\x->sqrt (r*r-x*x)) | ||
285 | |||
286 | epsTol = 1E-8::Double | ||
287 | |||
288 | integrateTest = test "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < epsTol) | ||
289 | |||
290 | --------------------------------------------------------------------- | ||
291 | |||
292 | besselTest = test "bessel_J0_e" ( abs (r-expected) < e ) | ||
293 | where (r,e) = bessel_J0_e 5.0 | ||
294 | expected = -0.17759677131433830434739701 | ||
295 | |||
296 | exponentialTest = test "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 ) | ||
297 | where (v,e,err) = exp_e10_e 30.0 | ||
298 | expected = exp 30.0 | ||
299 | |||
300 | gammaTest = test "gamma" (gamma 5 == 24.0) | ||
301 | |||
302 | --------------------------------------------------------------------- | ||
303 | |||
304 | cholRTest = chol ((2><2) [1,2,2,9::Double]) == (2><2) [1,2,0,2.23606797749979] | ||
305 | cholCTest = chol ((2><2) [1,2,2,9::Complex Double]) == (2><2) [1,2,0,2.23606797749979] | ||
306 | |||
307 | --------------------------------------------------------------------- | ||
308 | |||
309 | qrTest qr m = q <> r |~| m && unitary q && upperTriang r | ||
310 | where (q,r) = qr m | ||
311 | |||
312 | --------------------------------------------------------------------- | ||
313 | |||
314 | hessTest m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h | ||
315 | where (p,h) = hess m | ||
316 | |||
317 | --------------------------------------------------------------------- | ||
318 | |||
319 | schurTest1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s | ||
320 | where (u,s) = schur m | ||
321 | |||
322 | schurTest2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme | ||
323 | where (u,s) = schur m | ||
324 | |||
325 | --------------------------------------------------------------------- | ||
326 | |||
327 | nd1 = (3><3) [ 1/2, 1/4, 1/4 | ||
328 | , 0/1, 1/2, 1/4 | ||
329 | , 1/2, 1/4, 1/2 :: Double] | ||
330 | |||
331 | nd2 = (2><2) [1, 0, 1, 1:: Complex Double] | ||
332 | |||
333 | expmTest1 = expm nd1 :~14~: (3><3) | ||
334 | [ 1.762110887278176 | ||
335 | , 0.478085470590435 | ||
336 | , 0.478085470590435 | ||
337 | , 0.104719410945666 | ||
338 | , 1.709751181805343 | ||
339 | , 0.425725765117601 | ||
340 | , 0.851451530235203 | ||
341 | , 0.530445176063267 | ||
342 | , 1.814470592751009 ] | ||
343 | |||
344 | expmTest2 = expm nd2 :~15~: (2><2) | ||
345 | [ 2.718281828459045 | ||
346 | , 0.000000000000000 | ||
347 | , 2.718281828459045 | ||
348 | , 2.718281828459045 ] | ||
349 | |||
350 | expmTestDiag m = expm (logm m) |~| complex m | ||
351 | where logm m = matFunc Prelude.log m | ||
352 | |||
353 | |||
354 | |||
355 | --------------------------------------------------------------------- | ||
356 | |||
357 | asFortran m = (rows m >|< cols m) $ toList (flatten $ trans m) | ||
358 | asC m = (rows m >< cols m) $ toList (flatten m) | ||
359 | |||
360 | mulC a b = a <> b | ||
361 | mulF a b = trans $ trans b <> trans a | ||
362 | |||
363 | ------------------------------------------------------------------------- | ||
364 | |||
365 | multiplyG a b = reshape (cols b) $ fromList $ concat $ multiplyL (toLists a) (toLists b) | ||
366 | where multiplyL a b = [[dotL x y | y <- transpose b] | x <- a] | ||
367 | dotL a b = sum (zipWith (*) a b) | ||
368 | |||
369 | r >|< c = f where | ||
370 | f l | dim v == r*c = reshapeF r v | ||
371 | | otherwise = error "(>|<)" | ||
372 | where v = fromList l | ||
373 | reshapeF r = trans . reshape r | ||
374 | |||
375 | --------------------------------------------------------------------- | ||
376 | |||
377 | rot :: Double -> Matrix Double | ||
378 | rot a = (3><3) [ c,0,s | ||
379 | , 0,1,0 | ||
380 | ,-s,0,c ] | ||
381 | where c = cos a | ||
382 | s = sin a | ||
383 | |||
384 | fun n = foldl1' (<>) (map rot angles) | ||
385 | where angles = toList $ linspace n (0,1) | ||
386 | |||
387 | rotTest = fun (10^5) :~12~: rot 5E4 | ||
388 | |||
389 | --------------------------------------------------------------------- | ||
390 | |||
391 | tests = do | ||
392 | setErrorHandlerOff | ||
393 | putStrLn "--------- internal -----" | ||
394 | quickCheck ((\m -> m == trans m).sym :: Sym Double -> Bool) | ||
395 | quickCheck ((\m -> m == trans m).sym :: Sym (Complex Double) -> Bool) | ||
396 | quickCheck $ \l -> null l || (toList . fromList) l == (l :: [Double]) | ||
397 | quickCheck $ \l -> null l || (toList . fromList) l == (l :: [Complex Double]) | ||
398 | quickCheck $ \m -> m == asC (m :: RM) | ||
399 | quickCheck $ \m -> m == asC (m :: CM) | ||
400 | quickCheck $ \m -> m == asFortran (m :: RM) | ||
401 | quickCheck $ \m -> m == asFortran (m :: CM) | ||
402 | quickCheck $ \m -> m == (asC . asFortran) (m :: RM) | ||
403 | quickCheck $ \m -> m == (asC . asFortran) (m :: CM) | ||
404 | runTestTT $ TestList | ||
405 | [ test "1E5 rots" rotTest | ||
406 | ] | ||
407 | putStrLn "--------- multiply ----" | ||
408 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: RM) | ||
409 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: CM) | ||
410 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: RM)) | ||
411 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: CM)) | ||
412 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: RM) | ||
413 | quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: CM) | ||
414 | putStrLn "--------- lu ---------" | ||
415 | quickCheck (luTest1 :: RM->Bool) | ||
416 | quickCheck (luTest1 :: CM->Bool) | ||
417 | quickCheck (detTest2 . sqm :: SqM Double -> Bool) | ||
418 | quickCheck (detTest2 . sqm :: SqM (Complex Double) -> Bool) | ||
419 | runTestTT $ TestList | ||
420 | [ test "det1" detTest1 | ||
421 | ] | ||
422 | putStrLn "--------- svd ---------" | ||
423 | quickCheck (svdTest svdR) | ||
424 | quickCheck (svdTest svdRdd) | ||
425 | quickCheck (svdTest svdC) | ||
426 | quickCheck (svdTest' svdR) | ||
427 | quickCheck (svdTest' svdRdd) | ||
428 | quickCheck (svdTest' svdC) | ||
429 | quickCheck (svdTest' GSL.svdg) | ||
430 | putStrLn "--------- eig ---------" | ||
431 | quickCheck (eigTest . sqm :: SqM Double -> Bool) | ||
432 | quickCheck (eigTest . sqm :: SqM (Complex Double) -> Bool) | ||
433 | quickCheck (eigTestSH . sym :: Sym Double -> Bool) | ||
434 | quickCheck (eigTestSH . her :: Her -> Bool) | ||
435 | putStrLn "--------- inv ------" | ||
436 | quickCheck (invTest . sqm :: SqM Double -> Bool) | ||
437 | quickCheck (invTest . sqm :: SqM (Complex Double) -> Bool) | ||
438 | putStrLn "--------- pinv ------" | ||
439 | quickCheck (pinvTest ::RM->Bool) | ||
440 | if os == "mingw32" | ||
441 | then putStrLn "complex pinvTest skipped in this OS" | ||
442 | else quickCheck (pinvTest ::CM->Bool) | ||
443 | putStrLn "--------- chol ------" | ||
444 | runTestTT $ TestList | ||
445 | [ test "cholR" cholRTest | ||
446 | , test "cholC" cholCTest | ||
447 | ] | ||
448 | putStrLn "--------- qr ---------" | ||
449 | quickCheck (qrTest GSL.qr) | ||
450 | quickCheck (qrTest (GSL.unpackQR . GSL.qrPacked)) | ||
451 | quickCheck (qrTest ( unpackQR . GSL.qrPacked)) | ||
452 | quickCheck (qrTest qr ::RM->Bool) | ||
453 | quickCheck (qrTest qr ::CM->Bool) | ||
454 | putStrLn "--------- hess --------" | ||
455 | quickCheck (hessTest . sqm ::SqM Double->Bool) | ||
456 | quickCheck (hessTest . sqm ::SqM (Complex Double) -> Bool) | ||
457 | putStrLn "--------- schur --------" | ||
458 | quickCheck (schurTest2 . sqm ::SqM Double->Bool) | ||
459 | if os == "mingw32" | ||
460 | then putStrLn "complex schur skipped in this OS" | ||
461 | else quickCheck (schurTest1 . sqm ::SqM (Complex Double) -> Bool) | ||
462 | putStrLn "--------- expm --------" | ||
463 | runTestTT $ TestList | ||
464 | [ test "expmd" (expmTestDiag $ (2><2) [1,2,3,5 :: Double]) | ||
465 | , test "expm1" (expmTest1) | ||
466 | , test "expm2" (expmTest2) | ||
467 | ] | ||
468 | putStrLn "--------- nullspace ------" | ||
469 | quickCheck (nullspaceTest :: RM -> Bool) | ||
470 | quickCheck (nullspaceTest :: CM -> Bool) | ||
471 | putStrLn "--------- vector operations ------" | ||
472 | quickCheck $ (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM)) | ||
473 | quickCheck $ (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM)) | ||
474 | quickCheck $ (\u -> cos u * tan u |~| sin (u::RM)) | ||
475 | quickCheck $ (\u -> (cos u * tan u) :~6~: sin (u::CM)) | ||
476 | runTestTT $ TestList | ||
477 | [ test "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM) | ||
478 | , test "arith2" $ (((1+i) .* ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( (140*i-51).*1 :: CM) | ||
479 | , test "arith3" $ exp (i.*ones(10,10)*pi) + 1 |~| 0 | ||
480 | , test "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3] | ||
481 | ] | ||
482 | putStrLn "--------- GSL ------" | ||
483 | quickCheck $ \v -> ifft (fft v) |~| v | ||
484 | runTestTT $ TestList | ||
485 | [ gammaTest | ||
486 | , besselTest | ||
487 | , exponentialTest | ||
488 | , integrateTest | ||
489 | , polySolveTest | ||
490 | ] | ||
491 | |||
492 | bigtests = do | 23 | bigtests = do |
493 | putStrLn "--------- big matrices -----" | 24 | putStrLn "--------- big matrices -----" |
494 | runTestTT $ TestList | 25 | runTestTT $ TestList |
495 | [ test "eigS" $ eigTestSH bigmat | 26 | [ utest "eigS" $ eigSHProp bigmat |
496 | , test "eigH" $ eigTestSH bigmatc | 27 | , utest "eigH" $ eigSHProp bigmatc |
497 | , test "eigR" $ eigTest bigmat | 28 | , utest "eigR" $ eigProp bigmat |
498 | , test "eigC" $ eigTest bigmatc | 29 | , utest "eigC" $ eigProp bigmatc |
499 | , test "det" $ det (feye 1000) == 1 && det (feye 1002) == -1 | 30 | , utest "det" $ det (feye 1000) == 1 && det (feye 1002) == -1 |
500 | ] | 31 | ] |
32 | return () | ||
501 | 33 | ||
502 | main = do | 34 | main = do |
503 | args <- getArgs | 35 | args <- getArgs |
504 | if "--big" `elem` args | 36 | if "--big" `elem` args |
505 | then bigtests | 37 | then bigtests |
506 | else tests | 38 | else runTests 20 |
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs index 7ecf812..96280fb 100644 --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ b/lib/Numeric/LinearAlgebra/Tests.hs | |||
@@ -23,42 +23,174 @@ import Numeric.LinearAlgebra | |||
23 | import Numeric.LinearAlgebra.Tests.Instances | 23 | import Numeric.LinearAlgebra.Tests.Instances |
24 | import Numeric.LinearAlgebra.Tests.Properties | 24 | import Numeric.LinearAlgebra.Tests.Properties |
25 | import Test.QuickCheck | 25 | import Test.QuickCheck |
26 | import Numeric.GSL(setErrorHandlerOff) | 26 | import Test.HUnit hiding ((~:),test) |
27 | import System.Info | ||
28 | import Data.List(foldl1') | ||
29 | import Numeric.GSL hiding (sin,cos,exp,choose) | ||
27 | 30 | ||
28 | qCheck n = check defaultConfig {configSize = const n} | 31 | qCheck n = check defaultConfig {configSize = const n} |
29 | 32 | ||
33 | utest str b = TestCase $ assertBool str b | ||
34 | |||
35 | feye n = flipud (ident n) :: Matrix Double | ||
36 | |||
37 | detTest1 = det m == 26 | ||
38 | && det mc == 38 :+ (-3) | ||
39 | && det (feye 2) == -1 | ||
40 | where | ||
41 | m = (3><3) | ||
42 | [ 1, 2, 3 | ||
43 | , 4, 5, 7 | ||
44 | , 2, 8, 4 :: Double | ||
45 | ] | ||
46 | mc = (3><3) | ||
47 | [ 1, 2, 3 | ||
48 | , 4, 5, 7 | ||
49 | , 2, 8, i | ||
50 | ] | ||
51 | |||
52 | -------------------------------------------------------------------- | ||
53 | |||
54 | polyEval cs x = foldr (\c ac->ac*x+c) 0 cs | ||
55 | |||
56 | polySolveProp p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p)) | ||
57 | |||
58 | --------------------------------------------------------------------- | ||
59 | |||
60 | quad f a b = fst $ integrateQAGS 1E-9 100 f a b | ||
61 | |||
62 | -- A multiple integral can be easily defined using partial application | ||
63 | quad2 f a b g1 g2 = quad h a b | ||
64 | where h x = quad (f x) (g1 x) (g2 x) | ||
65 | |||
66 | volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) | ||
67 | 0 r (const 0) (\x->sqrt (r*r-x*x)) | ||
68 | |||
69 | --------------------------------------------------------------------- | ||
70 | |||
71 | besselTest = utest "bessel_J0_e" ( abs (r-expected) < e ) | ||
72 | where (r,e) = bessel_J0_e 5.0 | ||
73 | expected = -0.17759677131433830434739701 | ||
74 | |||
75 | exponentialTest = utest "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 ) | ||
76 | where (v,e,err) = exp_e10_e 30.0 | ||
77 | expected = exp 30.0 | ||
78 | |||
79 | --------------------------------------------------------------------- | ||
80 | |||
81 | nd1 = (3><3) [ 1/2, 1/4, 1/4 | ||
82 | , 0/1, 1/2, 1/4 | ||
83 | , 1/2, 1/4, 1/2 :: Double] | ||
84 | |||
85 | nd2 = (2><2) [1, 0, 1, 1:: Complex Double] | ||
86 | |||
87 | expmTest1 = expm nd1 :~14~: (3><3) | ||
88 | [ 1.762110887278176 | ||
89 | , 0.478085470590435 | ||
90 | , 0.478085470590435 | ||
91 | , 0.104719410945666 | ||
92 | , 1.709751181805343 | ||
93 | , 0.425725765117601 | ||
94 | , 0.851451530235203 | ||
95 | , 0.530445176063267 | ||
96 | , 1.814470592751009 ] | ||
97 | |||
98 | expmTest2 = expm nd2 :~15~: (2><2) | ||
99 | [ 2.718281828459045 | ||
100 | , 0.000000000000000 | ||
101 | , 2.718281828459045 | ||
102 | , 2.718281828459045 ] | ||
103 | |||
104 | |||
105 | --------------------------------------------------------------------- | ||
106 | |||
107 | rot :: Double -> Matrix Double | ||
108 | rot a = (3><3) [ c,0,s | ||
109 | , 0,1,0 | ||
110 | ,-s,0,c ] | ||
111 | where c = cos a | ||
112 | s = sin a | ||
113 | |||
114 | rotTest = fun (10^5) :~12~: rot 5E4 | ||
115 | where fun n = foldl1' (<>) (map rot angles) | ||
116 | where angles = toList $ linspace n (0,1) | ||
117 | |||
118 | |||
30 | -- | It runs all the tests. | 119 | -- | It runs all the tests. |
31 | runTests :: Int -- ^ maximum dimension | 120 | runTests :: Int -- ^ maximum dimension |
32 | -> IO () | 121 | -> IO () |
33 | runTests n = do | 122 | runTests n = do |
34 | setErrorHandlerOff | 123 | setErrorHandlerOff |
35 | let test p = qCheck n p | 124 | let test p = qCheck n p |
125 | putStrLn "------ lu" | ||
36 | test (luProp . rM) | 126 | test (luProp . rM) |
37 | test (luProp . cM) | 127 | test (luProp . cM) |
128 | putStrLn "------ inv" | ||
38 | test (invProp . rSqWC) | 129 | test (invProp . rSqWC) |
39 | test (invProp . cSqWC) | 130 | test (invProp . cSqWC) |
131 | putStrLn "------ pinv" | ||
40 | test (pinvProp . rM) | 132 | test (pinvProp . rM) |
41 | test (pinvProp . cM) | 133 | if os == "mingw32" |
134 | then putStrLn "complex pinvTest skipped in this OS" | ||
135 | else test (pinvProp . cM) | ||
136 | putStrLn "------ det" | ||
137 | test (detProp . rSqWC) | ||
42 | test (detProp . cSqWC) | 138 | test (detProp . cSqWC) |
139 | putStrLn "------ svd" | ||
43 | test (svdProp1 . rM) | 140 | test (svdProp1 . rM) |
44 | test (svdProp1 . cM) | 141 | test (svdProp1 . cM) |
45 | test (svdProp2 . rM) | 142 | test (svdProp2 . rM) |
46 | test (svdProp2 . cM) | 143 | test (svdProp2 . cM) |
144 | putStrLn "------ eig" | ||
47 | test (eigSHProp . rHer) | 145 | test (eigSHProp . rHer) |
48 | test (eigSHProp . cHer) | 146 | test (eigSHProp . cHer) |
49 | test (eigProp . rSq) | 147 | test (eigProp . rSq) |
50 | test (eigProp . cSq) | 148 | test (eigProp . cSq) |
149 | putStrLn "------ nullSpace" | ||
51 | test (nullspaceProp . rM) | 150 | test (nullspaceProp . rM) |
52 | test (nullspaceProp . cM) | 151 | test (nullspaceProp . cM) |
152 | putStrLn "------ qr" | ||
53 | test (qrProp . rM) | 153 | test (qrProp . rM) |
54 | test (qrProp . cM) | 154 | test (qrProp . cM) |
155 | putStrLn "------ hess" | ||
55 | test (hessProp . rSq) | 156 | test (hessProp . rSq) |
56 | test (hessProp . cSq) | 157 | test (hessProp . cSq) |
158 | putStrLn "------ schur" | ||
57 | test (schurProp2 . rSq) | 159 | test (schurProp2 . rSq) |
58 | test (schurProp1 . cSq) | 160 | if os == "mingw32" |
161 | then putStrLn "complex schur skipped in this OS" | ||
162 | else test (schurProp1 . cSq) | ||
163 | putStrLn "------ chol" | ||
59 | test (cholProp . rPosDef) | 164 | test (cholProp . rPosDef) |
60 | test (cholProp . cPosDef) | 165 | test (cholProp . cPosDef) |
166 | putStrLn "------ expm" | ||
167 | test (expmDiagProp . rSqWC) | ||
168 | test (expmDiagProp . cSqWC) | ||
169 | putStrLn "------ fft" | ||
170 | test (\v -> ifft (fft v) |~| v) | ||
171 | putStrLn "------ vector operations" | ||
172 | test (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM)) | ||
173 | test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM)) | ||
174 | test (\u -> cos u * tan u |~| sin (u::RM)) | ||
175 | test (\u -> (cos u * tan u) |~| sin (u::CM)) | ||
176 | putStrLn "------ some unit tests" | ||
177 | runTestTT $ TestList | ||
178 | [ utest "1E5 rots" rotTest | ||
179 | , utest "det1" detTest1 | ||
180 | , utest "expm1" (expmTest1) | ||
181 | , utest "expm2" (expmTest2) | ||
182 | , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM) | ||
183 | , utest "arith2" $ (((1+i) .* ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( (140*i-51).*1 :: CM) | ||
184 | , utest "arith3" $ exp (i.*ones(10,10)*pi) + 1 |~| 0 | ||
185 | , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3] | ||
186 | , utest "gamma" (gamma 5 == 24.0) | ||
187 | , besselTest | ||
188 | , exponentialTest | ||
189 | , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < 1E-8) | ||
190 | , utest "polySolve" (polySolveProp [1,2,3,4]) | ||
191 | ] | ||
192 | return () | ||
61 | 193 | ||
62 | -- | Some additional tests on big matrices. They take a few minutes. | 194 | -- -- | Some additional tests on big matrices. They take a few minutes. |
63 | runBigTests :: IO () | 195 | -- runBigTests :: IO () |
64 | runBigTests = undefined | 196 | -- runBigTests = undefined |
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs index 0317469..0563e62 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs | |||
@@ -9,13 +9,33 @@ Maintainer : Alberto Ruiz (aruiz at um dot es) | |||
9 | Stability : provisional | 9 | Stability : provisional |
10 | Portability : portable | 10 | Portability : portable |
11 | 11 | ||
12 | Arbitrary instances for vectors, matrices. | 12 | Testing properties. |
13 | 13 | ||
14 | -} | 14 | -} |
15 | 15 | ||
16 | module Numeric.LinearAlgebra.Tests.Properties | 16 | module Numeric.LinearAlgebra.Tests.Properties ( |
17 | 17 | dist, (|~|), (~:), Aprox((:~)), | |
18 | where | 18 | zeros, ones, |
19 | square, | ||
20 | unitary, | ||
21 | hermitian, | ||
22 | wellCond, | ||
23 | positiveDefinite, | ||
24 | upperTriang, | ||
25 | upperHessenberg, | ||
26 | luProp, | ||
27 | invProp, | ||
28 | pinvProp, | ||
29 | detProp, | ||
30 | nullspaceProp, | ||
31 | svdProp1, svdProp2, | ||
32 | eigProp, eigSHProp, | ||
33 | qrProp, | ||
34 | hessProp, | ||
35 | schurProp1, schurProp2, | ||
36 | cholProp, | ||
37 | expmDiagProp | ||
38 | ) where | ||
19 | 39 | ||
20 | import Numeric.LinearAlgebra | 40 | import Numeric.LinearAlgebra |
21 | import Numeric.LinearAlgebra.Tests.Instances(Sq(..),Her(..),Rot(..)) | 41 | import Numeric.LinearAlgebra.Tests.Instances(Sq(..),Her(..),Rot(..)) |
@@ -50,8 +70,6 @@ unitary m = square m && m <> ctrans m |~| ident (rows m) | |||
50 | 70 | ||
51 | hermitian m = square m && m |~| ctrans m | 71 | hermitian m = square m && m |~| ctrans m |
52 | 72 | ||
53 | degenerate m = rank m < min (rows m) (cols m) | ||
54 | |||
55 | wellCond m = rcond m > 1/100 | 73 | wellCond m = rcond m > 1/100 |
56 | 74 | ||
57 | positiveDefinite m = minimum (toList e) > 0 | 75 | positiveDefinite m = minimum (toList e) > 0 |
@@ -125,3 +143,7 @@ schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fix | |||
125 | cholProp m = m |~| ctrans c <> c && upperTriang c | 143 | cholProp m = m |~| ctrans c <> c && upperTriang c |
126 | where c = chol m | 144 | where c = chol m |
127 | pos = positiveDefinite m | 145 | pos = positiveDefinite m |
146 | |||
147 | expmDiagProp m = expm (logm m) |~| complex m | ||
148 | where logm m = matFunc log m | ||
149 | |||