module Numeric.LinearAlgebra.Tests(
qCheck, runTests
) where
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Tests.Instances
import Numeric.LinearAlgebra.Tests.Properties
import Test.HUnit hiding ((~:),test,Testable)
import System.Info
import Data.List(foldl1')
import Numeric.GSL hiding (sin,cos,exp,choose)
import Prelude hiding ((^))
import qualified Prelude
#include "Tests/quickCheckCompat.h"
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
]
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
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*rx*xy*y))
0 r (const 0) (\x->sqrt (r*rx*x))
besselTest = utest "bessel_J0_e" ( abs (rexpected) < 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*(x1)^2 + 20*(y2)^2 + 30
df [x,y] = [20*(x1), 40*(y2)]
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*(1x), b*(yx^2) ]
jacobian a b [x,_y] = [ [a , 0]
, [2*b*x, b] ]
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) :~12~: rot 5E4
where fun n = foldl1' (<>) (map rot angles)
where angles = toList $ linspace n (0,1)
runTests :: Int
-> IO ()
runTests n = do
setErrorHandlerOff
let test p = qCheck n p
putStrLn "------ mult"
test (multProp1 . rConsist)
test (multProp1 . cConsist)
test (multProp2 . rConsist)
test (multProp2 . cConsist)
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 "------ 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 (svdProp2 . rM)
test (svdProp2 . cM)
putStrLn "------ eig"
test (eigSHProp . rHer)
test (eigSHProp . cHer)
test (eigProp . rSq)
test (eigProp . cSq)
putStrLn "------ nullSpace"
test (nullspaceProp . rM)
test (nullspaceProp . cM)
putStrLn "------ qr"
test (qrProp . rM)
test (qrProp . 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)
putStrLn "------ expm"
test (expmDiagProp . rSqWC)
test (expmDiagProp . cSqWC)
putStrLn "------ fft"
test (\v -> ifft (fft v) |~| v)
putStrLn "------ vector operations"
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 "------ 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)))
putStrLn "------ some unit tests"
runTestTT $ TestList
[ utest "1E5 rots" rotTest
, utest "det1" detTest1
, utest "expm1" (expmTest1)
, utest "expm2" (expmTest2)
, utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 7)**2 |~| (49 :: RM)
, utest "arith2" $ (((1+i) .* ones (100,100) * 5 + 2)/0.5 7)**2 |~| ( (140*i51).*1 :: CM)
, utest "arith3" $ exp (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 "integrate" (abs (volSphere 2.5 4/3*pi*2.5^3) < 1E-8)
, utest "polySolve" (polySolveProp [1,2,3,4])
, minimizationTest
, rootFindingTest
]
return ()
makeUnitary v | realPart n > 1 = v */ n
| otherwise = v
where n = sqrt (conj v <.> v)