From 500f5fca244dadab494655aa73d7183df1c87c50 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Sat, 23 Feb 2008 19:35:51 +0000 Subject: working on tests --- lib/Numeric/LinearAlgebra/Algorithms.hs | 2 +- lib/Numeric/LinearAlgebra/Tests.hs | 61 +++++++++++------------ lib/Numeric/LinearAlgebra/Tests/Instances.hs | 63 ++++++++++++++++++++++-- lib/Numeric/LinearAlgebra/Tests/Properties.hs | 70 ++++++++++++++++++++++++++- 4 files changed, 157 insertions(+), 39 deletions(-) (limited to 'lib') diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index 79cc64d..1de018c 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs @@ -142,7 +142,7 @@ eigSH m | m `equal` ctrans m = eigSH' m -- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf. -- --- If @c = chol m@ then @m == c \<> ctrans c@. +-- If @c = chol m@ then @m == ctrans c \<> c@. chol :: Field t => Matrix t -> Matrix t chol m | m `equal` ctrans m = cholSH m | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs index 534eb04..7ecf812 100644 --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ b/lib/Numeric/LinearAlgebra/Tests.hs @@ -15,8 +15,7 @@ Some tests. module Numeric.LinearAlgebra.Tests( module Numeric.LinearAlgebra.Tests.Instances, module Numeric.LinearAlgebra.Tests.Properties, - qCheck,RM,CM, rM,cM, rHer,cHer,rRot,cRot,rSq,cSq, - runTests + qCheck, runTests --, runBigTests ) where @@ -24,43 +23,41 @@ import Numeric.LinearAlgebra import Numeric.LinearAlgebra.Tests.Instances import Numeric.LinearAlgebra.Tests.Properties import Test.QuickCheck -import Debug.Trace +import Numeric.GSL(setErrorHandlerOff) qCheck n = check defaultConfig {configSize = const n} -debug x = trace (show x) x - -type RM = Matrix Double -type CM = Matrix (Complex Double) - -rM m = m :: RM -cM m = m :: CM - -rHer (Her m) = m :: RM -cHer (Her m) = m :: CM - -rRot (Rot m) = m :: RM -cRot (Rot m) = m :: CM - -rSq (Sq m) = m :: RM -cSq (Sq m) = m :: CM - -rWC (WC m) = m :: RM -cWC (WC m) = m :: CM - -- | It runs all the tests. runTests :: Int -- ^ maximum dimension -> IO () runTests n = do - qCheck n (hermitian . rHer) - qCheck n (hermitian . cHer) - qCheck n (unitary . rRot) - qCheck n (unitary . cRot) - qCheck n (wellCond . rWC) - qCheck n (wellCond . cWC) - -------------------------------- - qCheck n (luTest . rM) - qCheck n (luTest . cM) + setErrorHandlerOff + let test p = qCheck n p + test (luProp . rM) + test (luProp . cM) + test (invProp . rSqWC) + test (invProp . cSqWC) + test (pinvProp . rM) + test (pinvProp . cM) + test (detProp . cSqWC) + test (svdProp1 . rM) + test (svdProp1 . cM) + test (svdProp2 . rM) + test (svdProp2 . cM) + test (eigSHProp . rHer) + test (eigSHProp . cHer) + test (eigProp . rSq) + test (eigProp . cSq) + test (nullspaceProp . rM) + test (nullspaceProp . cM) + test (qrProp . rM) + test (qrProp . cM) + test (hessProp . rSq) + test (hessProp . cSq) + test (schurProp2 . rSq) + test (schurProp1 . cSq) + test (cholProp . rPosDef) + test (cholProp . cPosDef) -- | Some additional tests on big matrices. They take a few minutes. runBigTests :: IO () diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs index 583143a..af486c8 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs @@ -14,10 +14,13 @@ Arbitrary instances for vectors, matrices. -} module Numeric.LinearAlgebra.Tests.Instances( - Sq(..), - Rot(..), - Her(..), - WC(..) + Sq(..), rSq,cSq, + Rot(..), rRot,cRot, + Her(..), rHer,cHer, + WC(..), rWC,cWC, + SqWC(..), rSqWC, cSqWC, + PosDef(..), rPosDef, cPosDef, + RM,CM, rM,cM ) where import Numeric.LinearAlgebra @@ -74,7 +77,7 @@ instance (Field a, Arbitrary a) => Arbitrary (Her a) where return $ Her (m' + ctrans m') coarbitrary = undefined --- a well-conditioned matrix (the singular values are between 1 and 100) +-- a well-conditioned general matrix (the singular values are between 1 and 100) newtype (WC a) = WC (Matrix a) deriving Show instance (Field a, Arbitrary a) => Arbitrary (WC a) where arbitrary = do @@ -87,3 +90,53 @@ instance (Field a, Arbitrary a) => Arbitrary (WC a) where let s = diagRect (fromList sv) r c return $ WC (u <> real s <> trans v) coarbitrary = undefined + +-- a well-conditioned square matrix (the singular values are between 1 and 100) +newtype (SqWC a) = SqWC (Matrix a) deriving Show +instance (Field a, Arbitrary a) => Arbitrary (SqWC a) where + arbitrary = do + Sq m <- arbitrary + let (u,_,v) = svd m + n = rows m + sv <- replicateM n (choose (1,100)) + let s = diag (fromList sv) + return $ SqWC (u <> real s <> trans v) + coarbitrary = undefined + +-- a positive definite square matrix (the eigenvalues are between 0 and 100) +newtype (PosDef a) = PosDef (Matrix a) deriving Show +instance (Field a, Arbitrary a) => Arbitrary (PosDef a) where + arbitrary = do + Her m <- arbitrary + let (_,v) = eigSH m + n = rows m + l <- replicateM n (choose (0,100)) + let s = diag (fromList l) + p = v <> real s <> ctrans v + return $ PosDef (0.5 .* p + 0.5 .* ctrans p) + coarbitrary = undefined + +type RM = Matrix Double +type CM = Matrix (Complex Double) + +rM m = m :: RM +cM m = m :: CM + +rHer (Her m) = m :: RM +cHer (Her m) = m :: CM + +rRot (Rot m) = m :: RM +cRot (Rot m) = m :: CM + +rSq (Sq m) = m :: RM +cSq (Sq m) = m :: CM + +rWC (WC m) = m :: RM +cWC (WC m) = m :: CM + +rSqWC (SqWC m) = m :: RM +cSqWC (SqWC m) = m :: CM + +rPosDef (PosDef m) = m :: RM +cPosDef (PosDef m) = m :: CM + diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs index 351615b..0317469 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs @@ -19,6 +19,7 @@ where import Numeric.LinearAlgebra import Numeric.LinearAlgebra.Tests.Instances(Sq(..),Her(..),Rot(..)) +import Test.QuickCheck -- relative error dist :: (Normed t, Num t) => t -> t -> Double @@ -53,7 +54,74 @@ degenerate m = rank m < min (rows m) (cols m) wellCond m = rcond m > 1/100 +positiveDefinite m = minimum (toList e) > 0 + where (e,v) = eigSH m + +upperTriang m = rows m == 1 || down == z + where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m)) + z = constant 0 (dim down) + +upperHessenberg m = rows m < 3 || down == z + where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m)) + z = constant 0 (dim down) + +zeros (r,c) = reshape c (constant 0 (r*c)) + +ones (r,c) = zeros (r,c) + 1 + ----------------------------------------------------- -luTest m = m |~| p <> l <> u && det p == s +luProp m = m |~| p <> l <> u && det p == s where (l,u,p,s) = lu m + +invProp m = m <> inv m |~| ident (rows m) + +pinvProp m = m <> p <> m |~| m + && p <> m <> p |~| p + && hermitian (m<>p) + && hermitian (p<>m) + where p = pinv m + +detProp m = s d1 |~| s d2 + where d1 = det m + d2 = det' m * det q + det' m = product $ toList $ takeDiag r + (q,r) = qr m + s x = fromList [x] + +nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c)) + where nl = nullspacePrec 1 m + n = fromColumns nl + r = rows m + c = cols m - rank m + +svdProp1 m = u <> real d <> trans v |~| m + && unitary u && unitary v + where (u,d,v) = full svd m + +svdProp2 m = (m |~| 0) `trivial` ((m |~| 0) || u <> real (diag s) <> trans v |~| m) + where (u,s,v) = economy svd m + +eigProp m = complex m <> v |~| v <> diag s + where (s, v) = eig m + +eigSHProp m = m <> v |~| v <> real (diag s) + && unitary v + && m |~| v <> real (diag s) <> ctrans v + where (s, v) = eigSH m + +qrProp m = q <> r |~| m && unitary q && upperTriang r + where (q,r) = qr m + +hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h + where (p,h) = hess m + +schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s + where (u,s) = schur m + +schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme + where (u,s) = schur m + +cholProp m = m |~| ctrans c <> c && upperTriang c + where c = chol m + pos = positiveDefinite m -- cgit v1.2.3