From 77552d080e88fc70312f55fd3303fac3464ab46e Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 14 Dec 2011 13:08:43 +0100 Subject: new package hmatrix-tests --- lib/Numeric/LinearAlgebra/Tests.hs | 723 --------------------- lib/Numeric/LinearAlgebra/Tests/Instances.hs | 249 ------- lib/Numeric/LinearAlgebra/Tests/Properties.hs | 272 -------- lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h | 33 - 4 files changed, 1277 deletions(-) delete mode 100644 lib/Numeric/LinearAlgebra/Tests.hs delete mode 100644 lib/Numeric/LinearAlgebra/Tests/Instances.hs delete mode 100644 lib/Numeric/LinearAlgebra/Tests/Properties.hs delete mode 100644 lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h (limited to 'lib') diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs deleted file mode 100644 index e859450..0000000 --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ /dev/null @@ -1,723 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns #-} ------------------------------------------------------------------------------ -{- | -Module : Numeric.LinearAlgebra.Tests -Copyright : (c) Alberto Ruiz 2007-9 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) -Stability : provisional -Portability : portable - -Some tests. - --} - -module Numeric.LinearAlgebra.Tests( --- module Numeric.LinearAlgebra.Tests.Instances, --- module Numeric.LinearAlgebra.Tests.Properties, - qCheck, runTests, runBenchmarks, findNaN ---, runBigTests -) where - -import Data.Packed.Random -import Numeric.LinearAlgebra -import Numeric.LinearAlgebra.LAPACK -import Numeric.LinearAlgebra.Tests.Instances -import Numeric.LinearAlgebra.Tests.Properties -import Test.HUnit hiding ((~:),test,Testable,State) -import System.Info -import Data.List(foldl1') -import Numeric.GSL -import Prelude hiding ((^)) -import qualified Prelude -import System.CPUTime -import Text.Printf -import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr) -import Control.Arrow((***)) -import Debug.Trace - -#include "Tests/quickCheckCompat.h" - -debug x = trace (show x) x - -a ^ b = a Prelude.^ (b :: Int) - -utest str b = TestCase $ assertBool str b - -a ~~ b = fromList a |~| fromList b - -feye n = flipud (ident n) :: Matrix Double - ------------------------------------------------------------ - -detTest1 = det m == 26 - && det mc == 38 :+ (-3) - && det (feye 2) == -1 - where - m = (3><3) - [ 1, 2, 3 - , 4, 5, 7 - , 2, 8, 4 :: Double - ] - mc = (3><3) - [ 1, 2, 3 - , 4, 5, 7 - , 2, 8, i - ] - -detTest2 = inv1 |~| inv2 && [det1] ~~ [det2] - where - m = complex (feye 6) - inv1 = inv m - det1 = det m - (inv2,(lda,sa)) = invlndet m - det2 = sa * exp lda - --------------------------------------------------------------------- - -polyEval cs x = foldr (\c ac->ac*x+c) 0 cs - -polySolveProp p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p)) - ---------------------------------------------------------------------- - -quad f a b = fst $ integrateQAGS 1E-9 100 f a b - --- A multiple integral can be easily defined using partial application -quad2 f a b g1 g2 = quad h a b - where h x = quad (f x) (g1 x) (g2 x) - -volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) - 0 r (const 0) (\x->sqrt (r*r-x*x)) - ---------------------------------------------------------------------- - -derivTest = abs (d (\x-> x * d (\y-> x+y) 1) 1 - 1) < 1E-10 - where d f x = fst $ derivCentral 0.01 f x - ---------------------------------------------------------------------- - --- besselTest = utest "bessel_J0_e" ( abs (r-expected) < e ) --- where (r,e) = bessel_J0_e 5.0 --- expected = -0.17759677131433830434739701 - --- exponentialTest = utest "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 ) --- where (v,e,_err) = exp_e10_e 30.0 --- expected = exp 30.0 - ---------------------------------------------------------------------- - -nd1 = (3><3) [ 1/2, 1/4, 1/4 - , 0/1, 1/2, 1/4 - , 1/2, 1/4, 1/2 :: Double] - -nd2 = (2><2) [1, 0, 1, 1:: Complex Double] - -expmTest1 = expm nd1 :~14~: (3><3) - [ 1.762110887278176 - , 0.478085470590435 - , 0.478085470590435 - , 0.104719410945666 - , 1.709751181805343 - , 0.425725765117601 - , 0.851451530235203 - , 0.530445176063267 - , 1.814470592751009 ] - -expmTest2 = expm nd2 :~15~: (2><2) - [ 2.718281828459045 - , 0.000000000000000 - , 2.718281828459045 - , 2.718281828459045 ] - ---------------------------------------------------------------------- - -minimizationTest = TestList - [ utest "minimization conjugatefr" (minim1 f df [5,7] ~~ [1,2]) - , utest "minimization nmsimplex2" (minim2 f [5,7] `elem` [24,25]) - ] - where f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30 - df [x,y] = [20*(x-1), 40*(y-2)] - minim1 g dg ini = fst $ minimizeD ConjugateFR 1E-3 30 1E-2 1E-4 g dg ini - minim2 g ini = rows $ snd $ minimize NMSimplex2 1E-2 30 [1,1] g ini - ---------------------------------------------------------------------- - -rootFindingTest = TestList [ utest "root Hybrids" (fst sol1 ~~ [1,1]) - , utest "root Newton" (rows (snd sol2) == 2) - ] - where sol1 = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5] - sol2 = rootJ Newton 1E-7 30 (rosenbrock 1 10) (jacobian 1 10) [-10,-5] - rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ] - jacobian a b [x,_y] = [ [-a , 0] - , [-2*b*x, b] ] - ---------------------------------------------------------------------- - -odeTest = utest "ode" (last (toLists sol) ~~ [-1.7588880332411019, 8.364348908711941e-2]) - where sol = odeSolveV RK8pd 1E-6 1E-6 0 (l2v $ vanderpol 10) Nothing (fromList [1,0]) ts - ts = linspace 101 (0,100) - l2v f = \t -> fromList . f t . toList - vanderpol mu _t [x,y] = [y, -x + mu * y * (1-x^2) ] - ---------------------------------------------------------------------- - -fittingTest = utest "levmar" (ok1 && ok2) - where - xs = map return [0 .. 39] - sigma = 0.1 - ys = map return $ toList $ fromList (map (head . expModel [5,0.1,1]) xs) - + scalar sigma * (randomVector 0 Gaussian 40) - dats = zip xs (zip ys (repeat sigma)) - dat = zip xs ys - - expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b] - expModelDer [a,lambda,_b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]] - - sols = fst $ fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dats [1,0,0] - sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] - - ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d - ok2 = norm2 (fromList (map fst sols) - fromList sol) < 1E-5 - ------------------------------------------------------ - -mbCholTest = utest "mbCholTest" (ok1 && ok2) where - m1 = (2><2) [2,5,5,8 :: Double] - m2 = (2><2) [3,5,5,9 :: Complex Double] - ok1 = mbCholSH m1 == Nothing - ok2 = mbCholSH m2 == Just (chol m2) - ---------------------------------------------------------------------- - -randomTestGaussian = c :~1~: snd (meanCov dat) where - a = (3><3) [1,2,3, - 2,4,0, - -2,2,1] - m = 3 |> [1,2,3] - c = a <> trans a - dat = gaussianSample 7 (10^6) m c - -randomTestUniform = c :~1~: snd (meanCov dat) where - c = diag $ 3 |> map ((/12).(^2)) [1,2,3] - dat = uniformSample 7 (10^6) [(0,1),(1,3),(3,6)] - ---------------------------------------------------------------------- - -rot :: Double -> Matrix Double -rot a = (3><3) [ c,0,s - , 0,1,0 - ,-s,0,c ] - where c = cos a - s = sin a - -rotTest = fun (10^5) :~11~: rot 5E4 - where fun n = foldl1' (<>) (map rot angles) - where angles = toList $ linspace n (0,1) - ---------------------------------------------------------------------- --- vector <= 0.6.0.2 bug discovered by Patrick Perry --- http://trac.haskell.org/vector/ticket/31 - -offsetTest = y == y' where - x = fromList [0..3 :: Double] - y = subVector 1 3 x - (f,o,n) = unsafeToForeignPtr y - y' = unsafeFromForeignPtr f o n - ---------------------------------------------------------------------- - -normsVTest = TestList [ - utest "normv2CD" $ norm2PropC v - , utest "normv2CF" $ norm2PropC (single v) -#ifndef NONORMVTEST - , utest "normv2D" $ norm2PropR x - , utest "normv2F" $ norm2PropR (single x) -#endif - , utest "normv1CD" $ norm1 v == 8 - , utest "normv1CF" $ norm1 (single v) == 8 - , utest "normv1D" $ norm1 x == 6 - , utest "normv1F" $ norm1 (single x) == 6 - - , utest "normvInfCD" $ normInf v == 5 - , utest "normvInfCF" $ normInf (single v) == 5 - , utest "normvInfD" $ normInf x == 3 - , utest "normvInfF" $ normInf (single x) == 3 - - ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double) - x = fromList [1,2,-3] :: Vector Double -#ifndef NONORMVTEST - norm2PropR a = norm2 a =~= sqrt (dot a a) -#endif - norm2PropC a = norm2 a =~= realPart (sqrt (dot a (conj a))) - a =~= b = fromList [a] |~| fromList [b] - -normsMTest = TestList [ - utest "norm2mCD" $ pnorm PNorm2 v =~= 8.86164970498005 - , utest "norm2mCF" $ pnorm PNorm2 (single v) =~= 8.86164970498005 - , utest "norm2mD" $ pnorm PNorm2 x =~= 5.96667765076216 - , utest "norm2mF" $ pnorm PNorm2 (single x) =~= 5.96667765076216 - - , utest "norm1mCD" $ pnorm PNorm1 v == 9 - , utest "norm1mCF" $ pnorm PNorm1 (single v) == 9 - , utest "norm1mD" $ pnorm PNorm1 x == 7 - , utest "norm1mF" $ pnorm PNorm1 (single x) == 7 - - , utest "normmInfCD" $ pnorm Infinity v == 12 - , utest "normmInfCF" $ pnorm Infinity (single v) == 12 - , utest "normmInfD" $ pnorm Infinity x == 8 - , utest "normmInfF" $ pnorm Infinity (single x) == 8 - - , utest "normmFroCD" $ pnorm Frobenius v =~= 8.88819441731559 - , utest "normmFroCF" $ pnorm Frobenius (single v) =~~= 8.88819441731559 - , utest "normmFroD" $ pnorm Frobenius x =~= 6.24499799839840 - , utest "normmFroF" $ pnorm Frobenius (single x) =~~= 6.24499799839840 - - ] where v = (2><2) [1,-2*i,3:+4,7] :: Matrix (Complex Double) - x = (2><2) [1,2,-3,5] :: Matrix Double - a =~= b = fromList [a] :~10~: fromList [b] - a =~~= b = fromList [a] :~5~: fromList [b] - ---------------------------------------------------------------------- - -sumprodTest = TestList [ - utest "sumCD" $ sumElements z == 6 - , utest "sumCF" $ sumElements (single z) == 6 - , utest "sumD" $ sumElements v == 6 - , utest "sumF" $ sumElements (single v) == 6 - - , utest "prodCD" $ prodProp z - , utest "prodCF" $ prodProp (single z) - , utest "prodD" $ prodProp v - , utest "prodF" $ prodProp (single v) - ] where v = fromList [1,2,3] :: Vector Double - z = fromList [1,2-i,3+i] - prodProp x = prodElements x == product (toList x) - ---------------------------------------------------------------------- - -chainTest = utest "chain" $ foldl1' (<>) ms |~| optimiseMult ms where - ms = [ diag (fromList [1,2,3 :: Double]) - , konst 3 (3,5) - , (5><10) [1 .. ] - , konst 5 (10,2) - ] - ---------------------------------------------------------------------- - -conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m) - ---------------------------------------------------------------------- - -newtype State s a = State { runState :: s -> (a,s) } - -instance Monad (State s) where - return a = State $ \s -> (a,s) - m >>= f = State $ \s -> let (a,s') = runState m s - in runState (f a) s' - -state_get :: State s s -state_get = State $ \s -> (s,s) - -state_put :: s -> State s () -state_put s = State $ \_ -> ((),s) - -evalState :: State s a -> s -> a -evalState m s = let (a,s') = runState m s - in seq s' a - -newtype MaybeT m a = MaybeT { runMaybeT :: m (Maybe a) } - -instance Monad m => Monad (MaybeT m) where - return a = MaybeT $ return $ Just a - m >>= f = MaybeT $ do - res <- runMaybeT m - case res of - Nothing -> return Nothing - Just r -> runMaybeT (f r) - fail _ = MaybeT $ return Nothing - -lift_maybe m = MaybeT $ do - res <- m - return $ Just res - --- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs ---successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool -successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (dim v - 1) v))) (v @> 0) - where stp e = do - ep <- lift_maybe $ state_get - if t e ep - then lift_maybe $ state_put e - else (fail "successive_ test failed") - --- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input ---successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b -successive f v = evalState (mapVectorM stp (subVector 1 (dim v - 1) v)) (v @> 0) - where stp e = do - ep <- state_get - state_put e - return $ f ep e - - -succTest = utest "successive" $ - successive_ (>) (fromList [1 :: Double,2,3,4]) == True - && successive_ (>) (fromList [1 :: Double,3,2,4]) == False - && successive (+) (fromList [1..10 :: Double]) == 9 |> [3,5,7,9,11,13,15,17,19] - ---------------------------------------------------------------------- - -findAssocTest = utest "findAssoc" ok - where - ok = m1 == m2 - m1 = assoc (6,6) 7 $ zip (find (>0) (ident 5 :: Matrix Float)) [10 ..] :: Matrix Double - m2 = diagRect 7 (fromList[10..14]) 6 6 - ---------------------------------------------------------------------- - -condTest = utest "cond" ok - where - ok = step v * v == cond v 0 0 0 v - v = fromList [-7 .. 7 ] :: Vector Float - ---------------------------------------------------------------------- - -conformTest = utest "conform" ok - where - ok = 1 + row [1,2,3] + col [10,20,30,40] + (4><3) [1..] - == (4><3) [13,15,17 - ,26,28,30 - ,39,41,43 - ,52,54,56] - row = asRow . fromList - col = asColumn . fromList :: [Double] -> Matrix Double - ---------------------------------------------------------------------- - -accumTest = utest "accum" ok - where - x = ident 3 :: Matrix Double - ok = accum x (+) [((1,2),7), ((2,2),3)] - == (3><3) [1,0,0 - ,0,1,7 - ,0,0,4] - && - toList (flatten x) == [1,0,0,0,1,0,0,0,1] - ---------------------------------------------------------------------- - --- | All tests must pass with a maximum dimension of about 20 --- (some tests may fail with bigger sizes due to precision loss). -runTests :: Int -- ^ maximum dimension - -> IO () -runTests n = do - setErrorHandlerOff - let test p = qCheck n p - putStrLn "------ mult Double" - test (multProp1 10 . rConsist) - test (multProp1 10 . cConsist) - test (multProp2 10 . rConsist) - test (multProp2 10 . cConsist) - putStrLn "------ mult Float" - test (multProp1 6 . (single *** single) . rConsist) - test (multProp1 6 . (single *** single) . cConsist) - test (multProp2 6 . (single *** single) . rConsist) - test (multProp2 6 . (single *** single) . cConsist) - putStrLn "------ sub-trans" - test (subProp . rM) - test (subProp . cM) - putStrLn "------ ctrans" - test (conjuTest . cM) - test (conjuTest . zM) - putStrLn "------ lu" - test (luProp . rM) - test (luProp . cM) - putStrLn "------ inv (linearSolve)" - test (invProp . rSqWC) - test (invProp . cSqWC) - putStrLn "------ luSolve" - test (linearSolveProp (luSolve.luPacked) . rSqWC) - test (linearSolveProp (luSolve.luPacked) . cSqWC) - putStrLn "------ cholSolve" - test (linearSolveProp (cholSolve.chol) . rPosDef) - test (linearSolveProp (cholSolve.chol) . cPosDef) - putStrLn "------ luSolveLS" - test (linearSolveProp linearSolveLS . rSqWC) - test (linearSolveProp linearSolveLS . cSqWC) - test (linearSolveProp2 linearSolveLS . rConsist) - test (linearSolveProp2 linearSolveLS . cConsist) - putStrLn "------ pinv (linearSolveSVD)" - test (pinvProp . rM) - test (pinvProp . cM) - putStrLn "------ det" - test (detProp . rSqWC) - test (detProp . cSqWC) - putStrLn "------ svd" - test (svdProp1 . rM) - test (svdProp1 . cM) - test (svdProp1a svdR) - test (svdProp1a svdC) - test (svdProp1a svdRd) - test (svdProp1b svdR) - test (svdProp1b svdC) - test (svdProp1b svdRd) - test (svdProp2 thinSVDR) - test (svdProp2 thinSVDC) - test (svdProp2 thinSVDRd) - test (svdProp2 thinSVDCd) - test (svdProp3 . rM) - test (svdProp3 . cM) - test (svdProp4 . rM) - test (svdProp4 . cM) - test (svdProp5a) - test (svdProp5b) - test (svdProp6a) - test (svdProp6b) - test (svdProp7 . rM) - test (svdProp7 . cM) - putStrLn "------ svdCd" -#ifdef NOZGESDD - putStrLn "Omitted" -#else - test (svdProp1a svdCd) - test (svdProp1b svdCd) -#endif - putStrLn "------ eig" - test (eigSHProp . rHer) - test (eigSHProp . cHer) - test (eigProp . rSq) - test (eigProp . cSq) - test (eigSHProp2 . rHer) - test (eigSHProp2 . cHer) - test (eigProp2 . rSq) - test (eigProp2 . cSq) - putStrLn "------ nullSpace" - test (nullspaceProp . rM) - test (nullspaceProp . cM) - putStrLn "------ qr" - test (qrProp . rM) - test (qrProp . cM) - test (rqProp . rM) - test (rqProp . cM) - test (rqProp1 . cM) - test (rqProp2 . cM) - test (rqProp3 . cM) - putStrLn "------ hess" - test (hessProp . rSq) - test (hessProp . cSq) - putStrLn "------ schur" - test (schurProp2 . rSq) - test (schurProp1 . cSq) - putStrLn "------ chol" - test (cholProp . rPosDef) - test (cholProp . cPosDef) - test (exactProp . rPosDef) - test (exactProp . cPosDef) - putStrLn "------ expm" - test (expmDiagProp . complex. rSqWC) - test (expmDiagProp . cSqWC) - putStrLn "------ fft" - test (\v -> ifft (fft v) |~| v) - putStrLn "------ vector operations - Double" - test (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM)) - test $ (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::CM)) . liftMatrix makeUnitary - test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM)) - test (\u -> cos u * tan u |~| sin (u::RM)) - test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary - putStrLn "------ vector operations - Float" - test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM)) - test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary - test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM)) - test (\u -> cos u * tan u |~~| sin (u::FM)) - test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary - putStrLn "------ read . show" - test (\m -> (m::RM) == read (show m)) - test (\m -> (m::CM) == read (show m)) - test (\m -> toRows (m::RM) == read (show (toRows m))) - test (\m -> toRows (m::CM) == read (show (toRows m))) - test (\m -> (m::FM) == read (show m)) - test (\m -> (m::ZM) == read (show m)) - test (\m -> toRows (m::FM) == read (show (toRows m))) - test (\m -> toRows (m::ZM) == read (show (toRows m))) - putStrLn "------ some unit tests" - _ <- runTestTT $ TestList - [ utest "1E5 rots" rotTest - , utest "det1" detTest1 - , utest "invlndet" detTest2 - , utest "expm1" (expmTest1) - , utest "expm2" (expmTest2) - , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM) - , utest "arith2" $ ((scalar (1+i) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*i-51) :: CM) - , utest "arith3" $ exp (scalar i * ones(10,10)*pi) + 1 |~| 0 - , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3] --- , utest "gamma" (gamma 5 == 24.0) --- , besselTest --- , exponentialTest - , utest "deriv" derivTest - , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < 1E-8) - , utest "polySolve" (polySolveProp [1,2,3,4]) - , minimizationTest - , rootFindingTest - , utest "randomGaussian" randomTestGaussian - , utest "randomUniform" randomTestUniform - , utest "buildVector/Matrix" $ - complex (10 |> [0::Double ..]) == buildVector 10 fromIntegral - && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0) - , utest "rank" $ rank ((2><3)[1,0,0,1,6*eps,0]) == 1 - && rank ((2><3)[1,0,0,1,7*eps,0]) == 2 - , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM) - , odeTest - , fittingTest - , mbCholTest - , utest "offset" offsetTest - , normsVTest - , normsMTest - , sumprodTest - , chainTest - , succTest - , findAssocTest - , condTest - , conformTest - , accumTest - ] - return () - - --- single precision approximate equality -infixl 4 |~~| -a |~~| b = a :~6~: b - -makeUnitary v | realPart n > 1 = v / scalar n - | otherwise = v - where n = sqrt (conj v <.> v) - --- -- | Some additional tests on big matrices. They take a few minutes. --- runBigTests :: IO () --- runBigTests = undefined - --- testcase for nonempty fpu stack -findNaN :: Int -> Bool -findNaN n = all (bugProp . eye) (take n $ cycle [1..20]) - where eye m = ident m :: Matrix ( Double) - --------------------------------------------------------------------------------- - --- | Performance measurements. -runBenchmarks :: IO () -runBenchmarks = do - --cholBench - solveBench - subBench - multBench - svdBench - eigBench - putStrLn "" - --------------------------------- - -time msg act = do - putStr (msg++" ") - t0 <- getCPUTime - act `seq` putStr " " - t1 <- getCPUTime - printf "%6.2f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Double)) :: IO () - return () - --------------------------------- - -manymult n = foldl1' (<>) (map rot2 angles) where - angles = toList $ linspace n (0,1) - rot2 :: Double -> Matrix Double - rot2 a = (3><3) [ c,0,s - , 0,1,0 - ,-s,0,c ] - where c = cos a - s = sin a - -multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Double)) - --------------------------------- - -subBench = do - putStrLn "" - let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v)) - time "0.1M subVector " (g (constant 1 (1+10^5) :: Vector Double) @> 0) - let f = foldl1' (.) (replicate (10^5) (fromRows.toRows)) - time "subVector-join 3" (f (ident 3 :: Matrix Double) @@>(0,0)) - time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0)) - --------------------------------- - -multBench = do - let a = ident 1000 :: Matrix Double - let b = ident 2000 :: Matrix Double - a `seq` b `seq` putStrLn "" - time "product of 1M different 3x3 matrices" (manymult (10^6)) - putStrLn "" - time "product of 1M constant 1x1 matrices" (multb 1) - time "product of 1M constant 3x3 matrices" (multb 3) - --time "product of 1M constant 5x5 matrices" (multb 5) - time "product of 1M const. 10x10 matrices" (multb 10) - --time "product of 1M const. 15x15 matrices" (multb 15) - time "product of 1M const. 20x20 matrices" (multb 20) - --time "product of 1M const. 25x25 matrices" (multb 25) - putStrLn "" - time "product (1000 x 1000)<>(1000 x 1000)" (a<>a) - time "product (2000 x 2000)<>(2000 x 2000)" (b<>b) - --------------------------------- - -eigBench = do - let m = reshape 1000 (randomVector 777 Uniform (1000*1000)) - s = m + trans m - m `seq` s `seq` putStrLn "" - time "eigenvalues symmetric 1000x1000" (eigenvaluesSH' m) - time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m) - time "eigenvalues general 1000x1000" (eigenvalues m) - time "eigenvectors general 1000x1000" (snd $ eig m) - --------------------------------- - -svdBench = do - let a = reshape 500 (randomVector 777 Uniform (3000*500)) - b = reshape 1000 (randomVector 777 Uniform (1000*1000)) - fv (_,_,v) = v@@>(0,0) - a `seq` b `seq` putStrLn "" - time "singular values 3000x500" (singularValues a) - time "thin svd 3000x500" (fv $ thinSVD a) - time "full svd 3000x500" (fv $ svd a) - time "singular values 1000x1000" (singularValues b) - time "full svd 1000x1000" (fv $ svd b) - --------------------------------- - -solveBenchN n = do - let x = uniformSample 777 (2*n) (replicate n (-1,1)) - a = trans x <> x - b = asColumn $ randomVector 666 Uniform n - a `seq` b `seq` putStrLn "" - time ("svd solve " ++ show n) (linearSolveSVD a b) - time (" ls solve " ++ show n) (linearSolveLS a b) - time (" solve " ++ show n) (linearSolve a b) - time ("cholSolve " ++ show n) (cholSolve (chol a) b) - -solveBench = do - solveBenchN 500 - solveBenchN 1000 - -- solveBenchN 1500 - --------------------------------- - -cholBenchN n = do - let x = uniformSample 777 (2*n) (replicate n (-1,1)) - a = trans x <> x - a `seq` putStrLn "" - time ("chol " ++ show n) (chol a) - -cholBench = do - cholBenchN 1200 - cholBenchN 600 - cholBenchN 300 --- cholBenchN 150 --- cholBenchN 50 diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs deleted file mode 100644 index 6dd9cfe..0000000 --- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs +++ /dev/null @@ -1,249 +0,0 @@ -{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-} -{-# OPTIONS_GHC -fno-warn-unused-imports #-} ------------------------------------------------------------------------------ -{- | -Module : Numeric.LinearAlgebra.Tests.Instances -Copyright : (c) Alberto Ruiz 2008 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) -Stability : provisional -Portability : portable - -Arbitrary instances for vectors, matrices. - --} - -module Numeric.LinearAlgebra.Tests.Instances( - Sq(..), rSq,cSq, - Rot(..), rRot,cRot, - Her(..), rHer,cHer, - WC(..), rWC,cWC, - SqWC(..), rSqWC, cSqWC, - PosDef(..), rPosDef, cPosDef, - Consistent(..), rConsist, cConsist, - RM,CM, rM,cM, - FM,ZM, fM,zM -) where - -import System.Random - -import Numeric.LinearAlgebra -import Control.Monad(replicateM) -#include "quickCheckCompat.h" - -#if MIN_VERSION_QuickCheck(2,0,0) -shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]] -shrinkListElementwise [] = [] -shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ] - ++ [ x:ys | ys <- shrinkListElementwise xs ] - -shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)] -shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ] -#endif - -#if MIN_VERSION_QuickCheck(2,1,1) -#else -instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where - arbitrary = do - re <- arbitrary - im <- arbitrary - return (re :+ im) - -#if MIN_VERSION_QuickCheck(2,0,0) - shrink (re :+ im) = - [ u :+ v | (u,v) <- shrinkPair (re,im) ] -#else - -- this has been moved to the 'Coarbitrary' class in QuickCheck 2 - coarbitrary = undefined -#endif - -#endif - -chooseDim = sized $ \m -> choose (1,max 1 m) - -instance (Field a, Arbitrary a) => Arbitrary (Vector a) where - arbitrary = do m <- chooseDim - l <- vector m - return $ fromList l - -#if MIN_VERSION_QuickCheck(2,0,0) - -- shrink any one of the components - shrink = map fromList . shrinkListElementwise . toList - -#else - coarbitrary = undefined -#endif - -instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where - arbitrary = do - m <- chooseDim - n <- chooseDim - l <- vector (m*n) - return $ (m>< cols a) - . shrinkListElementwise - . concat . toLists - $ a -#else - coarbitrary = undefined -#endif - - --- a square matrix -newtype (Sq a) = Sq (Matrix a) deriving Show -instance (Element a, Arbitrary a) => Arbitrary (Sq a) where - arbitrary = do - n <- chooseDim - l <- vector (n*n) - return $ Sq $ (n> Arbitrary (Rot a) where - arbitrary = do - Sq m <- arbitrary - let (q,_) = qr m - return (Rot q) - -#if MIN_VERSION_QuickCheck(2,0,0) -#else - coarbitrary = undefined -#endif - - --- a complex hermitian or real symmetric matrix -newtype (Her a) = Her (Matrix a) deriving Show -instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where - arbitrary = do - Sq m <- arbitrary - let m' = m/2 - return $ Her (m' + ctrans m') - -#if MIN_VERSION_QuickCheck(2,0,0) -#else - coarbitrary = undefined -#endif - -class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a -instance ArbitraryField Double -instance ArbitraryField (Complex Double) - - --- a well-conditioned general matrix (the singular values are between 1 and 100) -newtype (WC a) = WC (Matrix a) deriving Show -instance (ArbitraryField a) => Arbitrary (WC a) where - arbitrary = do - m <- arbitrary - let (u,_,v) = svd m - r = rows m - c = cols m - n = min r c - sv' <- replicateM n (choose (1,100)) - let s = diagRect 0 (fromList sv') r c - return $ WC (u <> real s <> trans v) - -#if MIN_VERSION_QuickCheck(2,0,0) -#else - coarbitrary = undefined -#endif - - --- a well-conditioned square matrix (the singular values are between 1 and 100) -newtype (SqWC a) = SqWC (Matrix a) deriving Show -instance (ArbitraryField a) => Arbitrary (SqWC a) where - arbitrary = do - Sq m <- arbitrary - let (u,_,v) = svd m - n = rows m - sv' <- replicateM n (choose (1,100)) - let s = diag (fromList sv') - return $ SqWC (u <> real s <> trans v) - -#if MIN_VERSION_QuickCheck(2,0,0) -#else - coarbitrary = undefined -#endif - - --- a positive definite square matrix (the eigenvalues are between 0 and 100) -newtype (PosDef a) = PosDef (Matrix a) deriving Show -instance (ArbitraryField a, Num (Vector a)) - => Arbitrary (PosDef a) where - arbitrary = do - Her m <- arbitrary - let (_,v) = eigSH m - n = rows m - l <- replicateM n (choose (0,100)) - let s = diag (fromList l) - p = v <> real s <> ctrans v - return $ PosDef (0.5 * p + 0.5 * ctrans p) - -#if MIN_VERSION_QuickCheck(2,0,0) -#else - coarbitrary = undefined -#endif - - --- a pair of matrices that can be multiplied -newtype (Consistent a) = Consistent (Matrix a, Matrix a) deriving Show -instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where - arbitrary = do - n <- chooseDim - k <- chooseDim - m <- chooseDim - la <- vector (n*k) - lb <- vector (k*m) - return $ Consistent ((n> c t -> c t -> Double -dist a b = realToFrac r - where norm = pnorm Infinity - na = norm a - nb = norm b - nab = norm (a-b) - mx = max na nb - mn = min na nb - r = if mn < peps - then mx - else nab/mx - -infixl 4 |~| -a |~| b = a :~10~: b ---a |~| b = dist a b < 10^^(-10) - -data Aprox a = (:~) a Int --- (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool -a :~n~: b = dist a b < 10^^(-n) - ------------------------------------------------------- - -square m = rows m == cols m - --- orthonormal columns -orthonormal m = ctrans m <> m |~| ident (cols m) - -unitary m = square m && orthonormal m - -hermitian m = square m && m |~| ctrans m - -wellCond m = rcond m > 1/100 - -positiveDefinite m = minimum (toList e) > 0 - where (e,_v) = eigSH m - -upperTriang m = rows m == 1 || down == z - where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m)) - z = constant 0 (dim down) - -upperHessenberg m = rows m < 3 || down == z - where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m)) - z = constant 0 (dim down) - -zeros (r,c) = reshape c (constant 0 (r*c)) - -ones (r,c) = zeros (r,c) + 1 - ------------------------------------------------------ - -luProp m = m |~| p <> l <> u && f (det p) |~| f s - where (l,u,p,s) = lu m - f x = fromList [x] - -invProp m = m <> inv m |~| ident (rows m) - -pinvProp m = m <> p <> m |~| m - && p <> m <> p |~| p - && hermitian (m<>p) - && hermitian (p<>m) - where p = pinv m - -detProp m = s d1 |~| s d2 - where d1 = det m - d2 = det' * det q - det' = product $ toList $ takeDiag r - (q,r) = qr m - s x = fromList [x] - -nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c) - && orthonormal (fromColumns nl)) - where nl = nullspacePrec 1 m - n = fromColumns nl - r = rows m - c = cols m - rank m - ------------------------------------------------------------------- - --- testcase for nonempty fpu stack --- uncommenting unitary' signature eliminates the problem -bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v - where (u,d,v) = fullSVD m - -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool - unitary' a = unitary a - ------------------------------------------------------------------- - --- fullSVD -svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v - where (u,d,v) = fullSVD m - -svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where - (u,s,v) = svdfun m - d = diagRect 0 s (rows m) (cols m) - -svdProp1b svdfun m = unitary u && unitary v where - (u,_,v) = svdfun m - --- thinSVD -svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m) - where (u,s,v) = thinSVDfun m - --- compactSVD -svdProp3 m = (m |~| u <> real (diag s) <> trans v - && orthonormal u && orthonormal v) - where (u,s,v) = compactSVD m - -svdProp4 m' = m |~| u <> real (diag s) <> trans v - && orthonormal u && orthonormal v - && (dim s == r || r == 0 && dim s == 1) - where (u,s,v) = compactSVD m - m = fromBlocks [[m'],[m']] - r = rank m' - -svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where - s1 = svR m - s2 = svRd m - (_,s3,_) = svdR m - (_,s4,_) = svdRd m - (_,s5,_) = thinSVDR m - (_,s6,_) = thinSVDRd m - -svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where - s1 = svC m - s2 = svCd m - (_,s3,_) = svdC m - (_,s4,_) = svdCd m - (_,s5,_) = thinSVDC m - (_,s6,_) = thinSVDCd m - -svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' - where (u,s,v) = svdR m - (s',v') = rightSVR m - (u',s'') = leftSVR m - -svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' - where (u,s,v) = svdC m - (s',v') = rightSVC m - (u',s'') = leftSVC m - -svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s''' - where (u,s,v) = svd m - (s',v') = rightSV m - (u',_s'') = leftSV m - s''' = singularValues m - ------------------------------------------------------------------- - -eigProp m = complex m <> v |~| v <> diag s - where (s, v) = eig m - -eigSHProp m = m <> v |~| v <> real (diag s) - && unitary v - && m |~| v <> real (diag s) <> ctrans v - where (s, v) = eigSH m - -eigProp2 m = fst (eig m) |~| eigenvalues m - -eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m - ------------------------------------------------------------------- - -qrProp m = q <> r |~| m && unitary q && upperTriang r - where (q,r) = qr m - -rqProp m = r <> q |~| m && unitary q && upperTriang' r - where (r,q) = rq m - -rqProp1 m = r <> q |~| m - where (r,q) = rq m - -rqProp2 m = unitary q - where (r,q) = rq m - -rqProp3 m = upperTriang' r - where (r,q) = rq m - -upperTriang' r = upptr (rows r) (cols r) * r |~| r - where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1 - where t = f-c - -hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h - where (p,h) = hess m - -schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s - where (u,s) = schur m - -schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme - where (u,s) = schur m - -cholProp m = m |~| ctrans c <> c && upperTriang c - where c = chol m - -exactProp m = chol m == chol (m+0) - -expmDiagProp m = expm (logm m) :~ 7 ~: complex m - where logm = matFunc log - --- reference multiply -mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ] - where doth u v = sum $ zipWith (*) (toList u) (toList v) - -multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) - -multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a) - -linearSolveProp f m = f m m |~| ident (rows m) - -linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b) - where q = min (rows a) (cols a) - b = a <> x - wc = rank a == q - -subProp m = m == (trans . fromColumns . toRows) m - diff --git a/lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h b/lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h deleted file mode 100644 index 714587b..0000000 --- a/lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h +++ /dev/null @@ -1,33 +0,0 @@ -#ifndef MIN_VERSION_QuickCheck -#define MIN_VERSION_QuickCheck(A,B,C) 1 -#endif - -#if MIN_VERSION_QuickCheck(2,0,0) -import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector - ,sized,classify,Testable,Property - - ,quickCheckWith,maxSize,stdArgs,shrink) - -#else -import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector - ,sized,classify,Testable,Property - - ,check,configSize,defaultConfig,trivial) -#endif - - - -#if MIN_VERSION_QuickCheck(2,0,0) -trivial :: Testable a => Bool -> a -> Property -trivial = (`classify` "trivial") -#else -#endif - - --- define qCheck, which used to be in Tests.hs -#if MIN_VERSION_QuickCheck(2,0,0) -qCheck n = quickCheckWith stdArgs {maxSize = n} -#else -qCheck n = check defaultConfig {configSize = const n} -#endif - -- cgit v1.2.3