summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--examples/tests.hs490
-rw-r--r--lib/Numeric/LinearAlgebra/Tests.hs144
-rw-r--r--lib/Numeric/LinearAlgebra/Tests/Properties.hs34
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
3module Main where 1module Main where
4 2
5import Numeric.GSL hiding (sin,cos,exp,choose)
6import Numeric.LinearAlgebra 3import Numeric.LinearAlgebra
7import Numeric.LinearAlgebra.LAPACK 4import Numeric.LinearAlgebra.Tests
8import qualified Numeric.GSL.Matrix as GSL
9import Test.QuickCheck hiding (test)
10import Test.HUnit hiding ((~:),test)
11import System.Random(randomRs,mkStdGen) 5import System.Random(randomRs,mkStdGen)
12import System.Info 6import Test.HUnit hiding (test)
13import Data.List(foldl1', transpose)
14import System(getArgs) 7import System(getArgs)
15import Debug.Trace(trace)
16
17debug x = trace (show x) x
18
19type RM = Matrix Double
20type CM = Matrix (Complex Double)
21
22-- relative error
23dist :: (Normed t, Num t) => t -> t -> Double
24dist 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
35infixl 4 |~|
36a |~| b = a :~10~: b
37
38data Aprox a = (:~) a Int
39
40(~:) :: (Normed a, Num a) => Aprox a -> a -> Bool
41a :~n~: b = dist a b < 10^^(-n)
42
43
44maxdim = 10
45
46instance (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
53instance (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
65data PairM a = PairM (Matrix a) (Matrix a) deriving Show
66instance (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
77data SqM a = SqM (Matrix a) deriving Show
78sqm (SqM a) = a
79instance (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
86data Sym a = Sym (Matrix a) deriving Show
87sym (Sym a) = a
88instance (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
94data Her = Her (Matrix (Complex Double)) deriving Show
95her (Her a) = a
96instance {-(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
102data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show 10pseudorandomR seed (n,m) = reshape m $ fromList $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed
103instance (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
113instance (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
120data PairV a = PairV (Vector a) (Vector a)
121instance (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
131test str b = TestCase $ assertBool str b
132
133----------------------------------------------------------------------
134
135pseudorandomR seed (n,m) = reshape m $ fromList $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed
136 11
137pseudorandomC seed (n,m) = toComplex (pseudorandomR seed (n,m), pseudorandomR (seed+1) (n,m)) 12pseudorandomC seed (n,m) = toComplex (pseudorandomR seed (n,m), pseudorandomR (seed+1) (n,m))
138 13
@@ -141,366 +16,23 @@ bigmat = m + trans m :: RM
141bigmatc = mc + ctrans mc ::CM 16bigmatc = mc + ctrans mc ::CM
142 where mc = pseudorandomC 19 (1000,1000) 17 where mc = pseudorandomC 19 (1000,1000)
143 18
144---------------------------------------------------------------------- 19utest str b = TestCase $ assertBool str b
145
146
147m = (3><3)
148 [ 1, 2, 3
149 , 4, 5, 7
150 , 2, 8, 4 :: Double
151 ]
152
153mc = (3><3)
154 [ 1, 2, 3
155 , 4, 5, 7
156 , 2, 8, i
157 ]
158
159
160mr = (3><4)
161 [ 1, 2, 3, 4,
162 2, 4, 6, 8,
163 1, 1, 1, 2:: Double
164 ]
165
166mrc = (3><4)
167 [ 1, 2, 3, 4,
168 2, 4, 6, 8,
169 i, i, i, 2
170 ]
171
172a = (3><4)
173 [ 1, 0, 0, 0
174 , 0, 2, 0, 0
175 , 0, 0, 0, 0 :: Double
176 ]
177
178b = (3><4)
179 [ 1, 0, 0, 0
180 , 0, 2, 3, 0
181 , 0, 0, 4, 0 :: Double
182 ]
183
184ac = (2><3) [1 .. 6::Double]
185bc = (3><4) [7 .. 18::Double]
186
187af = (2>|<3) [1,4,2,5,3,6::Double]
188bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double]
189
190-------------------------------------------------------
191 20
192feye n = flipud (ident n) :: Matrix Double 21feye n = flipud (ident n) :: Matrix Double
193 22
194
195luTest1 m = m |~| p <> l <> u
196 where (l,u,p,_) = lu m
197
198detTest1 = det m == 26
199 && det mc == 38 :+ (-3)
200 && det (feye 2) == -1
201
202detTest2 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
209invTest m = degenerate m || m <> inv m |~| ident (rows m)
210
211pinvTest m = m <> p <> m |~| m
212 && p <> m <> p |~| p
213 && hermitian (m<>p)
214 && hermitian (p<>m)
215 where p = pinv m
216
217square m = rows m == cols m
218
219unitary m = square m && m <> ctrans m |~| ident (rows m)
220
221hermitian m = m |~| ctrans m
222
223upperTriang m = rows m == 1 || down == z
224 where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
225 z = constant 0 (dim down)
226
227upperHessenberg m = rows m < 3 || down == z
228 where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))
229 z = constant 0 (dim down)
230
231svdTest svd m = u <> real d <> trans v |~| m
232 && unitary u && unitary v
233 where (u,d,v) = full svd m
234
235svdTest' svd m = m |~| 0 || u <> real (diag s) <> trans v |~| m
236 where (u,s,v) = economy svd m
237
238eigTest m = complex m <> v |~| v <> diag s
239 where (s, v) = eig m
240
241eigTestSH 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
246zeros (r,c) = reshape c (constant 0 (r*c))
247
248ones (r,c) = zeros (r,c) + 1
249
250degenerate m = rank m < min (rows m) (cols m)
251
252prec = 1E-15
253
254singular m = s1 < prec || s2/s1 < prec
255 where (_,ss,_) = svd m
256 s = toList ss
257 s1 = maximum s
258 s2 = minimum s
259
260nullspaceTest 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
268polyEval cs x = foldr (\c ac->ac*x+c) 0 cs
269
270polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))
271
272
273polySolveTest = test "polySolve" (polySolveTest' [1,2,3,4])
274
275---------------------------------------------------------------------
276
277quad f a b = fst $ integrateQAGS 1E-9 100 f a b
278
279-- A multiple integral can be easily defined using partial application
280quad2 f a b g1 g2 = quad h a b
281 where h x = quad (f x) (g1 x) (g2 x)
282
283volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y))
284 0 r (const 0) (\x->sqrt (r*r-x*x))
285
286epsTol = 1E-8::Double
287
288integrateTest = test "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < epsTol)
289
290---------------------------------------------------------------------
291
292besselTest = test "bessel_J0_e" ( abs (r-expected) < e )
293 where (r,e) = bessel_J0_e 5.0
294 expected = -0.17759677131433830434739701
295
296exponentialTest = 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
300gammaTest = test "gamma" (gamma 5 == 24.0)
301
302---------------------------------------------------------------------
303
304cholRTest = chol ((2><2) [1,2,2,9::Double]) == (2><2) [1,2,0,2.23606797749979]
305cholCTest = chol ((2><2) [1,2,2,9::Complex Double]) == (2><2) [1,2,0,2.23606797749979]
306
307---------------------------------------------------------------------
308
309qrTest qr m = q <> r |~| m && unitary q && upperTriang r
310 where (q,r) = qr m
311
312---------------------------------------------------------------------
313
314hessTest m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
315 where (p,h) = hess m
316
317---------------------------------------------------------------------
318
319schurTest1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s
320 where (u,s) = schur m
321
322schurTest2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme
323 where (u,s) = schur m
324
325---------------------------------------------------------------------
326
327nd1 = (3><3) [ 1/2, 1/4, 1/4
328 , 0/1, 1/2, 1/4
329 , 1/2, 1/4, 1/2 :: Double]
330
331nd2 = (2><2) [1, 0, 1, 1:: Complex Double]
332
333expmTest1 = 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
344expmTest2 = expm nd2 :~15~: (2><2)
345 [ 2.718281828459045
346 , 0.000000000000000
347 , 2.718281828459045
348 , 2.718281828459045 ]
349
350expmTestDiag m = expm (logm m) |~| complex m
351 where logm m = matFunc Prelude.log m
352
353
354
355---------------------------------------------------------------------
356
357asFortran m = (rows m >|< cols m) $ toList (flatten $ trans m)
358asC m = (rows m >< cols m) $ toList (flatten m)
359
360mulC a b = a <> b
361mulF a b = trans $ trans b <> trans a
362
363-------------------------------------------------------------------------
364
365multiplyG 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
369r >|< 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
377rot :: Double -> Matrix Double
378rot a = (3><3) [ c,0,s
379 , 0,1,0
380 ,-s,0,c ]
381 where c = cos a
382 s = sin a
383
384fun n = foldl1' (<>) (map rot angles)
385 where angles = toList $ linspace n (0,1)
386
387rotTest = fun (10^5) :~12~: rot 5E4
388
389---------------------------------------------------------------------
390
391tests = 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
492bigtests = do 23bigtests = 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
502main = do 34main = 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
23import Numeric.LinearAlgebra.Tests.Instances 23import Numeric.LinearAlgebra.Tests.Instances
24import Numeric.LinearAlgebra.Tests.Properties 24import Numeric.LinearAlgebra.Tests.Properties
25import Test.QuickCheck 25import Test.QuickCheck
26import Numeric.GSL(setErrorHandlerOff) 26import Test.HUnit hiding ((~:),test)
27import System.Info
28import Data.List(foldl1')
29import Numeric.GSL hiding (sin,cos,exp,choose)
27 30
28qCheck n = check defaultConfig {configSize = const n} 31qCheck n = check defaultConfig {configSize = const n}
29 32
33utest str b = TestCase $ assertBool str b
34
35feye n = flipud (ident n) :: Matrix Double
36
37detTest1 = 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
54polyEval cs x = foldr (\c ac->ac*x+c) 0 cs
55
56polySolveProp p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))
57
58---------------------------------------------------------------------
59
60quad f a b = fst $ integrateQAGS 1E-9 100 f a b
61
62-- A multiple integral can be easily defined using partial application
63quad2 f a b g1 g2 = quad h a b
64 where h x = quad (f x) (g1 x) (g2 x)
65
66volSphere 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
71besselTest = utest "bessel_J0_e" ( abs (r-expected) < e )
72 where (r,e) = bessel_J0_e 5.0
73 expected = -0.17759677131433830434739701
74
75exponentialTest = 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
81nd1 = (3><3) [ 1/2, 1/4, 1/4
82 , 0/1, 1/2, 1/4
83 , 1/2, 1/4, 1/2 :: Double]
84
85nd2 = (2><2) [1, 0, 1, 1:: Complex Double]
86
87expmTest1 = 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
98expmTest2 = expm nd2 :~15~: (2><2)
99 [ 2.718281828459045
100 , 0.000000000000000
101 , 2.718281828459045
102 , 2.718281828459045 ]
103
104
105---------------------------------------------------------------------
106
107rot :: Double -> Matrix Double
108rot a = (3><3) [ c,0,s
109 , 0,1,0
110 ,-s,0,c ]
111 where c = cos a
112 s = sin a
113
114rotTest = 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.
31runTests :: Int -- ^ maximum dimension 120runTests :: Int -- ^ maximum dimension
32 -> IO () 121 -> IO ()
33runTests n = do 122runTests 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.
63runBigTests :: IO () 195-- runBigTests :: IO ()
64runBigTests = 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)
9Stability : provisional 9Stability : provisional
10Portability : portable 10Portability : portable
11 11
12Arbitrary instances for vectors, matrices. 12Testing properties.
13 13
14-} 14-}
15 15
16module Numeric.LinearAlgebra.Tests.Properties 16module Numeric.LinearAlgebra.Tests.Properties (
17 17 dist, (|~|), (~:), Aprox((:~)),
18where 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
20import Numeric.LinearAlgebra 40import Numeric.LinearAlgebra
21import Numeric.LinearAlgebra.Tests.Instances(Sq(..),Her(..),Rot(..)) 41import Numeric.LinearAlgebra.Tests.Instances(Sq(..),Her(..),Rot(..))
@@ -50,8 +70,6 @@ unitary m = square m && m <> ctrans m |~| ident (rows m)
50 70
51hermitian m = square m && m |~| ctrans m 71hermitian m = square m && m |~| ctrans m
52 72
53degenerate m = rank m < min (rows m) (cols m)
54
55wellCond m = rcond m > 1/100 73wellCond m = rcond m > 1/100
56 74
57positiveDefinite m = minimum (toList e) > 0 75positiveDefinite m = minimum (toList e) > 0
@@ -125,3 +143,7 @@ schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fix
125cholProp m = m |~| ctrans c <> c && upperTriang c 143cholProp 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
147expmDiagProp m = expm (logm m) |~| complex m
148 where logm m = matFunc log m
149