summaryrefslogtreecommitdiff
path: root/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
diff options
context:
space:
mode:
Diffstat (limited to 'packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs')
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs170
1 files changed, 91 insertions, 79 deletions
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
index a5c37f4..046644f 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-----------------------------------------------------------------------------
@@ -15,7 +14,7 @@ Testing properties.
15-} 14-}
16 15
17module Numeric.LinearAlgebra.Tests.Properties ( 16module Numeric.LinearAlgebra.Tests.Properties (
18 dist, (|~|), (~~), (~:), Aprox((:~)), 17 dist, (|~|), (~~), (~:), Aprox((:~)), (~=),
19 zeros, ones, 18 zeros, ones,
20 square, 19 square,
21 unitary, 20 unitary,
@@ -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,
@@ -40,23 +39,21 @@ module Numeric.LinearAlgebra.Tests.Properties (
40 expmDiagProp, 39 expmDiagProp,
41 multProp1, multProp2, 40 multProp1, multProp2,
42 subProp, 41 subProp,
43 linearSolveProp, linearSolveProp2 42 linearSolveProp, linearSolvePropH, 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 47
49import Debug.Trace 48(~=) :: Double -> Double -> Bool
50import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector 49a ~= b = abs (a - b) < 1e-10
51 ,sized,classify,Testable,Property
52 ,quickCheckWith,maxSize,stdArgs,shrink)
53 50
54trivial :: Testable a => Bool -> a -> Property 51trivial :: Testable a => Bool -> a -> Property
55trivial = (`classify` "trivial") 52trivial = (`classify` "trivial")
56 53
57-- relative error 54-- relative error
58dist :: (Normed c t, Num (c t)) => c t -> c t -> Double 55dist :: (Num a, Normed a) => a -> a -> Double
59dist = relativeError Infinity 56dist = relativeError norm_Inf
60 57
61infixl 4 |~| 58infixl 4 |~|
62a |~| b = a :~10~: b 59a |~| b = a :~10~: b
@@ -73,11 +70,11 @@ a :~n~: b = dist a b < 10^^(-n)
73square m = rows m == cols m 70square m = rows m == cols m
74 71
75-- orthonormal columns 72-- orthonormal columns
76orthonormal m = ctrans m <> m |~| ident (cols m) 73orthonormal m = tr m <> m |~| ident (cols m)
77 74
78unitary m = square m && orthonormal m 75unitary m = square m && orthonormal m
79 76
80hermitian m = square m && m |~| ctrans m 77hermitian m = square m && m |~| tr m
81 78
82wellCond m = rcond m > 1/100 79wellCond m = rcond m > 1/100
83 80
@@ -85,12 +82,12 @@ positiveDefinite m = minimum (toList e) > 0
85 where (e,_v) = eigSH m 82 where (e,_v) = eigSH m
86 83
87upperTriang m = rows m == 1 || down == z 84upperTriang m = rows m == 1 || down == z
88 where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m)) 85 where down = fromList $ concat $ zipWith drop [1..] (toLists (tr m))
89 z = konst 0 (dim down) 86 z = konst 0 (size down)
90 87
91upperHessenberg m = rows m < 3 || down == z 88upperHessenberg m = rows m < 3 || down == z
92 where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m)) 89 where down = fromList $ concat $ zipWith drop [2..] (toLists (tr m))
93 z = konst 0 (dim down) 90 z = konst 0 (size down)
94 91
95zeros (r,c) = reshape c (konst 0 (r*c)) 92zeros (r,c) = reshape c (konst 0 (r*c))
96 93
@@ -118,81 +115,94 @@ detProp m = s d1 |~| s d2
118 s x = fromList [x] 115 s x = fromList [x]
119 116
120nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c) 117nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c)
121 && orthonormal (fromColumns nl)) 118 && orthonormal n)
122 where nl = nullspacePrec 1 m 119 where n = nullspaceSVD (Left (1*peps)) m (rightSV m)
123 n = fromColumns nl 120 nl = toColumns n
124 r = rows m 121 r = rows m
125 c = cols m - rank m 122 c = cols m - rank m
126 123
127------------------------------------------------------------------ 124------------------------------------------------------------------
128 125{-
129-- testcase for nonempty fpu stack 126-- testcase for nonempty fpu stack
130-- uncommenting unitary' signature eliminates the problem 127-- uncommenting unitary' signature eliminates the problem
131bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v 128bugProp m = m |~| u <> real d <> tr v && unitary' u && unitary' v
132 where (u,d,v) = fullSVD m 129 where (u,d,v) = svd m
133 -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool 130 -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool
134 unitary' a = unitary a 131 unitary' a = unitary a
135 132-}
136------------------------------------------------------------------ 133------------------------------------------------------------------
137 134
138-- fullSVD 135-- fullSVD
139svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v 136svdProp1 m = m |~| u <> real d <> tr v && unitary u && unitary v
140 where (u,d,v) = fullSVD m 137 where
138 (u,s,v) = svd m
139 d = diagRect 0 s (rows m) (cols m)
141 140
142svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where 141svdProp1a svdfun m = m |~| u <> real d <> tr v && unitary u && unitary v
142 where
143 (u,s,v) = svdfun m 143 (u,s,v) = svdfun m
144 d = diagRect 0 s (rows m) (cols m) 144 d = diagRect 0 s (rows m) (cols m)
145 145
146svdProp1b svdfun m = unitary u && unitary v where 146svdProp1b svdfun m = unitary u && unitary v
147 where
147 (u,_,v) = svdfun m 148 (u,_,v) = svdfun m
148 149
149-- thinSVD 150-- thinSVD
150svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m) 151svdProp2 thinSVDfun m
151 where (u,s,v) = thinSVDfun m 152 = m |~| u <> diag (real s) <> tr v
153 && orthonormal u && orthonormal v
154 && size s == min (rows m) (cols m)
155 where
156 (u,s,v) = thinSVDfun m
152 157
153-- compactSVD 158-- compactSVD
154svdProp3 m = (m |~| u <> real (diag s) <> trans v 159svdProp3 m = (m |~| u <> real (diag s) <> tr v
155 && orthonormal u && orthonormal v) 160 && orthonormal u && orthonormal v)
156 where (u,s,v) = compactSVD m 161 where
162 (u,s,v) = compactSVD m
157 163
158svdProp4 m' = m |~| u <> real (diag s) <> trans v 164svdProp4 m' = m |~| u <> real (diag s) <> tr v
159 && orthonormal u && orthonormal v 165 && orthonormal u && orthonormal v
160 && (dim s == r || r == 0 && dim s == 1) 166 && (size s == r || r == 0 && size s == 1)
161 where (u,s,v) = compactSVD m 167 where
162 m = fromBlocks [[m'],[m']] 168 (u,s,v) = compactSVD m
163 r = rank m' 169 m = fromBlocks [[m'],[m']]
164 170 r = rank m'
165svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where 171
166 s1 = svR m 172svdProp5a m = all (s1|~|) [s3,s5] where
167 s2 = svRd m 173 s1 = singularValues (m :: Matrix Double)
168 (_,s3,_) = svdR m 174-- s2 = svRd m
169 (_,s4,_) = svdRd m 175 (_,s3,_) = svd m
170 (_,s5,_) = thinSVDR m 176-- (_,s4,_) = svdRd m
171 (_,s6,_) = thinSVDRd m 177 (_,s5,_) = thinSVD m
172 178-- (_,s6,_) = thinSVDRd m
173svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where 179
174 s1 = svC m 180svdProp5b m = all (s1|~|) [s3,s5] where
175 s2 = svCd m 181 s1 = singularValues (m :: Matrix (Complex Double))
176 (_,s3,_) = svdC m 182-- s2 = svCd m
177 (_,s4,_) = svdCd m 183 (_,s3,_) = svd m
178 (_,s5,_) = thinSVDC m 184-- (_,s4,_) = svdCd m
179 (_,s6,_) = thinSVDCd m 185 (_,s5,_) = thinSVD m
186-- (_,s6,_) = thinSVDCd m
180 187
181svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' 188svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
182 where (u,s,v) = svdR m 189 where
183 (s',v') = rightSVR m 190 (u,s,v) = svd (m :: Matrix Double)
184 (u',s'') = leftSVR m 191 (s',v') = rightSV m
192 (u',s'') = leftSV m
185 193
186svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' 194svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
187 where (u,s,v) = svdC m 195 where
188 (s',v') = rightSVC m 196 (u,s,v) = svd (m :: Matrix (Complex Double))
189 (u',s'') = leftSVC m 197 (s',v') = rightSV m
198 (u',s'') = leftSV m
190 199
191svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s''' 200svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''
192 where (u,s,v) = svd m 201 where
193 (s',v') = rightSV m 202 (u,s,v) = svd m
194 (u',_s'') = leftSV m 203 (s',v') = rightSV m
195 s''' = singularValues m 204 (u',_s'') = leftSV m
205 s''' = singularValues m
196 206
197------------------------------------------------------------------ 207------------------------------------------------------------------
198 208
@@ -201,12 +211,12 @@ eigProp m = complex m <> v |~| v <> diag s
201 211
202eigSHProp m = m <> v |~| v <> real (diag s) 212eigSHProp m = m <> v |~| v <> real (diag s)
203 && unitary v 213 && unitary v
204 && m |~| v <> real (diag s) <> ctrans v 214 && m |~| v <> real (diag s) <> tr v
205 where (s, v) = eigSH m 215 where (s, v) = eigSH' m
206 216
207eigProp2 m = fst (eig m) |~| eigenvalues m 217eigProp2 m = fst (eig m) |~| eigenvalues m
208 218
209eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m 219eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m
210 220
211------------------------------------------------------------------ 221------------------------------------------------------------------
212 222
@@ -226,22 +236,22 @@ rqProp3 m = upperTriang' r
226 where (r,_q) = rq m 236 where (r,_q) = rq m
227 237
228upperTriang' r = upptr (rows r) (cols r) * r |~| r 238upperTriang' 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 239 where upptr f c = build (f,c) $ \r' c' -> if r'-t > c' then 0 else 1
230 where t = f-c 240 where t = fromIntegral (f-c)
231 241
232hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h 242hessProp m = m |~| p <> h <> tr p && unitary p && upperHessenberg h
233 where (p,h) = hess m 243 where (p,h) = hess m
234 244
235schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s 245schurProp1 m = m |~| u <> s <> tr u && unitary u && upperTriang s
236 where (u,s) = schur m 246 where (u,s) = schur m
237 247
238schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme 248schurProp2 m = m |~| u <> s <> tr u && unitary u && upperHessenberg s -- fixme
239 where (u,s) = schur m 249 where (u,s) = schur m
240 250
241cholProp m = m |~| ctrans c <> c && upperTriang c 251cholProp m = m |~| tr c <> c && upperTriang c
242 where c = chol m 252 where c = chol (trustSym m)
243 253
244exactProp m = chol m == chol (m+0) 254exactProp m = chol (trustSym m) == chol (trustSym (m+0))
245 255
246expmDiagProp m = expm (logm m) :~ 7 ~: complex m 256expmDiagProp m = expm (logm m) :~ 7 ~: complex m
247 where logm = matFunc log 257 where logm = matFunc log
@@ -252,14 +262,16 @@ mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ]
252 262
253multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) 263multProp1 p (a,b) = (a <> b) :~p~: (mulH a b)
254 264
255multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a) 265multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a)
256 266
257linearSolveProp f m = f m m |~| ident (rows m) 267linearSolveProp f m = f m m |~| ident (rows m)
258 268
269linearSolvePropH f m = f m (unSym m) |~| ident (rows (unSym m))
270
259linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b) 271linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
260 where q = min (rows a) (cols a) 272 where q = min (rows a) (cols a)
261 b = a <> x 273 b = a <> x
262 wc = rank a == q 274 wc = rank a == q
263 275
264subProp m = m == (trans . fromColumns . toRows) m 276subProp m = m == (conj . tr . fromColumns . toRows) m
265 277