-- -- QuickCheck tests -- ----------------------------------------------------------------------------- import Data.Packed.Internal import Data.Packed.Vector import Data.Packed.Matrix import Data.Packed.Internal.Matrix import GSL.Vector import GSL.Integration import GSL.Differentiation import GSL.Special import GSL.Fourier import GSL.Polynomials import LAPACK import Test.QuickCheck import Test.HUnit import Complex {- -- Bravo por quickCheck! pinvProp1 tol m = (rank m == cols m) ==> pinv m <> m ~~ ident (cols m) where infix 2 ~~ (~~) = approxEqual tol pinvProp2 tol m = 0 < r && r <= c ==> (r==c) `trivial` (m <> pinv m <> m ~~ m) where r = rank m c = cols m infix 2 ~~ (~~) = approxEqual tol nullspaceProp tol m = cr > 0 ==> m <> nt ~~ zeros where nt = trans (nullspace m) cr = corank m r = rows m zeros = create [r,cr] $ replicate (r*cr) 0 -} ac = (2><3) [1 .. 6::Double] bc = (3><4) [7 .. 18::Double] mz = (2 >< 3) [1,2,3,4,5,6:+(1::Double)] af = (2>|<3) [1,4,2,5,3,6::Double] bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double] a |=| b = rows a == rows b && cols a == cols b && toList (cdat a) == toList (cdat b) && toList (fdat a) == toList (fdat b) aprox fun a b = rows a == rows b && cols a == cols b && eps > aproxL fun (toList (t a)) (toList (t b)) where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat aproxL fun v1 v2 = sum (zipWith (\a b-> fun (a-b)) v1 v2) / fromIntegral (length v1) normVR a b = toScalarR AbsSum (vectorZipR Sub a b) a |~| b = rows a == rows b && cols a == cols b && eps > normVR (t a) (t b) where t = if (order a == RowMajor) `xor` isTrans a then cdat else fdat (|~~|) = aprox magnitude v1 ~~ v2 = reshape 1 v1 |~~| reshape 1 v2 u ~|~ v = normVR u v < eps eps = 1E-8::Double asFortran m = (rows m >|< cols m) $ toList (fdat m) asC m = (rows m >< cols m) $ toList (cdat m) mulC a b = multiply RowMajor a b mulF a b = multiply ColumnMajor a b cc = mulC ac bf cf = mulF af bc r = mulC cc (trans cf) rd = (2><2) [ 27736.0, 65356.0 , 65356.0, 154006.0 ::Double] instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where arbitrary = do r <- arbitrary i <- arbitrary return (r:+i) coarbitrary = undefined instance (Field a, Arbitrary a) => Arbitrary (Matrix a) where arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) m <- choose (1,10) n <- choose (1,10) l <- vector (m*n) ctype <- arbitrary let h = if ctype then (m>| Arbitrary (PairM a) where arbitrary = do a <- choose (1,10) b <- choose (1,10) c <- choose (1,10) l1 <- vector (a*b) l2 <- vector (b*c) return $ PairM ((a> Arbitrary (SqM a) where arbitrary = do n <- choose (1,10) l <- vector (n*n) return $ SqM $ (n> Arbitrary (Sym a) where arbitrary = do SqM m <- arbitrary return $ Sym (m `addM` trans m) coarbitrary = undefined data Her = Her (Matrix (Complex Double)) deriving Show instance {-(Field a, Arbitrary a, Num a) =>-} Arbitrary Her where arbitrary = do SqM m <- arbitrary return $ Her (m `addM` (liftMatrix conj) (trans m)) coarbitrary = undefined data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show instance (Num a, Field a, Arbitrary a) => Arbitrary (PairSM a) where arbitrary = do a <- choose (1,10) c <- choose (1,10) l1 <- vector (a*a) l2 <- vector (a*c) return $ PairSM ((a> Arbitrary (Vector a) where arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) m <- choose (1,100) l <- vector m return $ fromList l coarbitrary = undefined data PairV a = PairV (Vector a) (Vector a) instance (Field a, Arbitrary a) => Arbitrary (PairV a) where arbitrary = do --m <- sized $ \max -> choose (1,1+3*max) m <- choose (1,100) l1 <- vector m l2 <- vector m return $ PairV (fromList l1) (fromList l2) coarbitrary = undefined addM m1 m2 = liftMatrix2 add m1 m2 type BaseType = Double svdTestR fun prod m = u <> s <> trans v |~| m && u <> trans u |~| ident (rows m) && v <> trans v |~| ident (cols m) where (u,s,v) = fun m (<>) = prod svdTestC prod m = u <> s' <> (trans v) |~~| m && u <> (liftMatrix conj) (trans u) |~~| ident (rows m) && v <> (liftMatrix conj) (trans v) |~~| ident (cols m) where (u,s,v) = svdC m (<>) = prod s' = liftMatrix comp s eigTestC prod (SqM m) = (m <> v) |~~| (v <> diag s) && takeDiag ((liftMatrix conj (trans v)) <> v) ~~ constant (rows m) 1 --normalized where (s,v) = eigC m (<>) = prod eigTestR prod (SqM m) = (liftMatrix comp m <> v) |~~| (v <> diag s) -- && takeDiag ((liftMatrix conj (trans v)) <> v) ~~ constant (rows m) 1 --normalized ??? where (s,v) = eigR m (<>) = prod eigTestS prod (Sym m) = (m <> v) |~| (v <> diag s) && v <> trans v |~| ident (cols m) where (s,v) = eigS m (<>) = prod eigTestH prod (Her m) = (m <> v) |~~| (v <> diag (comp s)) && v <> (liftMatrix conj) (trans v) |~~| ident (cols m) where (s,v) = eigH m (<>) = prod linearSolveSQTest fun eqfun singu prod (PairSM a b) = singu a || (a <> fun a b) ==== b where (<>) = prod (====) = eqfun prec = 1E-15 singular fun m = s1 < prec || s2/s1 < prec where (_,ss,v) = fun m s = toList ss s1 = maximum s s2 = minimum s {- invTest msg m = do assertBool msg $ m <> inv m =~= ident (rows m) invComplexTest msg m = do assertBool msg $ m <> invC m =~= identC (rows m) invC m = linearSolveC m (identC (rows m)) identC n = toComplex(ident n, (0::Double) <>ident n) -} -------------------------------------------------------------------- pinvTest f feq m = (m <> f m <> m) `feq` m where (<>) = mulF pinvR m = linearSolveLSR m (ident (rows m)) pinvC m = linearSolveLSC m (ident (rows m)) pinvSVDR m = linearSolveSVDR Nothing m (ident (rows m)) pinvSVDC m = linearSolveSVDC Nothing m (ident (rows m)) -------------------------------------------------------------------- polyEval cs x = foldr (\c ac->ac*x+c) 0 cs polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p)) where l1 |~~| l2 = eps > aproxL magnitude l1 l2 polySolveTest = assertBool "polySolve" (polySolveTest' [1,2,3,4]) --------------------------------------------------------------------- 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)) integrateTest = assertBool "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < eps) --------------------------------------------------------------------- arit1 u = vectorMapValR PowVS 2 (vectorMapR Sin u) `add` vectorMapValR PowVS 2 (vectorMapR Cos u) ~|~ constant (dim u) 1 arit2 u = (vectorMapR Cos u) `mul` (vectorMapR Tan u) ~|~ vectorMapR Sin u -- arit3 (PairV u v) = --------------------------------------------------------------------- besselTest = do let (r,e) = bessel_J0_e 5.0 let expected = -0.17759677131433830434739701 assertBool "bessel_J0_e" ( abs (r-expected) < e ) exponentialTest = do let (v,e,err) = exp_e10_e 30.0 let expected = exp 30.0 assertBool "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 ) tests = TestList [ TestCase $ besselTest , TestCase $ exponentialTest , TestCase $ polySolveTest , TestCase $ integrateTest ] ---------------------------------------------------------------------- main = do putStrLn "--------- general -----" quickCheck (\(Sym m) -> m |=| (trans m:: Matrix BaseType)) quickCheck $ \l -> null l || (toList . fromList) l == (l :: [BaseType]) quickCheck $ \m -> m |=| asC (m :: Matrix BaseType) quickCheck $ \m -> m |=| asFortran (m :: Matrix BaseType) quickCheck $ \m -> m |=| (asC . asFortran) (m :: Matrix BaseType) putStrLn "--------- MULTIPLY ----" quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| mulF m1 (m2 :: Matrix BaseType) quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| trans (mulF (trans m2) (trans m1 :: Matrix BaseType)) quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| multiplyG m1 (m2 :: Matrix BaseType) putStrLn "--------- SVD ---------" quickCheck (svdTestR svdR mulC) quickCheck (svdTestR svdR mulF) quickCheck (svdTestR svdRdd mulC) quickCheck (svdTestR svdRdd mulF) quickCheck (svdTestC mulC) quickCheck (svdTestC mulF) putStrLn "--------- EIG ---------" quickCheck (eigTestC mulC) quickCheck (eigTestC mulF) quickCheck (eigTestR mulC) quickCheck (eigTestR mulF) quickCheck (eigTestS mulC) quickCheck (eigTestS mulF) quickCheck (eigTestH mulC) quickCheck (eigTestH mulF) putStrLn "--------- SOLVE ---------" quickCheck (linearSolveSQTest linearSolveR (|~|) (singular svdR') mulC) quickCheck (linearSolveSQTest linearSolveC (|~~|) (singular svdC') mulF) quickCheck (pinvTest pinvR (|~|)) quickCheck (pinvTest pinvC (|~~|)) quickCheck (pinvTest pinvSVDR (|~|)) quickCheck (pinvTest pinvSVDC (|~~|)) putStrLn "--------- VEC OPER ------" quickCheck arit1 quickCheck arit2 putStrLn "--------- GSL ------" runTestTT tests quickCheck $ \v -> ifft (fft v) ~~ v kk = (2><2) [ 1.0, 0.0 , -1.5, 1.0 ::Double] v = 11 # [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0::Double] pol = [14.125,-7.666666666666667,-14.3,-13.0,-7.0,-9.6,4.666666666666666,13.0,0.5]