cl3-1.0.0.2: tests/TestCl3.hs
{-# LANGUAGE ViewPatterns #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}
-------------------------------------------------------------------
-- |
-- Copyright : (c) 2017 Nathan Waivio
-- License : BSD3
-- Maintainer : Nathan Waivio <nathan.waivio@gmail.com>
--
-- A program to test Algebra.Geometric.Cl3
-- The code runs tests on some standard test input and then
-- runs quckcheck for some trig identities.
--
-------------------------------------------------------------------
module Main (main) where
import Test.QuickCheck (Arbitrary, arbitrary, oneof, suchThat, quickCheckWith, stdArgs, maxSuccess)
import Algebra.Geometric.Cl3
import Control.Applicative ((<*>), (<$>))
------------------------------------------------------------------
-- |
-- This program verifies the approximate equality of various trig
-- identities to the with the following limitations:
--
-- * The magnitude of the cliffor is limited in some cases.
--
-- * The imaginary part of the eigenvalues are unwrapped, due to the cyclical nature of some of the results, in a few cases.
--
-- * The poles of the functions are excluded.
--
-- * The poles of the derivatives of the functions are excluded when the cliffor is has a nilpotent component.
--
-- * Approximate equivalence is tested due to limitations with respect to floating point math.
--
-- * The implementation of Arbitrary for Cl3 limits the arbitrary cliffor such that the absolute value of cliff is less than 15
--
-- The following properties are verified in this module:
--
-- * log.exp Identity
--
-- * exp.log Identity
--
-- * abs*signum law
--
-- * The definition of recip
--
-- * recip.recip Identity
--
-- * sin.asin Identity
--
-- * asin.sin Identity
--
-- * cos.acos Identity
--
-- * acos.cos Identity
--
-- * sinh.asinh Identity
--
-- * asinh.sinh Identity
--
-- * cosh.acosh Identity
--
-- * acosh.cosh Identity
--
-- * Double Sin Identity
--
-- * Double Cos Identity
--
-- * Double Tan Identity
--
-- * Double Sinh Identity
--
-- * Double Cosh Identity
--
-- * Double Tanh Identity
--
-- * Positive Sin Shift Identity
--
-- * Negative Sin Shift Identity
--
-- * sin^2+cos^2 Identity
--
-- * cosh^2-sinh^2 Identity
--
-- * Symmetry of Cosh
--
-- * Symmetry of Sinh
--
-- * Double I Sin
--
-- * Composition Algebra Tests
--
-------------------------------------------------------------------
main :: IO ()
main = do moduleTests
print "Testing log.exp Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_LogExp
print "Testing exp.log Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_ExpLog
print "Testing abs*signum law:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_AbsSignum
print "Testing the definition of recip:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_RecipDef
print "Testing recip.recip Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_RecipID
print "Testing sin.asin Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_SinAsin
print "Testing asin.sin Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_AsinSin
print "Testing cos.acos Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_CosAcos
print "Testing acos.cos Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_AcosCos
print "Testing sinh.asinh Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_SinhAsinh
print "Testing asinh.sinh Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_AsinhSinh
print "Testing cosh.acosh Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_CoshAcosh
print "Testing acosh.cosh Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_AcoshCosh
print "Testing acosh.cosh Identity2:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_AcoshCosh2
print "Testing Double Sin Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubSin
print "Testing Double Cos Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubCos
print "Testing Double Tan Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubTan
print "Testing Double Sinh Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubSinh
print "Testing Double Cosh Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubCosh
print "Testing Double Tanh Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubTanh
print "Testing Positive Sin Shift Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_PosSinShift
print "Testing Negative Sin Shift Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_NegSinShift
print "Testing sin^2+cos^2 Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_SinSqCosSq
print "Testing cosh^2-sinh^2 Identity:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_CoshSqmSinhSq
print "Testing Symmetry of Cosh:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_SymCosh
print "Testing Symmetry of Sinh:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_SymSinh
print "Testing Double I Sin:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_DoubleISin
print "Is has Composition Sub-Algebras:"
quickCheckWith stdArgs { maxSuccess = 30000 } prop_CompAlg
----------------------------------------------------------
-- |Start of Module Tests
moduleTests :: IO ()
moduleTests = sequence_ $ tests <*> inputs
inputs :: [Cl3]
inputs = [R 0
,APS 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
,PV 0.5 0 0 0.5
,PV 0.5 0 0 (-0.5)
,BPV 0.5 0 0 0 (-0.5) 0
,BPV 0.5 0 0 0 0.5 0
,R 1
,V3 1 0 0
,APS 1 0.5 0 0 0 0.5 0 0
,APS 1 0.5 0 0 0 (-0.5) 0 0
,PV 1 1 0 0
,R 1
,R (-1)
,V3 1 0 0
,V3 (-1) 0 0
,V3 0 1 0
,V3 0 (-1) 0
,V3 0 0 1
,V3 0 0 (-1)
,BV 1 0 0
,BV (-1) 0 0
,BV 0 1 0
,BV 0 (-1) 0
,BV 0 0 1
,BV 0 0 (-1)
,I 1
,I (-1)
]
-- | 'tests' is a list of tests
-- The out of bounds can be the poles of the function or if the cliffor has
-- nilpotent content then the poles of the derivative as well as the function
tests :: [Cl3 -> IO()]
tests = [runTest "Log.Exp Identity" (log.exp) id (const False)
,runTest "Exp.Log Identity" (exp.log) id (\z -> lsv z < tol) -- singular inputs are out of bounds
,runTest "Abs*Signum Identity" (\x->abs x * signum x) id (const False)
,runTest "Reciprical Identity" (recip.recip) id (\z -> lsv z < tol) -- singular inputs are out of bounds
,runTest "sin.asin" (sin.asin) id (\z -> hasNilpotent z && poles [R 1, R (-1)] z)
,runTest "asin.sin" (asin.sin) (\z -> negate (I 1) * log (0.5 * (exp (I 1 * z) - exp (I (-1) * z)) +
sqrt (1+0.25*(exp (I (-1) * z) - exp (I 1 * z))^2))) (const False)
,runTest "cos.acos" (cos.acos) id (\z -> hasNilpotent z && poles [R 1, R (-1)] z)
,runTest "acos.cos" (acos.cos) (\z -> 0.5 * (pi - 2 * asin(cos z))) (\z -> hasNilpotent z && poles [R 0, pi, negate pi] z)
,runTest "sinh.asinh" (sinh.asinh) id (const False)
,runTest "asinh.sinh" (asinh.sinh) (\z -> log (0.5*(exp z - exp (negate z)) + sqrt (0.25 * (exp z - exp (negate z))^2 + 1))) (const False)
,runTest "cosh.acosh" (cosh.acosh) id (\z -> hasNilpotent z && poles [R 1, R (-1)] z)
,runTest "acosh.cosh" (acosh.cosh) (\z -> log (0.5*(exp z + exp (negate z)) +
sqrt (0.5*(exp z + exp (negate z)) - 1) * sqrt (0.5*(exp z + exp (negate z)) + 1))) (const False)
,runTest "Double Angle sin" (\z -> sin (2 * z)) (\z -> 2 * sin z * cos z) (const False)
,runTest "Double Angle cos" (\z -> cos (2 * z)) (\z -> cos z ^ 2 - sin z ^ 2) (const False)
,runTest "Double Angle tan" (\z -> tan (2 * z)) (\z -> (2 * tan z) / (1 - tan z ^ 2)) (const False)
,runTest "+Sin Shift" (\z -> sin (pi/2 + z)) cos (const False)
,runTest "-Sin Shift" (\z -> sin (pi/2 - z)) cos (const False)
,runTest "Double Angle sinh" (\z -> sinh (2 * z)) (\z -> 2 * sinh z * cosh z) (const False)
,runTest "Double Angle cosh" (\z -> cosh (2 * z)) (\z -> 2 * cosh z ^ 2 - 1) (const False)
,runTest "Double Angle tanh" (\z -> tanh (2 * z)) (\z -> (2 * tanh z) / (1 + tanh z ^ 2)) (const False)
,runTest "sin^2+cos^2" (\z -> sin z ^ 2 + cos z ^ 2) (const $ R 1) (const False)
,runTest "cosh^2-sinh^2" (\z -> cosh z ^ 2 - sinh z ^ 2) (const $ R 1) (const False)
,runTest "Symetry of cosh" (cosh.negate) cosh (const False)
,runTest "Symetry of sinh" (sinh.negate) (negate.sinh) (const False)
,runTest "sin.acos" (sin.acos) (\z -> sqrt (1 - z^2)) (\z -> hasNilpotent z && poles [R 1, R (-1)] z)
,runTest "sin.atan" (sin.atan) (\z -> z / sqrt (1 + z^2)) (poles [I 1, I (-1)])
,runTest "cos.atan" (cos.atan) (\z -> recip.sqrt $ 1 + z^2) (poles [I 1, I (-1)])
,runTest "cos.asin" (cos.asin) (\z -> sqrt (1 - z^2)) (\z -> hasNilpotent z && poles [R 1, R (-1)] z)
,runTest "tan.asin" (tan.asin) (\z -> z / sqrt (1 - z^2)) (poles [R 1, R (-1)])
,runTest "tan.acos" (tan.acos) (\z -> sqrt (1 - z^2) / z) (\z -> if hasNilpotent z then poles [R 1, R 0, R (-1)] z else poles [R 0] z)
]
-- | The Properties
prop_LogExp :: ArbCl3 -> Bool
prop_LogExp (Arb cliffor) = (abs cliffor > 10) || (
let cliffor' = unWrapIPartEigs cliffor -- imaginary part of log.exp repeats
-- round off errors get large for exp larger than 5 use spectproj (log.exp) for accuracy
in log (exp cliffor') ≈≈ cliffor')
-- log 0 is -Inf, Infinite vectors don't play nice
-- spectproj (exp.log) doesn't have this issue
prop_ExpLog :: ArbCl3 -> Bool
prop_ExpLog (Arb cliffor) = (lsv cliffor < tol) || (exp (log cliffor) ≈≈ cliffor)
prop_AbsSignum :: ArbCl3 -> Bool
prop_AbsSignum (Arb cliffor) = abs cliffor * signum cliffor ≈≈ cliffor
prop_RecipDef :: ArbCl3 -> Bool
prop_RecipDef (Arb cliffor) = (lsv cliffor < tol) || (recip cliffor * cliffor ≈≈ 1)
-- singular inputs don't recip also suffers from roundoff errors at large values
prop_RecipID :: ArbCl3 -> Bool
prop_RecipID (Arb cliffor) = (lsv cliffor < tol) || (recip (recip cliffor) ≈≈ cliffor)
prop_SinAsin :: ArbCl3 -> Bool
prop_SinAsin (Arb cliffor) = if hasNilpotent cliffor
then poles [R 1, R (-1)] cliffor || (sin (asin cliffor) ≈≈ cliffor)
else sin (asin cliffor) ≈≈ cliffor
prop_AsinSin :: ArbCl3 -> Bool
prop_AsinSin (Arb cliffor) = (abs cliffor > 10) || (asin (sin cliffor) ≈≈ (I (-1) * log (0.5 * (exp (I 1 * cliffor) - exp (I (-1) * cliffor)) +
sqrt (1+0.25*(exp (I (-1) * cliffor) - exp (I 1 * cliffor))^2))))
prop_CosAcos :: ArbCl3 -> Bool
prop_CosAcos (Arb cliffor) = if hasNilpotent cliffor
then poles [R 1, R (-1)] cliffor || (cos (acos cliffor) ≈≈ cliffor)
else cos (acos cliffor) ≈≈ cliffor
prop_AcosCos :: ArbCl3 -> Bool
prop_AcosCos (Arb cliffor) = (abs cliffor > 10) || (if hasNilpotent cliffor
then poles [R 0, pi, negate pi] cliffor || (acos (cos cliffor) ≈≈ 0.5 * (pi - 2 * asin(cos cliffor)))
else acos (cos cliffor) ≈≈ 0.5 * (pi - 2 * asin(cos cliffor)))
prop_SinhAsinh :: ArbCl3 -> Bool
prop_SinhAsinh (Arb cliffor) = sinh (asinh cliffor) ≈≈ cliffor
prop_AsinhSinh :: ArbCl3 -> Bool
prop_AsinhSinh (Arb cliffor) = (abs cliffor > 10) || (asinh (sinh cliffor) ≈≈ log (0.5*(exp cliffor - exp (negate cliffor)) +
sqrt (0.25 * (exp cliffor - exp (negate cliffor))^2 + 1)))
prop_CoshAcosh :: ArbCl3 -> Bool
prop_CoshAcosh (Arb cliffor) = if hasNilpotent cliffor
then poles [R 1, R (-1)] cliffor || (cosh (acosh cliffor) ≈≈ cliffor)
else cosh (acosh cliffor) ≈≈ cliffor
prop_AcoshCosh :: ArbCl3 -> Bool
prop_AcoshCosh (Arb cliffor) = acosh (cosh cliffor) ≈≈ log (0.5*(exp cliffor + exp (negate cliffor)) +
sqrt (0.5*(exp cliffor + exp (negate cliffor)) - 1) *
sqrt (0.5*(exp cliffor + exp (negate cliffor)) + 1))
prop_AcoshCosh2 :: ArbCl3 -> Bool
prop_AcoshCosh2 (Arb cliffor) = acosh (cosh cliffor) ≈≈ log (cosh cliffor + sqrt (cosh cliffor - 1) * sqrt (cosh cliffor + 1))
prop_DubSin :: ArbCl3 -> Bool
prop_DubSin (Arb cliffor) = sin (2 * cliffor) ≈≈ 2 * sin cliffor * cos cliffor
prop_DubCos :: ArbCl3 -> Bool
prop_DubCos (Arb cliffor) = cos (2 * cliffor) ≈≈ cos cliffor ^ 2 - sin cliffor ^ 2
prop_DubTan :: ArbCl3 -> Bool
prop_DubTan (Arb cliffor) = tan (2 * cliffor) ≈≈ (2 * tan cliffor) / (1 - tan cliffor ^ 2)
prop_DubSinh :: ArbCl3 -> Bool
prop_DubSinh (Arb cliffor) = sinh (2 * cliffor) ≈≈ 2 * sinh cliffor * cosh cliffor
prop_DubCosh :: ArbCl3 -> Bool
prop_DubCosh (Arb cliffor) = cosh (2 * cliffor) ≈≈ 2 * cosh cliffor ^ 2 - 1
prop_DubTanh :: ArbCl3 -> Bool
prop_DubTanh (Arb cliffor) = tanh (2 * cliffor) ≈≈ (2 * tanh cliffor) / (1 + tanh cliffor ^ 2)
prop_PosSinShift :: ArbCl3 -> Bool
prop_PosSinShift (Arb cliffor) = sin (pi/2 + cliffor) ≈≈ cos cliffor
prop_NegSinShift :: ArbCl3 -> Bool
prop_NegSinShift (Arb cliffor) = sin (pi/2 - cliffor) ≈≈ cos cliffor
prop_SinSqCosSq :: ArbCl3 -> Bool
prop_SinSqCosSq (Arb cliffor) = (abs cliffor > 10) || (sin cliffor ^ 2 + cos cliffor ^ 2 ≈≈ 1)
prop_CoshSqmSinhSq :: ArbCl3 -> Bool
prop_CoshSqmSinhSq (Arb cliffor) = (abs cliffor > 10) || (cosh cliffor ^ 2 - sinh cliffor ^ 2 ≈≈ 1)
prop_SymCosh :: ArbCl3 -> Bool
prop_SymCosh (Arb cliffor) = cosh (negate cliffor) ≈≈ cosh cliffor
prop_SymSinh :: ArbCl3 -> Bool
prop_SymSinh (Arb cliffor) = sinh (negate cliffor) ≈≈ negate (sinh cliffor)
prop_DoubleISin :: ArbCl3 -> Bool
prop_DoubleISin (Arb cliffor) = 2 * I 1 * sin cliffor ≈≈ exp(I 1 * cliffor) - exp (I (-1) * cliffor)
-- | Composition Sub-Algebras have a distributive norm over multiplication, like this:
--
-- > norm $ clif * clif' = norm clif * norm clif'
--
-- Strangly the constructor combinations with the "= True" don't play nice with 'abs'
-- they are the constructors with non-zero zero-divisors.
prop_CompAlg :: (ArbCl3, ArbCl3) -> Bool
prop_CompAlg (Arb PV{}, Arb PV{}) = True
prop_CompAlg (Arb PV{}, Arb BPV{}) = True
prop_CompAlg (Arb PV{}, Arb TPV{}) = True
prop_CompAlg (Arb PV{}, Arb APS{}) = True
prop_CompAlg (Arb BPV{}, Arb PV{}) = True
prop_CompAlg (Arb TPV{}, Arb PV{}) = True
prop_CompAlg (Arb APS{}, Arb PV{}) = True
prop_CompAlg (Arb BPV{}, Arb BPV{}) = True
prop_CompAlg (Arb BPV{}, Arb TPV{}) = True
prop_CompAlg (Arb BPV{}, Arb APS{}) = True
prop_CompAlg (Arb TPV{}, Arb BPV{}) = True
prop_CompAlg (Arb APS{}, Arb BPV{}) = True
prop_CompAlg (Arb TPV{}, Arb TPV{}) = True
prop_CompAlg (Arb TPV{}, Arb APS{}) = True
prop_CompAlg (Arb APS{}, Arb TPV{}) = True
prop_CompAlg (Arb APS{}, Arb APS{}) = True
prop_CompAlg (Arb cliffor, Arb cliffor') = abs ( cliffor * cliffor') ≈≈ abs cliffor * abs cliffor'
-- Run the test
-- compare the function under test (fUT) to a golden test funcion (gTF)
-- if the input is within bounds
runTest :: String -> (Cl3 -> Cl3) -> (Cl3 -> Cl3) -> (Cl3 -> Bool) -> Cl3 -> IO()
runTest testName fUT gTF outOB iVal =
let f = fUT iVal
g = gTF iVal
in if outOB iVal
then putStr (unlines [testName ++ ": Input Out of Bounds"])
else if f ≈≈ g
then putStr (unlines [testName ++ ": Passed"])
else putStr (unlines [testName ++ ": Failed"
,"Expected: " ++ show g
," got: " ++ show f
,"on input: " ++ show iVal
])
-- | '≈≈' aproximately equal
(≈≈) :: Cl3 -> Cl3 -> Bool
(reduce -> clifforA) ≈≈ (reduce -> clifforB) =
let ave = (abs clifforA + abs clifforB) / 2
in abs (clifforA - clifforB) <= 1e-5*ave + tol
infix 4 ≈≈
-- | 'poles' a function that tests if a cliffor is one of the defined poles
poles :: [Cl3] -> Cl3 -> Bool
poles [] _ = False
poles [p] cliffor = eig1 ≈≈ p || eig2 ≈≈ p
where (eig1,eig2) = eigvals cliffor
poles (p:ps) cliffor = (eig1 ≈≈ p || eig2 ≈≈ p) || poles ps cliffor
where (eig1,eig2) = eigvals cliffor
unWrapIPartEigs :: Cl3 -> Cl3
unWrapIPartEigs cliffor = reduce $ spectraldcmp unWrapI id cliffor
where unWrapI (R a0) = R a0
unWrapI (I a123) | a123 > pi = unWrapI $ I (a123 - 2*pi)
| a123 < (-pi) = unWrapI $ I (a123 + 2*pi)
| otherwise = I a123
unWrapI (C a0 a123) | a123 > pi = unWrapI $ C a0 (a123 - 2*pi)
| a123 < (-pi) = unWrapI $ C a0 (a123 + 2*pi)
| otherwise = C a0 a123
unWrapI _ = error "unWrapI should only be unWrapping R I and C"
----------------------------------------------------------
-------------------------------------------------------------------
--
-- Arbitrary Instance of Cl3 types, typically for use with the
-- "Test.QuickCheck" library.
--
-------------------------------------------------------------------
-- | 'ArbCl3' to provide a newtype wrapper to avoid the orphan instance
newtype ArbCl3 = Arb Cl3 deriving (Show)
-- | 'Arbitrary' instance that has its largest singular value less than or equal to 15
instance Arbitrary ArbCl3 where
arbitrary =
oneof [(Arb.)R <$> arbitrary,
((Arb.).).V3 <$> arbitrary <*> arbitrary <*> arbitrary,
((Arb.).).BV <$> arbitrary <*> arbitrary <*> arbitrary,
(Arb.)I <$> arbitrary,
(((Arb.).).).PV <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary,
(((Arb.).).).H <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary,
(Arb.).C <$> arbitrary <*> arbitrary,
(((((Arb.).).).).).BPV <$> arbitrary <*> arbitrary <*> arbitrary
<*> arbitrary <*> arbitrary <*> arbitrary,
(((Arb.).).).ODD <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary,
(((Arb.).).).TPV <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary,
(((((((Arb.).).).).).).).APS <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary
<*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary
] `suchThat` lessThan15
where
lessThan15 (Arb cliffor) = abs cliffor <= 15