diff --git a/CHANGES b/CHANGES new file mode 100644 --- /dev/null +++ b/CHANGES @@ -0,0 +1,5 @@ +0.1 +=== + +Created a separate testing package. + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,2 @@ +Copyright Alberto Ruiz 2010 +GPL license diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,5 @@ +#! /usr/bin/env runhaskell + +> import Distribution.Simple +> main = defaultMain + diff --git a/hmatrix-tests.cabal b/hmatrix-tests.cabal new file mode 100644 --- /dev/null +++ b/hmatrix-tests.cabal @@ -0,0 +1,45 @@ +Name: hmatrix-tests +Version: 0.1.0.0 +License: GPL +License-file: LICENSE +Author: Alberto Ruiz +Maintainer: Alberto Ruiz +Stability: provisional +Homepage: http://perception.inf.um.es/hmatrix +Synopsis: Tests for hmatrix +Description: Tests for hmatrix +Category: Math +tested-with: GHC==7.0.4 + +cabal-version: >=1.8 + +build-type: Simple + +extra-source-files: CHANGES + src/tests.hs + +library + + Build-Depends: base >= 4 && < 5, + hmatrix >= 0.13, + QuickCheck >= 2, HUnit, random + + hs-source-dirs: src + + exposed-modules: Numeric.LinearAlgebra.Tests + + other-modules: Numeric.LinearAlgebra.Tests.Instances, + Numeric.LinearAlgebra.Tests.Properties + + ghc-options: -Wall -fno-warn-missing-signatures -fno-warn-orphans + + +source-repository head + type: git + location: https://github.com/AlbertoRuiz/hmatrix + +Test-Suite basic + Build-Depends: base, hmatrix-tests + type: exitcode-stdio-1.0 + main-is: src/tests.hs + diff --git a/src/Numeric/LinearAlgebra/Tests.hs b/src/Numeric/LinearAlgebra/Tests.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/LinearAlgebra/Tests.hs @@ -0,0 +1,738 @@ +{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.LinearAlgebra.Tests +Copyright : (c) Alberto Ruiz 2007-11 +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 System.Exit +import Text.Printf +import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr) +import Control.Arrow((***)) +import Debug.Trace +import Control.Monad(when) + +import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector + ,sized,classify,Testable,Property + ,quickCheckWithResult,maxSize,stdArgs,shrink) + +import Test.QuickCheck.Test(isSuccess) + +qCheck n x = do + r <- quickCheckWithResult stdArgs {maxSize = n} x + when (not $ isSuccess r) (exitFailure) + +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" + c <- 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 + ] + when (errors c + failures c > 0) exitFailure + 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 + solveBench + subBench + multBench + cholBench + 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` putStr "" + time ("chol " ++ show n) (chol a) + +cholBench = do + putStrLn "" + cholBenchN 1200 + cholBenchN 600 + cholBenchN 300 +-- cholBenchN 150 +-- cholBenchN 50 diff --git a/src/Numeric/LinearAlgebra/Tests/Instances.hs b/src/Numeric/LinearAlgebra/Tests/Instances.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/LinearAlgebra/Tests/Instances.hs @@ -0,0 +1,251 @@ +{-# 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) +import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector + ,sized,classify,Testable,Property + ,quickCheckWith,maxSize,stdArgs,shrink) + +#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> Bool -> a -> Property +trivial = (`classify` "trivial") + + +-- relative error +dist :: (Normed c t, Num (c t)) => 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/src/tests.hs b/src/tests.hs new file mode 100644 --- /dev/null +++ b/src/tests.hs @@ -0,0 +1,3 @@ +import Numeric.LinearAlgebra.Tests + +main = runTests 20