diff options
Diffstat (limited to 'packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs')
-rw-r--r-- | packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs | 170 |
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 | ||
17 | module Numeric.LinearAlgebra.Tests.Properties ( | 16 | module 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 | ||
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 | 47 | |
49 | import Debug.Trace | 48 | (~=) :: Double -> Double -> Bool |
50 | import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector | 49 | a ~= b = abs (a - b) < 1e-10 |
51 | ,sized,classify,Testable,Property | ||
52 | ,quickCheckWith,maxSize,stdArgs,shrink) | ||
53 | 50 | ||
54 | trivial :: Testable a => Bool -> a -> Property | 51 | trivial :: Testable a => Bool -> a -> Property |
55 | trivial = (`classify` "trivial") | 52 | trivial = (`classify` "trivial") |
56 | 53 | ||
57 | -- relative error | 54 | -- relative error |
58 | dist :: (Normed c t, Num (c t)) => c t -> c t -> Double | 55 | dist :: (Num a, Normed a) => a -> a -> Double |
59 | dist = relativeError Infinity | 56 | dist = relativeError norm_Inf |
60 | 57 | ||
61 | infixl 4 |~| | 58 | infixl 4 |~| |
62 | a |~| b = a :~10~: b | 59 | a |~| b = a :~10~: b |
@@ -73,11 +70,11 @@ a :~n~: b = dist a b < 10^^(-n) | |||
73 | square m = rows m == cols m | 70 | square m = rows m == cols m |
74 | 71 | ||
75 | -- orthonormal columns | 72 | -- orthonormal columns |
76 | orthonormal m = ctrans m <> m |~| ident (cols m) | 73 | orthonormal m = tr m <> m |~| ident (cols m) |
77 | 74 | ||
78 | unitary m = square m && orthonormal m | 75 | unitary m = square m && orthonormal m |
79 | 76 | ||
80 | hermitian m = square m && m |~| ctrans m | 77 | hermitian m = square m && m |~| tr m |
81 | 78 | ||
82 | wellCond m = rcond m > 1/100 | 79 | wellCond 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 | ||
87 | upperTriang m = rows m == 1 || down == z | 84 | upperTriang 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 | ||
91 | upperHessenberg m = rows m < 3 || down == z | 88 | upperHessenberg 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 | ||
95 | zeros (r,c) = reshape c (konst 0 (r*c)) | 92 | zeros (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 | ||
120 | nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c) | 117 | nullspaceProp 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 |
131 | bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v | 128 | bugProp 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 |
139 | svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v | 136 | svdProp1 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 | ||
142 | svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where | 141 | svdProp1a 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 | ||
146 | svdProp1b svdfun m = unitary u && unitary v where | 146 | svdProp1b svdfun m = unitary u && unitary v |
147 | where | ||
147 | (u,_,v) = svdfun m | 148 | (u,_,v) = svdfun m |
148 | 149 | ||
149 | -- thinSVD | 150 | -- thinSVD |
150 | svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m) | 151 | svdProp2 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 |
154 | svdProp3 m = (m |~| u <> real (diag s) <> trans v | 159 | svdProp3 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 | ||
158 | svdProp4 m' = m |~| u <> real (diag s) <> trans v | 164 | svdProp4 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' | |
165 | svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where | 171 | |
166 | s1 = svR m | 172 | svdProp5a 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 | |
173 | svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where | 179 | |
174 | s1 = svC m | 180 | svdProp5b 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 | ||
181 | svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' | 188 | svdProp6a 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 | ||
186 | svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' | 194 | svdProp6b 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 | ||
191 | svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s''' | 200 | svdProp7 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 | ||
202 | eigSHProp m = m <> v |~| v <> real (diag s) | 212 | eigSHProp 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 | ||
207 | eigProp2 m = fst (eig m) |~| eigenvalues m | 217 | eigProp2 m = fst (eig m) |~| eigenvalues m |
208 | 218 | ||
209 | eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m | 219 | eigSHProp2 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 | ||
228 | upperTriang' r = upptr (rows r) (cols r) * r |~| r | 238 | 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 | 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 | ||
232 | hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h | 242 | hessProp m = m |~| p <> h <> tr p && unitary p && upperHessenberg h |
233 | where (p,h) = hess m | 243 | where (p,h) = hess m |
234 | 244 | ||
235 | schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s | 245 | schurProp1 m = m |~| u <> s <> tr u && unitary u && upperTriang s |
236 | where (u,s) = schur m | 246 | where (u,s) = schur m |
237 | 247 | ||
238 | schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme | 248 | schurProp2 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 | ||
241 | cholProp m = m |~| ctrans c <> c && upperTriang c | 251 | cholProp m = m |~| tr c <> c && upperTriang c |
242 | where c = chol m | 252 | where c = chol (trustSym m) |
243 | 253 | ||
244 | exactProp m = chol m == chol (m+0) | 254 | exactProp m = chol (trustSym m) == chol (trustSym (m+0)) |
245 | 255 | ||
246 | expmDiagProp m = expm (logm m) :~ 7 ~: complex m | 256 | expmDiagProp 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 | ||
253 | multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) | 263 | multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) |
254 | 264 | ||
255 | multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a) | 265 | multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a) |
256 | 266 | ||
257 | linearSolveProp f m = f m m |~| ident (rows m) | 267 | linearSolveProp f m = f m m |~| ident (rows m) |
258 | 268 | ||
269 | linearSolvePropH f m = f m (unSym m) |~| ident (rows (unSym m)) | ||
270 | |||
259 | linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b) | 271 | linearSolveProp2 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 | ||
264 | subProp m = m == (trans . fromColumns . toRows) m | 276 | subProp m = m == (conj . tr . fromColumns . toRows) m |
265 | 277 | ||