easytensor-2.1.1.0: src/Numeric/Matrix/Internal.hs
{-# LANGUAGE CPP #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE UnboxedTuples #-}
{-# LANGUAGE UndecidableInstances #-}
module Numeric.Matrix.Internal
( MatrixTranspose (..)
, SquareMatrix (..)
, MatrixDeterminant (..)
, MatrixInverse (..)
, Matrix
, HomTransform4 (..)
, Mat22f, Mat23f, Mat24f
, Mat32f, Mat33f, Mat34f
, Mat42f, Mat43f, Mat44f
, Mat22d, Mat23d, Mat24d
, Mat32d, Mat33d, Mat34d
, Mat42d, Mat43d, Mat44d
, mat22, mat33, mat44
, (%*)
) where
import GHC.Base
import Numeric.DataFrame.Contraction ((%*))
import Numeric.DataFrame.Internal.PrimArray
import Numeric.DataFrame.Type
import Numeric.Dimensions
import Numeric.Scalar.Internal
import Numeric.Vector.Internal
-- | Alias for DataFrames of rank 2
type Matrix (t :: l) (n :: k) (m :: k) = DataFrame t '[n,m]
class MatrixTranspose t (n :: k) (m :: k) where
-- | Transpose Mat
transpose :: Matrix t n m -> Matrix t m n
class SquareMatrix t (n :: Nat) where
-- | Mat with 1 on diagonal and 0 elsewhere
eye :: Matrix t n n
-- | Put the same value on the Mat diagonal, 0 otherwise
diag :: Scalar t -> Matrix t n n
-- | Sum of diagonal elements
trace :: Matrix t n n -> Scalar t
class MatrixDeterminant t (n :: Nat) where
-- | Determinant of Mat
det :: Matrix t n n -> Scalar t
class MatrixInverse t (n :: Nat) where
-- | Matrix inverse
inverse :: Matrix t n n -> Matrix t n n
-- | Operations on 4x4 transformation matrices and vectors in homogeneous coordinates.
-- All angles are specified in radians.
--
-- Note: since version 2 of @easytensor@, DataFrames and matrices are row-major.
-- A good SIMD implementation may drastically improve performance
-- of 4D vector-matrix products of the form @v %* m@, but not so much
-- for products of the form @m %* v@ (due to memory layout).
-- Thus, all operations here assume the former form to benefit more from
-- SIMD in future.
class HomTransform4 t where
-- | Create a translation matrix from a vector. The 4th coordinate is ignored.
--
-- If @p ! 3 == 1@ and @v ! 3 == 0@, then
--
-- > p %* translate4 v == p + v
--
translate4 :: Vector t 4 -> Matrix t 4 4
-- | Create a translation matrix from a vector.
--
-- If @p ! 3 == 1@, then
--
-- > p %* translate3 v == p + toHomVector v
--
translate3 :: Vector t 3 -> Matrix t 4 4
-- | Rotation matrix for a rotation around the X axis, angle is given in radians.
-- e.g. @p %* rotateX (pi/2)@ rotates point @p@ around @Ox@ by 90 degrees.
rotateX :: t -> Matrix t 4 4
-- | Rotation matrix for a rotation around the Y axis, angle is given in radians.
-- e.g. @p %* rotateY (pi/2)@ rotates point @p@ around @Oy@ by 90 degrees.
rotateY :: t -> Matrix t 4 4
-- | Rotation matrix for a rotation around the Z axis, angle is given in radians.
-- e.g. @p %* rotateZ (pi/2)@ rotates point @p@ around @Oz@ by 90 degrees.
rotateZ :: t -> Matrix t 4 4
-- | Rotation matrix for a rotation around an arbitrary normalized vector
-- e.g. @p %* rotate (pi/2) v@ rotates point @p@ around @v@ by 90 degrees.
rotate :: Vector t 3 -> 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 == rotateZ roll %* rotateY yaw %* rotateX pitch
--
-- 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
-> Vector t 3 -- ^ The viewers position
-> Vector t 3 -- ^ The point to look at
-> Matrix t 4 4
-- | A perspective symmetric projection matrix. Right-handed coordinate system. (@x@ - right, @y@ - top)
-- http://en.wikibooks.org/wiki/GLSL_Programming/Vertex_Transformations
perspective :: t -- ^ Near plane clipping distance (always positive)
-> t -- ^ Far plane clipping distance (always positive)
-> t -- ^ Field of view of the y axis, in radians
-> t -- ^ Aspect ratio, i.e. screen's width\/height
-> Matrix t 4 4
-- | An orthogonal symmetric projection matrix. Right-handed coordinate system. (@x@ - right, @y@ - top)
-- http://en.wikibooks.org/wiki/GLSL_Programming/Vertex_Transformations
orthogonal :: t -- ^ Near plane clipping distance
-> t -- ^ Far plane clipping distance
-> t -- ^ width
-> t -- ^ height
-> Matrix t 4 4
-- | Add one more dimension and set it to 1.
toHomPoint :: Vector t 3 -> Vector t 4
-- | Add one more dimension and set it to 0.
toHomVector :: Vector t 3 -> Vector t 4
-- | Transform a homogenous vector or point into a normal 3D vector.
-- If the last coordinate is not zero, divide the rest by it.
fromHom :: Vector t 4 -> Vector t 3
-- Type abbreviations
type Mat22f = Matrix Float 2 2
type Mat32f = Matrix Float 3 2
type Mat42f = Matrix Float 4 2
type Mat23f = Matrix Float 2 3
type Mat33f = Matrix Float 3 3
type Mat43f = Matrix Float 4 3
type Mat24f = Matrix Float 2 4
type Mat34f = Matrix Float 3 4
type Mat44f = Matrix Float 4 4
type Mat22d = Matrix Double 2 2
type Mat32d = Matrix Double 3 2
type Mat42d = Matrix Double 4 2
type Mat23d = Matrix Double 2 3
type Mat33d = Matrix Double 3 3
type Mat43d = Matrix Double 4 3
type Mat24d = Matrix Double 2 4
type Mat34d = Matrix Double 3 4
type Mat44d = Matrix Double 4 4
-- | Compose a 2x2D matrix
mat22 :: PrimBytes (t :: Type)
=> Vector t 2 -> Vector t 2 -> Matrix t 2 2
mat22 = packDF @_ @2 @'[2]
{-# INLINE mat22 #-}
-- | Compose a 3x3D matrix
mat33 :: PrimBytes (t :: Type)
=> Vector t 3 -> Vector t 3 -> Vector t 3 -> Matrix t 3 3
mat33 = packDF @_ @3 @'[3]
{-# INLINE mat33 #-}
-- | Compose a 4x4D matrix
mat44 :: PrimBytes (t :: Type)
=> Vector t 4
-> Vector t 4
-> Vector t 4
-> Vector t 4
-> Matrix t 4 4
mat44 = packDF @_ @4 @'[4]
{-# INLINE mat44 #-}
instance ( KnownDim n, KnownDim m
, PrimArray t (Matrix t n m)
, PrimArray t (Matrix t m n)
) => MatrixTranspose t (n :: Nat) (m :: Nat) where
transpose df = case uniqueOrCumulDims df of
Left a -> broadcast a
Right _
| wm <- dimVal' @m
, wn <- dimVal' @n
, m <- case wm of W# w -> word2Int# w
, n <- case wn of W# w -> word2Int# w
-> let f ( I# i, I# j )
| isTrue# (i ==# n) = f ( 0, I# (j +# 1#) )
| otherwise = (# ( I# (i +# 1#), I# j )
, ix# (i *# m +# j) df
#)
in case gen# (CumulDims [wm*wn, wn, 1]) f (0,0) of (# _, r #) -> r
instance MatrixTranspose (t :: Type) (xn :: XNat) (xm :: XNat) where
transpose (XFrame (df :: DataFrame t ns))
| ((D :: Dim n) :* (D :: Dim m) :* U) <- dims @ns
, Dict <- inferPrimElem df
= XFrame (transpose df :: Matrix t m n)
#if !MIN_VERSION_GLASGOW_HASKELL(8,10,0,0)
transpose _ = error "MatrixTranspose/transpose: impossible argument"
#endif
instance (KnownDim n, PrimArray t (Matrix t n n), Num t)
=> SquareMatrix t n where
eye
| n <- dimVal' @n
= let f 0 = (# n , 1 #)
f k = (# k - 1, 0 #)
in case gen# (CumulDims [n*n, n, 1]) f 0 of
(# _, r #) -> r
diag se
| n <- dimVal' @n
, e <- unScalar se
= let f 0 = (# n , e #)
f k = (# k - 1, 0 #)
in case gen# (CumulDims [n*n, n, 1]) f 0 of
(# _, r #) -> r
trace df
| I# n <- fromIntegral $ dimVal' @n
, n1 <- n +# 1#
= let f 0# = ix# 0# df
f k = ix# k df + f (k -# n1)
in scalar $ f (n *# n -# 1#)