easytensor 1.0.0.1 → 1.0.1.0

raw patch · 10 files changed

+919/−136 lines, 10 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• easytensor.cabal +18/−3
• src/Numeric/DataFrame/Internal/Array/Family/ArrayBase.hs +66/−68
• src/Numeric/Matrix.hs +3/−0
• src/Numeric/Matrix/Class.hs +12/−4
• src/Numeric/Matrix/Mat44d.hs +196/−0
• src/Numeric/Matrix/Mat44f.hs +196/−0
• test/Numeric/MatrixDoubleTest.hs +212/−0
• test/Numeric/MatrixFloatTest.hs +212/−0
• test/Numeric/MatrixTest.hs +0/−59
• test/Spec.hs +4/−2

@@ -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

@@ -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 #-} 

@@ -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

@@ -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

@@ -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#))

@@ -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#))

@@ -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

@@ -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

@@ -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

@@ -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   ]