diff --git a/bench/misc.hs b/bench/misc.hs --- a/bench/misc.hs +++ b/bench/misc.hs @@ -9,8 +9,12 @@ import Numeric.DataFrame import Numeric.Dimensions +import qualified Control.Monad.ST as ST +import qualified Numeric.DataFrame.ST as ST +-- import qualified Numeric.Dimensions.Traverse.ST as ST + main :: IO () main = do putStrLn "Hello world!" @@ -75,3 +79,19 @@ -- Updating existing frames print $ update (2:!Z) (scalar 777) rVec print $ update (2:!3:!Z) (vec2 999 999) x + + let matX = iwgen (scalar . fromEnum) :: DataFrame Int '[2,5,4] + matY = iwgen (scalar . fromEnum) :: DataFrame Int '[5,4] + putStrLn "Check carefully that this returns no garbage" + print matX + print (ewmap (<+:> scalar 111) matX :: DataFrame Int '[3,5,4]) + print matY + print (ewmap fromScalar matY :: DataFrame Int '[3,5,4]) + + -- Working with mutable frames + print $ ST.runST $ do + sdf <- ST.thawDataFrame matY + ST.writeDataFrame sdf (1:!1:!Z) 900101 + ST.writeDataFrame sdf (3:!3:!Z) 900303 + ST.writeDataFrame sdf (5:!3:!Z) 900503 + ST.unsafeFreezeDataFrame sdf diff --git a/easytensor.cabal b/easytensor.cabal --- a/easytensor.cabal +++ b/easytensor.cabal @@ -1,8 +1,8 @@ name: easytensor -version: 0.3.0.0 +version: 0.3.1.0 cabal-version: >=1.20 build-type: Simple -license: MIT +license: BSD3 license-file: LICENSE copyright: (c) Artem Chirkin maintainer: chirkin@arch.ethz.ch @@ -42,6 +42,7 @@ Numeric.Matrix Numeric.Vector Numeric.Scalar + Numeric.Quaternion build-depends: base >=4.9 && <5, ghc-prim >=0.5, @@ -65,7 +66,12 @@ Numeric.DataFrame.Type Numeric.DataFrame.Inference Numeric.DataFrame.Shape - Numeric.Matrix.Type + Numeric.Matrix.Class + Numeric.Matrix.Mat44d + Numeric.Matrix.Mat44f + Numeric.Quaternion.Class + Numeric.Quaternion.QDouble + Numeric.Quaternion.QFloat if impl(ghcjs) other-modules: Numeric.Array.Family.ArrayT @@ -86,18 +92,24 @@ Numeric.Array.Family.FloatX2 Numeric.Array.Family.FloatX3 Numeric.Array.Family.FloatX4 + Numeric.Array.Family.DoubleX2 + Numeric.Array.Family.DoubleX3 + Numeric.Array.Family.DoubleX4 js-sources: src-ghcjs/Numeric/Array/Family/ArrayT.js - ghc-options: -Wall -fwarn-tabs -O2 + src-ghcjs/Numeric/Quaternion/Quaternion.js + src-ghcjs/Numeric/Matrix/Mat44.js + ghc-options: -Wall -fwarn-tabs -fwarn-unused-do-bind -fwarn-monomorphism-restriction -O2 test-suite et-test - type: detailed-0.9 - test-module: Spec + type: exitcode-stdio-1.0 + main-is: Spec.hs other-modules: Numeric.DataFrame.Arbitraries Numeric.DataFrame.SubSpaceTest Numeric.DataFrame.BasicTest + Numeric.QuaternionTest build-depends: base -any, Cabal >=1.20, diff --git a/src-base/Numeric/Array.hs b/src-base/Numeric/Array.hs --- a/src-base/Numeric/Array.hs +++ b/src-base/Numeric/Array.hs @@ -31,3 +31,7 @@ import Numeric.Array.Family.FloatX2 () import Numeric.Array.Family.FloatX3 () import Numeric.Array.Family.FloatX4 () + +import Numeric.Array.Family.DoubleX2 () +import Numeric.Array.Family.DoubleX3 () +import Numeric.Array.Family.DoubleX4 () diff --git a/src-base/Numeric/Array/Family.hs b/src-base/Numeric/Array/Family.hs --- a/src-base/Numeric/Array/Family.hs +++ b/src-base/Numeric/Array/Family.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE CPP #-} {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE ExistentialQuantification #-} @@ -31,17 +32,23 @@ , ArrayW (..), ArrayW8 (..), ArrayW16 (..), ArrayW32 (..), ArrayW64 (..) , Scalar (..) , FloatX2 (..), FloatX3 (..), FloatX4 (..) + , DoubleX2 (..), DoubleX3 (..), DoubleX4 (..) , ArrayInstanceInference, ElemType (..), ArraySize (..) , ElemTypeInference (..), ArraySizeInference (..), ArrayInstanceEvidence , getArrayInstance, ArrayInstance (..), inferArrayInstance ) where +#include "MachDeps.h" import Data.Int (Int16, Int32, Int64, Int8) import Data.Type.Equality ((:~:) (..)) import Data.Word (Word16, Word32, Word64, Word8) -import GHC.Prim (ByteArray#, Double#, Float#, Int#, - Word#, unsafeCoerce#) +import GHC.Prim ( ByteArray#, Double#, Float# +#if WORD_SIZE_IN_BITS < 64 + , Int64#, Word64# +#endif + , Int#, Word#, unsafeCoerce#) +import GHC.Exts (RuntimeRep(..)) import Numeric.Array.ElementWise import Numeric.Commons @@ -54,6 +61,9 @@ Array Float '[2] = FloatX2 Array Float '[3] = FloatX3 Array Float '[4] = FloatX4 + Array Double '[2] = DoubleX2 + Array Double '[3] = DoubleX3 + Array Double '[4] = DoubleX4 Array Float (d ': ds) = ArrayF (d ': ds) Array Double (d ': ds) = ArrayD (d ': ds) Array Int (d ': ds) = ArrayI (d ': ds) @@ -75,19 +85,47 @@ instance Show t => Show (Scalar t) where show (Scalar t) = "{ " ++ show t ++ " }" -type instance ElemRep (Scalar t) = ElemRep t +type instance ElemRep (Scalar Float ) = 'FloatRep +type instance ElemRep (Scalar Double) = 'DoubleRep +type instance ElemRep (Scalar Int ) = 'IntRep +type instance ElemRep (Scalar Int8 ) = 'IntRep +type instance ElemRep (Scalar Int16 ) = 'IntRep +type instance ElemRep (Scalar Int32 ) = 'IntRep +#if WORD_SIZE_IN_BITS < 64 +type instance ElemRep (Scalar Int64 ) = 'Int64Rep +#else +type instance ElemRep (Scalar Int64 ) = 'IntRep +#endif +type instance ElemRep (Scalar Word ) = 'WordRep +type instance ElemRep (Scalar Word8 ) = 'WordRep +type instance ElemRep (Scalar Word16) = 'WordRep +type instance ElemRep (Scalar Word32) = 'WordRep +#if WORD_SIZE_IN_BITS < 64 +type instance ElemRep (Scalar Word64) = 'Word64Rep +#else +type instance ElemRep (Scalar Word64) = 'WordRep +#endif + type instance ElemPrim (Scalar Float ) = Float# type instance ElemPrim (Scalar Double) = Double# type instance ElemPrim (Scalar Int ) = Int# type instance ElemPrim (Scalar Int8 ) = Int# type instance ElemPrim (Scalar Int16 ) = Int# type instance ElemPrim (Scalar Int32 ) = Int# +#if WORD_SIZE_IN_BITS < 64 +type instance ElemPrim (Scalar Int64 ) = Int64# +#else type instance ElemPrim (Scalar Int64 ) = Int# +#endif type instance ElemPrim (Scalar Word ) = Word# type instance ElemPrim (Scalar Word8 ) = Word# type instance ElemPrim (Scalar Word16) = Word# type instance ElemPrim (Scalar Word32) = Word# +#if WORD_SIZE_IN_BITS < 64 +type instance ElemPrim (Scalar Word64) = Word64# +#else type instance ElemPrim (Scalar Word64) = Word# +#endif deriving instance PrimBytes (Scalar Float) deriving instance PrimBytes (Scalar Double) @@ -146,8 +184,13 @@ | FromScalarI16# Int# data ArrayI32 (ds :: [Nat]) = ArrayI32# Int# Int# ByteArray# | FromScalarI32# Int# +#if WORD_SIZE_IN_BITS < 64 data ArrayI64 (ds :: [Nat]) = ArrayI64# Int# Int# ByteArray# + | FromScalarI64# Int64# +#else +data ArrayI64 (ds :: [Nat]) = ArrayI64# Int# Int# ByteArray# | FromScalarI64# Int# +#endif data ArrayW (ds :: [Nat]) = ArrayW# Int# Int# ByteArray# | FromScalarW# Word# data ArrayW8 (ds :: [Nat]) = ArrayW8# Int# Int# ByteArray# @@ -156,8 +199,13 @@ | FromScalarW16# Word# data ArrayW32 (ds :: [Nat]) = ArrayW32# Int# Int# ByteArray# | FromScalarW32# Word# +#if WORD_SIZE_IN_BITS < 64 data ArrayW64 (ds :: [Nat]) = ArrayW64# Int# Int# ByteArray# + | FromScalarW64# Word64# +#else +data ArrayW64 (ds :: [Nat]) = ArrayW64# Int# Int# ByteArray# | FromScalarW64# Word# +#endif -- * Specialized types -- More efficient data types for small fixed-size tensors @@ -165,6 +213,10 @@ data FloatX3 = FloatX3# Float# Float# Float# data FloatX4 = FloatX4# Float# Float# Float# Float# +data DoubleX2 = DoubleX2# Double# Double# +data DoubleX3 = DoubleX3# Double# Double# Double# +data DoubleX4 = DoubleX4# Double# Double# Double# Double# + -- * Recovering type instances at runtime -- A combination of `ElemType t` and `ArraySize ds` should -- define an instance of `Array t ds` unambiguously. @@ -218,6 +270,9 @@ | ( Array t ds ~ FloatX2, ds ~ '[2], t ~ Float) => AIFloatX2 | ( Array t ds ~ FloatX3, ds ~ '[3], t ~ Float) => AIFloatX3 | ( Array t ds ~ FloatX4, ds ~ '[4], t ~ Float) => AIFloatX4 + | ( Array t ds ~ DoubleX2, ds ~ '[2], t ~ Double) => AIDoubleX2 + | ( Array t ds ~ DoubleX3, ds ~ '[3], t ~ Double) => AIDoubleX3 + | ( Array t ds ~ DoubleX4, ds ~ '[4], t ~ Double) => AIDoubleX4 -- | A singleton type used to prove that the given Array family instance -- has a known instance @@ -346,7 +401,7 @@ (ETWord64 , ASScalar) -> AIScalar (ETFloat , ASX2) -> AIFloatX2 - (ETDouble , ASX2) -> AIArrayD + (ETDouble , ASX2) -> AIDoubleX2 (ETInt , ASX2) -> AIArrayI (ETInt8 , ASX2) -> AIArrayI8 (ETInt16 , ASX2) -> AIArrayI16 @@ -359,7 +414,7 @@ (ETWord64 , ASX2) -> AIArrayW64 (ETFloat , ASX3) -> AIFloatX3 - (ETDouble , ASX3) -> AIArrayD + (ETDouble , ASX3) -> AIDoubleX3 (ETInt , ASX3) -> AIArrayI (ETInt8 , ASX3) -> AIArrayI8 (ETInt16 , ASX3) -> AIArrayI16 @@ -372,7 +427,7 @@ (ETWord64 , ASX3) -> AIArrayW64 (ETFloat , ASX4) -> AIFloatX4 - (ETDouble , ASX4) -> AIArrayD + (ETDouble , ASX4) -> AIDoubleX4 (ETInt , ASX4) -> AIArrayI (ETInt8 , ASX4) -> AIArrayI8 (ETInt16 , ASX4) -> AIArrayI16 @@ -385,7 +440,7 @@ (ETWord64 , ASX4) -> AIArrayW64 (ETFloat , ASXN) -> unsafeCoerce# (AIArrayF :: ArrayInstance Float '[5]) - (ETDouble , ASXN) -> AIArrayD + (ETDouble , ASXN) -> unsafeCoerce# (AIArrayD :: ArrayInstance Double '[5]) (ETInt , ASXN) -> AIArrayI (ETInt8 , ASXN) -> AIArrayI8 (ETInt16 , ASXN) -> AIArrayI16 diff --git a/src-base/Numeric/Array/Family/ArrayD.hs b/src-base/Numeric/Array/Family/ArrayD.hs --- a/src-base/Numeric/Array/Family/ArrayD.hs +++ b/src-base/Numeric/Array/Family/ArrayD.hs @@ -40,7 +40,7 @@ import Numeric.Dimensions import Numeric.Dimensions.Traverse import Numeric.TypeLits -import Numeric.Matrix.Type +import Numeric.Matrix.Class #include "MachDeps.h" diff --git a/src-base/Numeric/Array/Family/ArrayF.hs b/src-base/Numeric/Array/Family/ArrayF.hs --- a/src-base/Numeric/Array/Family/ArrayF.hs +++ b/src-base/Numeric/Array/Family/ArrayF.hs @@ -42,7 +42,7 @@ import Numeric.Dimensions import Numeric.Dimensions.Traverse import Numeric.TypeLits -import Numeric.Matrix.Type +import Numeric.Matrix.Class #include "MachDeps.h" #define ARR_TYPE ArrayF diff --git a/src-base/Numeric/Array/Family/ArrayI64.hs b/src-base/Numeric/Array/Family/ArrayI64.hs --- a/src-base/Numeric/Array/Family/ArrayI64.hs +++ b/src-base/Numeric/Array/Family/ArrayI64.hs @@ -27,10 +27,14 @@ module Numeric.Array.Family.ArrayI64 () where +#include "MachDeps.h" import GHC.Base (runRW#) import GHC.Prim import GHC.Types (Int (..), RuntimeRep (..), isTrue#) import GHC.Int (Int64 (..)) +#if WORD_SIZE_IN_BITS < 64 +import GHC.IntWord64 +#endif import Numeric.Array.ElementWise import Numeric.Array.Family @@ -39,11 +43,37 @@ import Numeric.Dimensions.Traverse -#include "MachDeps.h" + +#if SIZEOF_HSWORD < 8 #define ARR_TYPE ArrayI64 #define ARR_FROMSCALAR FromScalarI64# #define ARR_CONSTR ArrayI64# #define EL_TYPE_BOXED Int64 +#define EL_TYPE_PRIM Int64# +#define EL_RUNTIME_REP 'Int64Rep +#define EL_CONSTR I64# +#define EL_SIZE SIZEOF_INT64# +#define EL_ALIGNMENT ALIGNMENT_INT64# +#define EL_ZERO (intToInt64# 0#) +#define EL_ONE (intToInt64# 1#) +#define EL_MINUS_ONE (intToInt64# (-1#)) +#define INDEX_ARRAY indexInt64Array# +#define WRITE_ARRAY writeInt64Array# +#define OP_EQ (eqInt64#) +#define OP_NE (neInt64#) +#define OP_GT (gtInt64#) +#define OP_GE (geInt64#) +#define OP_LT (ltInt64#) +#define OP_LE (leInt64#) +#define OP_PLUS (plusInt64#) +#define OP_MINUS (minusInt64#) +#define OP_TIMES (timesInt64#) +#define OP_NEGATE negateInt64# +#else +#define ARR_TYPE ArrayI64 +#define ARR_FROMSCALAR FromScalarI64# +#define ARR_CONSTR ArrayI64# +#define EL_TYPE_BOXED Int64 #define EL_TYPE_PRIM Int# #define EL_RUNTIME_REP 'IntRep #define EL_CONSTR I64# @@ -64,28 +94,29 @@ #define OP_MINUS (-#) #define OP_TIMES (*#) #define OP_NEGATE negateInt# +#endif #include "Array.h" instance Num (ArrayI64 ds) where - (+) = zipV (+#) + (+) = zipV OP_PLUS {-# INLINE (+) #-} - (-) = zipV (-#) + (-) = zipV OP_MINUS {-# INLINE (-) #-} - (*) = zipV (*#) + (*) = zipV OP_TIMES {-# INLINE (*) #-} - negate = mapV negateInt# + negate = mapV OP_NEGATE {-# INLINE negate #-} - abs = mapV (\x -> if isTrue# (x >=# 0#) + abs = mapV (\x -> if isTrue# (OP_GE x EL_ZERO) then x - else negateInt# x + else OP_NEGATE x ) {-# INLINE abs #-} - signum = mapV (\x -> if isTrue# (x ># 0#) - then 1# - else if isTrue# (x <# 0#) - then -1# - else 0# + signum = mapV (\x -> if isTrue# (OP_GT x EL_ZERO) + then EL_ONE + else if isTrue# (OP_LT x EL_ZERO) + then EL_MINUS_ONE + else EL_ZERO ) {-# INLINE signum #-} fromInteger = broadcastArray . fromInteger diff --git a/src-base/Numeric/Array/Family/ArrayW64.hs b/src-base/Numeric/Array/Family/ArrayW64.hs --- a/src-base/Numeric/Array/Family/ArrayW64.hs +++ b/src-base/Numeric/Array/Family/ArrayW64.hs @@ -27,10 +27,14 @@ module Numeric.Array.Family.ArrayW64 () where +#include "MachDeps.h" import GHC.Base (runRW#) import GHC.Prim import GHC.Types (Int (..), RuntimeRep (..), isTrue#) import GHC.Word (Word64 (..)) +#if WORD_SIZE_IN_BITS < 64 +import GHC.IntWord64 +#endif import Numeric.Array.ElementWise import Numeric.Array.Family @@ -38,12 +42,42 @@ import Numeric.Dimensions import Numeric.Dimensions.Traverse +#if SIZEOF_HSWORD < 8 -#include "MachDeps.h" +plusWord64#, minusWord64#, timesWord64# + :: Word64# -> Word64# -> Word64# +plusWord64# x y = int64ToWord64# (word64ToInt64# x `plusInt64#` word64ToInt64# y) +minusWord64# x y = int64ToWord64# (word64ToInt64# x `minusInt64#` word64ToInt64# y) +timesWord64# x y = int64ToWord64# (word64ToInt64# x `timesInt64#` word64ToInt64# y) + + #define ARR_TYPE ArrayW64 #define ARR_FROMSCALAR FromScalarW64# #define ARR_CONSTR ArrayW64# #define EL_TYPE_BOXED Word64 +#define EL_TYPE_PRIM Word64# +#define EL_RUNTIME_REP 'Word64Rep +#define EL_CONSTR W64# +#define EL_SIZE SIZEOF_WORD64# +#define EL_ALIGNMENT ALIGNMENT_WORD64# +#define EL_ZERO (wordToWord64# 0##) +#define EL_ONE (wordToWord64# 1##) +#define INDEX_ARRAY indexWord64Array# +#define WRITE_ARRAY writeWord64Array# +#define OP_EQ (eqWord64#) +#define OP_NE (neWord64#) +#define OP_GT (gtWord64#) +#define OP_GE (geWord64#) +#define OP_LT (ltWord64#) +#define OP_LE (leWord64#) +#define OP_PLUS (plusWord64#) +#define OP_MINUS (minusWord64#) +#define OP_TIMES (timesWord64#) +#else +#define ARR_TYPE ArrayW64 +#define ARR_FROMSCALAR FromScalarW64# +#define ARR_CONSTR ArrayW64# +#define EL_TYPE_BOXED Word64 #define EL_TYPE_PRIM Word# #define EL_RUNTIME_REP 'WordRep #define EL_CONSTR W64# @@ -51,7 +85,6 @@ #define EL_ALIGNMENT ALIGNMENT_WORD64# #define EL_ZERO 0## #define EL_ONE 1## -#define EL_MINUS_ONE -1# #define INDEX_ARRAY indexWord64Array# #define WRITE_ARRAY writeWord64Array# #define OP_EQ eqWord# @@ -63,22 +96,27 @@ #define OP_PLUS plusWord# #define OP_MINUS minusWord# #define OP_TIMES timesWord# +#endif #include "Array.h" instance Num (ArrayW64 ds) where - (+) = zipV plusWord# + (+) = zipV OP_PLUS {-# INLINE (+) #-} - (-) = zipV minusWord# + (-) = zipV OP_MINUS {-# INLINE (-) #-} - (*) = zipV timesWord# + (*) = zipV OP_TIMES {-# INLINE (*) #-} +#if SIZEOF_HSWORD < 8 + negate = mapV (\x -> int64ToWord64# (negateInt64# (word64ToInt64# x))) +#else negate = mapV (\x -> int2Word# (negateInt# (word2Int# x))) +#endif {-# INLINE negate #-} abs = id {-# INLINE abs #-} - signum = mapV (\x -> if isTrue# (gtWord# x 0##) - then 1## - else 0## + signum = mapV (\x -> if isTrue# (OP_GT x EL_ZERO) + then EL_ONE + else EL_ZERO ) {-# INLINE signum #-} fromInteger = broadcastArray . fromInteger diff --git a/src-base/Numeric/Array/Family/DoubleX2.hs b/src-base/Numeric/Array/Family/DoubleX2.hs new file mode 100644 --- /dev/null +++ b/src-base/Numeric/Array/Family/DoubleX2.hs @@ -0,0 +1,262 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UnboxedTuples #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Array.Family.DoubleX2 +-- Copyright : (c) Artem Chirkin +-- License : BSD3 +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Array.Family.DoubleX2 () where + + +#include "MachDeps.h" + +import GHC.Base (runRW#) +import GHC.Prim +import GHC.Types (Double (..), RuntimeRep (..), + isTrue#) + +import Numeric.Array.ElementWise +import Numeric.Array.Family +import Numeric.Commons +import Numeric.Dimensions + + + + + +instance Show DoubleX2 where + show (DoubleX2# a1 a2) = "{ " ++ show (D# a1) + ++ ", " ++ show (D# a2) + ++ " }" + + + +instance Eq DoubleX2 where + DoubleX2# a1 a2 == DoubleX2# b1 b2 = isTrue# ( (a1 ==## b1) + `andI#` (a2 ==## b2) + ) + {-# INLINE (==) #-} + DoubleX2# a1 a2 /= DoubleX2# b1 b2 = isTrue# ( (a1 /=## b1) + `orI#` (a2 /=## b2) + ) + {-# INLINE (/=) #-} + + + +-- | Implement partial ordering for `>`, `<`, `>=`, `<=` +-- and lexicographical ordering for `compare` +instance Ord DoubleX2 where + DoubleX2# a1 a2 > DoubleX2# b1 b2 = isTrue# ( (a1 >## b1) + `andI#` (a2 >## b2) + ) + {-# INLINE (>) #-} + DoubleX2# a1 a2 < DoubleX2# b1 b2 = isTrue# ( (a1 <## b1) + `andI#` (a2 <## b2) + ) + {-# INLINE (<) #-} + DoubleX2# a1 a2 >= DoubleX2# b1 b2 = isTrue# ( (a1 >=## b1) + `andI#` (a2 >=## b2) + ) + {-# INLINE (>=) #-} + DoubleX2# a1 a2 <= DoubleX2# b1 b2 = isTrue# ( (a1 <=## b1) + `andI#` (a2 <=## b2) + ) + {-# INLINE (<=) #-} + -- | Compare lexicographically + compare (DoubleX2# a1 a2) (DoubleX2# b1 b2) + | isTrue# (a1 >## b1) = GT + | isTrue# (a1 <## b1) = LT + | isTrue# (a2 >## b2) = GT + | isTrue# (a2 <## b2) = LT + | otherwise = EQ + {-# INLINE compare #-} + -- | Element-wise minimum + min (DoubleX2# a1 a2) (DoubleX2# b1 b2) = + DoubleX2# (if isTrue# (a1 >## b1) then b1 else a1) + (if isTrue# (a2 >## b2) then b2 else a2) + {-# INLINE min #-} + -- | Element-wise maximum + max (DoubleX2# a1 a2) (DoubleX2# b1 b2) = + DoubleX2# (if isTrue# (a1 >## b1) then a1 else b1) + (if isTrue# (a2 >## b2) then a2 else b2) + {-# INLINE max #-} + + + +-- | element-wise operations for vectors +instance Num DoubleX2 where + DoubleX2# a1 a2 + DoubleX2# b1 b2 + = DoubleX2# ((+##) a1 b1) ((+##) a2 b2) + {-# INLINE (+) #-} + DoubleX2# a1 a2 - DoubleX2# b1 b2 + = DoubleX2# ((-##) a1 b1) ((-##) a2 b2) + {-# INLINE (-) #-} + DoubleX2# a1 a2 * DoubleX2# b1 b2 + = DoubleX2# ((*##) a1 b1) ((*##) a2 b2) + {-# INLINE (*) #-} + negate (DoubleX2# a1 a2) + = DoubleX2# (negateDouble# a1) (negateDouble# a2) + {-# INLINE negate #-} + abs (DoubleX2# a1 a2) + = DoubleX2# (if isTrue# (a1 >=## 0.0##) then a1 else negateDouble# a1) + (if isTrue# (a2 >=## 0.0##) then a2 else negateDouble# a2) + {-# INLINE abs #-} + signum (DoubleX2# a1 a2) + = DoubleX2# (if isTrue# (a1 >## 0.0##) + then 1.0## + else if isTrue# (a1 <## 0.0##) then -1.0## else 0.0## ) + (if isTrue# (a2 >## 0.0##) + then 1.0## + else if isTrue# (a2 <## 0.0##) then -1.0## else 0.0## ) + {-# INLINE signum #-} + fromInteger n = case fromInteger n of D# x -> DoubleX2# x x + {-# INLINE fromInteger #-} + + + +instance Fractional DoubleX2 where + DoubleX2# a1 a2 / DoubleX2# b1 b2 = DoubleX2# ((/##) a1 b1) + ((/##) a2 b2) + {-# INLINE (/) #-} + recip (DoubleX2# a1 a2) = DoubleX2# ((/##) 1.0## a1) + ((/##) 1.0## a2) + {-# INLINE recip #-} + fromRational r = case fromRational r of D# x -> DoubleX2# x x + {-# INLINE fromRational #-} + + + +instance Floating DoubleX2 where + pi = DoubleX2# 3.141592653589793238## 3.141592653589793238## + {-# INLINE pi #-} + exp (DoubleX2# a1 a2) = DoubleX2# (expDouble# a1) + (expDouble# a2) + {-# INLINE exp #-} + log (DoubleX2# a1 a2) = DoubleX2# (logDouble# a1) + (logDouble# a2) + {-# INLINE log #-} + sqrt (DoubleX2# a1 a2) = DoubleX2# (sqrtDouble# a1) + (sqrtDouble# a2) + {-# INLINE sqrt #-} + sin (DoubleX2# a1 a2) = DoubleX2# (sinDouble# a1) + (sinDouble# a2) + {-# INLINE sin #-} + cos (DoubleX2# a1 a2) = DoubleX2# (cosDouble# a1) + (cosDouble# a2) + {-# INLINE cos #-} + tan (DoubleX2# a1 a2) = DoubleX2# (tanDouble# a1) + (tanDouble# a2) + {-# INLINE tan #-} + asin (DoubleX2# a1 a2) = DoubleX2# (asinDouble# a1) + (asinDouble# a2) + {-# INLINE asin #-} + acos (DoubleX2# a1 a2) = DoubleX2# (acosDouble# a1) + (acosDouble# a2) + {-# INLINE acos #-} + atan (DoubleX2# a1 a2) = DoubleX2# (atanDouble# a1) + (atanDouble# a2) + {-# INLINE atan #-} + sinh (DoubleX2# a1 a2) = DoubleX2# (sinDouble# a1) + (sinDouble# a2) + {-# INLINE sinh #-} + cosh (DoubleX2# a1 a2) = DoubleX2# (coshDouble# a1) + (coshDouble# a2) + {-# INLINE cosh #-} + tanh (DoubleX2# a1 a2) = DoubleX2# (tanhDouble# a1) + (tanhDouble# a2) + {-# INLINE tanh #-} + DoubleX2# a1 a2 ** DoubleX2# b1 b2 = DoubleX2# ((**##) a1 b1) + ((**##) a2 b2) + {-# INLINE (**) #-} + + logBase x y = log y / log x + {-# INLINE logBase #-} + asinh x = log (x + sqrt (1.0+x*x)) + {-# INLINE asinh #-} + acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0))) + {-# INLINE acosh #-} + atanh x = 0.5 * log ((1.0+x) / (1.0-x)) + {-# INLINE atanh #-} + + + +type instance ElemRep DoubleX2 = 'DoubleRep +type instance ElemPrim DoubleX2 = Double# +instance PrimBytes DoubleX2 where + toBytes (DoubleX2# a1 a2) = case runRW# + ( \s0 -> case newByteArray# (SIZEOF_HSFLOAT# *# 2#) s0 of + (# s1, marr #) -> case writeDoubleArray# marr 0# a1 s1 of + s2 -> case writeDoubleArray# marr 1# a2 s2 of + s3 -> unsafeFreezeByteArray# marr s3 + ) of (# _, a #) -> (# 0#, 2#, a #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = DoubleX2# + (indexDoubleArray# arr off) + (indexDoubleArray# arr (off +# 1#)) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_HSFLOAT# *# 2# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_HSFLOAT# + {-# INLINE byteAlign #-} + elementByteSize _ = SIZEOF_HSFLOAT# + {-# INLINE elementByteSize #-} + ix 0# (DoubleX2# a1 _) = a1 + ix 1# (DoubleX2# _ a2) = a2 + ix _ _ = undefined + {-# INLINE ix #-} + + +instance ElementWise (Idx '[2]) Double DoubleX2 where + indexOffset# (DoubleX2# a1 _) 0# = D# a1 + indexOffset# (DoubleX2# _ a2) 1# = D# a2 + indexOffset# _ _ = undefined + {-# INLINE indexOffset# #-} + + (!) (DoubleX2# a1 _) ( 1 :! Z) = D# a1 + (!) (DoubleX2# _ a2) ( 2 :! Z) = D# a2 + (!) _ ( _ :! Z) = undefined + {-# INLINE (!) #-} + + broadcast (D# x) = DoubleX2# x x + {-# INLINE broadcast #-} + + ewmap f (DoubleX2# x y) = case (f (1:!Z) (D# x), f (2:!Z) (D# y)) of + (D# r1, D# r2) -> DoubleX2# r1 r2 + {-# INLINE ewmap #-} + + ewgen f = case (f (1:!Z), f (2:!Z)) of (D# r1, D# r2) -> DoubleX2# r1 r2 + {-# INLINE ewgen #-} + + ewgenA f = (\(D# r1) (D# r2) -> DoubleX2# r1 r2) <$> f (1:!Z) <*> f (2:!Z) + {-# INLINE ewgenA #-} + + ewfoldl f x0 (DoubleX2# x y) = f (2:!Z) (f (1:!Z) x0 (D# x)) (D# y) + {-# INLINE ewfoldl #-} + + ewfoldr f x0 (DoubleX2# x y) = f (1:!Z) (D# x) (f (2:!Z) (D# y) x0) + {-# INLINE ewfoldr #-} + + elementWise f (DoubleX2# x y) = (\(D# a) (D# b) -> DoubleX2# a b) + <$> f (D# x) <*> f (D# y) + {-# INLINE elementWise #-} + + indexWise f (DoubleX2# x y) = (\(D# a) (D# b) -> DoubleX2# a b) + <$> f (1:!Z) (D# x) <*> f (2:!Z) (D# y) + {-# INLINE indexWise #-} + + update (1 :! Z) (D# q) (DoubleX2# _ y) = DoubleX2# q y + update (2 :! Z) (D# q) (DoubleX2# x _) = DoubleX2# x q + update (_ :! Z) _ x = x + {-# INLINE update #-} diff --git a/src-base/Numeric/Array/Family/DoubleX3.hs b/src-base/Numeric/Array/Family/DoubleX3.hs new file mode 100644 --- /dev/null +++ b/src-base/Numeric/Array/Family/DoubleX3.hs @@ -0,0 +1,299 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UnboxedTuples #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Array.Family.DoubleX3 +-- Copyright : (c) Artem Chirkin +-- License : BSD3 +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Array.Family.DoubleX3 () where + + +#include "MachDeps.h" + +import GHC.Base (runRW#) +import GHC.Prim +import GHC.Types (Double (..), RuntimeRep (..), + isTrue#) + +import Numeric.Array.ElementWise +import Numeric.Array.Family +import Numeric.Commons +import Numeric.Dimensions + + + + + +instance Show DoubleX3 where + show (DoubleX3# a1 a2 a3) = "{ " ++ show (D# a1) + ++ ", " ++ show (D# a2) + ++ ", " ++ show (D# a3) + ++ " }" + + + +instance Eq DoubleX3 where + DoubleX3# a1 a2 a3 == DoubleX3# b1 b2 b3 = isTrue# ( (a1 ==## b1) + `andI#` (a2 ==## b2) + `andI#` (a3 ==## b3) + ) + {-# INLINE (==) #-} + DoubleX3# a1 a2 a3 /= DoubleX3# b1 b2 b3 = isTrue# ( (a1 /=## b1) + `orI#` (a2 /=## b2) + `orI#` (a3 /=## b3) + ) + {-# INLINE (/=) #-} + + + +-- | Implement partial ordering for `>`, `<`, `>=`, `<=` +-- and lexicographical ordering for `compare` +instance Ord DoubleX3 where + DoubleX3# a1 a2 a3 > DoubleX3# b1 b2 b3 = isTrue# ( (a1 >## b1) + `andI#` (a2 >## b2) + `andI#` (a3 >## b3) + ) + {-# INLINE (>) #-} + DoubleX3# a1 a2 a3 < DoubleX3# b1 b2 b3 = isTrue# ( (a1 <## b1) + `andI#` (a2 <## b2) + `andI#` (a3 <## b3) + ) + {-# INLINE (<) #-} + DoubleX3# a1 a2 a3 >= DoubleX3# b1 b2 b3 = isTrue# ( (a1 >=## b1) + `andI#` (a2 >=## b2) + `andI#` (a3 >=## b3) + ) + {-# INLINE (>=) #-} + DoubleX3# a1 a2 a3 <= DoubleX3# b1 b2 b3 = isTrue# ( (a1 <=## b1) + `andI#` (a2 <=## b2) + `andI#` (a3 <=## b3) + ) + {-# INLINE (<=) #-} + -- | Compare lexicographically + compare (DoubleX3# a1 a2 a3) (DoubleX3# b1 b2 b3) + | isTrue# (a1 >## b1) = GT + | isTrue# (a1 <## b1) = LT + | isTrue# (a2 >## b2) = GT + | isTrue# (a2 <## b2) = LT + | isTrue# (a3 >## b3) = GT + | isTrue# (a3 <## b3) = LT + | otherwise = EQ + {-# INLINE compare #-} + -- | Element-wise minimum + min (DoubleX3# a1 a2 a3) (DoubleX3# b1 b2 b3) = + DoubleX3# (if isTrue# (a1 >## b1) then b1 else a1) + (if isTrue# (a2 >## b2) then b2 else a2) + (if isTrue# (a3 >## b3) then b3 else a3) + {-# INLINE min #-} + -- | Element-wise maximum + max (DoubleX3# a1 a2 a3) (DoubleX3# b1 b2 b3) = + DoubleX3# (if isTrue# (a1 >## b1) then a1 else b1) + (if isTrue# (a2 >## b2) then a2 else b2) + (if isTrue# (a3 >## b3) then a3 else b3) + {-# INLINE max #-} + + + +-- | element-wise operations for vectors +instance Num DoubleX3 where + DoubleX3# a1 a2 a3 + DoubleX3# b1 b2 b3 + = DoubleX3# ((+##) a1 b1) ((+##) a2 b2) ((+##) a3 b3) + {-# INLINE (+) #-} + DoubleX3# a1 a2 a3 - DoubleX3# b1 b2 b3 + = DoubleX3# ((-##) a1 b1) ((-##) a2 b2) ((-##) a3 b3) + {-# INLINE (-) #-} + DoubleX3# a1 a2 a3 * DoubleX3# b1 b2 b3 + = DoubleX3# ((*##) a1 b1) ((*##) a2 b2) ((*##) a3 b3) + {-# INLINE (*) #-} + negate (DoubleX3# a1 a2 a3) + = DoubleX3# (negateDouble# a1) (negateDouble# a2) (negateDouble# a3) + {-# INLINE negate #-} + abs (DoubleX3# a1 a2 a3) + = DoubleX3# (if isTrue# (a1 >=## 0.0##) then a1 else negateDouble# a1) + (if isTrue# (a2 >=## 0.0##) then a2 else negateDouble# a2) + (if isTrue# (a3 >=## 0.0##) then a3 else negateDouble# a3) + {-# INLINE abs #-} + signum (DoubleX3# a1 a2 a3) + = DoubleX3# (if isTrue# (a1 >## 0.0##) + then 1.0## + else if isTrue# (a1 <## 0.0##) then -1.0## else 0.0## ) + (if isTrue# (a2 >## 0.0##) + then 1.0## + else if isTrue# (a2 <## 0.0##) then -1.0## else 0.0## ) + (if isTrue# (a3 >## 0.0##) + then 1.0## + else if isTrue# (a3 <## 0.0##) then -1.0## else 0.0## ) + {-# INLINE signum #-} + fromInteger n = case fromInteger n of D# x -> DoubleX3# x x x + {-# INLINE fromInteger #-} + + + +instance Fractional DoubleX3 where + DoubleX3# a1 a2 a3 / DoubleX3# b1 b2 b3 = DoubleX3# ((/##) a1 b1) + ((/##) a2 b2) + ((/##) a3 b3) + {-# INLINE (/) #-} + recip (DoubleX3# a1 a2 a3) = DoubleX3# ((/##) 1.0## a1) + ((/##) 1.0## a2) + ((/##) 1.0## a3) + {-# INLINE recip #-} + fromRational r = case fromRational r of D# x -> DoubleX3# x x x + {-# INLINE fromRational #-} + + + +instance Floating DoubleX3 where + pi = DoubleX3# 3.141592653589793238## 3.141592653589793238## 3.141592653589793238## + {-# INLINE pi #-} + exp (DoubleX3# a1 a2 a3) = DoubleX3# (expDouble# a1) + (expDouble# a2) + (expDouble# a3) + {-# INLINE exp #-} + log (DoubleX3# a1 a2 a3) = DoubleX3# (logDouble# a1) + (logDouble# a2) + (logDouble# a3) + {-# INLINE log #-} + sqrt (DoubleX3# a1 a2 a3) = DoubleX3# (sqrtDouble# a1) + (sqrtDouble# a2) + (sqrtDouble# a3) + {-# INLINE sqrt #-} + sin (DoubleX3# a1 a2 a3) = DoubleX3# (sinDouble# a1) + (sinDouble# a2) + (sinDouble# a3) + {-# INLINE sin #-} + cos (DoubleX3# a1 a2 a3) = DoubleX3# (cosDouble# a1) + (cosDouble# a2) + (cosDouble# a3) + {-# INLINE cos #-} + tan (DoubleX3# a1 a2 a3) = DoubleX3# (tanDouble# a1) + (tanDouble# a2) + (tanDouble# a3) + {-# INLINE tan #-} + asin (DoubleX3# a1 a2 a3) = DoubleX3# (asinDouble# a1) + (asinDouble# a2) + (asinDouble# a3) + {-# INLINE asin #-} + acos (DoubleX3# a1 a2 a3) = DoubleX3# (acosDouble# a1) + (acosDouble# a2) + (acosDouble# a3) + {-# INLINE acos #-} + atan (DoubleX3# a1 a2 a3) = DoubleX3# (atanDouble# a1) + (atanDouble# a2) + (atanDouble# a3) + {-# INLINE atan #-} + sinh (DoubleX3# a1 a2 a3) = DoubleX3# (sinDouble# a1) + (sinDouble# a2) + (sinDouble# a3) + {-# INLINE sinh #-} + cosh (DoubleX3# a1 a2 a3) = DoubleX3# (coshDouble# a1) + (coshDouble# a2) + (coshDouble# a3) + {-# INLINE cosh #-} + tanh (DoubleX3# a1 a2 a3) = DoubleX3# (tanhDouble# a1) + (tanhDouble# a2) + (tanhDouble# a3) + {-# INLINE tanh #-} + DoubleX3# a1 a2 a3 ** DoubleX3# b1 b2 b3 = DoubleX3# ((**##) a1 b1) + ((**##) a2 b2) + ((**##) a3 b3) + {-# INLINE (**) #-} + + logBase x y = log y / log x + {-# INLINE logBase #-} + asinh x = log (x + sqrt (1.0+x*x)) + {-# INLINE asinh #-} + acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0))) + {-# INLINE acosh #-} + atanh x = 0.5 * log ((1.0+x) / (1.0-x)) + {-# INLINE atanh #-} + + + +type instance ElemRep DoubleX3 = 'DoubleRep +type instance ElemPrim DoubleX3 = Double# +instance PrimBytes DoubleX3 where + toBytes (DoubleX3# a1 a2 a3) = case runRW# + ( \s0 -> case newByteArray# (SIZEOF_HSFLOAT# *# 3#) s0 of + (# s1, marr #) -> case writeDoubleArray# marr 0# a1 s1 of + s2 -> case writeDoubleArray# marr 1# a2 s2 of + s3 -> case writeDoubleArray# marr 2# a3 s3 of + s4 -> unsafeFreezeByteArray# marr s4 + ) of (# _, a #) -> (# 0#, 3#, a #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = DoubleX3# + (indexDoubleArray# arr off) + (indexDoubleArray# arr (off +# 1#)) + (indexDoubleArray# arr (off +# 2#)) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_HSFLOAT# *# 3# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_HSFLOAT# + {-# INLINE byteAlign #-} + elementByteSize _ = SIZEOF_HSFLOAT# + {-# INLINE elementByteSize #-} + ix 0# (DoubleX3# a1 _ _) = a1 + ix 1# (DoubleX3# _ a2 _) = a2 + ix 2# (DoubleX3# _ _ a3) = a3 + ix _ _ = undefined + {-# INLINE ix #-} + + +instance ElementWise (Idx '[3]) Double DoubleX3 where + indexOffset# (DoubleX3# a1 _ _) 0# = D# a1 + indexOffset# (DoubleX3# _ a2 _) 1# = D# a2 + indexOffset# (DoubleX3# _ _ a3) 2# = D# a3 + indexOffset# _ _ = undefined + {-# INLINE indexOffset# #-} + + (!) (DoubleX3# a1 _ _) ( 1 :! Z) = D# a1 + (!) (DoubleX3# _ a2 _) ( 2 :! Z) = D# a2 + (!) (DoubleX3# _ _ a3) ( 3 :! Z) = D# a3 + (!) _ ( _ :! Z) = undefined + {-# INLINE (!) #-} + + broadcast (D# x) = DoubleX3# x x x + {-# INLINE broadcast #-} + + ewmap f (DoubleX3# x y z) = case (f (1:!Z) (D# x), f (2:!Z) (D# y), f (3:!Z) (D# z)) of + (D# r1, D# r2, D# r3) -> DoubleX3# r1 r2 r3 + {-# INLINE ewmap #-} + + ewgen f = case (f (1:!Z), f (2:!Z), f (3:!Z)) of (D# r1, D# r2, D# r3) -> DoubleX3# r1 r2 r3 + {-# INLINE ewgen #-} + + ewgenA f = (\(D# r1) (D# r2) (D# r3) -> DoubleX3# r1 r2 r3) + <$> f (1:!Z) <*> f (2:!Z) <*> f (3:!Z) + {-# INLINE ewgenA #-} + + ewfoldl f x0 (DoubleX3# x y z) = f (3:!Z) (f (2:!Z) (f (1:!Z) x0 (D# x)) (D# y)) (D# z) + {-# INLINE ewfoldl #-} + + ewfoldr f x0 (DoubleX3# x y z) = f (1:!Z) (D# x) (f (2:!Z) (D# y) (f (3:!Z) (D# z) x0)) + {-# INLINE ewfoldr #-} + + elementWise f (DoubleX3# x y z) = (\(D# a) (D# b) (D# c) -> DoubleX3# a b c) + <$> f (D# x) <*> f (D# y) <*> f (D# z) + {-# INLINE elementWise #-} + + indexWise f (DoubleX3# x y z) = (\(D# a) (D# b) (D# c) -> DoubleX3# a b c) + <$> f (1:!Z) (D# x) <*> f (2:!Z) (D# y) <*> f (3:!Z) (D# z) + {-# INLINE indexWise #-} + + update (1 :! Z) (D# q) (DoubleX3# _ y z) = DoubleX3# q y z + update (2 :! Z) (D# q) (DoubleX3# x _ z) = DoubleX3# x q z + update (3 :! Z) (D# q) (DoubleX3# x y _) = DoubleX3# x y q + update (_ :! Z) _ x = x + {-# INLINE update #-} diff --git a/src-base/Numeric/Array/Family/DoubleX4.hs b/src-base/Numeric/Array/Family/DoubleX4.hs new file mode 100644 --- /dev/null +++ b/src-base/Numeric/Array/Family/DoubleX4.hs @@ -0,0 +1,335 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UnboxedTuples #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Array.Family.DoubleX4 +-- Copyright : (c) Artem Chirkin +-- License : BSD3 +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Array.Family.DoubleX4 () where + + +#include "MachDeps.h" + +import GHC.Base (runRW#) +import GHC.Prim +import GHC.Types (Double (..), RuntimeRep (..), + isTrue#) + +import Numeric.Array.ElementWise +import Numeric.Array.Family +import Numeric.Commons +import Numeric.Dimensions + + + + + +instance Show DoubleX4 where + show (DoubleX4# a1 a2 a3 a4) = "{ " ++ show (D# a1) + ++ ", " ++ show (D# a2) + ++ ", " ++ show (D# a3) + ++ ", " ++ show (D# a4) + ++ " }" + + + +instance Eq DoubleX4 where + DoubleX4# a1 a2 a3 a4 == DoubleX4# b1 b2 b3 b4 = isTrue# ( (a1 ==## b1) + `andI#` (a2 ==## b2) + `andI#` (a3 ==## b3) + `andI#` (a4 ==## b4) + ) + {-# INLINE (==) #-} + DoubleX4# a1 a2 a3 a4 /= DoubleX4# b1 b2 b3 b4 = isTrue# ( (a1 /=## b1) + `orI#` (a2 /=## b2) + `orI#` (a3 /=## b3) + `orI#` (a4 /=## b4) + ) + {-# INLINE (/=) #-} + + + +-- | Implement partial ordering for `>`, `<`, `>=`, `<=` +-- and lexicographical ordering for `compare` +instance Ord DoubleX4 where + DoubleX4# a1 a2 a3 a4 > DoubleX4# b1 b2 b3 b4 = isTrue# ( (a1 >## b1) + `andI#` (a2 >## b2) + `andI#` (a3 >## b3) + `andI#` (a4 >## b4) + ) + {-# INLINE (>) #-} + DoubleX4# a1 a2 a3 a4 < DoubleX4# b1 b2 b3 b4 = isTrue# ( (a1 <## b1) + `andI#` (a2 <## b2) + `andI#` (a3 <## b3) + `andI#` (a4 <## b4) + ) + {-# INLINE (<) #-} + DoubleX4# a1 a2 a3 a4 >= DoubleX4# b1 b2 b3 b4 = isTrue# ( (a1 >=## b1) + `andI#` (a2 >=## b2) + `andI#` (a3 >=## b3) + `andI#` (a4 >=## b4) + ) + {-# INLINE (>=) #-} + DoubleX4# a1 a2 a3 a4 <= DoubleX4# b1 b2 b3 b4 = isTrue# ( (a1 <=## b1) + `andI#` (a2 <=## b2) + `andI#` (a3 <=## b3) + `andI#` (a4 <=## b4) + ) + {-# INLINE (<=) #-} + -- | Compare lexicographically + compare (DoubleX4# a1 a2 a3 a4) (DoubleX4# b1 b2 b3 b4) + | isTrue# (a1 >## b1) = GT + | isTrue# (a1 <## b1) = LT + | isTrue# (a2 >## b2) = GT + | isTrue# (a2 <## b2) = LT + | isTrue# (a3 >## b3) = GT + | isTrue# (a3 <## b3) = LT + | isTrue# (a4 >## b4) = GT + | isTrue# (a4 <## b4) = LT + | otherwise = EQ + {-# INLINE compare #-} + -- | Element-wise minimum + min (DoubleX4# a1 a2 a3 a4) (DoubleX4# b1 b2 b3 b4) = + DoubleX4# (if isTrue# (a1 >## b1) then b1 else a1) + (if isTrue# (a2 >## b2) then b2 else a2) + (if isTrue# (a3 >## b3) then b3 else a3) + (if isTrue# (a4 >## b4) then b4 else a4) + {-# INLINE min #-} + -- | Element-wise maximum + max (DoubleX4# a1 a2 a3 a4) (DoubleX4# b1 b2 b3 b4) = + DoubleX4# (if isTrue# (a1 >## b1) then a1 else b1) + (if isTrue# (a2 >## b2) then a2 else b2) + (if isTrue# (a3 >## b3) then a3 else b3) + (if isTrue# (a4 >## b4) then a4 else b4) + {-# INLINE max #-} + + + +-- | element-wise operations for vectors +instance Num DoubleX4 where + DoubleX4# a1 a2 a3 a4 + DoubleX4# b1 b2 b3 b4 + = DoubleX4# ((+##) a1 b1) ((+##) a2 b2) ((+##) a3 b3) ((+##) a4 b4) + {-# INLINE (+) #-} + DoubleX4# a1 a2 a3 a4 - DoubleX4# b1 b2 b3 b4 + = DoubleX4# ((-##) a1 b1) ((-##) a2 b2) ((-##) a3 b3) ((-##) a4 b4) + {-# INLINE (-) #-} + DoubleX4# a1 a2 a3 a4 * DoubleX4# b1 b2 b3 b4 + = DoubleX4# ((*##) a1 b1) ((*##) a2 b2) ((*##) a3 b3) ((*##) a4 b4) + {-# INLINE (*) #-} + negate (DoubleX4# a1 a2 a3 a4) + = DoubleX4# (negateDouble# a1) (negateDouble# a2) (negateDouble# a3) (negateDouble# a4) + {-# INLINE negate #-} + abs (DoubleX4# a1 a2 a3 a4) + = DoubleX4# (if isTrue# (a1 >=## 0.0##) then a1 else negateDouble# a1) + (if isTrue# (a2 >=## 0.0##) then a2 else negateDouble# a2) + (if isTrue# (a3 >=## 0.0##) then a3 else negateDouble# a3) + (if isTrue# (a4 >=## 0.0##) then a4 else negateDouble# a4) + {-# INLINE abs #-} + signum (DoubleX4# a1 a2 a3 a4) + = DoubleX4# (if isTrue# (a1 >## 0.0##) + then 1.0## + else if isTrue# (a1 <## 0.0##) then -1.0## else 0.0## ) + (if isTrue# (a2 >## 0.0##) + then 1.0## + else if isTrue# (a2 <## 0.0##) then -1.0## else 0.0## ) + (if isTrue# (a3 >## 0.0##) + then 1.0## + else if isTrue# (a3 <## 0.0##) then -1.0## else 0.0## ) + (if isTrue# (a4 >## 0.0##) + then 1.0## + else if isTrue# (a4 <## 0.0##) then -1.0## else 0.0## ) + {-# INLINE signum #-} + fromInteger n = case fromInteger n of D# x -> DoubleX4# x x x x + {-# INLINE fromInteger #-} + + + +instance Fractional DoubleX4 where + DoubleX4# a1 a2 a3 a4 / DoubleX4# b1 b2 b3 b4 = DoubleX4# ((/##) a1 b1) + ((/##) a2 b2) + ((/##) a3 b3) + ((/##) a4 b4) + {-# INLINE (/) #-} + recip (DoubleX4# a1 a2 a3 a4) = DoubleX4# ((/##) 1.0## a1) + ((/##) 1.0## a2) + ((/##) 1.0## a3) + ((/##) 1.0## a4) + {-# INLINE recip #-} + fromRational r = case fromRational r of D# x -> DoubleX4# x x x x + {-# INLINE fromRational #-} + + + +instance Floating DoubleX4 where + pi = DoubleX4# 3.141592653589793238## 3.141592653589793238## 3.141592653589793238## 3.141592653589793238## + {-# INLINE pi #-} + exp (DoubleX4# a1 a2 a3 a4) = DoubleX4# (expDouble# a1) + (expDouble# a2) + (expDouble# a3) + (expDouble# a4) + {-# INLINE exp #-} + log (DoubleX4# a1 a2 a3 a4) = DoubleX4# (logDouble# a1) + (logDouble# a2) + (logDouble# a3) + (logDouble# a4) + {-# INLINE log #-} + sqrt (DoubleX4# a1 a2 a3 a4) = DoubleX4# (sqrtDouble# a1) + (sqrtDouble# a2) + (sqrtDouble# a3) + (sqrtDouble# a4) + {-# INLINE sqrt #-} + sin (DoubleX4# a1 a2 a3 a4) = DoubleX4# (sinDouble# a1) + (sinDouble# a2) + (sinDouble# a3) + (sinDouble# a4) + {-# INLINE sin #-} + cos (DoubleX4# a1 a2 a3 a4) = DoubleX4# (cosDouble# a1) + (cosDouble# a2) + (cosDouble# a3) + (cosDouble# a4) + {-# INLINE cos #-} + tan (DoubleX4# a1 a2 a3 a4) = DoubleX4# (tanDouble# a1) + (tanDouble# a2) + (tanDouble# a3) + (tanDouble# a4) + {-# INLINE tan #-} + asin (DoubleX4# a1 a2 a3 a4) = DoubleX4# (asinDouble# a1) + (asinDouble# a2) + (asinDouble# a3) + (asinDouble# a4) + {-# INLINE asin #-} + acos (DoubleX4# a1 a2 a3 a4) = DoubleX4# (acosDouble# a1) + (acosDouble# a2) + (acosDouble# a3) + (acosDouble# a4) + {-# INLINE acos #-} + atan (DoubleX4# a1 a2 a3 a4) = DoubleX4# (atanDouble# a1) + (atanDouble# a2) + (atanDouble# a3) + (atanDouble# a4) + {-# INLINE atan #-} + sinh (DoubleX4# a1 a2 a3 a4) = DoubleX4# (sinDouble# a1) + (sinDouble# a2) + (sinDouble# a3) + (sinDouble# a4) + {-# INLINE sinh #-} + cosh (DoubleX4# a1 a2 a3 a4) = DoubleX4# (coshDouble# a1) + (coshDouble# a2) + (coshDouble# a3) + (coshDouble# a4) + {-# INLINE cosh #-} + tanh (DoubleX4# a1 a2 a3 a4) = DoubleX4# (tanhDouble# a1) + (tanhDouble# a2) + (tanhDouble# a3) + (tanhDouble# a4) + {-# INLINE tanh #-} + DoubleX4# a1 a2 a3 a4 ** DoubleX4# b1 b2 b3 b4 = DoubleX4# ((**##) a1 b1) + ((**##) a2 b2) + ((**##) a3 b3) + ((**##) a4 b4) + {-# INLINE (**) #-} + + logBase x y = log y / log x + {-# INLINE logBase #-} + asinh x = log (x + sqrt (1.0+x*x)) + {-# INLINE asinh #-} + acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0))) + {-# INLINE acosh #-} + atanh x = 0.5 * log ((1.0+x) / (1.0-x)) + {-# INLINE atanh #-} + + + +type instance ElemRep DoubleX4 = 'DoubleRep +type instance ElemPrim DoubleX4 = Double# +instance PrimBytes DoubleX4 where + toBytes (DoubleX4# a1 a2 a3 a4) = case runRW# + ( \s0 -> case newByteArray# (SIZEOF_HSFLOAT# *# 3#) s0 of + (# s1, marr #) -> case writeDoubleArray# marr 0# a1 s1 of + s2 -> case writeDoubleArray# marr 1# a2 s2 of + s3 -> case writeDoubleArray# marr 2# a3 s3 of + s4 -> case writeDoubleArray# marr 3# a4 s4 of + s5 -> unsafeFreezeByteArray# marr s5 + ) of (# _, a #) -> (# 0#, 4#, a #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = DoubleX4# + (indexDoubleArray# arr off) + (indexDoubleArray# arr (off +# 1#)) + (indexDoubleArray# arr (off +# 2#)) + (indexDoubleArray# arr (off +# 3#)) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_HSFLOAT# *# 4# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_HSFLOAT# + {-# INLINE byteAlign #-} + elementByteSize _ = SIZEOF_HSFLOAT# + {-# INLINE elementByteSize #-} + ix 0# (DoubleX4# a1 _ _ _) = a1 + ix 1# (DoubleX4# _ a2 _ _) = a2 + ix 2# (DoubleX4# _ _ a3 _) = a3 + ix 3# (DoubleX4# _ _ _ a4) = a4 + ix _ _ = undefined + {-# INLINE ix #-} + + +instance ElementWise (Idx '[4]) Double DoubleX4 where + indexOffset# (DoubleX4# a1 _ _ _) 0# = D# a1 + indexOffset# (DoubleX4# _ a2 _ _) 1# = D# a2 + indexOffset# (DoubleX4# _ _ a3 _) 2# = D# a3 + indexOffset# (DoubleX4# _ _ _ a4) 3# = D# a4 + indexOffset# _ _ = undefined + {-# INLINE indexOffset# #-} + + (!) (DoubleX4# a1 _ _ _) ( 1 :! Z) = D# a1 + (!) (DoubleX4# _ a2 _ _) ( 2 :! Z) = D# a2 + (!) (DoubleX4# _ _ a3 _) ( 3 :! Z) = D# a3 + (!) (DoubleX4# _ _ _ a4) ( 4 :! Z) = D# a4 + (!) _ ( _ :! Z) = undefined + {-# INLINE (!) #-} + + broadcast (D# x) = DoubleX4# x x x x + {-# INLINE broadcast #-} + + ewmap f (DoubleX4# x y z w) = case (f (1:!Z) (D# x), f (2:!Z) (D# y), f (3:!Z) (D# z), f (3:!Z) (D# w)) of + (D# r1, D# r2, D# r3, D# r4) -> DoubleX4# r1 r2 r3 r4 + {-# INLINE ewmap #-} + + ewgen f = case (f (1:!Z), f (2:!Z), f (3:!Z), f (4:!Z)) of (D# r1, D# r2, D# r3, D# r4) -> DoubleX4# r1 r2 r3 r4 + {-# INLINE ewgen #-} + + ewgenA f = (\(D# a) (D# b) (D# c) (D# d) -> DoubleX4# a b c d) + <$> f (1:!Z) <*> f (2:!Z) <*> f (3:!Z) <*> f (4:!Z) + {-# INLINE ewgenA #-} + + ewfoldl f x0 (DoubleX4# x y z w) = f (4:!Z) (f (3:!Z) (f (2:!Z) (f (1:!Z) x0 (D# x)) (D# y)) (D# z)) (D# w) + {-# INLINE ewfoldl #-} + + ewfoldr f x0 (DoubleX4# x y z w) = f (1:!Z) (D# x) (f (2:!Z) (D# y) (f (3:!Z) (D# z) (f (4:!Z) (D# w) x0))) + {-# INLINE ewfoldr #-} + + elementWise f (DoubleX4# x y z w) = (\(D# a) (D# b) (D# c) (D# d) -> DoubleX4# a b c d) + <$> f (D# x) <*> f (D# y) <*> f (D# z) <*> f (D# w) + {-# INLINE elementWise #-} + + indexWise f (DoubleX4# x y z w) = (\(D# a) (D# b) (D# c) (D# d) -> DoubleX4# a b c d) + <$> f (1:!Z) (D# x) <*> f (2:!Z) (D# y) <*> f (3:!Z) (D# z) <*> f (4:!Z) (D# w) + {-# INLINE indexWise #-} + + update (1 :! Z) (D# q) (DoubleX4# _ y z w) = DoubleX4# q y z w + update (2 :! Z) (D# q) (DoubleX4# x _ z w) = DoubleX4# x q z w + update (3 :! Z) (D# q) (DoubleX4# x y _ w) = DoubleX4# x y q w + update (4 :! Z) (D# q) (DoubleX4# x y z _) = DoubleX4# x y z q + update (_ :! Z) _ x = x + {-# INLINE update #-} diff --git a/src-base/Numeric/DataFrame/Contraction.hs b/src-base/Numeric/DataFrame/Contraction.hs --- a/src-base/Numeric/DataFrame/Contraction.hs +++ b/src-base/Numeric/DataFrame/Contraction.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} @@ -30,6 +31,7 @@ ( Contraction (..), (%*) ) where +#include "MachDeps.h" import Data.Int (Int16, Int32, Int64, Int8) import Data.Word (Word16, Word32, Word64, Word8) @@ -198,10 +200,17 @@ , I# m <- intNatVal pm , I# n <- totalDim (Proxy @as) , I# k <- totalDim (Proxy @bs) +#if WORD_SIZE_IN_BITS < 64 + , Refl <- unsafeCoerce Refl :: ElemRep (Array Int64 (m : bs) ) :~: 'Int64Rep + , Refl <- unsafeCoerce Refl :: ElemPrim (Array Int64 (m : bs) ) :~: Int64# + , Refl <- unsafeCoerce Refl :: ElemRep (Array Int64 (as +: m)) :~: 'Int64Rep + , Refl <- unsafeCoerce Refl :: ElemPrim (Array Int64 (as +: m)) :~: Int64# +#else , Refl <- unsafeCoerce Refl :: ElemRep (Array Int64 (m : bs) ) :~: 'IntRep , Refl <- unsafeCoerce Refl :: ElemPrim (Array Int64 (m : bs) ) :~: Int# , Refl <- unsafeCoerce Refl :: ElemRep (Array Int64 (as +: m)) :~: 'IntRep , Refl <- unsafeCoerce Refl :: ElemPrim (Array Int64 (as +: m)) :~: Int# +#endif = prodI64 n m k x y where getM :: forall m p . p (m ': bs) -> Proxy m @@ -290,18 +299,27 @@ , I# m <- intNatVal pm , I# n <- totalDim (Proxy @as) , I# k <- totalDim (Proxy @bs) +#if WORD_SIZE_IN_BITS < 64 + , Refl <- unsafeCoerce Refl :: ElemRep (Array Word64 (m : bs) ) :~: 'Word64Rep + , Refl <- unsafeCoerce Refl :: ElemPrim (Array Word64 (m : bs) ) :~: Word64# + , Refl <- unsafeCoerce Refl :: ElemRep (Array Word64 (as +: m)) :~: 'Word64Rep + , Refl <- unsafeCoerce Refl :: ElemPrim (Array Word64 (as +: m)) :~: Word64# +#else , Refl <- unsafeCoerce Refl :: ElemRep (Array Word64 (m : bs) ) :~: 'WordRep , Refl <- unsafeCoerce Refl :: ElemPrim (Array Word64 (m : bs) ) :~: Word# , Refl <- unsafeCoerce Refl :: ElemRep (Array Word64 (as +: m)) :~: 'WordRep , Refl <- unsafeCoerce Refl :: ElemPrim (Array Word64 (as +: m)) :~: Word# +#endif = prodW64 n m k x y where getM :: forall m p . p (m ': bs) -> Proxy m getM _ = Proxy - -prodF :: (FloatBytes a, FloatBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodF :: (PrimBytes a, PrimBytes b, PrimBytes c + , ElemPrim a ~ Float#, ElemRep a ~ 'FloatRep + , ElemPrim b ~ Float#, ElemRep b ~ 'FloatRep + ) => Int# -> Int# -> Int# -> a -> b -> c prodF n m k x y = case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -317,7 +335,10 @@ bs = n *# k *# elementByteSize x {-# INLINE prodF #-} -prodD :: (DoubleBytes a, DoubleBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodD :: (PrimBytes a, PrimBytes b, PrimBytes c + , ElemPrim a ~ Double#, ElemRep a ~ 'DoubleRep + , ElemPrim b ~ Double#, ElemRep b ~ 'DoubleRep + ) => Int# -> Int# -> Int# -> a -> b -> c prodD n m k x y= case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -333,7 +354,10 @@ bs = n *# k *# elementByteSize x {-# INLINE prodD #-} -prodI :: (IntBytes a, IntBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodI :: (PrimBytes a, PrimBytes b, PrimBytes c + , ElemPrim a ~ Int#, ElemRep a ~ 'IntRep + , ElemPrim b ~ Int#, ElemRep b ~ 'IntRep + ) => Int# -> Int# -> Int# -> a -> b -> c prodI n m k x y= case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -349,7 +373,10 @@ bs = n *# k *# elementByteSize x {-# INLINE prodI #-} -prodI8 :: (IntBytes a, IntBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodI8 :: (PrimBytes a, PrimBytes b, PrimBytes c + , ElemPrim a ~ Int#, ElemRep a ~ 'IntRep + , ElemPrim b ~ Int#, ElemRep b ~ 'IntRep + ) => Int# -> Int# -> Int# -> a -> b -> c prodI8 n m k x y= case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -366,7 +393,10 @@ {-# INLINE prodI8 #-} -prodI16 :: (IntBytes a, IntBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodI16 :: ( PrimBytes a, PrimBytes b, PrimBytes c + , ElemPrim a ~ Int#, ElemRep a ~ 'IntRep + , ElemPrim b ~ Int#, ElemRep b ~ 'IntRep + ) => Int# -> Int# -> Int# -> a -> b -> c prodI16 n m k x y= case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -383,7 +413,10 @@ {-# INLINE prodI16 #-} -prodI32 :: (IntBytes a, IntBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodI32 :: ( PrimBytes a, PrimBytes b, PrimBytes c + , ElemPrim a ~ Int#, ElemRep a ~ 'IntRep + , ElemPrim b ~ Int#, ElemRep b ~ 'IntRep + ) => Int# -> Int# -> Int# -> a -> b -> c prodI32 n m k x y= case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -399,7 +432,19 @@ bs = n *# k *# elementByteSize x {-# INLINE prodI32 #-} -prodI64 :: (IntBytes a, IntBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c + +prodI64 :: ( PrimBytes a, PrimBytes b, PrimBytes c +#if WORD_SIZE_IN_BITS < 64 + , ElemPrim a ~ Int64#, ElemRep a ~ 'Int64Rep + , ElemPrim b ~ Int64#, ElemRep b ~ 'Int64Rep +#else + , ElemPrim a ~ Int#, ElemRep a ~ 'IntRep + , ElemPrim b ~ Int#, ElemRep b ~ 'IntRep +#endif + ) => Int# -> Int# -> Int# -> a -> b -> c +#if WORD_SIZE_IN_BITS < 64 +prodI64 = undefined +#else prodI64 n m k x y= case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -414,9 +459,12 @@ where bs = n *# k *# elementByteSize x {-# INLINE prodI64 #-} - +#endif -prodW :: (WordBytes a, WordBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodW :: ( PrimBytes a, PrimBytes b, PrimBytes c + , ElemPrim a ~ Word#, ElemRep a ~ 'WordRep + , ElemPrim b ~ Word#, ElemRep b ~ 'WordRep + ) => Int# -> Int# -> Int# -> a -> b -> c prodW n m k x y = case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -432,7 +480,10 @@ bs = n *# k *# elementByteSize x {-# INLINE prodW #-} -prodW8 :: (WordBytes a, WordBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodW8 :: ( PrimBytes a, PrimBytes b, PrimBytes c + , ElemPrim a ~ Word#, ElemRep a ~ 'WordRep + , ElemPrim b ~ Word#, ElemRep b ~ 'WordRep + ) => Int# -> Int# -> Int# -> a -> b -> c prodW8 n m k x y = case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -449,7 +500,10 @@ {-# INLINE prodW8 #-} -prodW16 :: (WordBytes a, WordBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodW16 :: ( PrimBytes a, PrimBytes b, PrimBytes c + , ElemPrim a ~ Word#, ElemRep a ~ 'WordRep + , ElemPrim b ~ Word#, ElemRep b ~ 'WordRep + ) => Int# -> Int# -> Int# -> a -> b -> c prodW16 n m k x y = case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -465,7 +519,10 @@ bs = n *# k *# elementByteSize x {-# INLINE prodW16 #-} -prodW32 :: (WordBytes a, WordBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodW32 :: ( PrimBytes a, PrimBytes b, PrimBytes c + , ElemPrim a ~ Word#, ElemRep a ~ 'WordRep + , ElemPrim b ~ Word#, ElemRep b ~ 'WordRep + ) => Int# -> Int# -> Int# -> a -> b -> c prodW32 n m k x y = case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -481,7 +538,18 @@ bs = n *# k *# elementByteSize x {-# INLINE prodW32 #-} -prodW64 :: (WordBytes a, WordBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodW64 :: ( PrimBytes a, PrimBytes b, PrimBytes c +#if WORD_SIZE_IN_BITS < 64 + , ElemPrim a ~ Word64#, ElemRep a ~ 'Word64Rep + , ElemPrim b ~ Word64#, ElemRep b ~ 'Word64Rep +#else + , ElemPrim a ~ Word#, ElemRep a ~ 'WordRep + , ElemPrim b ~ Word#, ElemRep b ~ 'WordRep +#endif + ) => Int# -> Int# -> Int# -> a -> b -> c +#if WORD_SIZE_IN_BITS < 64 +prodW64 = undefined +#else prodW64 n m k x y = case runRW# ( \s0 -> case newByteArray# bs s0 of (# s1, marr #) -> @@ -496,6 +564,7 @@ where bs = n *# k *# elementByteSize x {-# INLINE prodW64 #-} +#endif -- | Do something in a loop for int i from 0 to n-1 and j from 0 to m-1 loop2# :: Int# -> Int# -> (Int# -> Int#-> State# s -> State# s) -> State# s -> State# s @@ -505,5 +574,3 @@ | isTrue# (i ==# n) = loop' 0# (j +# 1#) s | otherwise = case f i j s of s1 -> loop' (i +# 1#) j s1 {-# INLINE loop2# #-} - - diff --git a/src-base/Numeric/DataFrame/Inference.hs b/src-base/Numeric/DataFrame/Inference.hs --- a/src-base/Numeric/DataFrame/Inference.hs +++ b/src-base/Numeric/DataFrame/Inference.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE ExistentialQuantification #-} {-# LANGUAGE FlexibleContexts #-} @@ -61,7 +62,7 @@ ETWord8 -> Evidence ETWord16 -> Evidence ETWord32 -> Evidence - ETWord64 -> Evidence + ETWord64 -> Evidence AIArrayF -> Evidence AIArrayD -> Evidence AIArrayI -> Evidence @@ -77,6 +78,9 @@ AIFloatX2 -> Evidence AIFloatX3 -> Evidence AIFloatX4 -> Evidence + AIDoubleX2 -> Evidence + AIDoubleX3 -> Evidence + AIDoubleX4 -> Evidence inferElementWise :: forall t (ds :: [Nat]) . ( ArrayInstanceInference t ds @@ -100,6 +104,9 @@ AIFloatX2 -> Evidence AIFloatX3 -> Evidence AIFloatX4 -> Evidence + AIDoubleX2 -> Evidence + AIDoubleX3 -> Evidence + AIDoubleX4 -> Evidence inferNumericFrame :: forall t (ds :: [Nat]) @@ -113,6 +120,9 @@ AIFloatX2 -> Evidence AIFloatX3 -> Evidence AIFloatX4 -> Evidence + AIDoubleX2 -> Evidence + AIDoubleX3 -> Evidence + AIDoubleX4 -> Evidence AIScalar -> case elemTypeInstance @t of ETFloat -> Evidence ETDouble -> Evidence @@ -125,7 +135,7 @@ ETWord8 -> Evidence ETWord16 -> Evidence ETWord32 -> Evidence - ETWord64 -> Evidence + ETWord64 -> Evidence AIArrayF -> Evidence AIArrayD -> Evidence AIArrayI -> Evidence diff --git a/src-base/Numeric/DataFrame/Mutable.hs b/src-base/Numeric/DataFrame/Mutable.hs new file mode 100644 --- /dev/null +++ b/src-base/Numeric/DataFrame/Mutable.hs @@ -0,0 +1,224 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UnboxedTuples #-} +{-# LANGUAGE TypeOperators #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.DataFrame.Mutable +-- Copyright : (c) Artem Chirkin +-- License : BSD3 +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- Interfrace to perform primitive stateful operations on mutable frames. +-- +----------------------------------------------------------------------------- + +module Numeric.DataFrame.Mutable + ( MutableFrame (..), MDataFrame () + , newDataFrame#, copyDataFrame#, copyMDataFrame#, unsafeFreezeDataFrame# + , freezeDataFrame#, thawDataFrame# + , writeDataFrame#, readDataFrame# + ) where + + +import GHC.Int (Int16 (..), Int32 (..), Int64 (..), + Int8 (..)) +import GHC.Prim +import GHC.Types (Double (..), Float (..), Int (..), + Word (..)) +import GHC.Word (Word16 (..), Word32 (..), Word64 (..), + Word8 (..)) + +import Numeric.Commons +import Numeric.DataFrame.Type +import Numeric.Dimensions +import Numeric.TypeLits + + +-- | Mutable DataFrame type +data MDataFrame s t (ns :: [Nat]) = MDataFrame# Int# Int# (MutableByteArray# s) + +-- | Create a new mutable DataFrame. +newDataFrame# :: forall t (ns :: [Nat]) s + . ( PrimBytes t, Dimensions ns) + => State# s -> (# State# s, MDataFrame s t ns #) +newDataFrame# s0 + | elS <- elementByteSize (undefined :: t) + , I# n <- totalDim (Proxy @ns) + , (# s1, mba #) <- newByteArray# (n *# elS) s0 + = (# s1, MDataFrame# 0# n mba #) +{-# INLINE newDataFrame# #-} + +-- | Copy one DataFrame into another mutable DataFrame at specified position. +copyDataFrame# :: forall t (as :: [Nat]) (b' :: Nat) (b :: Nat) (bs :: [Nat]) (asbs :: [Nat]) s + . ( PrimBytes (DataFrame t (as +: b')) + , ConcatList as (b :+ bs) asbs + , Dimensions (b :+ bs) + ) + => DataFrame t (as +: b') -> Idx (b :+ bs) -> MDataFrame s t asbs -> State# s -> (# State# s, () #) +copyDataFrame# df ei (MDataFrame# offM _ arrM) s + | (# offA, lenA, arrA #) <- toBytes df + , elS <- elementByteSize df + , I# i <- fromEnum ei + = (# copyByteArray# arrA (offA *# elS) arrM ((offM +# i) *# elS) (lenA *# elS) s, () #) +{-# INLINE copyDataFrame# #-} + +-- | Copy one mutable DataFrame into another mutable DataFrame at specified position. +copyMDataFrame# :: forall t (as :: [Nat]) (b' :: Nat) (b :: Nat) (bs :: [Nat]) (asbs :: [Nat]) s + . ( PrimBytes t + , ConcatList as (b :+ bs) asbs + , Dimensions (b :+ bs) + ) + => MDataFrame s t (as +: b') -> Idx (b :+ bs) -> MDataFrame s t asbs -> State# s -> (# State# s, () #) +copyMDataFrame# (MDataFrame# offA lenA arrA) ei (MDataFrame# offM _ arrM) s + | elS <- elementByteSize (undefined :: t) + , I# i <- fromEnum ei + = (# copyMutableByteArray# arrA (offA *# elS) arrM ((offM +# i) *# elS) (lenA *# elS) s, () #) +{-# INLINE copyMDataFrame# #-} + +-- | Make a mutable DataFrame immutable, without copying. +unsafeFreezeDataFrame# :: forall t (ns :: [Nat]) s + . PrimBytes (DataFrame t ns) + => MDataFrame s t ns -> State# s -> (# State# s, DataFrame t ns #) +unsafeFreezeDataFrame# (MDataFrame# offM lenM arrM) s1 + | (# s2, arrA #) <- unsafeFreezeByteArray# arrM s1 + = (# s2, fromBytes (# offM, lenM, arrA #) #) +{-# INLINE unsafeFreezeDataFrame# #-} + +-- | Copy content of a mutable DataFrame into a new immutable DataFrame. +freezeDataFrame# :: forall t (ns :: [Nat]) s + . PrimBytes (DataFrame t ns) + => MDataFrame s t ns -> State# s -> (# State# s, DataFrame t ns #) +freezeDataFrame# (MDataFrame# offM n arrM) s0 + | elS <- elementByteSize (undefined :: DataFrame t ns) + , (# s1, mba #) <- newByteArray# (n *# elS) s0 + , s2 <- copyMutableByteArray# arrM (offM *# elS) mba 0# (n *# elS) s1 + , (# s3, arrA #) <- unsafeFreezeByteArray# mba s2 + = (# s3, fromBytes (# 0#, n, arrA #) #) +{-# INLINE freezeDataFrame# #-} + +-- | Create a new mutable DataFrame and copy content of immutable one in there. +thawDataFrame# :: forall t (ns :: [Nat]) s + . PrimBytes (DataFrame t ns) + => DataFrame t ns -> State# s -> (# State# s, MDataFrame s t ns #) +thawDataFrame# df s0 + | elS <- elementByteSize (undefined :: DataFrame t ns) + , (# offA, n, arrA #) <- toBytes df + , (# s1, arrM #) <- newByteArray# (n *# elS) s0 + , s2 <- copyByteArray# arrA (offA *# elS) arrM 0# (n *# elS) s1 + = (# s2, MDataFrame# 0# n arrM #) +{-# INLINE thawDataFrame# #-} + +-- | Write a single element at the specified index +writeDataFrame# :: forall t (ns :: [Nat]) s + . ( MutableFrame t ns, Dimensions ns ) + => MDataFrame s t ns -> Idx ns -> t -> State# s -> (# State# s, () #) +writeDataFrame# mdf ei x s | I# i <- fromEnum ei = (# writeDataFrameOff# mdf i x s, () #) +{-# INLINE writeDataFrame# #-} + +-- | Read a single element at the specified index +readDataFrame# :: forall t (ns :: [Nat]) s + . ( MutableFrame t ns, Dimensions ns ) + => MDataFrame s t ns -> Idx ns -> State# s -> (# State# s, t #) +readDataFrame# mdf ei | I# i <- fromEnum ei = readDataFrameOff# mdf i +{-# INLINE readDataFrame# #-} + +class MutableFrame t (ns :: [Nat]) where + -- | Write a single element at the specified element offset + writeDataFrameOff# :: MDataFrame s t ns -> Int# -> t -> State# s -> State# s + -- | Read a single element at the specified element offset + readDataFrameOff# :: MDataFrame s t ns -> Int# -> State# s -> (# State# s, t #) + +instance MutableFrame Float (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (F# x) = writeFloatArray# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readFloatArray# arr (off +# i) s0 = (# s1, F# r #) + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Double (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (D# x) = writeDoubleArray# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readDoubleArray# arr (off +# i) s0 = (# s1, D# r #) + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Int (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (I# x) = writeIntArray# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readIntArray# arr (off +# i) s0 = (# s1, I# r #) + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Int8 (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (I8# x) = writeInt8Array# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readInt8Array# arr (off +# i) s0 = (# s1, I8# r #) + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Int16 (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (I16# x) = writeInt16Array# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readInt16Array# arr (off +# i) s0 = (# s1, I16# r #) + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Int32 (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (I32# x) = writeInt32Array# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readInt32Array# arr (off +# i) s0 = (# s1, I32# r #) + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Int64 (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (I64# x) = writeInt64Array# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readInt64Array# arr (off +# i) s0 = (# s1, I64# r #) + {-# INLINE readDataFrameOff# #-} + + +instance MutableFrame Word (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (W# x) = writeWordArray# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readWordArray# arr (off +# i) s0 = (# s1, W# r #) + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Word8 (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (W8# x) = writeWord8Array# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readWord8Array# arr (off +# i) s0 = (# s1, W8# r #) + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Word16 (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (W16# x) = writeWord16Array# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readWord16Array# arr (off +# i) s0 = (# s1, W16# r #) + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Word32 (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (W32# x) = writeWord32Array# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readWord32Array# arr (off +# i) s0 = (# s1, W32# r #) + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Word64 (ns :: [Nat]) where + writeDataFrameOff# (MDataFrame# off _ arr) i (W64# x) = writeWord64Array# arr (off +# i) x + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# (MDataFrame# off _ arr) i s0 + | (# s1, r #) <- readWord64Array# arr (off +# i) s0 = (# s1, W64# r #) + {-# INLINE readDataFrameOff# #-} diff --git a/src-base/Numeric/Matrix/Mat44d.hs b/src-base/Numeric/Matrix/Mat44d.hs new file mode 100644 --- /dev/null +++ b/src-base/Numeric/Matrix/Mat44d.hs @@ -0,0 +1,37 @@ +{-# OPTIONS_GHC -fno-warn-orphans #-} +module Numeric.Matrix.Mat44d () where + + +import Numeric.Matrix.Class + +notYet :: a +notYet = error "Sorry, this function is not implemented for current platform yet." + +instance HomTransform4 Double where + translate4 = notYet + {-# INLINE translate4 #-} + translate3 = notYet + {-# INLINE translate3 #-} + rotateX = notYet + {-# INLINE rotateX #-} + rotateY = notYet + {-# INLINE rotateY #-} + rotateZ = notYet + {-# INLINE rotateZ #-} + rotate = notYet + {-# INLINE rotate #-} + rotateEuler = notYet + {-# INLINE rotateEuler #-} + lookAt = notYet + {-# INLINE lookAt #-} + perspective = notYet + {-# INLINE perspective #-} + orthogonal = notYet + {-# INLINE orthogonal #-} + toHomPoint = notYet + {-# INLINE toHomPoint #-} + toHomVector = notYet + {-# INLINE toHomVector #-} + fromHom = notYet + {-# INLINE fromHom #-} + diff --git a/src-base/Numeric/Matrix/Mat44f.hs b/src-base/Numeric/Matrix/Mat44f.hs new file mode 100644 --- /dev/null +++ b/src-base/Numeric/Matrix/Mat44f.hs @@ -0,0 +1,37 @@ +{-# OPTIONS_GHC -fno-warn-orphans #-} +module Numeric.Matrix.Mat44f () where + + +import Numeric.Matrix.Class + +notYet :: a +notYet = error "Sorry, this function is not implemented for current platform yet." + +instance HomTransform4 Float where + translate4 = notYet + {-# INLINE translate4 #-} + translate3 = notYet + {-# INLINE translate3 #-} + rotateX = notYet + {-# INLINE rotateX #-} + rotateY = notYet + {-# INLINE rotateY #-} + rotateZ = notYet + {-# INLINE rotateZ #-} + rotate = notYet + {-# INLINE rotate #-} + rotateEuler = notYet + {-# INLINE rotateEuler #-} + lookAt = notYet + {-# INLINE lookAt #-} + perspective = notYet + {-# INLINE perspective #-} + orthogonal = notYet + {-# INLINE orthogonal #-} + toHomPoint = notYet + {-# INLINE toHomPoint #-} + toHomVector = notYet + {-# INLINE toHomVector #-} + fromHom = notYet + {-# INLINE fromHom #-} + diff --git a/src-base/Numeric/Quaternion/QDouble.hs b/src-base/Numeric/Quaternion/QDouble.hs new file mode 100644 --- /dev/null +++ b/src-base/Numeric/Quaternion/QDouble.hs @@ -0,0 +1,542 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE UnboxedTuples #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +module Numeric.Quaternion.QDouble + ( QDouble, Quater (..) + ) where + +import GHC.Exts +import Data.Coerce (coerce) + +import Numeric.Array +import Numeric.DataFrame.Type +import Numeric.Commons +import Numeric.Dimensions +import Numeric.Scalar +import Numeric.Vector +import Numeric.Matrix +import qualified Numeric.DataFrame.ST as ST +import qualified Numeric.Dimensions.Traverse.ST as ST +import qualified Control.Monad.ST as ST + +import Numeric.Quaternion.Class + + +type QDouble = Quater Double + +instance Quaternion Double where + newtype Quater Double = QDouble DoubleX4 + {-# INLINE packQ #-} + packQ (D# x) (D# y) (D# z) (D# w) = QDouble (DoubleX4# x y z w) + {-# INLINE unpackQ #-} + unpackQ (QDouble (DoubleX4# x y z w)) = (D# x, D# y, D# z, D# w) + {-# INLINE fromVecNum #-} + fromVecNum (KnownDataFrame (DoubleX3# x y z)) (D# w) = QDouble (DoubleX4# x y z w) + {-# INLINE fromVec4 #-} + fromVec4 = coerce + {-# INLINE toVec4 #-} + toVec4 = coerce + {-# INLINE square #-} + square q = D# (qdot q) + {-# INLINE im #-} + im (QDouble (DoubleX4# x y z _)) = QDouble (DoubleX4# x y z 0.0##) + {-# INLINE re #-} + re (QDouble (DoubleX4# _ _ _ w)) = QDouble (DoubleX4# 0.0## 0.0## 0.0## w) + {-# INLINE imVec #-} + imVec (QDouble (DoubleX4# x y z _)) = KnownDataFrame (DoubleX3# x y z) + {-# INLINE taker #-} + taker (QDouble (DoubleX4# _ _ _ w)) = D# w + {-# INLINE takei #-} + takei (QDouble (DoubleX4# x _ _ _)) = D# x + {-# INLINE takej #-} + takej (QDouble (DoubleX4# _ y _ _)) = D# y + {-# INLINE takek #-} + takek (QDouble (DoubleX4# _ _ z _)) = D# z + {-# INLINE conjugate #-} + conjugate (QDouble (DoubleX4# x y z w)) = QDouble (DoubleX4# + (negateDouble# x) + (negateDouble# y) + (negateDouble# z) w) + {-# INLINE rotScale #-} + rotScale (QDouble (DoubleX4# i j k t)) + (KnownDataFrame (DoubleX3# x y z)) + = let l = t*##t -## i*##i -## j*##j -## k*##k + d = 2.0## *## ( i*##x +## j*##y +## k*##z) + t2 = t *## 2.0## + in KnownDataFrame + ( DoubleX3# + (l*##x +## d*##i +## t2 *## (z*##j -## y*##k)) + (l*##y +## d*##j +## t2 *## (x*##k -## z*##i)) + (l*##z +## d*##k +## t2 *## (y*##i -## x*##j)) + ) + {-# INLINE getRotScale #-} + getRotScale _ (KnownDataFrame (DoubleX3# 0.0## 0.0## 0.0##)) + = QDouble (DoubleX4# 0.0## 0.0## 0.0## 0.0##) + getRotScale (KnownDataFrame (DoubleX3# 0.0## 0.0## 0.0##)) _ + = case infty of D# x -> QDouble (DoubleX4# x x x x) + getRotScale a@(KnownDataFrame (DoubleX3# a1 a2 a3)) + b@(KnownDataFrame (DoubleX3# b1 b2 b3)) + = let ma = sqrtDouble# (a1*##a1 +## a2*##a2 +## a3*##a3) + mb = sqrtDouble# (b1*##b1 +## b2*##b2 +## b3*##b3) + d = a1*##b1 +## a2*##b2 +## a3*##b3 + c = sqrtDouble# (ma*##mb +## d) + ma2 = ma *## sqrtDouble# 2.0## + r = 1.0## /## (ma2 *## c) + in case cross a b of + KnownDataFrame (DoubleX3# 0.0## 0.0## 0.0##) -> + if isTrue# (d >## 0.0##) + then QDouble (DoubleX4# 0.0## 0.0## 0.0## (sqrtDouble# (mb /## ma))) + -- Shall we move result from k to i component? + else QDouble (DoubleX4# 0.0## 0.0## (sqrtDouble# (mb /## ma)) 0.0##) + KnownDataFrame (DoubleX3# t1 t2 t3) -> QDouble + ( DoubleX4# + (t1 *## r) + (t2 *## r) + (t3 *## r) + (c /## ma2) + ) + {-# INLINE axisRotation #-} + axisRotation (KnownDataFrame (DoubleX3# 0.0## 0.0## 0.0##)) _ + = QDouble (DoubleX4# 0.0## 0.0## 0.0## 1.0##) + axisRotation (KnownDataFrame (DoubleX3# x y z)) (D# a) + = let c = cosDouble# (a *## 0.5##) + s = sinDouble# (a *## 0.5##) + /## sqrtDouble# (x*##x +## y*##y +## z*##z) + in QDouble + ( DoubleX4# + (x *## s) + (y *## s) + (z *## s) + c + ) + {-# INLINE qArg #-} + qArg (QDouble (DoubleX4# x y z w)) + = case atan2 (D# (sqrtDouble# (x*##x +## y*##y +## z*##z))) + (D# w) of + D# a -> D# (a *## 2.0##) + {-# INLINE fromMatrix33 #-} + fromMatrix33 m + = let d = + ( ix 0# m *## ( ix 4# m *## ix 8# m -## ix 5# m *## ix 7# m ) + -## ix 1# m *## ( ix 3# m *## ix 8# m -## ix 5# m *## ix 6# m ) + +## ix 2# m *## ( ix 3# m *## ix 7# m -## ix 4# m *## ix 6# m ) + ) **## 0.33333333333333333333333333333333## + in QDouble + ( DoubleX4# + (sqrtDouble# (max# 0.0## (d +## ix 0# m -## ix 4# m -## ix 8# m )) *## sign# (ix 5# m -## ix 7# m) *## 0.5##) + (sqrtDouble# (max# 0.0## (d -## ix 0# m +## ix 4# m -## ix 8# m )) *## sign# (ix 6# m -## ix 2# m) *## 0.5##) + (sqrtDouble# (max# 0.0## (d -## ix 0# m -## ix 4# m +## ix 8# m )) *## sign# (ix 1# m -## ix 3# m) *## 0.5##) + (sqrtDouble# (max# 0.0## (d +## ix 0# m +## ix 4# m +## ix 8# m )) *## 0.5##) + ) + {-# INLINE fromMatrix44 #-} + fromMatrix44 m + = let d = + ( ix 0# m *## ( ix 5# m *## ix 10# m -## ix 6# m *## ix 9# m ) + -## ix 1# m *## ( ix 4# m *## ix 10# m -## ix 6# m *## ix 8# m ) + +## ix 2# m *## ( ix 4# m *## ix 9# m -## ix 5# m *## ix 8# m ) + ) **## 0.33333333333333333333333333333333## + c = 0.5## /## ix 15# m + in QDouble + ( DoubleX4# + (sqrtDouble# (max# 0.0## (d +## ix 0# m -## ix 5# m -## ix 10# m )) *## sign# (ix 6# m -## ix 9# m) *## c) + (sqrtDouble# (max# 0.0## (d -## ix 0# m +## ix 5# m -## ix 10# m )) *## sign# (ix 8# m -## ix 2# m) *## c) + (sqrtDouble# (max# 0.0## (d -## ix 0# m -## ix 5# m +## ix 10# m )) *## sign# (ix 1# m -## ix 4# m) *## c) + (sqrtDouble# (max# 0.0## (d +## ix 0# m +## ix 5# m +## ix 10# m )) *## c) + ) + {-# INLINE toMatrix33 #-} + toMatrix33 (QDouble (DoubleX4# 0.0## 0.0## 0.0## w)) = diag (scalar (D# (w *## w))) + toMatrix33 (QDouble (DoubleX4# x' y' z' w')) = + let x = scalar (D# x') + y = scalar (D# y') + z = scalar (D# z') + w = scalar (D# w') + x2 = x * x + y2 = y * y + z2 = z * z + w2 = w * w + l2 = x2 + y2 + z2 + w2 + in ST.runST $ do + df <- ST.newDataFrame + ST.writeDataFrameOff df 0 $ l2 - 2*(z2 + y2) + ST.writeDataFrameOff df 1 $ 2*(x*y + z*w) + ST.writeDataFrameOff df 2 $ 2*(x*z - y*w) + ST.writeDataFrameOff df 3 $ 2*(x*y - z*w) + ST.writeDataFrameOff df 4 $ l2 - 2*(z2 + x2) + ST.writeDataFrameOff df 5 $ 2*(y*z + x*w) + ST.writeDataFrameOff df 6 $ 2*(x*z + y*w) + ST.writeDataFrameOff df 7 $ 2*(y*z - x*w) + ST.writeDataFrameOff df 8 $ l2 - 2*(y2 + x2) + ST.unsafeFreezeDataFrame df + {-# INLINE toMatrix44 #-} + toMatrix44 (QDouble (DoubleX4# 0.0## 0.0## 0.0## w)) = ST.runST $ do + df <- ST.newDataFrame + ST.overDimOff_ (dim :: Dim '[4,4]) (\i -> ST.writeDataFrameOff df (I# i) 0) 0# 1# + let w2 = scalar (D# (w *## w)) + ST.writeDataFrameOff df 0 w2 + ST.writeDataFrameOff df 5 w2 + ST.writeDataFrameOff df 10 w2 + ST.writeDataFrameOff df 15 1 + ST.unsafeFreezeDataFrame df + toMatrix44 (QDouble (DoubleX4# x' y' z' w')) = + let x = scalar (D# x') + y = scalar (D# y') + z = scalar (D# z') + w = scalar (D# w') + x2 = x * x + y2 = y * y + z2 = z * z + w2 = w * w + l2 = x2 + y2 + z2 + w2 + in ST.runST $ do + df <- ST.newDataFrame + ST.writeDataFrameOff df 0 $ l2 - 2*(z2 + y2) + ST.writeDataFrameOff df 1 $ 2*(x*y + z*w) + ST.writeDataFrameOff df 2 $ 2*(x*z - y*w) + ST.writeDataFrameOff df 3 0 + ST.writeDataFrameOff df 4 $ 2*(x*y - z*w) + ST.writeDataFrameOff df 5 $ l2 - 2*(z2 + x2) + ST.writeDataFrameOff df 6 $ 2*(y*z + x*w) + ST.writeDataFrameOff df 7 0 + ST.writeDataFrameOff df 8 $ 2*(x*z + y*w) + ST.writeDataFrameOff df 9 $ 2*(y*z - x*w) + ST.writeDataFrameOff df 10 $ l2 - 2*(y2 + x2) + ST.writeDataFrameOff df 11 0 + ST.writeDataFrameOff df 12 0 + ST.writeDataFrameOff df 13 0 + ST.writeDataFrameOff df 14 0 + ST.writeDataFrameOff df 15 1 + ST.unsafeFreezeDataFrame df + +qdot :: QDouble -> Double# +qdot (QDouble (DoubleX4# x y z w)) = (x *## x) +## + (y *## y) +## + (z *## z) +## + (w *## w) +{-# INLINE qdot #-} + +infty :: Double +infty = read "Infinity" + +max# :: Double# -> Double# -> Double# +max# a b | isTrue# (a >## b) = a + | otherwise = b +{-# INLINE max# #-} + +sign# :: Double# -> Double# +sign# a | isTrue# (a >## 0.0##) = 1.0## + | isTrue# (a <## 0.0##) = negateDouble# 1.0## + | otherwise = 0.0## +{-# INLINE sign# #-} + +-------------------------------------------------------------------------- +-- Num +-------------------------------------------------------------------------- + +instance Num QDouble where + QDouble a + QDouble b + = QDouble (a + b) + {-# INLINE (+) #-} + QDouble a - QDouble b + = QDouble (a - b) + {-# INLINE (-) #-} + QDouble (DoubleX4# a1 a2 a3 a4) * QDouble (DoubleX4# b1 b2 b3 b4) + = QDouble + ( DoubleX4# + ((a4 *## b1) +## + (a1 *## b4) +## + (a2 *## b3) -## + (a3 *## b2) + ) + ((a4 *## b2) -## + (a1 *## b3) +## + (a2 *## b4) +## + (a3 *## b1) + ) + ((a4 *## b3) +## + (a1 *## b2) -## + (a2 *## b1) +## + (a3 *## b4) + ) + ((a4 *## b4) -## + (a1 *## b1) -## + (a2 *## b2) -## + (a3 *## b3) + ) + ) + {-# INLINE (*) #-} + negate (QDouble a) = QDouble (negate a) + {-# INLINE negate #-} + abs q = QDouble (DoubleX4# 0.0## 0.0## 0.0## (sqrtDouble# (qdot q))) + {-# INLINE abs #-} + signum q@(QDouble (DoubleX4# x y z w)) + = case qdot q of + 0.0## -> QDouble (DoubleX4# 0.0## 0.0## 0.0## 0.0##) + qd -> case 1.0## /## sqrtDouble# qd of + s -> QDouble + ( DoubleX4# + (x *## s) + (y *## s) + (z *## s) + (w *## s) + ) + {-# INLINE signum #-} + fromInteger n = case fromInteger n of + D# x -> QDouble (DoubleX4# 0.0## 0.0## 0.0## x) + {-# INLINE fromInteger #-} + + + +-------------------------------------------------------------------------- +-- Fractional +-------------------------------------------------------------------------- + +instance Fractional QDouble where + {-# INLINE recip #-} + recip q@(QDouble (DoubleX4# x y z w)) = case -1.0## /## qdot q of + c -> QDouble + ( DoubleX4# + (x *## c) + (y *## c) + (z *## c) + (negateDouble# (w *## c)) + ) + {-# INLINE (/) #-} + a / b = a * recip b + {-# INLINE fromRational #-} + fromRational q = case fromRational q of + D# x -> QDouble (DoubleX4# 0.0## 0.0## 0.0## x) + +-------------------------------------------------------------------------- +-- Doubleing +-------------------------------------------------------------------------- + +instance Floating QDouble where + {-# INLINE pi #-} + pi = QDouble (DoubleX4# 0.0## 0.0## 0.0## + 3.141592653589793##) + {-# INLINE exp #-} + exp (QDouble (DoubleX4# x y z w)) + = case (# (x *## x) +## + (y *## y) +## + (z *## z) + , expDouble# w + #) of + (# 0.0##, et #) -> QDouble (DoubleX4# 0.0## 0.0## 0.0## et) + (# mv2, et #) -> case sqrtDouble# mv2 of + mv -> case et *## sinDouble# mv + /## mv of + l -> QDouble + ( DoubleX4# + (x *## l) + (y *## l) + (z *## l) + (et *## cosDouble# mv) + ) + {-# INLINE log #-} + log (QDouble (DoubleX4# x y z w)) + = case (x *## x) +## + (y *## y) +## + (z *## z) of + 0.0## -> if isTrue# (w >=## 0.0##) + then QDouble (DoubleX4# 0.0## 0.0## 0.0## (logDouble# w)) + else QDouble (DoubleX4# 3.141592653589793## 0.0## 0.0## + (logDouble# (negateDouble# w))) + mv2 -> case (# sqrtDouble# (mv2 +## (w *## w)) + , sqrtDouble# mv2 + #) of + (# mq, mv #) -> case atan2 (D# mv) (D# w) / D# mv of + D# l -> QDouble + ( DoubleX4# + (x *## l) + (y *## l) + (z *## l) + (logDouble# mq) + ) + {-# INLINE sqrt #-} + sqrt (QDouble (DoubleX4# x y z w)) + = case (x *## x) +## + (y *## y) +## + (z *## z) of + 0.0## -> if isTrue# (w >=## 0.0##) + then QDouble (DoubleX4# 0.0## 0.0## 0.0## (sqrtDouble# w)) + else QDouble (DoubleX4# (sqrtDouble# (negateDouble# w)) 0.0## 0.0## 0.0##) + mv2 -> + let mq = sqrtDouble# (mv2 +## w *## w) + l2 = sqrtDouble# mq + tq = w /## (mq *## 2.0##) + sina = sqrtDouble# (0.5## -## tq) *## l2 /## sqrtDouble# mv2 + in QDouble + ( DoubleX4# + (x *## sina) + (y *## sina) + (z *## sina) + (sqrtDouble# (0.5## +## tq) *## l2) + ) + {-# INLINE sin #-} + sin (QDouble (DoubleX4# x y z w)) + = case (x *## x) +## + (y *## y) +## + (z *## z) of + 0.0## -> QDouble (DoubleX4# 0.0## 0.0## 0.0## (sinDouble# w)) + mv2 -> case sqrtDouble# mv2 of + mv -> case cosDouble# w *## sinhDouble# mv + /## mv of + l -> QDouble + ( DoubleX4# + (x *## l) + (y *## l) + (z *## l) + (sinDouble# w *## coshDouble# mv) + ) + {-# INLINE cos #-} + cos (QDouble (DoubleX4# x y z w)) + = case (x *## x) +## + (y *## y) +## + (z *## z) of + 0.0## -> QDouble (DoubleX4# 0.0## 0.0## 0.0## (cosDouble# w)) + mv2 -> case sqrtDouble# mv2 of + mv -> case sinDouble# w *## sinhDouble# mv + /## negateDouble# mv of + l -> QDouble + ( DoubleX4# + (x *## l) + (y *## l) + (z *## l) + (cosDouble# w *## coshDouble# mv) + ) + {-# INLINE tan #-} + tan (QDouble (DoubleX4# x y z w)) + = case (x *## x) +## + (y *## y) +## + (z *## z) of + 0.0## -> QDouble (DoubleX4# 0.0## 0.0## 0.0## (tanDouble# w)) + mv2 -> + let mv = sqrtDouble# mv2 + chv = coshDouble# mv + shv = sinhDouble# mv + ct = cosDouble# w + st = sinDouble# w + cq = 1.0## /## + ( (ct *## ct *## chv *## chv) + +## + (st *## st *## shv *## shv) + ) + l = chv *## shv *## cq + /## mv + in QDouble + ( DoubleX4# + (x *## l) + (y *## l) + (z *## l) + (ct *## st *## cq) + ) + {-# INLINE sinh #-} + sinh (QDouble (DoubleX4# x y z w)) + = case (x *## x) +## + (y *## y) +## + (z *## z) of + 0.0## -> QDouble (DoubleX4# 0.0## 0.0## 0.0## (sinhDouble# w)) + mv2 -> case sqrtDouble# mv2 of + mv -> case coshDouble# w *## sinDouble# mv + /## mv of + l -> QDouble + ( DoubleX4# + (x *## l) + (y *## l) + (z *## l) + (sinhDouble# w *## cosDouble# mv) + ) + {-# INLINE cosh #-} + cosh (QDouble (DoubleX4# x y z w)) + = case (x *## x) +## + (y *## y) +## + (z *## z) of + 0.0## -> QDouble (DoubleX4# 0.0## 0.0## 0.0## (coshDouble# w)) + mv2 -> case sqrtDouble# mv2 of + mv -> case sinhDouble# w *## sinDouble# mv + /## mv of + l -> QDouble + ( DoubleX4# + (x *## l) + (y *## l) + (z *## l) + (coshDouble# w *## cosDouble# mv) + ) + {-# INLINE tanh #-} + tanh (QDouble (DoubleX4# x y z w)) + = case (x *## x) +## + (y *## y) +## + (z *## z) of + 0.0## -> QDouble (DoubleX4# 0.0## 0.0## 0.0## (tanhDouble# w)) + mv2 -> + let mv = sqrtDouble# mv2 + cv = cosDouble# mv + sv = sinDouble# mv + cht = coshDouble# w + sht = sinhDouble# w + cq = 1.0## /## + ( (cht *## cht *## cv *## cv) + +## + (sht *## sht *## sv *## sv) + ) + l = cv *## sv *## cq + /## mv + in QDouble + ( DoubleX4# + (x *## l) + (y *## l) + (z *## l) + (cht *## sht *## cq) + ) + {-# INLINE asin #-} + asin q = -i * log (i*q + sqrt (1 - q*q)) + where + i = case signum . im $ q of + 0 -> QDouble (DoubleX4# 1.0## 0.0## 0.0## 0.0##) + i' -> i' + {-# INLINE acos #-} + acos q = pi/2 - asin q + {-# INLINE atan #-} + atan q@(QDouble (DoubleX4# _ _ _ w)) + = if square imq == 0 + then QDouble (DoubleX4# 0.0## 0.0## 0.0## (atanDouble# w)) + else i / 2 * log ( (i + q) / (i - q) ) + where + i = signum imq + imq = im q + {-# INLINE asinh #-} + asinh q = log (q + sqrt (q*q + 1)) + {-# INLINE acosh #-} + acosh q = log (q + sqrt (q*q - 1)) + {-# INLINE atanh #-} + atanh q = 0.5 * log ((1+q)/(1-q)) + +-------------------------------------------------------------------------- +-- Eq +-------------------------------------------------------------------------- + +instance Eq QDouble where + {-# INLINE (==) #-} + QDouble a == QDouble b = a == b + {-# INLINE (/=) #-} + QDouble a /= QDouble b = a /= b + + + +-------------------------------------------------------------------------- +-- Show +-------------------------------------------------------------------------- + +instance Show QDouble where + show (QDouble (DoubleX4# x y z w)) = + show (D# w) ++ ss x ++ "i" + ++ ss y ++ "j" + ++ ss z ++ "k" + where + ss a# = case D# a# of + a -> if a >= 0 then " + " ++ show a + else " - " ++ show (negate a) diff --git a/src-base/Numeric/Quaternion/QFloat.hs b/src-base/Numeric/Quaternion/QFloat.hs new file mode 100644 --- /dev/null +++ b/src-base/Numeric/Quaternion/QFloat.hs @@ -0,0 +1,562 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE UnboxedTuples #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +module Numeric.Quaternion.QFloat + ( QFloat, Quater (..) + ) where + +import GHC.Exts +import Data.Coerce (coerce) + +import Numeric.Array +import Numeric.DataFrame.Type +import Numeric.Commons +import Numeric.Dimensions +import Numeric.Scalar +import Numeric.Vector +import Numeric.Matrix +import qualified Numeric.DataFrame.ST as ST +import qualified Numeric.Dimensions.Traverse.ST as ST +import qualified Control.Monad.ST as ST + +import Numeric.Quaternion.Class + + +type QFloat = Quater Float + +instance Quaternion Float where + newtype Quater Float = QFloat FloatX4 + {-# INLINE packQ #-} + packQ (F# x) (F# y) (F# z) (F# w) = QFloat (FloatX4# x y z w) + {-# INLINE unpackQ #-} + unpackQ (QFloat (FloatX4# x y z w)) = (F# x, F# y, F# z, F# w) + {-# INLINE fromVecNum #-} + fromVecNum (KnownDataFrame (FloatX3# x y z)) (F# w) = QFloat (FloatX4# x y z w) + {-# INLINE fromVec4 #-} + fromVec4 = coerce + {-# INLINE toVec4 #-} + toVec4 = coerce + {-# INLINE square #-} + square q = F# (qdot q) + {-# INLINE im #-} + im (QFloat (FloatX4# x y z _)) = QFloat (FloatX4# x y z 0.0#) + {-# INLINE re #-} + re (QFloat (FloatX4# _ _ _ w)) = QFloat (FloatX4# 0.0# 0.0# 0.0# w) + {-# INLINE imVec #-} + imVec (QFloat (FloatX4# x y z _)) = KnownDataFrame (FloatX3# x y z) + {-# INLINE taker #-} + taker (QFloat (FloatX4# _ _ _ w)) = F# w + {-# INLINE takei #-} + takei (QFloat (FloatX4# x _ _ _)) = F# x + {-# INLINE takej #-} + takej (QFloat (FloatX4# _ y _ _)) = F# y + {-# INLINE takek #-} + takek (QFloat (FloatX4# _ _ z _)) = F# z + {-# INLINE conjugate #-} + conjugate (QFloat (FloatX4# x y z w)) = QFloat (FloatX4# + (negateFloat# x) + (negateFloat# y) + (negateFloat# z) w) + {-# INLINE rotScale #-} + rotScale (QFloat (FloatX4# i j k t)) + (KnownDataFrame (FloatX3# x y z)) + = let l = t*%t -% i*%i -% j*%j -% k*%k + d = 2.0# *% ( i*%x +% j*%y +% k*%z) + t2 = t *% 2.0# + in KnownDataFrame + ( FloatX3# + (l*%x +% d*%i +% t2 *% (z*%j -% y*%k)) + (l*%y +% d*%j +% t2 *% (x*%k -% z*%i)) + (l*%z +% d*%k +% t2 *% (y*%i -% x*%j)) + ) + {-# INLINE getRotScale #-} + getRotScale _ (KnownDataFrame (FloatX3# 0.0# 0.0# 0.0#)) + = QFloat (FloatX4# 0.0# 0.0# 0.0# 0.0#) + getRotScale (KnownDataFrame (FloatX3# 0.0# 0.0# 0.0#)) _ + = case infty of F# x -> QFloat (FloatX4# x x x x) + getRotScale a@(KnownDataFrame (FloatX3# a1 a2 a3)) + b@(KnownDataFrame (FloatX3# b1 b2 b3)) + = let ma = sqrtFloat# (a1*%a1 +% a2*%a2 +% a3*%a3) + mb = sqrtFloat# (b1*%b1 +% b2*%b2 +% b3*%b3) + d = a1*%b1 +% a2*%b2 +% a3*%b3 + c = sqrtFloat# (ma*%mb +% d) + ma2 = ma *% sqrtFloat# 2.0# + r = 1.0# /% (ma2 *% c) + in case cross a b of + KnownDataFrame (FloatX3# 0.0# 0.0# 0.0#) -> + if isTrue# (gtFloat# d 0.0#) + then QFloat (FloatX4# 0.0# 0.0# 0.0# (sqrtFloat# (mb /% ma))) + -- Shall we move result from k to i component? + else QFloat (FloatX4# 0.0# 0.0# (sqrtFloat# (mb /% ma)) 0.0#) + KnownDataFrame (FloatX3# t1 t2 t3) -> QFloat + ( FloatX4# + (t1 *% r) + (t2 *% r) + (t3 *% r) + (c /% ma2) + ) + {-# INLINE axisRotation #-} + axisRotation (KnownDataFrame (FloatX3# 0.0# 0.0# 0.0#)) _ + = QFloat (FloatX4# 0.0# 0.0# 0.0# 1.0#) + axisRotation (KnownDataFrame (FloatX3# x y z)) (F# a) + = let c = cosFloat# (a *% 0.5#) + s = sinFloat# (a *% 0.5#) + /% sqrtFloat# (x*%x +% y*%y +% z*%z) + in QFloat + ( FloatX4# + (x *% s) + (y *% s) + (z *% s) + c + ) + {-# INLINE qArg #-} + qArg (QFloat (FloatX4# x y z w)) + = case atan2 (F# (sqrtFloat# (x*%x +% y*%y +% z*%z))) + (F# w) of + F# a -> F# (a *% 2.0#) + {-# INLINE fromMatrix33 #-} + fromMatrix33 m + = let d = powerFloat# + ( ix 0# m *% ( ix 4# m *% ix 8# m -% ix 5# m *% ix 7# m ) + -% ix 1# m *% ( ix 3# m *% ix 8# m -% ix 5# m *% ix 6# m ) + +% ix 2# m *% ( ix 3# m *% ix 7# m -% ix 4# m *% ix 6# m ) + ) 0.33333333333333333333333333333333# + in QFloat + ( FloatX4# + (sqrtFloat# (max# 0.0# (d +% ix 0# m -% ix 4# m -% ix 8# m )) *% sign# (ix 5# m -% ix 7# m) *% 0.5#) + (sqrtFloat# (max# 0.0# (d -% ix 0# m +% ix 4# m -% ix 8# m )) *% sign# (ix 6# m -% ix 2# m) *% 0.5#) + (sqrtFloat# (max# 0.0# (d -% ix 0# m -% ix 4# m +% ix 8# m )) *% sign# (ix 1# m -% ix 3# m) *% 0.5#) + (sqrtFloat# (max# 0.0# (d +% ix 0# m +% ix 4# m +% ix 8# m )) *% 0.5#) + ) + {-# INLINE fromMatrix44 #-} + fromMatrix44 m + = let d = powerFloat# + ( ix 0# m *% ( ix 5# m *% ix 10# m -% ix 6# m *% ix 9# m ) + -% ix 1# m *% ( ix 4# m *% ix 10# m -% ix 6# m *% ix 8# m ) + +% ix 2# m *% ( ix 4# m *% ix 9# m -% ix 5# m *% ix 8# m ) + ) 0.33333333333333333333333333333333# + c = 0.5# /% ix 15# m + in QFloat + ( FloatX4# + (sqrtFloat# (max# 0.0# (d +% ix 0# m -% ix 5# m -% ix 10# m )) *% sign# (ix 6# m -% ix 9# m) *% c) + (sqrtFloat# (max# 0.0# (d -% ix 0# m +% ix 5# m -% ix 10# m )) *% sign# (ix 8# m -% ix 2# m) *% c) + (sqrtFloat# (max# 0.0# (d -% ix 0# m -% ix 5# m +% ix 10# m )) *% sign# (ix 1# m -% ix 4# m) *% c) + (sqrtFloat# (max# 0.0# (d +% ix 0# m +% ix 5# m +% ix 10# m )) *% c) + ) + {-# INLINE toMatrix33 #-} + toMatrix33 (QFloat (FloatX4# 0.0# 0.0# 0.0# w)) = diag (scalar (F# (w *% w))) + toMatrix33 (QFloat (FloatX4# x' y' z' w')) = + let x = scalar (F# x') + y = scalar (F# y') + z = scalar (F# z') + w = scalar (F# w') + x2 = x * x + y2 = y * y + z2 = z * z + w2 = w * w + l2 = x2 + y2 + z2 + w2 + in ST.runST $ do + df <- ST.newDataFrame + ST.writeDataFrameOff df 0 $ l2 - 2*(z2 + y2) + ST.writeDataFrameOff df 1 $ 2*(x*y + z*w) + ST.writeDataFrameOff df 2 $ 2*(x*z - y*w) + ST.writeDataFrameOff df 3 $ 2*(x*y - z*w) + ST.writeDataFrameOff df 4 $ l2 - 2*(z2 + x2) + ST.writeDataFrameOff df 5 $ 2*(y*z + x*w) + ST.writeDataFrameOff df 6 $ 2*(x*z + y*w) + ST.writeDataFrameOff df 7 $ 2*(y*z - x*w) + ST.writeDataFrameOff df 8 $ l2 - 2*(y2 + x2) + ST.unsafeFreezeDataFrame df + {-# INLINE toMatrix44 #-} + toMatrix44 (QFloat (FloatX4# 0.0# 0.0# 0.0# w)) = ST.runST $ do + df <- ST.newDataFrame + ST.overDimOff_ (dim :: Dim '[4,4]) (\i -> ST.writeDataFrameOff df (I# i) 0) 0# 1# + let w2 = scalar (F# (w *% w)) + ST.writeDataFrameOff df 0 w2 + ST.writeDataFrameOff df 5 w2 + ST.writeDataFrameOff df 10 w2 + ST.writeDataFrameOff df 15 1 + ST.unsafeFreezeDataFrame df + toMatrix44 (QFloat (FloatX4# x' y' z' w')) = + let x = scalar (F# x') + y = scalar (F# y') + z = scalar (F# z') + w = scalar (F# w') + x2 = x * x + y2 = y * y + z2 = z * z + w2 = w * w + l2 = x2 + y2 + z2 + w2 + in ST.runST $ do + df <- ST.newDataFrame + ST.writeDataFrameOff df 0 $ l2 - 2*(z2 + y2) + ST.writeDataFrameOff df 1 $ 2*(x*y + z*w) + ST.writeDataFrameOff df 2 $ 2*(x*z - y*w) + ST.writeDataFrameOff df 3 0 + ST.writeDataFrameOff df 4 $ 2*(x*y - z*w) + ST.writeDataFrameOff df 5 $ l2 - 2*(z2 + x2) + ST.writeDataFrameOff df 6 $ 2*(y*z + x*w) + ST.writeDataFrameOff df 7 0 + ST.writeDataFrameOff df 8 $ 2*(x*z + y*w) + ST.writeDataFrameOff df 9 $ 2*(y*z - x*w) + ST.writeDataFrameOff df 10 $ l2 - 2*(y2 + x2) + ST.writeDataFrameOff df 11 0 + ST.writeDataFrameOff df 12 0 + ST.writeDataFrameOff df 13 0 + ST.writeDataFrameOff df 14 0 + ST.writeDataFrameOff df 15 1 + ST.unsafeFreezeDataFrame df + +qdot :: QFloat -> Float# +qdot (QFloat (FloatX4# x y z w)) = (x *% x) +% + (y *% y) +% + (z *% z) +% + (w *% w) +{-# INLINE qdot #-} + +(*%) :: Float# -> Float# -> Float# +(*%) = timesFloat# +{-# INLINE (*%) #-} +infixl 7 *% + +(-%) :: Float# -> Float# -> Float# +(-%) = minusFloat# +{-# INLINE (-%) #-} +infixl 6 -% + +(+%) :: Float# -> Float# -> Float# +(+%) = plusFloat# +{-# INLINE (+%) #-} +infixl 6 +% + +(/%) :: Float# -> Float# -> Float# +(/%) = divideFloat# +{-# INLINE (/%) #-} +infixl 7 /% + +infty :: Float +infty = read "Infinity" + +max# :: Float# -> Float# -> Float# +max# a b | isTrue# (gtFloat# a b) = a + | otherwise = b +{-# INLINE max# #-} + +sign# :: Float# -> Float# +sign# a | isTrue# (gtFloat# a 0.0#) = 1.0# + | isTrue# (ltFloat# a 0.0#) = negateFloat# 1.0# + | otherwise = 0.0# +{-# INLINE sign# #-} + +-------------------------------------------------------------------------- +-- Num +-------------------------------------------------------------------------- + +instance Num QFloat where + QFloat a + QFloat b + = QFloat (a + b) + {-# INLINE (+) #-} + QFloat a - QFloat b + = QFloat (a - b) + {-# INLINE (-) #-} + QFloat (FloatX4# a1 a2 a3 a4) * QFloat (FloatX4# b1 b2 b3 b4) + = QFloat + ( FloatX4# + ((a4 *% b1) +% + (a1 *% b4) +% + (a2 *% b3) -% + (a3 *% b2) + ) + ((a4 *% b2) -% + (a1 *% b3) +% + (a2 *% b4) +% + (a3 *% b1) + ) + ((a4 *% b3) +% + (a1 *% b2) -% + (a2 *% b1) +% + (a3 *% b4) + ) + ((a4 *% b4) -% + (a1 *% b1) -% + (a2 *% b2) -% + (a3 *% b3) + ) + ) + {-# INLINE (*) #-} + negate (QFloat a) = QFloat (negate a) + {-# INLINE negate #-} + abs q = QFloat (FloatX4# 0.0# 0.0# 0.0# (sqrtFloat# (qdot q))) + {-# INLINE abs #-} + signum q@(QFloat (FloatX4# x y z w)) + = case qdot q of + 0.0# -> QFloat (FloatX4# 0.0# 0.0# 0.0# 0.0#) + qd -> case 1.0# /% sqrtFloat# qd of + s -> QFloat + ( FloatX4# + (x *% s) + (y *% s) + (z *% s) + (w *% s) + ) + {-# INLINE signum #-} + fromInteger n = case fromInteger n of + F# x -> QFloat (FloatX4# 0.0# 0.0# 0.0# x) + {-# INLINE fromInteger #-} + + + +-------------------------------------------------------------------------- +-- Fractional +-------------------------------------------------------------------------- + +instance Fractional QFloat where + {-# INLINE recip #-} + recip q@(QFloat (FloatX4# x y z w)) = case -1.0# /% qdot q of + c -> QFloat + ( FloatX4# + (x *% c) + (y *% c) + (z *% c) + (negateFloat# (w *% c)) + ) + {-# INLINE (/) #-} + a / b = a * recip b + {-# INLINE fromRational #-} + fromRational q = case fromRational q of + F# x -> QFloat (FloatX4# 0.0# 0.0# 0.0# x) + +-------------------------------------------------------------------------- +-- Floating +-------------------------------------------------------------------------- + +instance Floating QFloat where + {-# INLINE pi #-} + pi = QFloat (FloatX4# 0.0# 0.0# 0.0# + 3.141592653589793#) + {-# INLINE exp #-} + exp (QFloat (FloatX4# x y z w)) + = case (# (x *% x) +% + (y *% y) +% + (z *% z) + , expFloat# w + #) of + (# 0.0#, et #) -> QFloat (FloatX4# 0.0# 0.0# 0.0# et) + (# mv2, et #) -> case sqrtFloat# mv2 of + mv -> case et *% sinFloat# mv + /% mv of + l -> QFloat + ( FloatX4# + (x *% l) + (y *% l) + (z *% l) + (et *% cosFloat# mv) + ) + {-# INLINE log #-} + log (QFloat (FloatX4# x y z w)) + = case (x *% x) +% + (y *% y) +% + (z *% z) of + 0.0# -> if isTrue# (w `geFloat#` 0.0#) + then QFloat (FloatX4# 0.0# 0.0# 0.0# (logFloat# w)) + else QFloat (FloatX4# 3.141592653589793# 0.0# 0.0# + (logFloat# (negateFloat# w))) + mv2 -> case (# sqrtFloat# (mv2 +% (w *% w)) + , sqrtFloat# mv2 + #) of + (# mq, mv #) -> case atan2 (F# mv) (F# w) / F# mv of + F# l -> QFloat + ( FloatX4# + (x *% l) + (y *% l) + (z *% l) + (logFloat# mq) + ) + {-# INLINE sqrt #-} + sqrt (QFloat (FloatX4# x y z w)) + = case (x *% x) +% + (y *% y) +% + (z *% z) of + 0.0# -> if isTrue# (w `geFloat#` 0.0#) + then QFloat (FloatX4# 0.0# 0.0# 0.0# (sqrtFloat# w)) + else QFloat (FloatX4# (sqrtFloat# (negateFloat# w)) 0.0# 0.0# 0.0#) + mv2 -> + let mq = sqrtFloat# (mv2 +% w *% w) + l2 = sqrtFloat# mq + tq = w /% (mq *% 2.0#) + sina = sqrtFloat# (0.5# -% tq) *% l2 /% sqrtFloat# mv2 + in QFloat + ( FloatX4# + (x *% sina) + (y *% sina) + (z *% sina) + (sqrtFloat# (0.5# +% tq) *% l2) + ) + {-# INLINE sin #-} + sin (QFloat (FloatX4# x y z w)) + = case (x *% x) +% + (y *% y) +% + (z *% z) of + 0.0# -> QFloat (FloatX4# 0.0# 0.0# 0.0# (sinFloat# w)) + mv2 -> case sqrtFloat# mv2 of + mv -> case cosFloat# w *% sinhFloat# mv + /% mv of + l -> QFloat + ( FloatX4# + (x *% l) + (y *% l) + (z *% l) + (sinFloat# w *% coshFloat# mv) + ) + {-# INLINE cos #-} + cos (QFloat (FloatX4# x y z w)) + = case (x *% x) +% + (y *% y) +% + (z *% z) of + 0.0# -> QFloat (FloatX4# 0.0# 0.0# 0.0# (cosFloat# w)) + mv2 -> case sqrtFloat# mv2 of + mv -> case sinFloat# w *% sinhFloat# mv + /% negateFloat# mv of + l -> QFloat + ( FloatX4# + (x *% l) + (y *% l) + (z *% l) + (cosFloat# w *% coshFloat# mv) + ) + {-# INLINE tan #-} + tan (QFloat (FloatX4# x y z w)) + = case (x *% x) +% + (y *% y) +% + (z *% z) of + 0.0# -> QFloat (FloatX4# 0.0# 0.0# 0.0# (tanFloat# w)) + mv2 -> + let mv = sqrtFloat# mv2 + chv = coshFloat# mv + shv = sinhFloat# mv + ct = cosFloat# w + st = sinFloat# w + cq = 1.0# /% + ( (ct *% ct *% chv *% chv) + +% + (st *% st *% shv *% shv) + ) + l = chv *% shv *% cq + /% mv + in QFloat + ( FloatX4# + (x *% l) + (y *% l) + (z *% l) + (ct *% st *% cq) + ) + {-# INLINE sinh #-} + sinh (QFloat (FloatX4# x y z w)) + = case (x *% x) +% + (y *% y) +% + (z *% z) of + 0.0# -> QFloat (FloatX4# 0.0# 0.0# 0.0# (sinhFloat# w)) + mv2 -> case sqrtFloat# mv2 of + mv -> case coshFloat# w *% sinFloat# mv + /% mv of + l -> QFloat + ( FloatX4# + (x *% l) + (y *% l) + (z *% l) + (sinhFloat# w *% cosFloat# mv) + ) + {-# INLINE cosh #-} + cosh (QFloat (FloatX4# x y z w)) + = case (x *% x) +% + (y *% y) +% + (z *% z) of + 0.0# -> QFloat (FloatX4# 0.0# 0.0# 0.0# (coshFloat# w)) + mv2 -> case sqrtFloat# mv2 of + mv -> case sinhFloat# w *% sinFloat# mv + /% mv of + l -> QFloat + ( FloatX4# + (x *% l) + (y *% l) + (z *% l) + (coshFloat# w *% cosFloat# mv) + ) + {-# INLINE tanh #-} + tanh (QFloat (FloatX4# x y z w)) + = case (x *% x) +% + (y *% y) +% + (z *% z) of + 0.0# -> QFloat (FloatX4# 0.0# 0.0# 0.0# (tanhFloat# w)) + mv2 -> + let mv = sqrtFloat# mv2 + cv = cosFloat# mv + sv = sinFloat# mv + cht = coshFloat# w + sht = sinhFloat# w + cq = 1.0# /% + ( (cht *% cht *% cv *% cv) + +% + (sht *% sht *% sv *% sv) + ) + l = cv *% sv *% cq + /% mv + in QFloat + ( FloatX4# + (x *% l) + (y *% l) + (z *% l) + (cht *% sht *% cq) + ) + {-# INLINE asin #-} + asin q = -i * log (i*q + sqrt (1 - q*q)) + where + i = case signum . im $ q of + 0 -> QFloat (FloatX4# 1.0# 0.0# 0.0# 0.0#) + i' -> i' + {-# INLINE acos #-} + acos q = pi/2 - asin q + {-# INLINE atan #-} + atan q@(QFloat (FloatX4# _ _ _ w)) + = if square imq == 0 + then QFloat (FloatX4# 0.0# 0.0# 0.0# (atanFloat# w)) + else i / 2 * log ( (i + q) / (i - q) ) + where + i = signum imq + imq = im q + {-# INLINE asinh #-} + asinh q = log (q + sqrt (q*q + 1)) + {-# INLINE acosh #-} + acosh q = log (q + sqrt (q*q - 1)) + {-# INLINE atanh #-} + atanh q = 0.5 * log ((1+q)/(1-q)) + +-------------------------------------------------------------------------- +-- Eq +-------------------------------------------------------------------------- + +instance Eq QFloat where + {-# INLINE (==) #-} + QFloat a == QFloat b = a == b + {-# INLINE (/=) #-} + QFloat a /= QFloat b = a /= b + + + +-------------------------------------------------------------------------- +-- Show +-------------------------------------------------------------------------- + +instance Show QFloat where + show (QFloat (FloatX4# x y z w)) = + show (F# w) ++ ss x ++ "i" + ++ ss y ++ "j" + ++ ss z ++ "k" + where + ss a# = case F# a# of + a -> if a >= 0 then " + " ++ show a + else " - " ++ show (negate a) diff --git a/src-ghcjs/Numeric/Array/Family.hs b/src-ghcjs/Numeric/Array/Family.hs --- a/src-ghcjs/Numeric/Array/Family.hs +++ b/src-ghcjs/Numeric/Array/Family.hs @@ -14,6 +14,8 @@ {-# LANGUAGE UnboxedTuples #-} {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE JavaScriptFFI #-} +{-# LANGUAGE UnliftedFFITypes #-} ----------------------------------------------------------------------------- -- | -- Module : Numeric.Array.Family @@ -38,7 +40,7 @@ import Data.Type.Equality ((:~:) (..)) import Data.Word (Word16, Word32, Word8) import GHC.Prim (Double#, Float#, Int#, - Word#, unsafeCoerce#) + Word#, unsafeCoerce#, ByteArray#) import GHC.Types (Int (..)) import GHCJS.Types @@ -74,9 +76,9 @@ type instance ElemRep Word8Clamped = ElemRep Int type instance ElemPrim Word8Clamped = Int# instance PrimBytes Word8Clamped where - toBytes (Clamped i) = toBytes (fromIntegral (min 0 (max 255 i)) :: Word8) + toBytes v = (# 0#, 1#, js_wrapWord8Clamped v #) {-# INLINE toBytes #-} - fromBytes bs = fromIntegral (fromBytes bs :: Word8) + fromBytes (# off, _, arr #) = js_unwrapWord8Clamped arr off {-# INLINE fromBytes #-} byteSize _ = 1# {-# INLINE byteSize #-} @@ -86,6 +88,8 @@ {-# INLINE elementByteSize #-} ix _ (Clamped (I# x)) = x {-# INLINE ix #-} +foreign import javascript unsafe "h$wrapBuffer((new Uint8ClampedArray([$1])).buffer)" js_wrapWord8Clamped :: Word8Clamped -> ByteArray# +foreign import javascript unsafe "($1.uc || new Uint8ClampedArray($1.buf))[$2]" js_unwrapWord8Clamped :: ByteArray# -> Int# -> Word8Clamped instance ElementWise (Idx ('[] :: [Nat])) Word8Clamped Word8Clamped where indexOffset# x _ = x diff --git a/src-ghcjs/Numeric/Array/Family/ArrayT.hs b/src-ghcjs/Numeric/Array/Family/ArrayT.hs --- a/src-ghcjs/Numeric/Array/Family/ArrayT.hs +++ b/src-ghcjs/Numeric/Array/Family/ArrayT.hs @@ -15,7 +15,6 @@ {-# LANGUAGE JavaScriptFFI #-} {-# LANGUAGE GHCForeignImportPrim #-} {-# LANGUAGE UnliftedFFITypes #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE Strict #-} {-# OPTIONS_GHC -fno-warn-orphans #-} module Numeric.Array.Family.ArrayT () where @@ -35,7 +34,7 @@ import Numeric.Dimensions import Numeric.Dimensions.Traverse import Numeric.TypeLits -import Numeric.Matrix.Type +import Numeric.Matrix.Class type instance ElemRep (ArrayT t ds) = ElemRep t @@ -255,11 +254,15 @@ 1 -> GT 0 -> EQ _ -> LT + max = js_arrayMax + min = js_arrayMin foreign import javascript unsafe "$1.every(function (e, i) { return e < $2[i]; })" js_arrayTLT :: ArrayT t ds -> ArrayT t ds -> Bool foreign import javascript unsafe "$1.every(function (e, i) { return e <= $2[i]; })" js_arrayTLE :: ArrayT t ds -> ArrayT t ds -> Bool foreign import javascript unsafe "$1.every(function (e, i) { return e > $2[i]; })" js_arrayTGT :: ArrayT t ds -> ArrayT t ds -> Bool foreign import javascript unsafe "$1.every(function (e, i) { return e >= $2[i]; })" js_arrayTGE :: ArrayT t ds -> ArrayT t ds -> Bool foreign import javascript unsafe "$1.reduce(function (r, e, i) { return r === 0 ? (e > $2[i] ? 1 : (e < $2[i] ? -1 : 0)) : r;}, 0)" js_arrayTCmp :: ArrayT t ds -> ArrayT t ds -> Int +foreign import javascript unsafe "$1.map(function (e, i) { return Math.max(e,$2[i]); })" js_arrayMax :: ArrayT t ds -> ArrayT t ds -> ArrayT t ds +foreign import javascript unsafe "$1.map(function (e, i) { return Math.min(e,$2[i]); })" js_arrayMin :: ArrayT t ds -> ArrayT t ds -> ArrayT t ds instance Dimensions ds => Num (ArrayT Float ds) where diff --git a/src-ghcjs/Numeric/Array/Family/ArrayT.js b/src-ghcjs/Numeric/Array/Family/ArrayT.js --- a/src-ghcjs/Numeric/Array/Family/ArrayT.js +++ b/src-ghcjs/Numeric/Array/Family/ArrayT.js @@ -170,6 +170,7 @@ function h$easytensor_transpose(n, mat) { + 'use strict'; var nmat = new mat.constructor(mat.length), m = Math.round(mat.length / n); for(var i = 0; i < n; i++) { @@ -181,99 +182,119 @@ } function h$easytensor_eyeFloat32(n) { - var mat = new Float32Array(n*n).fill(0); + 'use strict'; + var mat = new Float32Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=1;} return mat; } function h$easytensor_eyeFloat64(n) { - var mat = new Float64Array(n*n).fill(0); + 'use strict'; + var mat = new Float64Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=1;} return mat; } function h$easytensor_eyeInt8(n) { - var mat = new Int8Array(n*n).fill(0); + 'use strict'; + var mat = new Int8Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=1;} return mat; } function h$easytensor_eyeInt16(n) { - var mat = new Int16Array(n*n).fill(0); + 'use strict'; + var mat = new Int16Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=1;} return mat; } function h$easytensor_eyeInt32(n) { - var mat = new Int32Array(n*n).fill(0); + 'use strict'; + var mat = new Int32Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=1;} return mat; } function h$easytensor_eyeUint8(n) { - var mat = new Uint8Array(n*n).fill(0); + 'use strict'; + var mat = new Uint8Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=1;} return mat; } function h$easytensor_eyeUint8Clamped(n) { - var mat = new Uint8ClampedArray(n*n).fill(0); + 'use strict'; + var mat = new Uint8ClampedArray(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=1;} return mat; } function h$easytensor_eyeUint16(n) { - var mat = new Uint16Array(n*n).fill(0); + 'use strict'; + var mat = new Uint16Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=1;} return mat; } function h$easytensor_eyeUint32(n) { - var mat = new Uint32Array(n*n).fill(0); + 'use strict'; + var mat = new Uint32Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=1;} return mat; } function h$easytensor_diagFloat32(n,x) { - var mat = new Float32Array(n*n).fill(0); + 'use strict'; + 'use strict'; + var mat = new Float32Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=x;} return mat; } function h$easytensor_diagFloat64(n,x) { - var mat = new Float64Array(n*n).fill(0); + 'use strict'; + var mat = new Float64Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=x;} return mat; } function h$easytensor_diagInt8(n,x) { - var mat = new Int8Array(n*n).fill(0); + 'use strict'; + var mat = new Int8Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=x;} return mat; } function h$easytensor_diagInt16(n,x) { - var mat = new Int16Array(n*n).fill(0); + 'use strict'; + var mat = new Int16Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=x;} return mat; } function h$easytensor_diagInt32(n,x) { - var mat = new Int32Array(n*n).fill(0); + 'use strict'; + var mat = new Int32Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=x;} return mat; } function h$easytensor_diagUint8(n,x) { - var mat = new Uint8Array(n*n).fill(0); + 'use strict'; + var mat = new Uint8Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=x;} return mat; } function h$easytensor_diagUint8Clamped(n,x) { - var mat = new Uint8ClampedArray(n*n).fill(0); + 'use strict'; + var mat = new Uint8ClampedArray(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=x;} return mat; } function h$easytensor_diagUint16(n,x) { - var mat = new Uint16Array(n*n).fill(0); + 'use strict'; + var mat = new Uint16Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=x;} return mat; } function h$easytensor_diagUint32(n,x) { - var mat = new Uint32Array(n*n).fill(0); + 'use strict'; + var mat = new Uint32Array(n*n); for(var i = 0; i < n*n; i += n + 1){mat[i]=x;} return mat; } function h$easytensor_trace(mat, n) { + 'use strict'; var r = 0; for(var i = 0; i < n*n; i += n + 1){r+=mat[i];} return r; @@ -281,6 +302,7 @@ function h$easytensor_det(mat, n) { + 'use strict'; switch (n) { case 1: return mat[0]; @@ -296,10 +318,12 @@ } function h$easytensor_detJSMat2(mat) { + 'use strict'; return (mat[0]*mat[3] - mat[1]*mat[2]); } function h$easytensor_detJSMat3(mat) { + 'use strict'; return ( mat[0]*(mat[4]*mat[8]-mat[5]*mat[7]) - mat[1]*(mat[3]*mat[8]-mat[5]*mat[6]) @@ -308,6 +332,7 @@ } function h$easytensor_detJSMat4(mat) { + 'use strict'; var n11 = mat[ 0 ], n12 = mat[ 4 ], n13 = mat[ 8 ], n14 = mat[ 12 ]; var n21 = mat[ 1 ], n22 = mat[ 5 ], n23 = mat[ 9 ], n24 = mat[ 13 ]; var n31 = mat[ 2 ], n32 = mat[ 6 ], n33 = mat[ 10 ], n34 = mat[ 14 ]; @@ -353,6 +378,7 @@ function h$easytensor_inverse(mat, n) { + 'use strict'; switch (n) { case 1: return 1 / mat[0]; @@ -369,6 +395,7 @@ function h$easytensor_inverseJSM4(mat) { + 'use strict'; var rez = new mat.constructor(16); rez[0] = mat[13]*(mat[ 6]*mat[11]-mat[10]*mat[ 7])+mat[ 9]*(mat[14]*mat[ 7]-mat[ 6]*mat[15])+mat[ 5]*(mat[10]*mat[15]-mat[14]*mat[11]); rez[4] = mat[12]*(mat[10]*mat[ 7]-mat[ 6]*mat[11])+mat[ 8]*(mat[ 6]*mat[15]-mat[14]*mat[ 7])+mat[ 4]*(mat[14]*mat[11]-mat[10]*mat[15]); @@ -390,12 +417,13 @@ if (det === 0) { return undefined; } else { - for(var i = 0; i < 16; i++) {rez[i] !== det;} + for(var i = 0; i < 16; i++) {rez[i] /= det;} return rez; } } function h$easytensor_inverseJSM3(mat) { + 'use strict'; var rez = new mat.constructor(9); rez[0] = mat[4]*mat[8] - mat[7]*mat[5]; rez[3] = mat[6]*mat[5] - mat[3]*mat[8]; @@ -410,12 +438,13 @@ if (det === 0) { return undefined; } else { - for(var i = 0; i < 9; i++) {rez[i] !== det;} + for(var i = 0; i < 9; i++) {rez[i] /= det;} return rez; } } function h$easytensor_inverseJSM2(mat) { + 'use strict'; var det = mat[0]*mat[3] - mat[1]*mat[2]; if (det === 0) { return undefined; @@ -430,6 +459,7 @@ function h$easytensor_contract(n,m,k,lhs,rhs) { + 'use strict'; var t, rez = new lhs.constructor(n*k); for(var i = 0; i < n; i++) { for(var j = 0; j < k; j++) { @@ -441,4 +471,18 @@ } } return rez; +} + + +function h$easytensor_dot(lhs, rhs) { + 'use strict'; + return lhs.reduce(function (r, e, i) { return r + e*rhs[i];}, 0); +} + +function h$easytensor_cross(a, b) { + 'use strict'; + return [ a[1]*b[2]-a[2]*b[1] + , a[2]*b[0]-a[0]*b[2] + , a[0]*b[1]-a[1]*b[0] + ]; } diff --git a/src-ghcjs/Numeric/DataFrame/Contraction.hs b/src-ghcjs/Numeric/DataFrame/Contraction.hs --- a/src-ghcjs/Numeric/DataFrame/Contraction.hs +++ b/src-ghcjs/Numeric/DataFrame/Contraction.hs @@ -11,6 +11,9 @@ {-# LANGUAGE UnboxedTuples #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE InstanceSigs #-} +{-# LANGUAGE UnliftedFFITypes #-} +{-# LANGUAGE JavaScriptFFI #-} +{-# LANGUAGE GHCForeignImportPrim #-} ----------------------------------------------------------------------------- -- | -- Module : Numeric.DataFrame.Contraction @@ -62,17 +65,29 @@ {-# INLINE (%*) #-} infixl 7 %* -instance ( ConcatList as bs asbs +instance {-# OVERLAPPABLE #-} + ( ConcatList as bs asbs , Dimensions as , Dimensions bs + , asbs ~ (a' ': sbs') ) => Contraction t as bs asbs where contract :: forall m . KnownDim m => DataFrame t (as +: m) -> DataFrame t (m :+ bs) -> DataFrame t asbs contract dx dy - | Refl <- unsafeCoerce Refl :: Array t asbs :~: ArrayT t asbs - , Refl <- unsafeCoerce Refl :: Array t (as +: m) :~: ArrayT t (as +: m) + | Refl <- unsafeCoerce Refl :: Array t (as +: m) :~: ArrayT t (as +: m) , Evidence <- inferConcatDimensions @as @bs - = KnownDataFrame $ js_conctract @t @as @m @bs (dimVal (dim @as)) (dimVal' @m) (dimVal (dim @bs)) (coerce dx) (coerce dy) + = KnownDataFrame $ js_contract @t @as @m @bs (dimVal (dim @as)) (dimVal' @m) (dimVal (dim @bs)) (coerce dx) (coerce dy) foreign import javascript unsafe "h$easytensor_contract($1,$2,$3,$4,$5)" - js_conctract :: forall t as m bs . Int -> Int -> Int -> ArrayT t (as +: m) -> ArrayT t (m :+ bs) -> ArrayT t (as ++ bs) + js_contract :: forall t as m bs . Int -> Int -> Int -> ArrayT t (as +: m) -> ArrayT t (m :+ bs) -> ArrayT t (as ++ bs) + +instance {-# OVERLAPPING #-} + Contraction t '[] '[] '[] where + contract :: forall m . KnownDim m => DataFrame t '[m] -> DataFrame t '[m] -> DataFrame t ('[] :: [Nat]) + contract dx dy + = KnownDataFrame $ unsafeCoerce (js_contract0 (coerce dx) (coerce dy)) + +foreign import javascript unsafe "$1.reduce(function (r, e, i) { return e*$2[i] + r;}, 0)" + js_contract0 :: ArrayT t '[m] -> ArrayT t '[m] -> Any + + diff --git a/src-ghcjs/Numeric/DataFrame/Mutable.hs b/src-ghcjs/Numeric/DataFrame/Mutable.hs new file mode 100644 --- /dev/null +++ b/src-ghcjs/Numeric/DataFrame/Mutable.hs @@ -0,0 +1,296 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UnboxedTuples #-} +{-# LANGUAGE JavaScriptFFI #-} +{-# LANGUAGE GHCForeignImportPrim #-} +{-# LANGUAGE UnliftedFFITypes #-} +{-# LANGUAGE TypeOperators #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.DataFrame.Mutable +-- Copyright : (c) Artem Chirkin +-- License : BSD3 +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- Interfrace to perform primitive stateful operations on mutable frames. +-- +----------------------------------------------------------------------------- + +module Numeric.DataFrame.Mutable + ( MutableFrame (..), MDataFrame (..) + , newDataFrame#, copyDataFrame#, copyMDataFrame#, unsafeFreezeDataFrame# + , freezeDataFrame#, thawDataFrame# + , writeDataFrame#, readDataFrame# + , newArrayBuffer#, arrayBuffer#, viewFloatArray#, viewDoubleArray# + , viewIntArray#, viewInt32Array#, viewInt16Array#, viewInt8Array# + , viewWordArray#, viewWord32Array#, viewWord16Array#, viewWord8Array#, viewWord8ClampedArray# + ) where + + + +import GHCJS.Types (IsJSVal(), JSVal) +import GHC.Int (Int16 (..), Int32 (..),Int8 (..)) +import GHC.Prim +import GHC.Types (Double (..), Float (..), Int (..), Word (..)) +import GHC.Word (Word16 (..), Word32 (..), Word8 (..)) +import Unsafe.Coerce (unsafeCoerce) + +import Numeric.DataFrame.Type +import Numeric.Array.Family +import Numeric.Dimensions + + +-- | Mutable DataFrame type +newtype MDataFrame s t (ns :: [Nat]) = MDataFrame (MutableArrayT s t ns) +instance IsJSVal (MDataFrame s t ds) + + +-- | Create a new mutable DataFrame. +newDataFrame# :: forall t (ns :: [Nat]) s + . (ElemTypeInference t, Dimensions ns) + => State# s -> (# State# s, MDataFrame s t ns #) +newDataFrame# = case elemTypeInstance @t of + ETFloat -> js_createFloatArray n + ETDouble -> js_createDoubleArray n + ETInt -> js_createIntArray n + ETInt8 -> js_createInt8Array n + ETInt16 -> js_createInt16Array n + ETInt32 -> js_createInt32Array n + ETWord -> js_createWordArray n + ETWord8 -> js_createWord8Array n + ETWord16 -> js_createWord16Array n + ETWord32 -> js_createWord32Array n + ETWord8C -> js_createWord8ClampedArray n + where + n = dimVal (dim @ns) +{-# INLINE newDataFrame# #-} + + +-- | Copy one DataFrame into another mutable DataFrame at specified position. +copyDataFrame# :: forall t (as :: [Nat]) (b' :: Nat) (b :: Nat) (bs :: [Nat]) (asbs :: [Nat]) s + . ( ArraySizeInference (as +: b') + , ConcatList as (b :+ bs) asbs + , Dimensions as + , Dimensions (b :+ bs) + ) + => DataFrame t (as +: b') -> Idx (b :+ bs) -> MDataFrame s t asbs -> State# s -> (# State# s, () #) +copyDataFrame# df i mdf s0 = case arraySizeInstance @(as +: b') of + ASScalar -> df `seq` (# js_writeArrayOffsetJSVal# mdf (fromEnum i) (unsafeCoerce df) s0, () #) + ASArray -> js_copyDataFrame (coerce df) (fromEnum i * dimVal (dim @as)) mdf s0 +{-# INLINE copyDataFrame# #-} + + + +-- | Copy one mutable DataFrame into another mutable DataFrame at specified position. +copyMDataFrame# :: forall t (as :: [Nat]) (b' :: Nat) (b :: Nat) (bs :: [Nat]) (asbs :: [Nat]) s + . ( ConcatList as (b :+ bs) asbs + , Dimensions as + , Dimensions (b :+ bs) + ) + => MDataFrame s t (as +: b') -> Idx (b :+ bs) -> MDataFrame s t asbs -> State# s -> (# State# s, () #) +copyMDataFrame# d i = js_copyMDataFrame d (fromEnum i * dimVal (dim @as)) +{-# INLINE copyMDataFrame# #-} + + +-- | Make a mutable DataFrame immutable, without copying. +unsafeFreezeDataFrame# :: forall t (ns :: [Nat]) s + . (MutableFrame t ns, ArraySizeInference ns) + => MDataFrame s t ns -> State# s -> (# State# s, DataFrame t ns #) +unsafeFreezeDataFrame# a s = case arraySizeInstance @ns of + ASScalar -> case readDataFrameOff# a 0# s of + (# s1, v #) -> (# s1, coerce v #) + ASArray -> (# s, coerce a #) +{-# INLINE unsafeFreezeDataFrame# #-} + + +--unsafeThawArrayT# :: ArrayT t ds -> State# s -> (#State# s, MutableArrayT s t ds #) +--unsafeThawArrayT# a s = (# s, coerce a #) +--{-# INLINE unsafeThawArrayT# #-} + + +-- | Copy content of a mutable DataFrame into a new immutable DataFrame. +freezeDataFrame# :: forall t (ns :: [Nat]) s + . (MutableFrame t ns, ArraySizeInference ns) + => MDataFrame s t ns -> State# s -> (# State# s, DataFrame t ns #) +freezeDataFrame# a s = case arraySizeInstance @ns of + ASScalar -> case readDataFrameOff# a 0# s of + (# s1, v #) -> (# s1, coerce v #) + ASArray -> case js_freeze a s of + (# s1, v #) -> (# s1, coerce v #) +{-# INLINE freezeDataFrame# #-} + + + +-- | Create a new mutable DataFrame and copy content of immutable one in there. +thawDataFrame# :: forall t (ns :: [Nat]) s + . (MutableFrame t ns, ArrayInstanceInference t ns) + => DataFrame t ns -> State# s -> (# State# s, MDataFrame s t ns #) +thawDataFrame# a s = case arraySizeInstance @ns of + ASScalar -> case newDataFrame# @t @'[] s of + (# s1, df #) -> (# writeDataFrameOff# df 0# (coerce a) s1, df #) + ASArray -> js_thaw (coerce a) s +{-# INLINE thawDataFrame# #-} + +-- | Write a single element at the specified index +writeDataFrame# :: forall t (ns :: [Nat]) s + . ( MutableFrame t ns, Dimensions ns ) + => MDataFrame s t ns -> Idx ns -> t -> State# s -> (# State# s, () #) +writeDataFrame# mdf ei x s | I# i <- fromEnum ei = (# writeDataFrameOff# mdf i x s, () #) +{-# INLINE writeDataFrame# #-} + +-- | Read a single element at the specified index +readDataFrame# :: forall t (ns :: [Nat]) s + . ( MutableFrame t ns, Dimensions ns ) + => MDataFrame s t ns -> Idx ns -> State# s -> (# State# s, t #) +readDataFrame# mdf ei | I# i <- fromEnum ei = readDataFrameOff# mdf i +{-# INLINE readDataFrame# #-} + +class MutableFrame t (ns :: [Nat]) where + -- | Write a single element at the specified element offset + writeDataFrameOff# :: MDataFrame s t ns -> Int# -> t -> State# s -> State# s + -- | Read a single element at the specified element offset + readDataFrameOff# :: MDataFrame s t ns -> Int# -> State# s -> (# State# s, t #) + +instance MutableFrame Float (ns :: [Nat]) where + writeDataFrameOff# = js_writeArrayOffsetFloat# + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# = js_readArrayOffsetFloat# + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Double (ns :: [Nat]) where + writeDataFrameOff# = js_writeArrayOffsetDouble# + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# = js_readArrayOffsetDouble# + {-# INLINE readDataFrameOff# #-} + + +instance MutableFrame Int (ns :: [Nat]) where + writeDataFrameOff# = js_writeArrayOffsetInt# + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# = js_readArrayOffsetInt# + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Int8 (ns :: [Nat]) where + writeDataFrameOff# = js_writeArrayOffsetInt8# + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# = js_readArrayOffsetInt8# + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Int16 (ns :: [Nat]) where + writeDataFrameOff# = js_writeArrayOffsetInt16# + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# = js_readArrayOffsetInt16# + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Int32 (ns :: [Nat]) where + writeDataFrameOff# = js_writeArrayOffsetInt32# + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# = js_readArrayOffsetInt32# + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Word (ns :: [Nat]) where + writeDataFrameOff# = js_writeArrayOffsetWord# + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# = js_readArrayOffsetWord# + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Word8 (ns :: [Nat]) where + writeDataFrameOff# = js_writeArrayOffsetWord8# + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# = js_readArrayOffsetWord8# + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Word16 (ns :: [Nat]) where + writeDataFrameOff# = js_writeArrayOffsetWord16# + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# = js_readArrayOffsetWord16# + {-# INLINE readDataFrameOff# #-} + +instance MutableFrame Word32 (ns :: [Nat]) where + writeDataFrameOff# = js_writeArrayOffsetWord32# + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# = js_readArrayOffsetWord32# + {-# INLINE readDataFrameOff# #-} + + +instance MutableFrame Word8Clamped (ns :: [Nat]) where + writeDataFrameOff# = js_writeArrayOffsetWord8Clamped# + {-# INLINE writeDataFrameOff# #-} + readDataFrameOff# m o s = case js_readArrayOffsetWord8Clamped# m o s of + (# s1, y #) -> (# s1, coerce y #) + {-# INLINE readDataFrameOff# #-} + + + +foreign import javascript unsafe "new Float32Array($1)" js_createFloatArray :: Int -> State# s -> (# State# s, MDataFrame s Float ds #) +foreign import javascript unsafe "new Float64Array($1)" js_createDoubleArray :: Int -> State# s -> (# State# s, MDataFrame s Double ds #) +foreign import javascript unsafe "new Int32Array($1)" js_createIntArray :: Int -> State# s -> (# State# s, MDataFrame s Int ds #) +foreign import javascript unsafe "new Int32Array($1)" js_createInt32Array :: Int -> State# s -> (# State# s, MDataFrame s Int32 ds #) +foreign import javascript unsafe "new Int16Array($1)" js_createInt16Array :: Int -> State# s -> (# State# s, MDataFrame s Int16 ds #) +foreign import javascript unsafe "new Int8Array($1)" js_createInt8Array :: Int -> State# s -> (# State# s, MDataFrame s Int8 ds #) +foreign import javascript unsafe "new Uint32Array($1)" js_createWordArray :: Int -> State# s -> (# State# s, MDataFrame s Word ds #) +foreign import javascript unsafe "new Uint32Array($1)" js_createWord32Array :: Int -> State# s -> (# State# s, MDataFrame s Word32 ds #) +foreign import javascript unsafe "new Uint16Array($1)" js_createWord16Array :: Int -> State# s -> (# State# s, MDataFrame s Word16 ds #) +foreign import javascript unsafe "new Uint8Array($1)" js_createWord8Array :: Int -> State# s -> (# State# s, MDataFrame s Word8 ds #) +foreign import javascript unsafe "new Uint8ClampedArray($1)" js_createWord8ClampedArray :: Int -> State# s -> (# State# s, MDataFrame s Word8Clamped ds #) + + +foreign import javascript unsafe "$1[$2]" js_readArrayOffsetFloat# :: MDataFrame s Float ds -> Int# -> State# s -> (# State# s, Float #) +foreign import javascript unsafe "$1[$2]" js_readArrayOffsetDouble# :: MDataFrame s Double ds -> Int# -> State# s -> (# State# s, Double #) +foreign import javascript unsafe "$1[$2]" js_readArrayOffsetInt# :: MDataFrame s Int ds -> Int# -> State# s -> (# State# s, Int #) +foreign import javascript unsafe "$1[$2]" js_readArrayOffsetInt8# :: MDataFrame s Int8 ds -> Int# -> State# s -> (# State# s, Int8 #) +foreign import javascript unsafe "$1[$2]" js_readArrayOffsetInt16# :: MDataFrame s Int16 ds -> Int# -> State# s -> (# State# s, Int16 #) +foreign import javascript unsafe "$1[$2]" js_readArrayOffsetInt32# :: MDataFrame s Int32 ds -> Int# -> State# s -> (# State# s, Int32 #) +foreign import javascript unsafe "$1[$2]" js_readArrayOffsetWord# :: MDataFrame s Word ds -> Int# -> State# s -> (# State# s, Word #) +foreign import javascript unsafe "$1[$2]" js_readArrayOffsetWord8# :: MDataFrame s Word8 ds -> Int# -> State# s -> (# State# s, Word8 #) +foreign import javascript unsafe "$1[$2]" js_readArrayOffsetWord8Clamped# :: MDataFrame s Word8Clamped ds -> Int# -> State# s -> (# State# s, Int #) +foreign import javascript unsafe "$1[$2]" js_readArrayOffsetWord16# :: MDataFrame s Word16 ds -> Int# -> State# s -> (# State# s, Word16 #) +foreign import javascript unsafe "$1[$2]" js_readArrayOffsetWord32# :: MDataFrame s Word32 ds -> Int# -> State# s -> (# State# s, Word32 #) + + + +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetFloat# :: MDataFrame s Float ds -> Int# -> Float -> State# s -> State# s +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetDouble# :: MDataFrame s Double ds -> Int# -> Double -> State# s -> State# s +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetInt# :: MDataFrame s Int ds -> Int# -> Int -> State# s -> State# s +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetInt8# :: MDataFrame s Int8 ds -> Int# -> Int8 -> State# s -> State# s +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetInt16# :: MDataFrame s Int16 ds -> Int# -> Int16 -> State# s -> State# s +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetInt32# :: MDataFrame s Int32 ds -> Int# -> Int32 -> State# s -> State# s +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetWord# :: MDataFrame s Word ds -> Int# -> Word -> State# s -> State# s +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetWord8# :: MDataFrame s Word8 ds -> Int# -> Word8 -> State# s -> State# s +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetWord8Clamped# :: MDataFrame s Word8Clamped ds -> Int# -> Word8Clamped -> State# s -> State# s +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetWord16# :: MDataFrame s Word16 ds -> Int# -> Word16 -> State# s -> State# s +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetWord32# :: MDataFrame s Word32 ds -> Int# -> Word32 -> State# s -> State# s +foreign import javascript unsafe "$1[$2] = $3;" js_writeArrayOffsetJSVal# :: MDataFrame s t ds -> Int -> JSVal -> State# s -> State# s + + +foreign import javascript unsafe "$3.set($1, $2);" js_copyDataFrame :: ArrayT t as -> Int -> MDataFrame s t asbs -> State# s -> (# State# s, () #) +foreign import javascript unsafe "$3.set($1, $2);" js_copyMDataFrame :: MDataFrame s t as -> Int -> MDataFrame s t asbs -> State# s -> (# State# s, () #) + + +foreign import javascript unsafe "$1.slice()" js_freeze :: MDataFrame s t as -> State# s -> (# State# s, ArrayT t ds #) +foreign import javascript unsafe "$1.slice()" js_thaw :: ArrayT t as -> State# s -> (# State# s, MDataFrame s t ds #) + + + +foreign import javascript unsafe "new ArrayBuffer($1)" newArrayBuffer# :: Int -> State# s -> (# State# s, JSVal #) +foreign import javascript unsafe "new Float32Array($1)" viewFloatArray# :: JSVal -> State# s -> (# State# s, MDataFrame s Float ds #) +foreign import javascript unsafe "new Float64Array($1)" viewDoubleArray# :: JSVal -> State# s -> (# State# s, MDataFrame s Double ds #) +foreign import javascript unsafe "new Int32Array($1)" viewIntArray# :: JSVal -> State# s -> (# State# s, MDataFrame s Int ds #) +foreign import javascript unsafe "new Int32Array($1)" viewInt32Array# :: JSVal -> State# s -> (# State# s, MDataFrame s Int32 ds #) +foreign import javascript unsafe "new Int16Array($1)" viewInt16Array# :: JSVal -> State# s -> (# State# s, MDataFrame s Int16 ds #) +foreign import javascript unsafe "new Int8Array($1)" viewInt8Array# :: JSVal -> State# s -> (# State# s, MDataFrame s Int8 ds #) +foreign import javascript unsafe "new Uint32Array($1)" viewWordArray# :: JSVal -> State# s -> (# State# s, MDataFrame s Word ds #) +foreign import javascript unsafe "new Uint32Array($1)" viewWord32Array# :: JSVal -> State# s -> (# State# s, MDataFrame s Word32 ds #) +foreign import javascript unsafe "new Uint16Array($1)" viewWord16Array# :: JSVal -> State# s -> (# State# s, MDataFrame s Word16 ds #) +foreign import javascript unsafe "new Uint8Array($1)" viewWord8Array# :: JSVal -> State# s -> (# State# s, MDataFrame s Word8 ds #) +foreign import javascript unsafe "new Uint8ClampedArray($1)" viewWord8ClampedArray# :: JSVal -> State# s -> (# State# s, MDataFrame s Word8Clamped ds #) +foreign import javascript unsafe "$1.buffer" arrayBuffer# :: MDataFrame s t ds -> State# s -> (# State# s, JSVal #) diff --git a/src-ghcjs/Numeric/Matrix/Mat44.js b/src-ghcjs/Numeric/Matrix/Mat44.js new file mode 100644 --- /dev/null +++ b/src-ghcjs/Numeric/Matrix/Mat44.js @@ -0,0 +1,103 @@ +function h$easytensor_m4fromHom(v) { + 'use strict'; + var r = v.slice(0,3), t = v[3]; + if (t !== 0) { + r[0] /= t; r[1] /= t; r[2] /= t; + } + return r; +} + + +function h$easytensor_m4translate(v) { + 'use strict'; + var m = new v.constructor(16); + m.set(v, 12); + m[0] = 1; + m[5] = 1; + m[10] = 1; + m[15] = 1; + return m; +} + +function h$easytensor_m4rotateX(a) { + 'use strict'; + var c = Math.cos(a), s = Math.sin(a); + return [ 1, 0, 0, 0 + , 0, c, s, 0 + , 0,-s, c, 0 + , 0, 0, 0, 1]; +} + +function h$easytensor_m4rotateY(a) { + 'use strict'; + var c = Math.cos(a), s = Math.sin(a); + return [ c, 0,-s, 0 + , 0, 1, 0, 0 + , s, 0, c, 0 + , 0, 0, 0, 1]; +} + +function h$easytensor_m4rotateZ(a) { + 'use strict'; + var c = Math.cos(a), s = Math.sin(a); + return [ c, s, 0, 0 + ,-s, c, 0, 0 + , 0, 0, 1, 0 + , 0, 0, 0, 1]; +} + +function h$easytensor_m4rotate(vec, a) { + 'use strict'; + var c = Math.cos(a); + var s = Math.sin(a); + var c1 = 1 - c; + var x = vec[0], y = vec[1], z = vec[2]; + return [ c + c1*x*x, c1*x*y + s*z, c1*x*z - s*y, 0 + , c1*x*y - s*z, c + c1*y*y, c1*y*z + s*x, 0 + , c1*x*z + s*y, c1*y*z - s*x, c + c1*z*z, 0 + , 0, 0, 0, 1]; +} + +function h$easytensor_m4rotateEuler(x, y, z) { + 'use strict'; + var cx = Math.cos(x), sx = Math.sin(x), cy = Math.cos(y), sy = Math.sin(y), cz = Math.cos(z), sz = Math.sin(z); + return [ 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]; +} + +function h$easytensor_m4lookAt(up,camera,point) { + 'use strict'; + var zDir = [ camera[0] - point[0], camera[1] - point[1], camera[2] - point[2] ]; + var t = Math.hypot.apply(null,zDir); + zDir = zDir.map(function (e){return e / t;}); + var xDir = h$easytensor_cross(up,zDir); + t = Math.hypot.apply(null,xDir); + xDir = xDir.map(function (e){return e / t;}); + var yDir = h$easytensor_cross(zDir,xDir); + return [ xDir[0], yDir[0], zDir[0], 0 + , xDir[1], yDir[1], zDir[1], 0 + , xDir[2], yDir[2], zDir[2], 0 + , - h$easytensor_dot(xDir,camera), - h$easytensor_dot(yDir,camera), - h$easytensor_dot(zDir,camera), 1 + ]; +} + +function h$easytensor_m4perspective(n, f, fovy, aspect) { + 'use strict'; + var h2 = n*Math.tan(fovy/2); + var w2 = aspect*h2; + return [ n/w2, 0, 0, 0 + , 0, n/h2, 0, 0 + , 0, 0, (n+f)/(n-f),-1 + , 0, 0, 2*n*f/(n-f), 0 ]; +} + +function h$easytensor_m4orthogonal(n, f, w, h) { + 'use strict'; + return [ 2/w, 0, 0, 0 + , 0, 2/h, 0, 0 + , 0, 0, 2/(n-f), 0 + , 0, 0, (n+f)/(n-f), 1 ]; +} + diff --git a/src-ghcjs/Numeric/Matrix/Mat44d.hs b/src-ghcjs/Numeric/Matrix/Mat44d.hs new file mode 100644 --- /dev/null +++ b/src-ghcjs/Numeric/Matrix/Mat44d.hs @@ -0,0 +1,69 @@ +{-# LANGUAGE JavaScriptFFI #-} +{-# LANGUAGE DataKinds #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +module Numeric.Matrix.Mat44d () where + + +import Numeric.DataFrame.Type +import Numeric.Vector +import Numeric.Matrix.Class + + +instance HomTransform4 Double where + translate4 = js_translate + {-# INLINE translate4 #-} + translate3 = js_translate + {-# INLINE translate3 #-} + rotateX = js_rotateX + {-# INLINE rotateX #-} + rotateY = js_rotateY + {-# INLINE rotateY #-} + rotateZ = js_rotateZ + {-# INLINE rotateZ #-} + rotate = js_rotate + {-# INLINE rotate #-} + rotateEuler = js_rotateEuler + {-# INLINE rotateEuler #-} + lookAt = js_lookAt + {-# INLINE lookAt #-} + perspective = js_perspective + {-# INLINE perspective #-} + orthogonal = js_orthogonal + {-# INLINE orthogonal #-} + toHomPoint = js_toHomPoint + {-# INLINE toHomPoint #-} + toHomVector = js_toHomVector + {-# INLINE toHomVector #-} + fromHom = js_fromHom + {-# INLINE fromHom #-} + + + +foreign import javascript unsafe "h$easytensor_m4translate($1)" + js_translate :: Vector Double n -> Matrix Double 4 4 +foreign import javascript unsafe "new Float64Array(h$easytensor_m4rotateX($1))" + js_rotateX :: Double -> Matrix Double 4 4 +foreign import javascript unsafe "new Float64Array(h$easytensor_m4rotateY($1))" + js_rotateY :: Double -> Matrix Double 4 4 +foreign import javascript unsafe "new Float64Array(h$easytensor_m4rotateZ($1))" + js_rotateZ :: Double -> Matrix Double 4 4 +foreign import javascript unsafe "new Float64Array(h$easytensor_m4rotate($1, $2))" + js_rotate :: Vector Double 3 -> Double -> Matrix Double 4 4 +foreign import javascript unsafe "new Float64Array(h$easytensor_m4rotateEuler($1, $2, $3))" + js_rotateEuler :: Double -> Double -> Double -> Matrix Double 4 4 +foreign import javascript unsafe "new Float64Array(h$easytensor_m4lookAt($1,$2,$3))" + js_lookAt :: Vector Double 3 -> Vector Double 3 -> Vector Double 3 -> Matrix Double 4 4 +foreign import javascript unsafe "new Float64Array(h$easytensor_m4perspective($1, $2, $3, $4))" + js_perspective :: Double -> Double -> Double -> Double -> Matrix Double 4 4 +foreign import javascript unsafe "new Float64Array(h$easytensor_m4orthogonal($1, $2, $3, $4))" + js_orthogonal :: Double -> Double -> Double -> Double -> Matrix Double 4 4 + +foreign import javascript unsafe "$r = new $1.constructor(4); $r.set($1); $r[3] = 1;" + js_toHomPoint :: Vector Double 3 -> Vector Double 4 +foreign import javascript unsafe "$r = new $1.constructor(4); $r.set($1); $r[3] = 0;" + js_toHomVector :: Vector Double 3 -> Vector Double 4 +foreign import javascript unsafe "h$easytensor_m4fromHom($1)" + js_fromHom :: Vector Double 4 -> Vector Double 3 + + + diff --git a/src-ghcjs/Numeric/Matrix/Mat44f.hs b/src-ghcjs/Numeric/Matrix/Mat44f.hs new file mode 100644 --- /dev/null +++ b/src-ghcjs/Numeric/Matrix/Mat44f.hs @@ -0,0 +1,69 @@ +{-# LANGUAGE JavaScriptFFI #-} +{-# LANGUAGE DataKinds #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +module Numeric.Matrix.Mat44f () where + + +import Numeric.DataFrame.Type +import Numeric.Vector +import Numeric.Matrix.Class + + +instance HomTransform4 Float where + translate4 = js_translate + {-# INLINE translate4 #-} + translate3 = js_translate + {-# INLINE translate3 #-} + rotateX = js_rotateX + {-# INLINE rotateX #-} + rotateY = js_rotateY + {-# INLINE rotateY #-} + rotateZ = js_rotateZ + {-# INLINE rotateZ #-} + rotate = js_rotate + {-# INLINE rotate #-} + rotateEuler = js_rotateEuler + {-# INLINE rotateEuler #-} + lookAt = js_lookAt + {-# INLINE lookAt #-} + perspective = js_perspective + {-# INLINE perspective #-} + orthogonal = js_orthogonal + {-# INLINE orthogonal #-} + toHomPoint = js_toHomPoint + {-# INLINE toHomPoint #-} + toHomVector = js_toHomVector + {-# INLINE toHomVector #-} + fromHom = js_fromHom + {-# INLINE fromHom #-} + + + +foreign import javascript unsafe "h$easytensor_m4translate($1)" + js_translate :: Vector Float n -> Matrix Float 4 4 +foreign import javascript unsafe "new Float32Array(h$easytensor_m4rotateX($1))" + js_rotateX :: Float -> Matrix Float 4 4 +foreign import javascript unsafe "new Float32Array(h$easytensor_m4rotateY($1))" + js_rotateY :: Float -> Matrix Float 4 4 +foreign import javascript unsafe "new Float32Array(h$easytensor_m4rotateZ($1))" + js_rotateZ :: Float -> Matrix Float 4 4 +foreign import javascript unsafe "new Float32Array(h$easytensor_m4rotate($1, $2))" + js_rotate :: Vector Float 3 -> Float -> Matrix Float 4 4 +foreign import javascript unsafe "new Float32Array(h$easytensor_m4rotateEuler($1, $2, $3))" + js_rotateEuler :: Float -> Float -> Float -> Matrix Float 4 4 +foreign import javascript unsafe "new Float32Array(h$easytensor_m4lookAt($1,$2,$3))" + js_lookAt :: Vector Float 3 -> Vector Float 3 -> Vector Float 3 -> Matrix Float 4 4 +foreign import javascript unsafe "new Float32Array(h$easytensor_m4perspective($1, $2, $3, $4))" + js_perspective :: Float -> Float -> Float -> Float -> Matrix Float 4 4 +foreign import javascript unsafe "new Float32Array(h$easytensor_m4orthogonal($1, $2, $3, $4))" + js_orthogonal :: Float -> Float -> Float -> Float -> Matrix Float 4 4 + +foreign import javascript unsafe "$r = new $1.constructor(4); $r.set($1); $r[3] = 1;" + js_toHomPoint :: Vector Float 3 -> Vector Float 4 +foreign import javascript unsafe "$r = new $1.constructor(4); $r.set($1); $r[3] = 0;" + js_toHomVector :: Vector Float 3 -> Vector Float 4 +foreign import javascript unsafe "h$easytensor_m4fromHom($1)" + js_fromHom :: Vector Float 4 -> Vector Float 3 + + + diff --git a/src-ghcjs/Numeric/Quaternion/QDouble.hs b/src-ghcjs/Numeric/Quaternion/QDouble.hs new file mode 100644 --- /dev/null +++ b/src-ghcjs/Numeric/Quaternion/QDouble.hs @@ -0,0 +1,285 @@ +{-# LANGUAGE JavaScriptFFI #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE FlexibleInstances #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +module Numeric.Quaternion.QDouble + ( QDouble, Quater (..) + ) where + +import Data.JSString (JSString, unpack') +import Data.Coerce (coerce) + +import Numeric.Array +import Numeric.DataFrame.Type +import Numeric.Vector +import Numeric.Matrix + +import Numeric.Quaternion.Class + + +type QDouble = Quater Double + +instance Quaternion Double where + newtype Quater Double = QDouble (ArrayT Double '[4]) + {-# INLINE packQ #-} + packQ = js_packQ + {-# INLINE fromVecNum #-} + fromVecNum = js_fromVecNum + {-# INLINE fromVec4 #-} + fromVec4 = coerce + {-# INLINE toVec4 #-} + toVec4 = coerce + {-# INLINE unpackQ #-} + unpackQ = js_unpackQ + {-# INLINE square #-} + square = js_square + {-# INLINE im #-} + im = js_im + {-# INLINE re #-} + re = js_re + {-# INLINE imVec #-} + imVec = js_imVec + {-# INLINE taker #-} + taker = js_taker + {-# INLINE takei #-} + takei = js_takei + {-# INLINE takej #-} + takej = js_takej + {-# INLINE takek #-} + takek = js_takek + {-# INLINE conjugate #-} + conjugate = js_conjugate + {-# INLINE rotScale #-} + rotScale = js_rotScale + {-# INLINE getRotScale #-} + getRotScale = js_getRotScale + {-# INLINE axisRotation #-} + axisRotation = js_axisRotation + {-# INLINE qArg #-} + qArg = js_qArg + {-# INLINE fromMatrix33 #-} + fromMatrix33 = js_fromMatrix33 + {-# INLINE fromMatrix44 #-} + fromMatrix44 = js_fromMatrix44 + {-# INLINE toMatrix33 #-} + toMatrix33 = js_toMatrix33 + {-# INLINE toMatrix44 #-} + toMatrix44 = js_toMatrix44 + + +foreign import javascript unsafe "new Float64Array([$1,$2,$3,$4])" + js_packQ :: Double -> Double -> Double -> Double -> QDouble + +foreign import javascript unsafe "new Float64Array([$1[0],$1[1],$1[2],$2])" + js_fromVecNum :: Vector Double 3 -> Double -> QDouble + +foreign import javascript unsafe "$r1 = $1[0];$r2 = $1[1];$r3 = $1[2];$r4 = $1[3];" + js_unpackQ :: QDouble -> (Double,Double,Double,Double) + + +foreign import javascript unsafe "$1[0]*$1[0] + $1[1]*$1[1] + $1[2]*$1[2] + $1[3]*$1[3]" + js_square :: QDouble -> Double + +foreign import javascript unsafe "new Float64Array([$1[0],$1[1],$1[2],0])" + js_im :: QDouble -> QDouble + +foreign import javascript unsafe "$1.slice(0,3)" + js_imVec :: QDouble -> Vector Double 3 + +foreign import javascript unsafe "$r = new Float64Array(4); $r[3] = $1[3];" + js_re :: QDouble -> QDouble + +foreign import javascript unsafe "$1[3]" + js_taker :: QDouble -> Double + +foreign import javascript unsafe "$1[0]" + js_takei :: QDouble -> Double + +foreign import javascript unsafe "$1[1]" + js_takej :: QDouble -> Double + +foreign import javascript unsafe "$1[2]" + js_takek :: QDouble -> Double + +foreign import javascript unsafe "new Float64Array([-$1[0],-$1[1],-$1[2],$1[3]])" + js_conjugate :: QDouble -> QDouble + + +foreign import javascript unsafe "new Float64Array(h$easytensor_rotScale($1,$2))" + js_rotScale :: QDouble -> Vector Double 3 -> Vector Double 3 + +foreign import javascript unsafe "new Float64Array(h$easytensor_getRotScale($1,$2))" + js_getRotScale :: Vector Double 3 -> Vector Double 3 -> QDouble + +foreign import javascript unsafe "new Float64Array(h$easytensor_axisRotation($1,$2))" + js_axisRotation :: Vector Double 3 -> Double -> QDouble + +foreign import javascript unsafe "h$easytensor_qArg($1)" + js_qArg :: QDouble -> Double + +foreign import javascript unsafe "new Float64Array(h$easytensor_qfromMatrix33($1))" + js_fromMatrix33 :: Matrix Double 3 3 -> QDouble + +foreign import javascript unsafe "new Float64Array(h$easytensor_qfromMatrix44($1))" + js_fromMatrix44 :: Matrix Double 4 4 -> QDouble + +foreign import javascript unsafe "new Float64Array(h$easytensor_qtoMatrix33($1))" + js_toMatrix33 :: QDouble -> Matrix Double 3 3 + +foreign import javascript unsafe "new Float64Array(h$easytensor_qtoMatrix44($1))" + js_toMatrix44 :: QDouble -> Matrix Double 4 4 + + +-------------------------------------------------------------------------- +-- Num +-------------------------------------------------------------------------- + +instance Num QDouble where + {-# INLINE (+) #-} + (+) = js_plus + {-# INLINE (-) #-} + (-) = js_minus + {-# INLINE (*) #-} + (*) = js_times + {-# INLINE abs #-} + abs = js_abs + {-# INLINE signum #-} + signum = js_signum + {-# INLINE negate #-} + negate = fromVec4 . negate . toVec4 + {-# INLINE fromInteger #-} + fromInteger = js_toQuaternion . fromInteger + +foreign import javascript unsafe "$1.map(function (e, i) {return e + $2[i];})" + js_plus :: QDouble -> QDouble -> QDouble + + +foreign import javascript unsafe "$1.map(function (e, i) {return e - $2[i];})" + js_minus :: QDouble -> QDouble -> QDouble + +foreign import javascript unsafe "new Float64Array(\ + \[ $1[3]*$2[0] + $1[0]*$2[3] + $1[1]*$2[2] - $1[2]*$2[1]\ + \, $1[3]*$2[1] - $1[0]*$2[2] + $1[1]*$2[3] + $1[2]*$2[0]\ + \, $1[3]*$2[2] + $1[0]*$2[1] - $1[1]*$2[0] + $1[2]*$2[3]\ + \, $1[3]*$2[3] - $1[0]*$2[0] - $1[1]*$2[1] - $1[2]*$2[2] ])" + js_times :: QDouble -> QDouble -> QDouble + +foreign import javascript unsafe "new Float64Array([0,0,0,Math.hypot($1[0],$1[1],$1[2],$1[3])])" + js_abs :: QDouble -> QDouble + +foreign import javascript unsafe "var l = Math.hypot($1[0],$1[1],$1[2],$1[3]); $r = (l == 0) ? (new Float64Array(4)) : $1.map(function (e) {return e/l;})" + js_signum :: QDouble -> QDouble + +foreign import javascript unsafe "$r = new Float64Array(4); $r[3] = $1;" + js_toQuaternion :: Double -> QDouble + +{-# RULES +"realToFrac/DoubleQDouble" realToFrac = js_toQuaternion +"realToFrac/aQDouble" realToFrac = js_toQuaternion . realToFrac +"fromIntegral/aQDouble" fromIntegral = js_toQuaternion . fromIntegral + #-} + + +-------------------------------------------------------------------------- +-- Fractional +-------------------------------------------------------------------------- + +instance Fractional QDouble where + {-# INLINE recip #-} + recip = js_recip + {-# INLINE (/) #-} + (/) p = js_times p . js_recip + {-# INLINE fromRational #-} + fromRational = js_toQuaternion . fromRational + +foreign import javascript unsafe "new Float64Array(h$easytensor_qrecip($1))" + js_recip :: QDouble -> QDouble + + +-------------------------------------------------------------------------- +-- Floating +-------------------------------------------------------------------------- + +instance Floating QDouble where + {-# INLINE pi #-} + pi = js_toQuaternion pi + {-# INLINE exp #-} + exp = js_exp + {-# INLINE log #-} + log = js_log + {-# INLINE sqrt #-} + sqrt = js_sqrt + {-# INLINE sin #-} + sin = js_sin + {-# INLINE cos #-} + cos = js_cos + {-# INLINE tan #-} + tan = js_tan + {-# INLINE sinh #-} + sinh = js_sinh + {-# INLINE cosh #-} + cosh = js_cosh + {-# INLINE tanh #-} + tanh = js_tanh + {-# INLINE asin #-} + asin q = -i * log (i*q + sqrt (1 - q*q)) + where i = signum . im $ q + {-# INLINE acos #-} + acos q = pi/2 - asin q + {-# INLINE atan #-} + atan q = if square imq == 0 + then js_toQuaternion (atan $ taker q) + else i / 2 * log ( (i + q) / (i - q) ) + where i = signum imq + imq = im q + {-# INLINE asinh #-} + asinh q = log (q + sqrt (q*q + 1)) + {-# INLINE acosh #-} + acosh q = log (q + sqrt (q*q - 1)) + {-# INLINE atanh #-} + atanh q = 0.5 * log ((1+q)/(1-q)) + +foreign import javascript unsafe "new Float64Array(h$easytensor_qexp($1))" js_exp :: QDouble -> QDouble +foreign import javascript unsafe "new Float64Array(h$easytensor_qlog($1))" js_log :: QDouble -> QDouble +foreign import javascript unsafe "new Float64Array(h$easytensor_qsqrt($1))" js_sqrt :: QDouble -> QDouble +foreign import javascript unsafe "new Float64Array(h$easytensor_qsin($1))" js_sin :: QDouble -> QDouble +foreign import javascript unsafe "new Float64Array(h$easytensor_qcos($1))" js_cos :: QDouble -> QDouble +foreign import javascript unsafe "new Float64Array(h$easytensor_qtan($1))" js_tan :: QDouble -> QDouble +foreign import javascript unsafe "new Float64Array(h$easytensor_qsinh($1))" js_sinh :: QDouble -> QDouble +foreign import javascript unsafe "new Float64Array(h$easytensor_qcosh($1))" js_cosh :: QDouble -> QDouble +foreign import javascript unsafe "new Float64Array(h$easytensor_qtanh($1))" js_tanh :: QDouble -> QDouble + + +-------------------------------------------------------------------------- +-- Eq +-------------------------------------------------------------------------- + +instance Eq QDouble where + {-# INLINE (==) #-} + (==) = js_eq + {-# INLINE (/=) #-} + (/=) = js_neq + + + +foreign import javascript unsafe "$1[0] === $2[0] && $1[1] === $2[1] && $1[2] === $2[2] && $1[3] === $2[3]" + js_eq :: QDouble -> QDouble -> Bool +foreign import javascript unsafe "$1[0] !== $2[0] || $1[1] !== $2[1] || $1[2] !== $2[2] || $1[3] !== $2[3]" + js_neq :: QDouble -> QDouble -> Bool + + + +-------------------------------------------------------------------------- +-- Show +-------------------------------------------------------------------------- + +instance Show QDouble where + show = unpack' . js_show + +foreign import javascript unsafe "$1[3].toPrecision(8)\ + \ + ($1[0] >= 0 ? ' + ' : ' - ') + Math.abs($1[0]).toPrecision(8) + 'i'\ + \ + ($1[1] >= 0 ? ' + ' : ' - ') + Math.abs($1[1]).toPrecision(8) + 'j'\ + \ + ($1[2] >= 0 ? ' + ' : ' - ') + Math.abs($1[2]).toPrecision(8) + 'k'" + js_show:: QDouble -> JSString diff --git a/src-ghcjs/Numeric/Quaternion/QFloat.hs b/src-ghcjs/Numeric/Quaternion/QFloat.hs new file mode 100644 --- /dev/null +++ b/src-ghcjs/Numeric/Quaternion/QFloat.hs @@ -0,0 +1,295 @@ +-- Note, +-- +-- The whole module is made by copying Numeric.Quaternion.QDouble and replacing: +-- Double -> Float +-- Float64 -> Float32 +-- +-- If we are doing any refactoring of one of these modules, just do the same operation. + +{-# LANGUAGE JavaScriptFFI #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE FlexibleInstances #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +module Numeric.Quaternion.QFloat + ( QFloat, Quater (..) + ) where + +import Data.JSString (JSString, unpack') +import Data.Coerce (coerce) + +import Numeric.Array +import Numeric.DataFrame.Type +import Numeric.Vector +import Numeric.Matrix + +import Numeric.Quaternion.Class + + +type QFloat = Quater Float + +instance Quaternion Float where + newtype Quater Float = QFloat (ArrayT Float '[4]) + {-# INLINE packQ #-} + packQ = js_packQ + {-# INLINE fromVecNum #-} + fromVecNum = js_fromVecNum + {-# INLINE fromVec4 #-} + fromVec4 = coerce + {-# INLINE toVec4 #-} + toVec4 = coerce + {-# INLINE unpackQ #-} + unpackQ = js_unpackQ + {-# INLINE square #-} + square = js_square + {-# INLINE im #-} + im = js_im + {-# INLINE re #-} + re = js_re + {-# INLINE imVec #-} + imVec = js_imVec + {-# INLINE taker #-} + taker = js_taker + {-# INLINE takei #-} + takei = js_takei + {-# INLINE takej #-} + takej = js_takej + {-# INLINE takek #-} + takek = js_takek + {-# INLINE conjugate #-} + conjugate = js_conjugate + {-# INLINE rotScale #-} + rotScale = js_rotScale + {-# INLINE getRotScale #-} + getRotScale = js_getRotScale + {-# INLINE axisRotation #-} + axisRotation = js_axisRotation + {-# INLINE qArg #-} + qArg = js_qArg + {-# INLINE fromMatrix33 #-} + fromMatrix33 = js_fromMatrix33 + {-# INLINE fromMatrix44 #-} + fromMatrix44 = js_fromMatrix44 + {-# INLINE toMatrix33 #-} + toMatrix33 = js_toMatrix33 + {-# INLINE toMatrix44 #-} + toMatrix44 = js_toMatrix44 + + +foreign import javascript unsafe "new Float32Array([$1,$2,$3,$4])" + js_packQ :: Float -> Float -> Float -> Float -> QFloat + +foreign import javascript unsafe "new Float32Array([$1[0],$1[1],$1[2],$2])" + js_fromVecNum :: Vector Float 3 -> Float -> QFloat + +foreign import javascript unsafe "$r1 = $1[0];$r2 = $1[1];$r3 = $1[2];$r4 = $1[3];" + js_unpackQ :: QFloat -> (Float,Float,Float,Float) + + +foreign import javascript unsafe "$1[0]*$1[0] + $1[1]*$1[1] + $1[2]*$1[2] + $1[3]*$1[3]" + js_square :: QFloat -> Float + +foreign import javascript unsafe "new Float32Array([$1[0],$1[1],$1[2],0])" + js_im :: QFloat -> QFloat + +foreign import javascript unsafe "$1.slice(0,3)" + js_imVec :: QFloat -> Vector Float 3 + +foreign import javascript unsafe "$r = new Float32Array(4); $r[3] = $1[3];" + js_re :: QFloat -> QFloat + +foreign import javascript unsafe "$1[3]" + js_taker :: QFloat -> Float + +foreign import javascript unsafe "$1[0]" + js_takei :: QFloat -> Float + +foreign import javascript unsafe "$1[1]" + js_takej :: QFloat -> Float + +foreign import javascript unsafe "$1[2]" + js_takek :: QFloat -> Float + +foreign import javascript unsafe "new Float32Array([-$1[0],-$1[1],-$1[2],$1[3]])" + js_conjugate :: QFloat -> QFloat + + +foreign import javascript unsafe "new Float32Array(h$easytensor_rotScale($1,$2))" + js_rotScale :: QFloat -> Vector Float 3 -> Vector Float 3 + +foreign import javascript unsafe "new Float32Array(h$easytensor_getRotScale($1,$2))" + js_getRotScale :: Vector Float 3 -> Vector Float 3 -> QFloat + +foreign import javascript unsafe "new Float32Array(h$easytensor_axisRotation($1,$2))" + js_axisRotation :: Vector Float 3 -> Float -> QFloat + +foreign import javascript unsafe "h$easytensor_qArg($1)" + js_qArg :: QFloat -> Float + +foreign import javascript unsafe "new Float32Array(h$easytensor_qfromMatrix33($1))" + js_fromMatrix33 :: Matrix Float 3 3 -> QFloat + +foreign import javascript unsafe "new Float32Array(h$easytensor_qfromMatrix44($1))" + js_fromMatrix44 :: Matrix Float 4 4 -> QFloat + +foreign import javascript unsafe "new Float32Array(h$easytensor_qtoMatrix33($1))" + js_toMatrix33 :: QFloat -> Matrix Float 3 3 + +foreign import javascript unsafe "new Float32Array(h$easytensor_qtoMatrix44($1))" + js_toMatrix44 :: QFloat -> Matrix Float 4 4 + + +-------------------------------------------------------------------------- +-- Num +-------------------------------------------------------------------------- + +instance Num QFloat where + {-# INLINE (+) #-} + (+) = js_plus + {-# INLINE (-) #-} + (-) = js_minus + {-# INLINE (*) #-} + (*) = js_times + {-# INLINE abs #-} + abs = js_abs + {-# INLINE signum #-} + signum = js_signum + {-# INLINE negate #-} + negate = fromVec4 . negate . toVec4 + {-# INLINE fromInteger #-} + fromInteger = js_toQuaternion . fromInteger + +foreign import javascript unsafe "$1.map(function (e, i) {return e + $2[i];})" + js_plus :: QFloat -> QFloat -> QFloat + + +foreign import javascript unsafe "$1.map(function (e, i) {return e - $2[i];})" + js_minus :: QFloat -> QFloat -> QFloat + +foreign import javascript unsafe "new Float32Array(\ + \[ $1[3]*$2[0] + $1[0]*$2[3] + $1[1]*$2[2] - $1[2]*$2[1]\ + \, $1[3]*$2[1] - $1[0]*$2[2] + $1[1]*$2[3] + $1[2]*$2[0]\ + \, $1[3]*$2[2] + $1[0]*$2[1] - $1[1]*$2[0] + $1[2]*$2[3]\ + \, $1[3]*$2[3] - $1[0]*$2[0] - $1[1]*$2[1] - $1[2]*$2[2] ])" + js_times :: QFloat -> QFloat -> QFloat + +foreign import javascript unsafe "new Float32Array([0,0,0,Math.hypot($1[0],$1[1],$1[2],$1[3])])" + js_abs :: QFloat -> QFloat + +foreign import javascript unsafe "var l = Math.hypot($1[0],$1[1],$1[2],$1[3]); $r = (l == 0) ? (new Float32Array(4)) : $1.map(function (e) {return e/l;})" + js_signum :: QFloat -> QFloat + +foreign import javascript unsafe "$r = new Float32Array(4); $r[3] = $1;" + js_toQuaternion :: Float -> QFloat + +{-# RULES +"realToFrac/FloatQFloat" realToFrac = js_toQuaternion +"realToFrac/aQFloat" realToFrac = js_toQuaternion . realToFrac +"fromIntegral/aQFloat" fromIntegral = js_toQuaternion . fromIntegral + #-} + + +-------------------------------------------------------------------------- +-- Fractional +-------------------------------------------------------------------------- + +instance Fractional QFloat where + {-# INLINE recip #-} + recip = js_recip + {-# INLINE (/) #-} + (/) p = js_times p . js_recip + {-# INLINE fromRational #-} + fromRational = js_toQuaternion . fromRational + +foreign import javascript unsafe "new Float32Array(h$easytensor_qrecip($1))" + js_recip :: QFloat -> QFloat + + +-------------------------------------------------------------------------- +-- Floating +-------------------------------------------------------------------------- + +instance Floating QFloat where + {-# INLINE pi #-} + pi = js_toQuaternion pi + {-# INLINE exp #-} + exp = js_exp + {-# INLINE log #-} + log = js_log + {-# INLINE sqrt #-} + sqrt = js_sqrt + {-# INLINE sin #-} + sin = js_sin + {-# INLINE cos #-} + cos = js_cos + {-# INLINE tan #-} + tan = js_tan + {-# INLINE sinh #-} + sinh = js_sinh + {-# INLINE cosh #-} + cosh = js_cosh + {-# INLINE tanh #-} + tanh = js_tanh + {-# INLINE asin #-} + asin q = -i * log (i*q + sqrt (1 - q*q)) + where i = signum . im $ q + {-# INLINE acos #-} + acos q = pi/2 - asin q + {-# INLINE atan #-} + atan q = if square imq == 0 + then js_toQuaternion (atan $ taker q) + else i / 2 * log ( (i + q) / (i - q) ) + where i = signum imq + imq = im q + {-# INLINE asinh #-} + asinh q = log (q + sqrt (q*q + 1)) + {-# INLINE acosh #-} + acosh q = log (q + sqrt (q*q - 1)) + {-# INLINE atanh #-} + atanh q = 0.5 * log ((1+q)/(1-q)) + +foreign import javascript unsafe "new Float32Array(h$easytensor_qexp($1))" js_exp :: QFloat -> QFloat +foreign import javascript unsafe "new Float32Array(h$easytensor_qlog($1))" js_log :: QFloat -> QFloat +foreign import javascript unsafe "new Float32Array(h$easytensor_qsqrt($1))" js_sqrt :: QFloat -> QFloat +foreign import javascript unsafe "new Float32Array(h$easytensor_qsin($1))" js_sin :: QFloat -> QFloat +foreign import javascript unsafe "new Float32Array(h$easytensor_qcos($1))" js_cos :: QFloat -> QFloat +foreign import javascript unsafe "new Float32Array(h$easytensor_qtan($1))" js_tan :: QFloat -> QFloat +foreign import javascript unsafe "new Float32Array(h$easytensor_qsinh($1))" js_sinh :: QFloat -> QFloat +foreign import javascript unsafe "new Float32Array(h$easytensor_qcosh($1))" js_cosh :: QFloat -> QFloat +foreign import javascript unsafe "new Float32Array(h$easytensor_qtanh($1))" js_tanh :: QFloat -> QFloat + + +-------------------------------------------------------------------------- +-- Eq +-------------------------------------------------------------------------- + +instance Eq QFloat where + {-# INLINE (==) #-} + (==) = js_eq + {-# INLINE (/=) #-} + (/=) = js_neq + + + +foreign import javascript unsafe "$1[0] === $2[0] && $1[1] === $2[1] && $1[2] === $2[2] && $1[3] === $2[3]" + js_eq :: QFloat -> QFloat -> Bool +foreign import javascript unsafe "$1[0] !== $2[0] || $1[1] !== $2[1] || $1[2] !== $2[2] || $1[3] !== $2[3]" + js_neq :: QFloat -> QFloat -> Bool + + + +-------------------------------------------------------------------------- +-- Show +-------------------------------------------------------------------------- + +instance Show QFloat where + show = unpack' . js_show + +foreign import javascript unsafe "$1[3].toPrecision(8)\ + \ + ($1[0] >= 0 ? ' + ' : ' - ') + Math.abs($1[0]).toPrecision(8) + 'i'\ + \ + ($1[1] >= 0 ? ' + ' : ' - ') + Math.abs($1[1]).toPrecision(8) + 'j'\ + \ + ($1[2] >= 0 ? ' + ' : ' - ') + Math.abs($1[2]).toPrecision(8) + 'k'" + js_show:: QFloat -> JSString + + diff --git a/src-ghcjs/Numeric/Quaternion/Quaternion.js b/src-ghcjs/Numeric/Quaternion/Quaternion.js new file mode 100644 --- /dev/null +++ b/src-ghcjs/Numeric/Quaternion/Quaternion.js @@ -0,0 +1,205 @@ +function h$easytensor_rotScale(quat, vec) { + 'use strict'; + var i = quat[0], j = quat[1], k = quat[2], t = quat[3]; + var x = vec[0], y = vec[1], z = vec[2]; + var l = t*t - i*i - j*j - k*k; + var d = 2*(i*x + j*y + k*z); + t *= 2; + return [ l*x + d*i + t*(z*j - y*k) + , l*y + d*j + t*(x*k - z*i) + , l*z + d*k + t*(y*i - x*j) + ]; +} + +function h$easytensor_qArg(quat) { + 'use strict'; + return Math.atan2( Math.hypot(quat[0],quat[1],quat[2]) , quat[3] ) * 2 ; +} + +function h$easytensor_getRotScale(a, b) { + 'use strict'; + if (b[0] === 0 && b[1] === 0 && b[2] === 0) { return [0,0,0,0];} + if (a[0] === 0 && a[1] === 0 && a[2] === 0) { return [Infinity,Infinity,Infinity,Infinity];} + var t = h$easytensor_cross(a, b); + var ma = Math.hypot(a[0],a[1],a[2]); + var mb = Math.hypot(b[0],b[1],b[2]); + var dot = a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; + if (t[0] === 0 && t[1] === 0 && t[2] === 0) { + if (dot > 0) {return [0,0,0,Math.sqrt(mb/ma)];} + else {return [0,0,Math.sqrt(mb/ma),0];} + } + var c = Math.sqrt(ma*mb + dot); + ma *= Math.SQRT2; + return [ t[0]/(ma*c) + , t[1]/(ma*c) + , t[2]/(ma*c) + , c/ma + ]; +} + +function h$easytensor_axisRotation(axis, a) { + 'use strict'; + if (axis[0] === 0 && axis[1] === 0 && axis[2] === 0) { return [0,0,0,1];} + var c = Math.cos(a*0.5), s = Math.sin(a*0.5) / Math.hypot(axis[0],axis[1],axis[2]); + return [ axis[0]*s, axis[1]*s, axis[2]*s, c]; +} + +function h$easytensor_qfromMatrix33(m) { + 'use strict'; + var d = Math.cbrt( + m[0]*(m[4]*m[8]-m[5]*m[7]) + - m[1]*(m[3]*m[8]-m[5]*m[6]) + + m[2]*(m[3]*m[7]-m[4]*m[6])); + return [ Math.sqrt(Math.max( 0, d + m[0] - m[4] - m[8] )) * Math.sign(m[5] - m[7]) * 0.5 + , Math.sqrt(Math.max( 0, d - m[0] + m[4] - m[8] )) * Math.sign(m[6] - m[2]) * 0.5 + , Math.sqrt(Math.max( 0, d - m[0] - m[4] + m[8] )) * Math.sign(m[1] - m[3]) * 0.5 + , Math.sqrt(Math.max( 0, d + m[0] + m[4] + m[8] )) * 0.5 ]; +} + +function h$easytensor_qfromMatrix44(m) { + 'use strict'; + var d = Math.cbrt( + m[0]*(m[5]*m[10]-m[6]*m[9]) + - m[1]*(m[4]*m[10]-m[6]*m[8]) + + m[2]*(m[4]*m[ 9]-m[5]*m[8])); + return [ Math.sqrt(Math.max( 0, d + m[0] - m[5] - m[10] )) * Math.sign(m[6] - m[9]) * 0.5 / m[15] + , Math.sqrt(Math.max( 0, d - m[0] + m[5] - m[10] )) * Math.sign(m[8] - m[2]) * 0.5 / m[15] + , Math.sqrt(Math.max( 0, d - m[0] - m[5] + m[10] )) * Math.sign(m[1] - m[4]) * 0.5 / m[15] + , Math.sqrt(Math.max( 0, d + m[0] + m[5] + m[10] )) * 0.5 / m[15] ]; +} + +function h$easytensor_qtoMatrix33(quat) { + 'use strict'; + var x = quat[0], y = quat[1], z = quat[2], w = quat[3]; + var w2 = w*w; + if (x === 0 && y === 0 && z === 0) { + return [w2,0,0 + ,0,w2,0 + ,0,0,w2]; + } + var x2 = x*x, y2 = y*y, z2 = z*z; + var l2 = x2+y2+z2+w2; + return [ l2 - 2*(z2 + y2), 2*(x*y + z*w), 2*(x*z - y*w) + , 2*(x*y - z*w), l2 - 2*(z2 + x2), 2*(y*z + x*w) + , 2*(x*z + y*w), 2*(y*z - x*w), l2 - 2*(y2 + x2) ]; +} + +function h$easytensor_qtoMatrix44(quat) { + 'use strict'; + var x = quat[0], y = quat[1], z = quat[2], w = quat[3]; + var w2 = w*w; + if (x === 0 && y === 0 && z === 0) { + return [w2,0,0,0 + ,0,w2,0,0 + ,0,0,w2,0 + ,0,0,0,1]; + } + var x2 = x*x, y2 = y*y, z2 = z*z; + var l2 = x2+y2+z2+w2; + return [ l2 - 2*(z2 + y2), 2*(x*y + z*w), 2*(x*z - y*w), 0 + , 2*(x*y - z*w), l2 - 2*(z2 + x2), 2*(y*z + x*w), 0 + , 2*(x*z + y*w), 2*(y*z - x*w), l2 - 2*(y2 + x2), 0 + , 0, 0, 0, 1]; +} + +function h$easytensor_qrecip(q) { + 'use strict'; + var c = -1 / (q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]); + return [q[0]*c,q[1]*c,q[2]*c,-q[3]*c]; +} + +function h$easytensor_qexp(q) { + 'use strict'; + var mv = Math.hypot(q[0],q[1],q[2]), et = Math.exp(q[3]); + if(mv === 0) {return [0,0,0,et];} + var l = et * Math.sin(mv) / mv; + return [q[0]*l,q[1]*l,q[2]*l,et*Math.cos(mv)]; +} +function h$easytensor_qlog(q) { + 'use strict'; + var mv = q[0]*q[0] + q[1]*q[1] + q[2]*q[2]; + if(mv === 0) { + if(q[3] >= 0){ + return [0,0,0,Math.log(q[3])]; + } else { + return [Math.PI,0,0,Math.log(-q[3])]; + } + } + var mq = Math.sqrt(mv + q[3]*q[3]); + mv = Math.sqrt(mv); + var l = Math.atan2( mv, q[3] ) / mv; + return [q[0]*l,q[1]*l,q[2]*l,Math.log(mq)]; +} +function h$easytensor_qsqrt(q) { + 'use strict'; + var mv = q[0]*q[0] + q[1]*q[1] + q[2]*q[2]; + if(mv === 0) { + if(q[3] >= 0){ + return [0,0,0,Math.sqrt(q[3])]; + } else { + return [Math.sqrt(-q[3]),0,0,0]; + } + } + var l = Math.sqrt(mv + q[3]*q[3]); + var l2 = Math.sqrt(l); + var tq = q[3] / (l * 2); + var sina = Math.sqrt(0.5 - tq) * l2 / Math.sqrt(mv); + return [q[0]*sina,q[1]*sina,q[2]*sina,Math.sqrt(0.5 + tq) * l2]; +} + +// A good tutorial on complex number trigonometric functions is available here +// http://www.milefoot.com/math/complex/functionsofi.htm +// I extend it to complex numbers by replacing complex i with quaternion vector ijk + +function h$easytensor_qsin(q) { + 'use strict'; + var mv = q[0]*q[0] + q[1]*q[1] + q[2]*q[2]; + if(mv === 0) {return [0,0,0,Math.sin(q[3])];} + mv = Math.sqrt(mv); + var l = Math.cos(q[3]) * Math.sinh(mv) / mv; + return [q[0]*l,q[1]*l,q[2]*l, Math.sin(q[3])*Math.cosh(mv)]; +} +function h$easytensor_qcos(q) { + 'use strict'; + var mv = q[0]*q[0] + q[1]*q[1] + q[2]*q[2]; + if(mv === 0) {return [0,0,0,Math.cos(q[3])];} + mv = Math.sqrt(mv); + var l = - Math.sin(q[3]) * Math.sinh(mv) / mv; + return [q[0]*l,q[1]*l,q[2]*l, Math.cos(q[3])*Math.cosh(mv)]; +} +function h$easytensor_qtan(q) { + 'use strict'; + var mv = q[0]*q[0] + q[1]*q[1] + q[2]*q[2]; + if(mv === 0) {return [0,0,0,Math.tan(q[3])];} + mv = Math.sqrt(mv); + var chv = Math.cosh(mv), shv = Math.sinh(mv), ct = Math.cos(q[3]), st = Math.sin(q[3]); + var cq = 1 / (ct*ct*chv*chv + st*st*shv*shv); + var l = chv * shv * cq / mv; + return [q[0]*l,q[1]*l,q[2]*l, ct * st * cq]; +} +function h$easytensor_qsinh(q) { + 'use strict'; + var mv = q[0]*q[0] + q[1]*q[1] + q[2]*q[2]; + if(mv === 0) {return [0,0,0,Math.sinh(q[3])];} + mv = Math.sqrt(mv); + var l = Math.cosh(q[3]) * Math.sin(mv) / mv; + return [q[0]*l,q[1]*l,q[2]*l, Math.sinh(q[3])*Math.cos(mv)]; +} +function h$easytensor_qcosh(q) { + 'use strict'; + var mv = q[0]*q[0] + q[1]*q[1] + q[2]*q[2]; + if(mv === 0) {return [0,0,0,Math.cosh(q[3])];} + mv = Math.sqrt(mv); + var l = Math.sinh(q[3]) * Math.sin(mv) / mv; + return [q[0]*l,q[1]*l,q[2]*l, Math.cosh(q[3])*Math.cos(mv)]; +} +function h$easytensor_qtanh(q) { + 'use strict'; + var mv = q[0]*q[0] + q[1]*q[1] + q[2]*q[2]; + if(mv === 0) {return [0,0,0,Math.tanh(q[3])];} + mv = Math.sqrt(mv); + var cv = Math.cos(mv), sv = Math.sin(mv), cht = Math.cosh(q[3]), sht = Math.sinh(q[3]); + var cq = 1 / (cht*cht*cv*cv + sht*sht*sv*sv); + var l = cv * sv * cq / mv; + return [q[0]*l,q[1]*l,q[2]*l, cht * sht * cq]; +} diff --git a/src/Numeric/Commons.hs b/src/Numeric/Commons.hs --- a/src/Numeric/Commons.hs +++ b/src/Numeric/Commons.hs @@ -12,6 +12,10 @@ {-# LANGUAGE TypeInType #-} {-# LANGUAGE UnboxedTuples #-} {-# LANGUAGE UndecidableInstances #-} +#ifdef ghcjs_HOST_OS +{-# LANGUAGE JavaScriptFFI #-} +{-# LANGUAGE UnliftedFFITypes #-} +#endif ----------------------------------------------------------------------------- -- | -- Module : Numeric.Commons @@ -25,13 +29,14 @@ module Numeric.Commons ( ElemRep, ElemPrim - , PrimBytes (..), FloatBytes, DoubleBytes, IntBytes, WordBytes + , PrimBytes (..) ) where #include "MachDeps.h" -#include "HsBaseConfig.h" +#ifndef ghcjs_HOST_OS import GHC.Base (runRW#) +#endif import GHC.Int (Int16 (..), Int32 (..), Int64 (..), Int8 (..)) import GHC.Prim import GHC.Types (Double (..), Float (..), Int (..), RuntimeRep (..), @@ -47,7 +52,7 @@ type instance ElemRep Int8 = 'IntRep type instance ElemRep Int16 = 'IntRep type instance ElemRep Int32 = 'IntRep -#if SIZEOF_HSWORD < 8 +#if WORD_SIZE_IN_BITS < 64 type instance ElemRep Int64 = 'Int64Rep #else type instance ElemRep Int64 = 'IntRep @@ -56,7 +61,7 @@ type instance ElemRep Word8 = 'WordRep type instance ElemRep Word16 = 'WordRep type instance ElemRep Word32 = 'WordRep -#if SIZEOF_HSWORD < 8 +#if WORD_SIZE_IN_BITS < 64 type instance ElemRep Word64 = 'Word64Rep #else type instance ElemRep Word64 = 'WordRep @@ -69,7 +74,7 @@ type instance ElemPrim Int8 = Int# type instance ElemPrim Int16 = Int# type instance ElemPrim Int32 = Int# -#if SIZEOF_HSWORD < 8 +#if WORD_SIZE_IN_BITS < 64 type instance ElemPrim Int64 = Int64# #else type instance ElemPrim Int64 = Int# @@ -78,18 +83,13 @@ type instance ElemPrim Word8 = Word# type instance ElemPrim Word16 = Word# type instance ElemPrim Word32 = Word# -#if SIZEOF_HSWORD < 8 +#if WORD_SIZE_IN_BITS < 64 type instance ElemPrim Word64 = Word64# #else type instance ElemPrim Word64 = Word# #endif -type FloatBytes a = (PrimBytes a, ElemRep a ~ 'FloatRep , ElemPrim a ~ Float#) -type DoubleBytes a = (PrimBytes a, ElemRep a ~ 'DoubleRep, ElemPrim a ~ Double#) -type IntBytes a = (PrimBytes a, ElemRep a ~ 'IntRep , ElemPrim a ~ Int#) -type WordBytes a = (PrimBytes a, ElemRep a ~ 'WordRep , ElemPrim a ~ Word#) - -- | Facilities to convert to and from raw byte array. -- Warning! offsets and sizes are in elements, not in bytes! -- Therefore one must be really carefull if having a crazy idea of @@ -111,6 +111,12 @@ ix :: Int# -> a -> (ElemPrim a :: TYPE (ElemRep a)) instance PrimBytes Float where +#ifdef ghcjs_HOST_OS + toBytes v = (# 0#, 1#, js_wrapFloat v #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = js_unwrapFloat arr off + {-# INLINE fromBytes #-} +#else toBytes v@(F# x) = case runRW# ( \s0 -> case newByteArray# (byteSize v) s0 of (# s1, marr #) -> case writeFloatArray# marr 0# x s1 of @@ -119,6 +125,7 @@ {-# INLINE toBytes #-} fromBytes (# off, _, arr #) = F# (indexFloatArray# arr off) {-# INLINE fromBytes #-} +#endif byteSize _ = SIZEOF_HSFLOAT# {-# INLINE byteSize #-} byteAlign _ = ALIGNMENT_HSFLOAT# @@ -129,6 +136,12 @@ {-# INLINE ix #-} instance PrimBytes Double where +#ifdef ghcjs_HOST_OS + toBytes v = (# 0#, 1#, js_wrapDouble v #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = js_unwrapDouble arr off + {-# INLINE fromBytes #-} +#else toBytes v@(D# x) = case runRW# ( \s0 -> case newByteArray# (byteSize v) s0 of (# s1, marr #) -> case writeDoubleArray# marr 0# x s1 of @@ -137,6 +150,7 @@ {-# INLINE toBytes #-} fromBytes (# off, _, arr #) = D# (indexDoubleArray# arr off) {-# INLINE fromBytes #-} +#endif byteSize _ = SIZEOF_HSDOUBLE# {-# INLINE byteSize #-} byteAlign _ = ALIGNMENT_HSDOUBLE# @@ -147,6 +161,12 @@ {-# INLINE ix #-} instance PrimBytes Int where +#ifdef ghcjs_HOST_OS + toBytes v = (# 0#, 1#, js_wrapInt v #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = js_unwrapInt arr off + {-# INLINE fromBytes #-} +#else toBytes v@(I# x) = case runRW# ( \s0 -> case newByteArray# (byteSize v) s0 of (# s1, marr #) -> case writeIntArray# marr 0# x s1 of @@ -155,6 +175,7 @@ {-# INLINE toBytes #-} fromBytes (# off, _, arr #) = I# (indexIntArray# arr off) {-# INLINE fromBytes #-} +#endif byteSize _ = SIZEOF_HSINT# {-# INLINE byteSize #-} byteAlign _ = ALIGNMENT_HSINT# @@ -165,6 +186,12 @@ {-# INLINE ix #-} instance PrimBytes Int8 where +#ifdef ghcjs_HOST_OS + toBytes v = (# 0#, 1#, js_wrapInt8 v #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = js_unwrapInt8 arr off + {-# INLINE fromBytes #-} +#else toBytes v@(I8# x) = case runRW# ( \s0 -> case newByteArray# (byteSize v) s0 of (# s1, marr #) -> case writeInt8Array# marr 0# x s1 of @@ -173,6 +200,7 @@ {-# INLINE toBytes #-} fromBytes (# off, _, arr #) = I8# (indexInt8Array# arr off) {-# INLINE fromBytes #-} +#endif byteSize _ = SIZEOF_INT8# {-# INLINE byteSize #-} byteAlign _ = ALIGNMENT_INT8# @@ -183,6 +211,12 @@ {-# INLINE ix #-} instance PrimBytes Int16 where +#ifdef ghcjs_HOST_OS + toBytes v = (# 0#, 1#, js_wrapInt16 v #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = js_unwrapInt16 arr off + {-# INLINE fromBytes #-} +#else toBytes v@(I16# x) = case runRW# ( \s0 -> case newByteArray# (byteSize v) s0 of (# s1, marr #) -> case writeInt16Array# marr 0# x s1 of @@ -191,6 +225,7 @@ {-# INLINE toBytes #-} fromBytes (# off, _, arr #) = I16# (indexInt16Array# arr off) {-# INLINE fromBytes #-} +#endif byteSize _ = SIZEOF_INT16# {-# INLINE byteSize #-} byteAlign _ = ALIGNMENT_INT16# @@ -201,6 +236,12 @@ {-# INLINE ix #-} instance PrimBytes Int32 where +#ifdef ghcjs_HOST_OS + toBytes v = (# 0#, 1#, js_wrapInt32 v #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = js_unwrapInt32 arr off + {-# INLINE fromBytes #-} +#else toBytes v@(I32# x) = case runRW# ( \s0 -> case newByteArray# (byteSize v) s0 of (# s1, marr #) -> case writeInt32Array# marr 0# x s1 of @@ -209,6 +250,7 @@ {-# INLINE toBytes #-} fromBytes (# off, _, arr #) = I32# (indexInt32Array# arr off) {-# INLINE fromBytes #-} +#endif byteSize _ = SIZEOF_INT32# {-# INLINE byteSize #-} byteAlign _ = ALIGNMENT_INT32# @@ -239,6 +281,12 @@ #endif instance PrimBytes Word where +#ifdef ghcjs_HOST_OS + toBytes v = (# 0#, 1#, js_wrapWord v #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = js_unwrapWord arr off + {-# INLINE fromBytes #-} +#else toBytes v@(W# x) = case runRW# ( \s0 -> case newByteArray# (byteSize v) s0 of (# s1, marr #) -> case writeWordArray# marr 0# x s1 of @@ -247,6 +295,7 @@ {-# INLINE toBytes #-} fromBytes (# off, _, arr #) = W# (indexWordArray# arr off) {-# INLINE fromBytes #-} +#endif byteSize _ = SIZEOF_HSWORD# {-# INLINE byteSize #-} byteAlign _ = ALIGNMENT_HSWORD# @@ -257,6 +306,12 @@ {-# INLINE ix #-} instance PrimBytes Word8 where +#ifdef ghcjs_HOST_OS + toBytes v = (# 0#, 1#, js_wrapWord8 v #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = js_unwrapWord8 arr off + {-# INLINE fromBytes #-} +#else toBytes v@(W8# x) = case runRW# ( \s0 -> case newByteArray# (byteSize v) s0 of (# s1, marr #) -> case writeWord8Array# marr 0# x s1 of @@ -265,6 +320,7 @@ {-# INLINE toBytes #-} fromBytes (# off, _, arr #) = W8# (indexWord8Array# arr off) {-# INLINE fromBytes #-} +#endif byteSize _ = SIZEOF_WORD8# {-# INLINE byteSize #-} byteAlign _ = ALIGNMENT_WORD8# @@ -275,6 +331,12 @@ {-# INLINE ix #-} instance PrimBytes Word16 where +#ifdef ghcjs_HOST_OS + toBytes v = (# 0#, 1#, js_wrapWord16 v #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = js_unwrapWord16 arr off + {-# INLINE fromBytes #-} +#else toBytes v@(W16# x) = case runRW# ( \s0 -> case newByteArray# (byteSize v) s0 of (# s1, marr #) -> case writeWord16Array# marr 0# x s1 of @@ -283,6 +345,7 @@ {-# INLINE toBytes #-} fromBytes (# off, _, arr #) = W16# (indexWord16Array# arr off) {-# INLINE fromBytes #-} +#endif byteSize _ = SIZEOF_WORD16# {-# INLINE byteSize #-} byteAlign _ = ALIGNMENT_WORD16# @@ -293,6 +356,12 @@ {-# INLINE ix #-} instance PrimBytes Word32 where +#ifdef ghcjs_HOST_OS + toBytes v = (# 0#, 1#, js_wrapWord32 v #) + {-# INLINE toBytes #-} + fromBytes (# off, _, arr #) = js_unwrapWord32 arr off + {-# INLINE fromBytes #-} +#else toBytes v@(W32# x) = case runRW# ( \s0 -> case newByteArray# (byteSize v) s0 of (# s1, marr #) -> case writeWord32Array# marr 0# x s1 of @@ -301,6 +370,7 @@ {-# INLINE toBytes #-} fromBytes (# off, _, arr #) = W32# (indexWord32Array# arr off) {-# INLINE fromBytes #-} +#endif byteSize _ = SIZEOF_WORD32# {-# INLINE byteSize #-} byteAlign _ = ALIGNMENT_WORD32# @@ -331,3 +401,28 @@ {-# INLINE ix #-} #endif +#ifdef ghcjs_HOST_OS +foreign import javascript unsafe "h$wrapBuffer((new Float32Array([$1])).buffer)" js_wrapFloat :: Float -> ByteArray# +foreign import javascript unsafe "h$wrapBuffer((new Float64Array([$1])).buffer)" js_wrapDouble :: Double -> ByteArray# +foreign import javascript unsafe "h$wrapBuffer((new Int32Array([$1])).buffer)" js_wrapInt :: Int -> ByteArray# +foreign import javascript unsafe "h$wrapBuffer((new Int32Array([$1])).buffer)" js_wrapInt32 :: Int32 -> ByteArray# +foreign import javascript unsafe "h$wrapBuffer((new Int16Array([$1])).buffer)" js_wrapInt16 :: Int16 -> ByteArray# +foreign import javascript unsafe "h$wrapBuffer((new Int8Array([$1])).buffer)" js_wrapInt8 :: Int8 -> ByteArray# +foreign import javascript unsafe "h$wrapBuffer((new Uint32Array([$1])).buffer)" js_wrapWord :: Word -> ByteArray# +foreign import javascript unsafe "h$wrapBuffer((new Uint32Array([$1])).buffer)" js_wrapWord32 :: Word32 -> ByteArray# +foreign import javascript unsafe "h$wrapBuffer((new Uint16Array([$1])).buffer)" js_wrapWord16 :: Word16 -> ByteArray# +foreign import javascript unsafe "h$wrapBuffer((new Uint8Array([$1])).buffer)" js_wrapWord8 :: Word8 -> ByteArray# + + + +foreign import javascript unsafe "($1.f3 || new Float32Array($1.buf))[$2]" js_unwrapFloat :: ByteArray# -> Int# -> Float +foreign import javascript unsafe "($1.f6 || new Float64Array($1.buf))[$2]" js_unwrapDouble :: ByteArray# -> Int# -> Double +foreign import javascript unsafe "($1.i3 || new Int32Array($1.buf))[$2]" js_unwrapInt :: ByteArray# -> Int# -> Int +foreign import javascript unsafe "($1.i3 || new Int32Array($1.buf))[$2]" js_unwrapInt32 :: ByteArray# -> Int# -> Int32 +foreign import javascript unsafe "($1.i1 || new Int16Array($1.buf))[$2]" js_unwrapInt16 :: ByteArray# -> Int# -> Int16 +foreign import javascript unsafe "($1.i8 || new Int8Array($1.buf))[$2]" js_unwrapInt8 :: ByteArray# -> Int# -> Int8 +foreign import javascript unsafe "($1.u3 || new Uint32Array($1.buf))[$2]" js_unwrapWord :: ByteArray# -> Int# -> Word +foreign import javascript unsafe "($1.u3 || new Uint32Array($1.buf))[$2]" js_unwrapWord32 :: ByteArray# -> Int# -> Word32 +foreign import javascript unsafe "($1.u1 || new Uint16Array($1.buf))[$2]" js_unwrapWord16 :: ByteArray# -> Int# -> Word16 +foreign import javascript unsafe "($1.u8 || new Uint8Array($1.buf))[$2]" js_unwrapWord8 :: ByteArray# -> Int# -> Word8 +#endif diff --git a/src/Numeric/DataFrame.hs b/src/Numeric/DataFrame.hs --- a/src/Numeric/DataFrame.hs +++ b/src/Numeric/DataFrame.hs @@ -41,3 +41,5 @@ import Numeric.Matrix import Numeric.Scalar import Numeric.Vector + + diff --git a/src/Numeric/DataFrame/IO.hs b/src/Numeric/DataFrame/IO.hs --- a/src/Numeric/DataFrame/IO.hs +++ b/src/Numeric/DataFrame/IO.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE ExistentialQuantification #-} {-# LANGUAGE FlexibleContexts #-} @@ -9,6 +10,7 @@ {-# LANGUAGE TypeApplications #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE UnboxedTuples #-} +{-# LANGUAGE TypeOperators #-} ----------------------------------------------------------------------------- -- | -- Module : Numeric.DataFrame.IO @@ -17,22 +19,46 @@ -- -- Maintainer : chirkin@arch.ethz.ch -- --- Mutable DataFrames living in IO. +---- Mutable DataFrames living in IO. -- ----------------------------------------------------------------------------- module Numeric.DataFrame.IO - ( MutableFrame (), IODataFrame () + ( +#ifdef ghcjs_HOST_OS + MutableFrame (), IODataFrame (..), MDataFrame (..), MutableArrayT (..) +#else + MutableFrame (), IODataFrame () +#endif + , SomeIODataFrame (..) , newDataFrame, copyDataFrame, copyMutableDataFrame , unsafeFreezeDataFrame , freezeDataFrame, thawDataFrame , writeDataFrame, readDataFrame , writeDataFrameOff, readDataFrameOff +#ifdef ghcjs_HOST_OS + -- * JavaScript-specific functions + , MutableArrayBuffer + , newArrayBuffer, arrayBuffer, viewFloatArray, viewDoubleArray + , viewIntArray, viewInt32Array, viewInt16Array, viewInt8Array + , viewWordArray, viewWord32Array, viewWord16Array, viewWord8Array, viewWord8ClampedArray +#endif ) where import GHC.Prim (RealWorld) import GHC.Types (Int (..), IO (..)) + +#ifdef ghcjs_HOST_OS +import Numeric.Array.Family hiding (Scalar) +import JavaScript.TypedArray.ArrayBuffer +import GHC.Prim +import Data.Int +import Data.Word +import Data.Maybe +import GHCJS.Types +import Numeric.DataFrame.Inference +#endif import Numeric.Commons import Numeric.DataFrame.Type import Numeric.DataFrame.Mutable @@ -42,32 +68,53 @@ -- | Mutable DataFrame that lives in IO. -- Internal representation is always a ByteArray. newtype IODataFrame t (ns :: [Nat]) = IODataFrame (MDataFrame RealWorld t (ns :: [Nat])) - +#ifdef ghcjs_HOST_OS +instance IsJSVal (IODataFrame t ds) +#endif +-- | Mutable DataFrame of unknown dimensionality +data SomeIODataFrame t (xns :: [XNat]) + = forall (ns :: [Nat]) + . ( FixedDim xns ns ~ ns + , FixedXDim xns ns ~ xns + , NumericFrame t ns + ) + => SomeIODataFrame (IODataFrame t ns) -- | Create a new mutable DataFrame. newDataFrame :: forall t (ns :: [Nat]) +#ifdef ghcjs_HOST_OS + . ( ElemTypeInference t, Dimensions ns) +#else . ( PrimBytes t, Dimensions ns) +#endif => IO (IODataFrame t ns) newDataFrame = IODataFrame <$> IO (newDataFrame# @t @ns) {-# INLINE newDataFrame #-} -- | Copy one DataFrame into another mutable DataFrame at specified position. -copyDataFrame :: forall t (as :: [Nat]) (bs :: [Nat]) (asbs :: [Nat]) - . ( PrimBytes (DataFrame t as) - , ConcatList as bs asbs - , Dimensions bs - ) - => DataFrame t as -> Idx bs -> IODataFrame t asbs -> IO () +copyDataFrame :: forall t (as :: [Nat]) (b' :: Nat) (b :: Nat) (bs :: [Nat]) (asbs :: [Nat]) + . ( ConcatList as (b :+ bs) asbs, Dimensions (b :+ bs) +#ifdef ghcjs_HOST_OS + , ArraySizeInference (as +: b'), Dimensions as +#else + , PrimBytes (DataFrame t (as +: b')) +#endif + ) + => DataFrame t (as +: b') -> Idx (b :+ bs) -> IODataFrame t asbs -> IO () copyDataFrame df ei (IODataFrame mdf) = IO (copyDataFrame# df ei mdf) {-# INLINE copyDataFrame #-} + -- | Copy one mutable DataFrame into another mutable DataFrame at specified position. -copyMutableDataFrame :: forall t (as :: [Nat]) (bs :: [Nat]) (asbs :: [Nat]) +copyMutableDataFrame :: forall t (as :: [Nat]) (b' :: Nat) (b :: Nat) (bs :: [Nat]) (asbs :: [Nat]) . ( PrimBytes t - , ConcatList as bs asbs - , Dimensions bs + , ConcatList as (b :+ bs) asbs + , Dimensions (b :+ bs) +#ifdef ghcjs_HOST_OS + , Dimensions as +#endif ) - => IODataFrame t as -> Idx bs -> IODataFrame t asbs -> IO () + => IODataFrame t (as +: b') -> Idx (b :+ bs) -> IODataFrame t asbs -> IO () copyMutableDataFrame (IODataFrame mdfA) ei (IODataFrame mdfB) = IO (copyMDataFrame# mdfA ei mdfB) {-# INLINE copyMutableDataFrame #-} @@ -75,7 +122,11 @@ -- | Make a mutable DataFrame immutable, without copying. unsafeFreezeDataFrame :: forall t (ns :: [Nat]) - . PrimBytes (DataFrame t ns) +#ifdef ghcjs_HOST_OS + . (MutableFrame t ns, ArraySizeInference ns) +#else + . PrimBytes (DataFrame t ns) +#endif => IODataFrame t ns -> IO (DataFrame t ns) unsafeFreezeDataFrame (IODataFrame mdf) = IO (unsafeFreezeDataFrame# mdf) {-# INLINE unsafeFreezeDataFrame #-} @@ -83,14 +134,22 @@ -- | Copy content of a mutable DataFrame into a new immutable DataFrame. freezeDataFrame :: forall t (ns :: [Nat]) - . PrimBytes (DataFrame t ns) - => IODataFrame t ns -> IO (DataFrame t ns) +#ifdef ghcjs_HOST_OS + . (MutableFrame t ns, ArraySizeInference ns) +#else + . PrimBytes (DataFrame t ns) +#endif + => IODataFrame t ns -> IO (DataFrame t ns) freezeDataFrame (IODataFrame mdf) = IO (freezeDataFrame# mdf) {-# INLINE freezeDataFrame #-} -- | Create a new mutable DataFrame and copy content of immutable one in there. thawDataFrame :: forall t (ns :: [Nat]) - . PrimBytes (DataFrame t ns) +#ifdef ghcjs_HOST_OS + . (MutableFrame t ns, ArrayInstanceInference t ns) +#else + . PrimBytes (DataFrame t ns) +#endif => DataFrame t ns -> IO (IODataFrame t ns) thawDataFrame df = IODataFrame <$> IO (thawDataFrame# df) {-# INLINE thawDataFrame #-} @@ -126,3 +185,168 @@ => IODataFrame t ns -> Int -> IO (Scalar t) readDataFrameOff (IODataFrame mdf) (I# i) = scalar <$> IO (readDataFrameOff# mdf i) {-# INLINE readDataFrameOff #-} + + +#ifdef ghcjs_HOST_OS +newArrayBuffer :: Int -> IO MutableArrayBuffer +newArrayBuffer n = unsafeCoerce# <$> IO (newArrayBuffer# n) + +viewFloatArray :: forall ds + . ( Dimensions ds, ArraySizeInference ds) + => MutableArrayBuffer -> IO (SomeIODataFrame Float (AsXDims ds +: XN 0)) +viewFloatArray ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Float))) (dim @ds) ab + df <- fmap IODataFrame . IO $ viewFloatArray# (jsval ab) :: IO (IODataFrame Float (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Float ds) pn + of + Evidence -> case inferNumericFrame @Float @(ds +: n) of + Evidence -> SomeIODataFrame df + +viewDoubleArray :: forall ds + . ( Dimensions ds, ArraySizeInference ds) + => MutableArrayBuffer -> IO (SomeIODataFrame Double (AsXDims ds +: XN 0)) +viewDoubleArray ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Double))) (dim @ds) ab + df <- fmap IODataFrame . IO $ viewDoubleArray# (jsval ab) :: IO (IODataFrame Double (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Double ds) pn + of + Evidence -> case inferNumericFrame @Double @(ds +: n) of + Evidence -> SomeIODataFrame df + +viewIntArray :: forall ds + . ( Dimensions ds, ArraySizeInference ds) + => MutableArrayBuffer -> IO (SomeIODataFrame Int (AsXDims ds +: XN 0)) +viewIntArray ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Int))) (dim @ds) ab + df <- fmap IODataFrame . IO $ viewIntArray# (jsval ab) :: IO (IODataFrame Int (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Int ds) pn + of + Evidence -> case inferNumericFrame @Int @(ds +: n) of + Evidence -> SomeIODataFrame df + +viewInt32Array :: forall ds + . ( Dimensions ds, ArraySizeInference ds) + => MutableArrayBuffer -> IO (SomeIODataFrame Int32 (AsXDims ds +: XN 0)) +viewInt32Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Int32))) (dim @ds) ab + df <- fmap IODataFrame . IO $ viewInt32Array# (jsval ab) :: IO (IODataFrame Int32 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Int32 ds) pn + of + Evidence -> case inferNumericFrame @Int32 @(ds +: n) of + Evidence -> SomeIODataFrame df + +viewInt16Array :: forall ds + . ( Dimensions ds, ArraySizeInference ds) + => MutableArrayBuffer -> IO (SomeIODataFrame Int16 (AsXDims ds +: XN 0)) +viewInt16Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Int16))) (dim @ds) ab + df <- fmap IODataFrame . IO $ viewInt16Array# (jsval ab) :: IO (IODataFrame Int16 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Int16 ds) pn + of + Evidence -> case inferNumericFrame @Int16 @(ds +: n) of + Evidence -> SomeIODataFrame df + +viewInt8Array :: forall ds + . ( Dimensions ds, ArraySizeInference ds) + => MutableArrayBuffer -> IO (SomeIODataFrame Int8 (AsXDims ds +: XN 0)) +viewInt8Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Int8))) (dim @ds) ab + df <- fmap IODataFrame . IO $ viewInt8Array# (jsval ab) :: IO (IODataFrame Int8 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Int8 ds) pn + of + Evidence -> case inferNumericFrame @Int8 @(ds +: n) of + Evidence -> SomeIODataFrame df + +viewWordArray :: forall ds + . ( Dimensions ds, ArraySizeInference ds) + => MutableArrayBuffer -> IO (SomeIODataFrame Word (AsXDims ds +: XN 0)) +viewWordArray ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Word))) (dim @ds) ab + df <- fmap IODataFrame . IO $ viewWordArray# (jsval ab) :: IO (IODataFrame Word (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Word ds) pn + of + Evidence -> case inferNumericFrame @Word @(ds +: n) of + Evidence -> SomeIODataFrame df + +viewWord32Array :: forall ds + . ( Dimensions ds, ArraySizeInference ds) + => MutableArrayBuffer -> IO (SomeIODataFrame Word32 (AsXDims ds +: XN 0)) +viewWord32Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Word32))) (dim @ds) ab + df <- fmap IODataFrame . IO $ viewWord32Array# (jsval ab) :: IO (IODataFrame Word32 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Word32 ds) pn + of + Evidence -> case inferNumericFrame @Word32 @(ds +: n) of + Evidence -> SomeIODataFrame df + +viewWord16Array :: forall ds + . ( Dimensions ds, ArraySizeInference ds) + => MutableArrayBuffer -> IO (SomeIODataFrame Word16 (AsXDims ds +: XN 0)) +viewWord16Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Word16))) (dim @ds) ab + df <- fmap IODataFrame . IO $ viewWord16Array# (jsval ab) :: IO (IODataFrame Word16 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Word16 ds) pn + of + Evidence -> case inferNumericFrame @Word16 @(ds +: n) of + Evidence -> SomeIODataFrame df + +viewWord8Array :: forall ds + . ( Dimensions ds, ArraySizeInference ds) + => MutableArrayBuffer -> IO (SomeIODataFrame Word8 (AsXDims ds +: XN 0)) +viewWord8Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Word8))) (dim @ds) ab + df <- fmap IODataFrame . IO $ viewWord8Array# (jsval ab) :: IO (IODataFrame Word8 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Word8 ds) pn + of + Evidence -> case inferNumericFrame @Word8 @(ds +: n) of + Evidence -> SomeIODataFrame df + +viewWord8ClampedArray :: forall ds + . ( Dimensions ds, ArraySizeInference ds) + => MutableArrayBuffer -> IO (SomeIODataFrame Word8Clamped (AsXDims ds +: XN 0)) +viewWord8ClampedArray ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Word8Clamped))) (dim @ds) ab + df <- fmap IODataFrame . IO $ viewWord8ClampedArray# (jsval ab) :: IO (IODataFrame Word8Clamped (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Word8Clamped ds) pn + of + Evidence -> case inferNumericFrame @Word8Clamped @(ds +: n) of + Evidence -> SomeIODataFrame df + +arrayBuffer :: IODataFrame t ds -> IO MutableArrayBuffer +arrayBuffer (IODataFrame x) = unsafeCoerce# <$> IO (arrayBuffer# x) + + +foreign import javascript unsafe "$1.length" js_abLength :: MutableArrayBuffer -> IO Int + +abDim :: Int -> Dim (ds :: [Nat]) -> MutableArrayBuffer -> IO SomeDim +abDim elS d ab = fromMaybe (SomeDim (Dn :: Dim 0)) . someDimVal . (`quot` (elS * dimVal d)) <$> js_abLength ab + +unsafeForceFixedDims :: forall ds n + . Evidence ( FixedDim (AsXDims ds +: XN 0) (ds +: n) ~ (ds +: n) + , FixedXDim (AsXDims ds +: XN 0) (ds +: n) ~ (AsXDims ds +: XN 0) + ) +unsafeForceFixedDims = unsafeCoerce# (Evidence :: Evidence ( (ds +: n) ~ (ds +: n) , (ds +: n) ~ (ds +: n) )) + +#endif diff --git a/src/Numeric/DataFrame/Mutable.hs b/src/Numeric/DataFrame/Mutable.hs deleted file mode 100644 --- a/src/Numeric/DataFrame/Mutable.hs +++ /dev/null @@ -1,223 +0,0 @@ -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE ExistentialQuantification #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE KindSignatures #-} -{-# LANGUAGE MagicHash #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeApplications #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE UnboxedTuples #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.DataFrame.Mutable --- Copyright : (c) Artem Chirkin --- License : BSD3 --- --- Maintainer : chirkin@arch.ethz.ch --- --- Interfrace to perform primitive stateful operations on mutable frames. --- ------------------------------------------------------------------------------ - -module Numeric.DataFrame.Mutable - ( MutableFrame (..), MDataFrame () - , newDataFrame#, copyDataFrame#, copyMDataFrame#, unsafeFreezeDataFrame# - , freezeDataFrame#, thawDataFrame# - , writeDataFrame#, readDataFrame# - ) where - - -import GHC.Int (Int16 (..), Int32 (..), Int64 (..), - Int8 (..)) -import GHC.Prim -import GHC.Types (Double (..), Float (..), Int (..), - Word (..)) -import GHC.Word (Word16 (..), Word32 (..), Word64 (..), - Word8 (..)) - -import Numeric.Commons -import Numeric.DataFrame.Type -import Numeric.Dimensions -import Numeric.TypeLits - - --- | Mutable DataFrame type -data MDataFrame s t (ns :: [Nat]) = MDataFrame# Int# Int# (MutableByteArray# s) - --- | Create a new mutable DataFrame. -newDataFrame# :: forall t (ns :: [Nat]) s - . ( PrimBytes t, Dimensions ns) - => State# s -> (# State# s, MDataFrame s t ns #) -newDataFrame# s0 - | elS <- elementByteSize (undefined :: t) - , I# n <- totalDim (Proxy @ns) - , (# s1, mba #) <- newByteArray# (n *# elS) s0 - = (# s1, MDataFrame# 0# n mba #) -{-# INLINE newDataFrame# #-} - --- | Copy one DataFrame into another mutable DataFrame at specified position. -copyDataFrame# :: forall t (as :: [Nat]) (bs :: [Nat]) (asbs :: [Nat]) s - . ( PrimBytes (DataFrame t as) - , ConcatList as bs asbs - , Dimensions bs - ) - => DataFrame t as -> Idx bs -> MDataFrame s t asbs -> State# s -> (# State# s, () #) -copyDataFrame# df ei (MDataFrame# offM _ arrM) s - | (# offA, lenA, arrA #) <- toBytes df - , elS <- elementByteSize df - , I# i <- fromEnum ei - = (# copyByteArray# arrA (offA *# elS) arrM ((offM +# i) *# elS) (lenA *# elS) s, () #) -{-# INLINE copyDataFrame# #-} - --- | Copy one mutable DataFrame into another mutable DataFrame at specified position. -copyMDataFrame# :: forall t (as :: [Nat]) (bs :: [Nat]) (asbs :: [Nat]) s - . ( PrimBytes t - , ConcatList as bs asbs - , Dimensions bs - ) - => MDataFrame s t as -> Idx bs -> MDataFrame s t asbs -> State# s -> (# State# s, () #) -copyMDataFrame# (MDataFrame# offA lenA arrA) ei (MDataFrame# offM _ arrM) s - | elS <- elementByteSize (undefined :: t) - , I# i <- fromEnum ei - = (# copyMutableByteArray# arrA (offA *# elS) arrM ((offM +# i) *# elS) (lenA *# elS) s, () #) -{-# INLINE copyMDataFrame# #-} - --- | Make a mutable DataFrame immutable, without copying. -unsafeFreezeDataFrame# :: forall t (ns :: [Nat]) s - . PrimBytes (DataFrame t ns) - => MDataFrame s t ns -> State# s -> (# State# s, DataFrame t ns #) -unsafeFreezeDataFrame# (MDataFrame# offM lenM arrM) s1 - | (# s2, arrA #) <- unsafeFreezeByteArray# arrM s1 - = (# s2, fromBytes (# offM, lenM, arrA #) #) -{-# INLINE unsafeFreezeDataFrame# #-} - --- | Copy content of a mutable DataFrame into a new immutable DataFrame. -freezeDataFrame# :: forall t (ns :: [Nat]) s - . PrimBytes (DataFrame t ns) - => MDataFrame s t ns -> State# s -> (# State# s, DataFrame t ns #) -freezeDataFrame# (MDataFrame# offM n arrM) s0 - | elS <- elementByteSize (undefined :: DataFrame t ns) - , (# s1, mba #) <- newByteArray# (n *# elS) s0 - , s2 <- copyMutableByteArray# arrM (offM *# elS) mba 0# (n *# elS) s1 - , (# s3, arrA #) <- unsafeFreezeByteArray# mba s2 - = (# s3, fromBytes (# 0#, n, arrA #) #) -{-# INLINE freezeDataFrame# #-} - --- | Create a new mutable DataFrame and copy content of immutable one in there. -thawDataFrame# :: forall t (ns :: [Nat]) s - . PrimBytes (DataFrame t ns) - => DataFrame t ns -> State# s -> (# State# s, MDataFrame s t ns #) -thawDataFrame# df s0 - | elS <- elementByteSize (undefined :: DataFrame t ns) - , (# offA, n, arrA #) <- toBytes df - , (# s1, arrM #) <- newByteArray# (n *# elS) s0 - , s2 <- copyByteArray# arrA (offA *# elS) arrM 0# (n *# elS) s1 - = (# s2, MDataFrame# 0# n arrM #) -{-# INLINE thawDataFrame# #-} - --- | Write a single element at the specified index -writeDataFrame# :: forall t (ns :: [Nat]) s - . ( MutableFrame t ns, Dimensions ns ) - => MDataFrame s t ns -> Idx ns -> t -> State# s -> (# State# s, () #) -writeDataFrame# mdf ei x s | I# i <- fromEnum ei = (# writeDataFrameOff# mdf i x s, () #) -{-# INLINE writeDataFrame# #-} - --- | Read a single element at the specified index -readDataFrame# :: forall t (ns :: [Nat]) s - . ( MutableFrame t ns, Dimensions ns ) - => MDataFrame s t ns -> Idx ns -> State# s -> (# State# s, t #) -readDataFrame# mdf ei | I# i <- fromEnum ei = readDataFrameOff# mdf i -{-# INLINE readDataFrame# #-} - -class MutableFrame t (ns :: [Nat]) where - -- | Write a single element at the specified element offset - writeDataFrameOff# :: MDataFrame s t ns -> Int# -> t -> State# s -> State# s - -- | Read a single element at the specified element offset - readDataFrameOff# :: MDataFrame s t ns -> Int# -> State# s -> (# State# s, t #) - -instance MutableFrame Float (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (F# x) = writeFloatArray# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readFloatArray# arr (off +# i) s0 = (# s1, F# r #) - {-# INLINE readDataFrameOff# #-} - -instance MutableFrame Double (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (D# x) = writeDoubleArray# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readDoubleArray# arr (off +# i) s0 = (# s1, D# r #) - {-# INLINE readDataFrameOff# #-} - -instance MutableFrame Int (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (I# x) = writeIntArray# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readIntArray# arr (off +# i) s0 = (# s1, I# r #) - {-# INLINE readDataFrameOff# #-} - -instance MutableFrame Int8 (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (I8# x) = writeInt8Array# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readInt8Array# arr (off +# i) s0 = (# s1, I8# r #) - {-# INLINE readDataFrameOff# #-} - -instance MutableFrame Int16 (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (I16# x) = writeInt16Array# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readInt16Array# arr (off +# i) s0 = (# s1, I16# r #) - {-# INLINE readDataFrameOff# #-} - -instance MutableFrame Int32 (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (I32# x) = writeInt32Array# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readInt32Array# arr (off +# i) s0 = (# s1, I32# r #) - {-# INLINE readDataFrameOff# #-} - -instance MutableFrame Int64 (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (I64# x) = writeInt64Array# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readInt64Array# arr (off +# i) s0 = (# s1, I64# r #) - {-# INLINE readDataFrameOff# #-} - - -instance MutableFrame Word (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (W# x) = writeWordArray# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readWordArray# arr (off +# i) s0 = (# s1, W# r #) - {-# INLINE readDataFrameOff# #-} - -instance MutableFrame Word8 (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (W8# x) = writeWord8Array# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readWord8Array# arr (off +# i) s0 = (# s1, W8# r #) - {-# INLINE readDataFrameOff# #-} - -instance MutableFrame Word16 (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (W16# x) = writeWord16Array# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readWord16Array# arr (off +# i) s0 = (# s1, W16# r #) - {-# INLINE readDataFrameOff# #-} - -instance MutableFrame Word32 (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (W32# x) = writeWord32Array# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readWord32Array# arr (off +# i) s0 = (# s1, W32# r #) - {-# INLINE readDataFrameOff# #-} - -instance MutableFrame Word64 (ns :: [Nat]) where - writeDataFrameOff# (MDataFrame# off _ arr) i (W64# x) = writeWord64Array# arr (off +# i) x - {-# INLINE writeDataFrameOff# #-} - readDataFrameOff# (MDataFrame# off _ arr) i s0 - | (# s1, r #) <- readWord64Array# arr (off +# i) s0 = (# s1, W64# r #) - {-# INLINE readDataFrameOff# #-} diff --git a/src/Numeric/DataFrame/ST.hs b/src/Numeric/DataFrame/ST.hs --- a/src/Numeric/DataFrame/ST.hs +++ b/src/Numeric/DataFrame/ST.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE ExistentialQuantification #-} {-# LANGUAGE FlexibleContexts #-} @@ -9,6 +10,7 @@ {-# LANGUAGE TypeApplications #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE UnboxedTuples #-} +{-# LANGUAGE TypeOperators #-} ----------------------------------------------------------------------------- -- | -- Module : Numeric.DataFrame.ST @@ -22,18 +24,41 @@ ----------------------------------------------------------------------------- module Numeric.DataFrame.ST - ( MutableFrame (), STDataFrame () + ( +#ifdef ghcjs_HOST_OS + MutableFrame (), STDataFrame (..), MDataFrame (..), MutableArrayT (..) +#else + MutableFrame (), STDataFrame () +#endif + , SomeSTDataFrame (..) , newDataFrame, copyDataFrame, copyMutableDataFrame , unsafeFreezeDataFrame , freezeDataFrame, thawDataFrame , writeDataFrame, readDataFrame , writeDataFrameOff, readDataFrameOff +#ifdef ghcjs_HOST_OS + -- * JavaScript-specific functions + , STArrayBuffer + , newArrayBuffer, arrayBuffer, viewFloatArray, viewDoubleArray + , viewIntArray, viewInt32Array, viewInt16Array, viewInt8Array + , viewWordArray, viewWord32Array, viewWord16Array, viewWord8Array, viewWord8ClampedArray +#endif ) where import GHC.Types (Int (..)) import GHC.ST (ST(..)) +#ifdef ghcjs_HOST_OS +import Numeric.Array.Family hiding (Scalar) +import JavaScript.TypedArray.ArrayBuffer.ST +import GHC.Prim +import Data.Int +import Data.Word +import Data.Maybe +import GHCJS.Types +import Numeric.DataFrame.Inference +#endif import Numeric.Commons import Numeric.DataFrame.Type import Numeric.DataFrame.Mutable @@ -44,32 +69,52 @@ -- | Mutable DataFrame that lives in ST. -- Internal representation is always a ByteArray. newtype STDataFrame s t (ns :: [Nat]) = STDataFrame (MDataFrame s t (ns :: [Nat])) - +#ifdef ghcjs_HOST_OS +instance IsJSVal (STDataFrame s t ds) +#endif +-- | Mutable DataFrame of unknown dimensionality +data SomeSTDataFrame s t (xns :: [XNat]) + = forall (ns :: [Nat]) + . ( FixedDim xns ns ~ ns + , FixedXDim xns ns ~ xns + , NumericFrame t ns + ) + => SomeSTDataFrame (STDataFrame s t ns) -- | Create a new mutable DataFrame. newDataFrame :: forall t (ns :: [Nat]) s +#ifdef ghcjs_HOST_OS + . ( ElemTypeInference t, Dimensions ns) +#else . ( PrimBytes t, Dimensions ns) +#endif => ST s (STDataFrame s t ns) newDataFrame = STDataFrame <$> ST (newDataFrame# @t @ns) {-# INLINE newDataFrame #-} -- | Copy one DataFrame into another mutable DataFrame at specified position. -copyDataFrame :: forall t (as :: [Nat]) (bs :: [Nat]) (asbs :: [Nat]) s - . ( PrimBytes (DataFrame t as) - , ConcatList as bs asbs - , Dimensions bs - ) - => DataFrame t as -> Idx bs -> STDataFrame s t asbs -> ST s () +copyDataFrame :: forall t (as :: [Nat]) (b' :: Nat) (b :: Nat) (bs :: [Nat]) (asbs :: [Nat]) s + . ( ConcatList as (b :+ bs) asbs, Dimensions (b :+ bs) +#ifdef ghcjs_HOST_OS + , ArraySizeInference (as +: b'), Dimensions as +#else + , PrimBytes (DataFrame t (as +: b')) +#endif + ) + => DataFrame t (as +: b') -> Idx (b :+ bs) -> STDataFrame s t asbs -> ST s () copyDataFrame df ei (STDataFrame mdf) = ST (copyDataFrame# df ei mdf) {-# INLINE copyDataFrame #-} -- | Copy one mutable DataFrame into another mutable DataFrame at specified position. -copyMutableDataFrame :: forall t (as :: [Nat]) (bs :: [Nat]) (asbs :: [Nat]) s +copyMutableDataFrame :: forall t (as :: [Nat]) (b' :: Nat) (b :: Nat) (bs :: [Nat]) (asbs :: [Nat]) s . ( PrimBytes t - , ConcatList as bs asbs - , Dimensions bs + , ConcatList as (b :+ bs) asbs + , Dimensions (b :+ bs) +#ifdef ghcjs_HOST_OS + , Dimensions as +#endif ) - => STDataFrame s t as -> Idx bs -> STDataFrame s t asbs -> ST s () + => STDataFrame s t (as +: b') -> Idx (b :+ bs) -> STDataFrame s t asbs -> ST s () copyMutableDataFrame (STDataFrame mdfA) ei (STDataFrame mdfB) = ST (copyMDataFrame# mdfA ei mdfB) {-# INLINE copyMutableDataFrame #-} @@ -77,7 +122,11 @@ -- | Make a mutable DataFrame immutable, without copying. unsafeFreezeDataFrame :: forall t (ns :: [Nat]) s - . PrimBytes (DataFrame t ns) +#ifdef ghcjs_HOST_OS + . (MutableFrame t ns, ArraySizeInference ns) +#else + . PrimBytes (DataFrame t ns) +#endif => STDataFrame s t ns -> ST s (DataFrame t ns) unsafeFreezeDataFrame (STDataFrame mdf) = ST (unsafeFreezeDataFrame# mdf) {-# INLINE unsafeFreezeDataFrame #-} @@ -85,14 +134,22 @@ -- | Copy content of a mutable DataFrame into a new immutable DataFrame. freezeDataFrame :: forall t (ns :: [Nat]) s +#ifdef ghcjs_HOST_OS + . (MutableFrame t ns, ArraySizeInference ns) +#else . PrimBytes (DataFrame t ns) +#endif => STDataFrame s t ns -> ST s (DataFrame t ns) freezeDataFrame (STDataFrame mdf) = ST (freezeDataFrame# mdf) {-# INLINE freezeDataFrame #-} -- | Create a new mutable DataFrame and copy content of immutable one in there. thawDataFrame :: forall t (ns :: [Nat]) s - . PrimBytes (DataFrame t ns) +#ifdef ghcjs_HOST_OS + . (MutableFrame t ns, ArrayInstanceInference t ns) +#else + . PrimBytes (DataFrame t ns) +#endif => DataFrame t ns -> ST s (STDataFrame s t ns) thawDataFrame df = STDataFrame <$> ST (thawDataFrame# df) {-# INLINE thawDataFrame #-} @@ -128,3 +185,170 @@ => STDataFrame s t ns -> Int -> ST s (Scalar t) readDataFrameOff (STDataFrame mdf) (I# i) = scalar <$> ST (readDataFrameOff# mdf i) {-# INLINE readDataFrameOff #-} + + +#ifdef ghcjs_HOST_OS + +newArrayBuffer :: Int -> ST s (STArrayBuffer s) +newArrayBuffer n = unsafeCoerce# <$> ST (newArrayBuffer# n) + +viewFloatArray :: forall s ds + . ( Dimensions ds, ArraySizeInference ds) + => STArrayBuffer s -> ST s (SomeSTDataFrame s Float (AsXDims ds +: XN 0)) +viewFloatArray ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Float))) (dim @ds) ab + df <- fmap STDataFrame . ST $ viewFloatArray# (jsval ab) :: ST s (STDataFrame s Float (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Float ds) pn + of + Evidence -> case inferNumericFrame @Float @(ds +: n) of + Evidence -> SomeSTDataFrame df + +viewDoubleArray :: forall s ds + . ( Dimensions ds, ArraySizeInference ds) + => STArrayBuffer s -> ST s (SomeSTDataFrame s Double (AsXDims ds +: XN 0)) +viewDoubleArray ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Double))) (dim @ds) ab + df <- fmap STDataFrame . ST $ viewDoubleArray# (jsval ab) :: ST s (STDataFrame s Double (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Double ds) pn + of + Evidence -> case inferNumericFrame @Double @(ds +: n) of + Evidence -> SomeSTDataFrame df + +viewIntArray :: forall s ds + . ( Dimensions ds, ArraySizeInference ds) + => STArrayBuffer s -> ST s (SomeSTDataFrame s Int (AsXDims ds +: XN 0)) +viewIntArray ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Int))) (dim @ds) ab + df <- fmap STDataFrame . ST $ viewIntArray# (jsval ab) :: ST s (STDataFrame s Int (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Int ds) pn + of + Evidence -> case inferNumericFrame @Int @(ds +: n) of + Evidence -> SomeSTDataFrame df + +viewInt32Array :: forall s ds + . ( Dimensions ds, ArraySizeInference ds) + => STArrayBuffer s -> ST s (SomeSTDataFrame s Int32 (AsXDims ds +: XN 0)) +viewInt32Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Int32))) (dim @ds) ab + df <- fmap STDataFrame . ST $ viewInt32Array# (jsval ab) :: ST s (STDataFrame s Int32 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Int32 ds) pn + of + Evidence -> case inferNumericFrame @Int32 @(ds +: n) of + Evidence -> SomeSTDataFrame df + +viewInt16Array :: forall s ds + . ( Dimensions ds, ArraySizeInference ds) + => STArrayBuffer s -> ST s (SomeSTDataFrame s Int16 (AsXDims ds +: XN 0)) +viewInt16Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Int16))) (dim @ds) ab + df <- fmap STDataFrame . ST $ viewInt16Array# (jsval ab) :: ST s (STDataFrame s Int16 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Int16 ds) pn + of + Evidence -> case inferNumericFrame @Int16 @(ds +: n) of + Evidence -> SomeSTDataFrame df + +viewInt8Array :: forall s ds + . ( Dimensions ds, ArraySizeInference ds) + => STArrayBuffer s -> ST s (SomeSTDataFrame s Int8 (AsXDims ds +: XN 0)) +viewInt8Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Int8))) (dim @ds) ab + df <- fmap STDataFrame . ST $ viewInt8Array# (jsval ab) :: ST s (STDataFrame s Int8 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Int8 ds) pn + of + Evidence -> case inferNumericFrame @Int8 @(ds +: n) of + Evidence -> SomeSTDataFrame df + +viewWordArray :: forall s ds + . ( Dimensions ds, ArraySizeInference ds) + => STArrayBuffer s -> ST s (SomeSTDataFrame s Word (AsXDims ds +: XN 0)) +viewWordArray ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Word))) (dim @ds) ab + df <- fmap STDataFrame . ST $ viewWordArray# (jsval ab) :: ST s (STDataFrame s Word (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Word ds) pn + of + Evidence -> case inferNumericFrame @Word @(ds +: n) of + Evidence -> SomeSTDataFrame df + +viewWord32Array :: forall s ds + . ( Dimensions ds, ArraySizeInference ds) + => STArrayBuffer s -> ST s (SomeSTDataFrame s Word32 (AsXDims ds +: XN 0)) +viewWord32Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Word32))) (dim @ds) ab + df <- fmap STDataFrame . ST $ viewWord32Array# (jsval ab) :: ST s (STDataFrame s Word32 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Word32 ds) pn + of + Evidence -> case inferNumericFrame @Word32 @(ds +: n) of + Evidence -> SomeSTDataFrame df + +viewWord16Array :: forall s ds + . ( Dimensions ds, ArraySizeInference ds) + => STArrayBuffer s -> ST s (SomeSTDataFrame s Word16 (AsXDims ds +: XN 0)) +viewWord16Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Word16))) (dim @ds) ab + df <- fmap STDataFrame . ST $ viewWord16Array# (jsval ab) :: ST s (STDataFrame s Word16 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Word16 ds) pn + of + Evidence -> case inferNumericFrame @Word16 @(ds +: n) of + Evidence -> SomeSTDataFrame df + +viewWord8Array :: forall s ds + . ( Dimensions ds, ArraySizeInference ds) + => STArrayBuffer s -> ST s (SomeSTDataFrame s Word8 (AsXDims ds +: XN 0)) +viewWord8Array ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Word8))) (dim @ds) ab + df <- fmap STDataFrame . ST $ viewWord8Array# (jsval ab) :: ST s (STDataFrame s Word8 (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Word8 ds) pn + of + Evidence -> case inferNumericFrame @Word8 @(ds +: n) of + Evidence -> SomeSTDataFrame df + +viewWord8ClampedArray :: forall s ds + . ( Dimensions ds, ArraySizeInference ds) + => STArrayBuffer s -> ST s (SomeSTDataFrame s Word8Clamped (AsXDims ds +: XN 0)) +viewWord8ClampedArray ab = do + SomeDim (pn@Dn :: Dim (n :: Nat)) <- abDim (I# (byteSize (undefined :: Word8Clamped))) (dim @ds) ab + df <- fmap STDataFrame . ST $ viewWord8ClampedArray# (jsval ab) :: ST s (STDataFrame s Word8Clamped (ds +: n)) + return $ case unsafeForceFixedDims @ds @n + `sumEvs` inferSnocDimensions @ds @n + `sumEvs` inferSnocArrayInstance (undefined :: DataFrame Word8Clamped ds) pn + of + Evidence -> case inferNumericFrame @Word8Clamped @(ds +: n) of + Evidence -> SomeSTDataFrame df + +arrayBuffer :: STDataFrame s t ds -> ST s (STArrayBuffer s) +arrayBuffer (STDataFrame x) = unsafeCoerce# <$> ST (arrayBuffer# x) + + +foreign import javascript unsafe "$1.length" js_abLength :: STArrayBuffer s -> Int + +abDim :: Int -> Dim (ds :: [Nat]) -> STArrayBuffer s -> ST s SomeDim +abDim elS d ab = fromMaybe (SomeDim (Dn :: Dim 0)) . someDimVal . (`quot` (elS * dimVal d)) <$> pure (js_abLength ab) +{-# NOINLINE abDim #-} + +unsafeForceFixedDims :: forall ds n + . Evidence ( FixedDim (AsXDims ds +: XN 0) (ds +: n) ~ (ds +: n) + , FixedXDim (AsXDims ds +: XN 0) (ds +: n) ~ (AsXDims ds +: XN 0) + ) +unsafeForceFixedDims = unsafeCoerce# (Evidence :: Evidence ( (ds +: n) ~ (ds +: n) , (ds +: n) ~ (ds +: n) )) + +#endif diff --git a/src/Numeric/DataFrame/SubSpace.hs b/src/Numeric/DataFrame/SubSpace.hs --- a/src/Numeric/DataFrame/SubSpace.hs +++ b/src/Numeric/DataFrame/SubSpace.hs @@ -69,7 +69,7 @@ | asbs as -> bs, asbs bs -> as, as bs -> asbs where -- | Unsafely get a sub-dataframe by its primitive element subset. -- The offset is not checked to be aligned to the space structure or for bounds. - -- Arguments are zero-based element offset and element size (aka `totalDim` of sub dataframe) + -- Arguments are zero-based primitive element offset and subset ("as" element) size (aka `totalDim` of sub dataframe) -- -- Normal indexing can be expressed in terms of `indexOffset#`: -- @@ -152,13 +152,13 @@ ewfoldMap :: forall t (as :: [Nat]) (bs :: [Nat]) (asbs :: [Nat]) m . (Monoid m, SubSpace t as bs asbs) => (DataFrame t as -> m) -> DataFrame t asbs -> m -ewfoldMap f = ewfoldl (\b -> mappend b . f) mempty +ewfoldMap f = ewfoldl (\m b -> m `seq` (mappend m $! f b)) mempty {-# INLINE ewfoldMap #-} iwfoldMap :: forall t (as :: [Nat]) (bs :: [Nat]) (asbs :: [Nat]) m . ( Monoid m, SubSpace t as bs asbs) => (Idx bs -> DataFrame t as -> m) -> DataFrame t asbs -> m -iwfoldMap f = iwfoldl (\i b -> mappend b . f i) mempty +iwfoldMap f = iwfoldl (\i m b -> m `seq` (mappend m $! f i b)) mempty {-# INLINE iwfoldMap #-} @@ -218,34 +218,44 @@ #endif {-# INLINE indexOffset# #-} + ewmap :: forall s (as' :: [Nat]) (asbs' :: [Nat]) + . SubSpace s as' bs asbs' + => (DataFrame s as' -> DataFrame t as) + -> DataFrame s asbs' -> DataFrame t asbs ewmap f df | elS <- elementByteSize (undefined :: DataFrame t asbs) , I# lenBS <- totalDim (Proxy @bs) , I# lenAS <- totalDim (Proxy @as) + , I# lenAS' <- totalDim (Proxy @as') , lenASB <- lenAS *# elS = case runRW# ( \s0 -> case newByteArray# (lenAS *# lenBS *# elS) s0 of (# s1, marr #) -> case overDimOff_# (dim @bs) - ( \pos s -> case toBytes $ f (indexOffset# pos lenAS df) of - (# offX, _, arrX #) -> copyByteArray# arrX (offX *# elS) marr (pos *# elS) lenASB s - ) 0# lenAS s1 of + ( \pos s -> case toBytes $ f (indexOffset# (pos *# lenAS') lenAS' df) of + (# offX, _, arrX #) -> copyByteArray# arrX (offX *# elS) marr (pos *# lenASB) lenASB s + ) 0# 1# s1 of s2 -> unsafeFreezeByteArray# marr s2 ) of (# _, r #) -> fromBytes (# 0#, lenAS *# lenBS, r #) {-# INLINE ewmap #-} + iwmap :: forall s (as' :: [Nat]) (asbs' :: [Nat]) + . SubSpace s as' bs asbs' + => (Idx bs -> DataFrame s as' -> DataFrame t as) + -> DataFrame s asbs' -> DataFrame t asbs iwmap f df | elS <- elementByteSize (undefined :: DataFrame t asbs) , I# lenBS <- totalDim (Proxy @bs) , I# lenAS <- totalDim (Proxy @as) + , I# lenAS' <- totalDim (Proxy @as') , lenASB <- lenAS *# elS = case runRW# ( \s0 -> case newByteArray# (lenAS *# lenBS *# elS) s0 of (# s1, marr #) -> case overDim_# (dim @bs) - ( \i pos s -> case toBytes $ f i (indexOffset# pos lenAS df) of - (# offX, _, arrX #) -> copyByteArray# arrX (offX *# elS) marr (pos *# elS) lenASB s - ) 0# lenAS s1 of + ( \i pos s -> case toBytes $ f i (indexOffset# (pos *# lenAS') lenAS' df) of + (# offX, _, arrX #) -> copyByteArray# arrX (offX *# elS) marr (pos *# lenASB) lenASB s + ) 0# 1# s1 of s2 -> unsafeFreezeByteArray# marr s2 ) of (# _, r #) -> fromBytes (# 0#, lenAS *# lenBS, r #) diff --git a/src/Numeric/DataFrame/Type.hs b/src/Numeric/DataFrame/Type.hs --- a/src/Numeric/DataFrame/Type.hs +++ b/src/Numeric/DataFrame/Type.hs @@ -38,10 +38,17 @@ , FPFRame, IntegralFrame, NumericVariantFrame, CommonOpFrame ) where +#include "MachDeps.h" import Data.Int (Int16, Int32, Int64, Int8) import Data.Word (Word16, Word32, Word64, Word8) import Foreign.Storable (Storable (..)) -import GHC.Exts (Int (..), Ptr (..), Float#, Double#, Int#, Word#) +import GHC.Exts (Int (..), Ptr (..), Float#, Double#, Int#, Word# +#ifndef ghcjs_HOST_OS +#if WORD_SIZE_IN_BITS < 64 + , Int64#, Word64# +#endif +#endif + ) import GHC.Prim (copyAddrToByteArray#, copyByteArrayToAddr#, newByteArray#, plusAddr#, quotInt#, @@ -210,8 +217,12 @@ type instance ElemPrim (DataFrame Int16 ds) = Int# type instance ElemPrim (DataFrame Int32 ds) = Int# #ifndef ghcjs_HOST_OS +#if WORD_SIZE_IN_BITS < 64 +type instance ElemPrim (DataFrame Int64 ds) = Int64# +#else type instance ElemPrim (DataFrame Int64 ds) = Int# #endif +#endif type instance ElemPrim (DataFrame Word ds) = Word# type instance ElemPrim (DataFrame Word8 ds) = Word# type instance ElemPrim (DataFrame Word16 ds) = Word# @@ -219,8 +230,12 @@ #ifdef ghcjs_HOST_OS type instance ElemPrim (DataFrame Word8Clamped ds) = Int# #else +#if WORD_SIZE_IN_BITS < 64 +type instance ElemPrim (DataFrame Word64 ds) = Word64# +#else type instance ElemPrim (DataFrame Word64 ds) = Word# #endif +#endif deriving instance ( PrimBytes (Array Float ds) , ElemPrim (Array Float ds) ~ Float# , ElemRep (Array Float ds) ~ 'FloatRep) => PrimBytes (DataFrame Float ds) @@ -241,9 +256,15 @@ , ElemRep (Array Int32 ds) ~ 'IntRep) => PrimBytes (DataFrame Int32 ds) #ifndef ghcjs_HOST_OS deriving instance ( PrimBytes (Array Int64 ds) +#if WORD_SIZE_IN_BITS < 64 + , ElemPrim (Array Int64 ds) ~ Int64# + , ElemRep (Array Int64 ds) ~ 'Int64Rep +#else , ElemPrim (Array Int64 ds) ~ Int# - , ElemRep (Array Int64 ds) ~ 'IntRep) => PrimBytes (DataFrame Int64 ds) + , ElemRep (Array Int64 ds) ~ 'IntRep #endif + ) => PrimBytes (DataFrame Int64 ds) +#endif deriving instance ( PrimBytes (Array Word ds) , ElemPrim (Array Word ds) ~ Word# , ElemRep (Array Word ds) ~ 'WordRep) => PrimBytes (DataFrame Word ds) @@ -262,8 +283,14 @@ , ElemRep (Array Word8Clamped ds) ~ 'IntRep) => PrimBytes (DataFrame Word8Clamped ds) #else deriving instance ( PrimBytes (Array Word64 ds) +#if WORD_SIZE_IN_BITS < 64 + , ElemPrim (Array Word64 ds) ~ Word64# + , ElemRep (Array Word64 ds) ~ 'Word64Rep +#else , ElemPrim (Array Word64 ds) ~ Word# - , ElemRep (Array Word64 ds) ~ 'WordRep) => PrimBytes (DataFrame Word64 ds) + , ElemRep (Array Word64 ds) ~ 'WordRep +#endif + ) => PrimBytes (DataFrame Word64 ds) #endif diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs --- a/src/Numeric/Matrix.hs +++ b/src/Numeric/Matrix.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE ExistentialQuantification #-} {-# LANGUAGE FlexibleContexts #-} @@ -20,39 +21,40 @@ ( MatrixCalculus (..) , SquareMatrixCalculus (..) , MatrixInverse (..) + , HomTransform4 (..) , Matrix , Mat22f, Mat23f, Mat24f , Mat32f, Mat33f, Mat34f , Mat42f, Mat43f, Mat44f + , Mat22d, Mat23d, Mat24d + , Mat32d, Mat33d, Mat34d + , Mat42d, Mat43d, Mat44d , mat22, mat33, mat44 , (%*) ) where + + +#ifdef ghcjs_HOST_OS +import Numeric.Array.Family (ElemTypeInference) +#endif + import GHC.Types (Type) import Numeric.Commons import Numeric.DataFrame.Contraction ((%*)) import Numeric.DataFrame.Shape -import Numeric.Dimensions (Nat) -import Numeric.Matrix.Type +import Numeric.Dimensions (Nat, Idx (..)) +import Numeric.Matrix.Class +import Numeric.Matrix.Mat44d () +import Numeric.Matrix.Mat44f () import Numeric.Vector import Control.Monad.ST import Numeric.DataFrame.ST --- 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 - -- | Compose a 2x2D matrix mat22 :: ( PrimBytes (Vector t 2) , PrimBytes (Matrix t 2 2) @@ -61,30 +63,40 @@ mat22 = (<::>) -- | Compose a 3x3D matrix -mat33 :: ( PrimBytes t +mat33 :: ( +#ifdef ghcjs_HOST_OS + ElemTypeInference t, MutableFrame t '[3,3] +#else + PrimBytes t , PrimBytes (Vector t 3) , PrimBytes (Matrix t 3 3) +#endif ) => Vector t 3 -> Vector t 3 -> Vector t 3 -> Matrix t 3 3 mat33 a b c = runST $ do mmat <- newDataFrame - copyDataFrame a 1 mmat - copyDataFrame b 2 mmat - copyDataFrame c 3 mmat + copyDataFrame a (1:!1:!Z) mmat + copyDataFrame b (1:!2:!Z) mmat + copyDataFrame c (1:!3:!Z) mmat unsafeFreezeDataFrame mmat -- | Compose a 4x4D matrix mat44 :: forall (t :: Type) - . ( PrimBytes t + . ( +#ifdef ghcjs_HOST_OS + ElemTypeInference t, MutableFrame t '[4,4] +#else + PrimBytes t , PrimBytes (Vector t (4 :: Nat)) , PrimBytes (Matrix t (4 :: Nat) (4 :: Nat)) +#endif ) => Vector t (4 :: Nat) -> Vector t (4 :: Nat) -> Vector t (4 :: Nat) -> Vector t (4 :: Nat) -> Matrix t (4 :: Nat) (4 :: Nat) mat44 a b c d = runST $ do mmat <- newDataFrame - copyDataFrame a 1 mmat - copyDataFrame b 2 mmat - copyDataFrame c 3 mmat - copyDataFrame d 4 mmat + copyDataFrame a (1:!1:!Z) mmat + copyDataFrame b (1:!2:!Z) mmat + copyDataFrame c (1:!3:!Z) mmat + copyDataFrame d (1:!4:!Z) mmat unsafeFreezeDataFrame mmat diff --git a/src/Numeric/Matrix/Class.hs b/src/Numeric/Matrix/Class.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Matrix/Class.hs @@ -0,0 +1,114 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE TypeFamilies #-} +module Numeric.Matrix.Class + ( MatrixCalculus (..) + , SquareMatrixCalculus (..) + , MatrixInverse (..) + , Matrix + , HomTransform4 (..) + , Mat22f, Mat23f, Mat24f + , Mat32f, Mat33f, Mat34f + , Mat42f, Mat43f, Mat44f + , Mat22d, Mat23d, Mat24d + , Mat32d, Mat33d, Mat34d + , Mat42d, Mat43d, Mat44d + ) where + +import Numeric.Commons +import Numeric.DataFrame.Type +import Numeric.Dimensions (Nat) +import Numeric.Scalar +import Numeric.Vector + +-- | Alias for DataFrames of rank 2 +type Matrix t (n :: Nat) (m :: Nat) = DataFrame t '[n,m] + +class MatrixCalculus t (n :: Nat) (m :: Nat) where + -- | Transpose Mat + transpose :: (MatrixCalculus t m n, PrimBytes (Matrix t m n)) => Matrix t n m -> Matrix t m n + + +class SquareMatrixCalculus 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 + -- | Determinant of Mat + det :: Matrix t n n -> Scalar t + -- | Sum of diagonal elements + trace :: Matrix t n n -> Scalar t + +class MatrixInverse t (n :: Nat) where + -- | Matrix inverse + inverse :: DataFrame t '[n,n] -> DataFrame t '[n,n] + + +-- | 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 + translate4 :: Vector t 4 -> Matrix t 4 4 + -- | 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 + -- | Rotation matrix for a rotation around the Y axis, angle is given in radians. + rotateY :: t -> Matrix t 4 4 + -- | Rotation matrix for a rotation around the Z axis, angle is given in radians. + 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 + -- | 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 diff --git a/src/Numeric/Matrix/Type.hs b/src/Numeric/Matrix/Type.hs deleted file mode 100644 --- a/src/Numeric/Matrix/Type.hs +++ /dev/null @@ -1,49 +0,0 @@ -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE KindSignatures #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE TypeFamilies #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.Matrix.Type --- Copyright : (c) Artem Chirkin --- License : BSD3 --- --- Maintainer : chirkin@arch.ethz.ch --- --- ------------------------------------------------------------------------------ - -module Numeric.Matrix.Type - ( MatrixCalculus (..) - , SquareMatrixCalculus (..) - , MatrixInverse (..) - , Matrix - ) where - -import Numeric.Commons -import Numeric.DataFrame.Type -import Numeric.Dimensions (Nat) -import Numeric.Scalar - --- | Alias for DataFrames of rank 2 -type Matrix t (n :: Nat) (m :: Nat) = DataFrame t '[n,m] - -class MatrixCalculus t (n :: Nat) (m :: Nat) where - -- | Transpose Mat - transpose :: (MatrixCalculus t m n, PrimBytes (Matrix t m n)) => Matrix t n m -> Matrix t m n - - -class SquareMatrixCalculus 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 - -- | Determinant of Mat - det :: Matrix t n n -> Scalar t - -- | Sum of diagonal elements - trace :: Matrix t n n -> Scalar t - -class MatrixInverse t (n :: Nat) where - -- | Matrix inverse - inverse :: DataFrame t '[n,n] -> DataFrame t '[n,n] diff --git a/src/Numeric/Quaternion.hs b/src/Numeric/Quaternion.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Quaternion.hs @@ -0,0 +1,35 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Quaternion +-- Copyright : (c) Artem Chirkin +-- License : BSD3 +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- Quaternion operations implemented for Floats and Doubles. +-- +-- The types `QDouble` and `QFloat` have the same representation as corresponding `Vector t 4`. +-- This means, you can do a cheap conversion between the types. +-- +-- However, arithmetic instances, such as Num and Floating are implemented in substentially different ways. +-- For example, fromInteger fills a vector fully but sets only real part to a quaternion: +-- +-- >>> 1 = vec4 1 1 1 1 +-- >>> 1 = packQ 0 0 0 1 +-- +-- All other numeric operations for vectors are element-wise, but for quaternions I have implemented +-- the actual quaternion math. +-- Some of the operations (such as trigonometric operations) are ambiguous for quaternions; +-- the general rules I follow: +-- +-- 1. Preserve imaginary vector axis same if possible; +-- 2. If both @+q@ and @-q@ are possible, prefer real value positive (@re q >= 0@). +----------------------------------------------------------------------------- +module Numeric.Quaternion + ( Quaternion (..) + , QDouble, QFloat + ) where + +import Numeric.Quaternion.Class +import Numeric.Quaternion.QFloat +import Numeric.Quaternion.QDouble diff --git a/src/Numeric/Quaternion/Class.hs b/src/Numeric/Quaternion/Class.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Quaternion/Class.hs @@ -0,0 +1,85 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE TypeFamilies #-} +module Numeric.Quaternion.Class + ( Quaternion (..) + ) where + +import Numeric.Matrix (Matrix) +import Numeric.Vector (Vector) + +-- | Quaternion operations +class Quaternion t where + -- | Quaternion data type. The ordering of coordinates is @x y z w@, + -- where @w@ is an argument, and @x y z@ are components of a 3D vector + data Quater t + -- | Set the quaternion in format (x,y,z,w) + packQ :: t -> t -> t -> t -> Quater t + -- | Get the values of the quaternion in format (x,y,z,w) + unpackQ :: Quater t -> (t,t,t,t) + -- | Set the quaternion from 3D axis vector and argument + fromVecNum :: Vector t 3 -> t -> Quater t + -- | Set the quaternion from 4D vector in format (x,y,z,w) + fromVec4 :: Vector t 4 -> Quater t + -- | Transform the quaternion to 4D vector in format (x,y,z,w) + toVec4 :: Quater t -> Vector t 4 + -- | Get scalar square of the quaternion. + -- + -- >> realToFrac (square q) == q * conjugate q + square :: Quater t -> t + -- | Imagine part of quaternion (orientation vector) + im :: Quater t -> Quater t + -- | Real part of the quaternion + re :: Quater t -> Quater t + -- | Get imagenery 3D vector of the quaternion + imVec :: Quater t -> Vector t 3 + -- | Real part of the quaternion into number + taker :: Quater t -> t + -- | i-th component + takei :: Quater t -> t + -- | j-th component + takej :: Quater t -> t + -- | k-th component + takek :: Quater t -> t + -- | Conjugate quaternion (negate imaginary part) + conjugate :: Quater t -> Quater t + -- | Rotates and scales vector in 3D using quaternion. + -- Let @q = (cos (a\/2), sin (a\/2) * v)@; then the rotation angle is @a@, and the axis of rotation is @v@. + -- Scaling is proportional to @|v|^2@. + -- + -- >>> rotScale q x == q * x * (conjugate q) + rotScale :: Quater t -> Vector t 3 -> Vector t 3 + -- | Creates a quaternion @q@ from two vectors @a@ and @b@. + -- @rotScale q a == b@ + getRotScale :: Vector t 3 -> Vector t 3 -> Quater t + -- | Creates a rotation versor from an axis vector and an angle in radians. + -- Result is always a unit quaternion (versor). + -- If the argument vector is zero, then result is a real unit quaternion. + axisRotation :: Vector t 3 -> t -> Quater t + -- | Quaternion rotation angle + -- + -- >>> q /= 0 ==> axisRotation (imVec q) (qArg q) == signum q + qArg :: Quater t -> t + -- | Create a quaternion from a rotation matrix. + -- Note, that rotations of `q` and `-q` are equivalent, there result of this function may be + -- ambiguious. I decided to force its real part be positive: + -- + -- >>> taker (fromMatrix33 m) >= 0 + fromMatrix33 :: Matrix t 3 3 -> Quater t + -- | Create a quaternion from a homogenious coordinates trasform matrix. + -- Ignores matrix translation transform. + -- Note, that rotations of `q` and `-q` are equivalent, there result of this function may be + -- ambiguious. I decided to force its real part be positive: + -- + -- >>> taker (fromMatrix44 m) >= 0 + fromMatrix44 :: Matrix t 4 4 -> Quater t + -- | Create a rotation matrix from a quaternion. + -- Note, that rotations of `q` and `-q` are equivalent, so the following property holds: + -- + -- >>> toMatrix33 q == toMatrix33 (-q) + toMatrix33 :: Quater t -> Matrix t 3 3 + -- | Create a homogenious coordinates trasform matrix from a quaternion. + -- Translation of the output matrix is zero. + -- Note, that rotations of `q` and `-q` are equivalent, so the following property holds: + -- + -- >>> toMatrix44 q == toMatrix44 (-q) + toMatrix44 :: Quater t -> Matrix t 4 4 diff --git a/src/Numeric/Vector.hs b/src/Numeric/Vector.hs --- a/src/Numeric/Vector.hs +++ b/src/Numeric/Vector.hs @@ -20,8 +20,10 @@ -- * Common operations , (.*.), dot, (·) , normL1, normL2, normLPInf, normLNInf, normLP + , normalized , vec2, vec3, vec4 , det2, cross, (×) + , unpackV2, unpackV3, unpackV4 ) where import Numeric.Array.ElementWise @@ -86,6 +88,16 @@ => Vector t n -> Scalar t normL2 = scalar . sqrt . ewfoldr (const (\a -> (a*a +))) 0 +-- | Normalize vector w.r.t. Euclidean metric (L2). +normalized :: ( Floating t + , Fractional (Vector t n) + , ElementWise (Idx '[n]) t (Vector t n) + ) + => Vector t n -> Vector t n +normalized v = v / n + where + n = broadcast . sqrt $ ewfoldr (const (\a -> (a*a +))) 0 v + -- | Maximum of absolute values normLPInf :: ( Ord t, Num t , ElementWise (Idx '[n]) t (Vector t n) @@ -169,3 +181,21 @@ f (2 :! Z) = b f (3 :! Z) = c f _ = d + + +unpackV2 :: ElementWise (Idx '[2]) t (Vector t 2) + => Vector t 2 -> (t, t) +unpackV2 v = (v ! 1, v ! 2) +{-# INLINE unpackV2 #-} + + +unpackV3 :: ElementWise (Idx '[3]) t (Vector t 3) + => Vector t 3 -> (t, t, t) +unpackV3 v = (v ! 1, v ! 2, v ! 3) +{-# INLINE unpackV3 #-} + + +unpackV4 :: ElementWise (Idx '[4]) t (Vector t 4) + => Vector t 4 -> (t, t, t, t) +unpackV4 v = (v ! 1, v ! 2, v ! 3, v ! 4) +{-# INLINE unpackV4 #-} diff --git a/test/Numeric/QuaternionTest.hs b/test/Numeric/QuaternionTest.hs new file mode 100644 --- /dev/null +++ b/test/Numeric/QuaternionTest.hs @@ -0,0 +1,124 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE FlexibleInstances #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +module Numeric.QuaternionTest (runTests) where + +import Test.QuickCheck +import Numeric.Quaternion +import Numeric.Vector + + +instance (Quaternion t, Arbitrary t) => Arbitrary (Quater t) where + arbitrary = packQ <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary + shrink q | (x,y,z,t) <- unpackQ q = packQ <$> shrink x <*> shrink y <*> shrink z <*> shrink t + +instance Arbitrary Vec3d where + arbitrary = vec3 <$> arbitrary <*> arbitrary <*> arbitrary + shrink v | (x,y,z) <- unpackV3 v = vec3 <$> shrink x <*> shrink y <*> shrink z + + +(=~=) :: (Quaternion a, Num a, Ord a, Num (Quater a), Fractional a) => Quater a -> Quater a -> Bool +(=~=) a b = taker (abs (a - b)) <= eps + where + s = max 1e-6 $ max (taker (abs a)) (taker (abs b)) + eps = s * 1e-6 +infix 4 =~= + + +prop_Eq :: QDouble -> Bool +prop_Eq q = and + [ (\(x,y,z,t) -> packQ x y z t) (unpackQ q) == q + , packQ (takei q) (takej q) (takek q) (taker q) == q + , fromVecNum (imVec q) (taker q) == q + , im q + re q == q + ] + + +prop_DoubleConjugate :: QDouble -> Bool +prop_DoubleConjugate q = conjugate (conjugate q) == q + +prop_Square :: QDouble -> Bool +prop_Square q = q * conjugate q =~= realToFrac (square q) + + +prop_RotScale :: QDouble -> Vec3d -> Bool +prop_RotScale q v = fromVecNum (rotScale q v) 0 =~= q * fromVecNum v 0 * conjugate q + +prop_GetRotScale :: Vec3d -> Vec3d -> Bool +prop_GetRotScale a b = fromVecNum (rotScale q a) 0 =~= fromVecNum b 0 + where + q = getRotScale a b + + +prop_InverseRotScale :: QDouble -> Vec3d -> Bool +prop_InverseRotScale q v | q /= 0 = fromVecNum (rotScale (1/q) (rotScale q v)) 0 =~= fromVecNum v 0 + | otherwise = True + +prop_FromToMatrix33 :: QDouble -> Bool +prop_FromToMatrix33 q | q == 0 = True + | otherwise = fromMatrix33 (toMatrix33 (signum q)) =~= signum q * signum (re q) + +prop_FromToMatrix44 :: QDouble -> Bool +prop_FromToMatrix44 q | q == 0 = True + | otherwise = fromMatrix44 (toMatrix44 (signum q)) =~= signum q * signum (re q) + +prop_RotationArg :: QDouble -> Bool +prop_RotationArg q | q == 0 = axisRotation (imVec q) (qArg q) =~= 1 + | otherwise = axisRotation (imVec q) (qArg q) =~= signum q + + +prop_UnitQ :: QDouble -> Bool +prop_UnitQ q | q == 0 = True + | otherwise = realToFrac (square (q / q)) =~= (1 :: QDouble) + + +prop_ExpLog :: QDouble -> Bool +prop_ExpLog q | q == 0 = log (exp q) == q + | otherwise = exp (log q) =~= q + +prop_SinAsin :: QDouble -> Bool +prop_SinAsin q = sin (asin q') =~= q' + where + q' = signum q + +prop_CosAcos :: QDouble -> Bool +prop_CosAcos q = cos (acos q') =~= q' + where + q' = signum q + +prop_TanAtan :: QDouble -> Bool +prop_TanAtan q = tan (atan q') =~= q' + where -- protect agains big numbers and rounding errors + q' = if square q >= 1e6 then signum q else q + + +prop_SinhAsinh :: QDouble -> Bool +prop_SinhAsinh q = sinh (asinh q') =~= q' + where -- protect agains big numbers and rounding errors + q' = if square q >= 1e6 then signum q else q + +prop_CoshAcosh :: QDouble -> Bool +prop_CoshAcosh q = cosh (acosh q') =~= q' + where -- protect agains big numbers and rounding errors + q' = if square q >= 1e6 then signum q else q + +prop_TanhAtanh :: QDouble -> Bool +prop_TanhAtanh q = tanh (atanh q') =~= q' + where -- protect agains big numbers and rounding errors + q' = if square q >= 1e6 then signum q else q + +prop_SinCos :: QDouble -> Bool +prop_SinCos q = sin q' * sin q' + cos q' * cos q' =~= 1 + where + q' = signum q + +prop_SinhCosh :: QDouble -> Bool +prop_SinhCosh q = cosh q' * cosh q' - sinh q' * sinh q' =~= 1 + where + q' = signum q + +return [] +runTests :: IO Bool +runTests = $quickCheckAll diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -1,17 +1,40 @@ -module Spec (tests) where +module Main (tests, main) where +import System.Exit import Distribution.TestSuite import qualified Numeric.DataFrame.BasicTest import qualified Numeric.DataFrame.SubSpaceTest +import qualified Numeric.QuaternionTest + +-- | Collection of tests in detailed-0.9 format tests :: IO [Test] tests = return [ test "DataFrame.Basic" Numeric.DataFrame.BasicTest.runTests , test "DataFrame.SubSpace" Numeric.DataFrame.SubSpaceTest.runTests + , test "Quaternion" Numeric.QuaternionTest.runTests ] + + +-- | Run tests as exitcode-stdio-1.0 +main :: IO () +main = do + ts <- tests + trs <- mapM (\(Test ti) ->(,) (name ti) <$> run ti) ts + case filter (not . isGood) trs of + [] -> exitSuccess + xs -> do + putStrLn $ "Failed tests: " ++ unwords (fmap fst xs) + exitFailure + where + isGood (_, Finished Pass) = True + isGood _ = False + + +-- | Convert QuickCheck props into Cabal tests test :: String -> IO Bool -> Test test tName propOp = Test testI where