easytensor-2.1.0.0: test/Numeric/MatrixDoubleTest.hs
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeOperators #-}
module Numeric.MatrixDoubleTest (runTests) where
import Data.Fixed
import Numeric.Arbitraries
import Numeric.DataFrame
import Numeric.Dimensions
import Test.QuickCheck
type TestElem = Double
dropW :: (SubSpace t '[3] '[] '[3], SubSpace t '[4] '[] '[4])
=> Vector t 4 -> Vector t 3
dropW (Vec4 x y z _) = Vec3 x y z
prop_detTranspose :: SomeSquareMatrix AnyMatrix TestElem -> Property
prop_detTranspose (SSM m') = counterexample
(unlines
[ "Failed det/transpose:"
, "m: " ++ show m
, "mT: " ++ show (transpose m)
]
) $ approxEq mag a b
where
n2 = totalDim (dims `inSpaceOf` m)
mag = fromIntegral (n2*n2)
mm = recip $ maxElem m' `max` 1
m = ewmap (*mm) m'
a = det m
b = det $ transpose m
prop_inverse :: SomeSquareMatrix NonSingular TestElem -> Property
prop_inverse (SSM m) = check (m %* mi) .&&. check (mi %* m)
where
n2 = totalDim (dims `inSpaceOf` m)
mag = maxElem m `max` maxElem mi `max` 1
* fromIntegral (n2*n2)
mi = inverse m
check a = counterexample
( unlines
[ "Failed inverse:"
, " m: " ++ show m
, " mi: " ++ show mi
, " eye': " ++ show a
]
) $ approxEq mag a eye
prop_translate3vs4 :: Vector TestElem 4 -> Bool
prop_translate3vs4 v = translate4 v == translate3 (dropW v)
prop_translate4 :: Vector TestElem 4 -> Vector TestElem 3 -> Bool
prop_translate4 a b = toHomPoint b %* translate4 a == toHomPoint (dropW a + b)
prop_translate3 :: Vector TestElem 3 -> Vector TestElem 3 -> Bool
prop_translate3 a b = toHomPoint b %* translate3 a == toHomPoint (a + b)
prop_rotateX :: Vector TestElem 4 -> Property
prop_rotateX v@(Vec4 x y z w) =
conjoin [
v %* rotateX (-2 * pi) =~= v,
v %* rotateX (-1.5 * pi) =~= vec4 x (-z) y w,
v %* rotateX (-pi) =~= vec4 x (-y) (-z) w,
v %* rotateX (-0.5 * pi) =~= vec4 x z (-y) w,
v %* rotateX 0 =~= v,
v %* rotateX (0.5 * pi) =~= vec4 x (-z) y w,
v %* rotateX pi =~= vec4 x (-y) (-z) w,
v %* rotateX (1.5 * pi) =~= vec4 x z (-y) w,
v %* rotateX (2 * pi) =~= v
]
prop_rotateY :: Vector TestElem 4 -> Property
prop_rotateY v@(Vec4 x y z w) =
conjoin [
v %* rotateY (-2 * pi) =~= v,
v %* rotateY (-1.5 * pi) =~= vec4 z y (-x) w,
v %* rotateY (-pi) =~= vec4 (-x) y (-z) w,
v %* rotateY (-0.5 * pi) =~= vec4 (-z) y x w,
v %* rotateY 0 =~= v,
v %* rotateY (0.5 * pi) =~= vec4 z y (-x) w,
v %* rotateY pi =~= vec4 (-x) y (-z) w,
v %* rotateY (1.5 * pi) =~= vec4 (-z) y x w,
v %* rotateY (2 * pi) =~= v
]
prop_rotateZ :: Vector TestElem 4 -> Property
prop_rotateZ v@(Vec4 x y z w) =
conjoin [
v %* rotateZ (-2 * pi) =~= v,
v %* rotateZ (-1.5 * pi) =~= vec4 (-y) x z w,
v %* rotateZ (-pi) =~= vec4 (-x) (-y) z w,
v %* rotateZ (-0.5 * pi) =~= vec4 y (-x) z w,
v %* rotateZ 0 =~= v,
v %* rotateZ (0.5 * pi) =~= vec4 (-y) x z w,
v %* rotateZ pi =~= vec4 (-x) (-y) z w,
v %* rotateZ (1.5 * pi) =~= vec4 y (-x) z w,
v %* rotateZ (2 * pi) =~= v
]
prop_rotate :: TestElem -> Property
prop_rotate a =
conjoin [
rotate (vec3 1 0 0) a =~= rotateX a,
rotate (vec3 0 1 0) a =~= rotateY a,
rotate (vec3 0 0 1) a =~= rotateZ a
]
prop_rotateEuler :: TestElem -> TestElem -> TestElem -> Property
prop_rotateEuler pitch yaw roll =
rotateEuler pitch yaw roll =~= rotateZ roll %* rotateY yaw %* rotateX pitch
prop_lookAt :: Vector TestElem 3 -> Vector TestElem 3 -> Vector TestElem 3 -> Property
prop_lookAt up cam foc =
(apart cam foc && nonColinear up zb && nonColinear zb xb) ==>
conjoin [
(normalized . fromHom $ toHomPoint foc %* m) `aeq` vec3 0 0 (-1),
fromHom (toHomPoint cam %* m) `aeq` 0,
fromHom (toHomVector xb %* m) `aeq` vec3 1 0 0,
fromHom (toHomVector yb %* m) `aeq` vec3 0 1 0,
fromHom (toHomVector zb %* m) `aeq` vec3 0 0 1
]
where
aeq :: Vector TestElem 3 -> Vector TestElem 3 -> Property
aeq = approxEq 10
apart :: Vector TestElem 3 -> Vector TestElem 3 -> Bool
apart a b = normL2 (a - b) > max (2 * sqrt M_EPS) (0.1 * (normL2 a + normL2 b))
nonColinear :: Vector TestElem 3 -> Vector TestElem 3 -> Bool
nonColinear a b = normL2 (cross a b) > 0.1 * abs (dot a b)
m = lookAt up cam foc
zb = normalized $ cam - foc
xb = normalized $ up `cross` zb
yb = zb `cross` xb
prop_perspective :: TestElem -> TestElem -> TestElem -> TestElem -> Property
prop_perspective a b c d =
conjoin [
projectTo 0 0 n =~= vec3 0 0 (-1),
projectTo 0 0 f =~= vec3 0 0 1,
projectTo 1 1 n =~= vec3 1 1 (-1),
projectTo 1 (-1) n =~= vec3 1 (-1) (-1),
projectTo (-1) 1 n =~= vec3 (-1) 1 (-1),
projectTo (-1) (-1) n =~= vec3 (-1) (-1) (-1),
projectTo 1 1 f =~= vec3 1 1 1,
projectTo 1 (-1) f =~= vec3 1 (-1) 1,
projectTo (-1) 1 f =~= vec3 (-1) 1 1,
projectTo (-1) (-1) f =~= vec3 (-1) (-1) 1
]
where
n = 1.0 + mod' a 9.0 -- Near plane in range [1, 10)
f = 2*n + mod' b 99.0 -- Far plane in range [2*n, 2*n + 100)
fovy = (0.1 * pi) + mod' c (0.8 * pi) -- Y-axis field of view in range [0.1*pi, 0.9*pi)
aspect = 0.25 + mod' d 4.0 -- Aspect ration in range [1/4, 4/1]
hpd = tan (fovy * 0.5) -- height/distance
wpd = aspect * hpd -- width/distance
m = perspective n f fovy aspect
projectTo x' y' z = fromHom $ vec4 (x' * wpd * z) (y' * hpd * z) (-z) 1 %* m
prop_orthogonal :: TestElem -> TestElem -> TestElem -> TestElem -> Property
prop_orthogonal a b c d =
conjoin [
projectTo 0 0 n =~= vec3 0 0 (-1),
projectTo 0 0 f =~= vec3 0 0 1,
projectTo 1 1 n =~= vec3 1 1 (-1),
projectTo 1 (-1) n =~= vec3 1 (-1) (-1),
projectTo (-1) 1 n =~= vec3 (-1) 1 (-1),
projectTo (-1) (-1) n =~= vec3 (-1) (-1) (-1),
projectTo 1 1 f =~= vec3 1 1 1,
projectTo 1 (-1) f =~= vec3 1 (-1) 1,
projectTo (-1) 1 f =~= vec3 (-1) 1 1,
projectTo (-1) (-1) f =~= vec3 (-1) (-1) 1
]
where
n = 1.0 + mod' a 9.0 -- Near plane in range [1, 10)
f = 2*n + mod' b 99.0 -- Far plane in range [2*n, 2*n + 100)
w = 1.0 + mod' c 9999.0 -- Width in range [1, 10000)
h = 1.0 + mod' d 9999.0 -- Height in range [1, 10000)
m = orthogonal n f w h
projectTo x' y' z = fromHom $ vec4 (x' * w * 0.5) (y' * h * 0.5) (-z) 1 %* m
prop_toHomPoint :: Vector TestElem 3 -> Bool
prop_toHomPoint v@(Vec3 x y z) = toHomPoint v == vec4 x y z 1
prop_toHomVector :: Vector TestElem 3 -> Bool
prop_toHomVector v@(Vec3 x y z) = toHomVector v == vec4 x y z 0
prop_fromHom :: Vector TestElem 4 -> Property
prop_fromHom v@(Vec4 x y z w) =
case w of
0 -> fromHom v === vec3 x y z
_ -> fromHom v =~= vec3 (x/w) (y/w) (z/w)
return []
runTests :: Int -> IO Bool
runTests n = $forAllProperties
$ quickCheckWithResult stdArgs { maxSuccess = n }