diff options
Diffstat (limited to 'packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs')
-rw-r--r-- | packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs | 153 |
1 files changed, 80 insertions, 73 deletions
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 | ||
46 | import Numeric.Container | 45 | import Numeric.LinearAlgebra.HMatrix hiding (Testable,unitary) |
47 | import Numeric.LinearAlgebra --hiding (real,complex) | 46 | import Test.QuickCheck |
48 | import Numeric.LinearAlgebra.LAPACK | ||
49 | import Debug.Trace | ||
50 | import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector | ||
51 | ,sized,classify,Testable,Property | ||
52 | ,quickCheckWith,maxSize,stdArgs,shrink) | ||
53 | 47 | ||
54 | trivial :: Testable a => Bool -> a -> Property | 48 | trivial :: Testable a => Bool -> a -> Property |
55 | trivial = (`classify` "trivial") | 49 | trivial = (`classify` "trivial") |
56 | 50 | ||
57 | -- relative error | 51 | -- relative error |
58 | dist :: (Normed c t, Num (c t)) => c t -> c t -> Double | 52 | dist :: (Num a, Normed a) => a -> a -> Double |
59 | dist = relativeError Infinity | 53 | dist = relativeError norm_Inf |
60 | 54 | ||
61 | infixl 4 |~| | 55 | infixl 4 |~| |
62 | a |~| b = a :~10~: b | 56 | a |~| b = a :~10~: b |
@@ -73,11 +67,11 @@ a :~n~: b = dist a b < 10^^(-n) | |||
73 | square m = rows m == cols m | 67 | square m = rows m == cols m |
74 | 68 | ||
75 | -- orthonormal columns | 69 | -- orthonormal columns |
76 | orthonormal m = ctrans m <> m |~| ident (cols m) | 70 | orthonormal m = tr m <> m |~| ident (cols m) |
77 | 71 | ||
78 | unitary m = square m && orthonormal m | 72 | unitary m = square m && orthonormal m |
79 | 73 | ||
80 | hermitian m = square m && m |~| ctrans m | 74 | hermitian m = square m && m |~| tr m |
81 | 75 | ||
82 | wellCond m = rcond m > 1/100 | 76 | wellCond 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 | ||
87 | upperTriang m = rows m == 1 || down == z | 81 | upperTriang 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 | ||
91 | upperHessenberg m = rows m < 3 || down == z | 85 | upperHessenberg 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 | ||
95 | zeros (r,c) = reshape c (konst 0 (r*c)) | 89 | zeros (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 | ||
120 | nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c) | 114 | nullspaceProp 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 |
131 | bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v | 125 | bugProp 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 |
139 | svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v | 133 | svdProp1 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 | ||
142 | svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where | 138 | svdProp1a 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 | ||
146 | svdProp1b svdfun m = unitary u && unitary v where | 143 | svdProp1b svdfun m = unitary u && unitary v |
144 | where | ||
147 | (u,_,v) = svdfun m | 145 | (u,_,v) = svdfun m |
148 | 146 | ||
149 | -- thinSVD | 147 | -- thinSVD |
150 | svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m) | 148 | svdProp2 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 |
154 | svdProp3 m = (m |~| u <> real (diag s) <> trans v | 156 | svdProp3 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 | ||
158 | svdProp4 m' = m |~| u <> real (diag s) <> trans v | 161 | svdProp4 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' | |
165 | svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where | 168 | |
166 | s1 = svR m | 169 | svdProp5a 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 | |
173 | svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where | 176 | |
174 | s1 = svC m | 177 | svdProp5b 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 | ||
181 | svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' | 185 | svdProp6a 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 | ||
186 | svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' | 191 | svdProp6b 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 | ||
191 | svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s''' | 197 | svdProp7 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 | ||
202 | eigSHProp m = m <> v |~| v <> real (diag s) | 209 | eigSHProp 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 | ||
207 | eigProp2 m = fst (eig m) |~| eigenvalues m | 214 | eigProp2 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 | ||
228 | upperTriang' r = upptr (rows r) (cols r) * r |~| r | 235 | upperTriang' 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 | ||
232 | hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h | 239 | hessProp m = m |~| p <> h <> tr p && unitary p && upperHessenberg h |
233 | where (p,h) = hess m | 240 | where (p,h) = hess m |
234 | 241 | ||
235 | schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s | 242 | schurProp1 m = m |~| u <> s <> tr u && unitary u && upperTriang s |
236 | where (u,s) = schur m | 243 | where (u,s) = schur m |
237 | 244 | ||
238 | schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme | 245 | schurProp2 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 | ||
241 | cholProp m = m |~| ctrans c <> c && upperTriang c | 248 | cholProp m = m |~| tr c <> c && upperTriang c |
242 | where c = chol m | 249 | where c = chol m |
243 | 250 | ||
244 | exactProp m = chol m == chol (m+0) | 251 | exactProp 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 | ||
253 | multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) | 260 | multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) |
254 | 261 | ||
255 | multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a) | 262 | multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a) |
256 | 263 | ||
257 | linearSolveProp f m = f m m |~| ident (rows m) | 264 | linearSolveProp 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 | ||
264 | subProp m = m == (trans . fromColumns . toRows) m | 271 | subProp m = m == (conj . tr . fromColumns . toRows) m |
265 | 272 | ||