summaryrefslogtreecommitdiff
path: root/packages/tests/src/Numeric
diff options
context:
space:
mode:
Diffstat (limited to 'packages/tests/src/Numeric')
-rw-r--r--packages/tests/src/Numeric/GSL/Tests.hs4
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests.hs202
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs91
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs153
4 files changed, 183 insertions, 267 deletions
diff --git a/packages/tests/src/Numeric/GSL/Tests.hs b/packages/tests/src/Numeric/GSL/Tests.hs
index 9dff6f5..e5d205d 100644
--- a/packages/tests/src/Numeric/GSL/Tests.hs
+++ b/packages/tests/src/Numeric/GSL/Tests.hs
@@ -19,7 +19,7 @@ import System.Exit (exitFailure)
19 19
20import Test.HUnit (runTestTT, failures, Test(..), errors) 20import Test.HUnit (runTestTT, failures, Test(..), errors)
21 21
22import Numeric.LinearAlgebra 22import Numeric.LinearAlgebra.HMatrix
23import Numeric.GSL 23import Numeric.GSL
24import Numeric.LinearAlgebra.Tests (qCheck, utest) 24import Numeric.LinearAlgebra.Tests (qCheck, utest)
25import Numeric.LinearAlgebra.Tests.Properties ((|~|), (~~)) 25import Numeric.LinearAlgebra.Tests.Properties ((|~|), (~~))
@@ -42,7 +42,7 @@ fittingTest = utest "levmar" (ok1 && ok2)
42 sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] 42 sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
43 43
44 ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d 44 ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d
45 ok2 = norm2 (fromList (map fst sols) - fromList sol) < 1E-5 45 ok2 = norm_2 (fromList (map fst sols) - fromList sol) < 1E-5
46 46
47--------------------------------------------------------------------- 47---------------------------------------------------------------------
48 48
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 713af79..71c7c16 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -28,12 +28,9 @@ module Numeric.LinearAlgebra.Tests(
28--, runBigTests 28--, runBigTests
29) where 29) where
30 30
31import Numeric.LinearAlgebra 31import Numeric.LinearAlgebra hiding (unitary)
32import Numeric.LinearAlgebra.HMatrix hiding ((<>),linearSolve) 32import Numeric.LinearAlgebra.Devel hiding (vec)
33import Numeric.LinearAlgebra.Static(L) 33import Numeric.LinearAlgebra.Static(L)
34import Numeric.LinearAlgebra.Util(col,row)
35import Data.Packed
36import Numeric.LinearAlgebra.LAPACK
37import Numeric.LinearAlgebra.Tests.Instances 34import Numeric.LinearAlgebra.Tests.Instances
38import Numeric.LinearAlgebra.Tests.Properties 35import Numeric.LinearAlgebra.Tests.Properties
39import Test.HUnit hiding ((~:),test,Testable,State) 36import Test.HUnit hiding ((~:),test,Testable,State)
@@ -44,16 +41,13 @@ import qualified Prelude
44import System.CPUTime 41import System.CPUTime
45import System.Exit 42import System.Exit
46import Text.Printf 43import Text.Printf
47import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr) 44import Numeric.LinearAlgebra.Devel(unsafeFromForeignPtr,unsafeToForeignPtr)
48import Control.Arrow((***)) 45import Control.Arrow((***))
49import Debug.Trace 46import Debug.Trace
50import Control.Monad(when) 47import Control.Monad(when)
51import Numeric.LinearAlgebra.Util hiding (ones,row,col)
52import Control.Applicative 48import Control.Applicative
53import Control.Monad(ap) 49import Control.Monad(ap)
54 50
55import Data.Packed.ST
56
57import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector 51import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
58 ,sized,classify,Testable,Property 52 ,sized,classify,Testable,Property
59 ,quickCheckWithResult,maxSize,stdArgs,shrink) 53 ,quickCheckWithResult,maxSize,stdArgs,shrink)
@@ -89,7 +83,7 @@ detTest1 = det m == 26
89 mc = (3><3) 83 mc = (3><3)
90 [ 1, 2, 3 84 [ 1, 2, 3
91 , 4, 5, 7 85 , 4, 5, 7
92 , 2, 8, i 86 , 2, 8, iC
93 ] 87 ]
94 88
95detTest2 = inv1 |~| inv2 && [det1] ~~ [det2] 89detTest2 = inv1 |~| inv2 && [det1] ~~ [det2]
@@ -140,7 +134,7 @@ randomTestGaussian = c :~1~: snd (meanCov dat) where
140 2,4,0, 134 2,4,0,
141 -2,2,1] 135 -2,2,1]
142 m = 3 |> [1,2,3] 136 m = 3 |> [1,2,3]
143 c = a <> trans a 137 c = a <> tr a
144 dat = gaussianSample 7 (10^6) m c 138 dat = gaussianSample 7 (10^6) m c
145 139
146randomTestUniform = c :~1~: snd (meanCov dat) where 140randomTestUniform = c :~1~: snd (meanCov dat) where
@@ -174,54 +168,54 @@ offsetTest = y == y' where
174 168
175normsVTest = TestList [ 169normsVTest = TestList [
176 utest "normv2CD" $ norm2PropC v 170 utest "normv2CD" $ norm2PropC v
177 , utest "normv2CF" $ norm2PropC (single v) 171-- , utest "normv2CF" $ norm2PropC (single v)
178#ifndef NONORMVTEST 172#ifndef NONORMVTEST
179 , utest "normv2D" $ norm2PropR x 173 , utest "normv2D" $ norm2PropR x
180 , utest "normv2F" $ norm2PropR (single x) 174-- , utest "normv2F" $ norm2PropR (single x)
181#endif 175#endif
182 , utest "normv1CD" $ norm1 v == 8 176 , utest "normv1CD" $ norm_1 v == 8
183 , utest "normv1CF" $ norm1 (single v) == 8 177-- , utest "normv1CF" $ norm_1 (single v) == 8
184 , utest "normv1D" $ norm1 x == 6 178 , utest "normv1D" $ norm_1 x == 6
185 , utest "normv1F" $ norm1 (single x) == 6 179-- , utest "normv1F" $ norm_1 (single x) == 6
186 180
187 , utest "normvInfCD" $ normInf v == 5 181 , utest "normvInfCD" $ norm_Inf v == 5
188 , utest "normvInfCF" $ normInf (single v) == 5 182-- , utest "normvInfCF" $ norm_Inf (single v) == 5
189 , utest "normvInfD" $ normInf x == 3 183 , utest "normvInfD" $ norm_Inf x == 3
190 , utest "normvInfF" $ normInf (single x) == 3 184-- , utest "normvInfF" $ norm_Inf (single x) == 3
191 185
192 ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double) 186 ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double)
193 x = fromList [1,2,-3] :: Vector Double 187 x = fromList [1,2,-3] :: Vector Double
194#ifndef NONORMVTEST 188#ifndef NONORMVTEST
195 norm2PropR a = norm2 a =~= sqrt (udot a a) 189 norm2PropR a = norm_2 a =~= sqrt (udot a a)
196#endif 190#endif
197 norm2PropC a = norm2 a =~= realPart (sqrt (a <.> a)) 191 norm2PropC a = norm_2 a =~= realPart (sqrt (a `dot` a))
198 a =~= b = fromList [a] |~| fromList [b] 192 a =~= b = fromList [a] |~| fromList [b]
199 193
200normsMTest = TestList [ 194normsMTest = TestList [
201 utest "norm2mCD" $ pnorm PNorm2 v =~= 8.86164970498005 195 utest "norm2mCD" $ norm_2 v =~= 8.86164970498005
202 , utest "norm2mCF" $ pnorm PNorm2 (single v) =~= 8.86164970498005 196-- , utest "norm2mCF" $ norm_2 (single v) =~= 8.86164970498005
203 , utest "norm2mD" $ pnorm PNorm2 x =~= 5.96667765076216 197 , utest "norm2mD" $ norm_2 x =~= 5.96667765076216
204 , utest "norm2mF" $ pnorm PNorm2 (single x) =~= 5.96667765076216 198-- , utest "norm2mF" $ norm_2 (single x) =~= 5.96667765076216
205 199
206 , utest "norm1mCD" $ pnorm PNorm1 v == 9 200 , utest "norm1mCD" $ norm_1 v == 9
207 , utest "norm1mCF" $ pnorm PNorm1 (single v) == 9 201-- , utest "norm1mCF" $ norm_1 (single v) == 9
208 , utest "norm1mD" $ pnorm PNorm1 x == 7 202 , utest "norm1mD" $ norm_1 x == 7
209 , utest "norm1mF" $ pnorm PNorm1 (single x) == 7 203-- , utest "norm1mF" $ norm_1 (single x) == 7
210 204
211 , utest "normmInfCD" $ pnorm Infinity v == 12 205 , utest "normmInfCD" $ norm_Inf v == 12
212 , utest "normmInfCF" $ pnorm Infinity (single v) == 12 206-- , utest "normmInfCF" $ norm_Inf (single v) == 12
213 , utest "normmInfD" $ pnorm Infinity x == 8 207 , utest "normmInfD" $ norm_Inf x == 8
214 , utest "normmInfF" $ pnorm Infinity (single x) == 8 208-- , utest "normmInfF" $ norm_Inf (single x) == 8
215 209
216 , utest "normmFroCD" $ pnorm Frobenius v =~= 8.88819441731559 210 , utest "normmFroCD" $ norm_Frob v =~= 8.88819441731559
217 , utest "normmFroCF" $ pnorm Frobenius (single v) =~~= 8.88819441731559 211-- , utest "normmFroCF" $ norm_Frob (single v) =~~= 8.88819441731559
218 , utest "normmFroD" $ pnorm Frobenius x =~= 6.24499799839840 212 , utest "normmFroD" $ norm_Frob x =~= 6.24499799839840
219 , utest "normmFroF" $ pnorm Frobenius (single x) =~~= 6.24499799839840 213-- , utest "normmFroF" $ norm_Frob (single x) =~~= 6.24499799839840
220 214
221 ] where v = (2><2) [1,-2*i,3:+4,7] :: Matrix (Complex Double) 215 ] where v = (2><2) [1,-2*iC,3:+4,7] :: Matrix (Complex Double)
222 x = (2><2) [1,2,-3,5] :: Matrix Double 216 x = (2><2) [1,2,-3,5] :: Matrix Double
223 a =~= b = fromList [a] :~10~: fromList [b] 217 a =~= b = fromList [a] :~10~: fromList [b]
224 a =~~= b = fromList [a] :~5~: fromList [b] 218-- a =~~= b = fromList [a] :~5~: fromList [b]
225 219
226--------------------------------------------------------------------- 220---------------------------------------------------------------------
227 221
@@ -236,7 +230,7 @@ sumprodTest = TestList [
236 , utest "prodD" $ prodProp v 230 , utest "prodD" $ prodProp v
237 , utest "prodF" $ prodProp (single v) 231 , utest "prodF" $ prodProp (single v)
238 ] where v = fromList [1,2,3] :: Vector Double 232 ] where v = fromList [1,2,3] :: Vector Double
239 z = fromList [1,2-i,3+i] 233 z = fromList [1,2-iC,3+iC]
240 prodProp x = prodElements x == product (toList x) 234 prodProp x = prodElements x == product (toList x)
241 235
242--------------------------------------------------------------------- 236---------------------------------------------------------------------
@@ -250,7 +244,7 @@ chainTest = utest "chain" $ foldl1' (<>) ms |~| optimiseMult ms where
250 244
251--------------------------------------------------------------------- 245---------------------------------------------------------------------
252 246
253conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m) 247conjuTest m = cmap conjugate (flatten (conj (tr m))) == flatten (tr m)
254 248
255--------------------------------------------------------------------- 249---------------------------------------------------------------------
256 250
@@ -306,7 +300,7 @@ lift_maybe m = MaybeT $ do
306 300
307-- apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs 301-- apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs
308--successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool 302--successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool
309successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (dim v - 1) v))) (v @> 0) 303successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (size v - 1) v))) (v ! 0)
310 where stp e = do 304 where stp e = do
311 ep <- lift_maybe $ state_get 305 ep <- lift_maybe $ state_get
312 if t e ep 306 if t e ep
@@ -315,7 +309,7 @@ successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ s
315 309
316-- operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input 310-- operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input
317--successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b 311--successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b
318successive f v = evalState (mapVectorM stp (subVector 1 (dim v - 1) v)) (v @> 0) 312successive f v = evalState (mapVectorM stp (subVector 1 (size v - 1) v)) (v ! 0)
319 where stp e = do 313 where stp e = do
320 ep <- state_get 314 ep <- state_get
321 state_put e 315 state_put e
@@ -377,23 +371,6 @@ convolutionTest = utest "convolution" ok
377 371
378-------------------------------------------------------------------------------- 372--------------------------------------------------------------------------------
379 373
380kroneckerTest = utest "kronecker" ok
381 where
382 a,x,b :: Matrix Double
383 a = (3><4) [1..]
384 x = (4><2) [3,5..]
385 b = (2><5) [0,5..]
386 v1 = vec (a <> x <> b)
387 v2 = (trans b `kronecker` a) <> vec x
388 s = trans b <> b
389 v3 = vec s
390 v4 = (dup 5 :: Matrix Double) <> vech s
391 ok = v1 == v2 && v3 == v4
392 && vtrans 1 a == trans a
393 && vtrans (rows a) a == asColumn (vec a)
394
395--------------------------------------------------------------------------------
396
397sparseTest = utest "sparse" (fst $ checkT (undefined :: GMatrix)) 374sparseTest = utest "sparse" (fst $ checkT (undefined :: GMatrix))
398 375
399-------------------------------------------------------------------------------- 376--------------------------------------------------------------------------------
@@ -435,11 +412,11 @@ runTests n = do
435 test (multProp1 10 . cConsist) 412 test (multProp1 10 . cConsist)
436 test (multProp2 10 . rConsist) 413 test (multProp2 10 . rConsist)
437 test (multProp2 10 . cConsist) 414 test (multProp2 10 . cConsist)
438 putStrLn "------ mult Float" 415-- putStrLn "------ mult Float"
439 test (multProp1 6 . (single *** single) . rConsist) 416-- test (multProp1 6 . (single *** single) . rConsist)
440 test (multProp1 6 . (single *** single) . cConsist) 417-- test (multProp1 6 . (single *** single) . cConsist)
441 test (multProp2 6 . (single *** single) . rConsist) 418-- test (multProp2 6 . (single *** single) . rConsist)
442 test (multProp2 6 . (single *** single) . cConsist) 419-- test (multProp2 6 . (single *** single) . cConsist)
443 putStrLn "------ sub-trans" 420 putStrLn "------ sub-trans"
444 test (subProp . rM) 421 test (subProp . rM)
445 test (subProp . cM) 422 test (subProp . cM)
@@ -472,16 +449,16 @@ runTests n = do
472 putStrLn "------ svd" 449 putStrLn "------ svd"
473 test (svdProp1 . rM) 450 test (svdProp1 . rM)
474 test (svdProp1 . cM) 451 test (svdProp1 . cM)
475 test (svdProp1a svdR) 452 test (svdProp1a svd . rM)
476 test (svdProp1a svdC) 453 test (svdProp1a svd . cM)
477 test (svdProp1a svdRd) 454-- test (svdProp1a svdRd)
478 test (svdProp1b svdR) 455 test (svdProp1b svd . rM)
479 test (svdProp1b svdC) 456 test (svdProp1b svd . cM)
480 test (svdProp1b svdRd) 457-- test (svdProp1b svdRd)
481 test (svdProp2 thinSVDR) 458 test (svdProp2 thinSVD . rM)
482 test (svdProp2 thinSVDC) 459 test (svdProp2 thinSVD . cM)
483 test (svdProp2 thinSVDRd) 460-- test (svdProp2 thinSVDRd)
484 test (svdProp2 thinSVDCd) 461-- test (svdProp2 thinSVDCd)
485 test (svdProp3 . rM) 462 test (svdProp3 . rM)
486 test (svdProp3 . cM) 463 test (svdProp3 . cM)
487 test (svdProp4 . rM) 464 test (svdProp4 . rM)
@@ -492,12 +469,12 @@ runTests n = do
492 test (svdProp6b) 469 test (svdProp6b)
493 test (svdProp7 . rM) 470 test (svdProp7 . rM)
494 test (svdProp7 . cM) 471 test (svdProp7 . cM)
495 putStrLn "------ svdCd" 472-- putStrLn "------ svdCd"
496#ifdef NOZGESDD 473#ifdef NOZGESDD
497 putStrLn "Omitted" 474-- putStrLn "Omitted"
498#else 475#else
499 test (svdProp1a svdCd) 476-- test (svdProp1a svdCd)
500 test (svdProp1b svdCd) 477-- test (svdProp1b svdCd)
501#endif 478#endif
502 putStrLn "------ eig" 479 putStrLn "------ eig"
503 test (eigSHProp . rHer) 480 test (eigSHProp . rHer)
@@ -515,10 +492,10 @@ runTests n = do
515 test (qrProp . rM) 492 test (qrProp . rM)
516 test (qrProp . cM) 493 test (qrProp . cM)
517 test (rqProp . rM) 494 test (rqProp . rM)
518 test (rqProp . cM) 495-- test (rqProp . cM)
519 test (rqProp1 . cM) 496 test (rqProp1 . cM)
520 test (rqProp2 . cM) 497 test (rqProp2 . cM)
521 test (rqProp3 . cM) 498-- test (rqProp3 . cM)
522 putStrLn "------ hess" 499 putStrLn "------ hess"
523 test (hessProp . rSq) 500 test (hessProp . rSq)
524 test (hessProp . cSq) 501 test (hessProp . cSq)
@@ -539,12 +516,12 @@ runTests n = do
539 test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM)) 516 test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM))
540 test (\u -> cos u * tan u |~| sin (u::RM)) 517 test (\u -> cos u * tan u |~| sin (u::RM))
541 test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary 518 test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary
542 putStrLn "------ vector operations - Float" 519-- putStrLn "------ vector operations - Float"
543 test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM)) 520-- test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM))
544 test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary 521-- test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary
545 test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM)) 522-- test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM))
546 test (\u -> cos u * tan u |~~| sin (u::FM)) 523-- test (\u -> cos u * tan u |~~| sin (u::FM))
547 test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary 524-- test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary
548 putStrLn "------ read . show" 525 putStrLn "------ read . show"
549 test (\m -> (m::RM) == read (show m)) 526 test (\m -> (m::RM) == read (show m))
550 test (\m -> (m::CM) == read (show m)) 527 test (\m -> (m::CM) == read (show m))
@@ -562,8 +539,8 @@ runTests n = do
562 , utest "expm1" (expmTest1) 539 , utest "expm1" (expmTest1)
563 , utest "expm2" (expmTest2) 540 , utest "expm2" (expmTest2)
564 , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM) 541 , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)
565 , utest "arith2" $ ((scalar (1+i) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*i-51) :: CM) 542 , utest "arith2" $ ((scalar (1+iC) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*iC-51) :: CM)
566 , utest "arith3" $ exp (scalar i * ones(10,10)*pi) + 1 |~| 0 543 , utest "arith3" $ exp (scalar iC * ones(10,10)*pi) + 1 |~| 0
567 , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3] 544 , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3]
568-- , utest "gamma" (gamma 5 == 24.0) 545-- , utest "gamma" (gamma 5 == 24.0)
569-- , besselTest 546-- , besselTest
@@ -571,10 +548,10 @@ runTests n = do
571 , utest "randomGaussian" randomTestGaussian 548 , utest "randomGaussian" randomTestGaussian
572 , utest "randomUniform" randomTestUniform 549 , utest "randomUniform" randomTestUniform
573 , utest "buildVector/Matrix" $ 550 , utest "buildVector/Matrix" $
574 complex (10 |> [0::Double ..]) == buildVector 10 fromIntegral 551 complex (10 |> [0::Double ..]) == build 10 id
575 && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0) 552 && ident 5 == build (5,5) (\r c -> if r==c then 1::Double else 0)
576 , utest "rank" $ rank ((2><3)[1,0,0,1,5*eps,0]) == 1 553 , utest "rank" $ rank ((2><3)[1,0,0,1,5*peps,0::Double]) == 1
577 && rank ((2><3)[1,0,0,1,7*eps,0]) == 2 554 && rank ((2><3)[1,0,0,1,7*peps,0::Double]) == 2
578 , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM) 555 , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)
579 , mbCholTest 556 , mbCholTest
580 , utest "offset" offsetTest 557 , utest "offset" offsetTest
@@ -588,7 +565,6 @@ runTests n = do
588 , conformTest 565 , conformTest
589 , accumTest 566 , accumTest
590 , convolutionTest 567 , convolutionTest
591 , kroneckerTest
592 , sparseTest 568 , sparseTest
593 , staticTest 569 , staticTest
594 ] 570 ]
@@ -597,12 +573,12 @@ runTests n = do
597 573
598 574
599-- single precision approximate equality 575-- single precision approximate equality
600infixl 4 |~~| 576-- infixl 4 |~~|
601a |~~| b = a :~6~: b 577-- a |~~| b = a :~6~: b
602 578
603makeUnitary v | realPart n > 1 = v / scalar n 579makeUnitary v | realPart n > 1 = v / scalar n
604 | otherwise = v 580 | otherwise = v
605 where n = sqrt (v <.> v) 581 where n = sqrt (v `dot` v)
606 582
607-- -- | Some additional tests on big matrices. They take a few minutes. 583-- -- | Some additional tests on big matrices. They take a few minutes.
608-- runBigTests :: IO () 584-- runBigTests :: IO ()
@@ -668,9 +644,9 @@ manyvec5 xs = sumElements $ fromRows $ map (\x -> vec3 x (x**2) (x**3)) xs
668 644
669 645
670manyvec2 xs = sum $ map (\x -> sqrt(x^2 + (x**2)^2 +(x**3)^2)) xs 646manyvec2 xs = sum $ map (\x -> sqrt(x^2 + (x**2)^2 +(x**3)^2)) xs
671manyvec3 xs = sum $ map (pnorm PNorm2 . (\x -> fromList [x,x**2,x**3])) xs 647manyvec3 xs = sum $ map (norm_2 . (\x -> fromList [x,x**2,x**3])) xs
672 648
673manyvec4 xs = sum $ map (pnorm PNorm2 . (\x -> vec3 x (x**2) (x**3))) xs 649manyvec4 xs = sum $ map (norm_2 . (\x -> vec3 x (x**2) (x**3))) xs
674 650
675vec3 :: Double -> Double -> Double -> Vector Double 651vec3 :: Double -> Double -> Double -> Vector Double
676vec3 a b c = runSTVector $ do 652vec3 a b c = runSTVector $ do
@@ -695,11 +671,11 @@ mkVecBench = do
695 671
696subBench = do 672subBench = do
697 putStrLn "" 673 putStrLn ""
698 let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v)) 674 let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (size v -1) v))
699 time "0.1M subVector " (g (konst 1 (1+10^5) :: Vector Double) @> 0) 675 time "0.1M subVector " (g (konst 1 (1+10^5) :: Vector Double) ! 0)
700 let f = foldl1' (.) (replicate (10^5) (fromRows.toRows)) 676 let f = foldl1' (.) (replicate (10^5) (fromRows.toRows))
701 time "subVector-join 3" (f (ident 3 :: Matrix Double) @@>(0,0)) 677 time "subVector-join 3" (f (ident 3 :: Matrix Double) `atIndex` (0,0))
702 time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0)) 678 time "subVector-join 10" (f (ident 10 :: Matrix Double) `atIndex` (0,0))
703 679
704-------------------------------- 680--------------------------------
705 681
@@ -724,7 +700,7 @@ multBench = do
724 700
725eigBench = do 701eigBench = do
726 let m = reshape 1000 (randomVector 777 Uniform (1000*1000)) 702 let m = reshape 1000 (randomVector 777 Uniform (1000*1000))
727 s = m + trans m 703 s = m + tr m
728 m `seq` s `seq` putStrLn "" 704 m `seq` s `seq` putStrLn ""
729 time "eigenvalues symmetric 1000x1000" (eigenvaluesSH' m) 705 time "eigenvalues symmetric 1000x1000" (eigenvaluesSH' m)
730 time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m) 706 time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m)
@@ -736,7 +712,7 @@ eigBench = do
736svdBench = do 712svdBench = do
737 let a = reshape 500 (randomVector 777 Uniform (3000*500)) 713 let a = reshape 500 (randomVector 777 Uniform (3000*500))
738 b = reshape 1000 (randomVector 777 Uniform (1000*1000)) 714 b = reshape 1000 (randomVector 777 Uniform (1000*1000))
739 fv (_,_,v) = v@@>(0,0) 715 fv (_,_,v) = v `atIndex` (0,0)
740 a `seq` b `seq` putStrLn "" 716 a `seq` b `seq` putStrLn ""
741 time "singular values 3000x500" (singularValues a) 717 time "singular values 3000x500" (singularValues a)
742 time "thin svd 3000x500" (fv $ thinSVD a) 718 time "thin svd 3000x500" (fv $ thinSVD a)
@@ -748,7 +724,7 @@ svdBench = do
748 724
749solveBenchN n = do 725solveBenchN n = do
750 let x = uniformSample 777 (2*n) (replicate n (-1,1)) 726 let x = uniformSample 777 (2*n) (replicate n (-1,1))
751 a = trans x <> x 727 a = tr x <> x
752 b = asColumn $ randomVector 666 Uniform n 728 b = asColumn $ randomVector 666 Uniform n
753 a `seq` b `seq` putStrLn "" 729 a `seq` b `seq` putStrLn ""
754 time ("svd solve " ++ show n) (linearSolveSVD a b) 730 time ("svd solve " ++ show n) (linearSolveSVD a b)
@@ -765,7 +741,7 @@ solveBench = do
765 741
766cholBenchN n = do 742cholBenchN n = do
767 let x = uniformSample 777 (2*n) (replicate n (-1,1)) 743 let x = uniformSample 777 (2*n) (replicate n (-1,1))
768 a = trans x <> x 744 a = tr x <> x
769 a `seq` putStr "" 745 a `seq` putStr ""
770 time ("chol " ++ show n) (chol a) 746 time ("chol " ++ show n) (chol a)
771 747
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
index 53fc4d2..904ae05 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -1,5 +1,4 @@
1{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-} 1{-# LANGUAGE FlexibleContexts, UndecidableInstances, FlexibleInstances #-}
2{-# OPTIONS_GHC -fno-warn-unused-imports #-}
3----------------------------------------------------------------------------- 2-----------------------------------------------------------------------------
4{- | 3{- |
5Module : Numeric.LinearAlgebra.Tests.Instances 4Module : Numeric.LinearAlgebra.Tests.Instances
@@ -26,15 +25,11 @@ module Numeric.LinearAlgebra.Tests.Instances(
26 25
27import System.Random 26import System.Random
28 27
29import Numeric.LinearAlgebra 28import Numeric.LinearAlgebra.HMatrix hiding (vector)
30import Numeric.LinearAlgebra.Devel
31import Numeric.Container
32import Control.Monad(replicateM) 29import Control.Monad(replicateM)
33import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector 30import Test.QuickCheck(Arbitrary,arbitrary,choose,vector,sized,shrink)
34 ,sized,classify,Testable,Property 31
35 ,quickCheckWith,maxSize,stdArgs,shrink)
36 32
37#if MIN_VERSION_QuickCheck(2,0,0)
38shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]] 33shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]]
39shrinkListElementwise [] = [] 34shrinkListElementwise [] = []
40shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ] 35shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ]
@@ -42,25 +37,6 @@ shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ]
42 37
43shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)] 38shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)]
44shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ] 39shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ]
45#endif
46
47#if MIN_VERSION_QuickCheck(2,1,1)
48#else
49instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where
50 arbitrary = do
51 re <- arbitrary
52 im <- arbitrary
53 return (re :+ im)
54
55#if MIN_VERSION_QuickCheck(2,0,0)
56 shrink (re :+ im) =
57 [ u :+ v | (u,v) <- shrinkPair (re,im) ]
58#else
59 -- this has been moved to the 'Coarbitrary' class in QuickCheck 2
60 coarbitrary = undefined
61#endif
62
63#endif
64 40
65chooseDim = sized $ \m -> choose (1,max 1 m) 41chooseDim = sized $ \m -> choose (1,max 1 m)
66 42
@@ -68,15 +44,9 @@ instance (Field a, Arbitrary a) => Arbitrary (Vector a) where
68 arbitrary = do m <- chooseDim 44 arbitrary = do m <- chooseDim
69 l <- vector m 45 l <- vector m
70 return $ fromList l 46 return $ fromList l
71
72#if MIN_VERSION_QuickCheck(2,0,0)
73 -- shrink any one of the components 47 -- shrink any one of the components
74 shrink = map fromList . shrinkListElementwise . toList 48 shrink = map fromList . shrinkListElementwise . toList
75 49
76#else
77 coarbitrary = undefined
78#endif
79
80instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where 50instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where
81 arbitrary = do 51 arbitrary = do
82 m <- chooseDim 52 m <- chooseDim
@@ -84,16 +54,11 @@ instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where
84 l <- vector (m*n) 54 l <- vector (m*n)
85 return $ (m><n) l 55 return $ (m><n) l
86 56
87#if MIN_VERSION_QuickCheck(2,0,0)
88 -- shrink any one of the components 57 -- shrink any one of the components
89 shrink a = map (rows a >< cols a) 58 shrink a = map (rows a >< cols a)
90 . shrinkListElementwise 59 . shrinkListElementwise
91 . concat . toLists 60 . concat . toLists
92 $ a 61 $ a
93#else
94 coarbitrary = undefined
95#endif
96
97 62
98-- a square matrix 63-- a square matrix
99newtype (Sq a) = Sq (Matrix a) deriving Show 64newtype (Sq a) = Sq (Matrix a) deriving Show
@@ -103,11 +68,7 @@ instance (Element a, Arbitrary a) => Arbitrary (Sq a) where
103 l <- vector (n*n) 68 l <- vector (n*n)
104 return $ Sq $ (n><n) l 69 return $ Sq $ (n><n) l
105 70
106#if MIN_VERSION_QuickCheck(2,0,0)
107 shrink (Sq a) = [ Sq b | b <- shrink a ] 71 shrink (Sq a) = [ Sq b | b <- shrink a ]
108#else
109 coarbitrary = undefined
110#endif
111 72
112 73
113-- a unitary matrix 74-- a unitary matrix
@@ -118,11 +79,6 @@ instance (Field a, Arbitrary a) => Arbitrary (Rot a) where
118 let (q,_) = qr m 79 let (q,_) = qr m
119 return (Rot q) 80 return (Rot q)
120 81
121#if MIN_VERSION_QuickCheck(2,0,0)
122#else
123 coarbitrary = undefined
124#endif
125
126 82
127-- a complex hermitian or real symmetric matrix 83-- a complex hermitian or real symmetric matrix
128newtype (Her a) = Her (Matrix a) deriving Show 84newtype (Her a) = Her (Matrix a) deriving Show
@@ -130,12 +86,8 @@ instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where
130 arbitrary = do 86 arbitrary = do
131 Sq m <- arbitrary 87 Sq m <- arbitrary
132 let m' = m/2 88 let m' = m/2
133 return $ Her (m' + ctrans m') 89 return $ Her (m' + tr m')
134 90
135#if MIN_VERSION_QuickCheck(2,0,0)
136#else
137 coarbitrary = undefined
138#endif
139 91
140class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a 92class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a
141instance ArbitraryField Double 93instance ArbitraryField Double
@@ -144,7 +96,7 @@ instance ArbitraryField (Complex Double)
144 96
145-- a well-conditioned general matrix (the singular values are between 1 and 100) 97-- a well-conditioned general matrix (the singular values are between 1 and 100)
146newtype (WC a) = WC (Matrix a) deriving Show 98newtype (WC a) = WC (Matrix a) deriving Show
147instance (ArbitraryField a) => Arbitrary (WC a) where 99instance (Numeric a, ArbitraryField a) => Arbitrary (WC a) where
148 arbitrary = do 100 arbitrary = do
149 m <- arbitrary 101 m <- arbitrary
150 let (u,_,v) = svd m 102 let (u,_,v) = svd m
@@ -153,34 +105,24 @@ instance (ArbitraryField a) => Arbitrary (WC a) where
153 n = min r c 105 n = min r c
154 sv' <- replicateM n (choose (1,100)) 106 sv' <- replicateM n (choose (1,100))
155 let s = diagRect 0 (fromList sv') r c 107 let s = diagRect 0 (fromList sv') r c
156 return $ WC (u `mXm` real s `mXm` trans v) 108 return $ WC (u <> real s <> tr v)
157
158#if MIN_VERSION_QuickCheck(2,0,0)
159#else
160 coarbitrary = undefined
161#endif
162 109
163 110
164-- a well-conditioned square matrix (the singular values are between 1 and 100) 111-- a well-conditioned square matrix (the singular values are between 1 and 100)
165newtype (SqWC a) = SqWC (Matrix a) deriving Show 112newtype (SqWC a) = SqWC (Matrix a) deriving Show
166instance (ArbitraryField a) => Arbitrary (SqWC a) where 113instance (ArbitraryField a, Numeric a) => Arbitrary (SqWC a) where
167 arbitrary = do 114 arbitrary = do
168 Sq m <- arbitrary 115 Sq m <- arbitrary
169 let (u,_,v) = svd m 116 let (u,_,v) = svd m
170 n = rows m 117 n = rows m
171 sv' <- replicateM n (choose (1,100)) 118 sv' <- replicateM n (choose (1,100))
172 let s = diag (fromList sv') 119 let s = diag (fromList sv')
173 return $ SqWC (u `mXm` real s `mXm` trans v) 120 return $ SqWC (u <> real s <> tr v)
174
175#if MIN_VERSION_QuickCheck(2,0,0)
176#else
177 coarbitrary = undefined
178#endif
179 121
180 122
181-- a positive definite square matrix (the eigenvalues are between 0 and 100) 123-- a positive definite square matrix (the eigenvalues are between 0 and 100)
182newtype (PosDef a) = PosDef (Matrix a) deriving Show 124newtype (PosDef a) = PosDef (Matrix a) deriving Show
183instance (ArbitraryField a, Num (Vector a)) 125instance (Numeric a, ArbitraryField a, Num (Vector a))
184 => Arbitrary (PosDef a) where 126 => Arbitrary (PosDef a) where
185 arbitrary = do 127 arbitrary = do
186 Her m <- arbitrary 128 Her m <- arbitrary
@@ -188,13 +130,8 @@ instance (ArbitraryField a, Num (Vector a))
188 n = rows m 130 n = rows m
189 l <- replicateM n (choose (0,100)) 131 l <- replicateM n (choose (0,100))
190 let s = diag (fromList l) 132 let s = diag (fromList l)
191 p = v `mXm` real s `mXm` ctrans v 133 p = v <> real s <> tr v
192 return $ PosDef (0.5 * p + 0.5 * ctrans p) 134 return $ PosDef (0.5 * p + 0.5 * tr p)
193
194#if MIN_VERSION_QuickCheck(2,0,0)
195#else
196 coarbitrary = undefined
197#endif
198 135
199 136
200-- a pair of matrices that can be multiplied 137-- a pair of matrices that can be multiplied
@@ -208,11 +145,7 @@ instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where
208 lb <- vector (k*m) 145 lb <- vector (k*m)
209 return $ Consistent ((n><k) la, (k><m) lb) 146 return $ Consistent ((n><k) la, (k><m) lb)
210 147
211#if MIN_VERSION_QuickCheck(2,0,0)
212 shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ] 148 shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ]
213#else
214 coarbitrary = undefined
215#endif
216 149
217 150
218 151
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
index a5c37f4..207a303 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -1,5 +1,4 @@
1{-# LANGUAGE CPP, FlexibleContexts #-} 1{-# LANGUAGE FlexibleContexts #-}
2{-# OPTIONS_GHC -fno-warn-unused-imports #-}
3{-# LANGUAGE TypeFamilies #-} 2{-# LANGUAGE TypeFamilies #-}
4 3
5----------------------------------------------------------------------------- 4-----------------------------------------------------------------------------
@@ -29,7 +28,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
29 pinvProp, 28 pinvProp,
30 detProp, 29 detProp,
31 nullspaceProp, 30 nullspaceProp,
32 bugProp, 31-- bugProp,
33 svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, 32 svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,
34 svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, 33 svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
35 eigProp, eigSHProp, eigProp2, eigSHProp2, 34 eigProp, eigSHProp, eigProp2, eigSHProp2,
@@ -43,20 +42,15 @@ module Numeric.LinearAlgebra.Tests.Properties (
43 linearSolveProp, linearSolveProp2 42 linearSolveProp, linearSolveProp2
44) where 43) where
45 44
46import Numeric.Container 45import Numeric.LinearAlgebra.HMatrix hiding (Testable,unitary)
47import Numeric.LinearAlgebra --hiding (real,complex) 46import Test.QuickCheck
48import Numeric.LinearAlgebra.LAPACK
49import Debug.Trace
50import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
51 ,sized,classify,Testable,Property
52 ,quickCheckWith,maxSize,stdArgs,shrink)
53 47
54trivial :: Testable a => Bool -> a -> Property 48trivial :: Testable a => Bool -> a -> Property
55trivial = (`classify` "trivial") 49trivial = (`classify` "trivial")
56 50
57-- relative error 51-- relative error
58dist :: (Normed c t, Num (c t)) => c t -> c t -> Double 52dist :: (Num a, Normed a) => a -> a -> Double
59dist = relativeError Infinity 53dist = relativeError norm_Inf
60 54
61infixl 4 |~| 55infixl 4 |~|
62a |~| b = a :~10~: b 56a |~| b = a :~10~: b
@@ -73,11 +67,11 @@ a :~n~: b = dist a b < 10^^(-n)
73square m = rows m == cols m 67square m = rows m == cols m
74 68
75-- orthonormal columns 69-- orthonormal columns
76orthonormal m = ctrans m <> m |~| ident (cols m) 70orthonormal m = tr m <> m |~| ident (cols m)
77 71
78unitary m = square m && orthonormal m 72unitary m = square m && orthonormal m
79 73
80hermitian m = square m && m |~| ctrans m 74hermitian m = square m && m |~| tr m
81 75
82wellCond m = rcond m > 1/100 76wellCond m = rcond m > 1/100
83 77
@@ -85,12 +79,12 @@ positiveDefinite m = minimum (toList e) > 0
85 where (e,_v) = eigSH m 79 where (e,_v) = eigSH m
86 80
87upperTriang m = rows m == 1 || down == z 81upperTriang m = rows m == 1 || down == z
88 where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m)) 82 where down = fromList $ concat $ zipWith drop [1..] (toLists (tr m))
89 z = konst 0 (dim down) 83 z = konst 0 (size down)
90 84
91upperHessenberg m = rows m < 3 || down == z 85upperHessenberg m = rows m < 3 || down == z
92 where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m)) 86 where down = fromList $ concat $ zipWith drop [2..] (toLists (tr m))
93 z = konst 0 (dim down) 87 z = konst 0 (size down)
94 88
95zeros (r,c) = reshape c (konst 0 (r*c)) 89zeros (r,c) = reshape c (konst 0 (r*c))
96 90
@@ -118,81 +112,94 @@ detProp m = s d1 |~| s d2
118 s x = fromList [x] 112 s x = fromList [x]
119 113
120nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c) 114nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c)
121 && orthonormal (fromColumns nl)) 115 && orthonormal n)
122 where nl = nullspacePrec 1 m 116 where n = nullspaceSVD (Left (1*peps)) m (rightSV m)
123 n = fromColumns nl 117 nl = toColumns n
124 r = rows m 118 r = rows m
125 c = cols m - rank m 119 c = cols m - rank m
126 120
127------------------------------------------------------------------ 121------------------------------------------------------------------
128 122{-
129-- testcase for nonempty fpu stack 123-- testcase for nonempty fpu stack
130-- uncommenting unitary' signature eliminates the problem 124-- uncommenting unitary' signature eliminates the problem
131bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v 125bugProp m = m |~| u <> real d <> tr v && unitary' u && unitary' v
132 where (u,d,v) = fullSVD m 126 where (u,d,v) = svd m
133 -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool 127 -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool
134 unitary' a = unitary a 128 unitary' a = unitary a
135 129-}
136------------------------------------------------------------------ 130------------------------------------------------------------------
137 131
138-- fullSVD 132-- fullSVD
139svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v 133svdProp1 m = m |~| u <> real d <> tr v && unitary u && unitary v
140 where (u,d,v) = fullSVD m 134 where
135 (u,s,v) = svd m
136 d = diagRect 0 s (rows m) (cols m)
141 137
142svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where 138svdProp1a svdfun m = m |~| u <> real d <> tr v && unitary u && unitary v
139 where
143 (u,s,v) = svdfun m 140 (u,s,v) = svdfun m
144 d = diagRect 0 s (rows m) (cols m) 141 d = diagRect 0 s (rows m) (cols m)
145 142
146svdProp1b svdfun m = unitary u && unitary v where 143svdProp1b svdfun m = unitary u && unitary v
144 where
147 (u,_,v) = svdfun m 145 (u,_,v) = svdfun m
148 146
149-- thinSVD 147-- thinSVD
150svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m) 148svdProp2 thinSVDfun m
151 where (u,s,v) = thinSVDfun m 149 = m |~| u <> diag (real s) <> tr v
150 && orthonormal u && orthonormal v
151 && size s == min (rows m) (cols m)
152 where
153 (u,s,v) = thinSVDfun m
152 154
153-- compactSVD 155-- compactSVD
154svdProp3 m = (m |~| u <> real (diag s) <> trans v 156svdProp3 m = (m |~| u <> real (diag s) <> tr v
155 && orthonormal u && orthonormal v) 157 && orthonormal u && orthonormal v)
156 where (u,s,v) = compactSVD m 158 where
159 (u,s,v) = compactSVD m
157 160
158svdProp4 m' = m |~| u <> real (diag s) <> trans v 161svdProp4 m' = m |~| u <> real (diag s) <> tr v
159 && orthonormal u && orthonormal v 162 && orthonormal u && orthonormal v
160 && (dim s == r || r == 0 && dim s == 1) 163 && (size s == r || r == 0 && size s == 1)
161 where (u,s,v) = compactSVD m 164 where
162 m = fromBlocks [[m'],[m']] 165 (u,s,v) = compactSVD m
163 r = rank m' 166 m = fromBlocks [[m'],[m']]
164 167 r = rank m'
165svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where 168
166 s1 = svR m 169svdProp5a m = all (s1|~|) [s3,s5] where
167 s2 = svRd m 170 s1 = singularValues (m :: Matrix Double)
168 (_,s3,_) = svdR m 171-- s2 = svRd m
169 (_,s4,_) = svdRd m 172 (_,s3,_) = svd m
170 (_,s5,_) = thinSVDR m 173-- (_,s4,_) = svdRd m
171 (_,s6,_) = thinSVDRd m 174 (_,s5,_) = thinSVD m
172 175-- (_,s6,_) = thinSVDRd m
173svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where 176
174 s1 = svC m 177svdProp5b m = all (s1|~|) [s3,s5] where
175 s2 = svCd m 178 s1 = singularValues (m :: Matrix (Complex Double))
176 (_,s3,_) = svdC m 179-- s2 = svCd m
177 (_,s4,_) = svdCd m 180 (_,s3,_) = svd m
178 (_,s5,_) = thinSVDC m 181-- (_,s4,_) = svdCd m
179 (_,s6,_) = thinSVDCd m 182 (_,s5,_) = thinSVD m
183-- (_,s6,_) = thinSVDCd m
180 184
181svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' 185svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
182 where (u,s,v) = svdR m 186 where
183 (s',v') = rightSVR m 187 (u,s,v) = svd (m :: Matrix Double)
184 (u',s'') = leftSVR m 188 (s',v') = rightSV m
189 (u',s'') = leftSV m
185 190
186svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' 191svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
187 where (u,s,v) = svdC m 192 where
188 (s',v') = rightSVC m 193 (u,s,v) = svd (m :: Matrix (Complex Double))
189 (u',s'') = leftSVC m 194 (s',v') = rightSV m
195 (u',s'') = leftSV m
190 196
191svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s''' 197svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''
192 where (u,s,v) = svd m 198 where
193 (s',v') = rightSV m 199 (u,s,v) = svd m
194 (u',_s'') = leftSV m 200 (s',v') = rightSV m
195 s''' = singularValues m 201 (u',_s'') = leftSV m
202 s''' = singularValues m
196 203
197------------------------------------------------------------------ 204------------------------------------------------------------------
198 205
@@ -201,7 +208,7 @@ eigProp m = complex m <> v |~| v <> diag s
201 208
202eigSHProp m = m <> v |~| v <> real (diag s) 209eigSHProp m = m <> v |~| v <> real (diag s)
203 && unitary v 210 && unitary v
204 && m |~| v <> real (diag s) <> ctrans v 211 && m |~| v <> real (diag s) <> tr v
205 where (s, v) = eigSH m 212 where (s, v) = eigSH m
206 213
207eigProp2 m = fst (eig m) |~| eigenvalues m 214eigProp2 m = fst (eig m) |~| eigenvalues m
@@ -226,19 +233,19 @@ rqProp3 m = upperTriang' r
226 where (r,_q) = rq m 233 where (r,_q) = rq m
227 234
228upperTriang' r = upptr (rows r) (cols r) * r |~| r 235upperTriang' r = upptr (rows r) (cols r) * r |~| r
229 where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1 236 where upptr f c = build (f,c) $ \r' c' -> if r'-t > c' then 0 else 1
230 where t = f-c 237 where t = fromIntegral (f-c)
231 238
232hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h 239hessProp m = m |~| p <> h <> tr p && unitary p && upperHessenberg h
233 where (p,h) = hess m 240 where (p,h) = hess m
234 241
235schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s 242schurProp1 m = m |~| u <> s <> tr u && unitary u && upperTriang s
236 where (u,s) = schur m 243 where (u,s) = schur m
237 244
238schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme 245schurProp2 m = m |~| u <> s <> tr u && unitary u && upperHessenberg s -- fixme
239 where (u,s) = schur m 246 where (u,s) = schur m
240 247
241cholProp m = m |~| ctrans c <> c && upperTriang c 248cholProp m = m |~| tr c <> c && upperTriang c
242 where c = chol m 249 where c = chol m
243 250
244exactProp m = chol m == chol (m+0) 251exactProp m = chol m == chol (m+0)
@@ -252,7 +259,7 @@ mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ]
252 259
253multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) 260multProp1 p (a,b) = (a <> b) :~p~: (mulH a b)
254 261
255multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a) 262multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a)
256 263
257linearSolveProp f m = f m m |~| ident (rows m) 264linearSolveProp f m = f m m |~| ident (rows m)
258 265
@@ -261,5 +268,5 @@ linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
261 b = a <> x 268 b = a <> x
262 wc = rank a == q 269 wc = rank a == q
263 270
264subProp m = m == (trans . fromColumns . toRows) m 271subProp m = m == (conj . tr . fromColumns . toRows) m
265 272