summaryrefslogtreecommitdiff
path: root/packages/tests/src/Numeric/LinearAlgebra/Tests
diff options
context:
space:
mode:
Diffstat (limited to 'packages/tests/src/Numeric/LinearAlgebra/Tests')
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs19
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs147
2 files changed, 89 insertions, 77 deletions
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(
26 26
27import System.Random 27import System.Random
28 28
29import Numeric.LinearAlgebra 29import Numeric.LinearAlgebra.HMatrix hiding (vector)
30import Numeric.LinearAlgebra.Devel 30import Numeric.LinearAlgebra.Devel
31import Numeric.Container
32import Control.Monad(replicateM) 31import Control.Monad(replicateM)
33import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector 32import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
34 ,sized,classify,Testable,Property 33 ,sized,classify,Testable,Property
@@ -130,7 +129,7 @@ instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where
130 arbitrary = do 129 arbitrary = do
131 Sq m <- arbitrary 130 Sq m <- arbitrary
132 let m' = m/2 131 let m' = m/2
133 return $ Her (m' + ctrans m') 132 return $ Her (m' + tr m')
134 133
135#if MIN_VERSION_QuickCheck(2,0,0) 134#if MIN_VERSION_QuickCheck(2,0,0)
136#else 135#else
@@ -144,7 +143,7 @@ instance ArbitraryField (Complex Double)
144 143
145-- a well-conditioned general matrix (the singular values are between 1 and 100) 144-- a well-conditioned general matrix (the singular values are between 1 and 100)
146newtype (WC a) = WC (Matrix a) deriving Show 145newtype (WC a) = WC (Matrix a) deriving Show
147instance (ArbitraryField a) => Arbitrary (WC a) where 146instance (Numeric a, ArbitraryField a) => Arbitrary (WC a) where
148 arbitrary = do 147 arbitrary = do
149 m <- arbitrary 148 m <- arbitrary
150 let (u,_,v) = svd m 149 let (u,_,v) = svd m
@@ -153,7 +152,7 @@ instance (ArbitraryField a) => Arbitrary (WC a) where
153 n = min r c 152 n = min r c
154 sv' <- replicateM n (choose (1,100)) 153 sv' <- replicateM n (choose (1,100))
155 let s = diagRect 0 (fromList sv') r c 154 let s = diagRect 0 (fromList sv') r c
156 return $ WC (u `mXm` real s `mXm` trans v) 155 return $ WC (u <> real s <> tr v)
157 156
158#if MIN_VERSION_QuickCheck(2,0,0) 157#if MIN_VERSION_QuickCheck(2,0,0)
159#else 158#else
@@ -163,14 +162,14 @@ instance (ArbitraryField a) => Arbitrary (WC a) where
163 162
164-- a well-conditioned square matrix (the singular values are between 1 and 100) 163-- a well-conditioned square matrix (the singular values are between 1 and 100)
165newtype (SqWC a) = SqWC (Matrix a) deriving Show 164newtype (SqWC a) = SqWC (Matrix a) deriving Show
166instance (ArbitraryField a) => Arbitrary (SqWC a) where 165instance (ArbitraryField a, Numeric a) => Arbitrary (SqWC a) where
167 arbitrary = do 166 arbitrary = do
168 Sq m <- arbitrary 167 Sq m <- arbitrary
169 let (u,_,v) = svd m 168 let (u,_,v) = svd m
170 n = rows m 169 n = rows m
171 sv' <- replicateM n (choose (1,100)) 170 sv' <- replicateM n (choose (1,100))
172 let s = diag (fromList sv') 171 let s = diag (fromList sv')
173 return $ SqWC (u `mXm` real s `mXm` trans v) 172 return $ SqWC (u <> real s <> tr v)
174 173
175#if MIN_VERSION_QuickCheck(2,0,0) 174#if MIN_VERSION_QuickCheck(2,0,0)
176#else 175#else
@@ -180,7 +179,7 @@ instance (ArbitraryField a) => Arbitrary (SqWC a) where
180 179
181-- a positive definite square matrix (the eigenvalues are between 0 and 100) 180-- a positive definite square matrix (the eigenvalues are between 0 and 100)
182newtype (PosDef a) = PosDef (Matrix a) deriving Show 181newtype (PosDef a) = PosDef (Matrix a) deriving Show
183instance (ArbitraryField a, Num (Vector a)) 182instance (Numeric a, ArbitraryField a, Num (Vector a))
184 => Arbitrary (PosDef a) where 183 => Arbitrary (PosDef a) where
185 arbitrary = do 184 arbitrary = do
186 Her m <- arbitrary 185 Her m <- arbitrary
@@ -188,8 +187,8 @@ instance (ArbitraryField a, Num (Vector a))
188 n = rows m 187 n = rows m
189 l <- replicateM n (choose (0,100)) 188 l <- replicateM n (choose (0,100))
190 let s = diag (fromList l) 189 let s = diag (fromList l)
191 p = v `mXm` real s `mXm` ctrans v 190 p = v <> real s <> tr v
192 return $ PosDef (0.5 * p + 0.5 * ctrans p) 191 return $ PosDef (0.5 * p + 0.5 * tr p)
193 192
194#if MIN_VERSION_QuickCheck(2,0,0) 193#if MIN_VERSION_QuickCheck(2,0,0)
195#else 194#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 @@
1{-# LANGUAGE CPP, FlexibleContexts #-} 1{-# LANGUAGE CPP, FlexibleContexts #-}
2{-# OPTIONS_GHC -fno-warn-unused-imports #-} 2{-# OPTIONS_GHC -fno-warn-unused-imports #-}
3{-# LANGUAGE GADTs #-}
4
3----------------------------------------------------------------------------- 5-----------------------------------------------------------------------------
4{- | 6{- |
5Module : Numeric.LinearAlgebra.Tests.Properties 7Module : Numeric.LinearAlgebra.Tests.Properties
@@ -27,7 +29,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
27 pinvProp, 29 pinvProp,
28 detProp, 30 detProp,
29 nullspaceProp, 31 nullspaceProp,
30 bugProp, 32-- bugProp,
31 svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, 33 svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,
32 svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, 34 svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
33 eigProp, eigSHProp, eigProp2, eigSHProp2, 35 eigProp, eigSHProp, eigProp2, eigSHProp2,
@@ -41,9 +43,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
41 linearSolveProp, linearSolveProp2 43 linearSolveProp, linearSolveProp2
42) where 44) where
43 45
44import Numeric.Container 46import Numeric.LinearAlgebra.HMatrix hiding (Testable)--hiding (real,complex)
45import Numeric.LinearAlgebra --hiding (real,complex)
46import Numeric.LinearAlgebra.LAPACK
47import Debug.Trace 47import Debug.Trace
48import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector 48import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
49 ,sized,classify,Testable,Property 49 ,sized,classify,Testable,Property
@@ -53,8 +53,8 @@ trivial :: Testable a => Bool -> a -> Property
53trivial = (`classify` "trivial") 53trivial = (`classify` "trivial")
54 54
55-- relative error 55-- relative error
56dist :: (Normed c t, Num (c t)) => c t -> c t -> Double 56dist :: (Num a, Normed a) => a -> a -> Double
57dist = relativeError Infinity 57dist = relativeError norm_Inf
58 58
59infixl 4 |~| 59infixl 4 |~|
60a |~| b = a :~10~: b 60a |~| b = a :~10~: b
@@ -71,11 +71,11 @@ a :~n~: b = dist a b < 10^^(-n)
71square m = rows m == cols m 71square m = rows m == cols m
72 72
73-- orthonormal columns 73-- orthonormal columns
74orthonormal m = ctrans m <> m |~| ident (cols m) 74orthonormal m = tr m <> m |~| ident (cols m)
75 75
76unitary m = square m && orthonormal m 76unitary m = square m && orthonormal m
77 77
78hermitian m = square m && m |~| ctrans m 78hermitian m = square m && m |~| tr m
79 79
80wellCond m = rcond m > 1/100 80wellCond m = rcond m > 1/100
81 81
@@ -83,12 +83,12 @@ positiveDefinite m = minimum (toList e) > 0
83 where (e,_v) = eigSH m 83 where (e,_v) = eigSH m
84 84
85upperTriang m = rows m == 1 || down == z 85upperTriang m = rows m == 1 || down == z
86 where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m)) 86 where down = fromList $ concat $ zipWith drop [1..] (toLists (tr m))
87 z = konst 0 (dim down) 87 z = konst 0 (size down)
88 88
89upperHessenberg m = rows m < 3 || down == z 89upperHessenberg m = rows m < 3 || down == z
90 where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m)) 90 where down = fromList $ concat $ zipWith drop [2..] (toLists (tr m))
91 z = konst 0 (dim down) 91 z = konst 0 (size down)
92 92
93zeros (r,c) = reshape c (konst 0 (r*c)) 93zeros (r,c) = reshape c (konst 0 (r*c))
94 94
@@ -116,81 +116,94 @@ detProp m = s d1 |~| s d2
116 s x = fromList [x] 116 s x = fromList [x]
117 117
118nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c) 118nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c)
119 && orthonormal (fromColumns nl)) 119 && orthonormal n)
120 where nl = nullspacePrec 1 m 120 where n = nullspaceSVD (Left (1*peps)) m (rightSV m)
121 n = fromColumns nl 121 nl = toColumns n
122 r = rows m 122 r = rows m
123 c = cols m - rank m 123 c = cols m - rank m
124 124
125------------------------------------------------------------------ 125------------------------------------------------------------------
126 126{-
127-- testcase for nonempty fpu stack 127-- testcase for nonempty fpu stack
128-- uncommenting unitary' signature eliminates the problem 128-- uncommenting unitary' signature eliminates the problem
129bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v 129bugProp m = m |~| u <> real d <> tr v && unitary' u && unitary' v
130 where (u,d,v) = fullSVD m 130 where (u,d,v) = svd m
131 -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool 131 -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool
132 unitary' a = unitary a 132 unitary' a = unitary a
133 133-}
134------------------------------------------------------------------ 134------------------------------------------------------------------
135 135
136-- fullSVD 136-- fullSVD
137svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v 137svdProp1 m = m |~| u <> real d <> tr v && unitary u && unitary v
138 where (u,d,v) = fullSVD m 138 where
139 (u,s,v) = svd m
140 d = diagRect 0 s (rows m) (cols m)
139 141
140svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where 142svdProp1a svdfun m = m |~| u <> real d <> tr v && unitary u && unitary v
143 where
141 (u,s,v) = svdfun m 144 (u,s,v) = svdfun m
142 d = diagRect 0 s (rows m) (cols m) 145 d = diagRect 0 s (rows m) (cols m)
143 146
144svdProp1b svdfun m = unitary u && unitary v where 147svdProp1b svdfun m = unitary u && unitary v
148 where
145 (u,_,v) = svdfun m 149 (u,_,v) = svdfun m
146 150
147-- thinSVD 151-- thinSVD
148svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m) 152svdProp2 thinSVDfun m
149 where (u,s,v) = thinSVDfun m 153 = m |~| u <> diag (real s) <> tr v
154 && orthonormal u && orthonormal v
155 && size s == min (rows m) (cols m)
156 where
157 (u,s,v) = thinSVDfun m
150 158
151-- compactSVD 159-- compactSVD
152svdProp3 m = (m |~| u <> real (diag s) <> trans v 160svdProp3 m = (m |~| u <> real (diag s) <> tr v
153 && orthonormal u && orthonormal v) 161 && orthonormal u && orthonormal v)
154 where (u,s,v) = compactSVD m 162 where
163 (u,s,v) = compactSVD m
155 164
156svdProp4 m' = m |~| u <> real (diag s) <> trans v 165svdProp4 m' = m |~| u <> real (diag s) <> tr v
157 && orthonormal u && orthonormal v 166 && orthonormal u && orthonormal v
158 && (dim s == r || r == 0 && dim s == 1) 167 && (size s == r || r == 0 && size s == 1)
159 where (u,s,v) = compactSVD m 168 where
160 m = fromBlocks [[m'],[m']] 169 (u,s,v) = compactSVD m
161 r = rank m' 170 m = fromBlocks [[m'],[m']]
162 171 r = rank m'
163svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where 172
164 s1 = svR m 173svdProp5a m = all (s1|~|) [s3,s5] where
165 s2 = svRd m 174 s1 = singularValues (m :: Matrix Double)
166 (_,s3,_) = svdR m 175-- s2 = svRd m
167 (_,s4,_) = svdRd m 176 (_,s3,_) = svd m
168 (_,s5,_) = thinSVDR m 177-- (_,s4,_) = svdRd m
169 (_,s6,_) = thinSVDRd m 178 (_,s5,_) = thinSVD m
170 179-- (_,s6,_) = thinSVDRd m
171svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where 180
172 s1 = svC m 181svdProp5b m = all (s1|~|) [s3,s5] where
173 s2 = svCd m 182 s1 = singularValues (m :: Matrix (Complex Double))
174 (_,s3,_) = svdC m 183-- s2 = svCd m
175 (_,s4,_) = svdCd m 184 (_,s3,_) = svd m
176 (_,s5,_) = thinSVDC m 185-- (_,s4,_) = svdCd m
177 (_,s6,_) = thinSVDCd m 186 (_,s5,_) = thinSVD m
187-- (_,s6,_) = thinSVDCd m
178 188
179svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' 189svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
180 where (u,s,v) = svdR m 190 where
181 (s',v') = rightSVR m 191 (u,s,v) = svd (m :: Matrix Double)
182 (u',s'') = leftSVR m 192 (s',v') = rightSV m
193 (u',s'') = leftSV m
183 194
184svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' 195svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
185 where (u,s,v) = svdC m 196 where
186 (s',v') = rightSVC m 197 (u,s,v) = svd (m :: Matrix (Complex Double))
187 (u',s'') = leftSVC m 198 (s',v') = rightSV m
199 (u',s'') = leftSV m
188 200
189svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s''' 201svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''
190 where (u,s,v) = svd m 202 where
191 (s',v') = rightSV m 203 (u,s,v) = svd m
192 (u',_s'') = leftSV m 204 (s',v') = rightSV m
193 s''' = singularValues m 205 (u',_s'') = leftSV m
206 s''' = singularValues m
194 207
195------------------------------------------------------------------ 208------------------------------------------------------------------
196 209
@@ -199,7 +212,7 @@ eigProp m = complex m <> v |~| v <> diag s
199 212
200eigSHProp m = m <> v |~| v <> real (diag s) 213eigSHProp m = m <> v |~| v <> real (diag s)
201 && unitary v 214 && unitary v
202 && m |~| v <> real (diag s) <> ctrans v 215 && m |~| v <> real (diag s) <> tr v
203 where (s, v) = eigSH m 216 where (s, v) = eigSH m
204 217
205eigProp2 m = fst (eig m) |~| eigenvalues m 218eigProp2 m = fst (eig m) |~| eigenvalues m
@@ -224,19 +237,19 @@ rqProp3 m = upperTriang' r
224 where (r,_q) = rq m 237 where (r,_q) = rq m
225 238
226upperTriang' r = upptr (rows r) (cols r) * r |~| r 239upperTriang' r = upptr (rows r) (cols r) * r |~| r
227 where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1 240 where upptr f c = build (f,c) $ \r' c' -> if r'-t > c' then 0 else 1
228 where t = f-c 241 where t = fromIntegral (f-c)
229 242
230hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h 243hessProp m = m |~| p <> h <> tr p && unitary p && upperHessenberg h
231 where (p,h) = hess m 244 where (p,h) = hess m
232 245
233schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s 246schurProp1 m = m |~| u <> s <> tr u && unitary u && upperTriang s
234 where (u,s) = schur m 247 where (u,s) = schur m
235 248
236schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme 249schurProp2 m = m |~| u <> s <> tr u && unitary u && upperHessenberg s -- fixme
237 where (u,s) = schur m 250 where (u,s) = schur m
238 251
239cholProp m = m |~| ctrans c <> c && upperTriang c 252cholProp m = m |~| tr c <> c && upperTriang c
240 where c = chol m 253 where c = chol m
241 254
242exactProp m = chol m == chol (m+0) 255exactProp m = chol m == chol (m+0)
@@ -250,7 +263,7 @@ mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ]
250 263
251multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) 264multProp1 p (a,b) = (a <> b) :~p~: (mulH a b)
252 265
253multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a) 266multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a)
254 267
255linearSolveProp f m = f m m |~| ident (rows m) 268linearSolveProp f m = f m m |~| ident (rows m)
256 269
@@ -259,5 +272,5 @@ linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
259 b = a <> x 272 b = a <> x
260 wc = rank a == q 273 wc = rank a == q
261 274
262subProp m = m == (trans . fromColumns . toRows) m 275subProp m = m == (tr . fromColumns . toRows) m
263 276