diff --git a/easytensor.cabal b/easytensor.cabal --- a/easytensor.cabal +++ b/easytensor.cabal @@ -1,5 +1,5 @@ name: easytensor -version: 1.0.0.1 +version: 1.0.1.0 cabal-version: >=1.22 build-type: Simple license: BSD3 @@ -65,12 +65,20 @@ base >=4.10 && <5, dimensions >=1.0.0.0 default-language: Haskell2010 + if impl(ghc >= 8.6) + default-extensions: NoStarIsType hs-source-dirs: src other-modules: + Numeric.Matrix.Mat44d + Numeric.Matrix.Mat44f Numeric.Quaternion.Class Numeric.Quaternion.QDouble Numeric.Quaternion.QFloat - ghc-options: -Wall -fwarn-tabs -fwarn-unused-do-bind -fwarn-monomorphism-restriction + ghc-options: + -Wall + -fwarn-tabs + -fwarn-unused-do-bind + -fwarn-monomorphism-restriction test-suite et-test @@ -81,7 +89,8 @@ Numeric.DataFrame.Arbitraries Numeric.DataFrame.SubSpaceTest Numeric.DataFrame.BasicTest - Numeric.MatrixTest + Numeric.MatrixDoubleTest + Numeric.MatrixFloatTest Numeric.QuaternionTest build-depends: base -any, @@ -90,6 +99,8 @@ easytensor -any, dimensions -any default-language: Haskell2010 + if impl(ghc >= 8.6) + default-extensions: NoStarIsType, NoMonadFailDesugaring hs-source-dirs: test ghc-options: -Wall -fwarn-tabs -O2 @@ -103,6 +114,8 @@ easytensor -any, dimensions -any default-language: Haskell2010 + if impl(ghc >= 8.6) + default-extensions: NoStarIsType, NoMonadFailDesugaring hs-source-dirs: bench ghc-options: -Wall -fwarn-tabs -O2 @@ -117,5 +130,7 @@ dimensions -any, time -any default-language: Haskell2010 + if impl(ghc >= 8.6) + default-extensions: NoStarIsType, NoMonadFailDesugaring hs-source-dirs: bench ghc-options: -Wall -fwarn-tabs -O2 diff --git a/src/Numeric/DataFrame/Internal/Array/Family/ArrayBase.hs b/src/Numeric/DataFrame/Internal/Array/Family/ArrayBase.hs --- a/src/Numeric/DataFrame/Internal/Array/Family/ArrayBase.hs +++ b/src/Numeric/DataFrame/Internal/Array/Family/ArrayBase.hs @@ -1,23 +1,21 @@ -{-# LANGUAGE ConstraintKinds #-} -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE ExistentialQuantification #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE MagicHash #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE Rank2Types #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE TypeApplications #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE TypeFamilyDependencies #-} -{-# LANGUAGE TypeInType #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE UnboxedSums #-} -{-# LANGUAGE UnboxedTuples #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeFamilyDependencies #-} +{-# LANGUAGE TypeInType #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UnboxedSums #-} +{-# LANGUAGE UnboxedTuples #-} module Numeric.DataFrame.Internal.Array.Family.ArrayBase ( ArrayBase (..) @@ -25,7 +23,7 @@ import Data.Int import Data.Word -import GHC.Base hiding (foldr) +import GHC.Base hiding (foldr) import Numeric.DataFrame.Internal.Array.Class import Numeric.DataFrame.Internal.Array.PrimOps import Numeric.Dimensions @@ -44,6 +42,7 @@ | (# Int# -- Offset measured in elements. , Int# -- Number of elements. , ByteArray# -- Content. + , Evidence (PrimBytes t) #) #) @@ -67,7 +66,7 @@ | W# nw <- totalDim' @ds , n <- word2Int# nw , tbs <- byteSize t -> go tbs (tbs *# n) t - (# | (# _, _, arr #) #) -> + (# | (# _, _, arr, _ #) #) -> -- very weird trick with touch# allows to workaround GHC bug -- "internal error: ARR_WORDS object entered!" -- TODO: report this @@ -87,7 +86,7 @@ , tbs <- byteSize (undefined :: t) , (# offN, offRem #) <- quotRemInt# bOff tbs = case offRem of - 0# -> ArrayBase (# | (# offN, n , ba #) #) + 0# -> ArrayBase (# | (# offN, n , ba, E #) #) _ -> go n (tbs *# n) where go n bsize = case runRW# @@ -99,7 +98,7 @@ of (# s1, mba #) -> unsafeFreezeByteArray# mba (copyByteArray# ba bOff mba 0# bsize s1) - ) of (# _, r #) -> ArrayBase (# | (# 0# , n , r #) #) + ) of (# _, r #) -> ArrayBase (# | (# 0# , n , r, E #) #) {-# NOINLINE go #-} {-# INLINE fromBytes #-} @@ -111,14 +110,14 @@ = case newByteArray# bsize s0 of (# s1, mba1 #) -> case unsafeFreezeByteArray# mba1 (copyMutableByteArray# mba bOff mba1 0# bsize s1) of - (# s2, ba #) -> (# s2, ArrayBase (# | (# 0# , n , ba #) #) #) + (# s2, ba #) -> (# s2, ArrayBase (# | (# 0# , n , ba, E #) #) #) {-# INLINE readBytes #-} writeBytes mba bOff (ArrayBase c) | tbs <- byteSize (undefined :: t) = case c of (# t | #) | W# n <- totalDim' @ds -> loop# bOff tbs (bOff +# word2Int# n *# tbs) (\i -> writeBytes mba i t) - (# | (# offN, n, arr #) #) -> + (# | (# offN, n, arr, _ #) #) -> copyByteArray# arr (offN *# tbs) mba bOff (n *# tbs) {-# INLINE writeBytes #-} @@ -130,14 +129,14 @@ = case newByteArray# bsize s0 of (# s1, mba1 #) -> case unsafeFreezeByteArray# mba1 (copyAddrToByteArray# addr mba1 0# bsize s1) of - (# s2, ba #) -> (# s2, ArrayBase (# | (# 0# , n , ba #) #) #) + (# s2, ba #) -> (# s2, ArrayBase (# | (# 0# , n , ba, E #) #) #) {-# INLINE readAddr #-} writeAddr (ArrayBase c) addr | tbs <- byteSize (undefined :: t) = case c of (# t | #) | W# n <- totalDim' @ds -> loop# 0# tbs (word2Int# n *# tbs) (\i -> writeAddr t (plusAddr# addr i)) - (# | (# offN, n, arr #) #) -> + (# | (# offN, n, arr, _ #) #) -> copyByteArrayToAddr# arr (offN *# tbs) addr (n *# tbs) {-# INLINE writeAddr #-} @@ -150,25 +149,24 @@ {-# INLINE byteAlign #-} byteOffset (ArrayBase a) = case a of - (# _ | #) -> 0# - (# | (# off, _, _ #) #) -> off *# byteSize (undefined :: t) + (# _ | #) -> 0# + (# | (# off, _, _, _ #) #) -> off *# byteSize (undefined :: t) {-# INLINE byteOffset #-} indexArray ba off | W# nw <- totalDim' @ds , n <- word2Int# nw - = ArrayBase (# | (# off *# n, n, ba #) #) + = ArrayBase (# | (# off *# n, n, ba, E #) #) {-# INLINE indexArray #-} -- | Accumulates only idempotent operations! -- Being applied to FromScalars, executes only once! --- Here, idempotance means: assuming @f a b = g @, @g (g x) = g x@ +-- Here, idempotance means: assuming @f a b = g x@, @g (g x) = g x@ -- --- Also, I assume the size of arrays is the same -accumV2Idempotent :: PrimBytes t - => a +-- Also, I assume the sizes of arrays are the same +accumV2Idempotent :: a -> (t -> t -> a -> a) -> ArrayBase t ds -> ArrayBase t ds -> a accumV2Idempotent x f @@ -176,23 +174,23 @@ (ArrayBase (# b | #)) = f a b x accumV2Idempotent x f - a@(ArrayBase (# | (# _, nA, _ #) #)) - b@(ArrayBase (# | (# _, nB, _ #) #)) + a@(ArrayBase (# | (# _, nA, _, E #) #)) + b@(ArrayBase (# | (# _, nB, _, _ #) #)) = loop1a# (minInt# nA nB) (\i -> f (ix# i a) (ix# i b)) x accumV2Idempotent x f (ArrayBase (# a | #)) - b@(ArrayBase (# | (# _, n, _ #) #)) + b@(ArrayBase (# | (# _, n, _, E #) #)) = loop1a# n (\i -> f a (ix# i b)) x accumV2Idempotent x f - a@(ArrayBase (# | (# _, n, _ #) #)) + a@(ArrayBase (# | (# _, n, _, E #) #)) (ArrayBase (# b | #)) = loop1a# n (\i -> f (ix# i a) b) x {-# INLINE accumV2Idempotent #-} -mapV :: PrimBytes t => (t -> t) -> ArrayBase t ds -> ArrayBase t ds +mapV :: (t -> t) -> ArrayBase t ds -> ArrayBase t ds mapV f (ArrayBase (# t | #)) = ArrayBase (# f t | #) -mapV f x@(ArrayBase (# | (# offN, n, ba #) #)) +mapV f x@(ArrayBase (# | (# offN, n, ba, E #) #)) | tbs <- byteSize (undefEl x) = go (tbs *# n) where @@ -202,18 +200,18 @@ ( loop1# n (\i -> writeArray mba i (f (indexArray ba (offN +# i)))) s1 ) - ) of (# _, r #) -> ArrayBase (# | (# 0#, n, r #) #) + ) of (# _, r #) -> ArrayBase (# | (# 0#, n, r, E #) #) {-# NOINLINE go #-} {-# INLINE mapV #-} -zipV :: PrimBytes t => (t -> t -> t) +zipV :: (t -> t -> t) -> ArrayBase t ds -> ArrayBase t ds -> ArrayBase t ds zipV f (ArrayBase (# x | #)) b = mapV (f x) b -zipV f a (ArrayBase (# y | #)) = mapV (flip f y) a -zipV f a@(ArrayBase (# | (# oa, na, ba #) #)) - (ArrayBase (# | (# ob, nb, bb #) #)) - | n <- (minInt# na nb) +zipV f a (ArrayBase (# y | #)) = mapV (`f` y) a +zipV f a@(ArrayBase (# | (# oa, na, ba, E #) #)) + (ArrayBase (# | (# ob, nb, bb, _ #) #)) + | n <- minInt# na nb = go n (byteSize (undefEl a) *# n) where go n bsize = case runRW# @@ -226,7 +224,7 @@ ) ) s1 ) - ) of (# _, r #) -> ArrayBase (# | (# 0#, n, r #) #) + ) of (# _, r #) -> ArrayBase (# | (# 0#, n, r, E #) #) {-# NOINLINE go #-} {-# INLINE zipV #-} @@ -234,7 +232,7 @@ -- TODO: to improve performance, I can either compare bytearrays using memcmp -- or implement early termination if the first elements do not match. -- On the other hand, hopefully @(&&)@ and @(||)@ ops take care of that. -instance (Eq t, PrimBytes t) => Eq (ArrayBase t ds) where +instance Eq t => Eq (ArrayBase t ds) where {-# SPECIALIZE instance Eq (ArrayBase Float ds) #-} {-# SPECIALIZE instance Eq (ArrayBase Double ds) #-} {-# SPECIALIZE instance Eq (ArrayBase Int ds) #-} @@ -252,7 +250,7 @@ -- | Implement partial ordering for `>`, `<`, `>=`, `<=` -- and lexicographical ordering for `compare` -instance (Ord t, PrimBytes t) => Ord (ArrayBase t ds) where +instance Ord t => Ord (ArrayBase t ds) where {-# SPECIALIZE instance Ord (ArrayBase Float ds) #-} {-# SPECIALIZE instance Ord (ArrayBase Double ds) #-} {-# SPECIALIZE instance Ord (ArrayBase Int ds) #-} @@ -287,10 +285,10 @@ max = zipV max {-# INLINE max #-} -instance (Dimensions ds, PrimBytes t, Show t) +instance (Dimensions ds, Show t) => Show (ArrayBase t ds) where show x = case dims @_ @ds of - U -> "{ " ++ show (ix# 0# x) ++ " }" + U -> "{ " ++ show (ix U x) ++ " }" Dim :* U -> ('{' :) . drop 1 $ foldr (\i s -> ", " ++ show (ix i x) ++ s) " }" [minBound .. maxBound] @@ -330,7 +328,7 @@ maxBound = ArrayBase (# maxBound | #) minBound = ArrayBase (# minBound | #) -instance (Num t, PrimBytes t) => Num (ArrayBase t ds) where +instance Num t => Num (ArrayBase t ds) where {-# SPECIALIZE instance Num (ArrayBase Float ds) #-} {-# SPECIALIZE instance Num (ArrayBase Double ds) #-} {-# SPECIALIZE instance Num (ArrayBase Int ds) #-} @@ -358,7 +356,7 @@ fromInteger i = ArrayBase (# fromInteger i | #) {-# INLINE fromInteger #-} -instance (Fractional t, PrimBytes t) => Fractional (ArrayBase t ds) where +instance Fractional t => Fractional (ArrayBase t ds) where {-# SPECIALIZE instance Fractional (ArrayBase Float ds) #-} {-# SPECIALIZE instance Fractional (ArrayBase Double ds) #-} (/) = zipV (/) @@ -369,7 +367,7 @@ {-# INLINE fromRational #-} -instance (Floating t, PrimBytes t) => Floating (ArrayBase t ds) where +instance Floating t => Floating (ArrayBase t ds) where {-# SPECIALIZE instance Floating (ArrayBase Float ds) #-} {-# SPECIALIZE instance Floating (ArrayBase Double ds) #-} pi = ArrayBase (# pi | #) @@ -427,8 +425,8 @@ {-# INLINE broadcast #-} ix# i (ArrayBase a) = case a of - (# t | #) -> t - (# | (# off, _, arr #) #) -> indexArray arr (off +# i) + (# t | #) -> t + (# | (# off, _, arr, _ #) #) -> indexArray arr (off +# i) {-# INLINE ix# #-} gen# n f z0 = go (byteSize @t undefined *# n) @@ -438,7 +436,7 @@ (# s1, mba #) -> case loop0 mba 0# z0 s1 of (# s2, z1 #) -> case unsafeFreezeByteArray# mba s2 of (# s3, ba #) -> (# s3, (# z1, ba #) #) - ) of (# _, (# z1, ba #) #) -> (# z1, ArrayBase (# | (# 0# , n , ba #) #) #) + ) of (# _, (# z1, ba #) #) -> (# z1, ArrayBase (# | (# 0# , n , ba, E #) #) #) {-# NOINLINE go #-} loop0 mba i z s | isTrue# (i ==# n) = (# s, z #) @@ -454,9 +452,9 @@ (writeArray mba i x (loop1# n (\j -> writeArray mba j a) s1) ) - ) of (# _, r #) -> ArrayBase (# | (# 0# , n , r #) #) + ) of (# _, r #) -> ArrayBase (# | (# 0# , n , r, E #) #) {-# NOINLINE go #-} - upd# _ i x (ArrayBase (# | (# offN , n , ba #) #)) = go (byteSize x) + upd# _ i x (ArrayBase (# | (# offN , n , ba, E #) #)) = go (byteSize x) where go tbs = case runRW# ( \s0 -> case newByteArray# (tbs *# n) s0 of @@ -464,21 +462,21 @@ (writeArray mba i x (copyByteArray# ba (offN *# tbs) mba 0# (tbs *# n) s1) ) - ) of (# _, r #) -> ArrayBase (# | (# 0# , n , r #) #) + ) of (# _, r #) -> ArrayBase (# | (# 0# , n , r, E #) #) {-# NOINLINE go #-} {-# INLINE upd# #-} elemOffset (ArrayBase a) = case a of - (# _ | #) -> 0# - (# | (# off, _, _ #) #) -> off + (# _ | #) -> 0# + (# | (# off, _, _, _ #) #) -> off {-# INLINE elemOffset #-} elemSize0 (ArrayBase a) = case a of - (# _ | #) -> 0# - (# | (# _, n, _ #) #) -> n + (# _ | #) -> 0# + (# | (# _, n, _, _ #) #) -> n {-# INLINE elemSize0 #-} - fromElems off n ba = ArrayBase (# | (# off , n , ba #) #) + fromElems off n ba = ArrayBase (# | (# off , n , ba, E #) #) {-# INLINE fromElems #-} @@ -488,10 +486,10 @@ -------------------------------------------------------------------------------- -ix :: (PrimBytes t, Dimensions ds) => Idxs ds -> ArrayBase t ds -> t +ix :: Dimensions ds => Idxs ds -> ArrayBase t ds -> t ix i (ArrayBase a) = case a of (# t | #) -> t - (# | (# off, _, arr #) #) -> case fromEnum i of + (# | (# off, _, arr, E #) #) -> case fromEnum i of I# i# -> indexArray arr (off +# i#) {-# INLINE ix #-} diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs --- a/src/Numeric/Matrix.hs +++ b/src/Numeric/Matrix.hs @@ -34,6 +34,7 @@ , MatrixInverse (..) , MatrixLU (..), LUFact (..) , Matrix + , HomTransform4 (..) , Mat22f, Mat23f, Mat24f , Mat32f, Mat33f, Mat34f , Mat42f, Mat43f, Mat44f @@ -58,6 +59,8 @@ import Numeric.DataFrame.Type import Numeric.Dimensions import Numeric.Matrix.Class +import Numeric.Matrix.Mat44d () +import Numeric.Matrix.Mat44f () import Numeric.PrimBytes import Numeric.Scalar import Numeric.Vector diff --git a/src/Numeric/Matrix/Class.hs b/src/Numeric/Matrix/Class.hs --- a/src/Numeric/Matrix/Class.hs +++ b/src/Numeric/Matrix/Class.hs @@ -78,9 +78,9 @@ -- | Operations on 4x4 transformation matrices and vectors in homogeneous coordinates. -- All angles are specified in radians. class HomTransform4 t where - -- | Create a translation matrix from a vector + -- | Create a translation matrix from a vector. The 4th coordinate is ignored. translate4 :: Vector t 4 -> Matrix t 4 4 - -- | Create a translation matrix from a vector + -- | Create a translation matrix from a vector. translate3 :: Vector t 3 -> Matrix t 4 4 -- | Rotation matrix for a rotation around the X axis, angle is given in radians. rotateX :: t -> Matrix t 4 4 @@ -90,8 +90,16 @@ rotateZ :: t -> Matrix t 4 4 -- | Rotation matrix for a rotation around an arbitrary normalized vector rotate :: Vector t 3 -> t -> Matrix t 4 4 - -- | Rotation matrix from the Euler angles yaw pitch and roll - rotateEuler :: t -> t -> t -> Matrix t 4 4 + -- | Rotation matrix from the Euler angles roll (axis @Z@), yaw (axis @Y'@), and pitch (axis @X''@). + -- This order is known as Tait-Bryan angles (@Z-Y'-X''@ intrinsic rotations), or nautical angles, or Cardan angles. + -- + -- > rotateEuler pitch yaw roll == rotateX pitch %* rotateY yaw %* rotateZ roll + -- + -- https://en.wikipedia.org/wiki/Euler_angles#Conventions_2 + rotateEuler :: t -- ^ pitch (axis @X''@) + -> t -- ^ yaw (axis @Y'@) + -> t -- ^ roll (axis @Z@) + -> Matrix t 4 4 -- | Create a transform matrix using up direction, camera position and a point to look at. -- Just the same as GluLookAt. lookAt :: Vector t 3 -- ^ The up direction, not necessary unit length or perpendicular to the view vector diff --git a/src/Numeric/Matrix/Mat44d.hs b/src/Numeric/Matrix/Mat44d.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Matrix/Mat44d.hs @@ -0,0 +1,196 @@ +{-# LANGUAGE MagicHash #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} + +module Numeric.Matrix.Mat44d () where + +import qualified Control.Monad.ST as ST +import GHC.Exts +import Numeric.DataFrame.Internal.Array.Family.DoubleX3 +import Numeric.DataFrame.Internal.Array.Family.DoubleX4 +import qualified Numeric.DataFrame.ST as ST +import Numeric.DataFrame.SubSpace +import Numeric.DataFrame.Type +import Numeric.Matrix.Class +import Numeric.Scalar +import Numeric.Vector + +{-# INLINE mat44d #-} +mat44d :: + Scd -> Scd -> Scd -> Scd -> + Scd -> Scd -> Scd -> Scd -> + Scd -> Scd -> Scd -> Scd -> + Scd -> Scd -> Scd -> Scd -> + Mat44d +mat44d + _11 _12 _13 _14 + _21 _22 _23 _24 + _31 _32 _33 _34 + _41 _42 _43 _44 + = ST.runST $ do + df <- ST.newDataFrame + ST.writeDataFrameOff df 0 _11 + ST.writeDataFrameOff df 1 _21 + ST.writeDataFrameOff df 2 _31 + ST.writeDataFrameOff df 3 _41 + ST.writeDataFrameOff df 4 _12 + ST.writeDataFrameOff df 5 _22 + ST.writeDataFrameOff df 6 _32 + ST.writeDataFrameOff df 7 _42 + ST.writeDataFrameOff df 8 _13 + ST.writeDataFrameOff df 9 _23 + ST.writeDataFrameOff df 10 _33 + ST.writeDataFrameOff df 11 _43 + ST.writeDataFrameOff df 12 _14 + ST.writeDataFrameOff df 13 _24 + ST.writeDataFrameOff df 14 _34 + ST.writeDataFrameOff df 15 _44 + ST.unsafeFreezeDataFrame df + +instance HomTransform4 Double where + {-# INLINE translate4 #-} + translate4 (SingleFrame (DoubleX4# x# y# z# _)) = mat44d + 1 0 0 x + 0 1 0 y + 0 0 1 z + 0 0 0 1 + where + x = scalar $ D# x# + y = scalar $ D# y# + z = scalar $ D# z# + + {-# INLINE translate3 #-} + translate3 (SingleFrame (DoubleX3# x# y# z#)) = mat44d + 1 0 0 x + 0 1 0 y + 0 0 1 z + 0 0 0 1 + where + x = scalar $ D# x# + y = scalar $ D# y# + z = scalar $ D# z# + + {-# INLINE rotateX #-} + rotateX a = mat44d + 1 0 0 0 + 0 c n 0 + 0 s c 0 + 0 0 0 1 + where + c = scalar $ cos a + s = scalar $ sin a + n = -s + + {-# INLINE rotateY #-} + rotateY a = mat44d + c 0 s 0 + 0 1 0 0 + n 0 c 0 + 0 0 0 1 + where + c = scalar $ cos a + s = scalar $ sin a + n = -s + + {-# INLINE rotateZ #-} + rotateZ a = mat44d + c n 0 0 + s c 0 0 + 0 0 1 0 + 0 0 0 1 + where + c = scalar $ cos a + s = scalar $ sin a + n = -s + + {-# INLINE rotate #-} + rotate (SingleFrame (DoubleX3# x# y# z#)) a = mat44d + (c+xxv) (xyv-zs) (xzv+ys) 0 + (yxv+zs) (c+yyv) (yzv-xs) 0 + (zxv-ys) (zyv+xs) (c+zzv) 0 + 0 0 0 1 + where + c = scalar $ cos a + v = 1 - c -- v for versine + s = scalar $ sin a + x = scalar $ D# x# + y = scalar $ D# y# + z = scalar $ D# z# + xxv = x * x * v + xyv = x * y * v + xzv = x * z * v + yxv = xyv + yyv = y * y * v + yzv = y * z * v + zxv = xzv + zyv = yzv + zzv = z * z * v + xs = x * s + ys = y * s + zs = z * s + + {-# INLINE rotateEuler #-} + rotateEuler x y z = mat44d + (cy*cz) (-cy*sz) sy 0 + (cx*sz+sx*sy*cz) (cx*cz-sx*sy*sz) (-sx*cy) 0 + (sx*sz-cx*sy*cz) (sx*cz+cx*sy*sz) (cx*cy) 0 + 0 0 0 1 + where + cx = scalar $ cos x + sx = scalar $ sin x + cy = scalar $ cos y + sy = scalar $ sin y + cz = scalar $ cos z + sz = scalar $ sin z + + {-# INLINE lookAt #-} + lookAt up cam foc = mat44d + (xb!1) (xb!2) (xb!3) tx + (yb!1) (yb!2) (yb!3) ty + (zb!1) (zb!2) (zb!3) tz + 0 0 0 1 + where + zb = normalized $ cam - foc -- Basis vector for "backward", since +Z is behind the camera + xb = normalized $ up `cross` zb -- Basis vector for "right" + yb = zb `cross` xb -- Basis vector for "up" + ncam = -cam + tx = xb `dot` ncam + ty = yb `dot` ncam + tz = zb `dot` ncam + + {-# INLINE perspective #-} + perspective n f fovy aspect = mat44d + dpw 0 0 0 + 0 dph 0 0 + 0 0 a b + 0 0 (-1) 0 + where + hpd = tan (fovy * 0.5) -- height/distance + wpd = aspect * hpd; -- width/distance + dph = scalar $ 1 / hpd -- distance/height + dpw = scalar $ 1 / wpd -- distance/width + nmf = n - f + a = scalar $ (n + f) / nmf + b = scalar $ 2 * n * f / nmf + + {-# INLINE orthogonal #-} + orthogonal n f w h = mat44d + iw 0 0 0 + 0 ih 0 0 + 0 0 a b + 0 0 0 1 + where + ih = scalar $ 2 / h + iw = scalar $ 2 / w + nmf = n - f + a = scalar $ 2 / nmf + b = scalar $ (n + f) / nmf + + {-# INLINE toHomPoint #-} + toHomPoint (SingleFrame (DoubleX3# x# y# z#)) = SingleFrame (DoubleX4# x# y# z# 1.0##) + + {-# INLINE toHomVector #-} + toHomVector (SingleFrame (DoubleX3# x# y# z#)) = SingleFrame (DoubleX4# x# y# z# 0.0##) + + {-# INLINE fromHom #-} + fromHom (SingleFrame (DoubleX4# x# y# z# 0.0##)) = SingleFrame (DoubleX3# x# y# z#) + fromHom (SingleFrame (DoubleX4# x# y# z# w#)) = SingleFrame (DoubleX3# (x# /## w#) (y# /## w#) (z# /## w#)) diff --git a/src/Numeric/Matrix/Mat44f.hs b/src/Numeric/Matrix/Mat44f.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Matrix/Mat44f.hs @@ -0,0 +1,196 @@ +{-# LANGUAGE MagicHash #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} + +module Numeric.Matrix.Mat44f () where + +import qualified Control.Monad.ST as ST +import GHC.Exts +import Numeric.DataFrame.Internal.Array.Family.FloatX3 +import Numeric.DataFrame.Internal.Array.Family.FloatX4 +import qualified Numeric.DataFrame.ST as ST +import Numeric.DataFrame.SubSpace +import Numeric.DataFrame.Type +import Numeric.Matrix.Class +import Numeric.Scalar +import Numeric.Vector + +{-# INLINE mat44f #-} +mat44f :: + Scf -> Scf -> Scf -> Scf -> + Scf -> Scf -> Scf -> Scf -> + Scf -> Scf -> Scf -> Scf -> + Scf -> Scf -> Scf -> Scf -> + Mat44f +mat44f + _11 _12 _13 _14 + _21 _22 _23 _24 + _31 _32 _33 _34 + _41 _42 _43 _44 + = ST.runST $ do + df <- ST.newDataFrame + ST.writeDataFrameOff df 0 _11 + ST.writeDataFrameOff df 1 _21 + ST.writeDataFrameOff df 2 _31 + ST.writeDataFrameOff df 3 _41 + ST.writeDataFrameOff df 4 _12 + ST.writeDataFrameOff df 5 _22 + ST.writeDataFrameOff df 6 _32 + ST.writeDataFrameOff df 7 _42 + ST.writeDataFrameOff df 8 _13 + ST.writeDataFrameOff df 9 _23 + ST.writeDataFrameOff df 10 _33 + ST.writeDataFrameOff df 11 _43 + ST.writeDataFrameOff df 12 _14 + ST.writeDataFrameOff df 13 _24 + ST.writeDataFrameOff df 14 _34 + ST.writeDataFrameOff df 15 _44 + ST.unsafeFreezeDataFrame df + +instance HomTransform4 Float where + {-# INLINE translate4 #-} + translate4 (SingleFrame (FloatX4# x# y# z# _)) = mat44f + 1 0 0 x + 0 1 0 y + 0 0 1 z + 0 0 0 1 + where + x = scalar $ F# x# + y = scalar $ F# y# + z = scalar $ F# z# + + {-# INLINE translate3 #-} + translate3 (SingleFrame (FloatX3# x# y# z#)) = mat44f + 1 0 0 x + 0 1 0 y + 0 0 1 z + 0 0 0 1 + where + x = scalar $ F# x# + y = scalar $ F# y# + z = scalar $ F# z# + + {-# INLINE rotateX #-} + rotateX a = mat44f + 1 0 0 0 + 0 c n 0 + 0 s c 0 + 0 0 0 1 + where + c = scalar $ cos a + s = scalar $ sin a + n = -s + + {-# INLINE rotateY #-} + rotateY a = mat44f + c 0 s 0 + 0 1 0 0 + n 0 c 0 + 0 0 0 1 + where + c = scalar $ cos a + s = scalar $ sin a + n = -s + + {-# INLINE rotateZ #-} + rotateZ a = mat44f + c n 0 0 + s c 0 0 + 0 0 1 0 + 0 0 0 1 + where + c = scalar $ cos a + s = scalar $ sin a + n = -s + + {-# INLINE rotate #-} + rotate (SingleFrame (FloatX3# x# y# z#)) a = mat44f + (c+xxv) (xyv-zs) (xzv+ys) 0 + (yxv+zs) (c+yyv) (yzv-xs) 0 + (zxv-ys) (zyv+xs) (c+zzv) 0 + 0 0 0 1 + where + c = scalar $ cos a + v = 1 - c -- v for versine + s = scalar $ sin a + x = scalar $ F# x# + y = scalar $ F# y# + z = scalar $ F# z# + xxv = x * x * v + xyv = x * y * v + xzv = x * z * v + yxv = xyv + yyv = y * y * v + yzv = y * z * v + zxv = xzv + zyv = yzv + zzv = z * z * v + xs = x * s + ys = y * s + zs = z * s + + {-# INLINE rotateEuler #-} + rotateEuler x y z = mat44f + (cy*cz) (-cy*sz) sy 0 + (cx*sz+sx*sy*cz) (cx*cz-sx*sy*sz) (-sx*cy) 0 + (sx*sz-cx*sy*cz) (sx*cz+cx*sy*sz) (cx*cy) 0 + 0 0 0 1 + where + cx = scalar $ cos x + sx = scalar $ sin x + cy = scalar $ cos y + sy = scalar $ sin y + cz = scalar $ cos z + sz = scalar $ sin z + + {-# INLINE lookAt #-} + lookAt up cam foc = mat44f + (xb!1) (xb!2) (xb!3) tx + (yb!1) (yb!2) (yb!3) ty + (zb!1) (zb!2) (zb!3) tz + 0 0 0 1 + where + zb = normalized $ cam - foc -- Basis vector for "backward", since +Z is behind the camera + xb = normalized $ up `cross` zb -- Basis vector for "right" + yb = zb `cross` xb -- Basis vector for "up" + ncam = -cam + tx = xb `dot` ncam + ty = yb `dot` ncam + tz = zb `dot` ncam + + {-# INLINE perspective #-} + perspective n f fovy aspect = mat44f + dpw 0 0 0 + 0 dph 0 0 + 0 0 a b + 0 0 (-1) 0 + where + hpd = tan (fovy * 0.5) -- height/distance + wpd = aspect * hpd; -- width/distance + dph = scalar $ 1 / hpd -- distance/height + dpw = scalar $ 1 / wpd -- distance/width + nmf = n - f + a = scalar $ (n + f) / nmf + b = scalar $ 2 * n * f / nmf + + {-# INLINE orthogonal #-} + orthogonal n f w h = mat44f + iw 0 0 0 + 0 ih 0 0 + 0 0 a b + 0 0 0 1 + where + ih = scalar $ 2 / h + iw = scalar $ 2 / w + nmf = n - f + a = scalar $ 2 / nmf + b = scalar $ (n + f) / nmf + + {-# INLINE toHomPoint #-} + toHomPoint (SingleFrame (FloatX3# x# y# z#)) = SingleFrame (FloatX4# x# y# z# 1.0#) + + {-# INLINE toHomVector #-} + toHomVector (SingleFrame (FloatX3# x# y# z#)) = SingleFrame (FloatX4# x# y# z# 0.0#) + + {-# INLINE fromHom #-} + fromHom (SingleFrame (FloatX4# x# y# z# 0.0#)) = SingleFrame (FloatX3# x# y# z#) + fromHom (SingleFrame (FloatX4# x# y# z# w#)) = SingleFrame (FloatX3# (x# `divideFloat#` w#) (y# `divideFloat#` w#) (z# `divideFloat#` w#)) diff --git a/test/Numeric/MatrixDoubleTest.hs b/test/Numeric/MatrixDoubleTest.hs new file mode 100644 --- /dev/null +++ b/test/Numeric/MatrixDoubleTest.hs @@ -0,0 +1,212 @@ +{-# 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.DataFrame +import Numeric.DataFrame.Arbitraries () +import Numeric.DataFrame.Internal.Array.Class +import Numeric.Dimensions +import Numeric.PrimBytes +import Test.QuickCheck + +eps :: Scd +eps = 0.0000001 + +dropW :: (SubSpace t '[] '[3] '[3], SubSpace t '[] '[4] '[4]) => Vector t 4 -> Vector t 3 +dropW v | (x,y,z,_) <- unpackV4 v = vec3 x y z + +approxEq :: + forall (ds :: [Nat]). + ( + Dimensions ds, + Num (DataFrame Double ds), + PrimBytes (DataFrame Double ds), + PrimArray Double (DataFrame Double ds) + ) => + DataFrame Double ds -> DataFrame Double ds -> Bool +approxEq a b = (eps >=) . ewfoldl @_ @'[] max 0 . abs $ a - b +infix 4 `approxEq` + +prop_detTranspose :: Matrix '[Double, Double] (XN 2) (XN 2) -> Bool +prop_detTranspose (XFrame (x :*: y :*: Z)) + | -- infer KnownDim for both dimensions of matrix x (and y) + KnownDims <- dims `inSpaceOf` x + = let m = diag (ewfoldl max 0 $ abs x) + x %* transpose y + a = det m + b = det $ transpose m + in abs (a - b) / (abs a + abs b + 1) <= eps + +prop_inverse :: Matrix '[Double, Double] (XN 2) (XN 2) -> Bool +prop_inverse (XFrame (x :*: y :*: Z)) + | -- infer KnownDim for both dimensions of matrix x (and y) + (KnownDims :: Dims ns) <- dims `inSpaceOf` x + -- cumbersose inverse instance requires PrimBytes (Vector t n) + , E <- inferASing' @Double @'[Head ns] + , E <- inferPrim' @Double @'[Head ns] + = let m = diag base + x %* transpose y + mi = inverse m + err a b = ewfoldl max 0 (abs (b - a)) / base + base = ewfoldl max 0.5 (abs x) + ewfoldl max 0.5 (abs y) + in err eye (m %* mi) <= eps + && err eye (mi %* m) <= eps + +prop_LU :: Matrix '[Double, Double] (XN 2) (XN 2) -> Bool +prop_LU (XFrame (x :*: y :*: Z)) + | -- infer KnownDim for both dimensions of matrix x (and y) + (KnownDims :: Dims ns) <- dims `inSpaceOf` x + -- cumbersose inverse instance requires PrimBytes (Vector t n) + , E <- inferASing' @Double @'[Head ns] + , E <- inferPrim' @Double @'[Head ns] + = let m = diag base + x %* transpose y + f = lu m + err a b = ewfoldl max 0 (abs (b - a)) / base + base = ewfoldl max 0.5 (abs x) + ewfoldl max 0.5 (abs y) + in err (luPerm f %* m) (luLower f %* luUpper f) <= eps + +prop_translate3vs4 :: Vector Double 4 -> Bool +prop_translate3vs4 v = translate4 v == translate3 (dropW v) + +prop_translate4 :: Vector Double 4 -> Vector Double 3 -> Bool +prop_translate4 a b = translate4 a %* toHomPoint b == toHomPoint (dropW a + b) + +prop_translate3 :: Vector Double 3 -> Vector Double 3 -> Bool +prop_translate3 a b = translate3 a %* toHomPoint b == toHomPoint (a + b) + +prop_rotateX :: Vector Double 4 -> Bool +prop_rotateX v | (x,y,z,w) <- unpackV4 v = + and [ + rotateX (-2 * pi) %* v `approxEq` v, + rotateX (-1.5 * pi) %* v `approxEq` vec4 x (-z) y w, + rotateX (-pi) %* v `approxEq` vec4 x (-y) (-z) w, + rotateX (-0.5 * pi) %* v `approxEq` vec4 x z (-y) w, + rotateX 0 %* v `approxEq` v, + rotateX (0.5 * pi) %* v `approxEq` vec4 x (-z) y w, + rotateX pi %* v `approxEq` vec4 x (-y) (-z) w, + rotateX (1.5 * pi) %* v `approxEq` vec4 x z (-y) w, + rotateX (2 * pi) %* v `approxEq` v + ] + +prop_rotateY :: Vector Double 4 -> Bool +prop_rotateY v | (x,y,z,w) <- unpackV4 v = + and [ + rotateY (-2 * pi) %* v `approxEq` v, + rotateY (-1.5 * pi) %* v `approxEq` vec4 z y (-x) w, + rotateY (-pi) %* v `approxEq` vec4 (-x) y (-z) w, + rotateY (-0.5 * pi) %* v `approxEq` vec4 (-z) y x w, + rotateY 0 %* v `approxEq` v, + rotateY (0.5 * pi) %* v `approxEq` vec4 z y (-x) w, + rotateY pi %* v `approxEq` vec4 (-x) y (-z) w, + rotateY (1.5 * pi) %* v `approxEq` vec4 (-z) y x w, + rotateY (2 * pi) %* v `approxEq` v + ] + +prop_rotateZ :: Vector Double 4 -> Bool +prop_rotateZ v | (x,y,z,w) <- unpackV4 v = + and [ + rotateZ (-2 * pi) %* v `approxEq` v, + rotateZ (-1.5 * pi) %* v `approxEq` vec4 (-y) x z w, + rotateZ (-pi) %* v `approxEq` vec4 (-x) (-y) z w, + rotateZ (-0.5 * pi) %* v `approxEq` vec4 y (-x) z w, + rotateZ 0 %* v `approxEq` v, + rotateZ (0.5 * pi) %* v `approxEq` vec4 (-y) x z w, + rotateZ pi %* v `approxEq` vec4 (-x) (-y) z w, + rotateZ (1.5 * pi) %* v `approxEq` vec4 y (-x) z w, + rotateZ (2 * pi) %* v `approxEq` v + ] + +prop_rotate :: Double -> Bool +prop_rotate a = + and [ + rotate (vec3 1 0 0) a `approxEq` rotateX a, + rotate (vec3 0 1 0) a `approxEq` rotateY a, + rotate (vec3 0 0 1) a `approxEq` rotateZ a + ] + +prop_rotateEuler :: Double -> Double -> Double -> Bool +prop_rotateEuler pitch yaw roll = rotateEuler pitch yaw roll `approxEq` rotateX pitch %* rotateY yaw %* rotateZ roll + +prop_lookAt :: Vector Double 3 -> Vector Double 3 -> Vector Double 3 -> Bool +prop_lookAt up cam foc = + and [ + (normalized . fromHom $ m %* toHomPoint foc) `approxEq` vec3 0 0 (-1), + fromHom (m %* toHomPoint cam) `approxEq` 0, + fromHom (m %* toHomVector xb) `approxEq` vec3 1 0 0, + fromHom (m %* toHomVector yb) `approxEq` vec3 0 1 0, + fromHom (m %* toHomVector zb) `approxEq` vec3 0 0 1 + ] + where + m = lookAt up cam foc + zb = normalized $ cam - foc + xb = normalized $ up `cross` zb + yb = zb `cross` xb + +prop_perspective :: Double -> Double -> Double -> Double -> Bool +prop_perspective a b c d = + and [ + projectTo 0 0 n `approxEq` vec3 0 0 (-1), + projectTo 0 0 f `approxEq` vec3 0 0 1, + projectTo 1 1 n `approxEq` vec3 1 1 (-1), + projectTo 1 (-1) n `approxEq` vec3 1 (-1) (-1), + projectTo (-1) 1 n `approxEq` vec3 (-1) 1 (-1), + projectTo (-1) (-1) n `approxEq` vec3 (-1) (-1) (-1), + projectTo 1 1 f `approxEq` vec3 1 1 1, + projectTo 1 (-1) f `approxEq` vec3 1 (-1) 1, + projectTo (-1) 1 f `approxEq` vec3 (-1) 1 1, + projectTo (-1) (-1) f `approxEq` vec3 (-1) (-1) 1 + ] + where + n = 1.0 + mod' a 9.0 -- Near plane in range [1, 10) + f = n + 1.0 + mod' b 99.0 -- Far plane in range [n + 1, 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 $ m %* vec4 (x' * wpd * z) (y' * hpd * z) (-z) 1 + +prop_orthogonal :: Double -> Double -> Double -> Double -> Bool +prop_orthogonal a b c d = + and [ + projectTo 0 0 n `approxEq` vec3 0 0 (-1), + projectTo 0 0 f `approxEq` vec3 0 0 1, + projectTo 1 1 n `approxEq` vec3 1 1 (-1), + projectTo 1 (-1) n `approxEq` vec3 1 (-1) (-1), + projectTo (-1) 1 n `approxEq` vec3 (-1) 1 (-1), + projectTo (-1) (-1) n `approxEq` vec3 (-1) (-1) (-1), + projectTo 1 1 f `approxEq` vec3 1 1 1, + projectTo 1 (-1) f `approxEq` vec3 1 (-1) 1, + projectTo (-1) 1 f `approxEq` vec3 (-1) 1 1, + projectTo (-1) (-1) f `approxEq` vec3 (-1) (-1) 1 + ] + where + n = 1.0 + mod' a 9.0 -- Near plane in range [1, 10) + f = n + 1.0 + mod' b 99.0 -- Far plane in range [n + 1, 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 $ m %* vec4 (x' * w * 0.5) (y' * h * 0.5) (-z) 1 + +prop_toHomPoint :: Vector Double 3 -> Bool +prop_toHomPoint v | (x,y,z) <- unpackV3 v = toHomPoint v == vec4 x y z 1 + +prop_toHomVector :: Vector Double 3 -> Bool +prop_toHomVector v | (x,y,z) <- unpackV3 v = toHomVector v == vec4 x y z 0 + +prop_fromHom :: Vector Double 4 -> Bool +prop_fromHom v | (x,y,z,w) <- unpackV4 v = + case w of + 0 -> fromHom v == vec3 x y z + _ -> fromHom v `approxEq` vec3 (x/w) (y/w) (z/w) + +return [] +runTests :: IO Bool +runTests = $quickCheckAll diff --git a/test/Numeric/MatrixFloatTest.hs b/test/Numeric/MatrixFloatTest.hs new file mode 100644 --- /dev/null +++ b/test/Numeric/MatrixFloatTest.hs @@ -0,0 +1,212 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeOperators #-} + +module Numeric.MatrixFloatTest (runTests) where + + +import Data.Fixed +import Numeric.DataFrame +import Numeric.DataFrame.Arbitraries () +import Numeric.DataFrame.Internal.Array.Class +import Numeric.Dimensions +import Numeric.PrimBytes +import Test.QuickCheck + +eps :: Scf +eps = 0.01 + +dropW :: (SubSpace t '[] '[3] '[3], SubSpace t '[] '[4] '[4]) => Vector t 4 -> Vector t 3 +dropW v | (x,y,z,_) <- unpackV4 v = vec3 x y z + +approxEq :: + forall (ds :: [Nat]). + ( + Dimensions ds, + Num (DataFrame Float ds), + PrimBytes (DataFrame Float ds), + PrimArray Float (DataFrame Float ds) + ) => + DataFrame Float ds -> DataFrame Float ds -> Bool +approxEq a b = (eps >=) . ewfoldl @_ @'[] max 0 . abs $ a - b +infix 4 `approxEq` + +prop_detTranspose :: Matrix '[Float, Float] (XN 2) (XN 2) -> Bool +prop_detTranspose (XFrame (x :*: y :*: Z)) + | -- infer KnownDim for both dimensions of matrix x (and y) + KnownDims <- dims `inSpaceOf` x + = let m = diag (ewfoldl max 0 $ abs x) + x %* transpose y + a = det m + b = det $ transpose m + in abs (a - b) / (abs a + abs b + 1) <= eps + +prop_inverse :: Matrix '[Float, Float] (XN 2) (XN 2) -> Bool +prop_inverse (XFrame (x :*: y :*: Z)) + | -- infer KnownDim for both dimensions of matrix x (and y) + (KnownDims :: Dims ns) <- dims `inSpaceOf` x + -- cumbersose inverse instance requires PrimBytes (Vector t n) + , E <- inferASing' @Float @'[Head ns] + , E <- inferPrim' @Float @'[Head ns] + = let m = diag base + x %* transpose y + mi = inverse m + err a b = ewfoldl max 0 (abs (b - a)) / base + base = ewfoldl max 0.5 (abs x) + ewfoldl max 0.5 (abs y) + in err eye (m %* mi) <= eps + && err eye (mi %* m) <= eps + +prop_LU :: Matrix '[Float, Float] (XN 2) (XN 2) -> Bool +prop_LU (XFrame (x :*: y :*: Z)) + | -- infer KnownDim for both dimensions of matrix x (and y) + (KnownDims :: Dims ns) <- dims `inSpaceOf` x + -- cumbersose inverse instance requires PrimBytes (Vector t n) + , E <- inferASing' @Float @'[Head ns] + , E <- inferPrim' @Float @'[Head ns] + = let m = diag base + x %* transpose y + f = lu m + err a b = ewfoldl max 0 (abs (b - a)) / base + base = ewfoldl max 0.5 (abs x) + ewfoldl max 0.5 (abs y) + in err (luPerm f %* m) (luLower f %* luUpper f) <= eps + +prop_translate3vs4 :: Vector Float 4 -> Bool +prop_translate3vs4 v = translate4 v == translate3 (dropW v) + +prop_translate4 :: Vector Float 4 -> Vector Float 3 -> Bool +prop_translate4 a b = translate4 a %* toHomPoint b == toHomPoint (dropW a + b) + +prop_translate3 :: Vector Float 3 -> Vector Float 3 -> Bool +prop_translate3 a b = translate3 a %* toHomPoint b == toHomPoint (a + b) + +prop_rotateX :: Vector Float 4 -> Bool +prop_rotateX v | (x,y,z,w) <- unpackV4 v = + and [ + rotateX (-2 * pi) %* v `approxEq` v, + rotateX (-1.5 * pi) %* v `approxEq` vec4 x (-z) y w, + rotateX (-pi) %* v `approxEq` vec4 x (-y) (-z) w, + rotateX (-0.5 * pi) %* v `approxEq` vec4 x z (-y) w, + rotateX 0 %* v `approxEq` v, + rotateX (0.5 * pi) %* v `approxEq` vec4 x (-z) y w, + rotateX pi %* v `approxEq` vec4 x (-y) (-z) w, + rotateX (1.5 * pi) %* v `approxEq` vec4 x z (-y) w, + rotateX (2 * pi) %* v `approxEq` v + ] + +prop_rotateY :: Vector Float 4 -> Bool +prop_rotateY v | (x,y,z,w) <- unpackV4 v = + and [ + rotateY (-2 * pi) %* v `approxEq` v, + rotateY (-1.5 * pi) %* v `approxEq` vec4 z y (-x) w, + rotateY (-pi) %* v `approxEq` vec4 (-x) y (-z) w, + rotateY (-0.5 * pi) %* v `approxEq` vec4 (-z) y x w, + rotateY 0 %* v `approxEq` v, + rotateY (0.5 * pi) %* v `approxEq` vec4 z y (-x) w, + rotateY pi %* v `approxEq` vec4 (-x) y (-z) w, + rotateY (1.5 * pi) %* v `approxEq` vec4 (-z) y x w, + rotateY (2 * pi) %* v `approxEq` v + ] + +prop_rotateZ :: Vector Float 4 -> Bool +prop_rotateZ v | (x,y,z,w) <- unpackV4 v = + and [ + rotateZ (-2 * pi) %* v `approxEq` v, + rotateZ (-1.5 * pi) %* v `approxEq` vec4 (-y) x z w, + rotateZ (-pi) %* v `approxEq` vec4 (-x) (-y) z w, + rotateZ (-0.5 * pi) %* v `approxEq` vec4 y (-x) z w, + rotateZ 0 %* v `approxEq` v, + rotateZ (0.5 * pi) %* v `approxEq` vec4 (-y) x z w, + rotateZ pi %* v `approxEq` vec4 (-x) (-y) z w, + rotateZ (1.5 * pi) %* v `approxEq` vec4 y (-x) z w, + rotateZ (2 * pi) %* v `approxEq` v + ] + +prop_rotate :: Float -> Bool +prop_rotate a = + and [ + rotate (vec3 1 0 0) a `approxEq` rotateX a, + rotate (vec3 0 1 0) a `approxEq` rotateY a, + rotate (vec3 0 0 1) a `approxEq` rotateZ a + ] + +prop_rotateEuler :: Float -> Float -> Float -> Bool +prop_rotateEuler pitch yaw roll = rotateEuler pitch yaw roll `approxEq` rotateX pitch %* rotateY yaw %* rotateZ roll + +prop_lookAt :: Vector Float 3 -> Vector Float 3 -> Vector Float 3 -> Bool +prop_lookAt up cam foc = + and [ + (normalized . fromHom $ m %* toHomPoint foc) `approxEq` vec3 0 0 (-1), + fromHom (m %* toHomPoint cam) `approxEq` 0, + fromHom (m %* toHomVector xb) `approxEq` vec3 1 0 0, + fromHom (m %* toHomVector yb) `approxEq` vec3 0 1 0, + fromHom (m %* toHomVector zb) `approxEq` vec3 0 0 1 + ] + where + m = lookAt up cam foc + zb = normalized $ cam - foc + xb = normalized $ up `cross` zb + yb = zb `cross` xb + +prop_perspective :: Float -> Float -> Float -> Float -> Bool +prop_perspective a b c d = + and [ + projectTo 0 0 n `approxEq` vec3 0 0 (-1), + projectTo 0 0 f `approxEq` vec3 0 0 1, + projectTo 1 1 n `approxEq` vec3 1 1 (-1), + projectTo 1 (-1) n `approxEq` vec3 1 (-1) (-1), + projectTo (-1) 1 n `approxEq` vec3 (-1) 1 (-1), + projectTo (-1) (-1) n `approxEq` vec3 (-1) (-1) (-1), + projectTo 1 1 f `approxEq` vec3 1 1 1, + projectTo 1 (-1) f `approxEq` vec3 1 (-1) 1, + projectTo (-1) 1 f `approxEq` vec3 (-1) 1 1, + projectTo (-1) (-1) f `approxEq` vec3 (-1) (-1) 1 + ] + where + n = 1.0 + mod' a 9.0 -- Near plane in range [1, 10) + f = n + 1.0 + mod' b 99.0 -- Far plane in range [n + 1, 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 $ m %* vec4 (x' * wpd * z) (y' * hpd * z) (-z) 1 + +prop_orthogonal :: Float -> Float -> Float -> Float -> Bool +prop_orthogonal a b c d = + and [ + projectTo 0 0 n `approxEq` vec3 0 0 (-1), + projectTo 0 0 f `approxEq` vec3 0 0 1, + projectTo 1 1 n `approxEq` vec3 1 1 (-1), + projectTo 1 (-1) n `approxEq` vec3 1 (-1) (-1), + projectTo (-1) 1 n `approxEq` vec3 (-1) 1 (-1), + projectTo (-1) (-1) n `approxEq` vec3 (-1) (-1) (-1), + projectTo 1 1 f `approxEq` vec3 1 1 1, + projectTo 1 (-1) f `approxEq` vec3 1 (-1) 1, + projectTo (-1) 1 f `approxEq` vec3 (-1) 1 1, + projectTo (-1) (-1) f `approxEq` vec3 (-1) (-1) 1 + ] + where + n = 1.0 + mod' a 9.0 -- Near plane in range [1, 10) + f = n + 1.0 + mod' b 99.0 -- Far plane in range [n + 1, 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 $ m %* vec4 (x' * w * 0.5) (y' * h * 0.5) (-z) 1 + +prop_toHomPoint :: Vector Float 3 -> Bool +prop_toHomPoint v | (x,y,z) <- unpackV3 v = toHomPoint v == vec4 x y z 1 + +prop_toHomVector :: Vector Float 3 -> Bool +prop_toHomVector v | (x,y,z) <- unpackV3 v = toHomVector v == vec4 x y z 0 + +prop_fromHom :: Vector Float 4 -> Bool +prop_fromHom v | (x,y,z,w) <- unpackV4 v = + case w of + 0 -> fromHom v == vec3 x y z + _ -> fromHom v `approxEq` vec3 (x/w) (y/w) (z/w) + +return [] +runTests :: IO Bool +runTests = $quickCheckAll diff --git a/test/Numeric/MatrixTest.hs b/test/Numeric/MatrixTest.hs deleted file mode 100644 --- a/test/Numeric/MatrixTest.hs +++ /dev/null @@ -1,59 +0,0 @@ -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE TypeApplications #-} -{-# LANGUAGE TypeOperators #-} - -module Numeric.MatrixTest (runTests) where - - -import Numeric.DataFrame -import Numeric.DataFrame.Arbitraries () -import Numeric.Dimensions -import Test.QuickCheck - -eps :: Scd -eps = 0.0000001 - -prop_detTranspose :: Matrix '[Double, Double] (XN 2) (XN 2) -> Bool -prop_detTranspose (XFrame (x :*: y :*: Z)) - | -- infer KnownDim for both dimensions of matrix x (and y) - KnownDims <- dims `inSpaceOf` x - = let m = diag (ewfoldl max 0 $ abs x) + x %* transpose y - a = det m - b = det $ transpose m - in abs (a - b) / (abs a + abs b + 1) <= eps - -prop_inverse :: Matrix '[Double, Double] (XN 2) (XN 2) -> Bool -prop_inverse (XFrame (x :*: y :*: Z)) - | -- infer KnownDim for both dimensions of matrix x (and y) - (KnownDims :: Dims ns) <- dims `inSpaceOf` x - -- cumbersose inverse instance requires PrimBytes (Vector t n) - , E <- inferASing' @Double @'[Head ns] - , E <- inferPrim' @Double @'[Head ns] - = let m = diag base + x %* transpose y - mi = inverse m - err a b = ewfoldl max 0 (abs (b - a)) / base - base = ewfoldl max 0.5 (abs x) + ewfoldl max 0.5 (abs y) - in err eye (m %* mi) <= eps - && err eye (mi %* m) <= eps - -prop_LU :: Matrix '[Double, Double] (XN 2) (XN 2) -> Bool -prop_LU (XFrame (x :*: y :*: Z)) - | -- infer KnownDim for both dimensions of matrix x (and y) - (KnownDims :: Dims ns) <- dims `inSpaceOf` x - -- cumbersose inverse instance requires PrimBytes (Vector t n) - , E <- inferASing' @Double @'[Head ns] - , E <- inferPrim' @Double @'[Head ns] - = let m = diag base + x %* transpose y - f = lu m - err a b = ewfoldl max 0 (abs (b - a)) / base - base = ewfoldl max 0.5 (abs x) + ewfoldl max 0.5 (abs y) - in err (luPerm f %* m) (luLower f %* luUpper f) <= eps - - -return [] -runTests :: IO Bool -runTests = $quickCheckAll diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -5,7 +5,8 @@ import qualified Numeric.DataFrame.BasicTest import qualified Numeric.DataFrame.SubSpaceTest -import qualified Numeric.MatrixTest +import qualified Numeric.MatrixDoubleTest +import qualified Numeric.MatrixFloatTest import qualified Numeric.QuaternionTest @@ -14,7 +15,8 @@ tests = return [ test "DataFrame.Basic" Numeric.DataFrame.BasicTest.runTests , test "DataFrame.SubSpace" Numeric.DataFrame.SubSpaceTest.runTests - , test "Matrix" Numeric.MatrixTest.runTests + , test "MatrixDouble" Numeric.MatrixDoubleTest.runTests + , test "MatrixFloat" Numeric.MatrixFloatTest.runTests , test "Quaternion" Numeric.QuaternionTest.runTests ]