diff options
Diffstat (limited to 'packages/tests/src/Numeric')
-rw-r--r-- | packages/tests/src/Numeric/GSL/Tests.hs | 4 | ||||
-rw-r--r-- | packages/tests/src/Numeric/LinearAlgebra/Tests.hs | 202 | ||||
-rw-r--r-- | packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs | 91 | ||||
-rw-r--r-- | packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs | 153 |
4 files changed, 183 insertions, 267 deletions
diff --git a/packages/tests/src/Numeric/GSL/Tests.hs b/packages/tests/src/Numeric/GSL/Tests.hs index 9dff6f5..e5d205d 100644 --- a/packages/tests/src/Numeric/GSL/Tests.hs +++ b/packages/tests/src/Numeric/GSL/Tests.hs | |||
@@ -19,7 +19,7 @@ import System.Exit (exitFailure) | |||
19 | 19 | ||
20 | import Test.HUnit (runTestTT, failures, Test(..), errors) | 20 | import Test.HUnit (runTestTT, failures, Test(..), errors) |
21 | 21 | ||
22 | import Numeric.LinearAlgebra | 22 | import Numeric.LinearAlgebra.HMatrix |
23 | import Numeric.GSL | 23 | import Numeric.GSL |
24 | import Numeric.LinearAlgebra.Tests (qCheck, utest) | 24 | import Numeric.LinearAlgebra.Tests (qCheck, utest) |
25 | import Numeric.LinearAlgebra.Tests.Properties ((|~|), (~~)) | 25 | import Numeric.LinearAlgebra.Tests.Properties ((|~|), (~~)) |
@@ -42,7 +42,7 @@ fittingTest = utest "levmar" (ok1 && ok2) | |||
42 | sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] | 42 | sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] |
43 | 43 | ||
44 | ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d | 44 | ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d |
45 | ok2 = norm2 (fromList (map fst sols) - fromList sol) < 1E-5 | 45 | ok2 = norm_2 (fromList (map fst sols) - fromList sol) < 1E-5 |
46 | 46 | ||
47 | --------------------------------------------------------------------- | 47 | --------------------------------------------------------------------- |
48 | 48 | ||
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs index 713af79..71c7c16 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs | |||
@@ -28,12 +28,9 @@ module Numeric.LinearAlgebra.Tests( | |||
28 | --, runBigTests | 28 | --, runBigTests |
29 | ) where | 29 | ) where |
30 | 30 | ||
31 | import Numeric.LinearAlgebra | 31 | import Numeric.LinearAlgebra hiding (unitary) |
32 | import Numeric.LinearAlgebra.HMatrix hiding ((<>),linearSolve) | 32 | import Numeric.LinearAlgebra.Devel hiding (vec) |
33 | import Numeric.LinearAlgebra.Static(L) | 33 | import Numeric.LinearAlgebra.Static(L) |
34 | import Numeric.LinearAlgebra.Util(col,row) | ||
35 | import Data.Packed | ||
36 | import Numeric.LinearAlgebra.LAPACK | ||
37 | import Numeric.LinearAlgebra.Tests.Instances | 34 | import Numeric.LinearAlgebra.Tests.Instances |
38 | import Numeric.LinearAlgebra.Tests.Properties | 35 | import Numeric.LinearAlgebra.Tests.Properties |
39 | import Test.HUnit hiding ((~:),test,Testable,State) | 36 | import Test.HUnit hiding ((~:),test,Testable,State) |
@@ -44,16 +41,13 @@ import qualified Prelude | |||
44 | import System.CPUTime | 41 | import System.CPUTime |
45 | import System.Exit | 42 | import System.Exit |
46 | import Text.Printf | 43 | import Text.Printf |
47 | import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr) | 44 | import Numeric.LinearAlgebra.Devel(unsafeFromForeignPtr,unsafeToForeignPtr) |
48 | import Control.Arrow((***)) | 45 | import Control.Arrow((***)) |
49 | import Debug.Trace | 46 | import Debug.Trace |
50 | import Control.Monad(when) | 47 | import Control.Monad(when) |
51 | import Numeric.LinearAlgebra.Util hiding (ones,row,col) | ||
52 | import Control.Applicative | 48 | import Control.Applicative |
53 | import Control.Monad(ap) | 49 | import Control.Monad(ap) |
54 | 50 | ||
55 | import Data.Packed.ST | ||
56 | |||
57 | import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector | 51 | import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector |
58 | ,sized,classify,Testable,Property | 52 | ,sized,classify,Testable,Property |
59 | ,quickCheckWithResult,maxSize,stdArgs,shrink) | 53 | ,quickCheckWithResult,maxSize,stdArgs,shrink) |
@@ -89,7 +83,7 @@ detTest1 = det m == 26 | |||
89 | mc = (3><3) | 83 | mc = (3><3) |
90 | [ 1, 2, 3 | 84 | [ 1, 2, 3 |
91 | , 4, 5, 7 | 85 | , 4, 5, 7 |
92 | , 2, 8, i | 86 | , 2, 8, iC |
93 | ] | 87 | ] |
94 | 88 | ||
95 | detTest2 = inv1 |~| inv2 && [det1] ~~ [det2] | 89 | detTest2 = inv1 |~| inv2 && [det1] ~~ [det2] |
@@ -140,7 +134,7 @@ randomTestGaussian = c :~1~: snd (meanCov dat) where | |||
140 | 2,4,0, | 134 | 2,4,0, |
141 | -2,2,1] | 135 | -2,2,1] |
142 | m = 3 |> [1,2,3] | 136 | m = 3 |> [1,2,3] |
143 | c = a <> trans a | 137 | c = a <> tr a |
144 | dat = gaussianSample 7 (10^6) m c | 138 | dat = gaussianSample 7 (10^6) m c |
145 | 139 | ||
146 | randomTestUniform = c :~1~: snd (meanCov dat) where | 140 | randomTestUniform = c :~1~: snd (meanCov dat) where |
@@ -174,54 +168,54 @@ offsetTest = y == y' where | |||
174 | 168 | ||
175 | normsVTest = TestList [ | 169 | normsVTest = TestList [ |
176 | utest "normv2CD" $ norm2PropC v | 170 | utest "normv2CD" $ norm2PropC v |
177 | , utest "normv2CF" $ norm2PropC (single v) | 171 | -- , utest "normv2CF" $ norm2PropC (single v) |
178 | #ifndef NONORMVTEST | 172 | #ifndef NONORMVTEST |
179 | , utest "normv2D" $ norm2PropR x | 173 | , utest "normv2D" $ norm2PropR x |
180 | , utest "normv2F" $ norm2PropR (single x) | 174 | -- , utest "normv2F" $ norm2PropR (single x) |
181 | #endif | 175 | #endif |
182 | , utest "normv1CD" $ norm1 v == 8 | 176 | , utest "normv1CD" $ norm_1 v == 8 |
183 | , utest "normv1CF" $ norm1 (single v) == 8 | 177 | -- , utest "normv1CF" $ norm_1 (single v) == 8 |
184 | , utest "normv1D" $ norm1 x == 6 | 178 | , utest "normv1D" $ norm_1 x == 6 |
185 | , utest "normv1F" $ norm1 (single x) == 6 | 179 | -- , utest "normv1F" $ norm_1 (single x) == 6 |
186 | 180 | ||
187 | , utest "normvInfCD" $ normInf v == 5 | 181 | , utest "normvInfCD" $ norm_Inf v == 5 |
188 | , utest "normvInfCF" $ normInf (single v) == 5 | 182 | -- , utest "normvInfCF" $ norm_Inf (single v) == 5 |
189 | , utest "normvInfD" $ normInf x == 3 | 183 | , utest "normvInfD" $ norm_Inf x == 3 |
190 | , utest "normvInfF" $ normInf (single x) == 3 | 184 | -- , utest "normvInfF" $ norm_Inf (single x) == 3 |
191 | 185 | ||
192 | ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double) | 186 | ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double) |
193 | x = fromList [1,2,-3] :: Vector Double | 187 | x = fromList [1,2,-3] :: Vector Double |
194 | #ifndef NONORMVTEST | 188 | #ifndef NONORMVTEST |
195 | norm2PropR a = norm2 a =~= sqrt (udot a a) | 189 | norm2PropR a = norm_2 a =~= sqrt (udot a a) |
196 | #endif | 190 | #endif |
197 | norm2PropC a = norm2 a =~= realPart (sqrt (a <.> a)) | 191 | norm2PropC a = norm_2 a =~= realPart (sqrt (a `dot` a)) |
198 | a =~= b = fromList [a] |~| fromList [b] | 192 | a =~= b = fromList [a] |~| fromList [b] |
199 | 193 | ||
200 | normsMTest = TestList [ | 194 | normsMTest = TestList [ |
201 | utest "norm2mCD" $ pnorm PNorm2 v =~= 8.86164970498005 | 195 | utest "norm2mCD" $ norm_2 v =~= 8.86164970498005 |
202 | , utest "norm2mCF" $ pnorm PNorm2 (single v) =~= 8.86164970498005 | 196 | -- , utest "norm2mCF" $ norm_2 (single v) =~= 8.86164970498005 |
203 | , utest "norm2mD" $ pnorm PNorm2 x =~= 5.96667765076216 | 197 | , utest "norm2mD" $ norm_2 x =~= 5.96667765076216 |
204 | , utest "norm2mF" $ pnorm PNorm2 (single x) =~= 5.96667765076216 | 198 | -- , utest "norm2mF" $ norm_2 (single x) =~= 5.96667765076216 |
205 | 199 | ||
206 | , utest "norm1mCD" $ pnorm PNorm1 v == 9 | 200 | , utest "norm1mCD" $ norm_1 v == 9 |
207 | , utest "norm1mCF" $ pnorm PNorm1 (single v) == 9 | 201 | -- , utest "norm1mCF" $ norm_1 (single v) == 9 |
208 | , utest "norm1mD" $ pnorm PNorm1 x == 7 | 202 | , utest "norm1mD" $ norm_1 x == 7 |
209 | , utest "norm1mF" $ pnorm PNorm1 (single x) == 7 | 203 | -- , utest "norm1mF" $ norm_1 (single x) == 7 |
210 | 204 | ||
211 | , utest "normmInfCD" $ pnorm Infinity v == 12 | 205 | , utest "normmInfCD" $ norm_Inf v == 12 |
212 | , utest "normmInfCF" $ pnorm Infinity (single v) == 12 | 206 | -- , utest "normmInfCF" $ norm_Inf (single v) == 12 |
213 | , utest "normmInfD" $ pnorm Infinity x == 8 | 207 | , utest "normmInfD" $ norm_Inf x == 8 |
214 | , utest "normmInfF" $ pnorm Infinity (single x) == 8 | 208 | -- , utest "normmInfF" $ norm_Inf (single x) == 8 |
215 | 209 | ||
216 | , utest "normmFroCD" $ pnorm Frobenius v =~= 8.88819441731559 | 210 | , utest "normmFroCD" $ norm_Frob v =~= 8.88819441731559 |
217 | , utest "normmFroCF" $ pnorm Frobenius (single v) =~~= 8.88819441731559 | 211 | -- , utest "normmFroCF" $ norm_Frob (single v) =~~= 8.88819441731559 |
218 | , utest "normmFroD" $ pnorm Frobenius x =~= 6.24499799839840 | 212 | , utest "normmFroD" $ norm_Frob x =~= 6.24499799839840 |
219 | , utest "normmFroF" $ pnorm Frobenius (single x) =~~= 6.24499799839840 | 213 | -- , utest "normmFroF" $ norm_Frob (single x) =~~= 6.24499799839840 |
220 | 214 | ||
221 | ] where v = (2><2) [1,-2*i,3:+4,7] :: Matrix (Complex Double) | 215 | ] where v = (2><2) [1,-2*iC,3:+4,7] :: Matrix (Complex Double) |
222 | x = (2><2) [1,2,-3,5] :: Matrix Double | 216 | x = (2><2) [1,2,-3,5] :: Matrix Double |
223 | a =~= b = fromList [a] :~10~: fromList [b] | 217 | a =~= b = fromList [a] :~10~: fromList [b] |
224 | a =~~= b = fromList [a] :~5~: fromList [b] | 218 | -- a =~~= b = fromList [a] :~5~: fromList [b] |
225 | 219 | ||
226 | --------------------------------------------------------------------- | 220 | --------------------------------------------------------------------- |
227 | 221 | ||
@@ -236,7 +230,7 @@ sumprodTest = TestList [ | |||
236 | , utest "prodD" $ prodProp v | 230 | , utest "prodD" $ prodProp v |
237 | , utest "prodF" $ prodProp (single v) | 231 | , utest "prodF" $ prodProp (single v) |
238 | ] where v = fromList [1,2,3] :: Vector Double | 232 | ] where v = fromList [1,2,3] :: Vector Double |
239 | z = fromList [1,2-i,3+i] | 233 | z = fromList [1,2-iC,3+iC] |
240 | prodProp x = prodElements x == product (toList x) | 234 | prodProp x = prodElements x == product (toList x) |
241 | 235 | ||
242 | --------------------------------------------------------------------- | 236 | --------------------------------------------------------------------- |
@@ -250,7 +244,7 @@ chainTest = utest "chain" $ foldl1' (<>) ms |~| optimiseMult ms where | |||
250 | 244 | ||
251 | --------------------------------------------------------------------- | 245 | --------------------------------------------------------------------- |
252 | 246 | ||
253 | conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m) | 247 | conjuTest m = cmap conjugate (flatten (conj (tr m))) == flatten (tr m) |
254 | 248 | ||
255 | --------------------------------------------------------------------- | 249 | --------------------------------------------------------------------- |
256 | 250 | ||
@@ -306,7 +300,7 @@ lift_maybe m = MaybeT $ do | |||
306 | 300 | ||
307 | -- apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs | 301 | -- apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs |
308 | --successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool | 302 | --successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool |
309 | successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (dim v - 1) v))) (v @> 0) | 303 | successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (size v - 1) v))) (v ! 0) |
310 | where stp e = do | 304 | where stp e = do |
311 | ep <- lift_maybe $ state_get | 305 | ep <- lift_maybe $ state_get |
312 | if t e ep | 306 | if t e ep |
@@ -315,7 +309,7 @@ successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ s | |||
315 | 309 | ||
316 | -- operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input | 310 | -- operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input |
317 | --successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b | 311 | --successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b |
318 | successive f v = evalState (mapVectorM stp (subVector 1 (dim v - 1) v)) (v @> 0) | 312 | successive f v = evalState (mapVectorM stp (subVector 1 (size v - 1) v)) (v ! 0) |
319 | where stp e = do | 313 | where stp e = do |
320 | ep <- state_get | 314 | ep <- state_get |
321 | state_put e | 315 | state_put e |
@@ -377,23 +371,6 @@ convolutionTest = utest "convolution" ok | |||
377 | 371 | ||
378 | -------------------------------------------------------------------------------- | 372 | -------------------------------------------------------------------------------- |
379 | 373 | ||
380 | kroneckerTest = utest "kronecker" ok | ||
381 | where | ||
382 | a,x,b :: Matrix Double | ||
383 | a = (3><4) [1..] | ||
384 | x = (4><2) [3,5..] | ||
385 | b = (2><5) [0,5..] | ||
386 | v1 = vec (a <> x <> b) | ||
387 | v2 = (trans b `kronecker` a) <> vec x | ||
388 | s = trans b <> b | ||
389 | v3 = vec s | ||
390 | v4 = (dup 5 :: Matrix Double) <> vech s | ||
391 | ok = v1 == v2 && v3 == v4 | ||
392 | && vtrans 1 a == trans a | ||
393 | && vtrans (rows a) a == asColumn (vec a) | ||
394 | |||
395 | -------------------------------------------------------------------------------- | ||
396 | |||
397 | sparseTest = utest "sparse" (fst $ checkT (undefined :: GMatrix)) | 374 | sparseTest = utest "sparse" (fst $ checkT (undefined :: GMatrix)) |
398 | 375 | ||
399 | -------------------------------------------------------------------------------- | 376 | -------------------------------------------------------------------------------- |
@@ -435,11 +412,11 @@ runTests n = do | |||
435 | test (multProp1 10 . cConsist) | 412 | test (multProp1 10 . cConsist) |
436 | test (multProp2 10 . rConsist) | 413 | test (multProp2 10 . rConsist) |
437 | test (multProp2 10 . cConsist) | 414 | test (multProp2 10 . cConsist) |
438 | putStrLn "------ mult Float" | 415 | -- putStrLn "------ mult Float" |
439 | test (multProp1 6 . (single *** single) . rConsist) | 416 | -- test (multProp1 6 . (single *** single) . rConsist) |
440 | test (multProp1 6 . (single *** single) . cConsist) | 417 | -- test (multProp1 6 . (single *** single) . cConsist) |
441 | test (multProp2 6 . (single *** single) . rConsist) | 418 | -- test (multProp2 6 . (single *** single) . rConsist) |
442 | test (multProp2 6 . (single *** single) . cConsist) | 419 | -- test (multProp2 6 . (single *** single) . cConsist) |
443 | putStrLn "------ sub-trans" | 420 | putStrLn "------ sub-trans" |
444 | test (subProp . rM) | 421 | test (subProp . rM) |
445 | test (subProp . cM) | 422 | test (subProp . cM) |
@@ -472,16 +449,16 @@ runTests n = do | |||
472 | putStrLn "------ svd" | 449 | putStrLn "------ svd" |
473 | test (svdProp1 . rM) | 450 | test (svdProp1 . rM) |
474 | test (svdProp1 . cM) | 451 | test (svdProp1 . cM) |
475 | test (svdProp1a svdR) | 452 | test (svdProp1a svd . rM) |
476 | test (svdProp1a svdC) | 453 | test (svdProp1a svd . cM) |
477 | test (svdProp1a svdRd) | 454 | -- test (svdProp1a svdRd) |
478 | test (svdProp1b svdR) | 455 | test (svdProp1b svd . rM) |
479 | test (svdProp1b svdC) | 456 | test (svdProp1b svd . cM) |
480 | test (svdProp1b svdRd) | 457 | -- test (svdProp1b svdRd) |
481 | test (svdProp2 thinSVDR) | 458 | test (svdProp2 thinSVD . rM) |
482 | test (svdProp2 thinSVDC) | 459 | test (svdProp2 thinSVD . cM) |
483 | test (svdProp2 thinSVDRd) | 460 | -- test (svdProp2 thinSVDRd) |
484 | test (svdProp2 thinSVDCd) | 461 | -- test (svdProp2 thinSVDCd) |
485 | test (svdProp3 . rM) | 462 | test (svdProp3 . rM) |
486 | test (svdProp3 . cM) | 463 | test (svdProp3 . cM) |
487 | test (svdProp4 . rM) | 464 | test (svdProp4 . rM) |
@@ -492,12 +469,12 @@ runTests n = do | |||
492 | test (svdProp6b) | 469 | test (svdProp6b) |
493 | test (svdProp7 . rM) | 470 | test (svdProp7 . rM) |
494 | test (svdProp7 . cM) | 471 | test (svdProp7 . cM) |
495 | putStrLn "------ svdCd" | 472 | -- putStrLn "------ svdCd" |
496 | #ifdef NOZGESDD | 473 | #ifdef NOZGESDD |
497 | putStrLn "Omitted" | 474 | -- putStrLn "Omitted" |
498 | #else | 475 | #else |
499 | test (svdProp1a svdCd) | 476 | -- test (svdProp1a svdCd) |
500 | test (svdProp1b svdCd) | 477 | -- test (svdProp1b svdCd) |
501 | #endif | 478 | #endif |
502 | putStrLn "------ eig" | 479 | putStrLn "------ eig" |
503 | test (eigSHProp . rHer) | 480 | test (eigSHProp . rHer) |
@@ -515,10 +492,10 @@ runTests n = do | |||
515 | test (qrProp . rM) | 492 | test (qrProp . rM) |
516 | test (qrProp . cM) | 493 | test (qrProp . cM) |
517 | test (rqProp . rM) | 494 | test (rqProp . rM) |
518 | test (rqProp . cM) | 495 | -- test (rqProp . cM) |
519 | test (rqProp1 . cM) | 496 | test (rqProp1 . cM) |
520 | test (rqProp2 . cM) | 497 | test (rqProp2 . cM) |
521 | test (rqProp3 . cM) | 498 | -- test (rqProp3 . cM) |
522 | putStrLn "------ hess" | 499 | putStrLn "------ hess" |
523 | test (hessProp . rSq) | 500 | test (hessProp . rSq) |
524 | test (hessProp . cSq) | 501 | test (hessProp . cSq) |
@@ -539,12 +516,12 @@ runTests n = do | |||
539 | test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM)) | 516 | test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM)) |
540 | test (\u -> cos u * tan u |~| sin (u::RM)) | 517 | test (\u -> cos u * tan u |~| sin (u::RM)) |
541 | test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary | 518 | test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary |
542 | putStrLn "------ vector operations - Float" | 519 | -- putStrLn "------ vector operations - Float" |
543 | test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM)) | 520 | -- test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM)) |
544 | test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary | 521 | -- test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary |
545 | test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM)) | 522 | -- test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM)) |
546 | test (\u -> cos u * tan u |~~| sin (u::FM)) | 523 | -- test (\u -> cos u * tan u |~~| sin (u::FM)) |
547 | test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary | 524 | -- test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary |
548 | putStrLn "------ read . show" | 525 | putStrLn "------ read . show" |
549 | test (\m -> (m::RM) == read (show m)) | 526 | test (\m -> (m::RM) == read (show m)) |
550 | test (\m -> (m::CM) == read (show m)) | 527 | test (\m -> (m::CM) == read (show m)) |
@@ -562,8 +539,8 @@ runTests n = do | |||
562 | , utest "expm1" (expmTest1) | 539 | , utest "expm1" (expmTest1) |
563 | , utest "expm2" (expmTest2) | 540 | , utest "expm2" (expmTest2) |
564 | , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM) | 541 | , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM) |
565 | , utest "arith2" $ ((scalar (1+i) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*i-51) :: CM) | 542 | , utest "arith2" $ ((scalar (1+iC) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*iC-51) :: CM) |
566 | , utest "arith3" $ exp (scalar i * ones(10,10)*pi) + 1 |~| 0 | 543 | , utest "arith3" $ exp (scalar iC * ones(10,10)*pi) + 1 |~| 0 |
567 | , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3] | 544 | , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3] |
568 | -- , utest "gamma" (gamma 5 == 24.0) | 545 | -- , utest "gamma" (gamma 5 == 24.0) |
569 | -- , besselTest | 546 | -- , besselTest |
@@ -571,10 +548,10 @@ runTests n = do | |||
571 | , utest "randomGaussian" randomTestGaussian | 548 | , utest "randomGaussian" randomTestGaussian |
572 | , utest "randomUniform" randomTestUniform | 549 | , utest "randomUniform" randomTestUniform |
573 | , utest "buildVector/Matrix" $ | 550 | , utest "buildVector/Matrix" $ |
574 | complex (10 |> [0::Double ..]) == buildVector 10 fromIntegral | 551 | complex (10 |> [0::Double ..]) == build 10 id |
575 | && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0) | 552 | && ident 5 == build (5,5) (\r c -> if r==c then 1::Double else 0) |
576 | , utest "rank" $ rank ((2><3)[1,0,0,1,5*eps,0]) == 1 | 553 | , utest "rank" $ rank ((2><3)[1,0,0,1,5*peps,0::Double]) == 1 |
577 | && rank ((2><3)[1,0,0,1,7*eps,0]) == 2 | 554 | && rank ((2><3)[1,0,0,1,7*peps,0::Double]) == 2 |
578 | , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM) | 555 | , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM) |
579 | , mbCholTest | 556 | , mbCholTest |
580 | , utest "offset" offsetTest | 557 | , utest "offset" offsetTest |
@@ -588,7 +565,6 @@ runTests n = do | |||
588 | , conformTest | 565 | , conformTest |
589 | , accumTest | 566 | , accumTest |
590 | , convolutionTest | 567 | , convolutionTest |
591 | , kroneckerTest | ||
592 | , sparseTest | 568 | , sparseTest |
593 | , staticTest | 569 | , staticTest |
594 | ] | 570 | ] |
@@ -597,12 +573,12 @@ runTests n = do | |||
597 | 573 | ||
598 | 574 | ||
599 | -- single precision approximate equality | 575 | -- single precision approximate equality |
600 | infixl 4 |~~| | 576 | -- infixl 4 |~~| |
601 | a |~~| b = a :~6~: b | 577 | -- a |~~| b = a :~6~: b |
602 | 578 | ||
603 | makeUnitary v | realPart n > 1 = v / scalar n | 579 | makeUnitary v | realPart n > 1 = v / scalar n |
604 | | otherwise = v | 580 | | otherwise = v |
605 | where n = sqrt (v <.> v) | 581 | where n = sqrt (v `dot` v) |
606 | 582 | ||
607 | -- -- | Some additional tests on big matrices. They take a few minutes. | 583 | -- -- | Some additional tests on big matrices. They take a few minutes. |
608 | -- runBigTests :: IO () | 584 | -- runBigTests :: IO () |
@@ -668,9 +644,9 @@ manyvec5 xs = sumElements $ fromRows $ map (\x -> vec3 x (x**2) (x**3)) xs | |||
668 | 644 | ||
669 | 645 | ||
670 | manyvec2 xs = sum $ map (\x -> sqrt(x^2 + (x**2)^2 +(x**3)^2)) xs | 646 | manyvec2 xs = sum $ map (\x -> sqrt(x^2 + (x**2)^2 +(x**3)^2)) xs |
671 | manyvec3 xs = sum $ map (pnorm PNorm2 . (\x -> fromList [x,x**2,x**3])) xs | 647 | manyvec3 xs = sum $ map (norm_2 . (\x -> fromList [x,x**2,x**3])) xs |
672 | 648 | ||
673 | manyvec4 xs = sum $ map (pnorm PNorm2 . (\x -> vec3 x (x**2) (x**3))) xs | 649 | manyvec4 xs = sum $ map (norm_2 . (\x -> vec3 x (x**2) (x**3))) xs |
674 | 650 | ||
675 | vec3 :: Double -> Double -> Double -> Vector Double | 651 | vec3 :: Double -> Double -> Double -> Vector Double |
676 | vec3 a b c = runSTVector $ do | 652 | vec3 a b c = runSTVector $ do |
@@ -695,11 +671,11 @@ mkVecBench = do | |||
695 | 671 | ||
696 | subBench = do | 672 | subBench = do |
697 | putStrLn "" | 673 | putStrLn "" |
698 | let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v)) | 674 | let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (size v -1) v)) |
699 | time "0.1M subVector " (g (konst 1 (1+10^5) :: Vector Double) @> 0) | 675 | time "0.1M subVector " (g (konst 1 (1+10^5) :: Vector Double) ! 0) |
700 | let f = foldl1' (.) (replicate (10^5) (fromRows.toRows)) | 676 | let f = foldl1' (.) (replicate (10^5) (fromRows.toRows)) |
701 | time "subVector-join 3" (f (ident 3 :: Matrix Double) @@>(0,0)) | 677 | time "subVector-join 3" (f (ident 3 :: Matrix Double) `atIndex` (0,0)) |
702 | time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0)) | 678 | time "subVector-join 10" (f (ident 10 :: Matrix Double) `atIndex` (0,0)) |
703 | 679 | ||
704 | -------------------------------- | 680 | -------------------------------- |
705 | 681 | ||
@@ -724,7 +700,7 @@ multBench = do | |||
724 | 700 | ||
725 | eigBench = do | 701 | eigBench = do |
726 | let m = reshape 1000 (randomVector 777 Uniform (1000*1000)) | 702 | let m = reshape 1000 (randomVector 777 Uniform (1000*1000)) |
727 | s = m + trans m | 703 | s = m + tr m |
728 | m `seq` s `seq` putStrLn "" | 704 | m `seq` s `seq` putStrLn "" |
729 | time "eigenvalues symmetric 1000x1000" (eigenvaluesSH' m) | 705 | time "eigenvalues symmetric 1000x1000" (eigenvaluesSH' m) |
730 | time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m) | 706 | time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m) |
@@ -736,7 +712,7 @@ eigBench = do | |||
736 | svdBench = do | 712 | svdBench = do |
737 | let a = reshape 500 (randomVector 777 Uniform (3000*500)) | 713 | let a = reshape 500 (randomVector 777 Uniform (3000*500)) |
738 | b = reshape 1000 (randomVector 777 Uniform (1000*1000)) | 714 | b = reshape 1000 (randomVector 777 Uniform (1000*1000)) |
739 | fv (_,_,v) = v@@>(0,0) | 715 | fv (_,_,v) = v `atIndex` (0,0) |
740 | a `seq` b `seq` putStrLn "" | 716 | a `seq` b `seq` putStrLn "" |
741 | time "singular values 3000x500" (singularValues a) | 717 | time "singular values 3000x500" (singularValues a) |
742 | time "thin svd 3000x500" (fv $ thinSVD a) | 718 | time "thin svd 3000x500" (fv $ thinSVD a) |
@@ -748,7 +724,7 @@ svdBench = do | |||
748 | 724 | ||
749 | solveBenchN n = do | 725 | solveBenchN n = do |
750 | let x = uniformSample 777 (2*n) (replicate n (-1,1)) | 726 | let x = uniformSample 777 (2*n) (replicate n (-1,1)) |
751 | a = trans x <> x | 727 | a = tr x <> x |
752 | b = asColumn $ randomVector 666 Uniform n | 728 | b = asColumn $ randomVector 666 Uniform n |
753 | a `seq` b `seq` putStrLn "" | 729 | a `seq` b `seq` putStrLn "" |
754 | time ("svd solve " ++ show n) (linearSolveSVD a b) | 730 | time ("svd solve " ++ show n) (linearSolveSVD a b) |
@@ -765,7 +741,7 @@ solveBench = do | |||
765 | 741 | ||
766 | cholBenchN n = do | 742 | cholBenchN n = do |
767 | let x = uniformSample 777 (2*n) (replicate n (-1,1)) | 743 | let x = uniformSample 777 (2*n) (replicate n (-1,1)) |
768 | a = trans x <> x | 744 | a = tr x <> x |
769 | a `seq` putStr "" | 745 | a `seq` putStr "" |
770 | time ("chol " ++ show n) (chol a) | 746 | time ("chol " ++ show n) (chol a) |
771 | 747 | ||
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs index 53fc4d2..904ae05 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs | |||
@@ -1,5 +1,4 @@ | |||
1 | {-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-} | 1 | {-# LANGUAGE FlexibleContexts, UndecidableInstances, FlexibleInstances #-} |
2 | {-# OPTIONS_GHC -fno-warn-unused-imports #-} | ||
3 | ----------------------------------------------------------------------------- | 2 | ----------------------------------------------------------------------------- |
4 | {- | | 3 | {- | |
5 | Module : Numeric.LinearAlgebra.Tests.Instances | 4 | Module : Numeric.LinearAlgebra.Tests.Instances |
@@ -26,15 +25,11 @@ module Numeric.LinearAlgebra.Tests.Instances( | |||
26 | 25 | ||
27 | import System.Random | 26 | import System.Random |
28 | 27 | ||
29 | import Numeric.LinearAlgebra | 28 | import Numeric.LinearAlgebra.HMatrix hiding (vector) |
30 | import Numeric.LinearAlgebra.Devel | ||
31 | import Numeric.Container | ||
32 | import Control.Monad(replicateM) | 29 | import Control.Monad(replicateM) |
33 | import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector | 30 | import Test.QuickCheck(Arbitrary,arbitrary,choose,vector,sized,shrink) |
34 | ,sized,classify,Testable,Property | 31 | |
35 | ,quickCheckWith,maxSize,stdArgs,shrink) | ||
36 | 32 | ||
37 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
38 | shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]] | 33 | shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]] |
39 | shrinkListElementwise [] = [] | 34 | shrinkListElementwise [] = [] |
40 | shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ] | 35 | shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ] |
@@ -42,25 +37,6 @@ shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ] | |||
42 | 37 | ||
43 | shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)] | 38 | shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)] |
44 | shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ] | 39 | shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ] |
45 | #endif | ||
46 | |||
47 | #if MIN_VERSION_QuickCheck(2,1,1) | ||
48 | #else | ||
49 | instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where | ||
50 | arbitrary = do | ||
51 | re <- arbitrary | ||
52 | im <- arbitrary | ||
53 | return (re :+ im) | ||
54 | |||
55 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
56 | shrink (re :+ im) = | ||
57 | [ u :+ v | (u,v) <- shrinkPair (re,im) ] | ||
58 | #else | ||
59 | -- this has been moved to the 'Coarbitrary' class in QuickCheck 2 | ||
60 | coarbitrary = undefined | ||
61 | #endif | ||
62 | |||
63 | #endif | ||
64 | 40 | ||
65 | chooseDim = sized $ \m -> choose (1,max 1 m) | 41 | chooseDim = sized $ \m -> choose (1,max 1 m) |
66 | 42 | ||
@@ -68,15 +44,9 @@ instance (Field a, Arbitrary a) => Arbitrary (Vector a) where | |||
68 | arbitrary = do m <- chooseDim | 44 | arbitrary = do m <- chooseDim |
69 | l <- vector m | 45 | l <- vector m |
70 | return $ fromList l | 46 | return $ fromList l |
71 | |||
72 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
73 | -- shrink any one of the components | 47 | -- shrink any one of the components |
74 | shrink = map fromList . shrinkListElementwise . toList | 48 | shrink = map fromList . shrinkListElementwise . toList |
75 | 49 | ||
76 | #else | ||
77 | coarbitrary = undefined | ||
78 | #endif | ||
79 | |||
80 | instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where | 50 | instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where |
81 | arbitrary = do | 51 | arbitrary = do |
82 | m <- chooseDim | 52 | m <- chooseDim |
@@ -84,16 +54,11 @@ instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where | |||
84 | l <- vector (m*n) | 54 | l <- vector (m*n) |
85 | return $ (m><n) l | 55 | return $ (m><n) l |
86 | 56 | ||
87 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
88 | -- shrink any one of the components | 57 | -- shrink any one of the components |
89 | shrink a = map (rows a >< cols a) | 58 | shrink a = map (rows a >< cols a) |
90 | . shrinkListElementwise | 59 | . shrinkListElementwise |
91 | . concat . toLists | 60 | . concat . toLists |
92 | $ a | 61 | $ a |
93 | #else | ||
94 | coarbitrary = undefined | ||
95 | #endif | ||
96 | |||
97 | 62 | ||
98 | -- a square matrix | 63 | -- a square matrix |
99 | newtype (Sq a) = Sq (Matrix a) deriving Show | 64 | newtype (Sq a) = Sq (Matrix a) deriving Show |
@@ -103,11 +68,7 @@ instance (Element a, Arbitrary a) => Arbitrary (Sq a) where | |||
103 | l <- vector (n*n) | 68 | l <- vector (n*n) |
104 | return $ Sq $ (n><n) l | 69 | return $ Sq $ (n><n) l |
105 | 70 | ||
106 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
107 | shrink (Sq a) = [ Sq b | b <- shrink a ] | 71 | shrink (Sq a) = [ Sq b | b <- shrink a ] |
108 | #else | ||
109 | coarbitrary = undefined | ||
110 | #endif | ||
111 | 72 | ||
112 | 73 | ||
113 | -- a unitary matrix | 74 | -- a unitary matrix |
@@ -118,11 +79,6 @@ instance (Field a, Arbitrary a) => Arbitrary (Rot a) where | |||
118 | let (q,_) = qr m | 79 | let (q,_) = qr m |
119 | return (Rot q) | 80 | return (Rot q) |
120 | 81 | ||
121 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
122 | #else | ||
123 | coarbitrary = undefined | ||
124 | #endif | ||
125 | |||
126 | 82 | ||
127 | -- a complex hermitian or real symmetric matrix | 83 | -- a complex hermitian or real symmetric matrix |
128 | newtype (Her a) = Her (Matrix a) deriving Show | 84 | newtype (Her a) = Her (Matrix a) deriving Show |
@@ -130,12 +86,8 @@ instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where | |||
130 | arbitrary = do | 86 | arbitrary = do |
131 | Sq m <- arbitrary | 87 | Sq m <- arbitrary |
132 | let m' = m/2 | 88 | let m' = m/2 |
133 | return $ Her (m' + ctrans m') | 89 | return $ Her (m' + tr m') |
134 | 90 | ||
135 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
136 | #else | ||
137 | coarbitrary = undefined | ||
138 | #endif | ||
139 | 91 | ||
140 | class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a | 92 | class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a |
141 | instance ArbitraryField Double | 93 | instance ArbitraryField Double |
@@ -144,7 +96,7 @@ instance ArbitraryField (Complex Double) | |||
144 | 96 | ||
145 | -- a well-conditioned general matrix (the singular values are between 1 and 100) | 97 | -- a well-conditioned general matrix (the singular values are between 1 and 100) |
146 | newtype (WC a) = WC (Matrix a) deriving Show | 98 | newtype (WC a) = WC (Matrix a) deriving Show |
147 | instance (ArbitraryField a) => Arbitrary (WC a) where | 99 | instance (Numeric a, ArbitraryField a) => Arbitrary (WC a) where |
148 | arbitrary = do | 100 | arbitrary = do |
149 | m <- arbitrary | 101 | m <- arbitrary |
150 | let (u,_,v) = svd m | 102 | let (u,_,v) = svd m |
@@ -153,34 +105,24 @@ instance (ArbitraryField a) => Arbitrary (WC a) where | |||
153 | n = min r c | 105 | n = min r c |
154 | sv' <- replicateM n (choose (1,100)) | 106 | sv' <- replicateM n (choose (1,100)) |
155 | let s = diagRect 0 (fromList sv') r c | 107 | let s = diagRect 0 (fromList sv') r c |
156 | return $ WC (u `mXm` real s `mXm` trans v) | 108 | return $ WC (u <> real s <> tr v) |
157 | |||
158 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
159 | #else | ||
160 | coarbitrary = undefined | ||
161 | #endif | ||
162 | 109 | ||
163 | 110 | ||
164 | -- a well-conditioned square matrix (the singular values are between 1 and 100) | 111 | -- a well-conditioned square matrix (the singular values are between 1 and 100) |
165 | newtype (SqWC a) = SqWC (Matrix a) deriving Show | 112 | newtype (SqWC a) = SqWC (Matrix a) deriving Show |
166 | instance (ArbitraryField a) => Arbitrary (SqWC a) where | 113 | instance (ArbitraryField a, Numeric a) => Arbitrary (SqWC a) where |
167 | arbitrary = do | 114 | arbitrary = do |
168 | Sq m <- arbitrary | 115 | Sq m <- arbitrary |
169 | let (u,_,v) = svd m | 116 | let (u,_,v) = svd m |
170 | n = rows m | 117 | n = rows m |
171 | sv' <- replicateM n (choose (1,100)) | 118 | sv' <- replicateM n (choose (1,100)) |
172 | let s = diag (fromList sv') | 119 | let s = diag (fromList sv') |
173 | return $ SqWC (u `mXm` real s `mXm` trans v) | 120 | return $ SqWC (u <> real s <> tr v) |
174 | |||
175 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
176 | #else | ||
177 | coarbitrary = undefined | ||
178 | #endif | ||
179 | 121 | ||
180 | 122 | ||
181 | -- a positive definite square matrix (the eigenvalues are between 0 and 100) | 123 | -- a positive definite square matrix (the eigenvalues are between 0 and 100) |
182 | newtype (PosDef a) = PosDef (Matrix a) deriving Show | 124 | newtype (PosDef a) = PosDef (Matrix a) deriving Show |
183 | instance (ArbitraryField a, Num (Vector a)) | 125 | instance (Numeric a, ArbitraryField a, Num (Vector a)) |
184 | => Arbitrary (PosDef a) where | 126 | => Arbitrary (PosDef a) where |
185 | arbitrary = do | 127 | arbitrary = do |
186 | Her m <- arbitrary | 128 | Her m <- arbitrary |
@@ -188,13 +130,8 @@ instance (ArbitraryField a, Num (Vector a)) | |||
188 | n = rows m | 130 | n = rows m |
189 | l <- replicateM n (choose (0,100)) | 131 | l <- replicateM n (choose (0,100)) |
190 | let s = diag (fromList l) | 132 | let s = diag (fromList l) |
191 | p = v `mXm` real s `mXm` ctrans v | 133 | p = v <> real s <> tr v |
192 | return $ PosDef (0.5 * p + 0.5 * ctrans p) | 134 | return $ PosDef (0.5 * p + 0.5 * tr p) |
193 | |||
194 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
195 | #else | ||
196 | coarbitrary = undefined | ||
197 | #endif | ||
198 | 135 | ||
199 | 136 | ||
200 | -- a pair of matrices that can be multiplied | 137 | -- a pair of matrices that can be multiplied |
@@ -208,11 +145,7 @@ instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where | |||
208 | lb <- vector (k*m) | 145 | lb <- vector (k*m) |
209 | return $ Consistent ((n><k) la, (k><m) lb) | 146 | return $ Consistent ((n><k) la, (k><m) lb) |
210 | 147 | ||
211 | #if MIN_VERSION_QuickCheck(2,0,0) | ||
212 | shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ] | 148 | shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ] |
213 | #else | ||
214 | coarbitrary = undefined | ||
215 | #endif | ||
216 | 149 | ||
217 | 150 | ||
218 | 151 | ||
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 | ||