From dcc03a4a764cb8683b80758af97fcbcc9aadba73 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 8 Jan 2015 16:15:29 +0100 Subject: wip on tests --- .../src/Numeric/LinearAlgebra/Tests/Instances.hs | 19 ++- .../src/Numeric/LinearAlgebra/Tests/Properties.hs | 147 +++++++++++---------- 2 files changed, 89 insertions(+), 77 deletions(-) (limited to 'packages/tests/src/Numeric/LinearAlgebra/Tests') diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs index 53fc4d2..e2c3840 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs @@ -26,9 +26,8 @@ module Numeric.LinearAlgebra.Tests.Instances( import System.Random -import Numeric.LinearAlgebra +import Numeric.LinearAlgebra.HMatrix hiding (vector) import Numeric.LinearAlgebra.Devel -import Numeric.Container import Control.Monad(replicateM) import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector ,sized,classify,Testable,Property @@ -130,7 +129,7 @@ instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where arbitrary = do Sq m <- arbitrary let m' = m/2 - return $ Her (m' + ctrans m') + return $ Her (m' + tr m') #if MIN_VERSION_QuickCheck(2,0,0) #else @@ -144,7 +143,7 @@ instance ArbitraryField (Complex Double) -- a well-conditioned general matrix (the singular values are between 1 and 100) newtype (WC a) = WC (Matrix a) deriving Show -instance (ArbitraryField a) => Arbitrary (WC a) where +instance (Numeric a, ArbitraryField a) => Arbitrary (WC a) where arbitrary = do m <- arbitrary let (u,_,v) = svd m @@ -153,7 +152,7 @@ instance (ArbitraryField a) => Arbitrary (WC a) where n = min r c sv' <- replicateM n (choose (1,100)) let s = diagRect 0 (fromList sv') r c - return $ WC (u `mXm` real s `mXm` trans v) + return $ WC (u <> real s <> tr v) #if MIN_VERSION_QuickCheck(2,0,0) #else @@ -163,14 +162,14 @@ instance (ArbitraryField a) => Arbitrary (WC a) where -- a well-conditioned square matrix (the singular values are between 1 and 100) newtype (SqWC a) = SqWC (Matrix a) deriving Show -instance (ArbitraryField a) => Arbitrary (SqWC a) where +instance (ArbitraryField a, Numeric 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 `mXm` real s `mXm` trans v) + return $ SqWC (u <> real s <> tr v) #if MIN_VERSION_QuickCheck(2,0,0) #else @@ -180,7 +179,7 @@ instance (ArbitraryField a) => Arbitrary (SqWC a) where -- a positive definite square matrix (the eigenvalues are between 0 and 100) newtype (PosDef a) = PosDef (Matrix a) deriving Show -instance (ArbitraryField a, Num (Vector a)) +instance (Numeric a, ArbitraryField a, Num (Vector a)) => Arbitrary (PosDef a) where arbitrary = do Her m <- arbitrary @@ -188,8 +187,8 @@ instance (ArbitraryField a, Num (Vector a)) n = rows m l <- replicateM n (choose (0,100)) let s = diag (fromList l) - p = v `mXm` real s `mXm` ctrans v - return $ PosDef (0.5 * p + 0.5 * ctrans p) + p = v <> real s <> tr v + return $ PosDef (0.5 * p + 0.5 * tr p) #if MIN_VERSION_QuickCheck(2,0,0) #else diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs index 9bdf897..d9645c3 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs @@ -1,5 +1,7 @@ {-# LANGUAGE CPP, FlexibleContexts #-} {-# OPTIONS_GHC -fno-warn-unused-imports #-} +{-# LANGUAGE GADTs #-} + ----------------------------------------------------------------------------- {- | Module : Numeric.LinearAlgebra.Tests.Properties @@ -27,7 +29,7 @@ module Numeric.LinearAlgebra.Tests.Properties ( pinvProp, detProp, nullspaceProp, - bugProp, +-- bugProp, svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, eigProp, eigSHProp, eigProp2, eigSHProp2, @@ -41,9 +43,7 @@ module Numeric.LinearAlgebra.Tests.Properties ( linearSolveProp, linearSolveProp2 ) where -import Numeric.Container -import Numeric.LinearAlgebra --hiding (real,complex) -import Numeric.LinearAlgebra.LAPACK +import Numeric.LinearAlgebra.HMatrix hiding (Testable)--hiding (real,complex) import Debug.Trace import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector ,sized,classify,Testable,Property @@ -53,8 +53,8 @@ trivial :: Testable a => Bool -> a -> Property trivial = (`classify` "trivial") -- relative error -dist :: (Normed c t, Num (c t)) => c t -> c t -> Double -dist = relativeError Infinity +dist :: (Num a, Normed a) => a -> a -> Double +dist = relativeError norm_Inf infixl 4 |~| a |~| b = a :~10~: b @@ -71,11 +71,11 @@ a :~n~: b = dist a b < 10^^(-n) square m = rows m == cols m -- orthonormal columns -orthonormal m = ctrans m <> m |~| ident (cols m) +orthonormal m = tr m <> m |~| ident (cols m) unitary m = square m && orthonormal m -hermitian m = square m && m |~| ctrans m +hermitian m = square m && m |~| tr m wellCond m = rcond m > 1/100 @@ -83,12 +83,12 @@ 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 = konst 0 (dim down) + where down = fromList $ concat $ zipWith drop [1..] (toLists (tr m)) + z = konst 0 (size down) upperHessenberg m = rows m < 3 || down == z - where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m)) - z = konst 0 (dim down) + where down = fromList $ concat $ zipWith drop [2..] (toLists (tr m)) + z = konst 0 (size down) zeros (r,c) = reshape c (konst 0 (r*c)) @@ -116,81 +116,94 @@ detProp m = s d1 |~| s d2 s x = fromList [x] nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c) - && orthonormal (fromColumns nl)) - where nl = nullspacePrec 1 m - n = fromColumns nl + && orthonormal n) + where n = nullspaceSVD (Left (1*peps)) m (rightSV m) + nl = toColumns n r = rows m c = cols m - rank m ------------------------------------------------------------------ - +{- -- testcase for nonempty fpu stack -- uncommenting unitary' signature eliminates the problem -bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v - where (u,d,v) = fullSVD m +bugProp m = m |~| u <> real d <> tr v && unitary' u && unitary' v + where (u,d,v) = svd m -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool unitary' a = unitary a - +-} ------------------------------------------------------------------ -- fullSVD -svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v - where (u,d,v) = fullSVD m +svdProp1 m = m |~| u <> real d <> tr v && unitary u && unitary v + where + (u,s,v) = svd m + d = diagRect 0 s (rows m) (cols m) -svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where +svdProp1a svdfun m = m |~| u <> real d <> tr v && unitary u && unitary v + where (u,s,v) = svdfun m d = diagRect 0 s (rows m) (cols m) -svdProp1b svdfun m = unitary u && unitary v where +svdProp1b svdfun m = unitary u && unitary v + where (u,_,v) = svdfun m -- thinSVD -svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m) - where (u,s,v) = thinSVDfun m +svdProp2 thinSVDfun m + = m |~| u <> diag (real s) <> tr v + && orthonormal u && orthonormal v + && size s == min (rows m) (cols m) + where + (u,s,v) = thinSVDfun m -- compactSVD -svdProp3 m = (m |~| u <> real (diag s) <> trans v +svdProp3 m = (m |~| u <> real (diag s) <> tr v && orthonormal u && orthonormal v) - where (u,s,v) = compactSVD m + where + (u,s,v) = compactSVD m -svdProp4 m' = m |~| u <> real (diag s) <> trans v +svdProp4 m' = m |~| u <> real (diag s) <> tr v && orthonormal u && orthonormal v - && (dim s == r || r == 0 && dim s == 1) - where (u,s,v) = compactSVD m - m = fromBlocks [[m'],[m']] - r = rank m' - -svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where - s1 = svR m - s2 = svRd m - (_,s3,_) = svdR m - (_,s4,_) = svdRd m - (_,s5,_) = thinSVDR m - (_,s6,_) = thinSVDRd m - -svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where - s1 = svC m - s2 = svCd m - (_,s3,_) = svdC m - (_,s4,_) = svdCd m - (_,s5,_) = thinSVDC m - (_,s6,_) = thinSVDCd m + && (size s == r || r == 0 && size s == 1) + where + (u,s,v) = compactSVD m + m = fromBlocks [[m'],[m']] + r = rank m' + +svdProp5a m = all (s1|~|) [s3,s5] where + s1 = singularValues (m :: Matrix Double) +-- s2 = svRd m + (_,s3,_) = svd m +-- (_,s4,_) = svdRd m + (_,s5,_) = thinSVD m +-- (_,s6,_) = thinSVDRd m + +svdProp5b m = all (s1|~|) [s3,s5] where + s1 = singularValues (m :: Matrix (Complex Double)) +-- s2 = svCd m + (_,s3,_) = svd m +-- (_,s4,_) = svdCd m + (_,s5,_) = thinSVD m +-- (_,s6,_) = thinSVDCd m svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' - where (u,s,v) = svdR m - (s',v') = rightSVR m - (u',s'') = leftSVR m + where + (u,s,v) = svd (m :: Matrix Double) + (s',v') = rightSV m + (u',s'') = leftSV m svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' - where (u,s,v) = svdC m - (s',v') = rightSVC m - (u',s'') = leftSVC m + where + (u,s,v) = svd (m :: Matrix (Complex Double)) + (s',v') = rightSV m + (u',s'') = leftSV m svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s''' - where (u,s,v) = svd m - (s',v') = rightSV m - (u',_s'') = leftSV m - s''' = singularValues m + where + (u,s,v) = svd m + (s',v') = rightSV m + (u',_s'') = leftSV m + s''' = singularValues m ------------------------------------------------------------------ @@ -199,7 +212,7 @@ eigProp m = complex m <> v |~| v <> diag s eigSHProp m = m <> v |~| v <> real (diag s) && unitary v - && m |~| v <> real (diag s) <> ctrans v + && m |~| v <> real (diag s) <> tr v where (s, v) = eigSH m eigProp2 m = fst (eig m) |~| eigenvalues m @@ -224,19 +237,19 @@ rqProp3 m = upperTriang' r where (r,_q) = rq m upperTriang' r = upptr (rows r) (cols r) * r |~| r - where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1 - where t = f-c + where upptr f c = build (f,c) $ \r' c' -> if r'-t > c' then 0 else 1 + where t = fromIntegral (f-c) -hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h +hessProp m = m |~| p <> h <> tr p && unitary p && upperHessenberg h where (p,h) = hess m -schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s +schurProp1 m = m |~| u <> s <> tr u && unitary u && upperTriang s where (u,s) = schur m -schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme +schurProp2 m = m |~| u <> s <> tr u && unitary u && upperHessenberg s -- fixme where (u,s) = schur m -cholProp m = m |~| ctrans c <> c && upperTriang c +cholProp m = m |~| tr c <> c && upperTriang c where c = chol m exactProp m = chol m == chol (m+0) @@ -250,7 +263,7 @@ mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ] multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) -multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a) +multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a) linearSolveProp f m = f m m |~| ident (rows m) @@ -259,5 +272,5 @@ linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b) b = a <> x wc = rank a == q -subProp m = m == (trans . fromColumns . toRows) m +subProp m = m == (tr . fromColumns . toRows) m -- cgit v1.2.3