diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,7 @@ +Copyright (c) 2016 Artem M. Chirkin + +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/app/Main.hs b/app/Main.hs new file mode 100644 --- /dev/null +++ b/app/Main.hs @@ -0,0 +1,112 @@ + +{-# LANGUAGE DataKinds #-} +module Main where + + +--import qualified Numeric.Tensor as T +import Numeric.Commons +-- import Numeric.Vector ((<:>)) +-- import qualified Numeric.Vector as V +-- import qualified Numeric.Matrix as M +import Numeric.EasyTensor + +import Foreign.Storable +import Foreign.Ptr +import Foreign.Marshal + + +main :: IO () +main = do + putStrLn "Hello world!" + print (two + vec2 3 4) + print (two + vec2 3 4 + 5) + print (two <:> two == two <:> x) + print (two <:> two == two <:> two) + print (two <:> x) + print x + print (x < two) + print (fromBytes (toBytes x) + two / 7 + 5) + print ((x <:> two * 3 <:> two) / 4.2 <:> x) + putStrLn "Done pure!" + ptr <- mallocArray 3 + poke ptr x + pokeElemOff ptr 1 5 + pokeByteOff ptr 16 two + peek ptr >>= print + peekByteOff ptr 8 >>= (print :: Vec2f -> IO ()) + peek (plusPtr ptr 8 :: Ptr Vec2f) >>= print + peekElemOff ptr 2 >>= print + peek (plusPtr ptr 16 :: Ptr Vec2f) >>= print + peekByteOff ptr 0 >>= (print :: Vec Float 6 -> IO ()) + peekByteOff ptr 4 >>= (print :: Vec Float 3 -> IO ()) + putStrLn "Done IO!" + + + putStrLn "Matrices..." + print $ index 1 2 m1 + print m1 + putStrLn (show (index 1 1 m1) ++ " " ++ show (index 1 2 m1) ) + putStrLn (show (index 2 1 m1) ++ " " ++ show (index 2 2 m1) ) + print (indexCol 1 m1) + print (indexCol 2 m1) + print (indexRow 1 m1) + print (indexRow 2 m1) + print (indexRow 3 m1) +-- print (M.indexMat 2 4 m1) +-- print (M.indexCol 4 m1 :: V.Vec2f) + + putStrLn "Matrix products" + print y2 + + putStrLn "EasyTensor" +-- print a +-- print b +-- print c +-- print d +-- print (a %* c) + print (two <:> x <:> 7 / x) + print m32 + print m33 +-- print x3 +-- print v3 +-- print (v3 // 4) +-- print (10 \\ 4 :: Tensor Double 1 1) +-- print $ a %* c +-- print $ transpose x3 %* m32 + putStrLn "m33" + print m33 + putStrLn "determinants" + print (2 * eye :: Tensor Float 5 5) + print $ det (2 * eye :: Tensor Float 5 5) + print m33 + print $ det m33 + print $ two <:> x + print $ det (two <:> x) + print (x3 <:> (m33 %* x3) <:> (m33 %* x3)) + print $ det (x3 <:> (m33 %* x3) <:> (m33 %* x3)) + putStrLn "Inverse" + print m33 + print (inverse m33) + print $ inverse m33 %* m33 + print $ m33 %* inverse m33 + print (2 * eye :: Tensor Float 5 5) + print $ inverse (2 * eye :: Tensor Float 5 5) + print $ (2 * eye :: Tensor Float 5 5) * inverse (2 * eye :: Tensor Float 5 5) + print (x3 <:> v3 <:> v3) + print $ inverse (x3 <:> v3 <:> v3) + print $ (x3 <:> v3 <:> v3) %* inverse (x3 <:> v3 <:> v3) + where + two = vec2 2 2.001 :: Vec2f + x = two / vec2 3.2 (-2) + m1 = fromBytes (toBytes (two <:> x <:> 7 / x)) :: Mat Float 3 2 + m32 = m1 + x2 = 7 :: Vec2f + x3 = 9 :: Vec3f + y2 = m32 %* x2 + m33 = m32 <:> 17 + v3 = m33 %* x3 + +-- a = 1 :: Tensor Float 2 2 +-- b = 3 :: Tensor Float 1 1 +-- c = 4 :: Tensor Float 2 1 +-- d = 5 :: Tensor Float 1 2 diff --git a/easytensor.cabal b/easytensor.cabal new file mode 100644 --- /dev/null +++ b/easytensor.cabal @@ -0,0 +1,72 @@ +name: easytensor +version: 0.1.0.0 +cabal-version: >=1.20 +build-type: Simple +license: MIT +license-file: LICENSE +copyright: (c) Artem Chirkin +maintainer: chirkin@arch.ethz.ch +homepage: https://github.com/achirkin/easytensor#readme +synopsis: Initial project template from stack +description: + Pure haskell vector/matrix/tensor library. + Features dimensionality type-checking for all operations. + Generic n-dimensional versions are implemented using low-level prim ops. + Allows ad-hoc replacement with fixed low-dimensionality vectors and matrices without changing user interface. + Please see README.md +category: Math, Geometry +author: Artem Chirkin + +source-repository head + type: git + location: https://github.com/achirkin/easytensor.git + +flag unsafeindices + description: + Disable bound checks when accessing elements of a tensor/matrix/vector. + default: False + +library + + if flag(unsafeindices) + cpp-options: -DUNSAFE_INDICES + exposed-modules: + Numeric.EasyTensor + Numeric.Tensor + Numeric.Commons + build-depends: + base >=4.9 && <5, + ghc-prim >=0.5 + default-language: Haskell2010 + hs-source-dirs: src + other-modules: + Numeric.Matrix + Numeric.Matrix.Base.FloatXNM + Numeric.Matrix.Class + Numeric.Matrix.Family + Numeric.Vector + Numeric.Vector.Base.FloatX2 + Numeric.Vector.Base.FloatXN + Numeric.Vector.Class + Numeric.Vector.Family + ghc-options: -Wall -fwarn-tabs + +executable easytensor-exe + main-is: Main.hs + build-depends: + base -any, + easytensor -any, + ghc-prim >=0.5 + default-language: Haskell2010 + hs-source-dirs: app + ghc-options: -Wall -fwarn-tabs -threaded -rtsopts -with-rtsopts=-N + +test-suite easytensor-test + type: exitcode-stdio-1.0 + main-is: Spec.hs + build-depends: + base -any, + easytensor -any + default-language: Haskell2010 + hs-source-dirs: test + ghc-options: -Wall -fwarn-tabs -threaded -rtsopts -with-rtsopts=-N diff --git a/src/Numeric/Commons.hs b/src/Numeric/Commons.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Commons.hs @@ -0,0 +1,307 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MagicHash, UnboxedTuples #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Commons +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Commons + ( PrimBytes (..) + , FloatBytes (..) + , DoubleBytes (..) + , IntBytes (..) + , WordBytes (..) + ) where + +#include "MachDeps.h" +#include "HsBaseConfig.h" + +import GHC.Base (runRW#) +import GHC.Ptr +import GHC.Prim +import GHC.Types +import GHC.Int +import GHC.Word +import Foreign.Storable + + +class PrimBytes a where + -- | Store content of a data type in a primitive byte array + toBytes :: a -> ByteArray# + -- | Load content of a data type from a primitive byte array + fromBytes :: ByteArray# -> a + -- | Size of a data type in bytes + byteSize :: a -> Int# + -- | Alignment of a data type in bytes + byteAlign :: a -> Int# + +-- | Primitive indexing. No checks, no safety. +class FloatBytes a where + -- | Primitive get Float# (element offset) + ixF :: Int# -> a -> Float# + +-- | Primitive indexing. No checks, no safety. +class DoubleBytes a where + -- | Primitive get Double# (element offset) + ixD :: Int# -> a -> Double# + +-- | Primitive indexing. No checks, no safety. +class IntBytes a where + -- | Primitive get Int# (element offset) + ixI :: Int# -> a -> Int# + +-- | Primitive indexing. No checks, no safety. +class WordBytes a where + -- | Primitive get Word# (element offset) + ixW :: Int# -> a -> Word# + +instance PrimBytes a => Storable a where + sizeOf x = I# (byteSize x) + alignment x = I# (byteAlign x) + peekElemOff ptr (I# offset) = peekByteOff ptr (I# (offset *# byteSize (undefined :: a))) + pokeElemOff ptr (I# offset) = pokeByteOff ptr (I# (offset *# byteSize (undefined :: a))) + peekByteOff (Ptr addr) (I# offset) = IO $ \s0 -> case newByteArray# bsize s0 of + (# s1, marr #) -> case copyAddrToByteArray# (addr `plusAddr#` offset) marr 0# bsize s1 of + s2 -> case unsafeFreezeByteArray# marr s2 of + (# s3, arr #) -> (# s3, fromBytes arr #) + where + bsize = byteSize (undefined :: a) + pokeByteOff (Ptr addr) (I# offset) x = IO $ \s0 -> case copyByteArrayToAddr# (toBytes x) + 0# + (addr `plusAddr#` offset) + bsize s0 of + s2 -> (# s2, () #) + where + bsize = byteSize (undefined :: a) + peek ptr = peekByteOff ptr 0 + poke ptr = pokeByteOff ptr 0 + + +instance PrimBytes Float where + toBytes v@(F# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeFloatArray# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = F# (indexFloatArray# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_HSFLOAT# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_HSFLOAT# + {-# INLINE byteAlign #-} + +instance FloatBytes Float where + ixF _ (F# x) = x + {-# INLINE ixF #-} + +instance PrimBytes Double where + toBytes v@(D# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeDoubleArray# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = D# (indexDoubleArray# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_HSDOUBLE# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_HSDOUBLE# + {-# INLINE byteAlign #-} + +instance DoubleBytes Double where + ixD _ (D# x) = x + {-# INLINE ixD #-} + +instance PrimBytes Int where + toBytes v@(I# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeIntArray# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = I# (indexIntArray# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_HSINT# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_HSINT# + {-# INLINE byteAlign #-} + +instance IntBytes Int where + ixI _ (I# x) = x + {-# INLINE ixI #-} + +instance PrimBytes Int8 where + toBytes v@(I8# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeInt8Array# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = I8# (indexInt8Array# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_INT8# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_INT8# + {-# INLINE byteAlign #-} + +instance IntBytes Int8 where + ixI _ (I8# x) = x + {-# INLINE ixI #-} + +instance PrimBytes Int16 where + toBytes v@(I16# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeInt16Array# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = I16# (indexInt16Array# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_INT16# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_INT16# + {-# INLINE byteAlign #-} + +instance IntBytes Int16 where + ixI _ (I16# x) = x + {-# INLINE ixI #-} + +instance PrimBytes Int32 where + toBytes v@(I32# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeInt32Array# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = I32# (indexInt32Array# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_INT32# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_INT32# + {-# INLINE byteAlign #-} + +instance IntBytes Int32 where + ixI _ (I32# x) = x + {-# INLINE ixI #-} + +instance PrimBytes Int64 where + toBytes v@(I64# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeInt64Array# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = I64# (indexInt64Array# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_INT64# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_INT64# + {-# INLINE byteAlign #-} + +instance IntBytes Int64 where + ixI _ (I64# x) = x + {-# INLINE ixI #-} + +instance PrimBytes Word where + toBytes v@(W# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeWordArray# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = W# (indexWordArray# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_HSWORD# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_HSWORD# + {-# INLINE byteAlign #-} + +instance WordBytes Word where + ixW _ (W# x) = x + {-# INLINE ixW #-} + +instance PrimBytes Word8 where + toBytes v@(W8# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeWord8Array# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = W8# (indexWord8Array# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_WORD8# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_WORD8# + {-# INLINE byteAlign #-} + +instance WordBytes Word8 where + ixW _ (W8# x) = x + {-# INLINE ixW #-} + +instance PrimBytes Word16 where + toBytes v@(W16# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeWord16Array# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = W16# (indexWord16Array# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_WORD16# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_WORD16# + {-# INLINE byteAlign #-} + +instance WordBytes Word16 where + ixW _ (W16# x) = x + {-# INLINE ixW #-} + +instance PrimBytes Word32 where + toBytes v@(W32# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeWord32Array# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = W32# (indexWord32Array# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_WORD32# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_WORD32# + {-# INLINE byteAlign #-} + +instance WordBytes Word32 where + ixW _ (W32# x) = x + {-# INLINE ixW #-} + +instance PrimBytes Word64 where + toBytes v@(W64# x) = case runRW# + ( \s0 -> case newByteArray# (byteSize v) s0 of + (# s1, marr #) -> case writeWord64Array# marr 0# x s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = W64# (indexWord64Array# arr 0#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_WORD64# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_WORD64# + {-# INLINE byteAlign #-} + +instance WordBytes Word64 where + ixW _ (W64# x) = x + {-# INLINE ixW #-} + + + diff --git a/src/Numeric/EasyTensor.hs b/src/Numeric/EasyTensor.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/EasyTensor.hs @@ -0,0 +1,557 @@ +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeFamilyDependencies #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE MagicHash, UnboxedTuples #-} +{-# LANGUAGE KindSignatures, DataKinds #-} +{-# LANGUAGE TypeOperators, FlexibleInstances, ScopedTypeVariables #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.EasyTensor +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- This module generalizes matrices and vectors. +-- Yet it is limited to rank 2, allowing for a simple and nicely type-checked interface. +-- +-- +-- +----------------------------------------------------------------------------- + +module Numeric.EasyTensor + ( Tensor () + -- * Common operations + , fill + , prod, (%*) + , inverse, transpose + , (<:>), (//), (\\) + , index, indexCol, indexRow, dimN, dimM + , (V..*.), dot, (·), normL1, normL2, normLPInf, normLNInf, normLP + , eye, diag, det, trace + , toDiag, toDiag', fromDiag, fromDiag' + -- * Type abbreviations + , Mat, Vec, Vec' + , Vec2f, Vec3f, Vec4f + , Vec2f', Vec3f', Vec4f' + , Mat22f, Mat23f, Mat24f + , Mat32f, Mat33f, Mat34f + , Mat42f, Mat43f, Mat44f + -- * Simplified type constructors + , scalar, vec2, vec3, vec4, vec2', vec3', vec4' + , mat22, mat33, mat44 + -- * Some low-dimensional operations + , det2, det2', cross, (×) + ) where + +import GHC.Base (runRW#) +import GHC.Prim +import GHC.Types +import GHC.TypeLits +import Unsafe.Coerce (unsafeCoerce) + +import qualified Numeric.Vector as V +import qualified Numeric.Matrix as M +import qualified Numeric.Matrix.Class as M +import Numeric.Commons + + + +newtype Tensor t n m = Tensor { _unT :: TT t n m } +instance Show (TT t n m) => Show (Tensor t n m) where + show (Tensor t) = show t +deriving instance Eq (TT t n m) => Eq (Tensor t n m) +deriving instance Ord (TT t n m) => Ord (Tensor t n m) +deriving instance Num (TT t n m) => Num (Tensor t n m) +deriving instance Fractional (TT t n m) => Fractional (Tensor t n m) +deriving instance Floating (TT t n m) => Floating (Tensor t n m) +deriving instance V.VectorCalculus t n (TT t n 1) => V.VectorCalculus t n (Tensor t n 1) +deriving instance V.VectorCalculus t m (TT t 1 m) => V.VectorCalculus t m (Tensor t 1 m) +deriving instance M.MatrixCalculus t n m (TT t n m) => M.MatrixCalculus t n m (Tensor t n m) +deriving instance M.SquareMatrixCalculus t n (TT t n n) => M.SquareMatrixCalculus t n (Tensor t n n) +deriving instance M.MatrixInverse (TT t n n) => M.MatrixInverse (Tensor t n n) +deriving instance PrimBytes (TT t n m) => PrimBytes (Tensor t n m) +deriving instance FloatBytes (TT t n m) => FloatBytes (Tensor t n m) +deriving instance DoubleBytes (TT t n m) => DoubleBytes (Tensor t n m) +deriving instance IntBytes (TT t n m) => IntBytes (Tensor t n m) +deriving instance WordBytes (TT t n m) => WordBytes (Tensor t n m) + + + + +newtype Scalar t = Scalar { _unScalar :: t } + deriving ( Bounded, Enum, Eq, Integral, Num, Fractional, Floating, Ord, Read, Real, RealFrac, RealFloat + , PrimBytes, FloatBytes, DoubleBytes, IntBytes, WordBytes) +instance Show t => Show (Scalar t) where + show (Scalar t) = "{ " ++ show t ++ " }" + +newtype CoVector t n = CoVector {_unCoVec :: V.Vector t n} +instance Show (V.Vector t n) => Show (CoVector t n) where + show (CoVector t) = show t +deriving instance Eq (V.Vector t n) => Eq (CoVector t n) +deriving instance Ord (V.Vector t n) => Ord (CoVector t n) +deriving instance Num (V.Vector t n) => Num (CoVector t n) +deriving instance Fractional (V.Vector t n) => Fractional (CoVector t n) +deriving instance Floating (V.Vector t n) => Floating (CoVector t n) +deriving instance PrimBytes (V.Vector t n) => PrimBytes (CoVector t n) +deriving instance V.VectorCalculus t n (V.Vector t n) => V.VectorCalculus t n (CoVector t n) +deriving instance FloatBytes (V.Vector t n) => FloatBytes (CoVector t n) +deriving instance DoubleBytes (V.Vector t n) => DoubleBytes (CoVector t n) +deriving instance IntBytes (V.Vector t n) => IntBytes (CoVector t n) +deriving instance WordBytes (V.Vector t n) => WordBytes (CoVector t n) + +newtype ContraVector t n = ContraVector {_unContraVec :: V.Vector t n} +instance Show (V.Vector t n) => Show (ContraVector t n) where + show (ContraVector t) = show t ++ "'" +deriving instance Eq (V.Vector t n) => Eq (ContraVector t n) +deriving instance Ord (V.Vector t n) => Ord (ContraVector t n) +deriving instance Num (V.Vector t n) => Num (ContraVector t n) +deriving instance Fractional (V.Vector t n) => Fractional (ContraVector t n) +deriving instance Floating (V.Vector t n) => Floating (ContraVector t n) +deriving instance PrimBytes (V.Vector t n) => PrimBytes (ContraVector t n) +deriving instance V.VectorCalculus t n (V.Vector t n) => V.VectorCalculus t n (ContraVector t n) +deriving instance FloatBytes (V.Vector t n) => FloatBytes (ContraVector t n) +deriving instance DoubleBytes (V.Vector t n) => DoubleBytes (ContraVector t n) +deriving instance IntBytes (V.Vector t n) => IntBytes (ContraVector t n) +deriving instance WordBytes (V.Vector t n) => WordBytes (ContraVector t n) + +newtype Matrix t n m = Matrix (M.Matrix t n m) +instance Show (M.Matrix t n m) => Show (Matrix t n m) where + show (Matrix t) = show t +deriving instance Eq (M.Matrix t n m) => Eq (Matrix t n m) +deriving instance Ord (M.Matrix t n m) => Ord (Matrix t n m) +deriving instance Num (M.Matrix t n m) => Num (Matrix t n m) +deriving instance Fractional (M.Matrix t n m) => Fractional (Matrix t n m) +deriving instance Floating (M.Matrix t n m) => Floating (Matrix t n m) +deriving instance PrimBytes (M.Matrix t n m) => PrimBytes (Matrix t n m) +deriving instance M.MatrixCalculus t n m (M.Matrix t n m) => M.MatrixCalculus t n m (Matrix t n m) +deriving instance M.SquareMatrixCalculus t n (M.Matrix t n n) => M.SquareMatrixCalculus t n (Matrix t n n) +deriving instance M.MatrixInverse (M.Matrix t n n) => M.MatrixInverse (Matrix t n n) +deriving instance FloatBytes (M.Matrix t n m) => FloatBytes (Matrix t n m) +deriving instance DoubleBytes (M.Matrix t n m) => DoubleBytes (Matrix t n m) +deriving instance IntBytes (M.Matrix t n m) => IntBytes (Matrix t n m) +deriving instance WordBytes (M.Matrix t n m) => WordBytes (Matrix t n m) + + + + +type family TT t (n :: Nat) (m :: Nat) = v | v -> t n m where + TT t 1 1 = Scalar t + TT t n 1 = ContraVector t n + TT t 1 m = CoVector t m + TT t n m = Matrix t n m + +-- | Fill whole tensor with a single value +fill :: M.MatrixCalculus t n m (Tensor t n m) => Tensor t 1 1 -> Tensor t n m +fill = M.broadcastMat . _unScalar . _unT +{-# INLINE fill #-} + +-- | Get an element of a tensor +index :: M.MatrixCalculus t n m (Tensor t n m) => Int -> Int -> Tensor t n m -> Tensor t 1 1 +index i j = Tensor . Scalar . M.indexMat i j +{-# INLINE index #-} + +-- | Get a column vector of a matrix +indexCol :: ( M.MatrixCalculus t n m (Tensor t n m) + , V.VectorCalculus t n (Tensor t n 1) + , PrimBytes (Tensor t n 1) + ) + => Int -> Tensor t n m -> Tensor t n 1 +indexCol = M.indexCol +{-# INLINE indexCol #-} + +-- | Get a row vector of a matrix +indexRow :: ( M.MatrixCalculus t n m (Tensor t n m) + , V.VectorCalculus t m (Tensor t 1 m) + , PrimBytes (Tensor t 1 m) + ) + => Int -> Tensor t n m -> Tensor t 1 m +indexRow = M.indexRow +{-# INLINE indexRow #-} + +transpose :: ( M.MatrixCalculus t n m (Tensor t n m) + , M.MatrixCalculus t m n (Tensor t m n) + , PrimBytes (Tensor t m n) + ) => Tensor t n m -> Tensor t m n +transpose = M.transpose +{-# INLINE transpose #-} + + +dimN :: ( M.MatrixCalculus t n m (Tensor t n m)) => Tensor t n m -> Int +dimN = M.dimN +{-# INLINE dimN #-} + +dimM :: ( M.MatrixCalculus t n m (Tensor t n m)) => Tensor t n m -> Int +dimM = M.dimM +{-# INLINE dimM #-} + +-- | Matrix product for tensors rank 2, as well matrix-vector or vector-matrix products +infixl 7 %* +(%*) :: (M.MatrixProduct (Tensor t n m) (Tensor t m k) (Tensor t n k)) + => Tensor t n m -> Tensor t m k -> Tensor t n k +(%*) = prod +{-# INLINE (%*) #-} + +-- | Matrix product for tensors rank 2, as well matrix-vector or vector-matrix products +prod :: (M.MatrixProduct (Tensor t n m) (Tensor t m k) (Tensor t n k)) + => Tensor t n m -> Tensor t m k -> Tensor t n k +prod = M.prod +{-# INLINE prod #-} + +-- | Divide on the right: R = A * B^(-1) +(//) :: ( M.MatrixProduct (Tensor t n m) (Tensor t m m) (Tensor t n m) + , M.MatrixInverse (TT t m m)) + => Tensor t n m -> Tensor t m m -> Tensor t n m +(//) a b = prod a (inverse b) +{-# INLINE (//) #-} + +-- | Divide on the left: R = A^(-1) * b +(\\) :: ( M.MatrixProduct (Tensor t n n) (Tensor t n m) (Tensor t n m) + , M.MatrixInverse (TT t n n)) + => Tensor t n n -> Tensor t n m -> Tensor t n m +(\\) a b = prod (inverse a) b +{-# INLINE (\\) #-} + + +-- | Matrix inverse +inverse :: M.MatrixInverse (Tensor t n n) => Tensor t n n -> Tensor t n n +inverse = M.inverse +{-# INLINE inverse #-} + + +instance ( FloatBytes (Tensor Float n m) + , FloatBytes (Tensor Float m k) + , PrimBytes (Tensor Float n k) + , M.MatrixCalculus Float n m (Tensor Float n m) + , M.MatrixCalculus Float m k (Tensor Float m k) + ) + => M.MatrixProduct (Tensor Float n m) (Tensor Float m k) (Tensor Float n k) where + prod x y = case (dimN x, dimM x, dimM y) of + ( I# n, I# m, I# k) -> M.prodF n m k x y + {-# INLINE prod #-} + +instance ( DoubleBytes (Tensor Double n m) + , DoubleBytes (Tensor Double m k) + , PrimBytes (Tensor Double n k) + , M.MatrixCalculus Double n m (Tensor Double n m) + , M.MatrixCalculus Double m k (Tensor Double m k) + ) + => M.MatrixProduct (Tensor Double n m) (Tensor Double m k) (Tensor Double n k) where + prod x y = case (dimN x, dimM x, dimM y) of + ( I# n, I# m, I# k) -> M.prodD n m k x y + {-# INLINE prod #-} + + + +-- | Append one vector to another, adding up their dimensionality +(<:>) :: ( PrimBytes (Tensor t k n) + , PrimBytes (Tensor t k m) + , PrimBytes (Tensor t k (n + m)) + ) + => Tensor t k n -> Tensor t k m -> Tensor t k (n + m) +a <:> b = case (# toBytes a, toBytes b, byteSize a, byteSize b #) of + (# arr1, arr2, n, m #) -> case runRW# + ( \s0 -> case newByteArray# (n +# m) s0 of + (# s1, marr #) -> case copyByteArray# arr1 0# marr 0# n s1 of + s2 -> case copyByteArray# arr2 0# marr n m s2 of + s3 -> unsafeFreezeByteArray# marr s3 + ) of (# _, r #) -> fromBytes r +infixl 5 <:> + +-- Simple types + +type Mat t n m = Tensor t n m +type Vec t n = Tensor t n 1 +type Vec' t m = Tensor t 1 m + + +-- Even Sympler types + +type Vec2f = Tensor Float 2 1 +type Vec3f = Tensor Float 3 1 +type Vec4f = Tensor Float 4 1 +type Vec2f' = Tensor Float 1 2 +type Vec3f' = Tensor Float 1 3 +type Vec4f' = Tensor Float 1 4 + + +type Mat22f = Tensor Float 2 2 +type Mat32f = Tensor Float 3 2 +type Mat42f = Tensor Float 4 2 +type Mat23f = Tensor Float 2 3 +type Mat33f = Tensor Float 3 3 +type Mat43f = Tensor Float 4 3 +type Mat24f = Tensor Float 2 4 +type Mat34f = Tensor Float 3 4 +type Mat44f = Tensor Float 4 4 + + +-- construct tensors + +scalar :: t -> Tensor t 1 1 +scalar = Tensor . Scalar + +vec2 :: V.Vector2D t => t -> t -> Tensor t 2 1 +vec2 a b = Tensor . ContraVector $ V.vec2 a b + +vec2' :: V.Vector2D t => t -> t -> Tensor t 1 2 +vec2' a b = Tensor . CoVector $ V.vec2 a b + +vec3 :: V.Vector3D t => t -> t -> t -> Tensor t 3 1 +vec3 a b c = Tensor . ContraVector $ V.vec3 a b c + +vec3' :: V.Vector3D t => t -> t -> t -> Tensor t 1 3 +vec3' a b c = Tensor . CoVector $ V.vec3 a b c + +vec4 :: V.Vector4D t => t -> t -> t -> t -> Tensor t 4 1 +vec4 a b c d = Tensor . ContraVector $ V.vec4 a b c d + +vec4' :: V.Vector4D t => t -> t -> t -> t -> Tensor t 1 4 +vec4' a b c d = Tensor . CoVector $ V.vec4 a b c d + + + +-- | Compose a 2x2D matrix +mat22 :: M.Matrix2x2 t => Tensor t 2 1 -> Tensor t 2 1 -> Tensor t 2 2 +mat22 (Tensor (ContraVector a)) (Tensor (ContraVector b)) = Tensor . Matrix $ M.mat22 a b + +-- | Compose a 3x3D matrix +mat33 :: ( PrimBytes (Tensor t 3 3) + , PrimBytes (Tensor t 3 2) + , PrimBytes (Tensor t 3 1) + ) + => Tensor t 3 1 -> Tensor t 3 1 -> Tensor t 3 1 -> Tensor t 3 3 +mat33 a b c = a <:> b <:> c + +-- | Compose a 4x4D matrix +mat44 :: ( PrimBytes (Tensor t 4 4) + , PrimBytes (Tensor t 4 3) + , PrimBytes (Tensor t 4 2) + , PrimBytes (Tensor t 4 1) + ) + => Tensor t 4 1 -> Tensor t 4 1 -> Tensor t 4 1 -> Tensor t 4 1 -> Tensor t 4 4 +mat44 a b c d = a <:> b <:> c <:> d + +-- useful low-dimensional functions + +det2 :: V.Vector2D t => Tensor t 2 1 -> Tensor t 2 1 -> Tensor t 1 1 +det2 (Tensor (ContraVector a)) (Tensor (ContraVector b)) = Tensor . Scalar $ V.det2 a b + +det2' :: V.Vector2D t => Tensor t 1 2 -> Tensor t 1 2 -> Tensor t 1 1 +det2' (Tensor (CoVector a)) (Tensor (CoVector b)) = Tensor . Scalar $ V.det2 a b + +-- | Cross product for two vectors in 3D +infixl 7 × +(×) :: V.Vector3D t => Tensor t 3 1 -> Tensor t 3 1 -> Tensor t 3 1 +(×) = cross +{-# INLINE (×) #-} + +-- | Cross product for two vectors in 3D +cross :: V.Vector3D t => Tensor t 3 1 -> Tensor t 3 1 -> Tensor t 3 1 +cross (Tensor (ContraVector a)) (Tensor (ContraVector b)) = Tensor . ContraVector $ V.cross a b +{-# INLINE cross #-} + + +-- re-use functions provided by Vector and Matrix Calculus + +-- | Dot product of two vectors +infixl 7 · +(·) :: V.VectorCalculus t n v => v -> v -> Tensor t 1 1 +(·) = dot +{-# INLINE (·) #-} + + +-- | Dot product of two vectors +dot :: V.VectorCalculus t n v => v -> v -> Tensor t 1 1 +dot a b = Tensor . Scalar $ V.dot a b +{-# INLINE dot #-} + +-- | Sum of absolute values +normL1 :: V.VectorCalculus t n v => v -> Tensor t 1 1 +normL1 = Tensor . Scalar . V.normL1 +{-# INLINE normL1 #-} + +-- | hypot function (square root of squares) +normL2 :: V.VectorCalculus t n v => v -> Tensor t 1 1 +normL2 = Tensor . Scalar . V.normL2 +{-# INLINE normL2 #-} + +-- | Maximum of absolute values +normLPInf :: V.VectorCalculus t n v => v -> Tensor t 1 1 +normLPInf = Tensor . Scalar . V.normLPInf +{-# INLINE normLPInf #-} + +-- | Minimum of absolute values +normLNInf :: V.VectorCalculus t n v => v -> Tensor t 1 1 +normLNInf = Tensor . Scalar . V.normLNInf +{-# INLINE normLNInf #-} + +-- | Norm in Lp space +normLP :: V.VectorCalculus t n v => Int -> v -> Tensor t 1 1 +normLP p = Tensor . Scalar . V.normLP p +{-# INLINE normLP #-} + +-- | Identity matrix. Mat with 1 on diagonal and 0 elsewhere +eye :: M.SquareMatrixCalculus t n (Tensor t n n) => Tensor t n n +eye = M.eye +{-# INLINE eye #-} + +-- | Put the same value on the Mat diagonal, 0 otherwise +diag :: M.SquareMatrixCalculus t n (Tensor t n n) => Tensor t 1 1 -> Tensor t n n +diag = M.diag . _unScalar . _unT +{-# INLINE diag #-} + +-- | Determinant of Mat +det :: M.SquareMatrixCalculus t n (Tensor t n n) => Tensor t n n -> Tensor t 1 1 +det = Tensor . Scalar . M.det +{-# INLINE det #-} + +-- | Sum of diagonal elements +trace :: M.SquareMatrixCalculus t n (Tensor t n n) => Tensor t n n -> Tensor t 1 1 +trace = Tensor . Scalar . M.trace +{-# INLINE trace #-} + +-- | Get the diagonal elements from Mat into Vec +fromDiag :: ( M.SquareMatrixCalculus t n (Tensor t n n) + , V.VectorCalculus t n (Tensor t n 1) + , PrimBytes (Tensor t n 1) + ) + => Tensor t n n -> Tensor t n 1 +fromDiag = M.fromDiag +{-# INLINE fromDiag #-} + +-- | Get the diagonal elements from Mat into Vec +fromDiag' :: ( M.SquareMatrixCalculus t n (Tensor t n n) + , V.VectorCalculus t n (Tensor t 1 n) + , PrimBytes (Tensor t 1 n) + ) + => Tensor t n n -> Tensor t 1 n +fromDiag' = M.fromDiag +{-# INLINE fromDiag' #-} + +-- | Set Vec values into the diagonal elements of Mat +toDiag :: ( M.SquareMatrixCalculus t n (Tensor t n n) + , V.VectorCalculus t n (Tensor t n 1) + , PrimBytes (Tensor t n 1) + ) + => Tensor t n 1 -> Tensor t n n +toDiag = M.toDiag +{-# INLINE toDiag #-} + +-- | Set Vec values into the diagonal elements of Mat +toDiag' :: ( M.SquareMatrixCalculus t n (Tensor t n n) + , V.VectorCalculus t n (Tensor t 1 n) + , PrimBytes (Tensor t 1 n) + ) + => Tensor t 1 n -> Tensor t n n +toDiag' = M.toDiag +{-# INLINE toDiag' #-} + + +-- missing instances + + +instance Num t => V.VectorCalculus t 1 (Scalar t) where + broadcastVec = Scalar + {-# INLINE broadcastVec #-} + (.*.) = (*) + {-# INLINE (.*.) #-} + dot a = _unScalar . (a *) + {-# INLINE dot #-} + indexVec 1 = _unScalar + indexVec i = const . error $ "Bad index " ++ show i ++ " for a scalar" + {-# INLINE indexVec #-} + normL1 = _unScalar . abs + {-# INLINE normL1 #-} + normL2 = _unScalar . abs + {-# INLINE normL2 #-} + normLPInf = _unScalar + {-# INLINE normLPInf #-} + normLNInf = _unScalar + {-# INLINE normLNInf #-} + normLP _ = _unScalar . abs + {-# INLINE normLP #-} + dim _ = 1 + {-# INLINE dim #-} + +instance Num t => M.MatrixCalculus t 1 1 (Scalar t) where + broadcastMat = Scalar + {-# INLINE broadcastMat #-} + indexMat 1 1 = _unScalar + indexMat i j = const . error $ "Bad index (" ++ show i ++ ", " ++ show j ++ ") for a scalar" + {-# INLINE indexMat #-} + transpose = unsafeCoerce + {-# INLINE transpose #-} + dimN _ = 1 + {-# INLINE dimN #-} + dimM _ = 1 + {-# INLINE dimM #-} + indexCol 1 = unsafeCoerce + indexCol j = const . error $ "Bad column index " ++ show j ++ " for a scalar" + {-# INLINE indexCol #-} + indexRow 1 = unsafeCoerce + indexRow i = const . error $ "Bad row index " ++ show i ++ " for a scalar" + {-# INLINE indexRow #-} + + + +instance Num t => M.SquareMatrixCalculus t 1 (Scalar t) where + eye = Scalar 1 + {-# INLINE eye #-} + diag = Scalar + {-# INLINE diag #-} + det = _unScalar + {-# INLINE det #-} + trace = _unScalar + {-# INLINE trace #-} + fromDiag = unsafeCoerce + {-# INLINE fromDiag #-} + toDiag = unsafeCoerce + {-# INLINE toDiag #-} + +instance Fractional t => M.MatrixInverse (Scalar t) where + inverse = Scalar . recip . _unScalar + + +instance (KnownNat n, V.VectorCalculus t n (V.Vector t n)) => M.MatrixCalculus t n 1 (ContraVector t n) where + broadcastMat = ContraVector . V.broadcastVec + {-# INLINE broadcastMat #-} + indexMat i 1 (ContraVector v) = V.indexVec i v + indexMat i j (ContraVector v) = error $ "Bad index (" ++ show i ++ ", " ++ show j ++ ") for a " ++ show (V.dim v) ++ "x1D matrix" + {-# INLINE indexMat #-} + transpose (ContraVector v) = unsafeCoerce $ CoVector v + {-# INLINE transpose #-} + dimN = V.dim . _unContraVec + {-# INLINE dimN #-} + dimM _ = 1 + {-# INLINE dimM #-} + indexCol 1 (ContraVector v) = unsafeCoerce $ ContraVector v + indexCol j (ContraVector v) = error $ "Bad column index " ++ show j ++ " for a " ++ show (V.dim v) ++ "x1D matrix" + {-# INLINE indexCol #-} + indexRow i = unsafeCoerce . Scalar . V.indexVec i . _unContraVec + {-# INLINE indexRow #-} + + +instance (KnownNat m, V.VectorCalculus t m (V.Vector t m)) => M.MatrixCalculus t 1 m (CoVector t m) where + broadcastMat = CoVector . V.broadcastVec + {-# INLINE broadcastMat #-} + indexMat 1 i (CoVector v) = V.indexVec i v + indexMat j i (CoVector v) = error $ "Bad index (" ++ show j ++ ", " ++ show i ++ ") for a 1x" ++ show (V.dim v) ++ "D matrix" + {-# INLINE indexMat #-} + transpose (CoVector v) = unsafeCoerce $ ContraVector v + {-# INLINE transpose #-} + dimN _ = 1 + {-# INLINE dimN #-} + dimM = V.dim . _unCoVec + {-# INLINE dimM #-} + indexCol i = unsafeCoerce . Scalar . V.indexVec i . _unCoVec + {-# INLINE indexCol #-} + indexRow 1 x = unsafeCoerce x + indexRow j (CoVector v) = error $ "Bad column index " ++ show j ++ " for a 1x" ++ show (V.dim v) ++ "D matrix" + {-# INLINE indexRow #-} + + diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Matrix.hs @@ -0,0 +1,75 @@ +{-# LANGUAGE MagicHash, DataKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeOperators #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Matrix +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Matrix + ( Matrix + -- * Type classes + , MatrixCalculus (..), SquareMatrixCalculus (..) + , Matrix2x2 (..) + -- * Type abbreviations + , Mat22f, Mat23f, Mat24f + , Mat32f, Mat33f, Mat34f + , Mat42f, Mat43f, Mat44f + -- * Common functions on vectors +-- , (<:>) + , MatrixProduct (..) + , MatrixInverse (..) + ) where + +--import GHC.Base (runRW#) +--import GHC.Prim +--import GHC.TypeLits +-- +--import Numeric.Commons +import Numeric.Matrix.Class ( MatrixCalculus(..) + , SquareMatrixCalculus(..) + , Matrix2x2(..) + , MatrixProduct (..) + , MatrixInverse (..) + ) +import Numeric.Matrix.Family (Matrix) + + +-- Import instances +import Numeric.Matrix.Base.FloatXNM () + + +-- 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 + + + +---- | Append one vector to another, adding up their dimensionality +--(<:>) :: ( PrimBytes (Matrix t n1 m1) +-- , PrimBytes (Matrix t m) +-- , PrimBytes (Matrix t (n+m)) +-- ) +-- => Matrix t n -> Matrix t m -> Matrix t (n + m) +--a <:> b = case (# toBytes a, toBytes b, byteSize a, byteSize b #) of +-- (# arr1, arr2, n, m #) -> case runRW# +-- ( \s0 -> case newByteArray# (n +# m) s0 of +-- (# s1, marr #) -> case copyByteArray# arr1 0# marr 0# n s1 of +-- s2 -> case copyByteArray# arr2 0# marr n m s2 of +-- s3 -> unsafeFreezeByteArray# marr s3 +-- ) of (# _, r #) -> fromBytes r +--infixl 5 <:> diff --git a/src/Numeric/Matrix/Base/FloatXNM.hs b/src/Numeric/Matrix/Base/FloatXNM.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Matrix/Base/FloatXNM.hs @@ -0,0 +1,565 @@ +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MagicHash, UnboxedTuples, DataKinds #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Matrix.Base.FloatXNM +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Matrix.Base.FloatXNM () where + +#include "MachDeps.h" +#include "HsBaseConfig.h" + +import GHC.Base (runRW#) +import GHC.Prim +import GHC.Types +import GHC.TypeLits + + +import Numeric.Commons +import Numeric.Matrix.Class +import Numeric.Matrix.Family + + + + +instance (KnownNat n, KnownNat m) => Show (MFloatXNM n m) where + show x@(MFloatXNM arr) = "{" ++ drop 2 (loop' (n -# 1#) (m -# 1#) " }") + where + loop' i j acc | isTrue# (i ==# -1#) = acc + | isTrue# (j ==# -1#) = loop' (i -# 1#) (m -# 1#) ('\n':acc) + | otherwise = loop' i (j -# 1#) (", " ++ show (F# (indexFloatArray# arr (i +# n *# j))) ++ acc) + n = dimN# x + m = dimM# x + + + +instance (KnownNat n, KnownNat m) => Eq (MFloatXNM n m) where + a == b = accumV2 (\x y r -> r && isTrue# (x `eqFloat#` y)) a b True + {-# INLINE (==) #-} + a /= b = accumV2 (\x y r -> r || isTrue# (x `neFloat#` y)) a b False + {-# INLINE (/=) #-} + + + + +-- | Implement partial ordering for `>`, `<`, `>=`, `<=` and lexicographical ordering for `compare` +instance (KnownNat n, KnownNat m) => Ord (MFloatXNM n m) where + a > b = accumV2 (\x y r -> r && isTrue# (x `gtFloat#` y)) a b True + {-# INLINE (>) #-} + a < b = accumV2 (\x y r -> r && isTrue# (x `ltFloat#` y)) a b True + {-# INLINE (<) #-} + a >= b = accumV2 (\x y r -> r && isTrue# (x `geFloat#` y)) a b True + {-# INLINE (>=) #-} + a <= b = accumV2 (\x y r -> r && isTrue# (x `leFloat#` y)) a b True + {-# INLINE (<=) #-} + -- | Compare lexicographically + compare a b = accumV2 (\x y r -> r `mappend` + if isTrue# (x `gtFloat#` y) + then GT + else if isTrue# (x `ltFloat#` y) + then LT + else EQ + ) a b EQ + {-# INLINE compare #-} + + + + +instance (KnownNat n, KnownNat m) => Num (MFloatXNM n m) where + (+) = zipV plusFloat# + {-# INLINE (+) #-} + (-) = zipV minusFloat# + {-# INLINE (-) #-} + (*) = zipV timesFloat# + {-# INLINE (*) #-} + negate = mapV negateFloat# + {-# INLINE negate #-} + abs = mapV (\x -> if isTrue# (x `geFloat#` 0.0#) then x else negateFloat# x) + {-# INLINE abs #-} + signum = mapV (\x -> if isTrue# (x `gtFloat#` 0.0#) then 1.0# else if isTrue# (x `ltFloat#` 0.0#) then -1.0# else 0.0#) + {-# INLINE signum #-} + fromInteger = broadcastMat . fromInteger + {-# INLINE fromInteger #-} + + + +instance (KnownNat n, KnownNat m) => Fractional (MFloatXNM n m) where + (/) = zipV divideFloat# + {-# INLINE (/) #-} + recip = mapV (divideFloat# 1.0#) + {-# INLINE recip #-} + fromRational = broadcastMat . fromRational + {-# INLINE fromRational #-} + + + +instance (KnownNat n, KnownNat m) => Floating (MFloatXNM n m) where + pi = broadcastMat pi + {-# INLINE pi #-} + exp = mapV expFloat# + {-# INLINE exp #-} + log = mapV logFloat# + {-# INLINE log #-} + sqrt = mapV sqrtFloat# + {-# INLINE sqrt #-} + sin = mapV sinFloat# + {-# INLINE sin #-} + cos = mapV cosFloat# + {-# INLINE cos #-} + tan = mapV tanFloat# + {-# INLINE tan #-} + asin = mapV asinFloat# + {-# INLINE asin #-} + acos = mapV acosFloat# + {-# INLINE acos #-} + atan = mapV atanFloat# + {-# INLINE atan #-} + sinh = mapV sinFloat# + {-# INLINE sinh #-} + cosh = mapV coshFloat# + {-# INLINE cosh #-} + tanh = mapV tanhFloat# + {-# INLINE tanh #-} + (**) = zipV powerFloat# + {-# INLINE (**) #-} + + logBase = zipV (\x y -> logFloat# y `divideFloat#` logFloat# x) + {-# INLINE logBase #-} + asinh = mapV (\x -> logFloat# (x `plusFloat#` sqrtFloat# (1.0# `plusFloat#` timesFloat# x x))) + {-# INLINE asinh #-} + acosh = mapV (\x -> case plusFloat# x 1.0# of + y -> logFloat# ( x `plusFloat#` timesFloat# y (sqrtFloat# (minusFloat# x 1.0# `divideFloat#` y))) + ) + {-# INLINE acosh #-} + atanh = mapV (\x -> 0.5# `timesFloat#` logFloat# (plusFloat# 1.0# x `divideFloat#` minusFloat# 1.0# x)) + {-# INLINE atanh #-} + + + + + + + + +instance (KnownNat n, KnownNat m) => MatrixCalculus Float n m (MFloatXNM n m) where + broadcastMat (F# x) = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop# n + (\i s' -> writeFloatArray# marr i x s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> MFloatXNM r + where + n = dimN# (undefined :: MFloatXNM n m) *# dimM# (undefined :: MFloatXNM n m) + bs = n *# SIZEOF_HSFLOAT# + {-# INLINE broadcastMat #-} + indexMat (I# i) (I# j) x@(MFloatXNM arr) +#ifndef UNSAFE_INDICES + | isTrue# ( (i ># n) + `orI#` (i <=# 0#) + `orI#` (j ># _m) + `orI#` (j <=# 0#) + ) = error $ "Bad index (" ++ show (I# i) ++ ", " ++ show (I# j) ++ ") for " + ++ show (I# n) ++ "x" ++ show (I# _m) ++ "D matrix" + | otherwise +#endif + = F# (indexFloatArray# arr (i -# 1# +# n *# (j -# 1#))) + where + n = dimN# x + _m = dimM# x + {-# INLINE indexMat #-} + transpose x@(MFloatXNM arr) = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop2# n m + (\i j s' -> writeFloatArray# marr (i +# n *# j) (indexFloatArray# arr (i *# m +# j)) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + n = dimN# x + m = dimM# x + bs = n *# m *# SIZEOF_HSFLOAT# + dimN x = I# (dimN# x) + {-# INLINE dimN #-} + dimM x = I# (dimM# x) + {-# INLINE dimM #-} + indexCol (I# j) x@(MFloatXNM arr) +#ifndef UNSAFE_INDICES + | isTrue# ( (j ># dimM# x) + `orI#` (j <=# 0#) + ) = error $ "Bad column index " ++ show (I# j) ++ " for " + ++ show (I# n) ++ "x" ++ show (dimM x) ++ "D matrix" + | otherwise +#endif + = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case copyByteArray# arr (bs *# (j -# 1#)) marr 0# bs s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + n = dimN# x + bs = n *# SIZEOF_HSFLOAT# + indexRow (I# i) x@(MFloatXNM arr) +#ifndef UNSAFE_INDICES + | isTrue# ( (i ># n) + `orI#` (i <=# 0#) + ) = error $ "Bad row index " ++ show (I# i) ++ " for " + ++ show (I# n) ++ "x" ++ show (I# m) ++ "D matrix" + | otherwise +#endif + = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop# m + (\j s' -> writeFloatArray# marr j (indexFloatArray# arr (i -# 1# +# n *# j)) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + n = dimN# x + m = dimM# x + bs = m *# SIZEOF_HSFLOAT# + + +instance (KnownNat n, KnownNat m) => PrimBytes (MFloatXNM n m) where + toBytes (MFloatXNM a) = a + {-# INLINE toBytes #-} + fromBytes = MFloatXNM + {-# INLINE fromBytes #-} + byteSize x = SIZEOF_HSFLOAT# *# dimN# x *# dimM# x + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_HSFLOAT# + {-# INLINE byteAlign #-} + +instance FloatBytes (MFloatXNM n m) where + ixF i (MFloatXNM a) = indexFloatArray# a i + {-# INLINE ixF #-} + + +instance KnownNat n => SquareMatrixCalculus Float n (MFloatXNM n n) where + eye = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop# n + (\j s' -> writeFloatArray# marr (j *# n1) 1.0# s' + ) (setByteArray# marr 0# bs 0# s1) of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + n1 = n +# 1# + n = dimN# (undefined :: MFloatXNM n n) + bs = n *# n *# SIZEOF_HSFLOAT# + {-# INLINE eye #-} + diag (F# v) = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop# n + (\j s' -> writeFloatArray# marr (j *# n1) v s' + ) (setByteArray# marr 0# bs 0# s1) of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + n1 = n +# 1# + n = dimN# (undefined :: MFloatXNM n n) + bs = n *# n *# SIZEOF_HSFLOAT# + {-# INLINE diag #-} + + + det v@(MFloatXNM arr) = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, mat #) -> case newByteArray# (n *# SIZEOF_HSFLOAT#) (copyByteArray# arr 0# mat 0# bs s1) of + (# s2, vec #) -> + let f i x s | isTrue# (i >=# n) = (# s, x #) + | otherwise = + let (# s' , j #) = maxInRowRem# n n i mat s + (# s'', x' #) = if isTrue# (i /=# j) then (# swapCols# n i j vec mat s', negateFloat# x #) + else (# s', x #) + (# s''', y #) = clearRowEnd# n n i mat s'' + in if isTrue# (eqFloat# 0.0# y) then (# s''', 0.0# #) + else f (i +# 1#) (timesFloat# x' y) s''' + in f 0# 1.0# s2 + ) of (# _, r #) -> F# r + where + n = dimN# v + bs = n *# n *# SIZEOF_HSFLOAT# + + {-# INLINE det #-} + + + + trace x@(MFloatXNM a) = F# (loop' 0# 0.0#) + where + n1 = n +# 1# + n = dimN# x + nn = n *# n + loop' i acc | isTrue# (i ># nn) = acc + | otherwise = loop' (i +# n1) (indexFloatArray# a i `plusFloat#` acc) + {-# INLINE trace #-} + fromDiag x@(MFloatXNM a) = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop# n + (\j s' -> writeFloatArray# marr j (indexFloatArray# a (j *# n1)) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + n1 = n +# 1# + n = dimN# x + bs = n *# SIZEOF_HSFLOAT# + {-# INLINE fromDiag #-} + toDiag x = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop# n + (\j s' -> writeFloatArray# marr (j *# n1) (indexFloatArray# a j) s' + ) (setByteArray# marr 0# bs 0# s1) of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + a = toBytes x + n1 = n +# 1# + n = dimN# (undefined :: MFloatXNM n n) + bs = n *# n *# SIZEOF_HSFLOAT# + {-# INLINE toDiag #-} + + +instance KnownNat n => MatrixInverse (MFloatXNM n n) where + inverse v@(MFloatXNM arr) = case runRW# + ( \s0 -> case newByteArray# (bs *# 2#) s0 of + (# s1, mat #) -> case newByteArray# (n *# SIZEOF_HSFLOAT#) + -- copy original matrix to the top of an augmented matrix + (loop# n (\i s -> writeFloatArray# mat (i *# nn +# i +# n) 1.0# (copyByteArray# arr (i *# vs) mat (2# *# i *# vs) vs s)) + (setByteArray# mat 0# (bs *# 2#) 0# s1) + ) of + (# s2, vec #) -> + let f i s | isTrue# (i >=# n) = s + | otherwise = + let (# s' , j #) = maxInRowRem# nn n i mat s + s'' = if isTrue# (i /=# j) then swapCols# nn i j vec mat s' + else s' + (# s''', _ #) = clearRowAll# nn n i mat s'' + in f (i +# 1#) s''' + in unsafeFreezeByteArray# mat + ( shrinkMutableByteArray# mat bs + (-- copy inverse matrix from the augmented part + loop# n (\i s -> copyMutableByteArray# mat (2# *# i *# vs +# vs) mat (i *# vs) vs s) + (f 0# s2) + ) + ) + ) of (# _, r #) -> MFloatXNM r + where + nn = 2# *# n + n = dimN# v + vs = n *# SIZEOF_HSFLOAT# + bs = n *# n *# SIZEOF_HSFLOAT# + + +----------------------------------------------------------------------------- +-- Helpers +----------------------------------------------------------------------------- + + +proxyN# :: MFloatXNM n m -> Proxy# n +proxyN# _ = proxy# + +dimN# :: KnownNat n => MFloatXNM n m -> Int# +dimN# x = case fromInteger (natVal' (proxyN# x)) of I# n -> n +{-# INLINE dimN# #-} + +dimM# :: KnownNat m => MFloatXNM n m -> Int# +dimM# x = case fromInteger (natVal x) of I# n -> n +{-# INLINE dimM# #-} + +-- | Do something in a loop for int i from 0 to n +loop# :: Int# -> (Int# -> State# s -> State# s) -> State# s -> State# s +loop# n f = loop' 0# + where + loop' i s | isTrue# (i ==# n) = s + | otherwise = case f i s of s1 -> loop' (i +# 1#) s1 +{-# INLINE loop# #-} + +-- | 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 +loop2# n m f = loop' 0# 0# + where + loop' i j s | isTrue# (j ==# m) = s + | isTrue# (i ==# n) = loop' 0# (j +# 1#) s + | otherwise = case f i j s of s1 -> loop' (i +# 1#) j s1 +{-# INLINE loop2# #-} + +zipV :: (KnownNat n, KnownNat m) => (Float# -> Float# -> Float#) -> MFloatXNM n m -> MFloatXNM n m -> MFloatXNM n m +zipV f x@(MFloatXNM a) (MFloatXNM b) = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop# n + (\i s' -> case f (indexFloatArray# a i) (indexFloatArray# b i) of + r -> writeFloatArray# marr i r s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> MFloatXNM r + where + n = dimN# x *# dimM# x + bs = n *# SIZEOF_HSFLOAT# + +mapV :: (KnownNat n, KnownNat m) => (Float# -> Float#) -> MFloatXNM n m -> MFloatXNM n m +mapV f x@(MFloatXNM a) = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop# n + (\i s' -> case f (indexFloatArray# a i) of + r -> writeFloatArray# marr i r s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> MFloatXNM r + where + n = dimN# x *# dimM# x + bs = n *# SIZEOF_HSFLOAT# + + +--accumVFloat :: (KnownNat n, KnownNat m) => (Float# -> Float# -> Float#) -> MFloatXNM n m -> Float# -> Float# +--accumVFloat f x@(MFloatXNM a) = loop' 0# +-- where +-- loop' i acc | isTrue# (i ==# n) = acc +-- | otherwise = loop' (i +# 1#) (f (indexFloatArray# a i) acc) +-- n = dimN# x *# dimM# x + +accumV2 :: (KnownNat n, KnownNat m) => (Float# -> Float# -> a -> a) -> MFloatXNM n m -> MFloatXNM n m -> a -> a +accumV2 f x@(MFloatXNM a) (MFloatXNM b) = loop' 0# + where + loop' i acc | isTrue# (i ==# n) = acc + | otherwise = loop' (i +# 1#) (f (indexFloatArray# a i) (indexFloatArray# b i) acc) + n = dimN# x *# dimM# x + + + +-- | Swap columns i and j. Does not check if i or j is larger than matrix width m +swapCols# :: Int# -- n + -> Int# -- ith column to swap + -> Int# -- jth column to swap + -> MutableByteArray# s -- buffer byte array of length of n elems + -> MutableByteArray# s -- byte array of matrix + -> State# s -- previous state + -> State# s -- next state +swapCols# n i j vec mat s0 = + -- copy ith column to bugger vec + case copyMutableByteArray# mat (i *# bs) vec 0# bs s0 of + s1 -> case copyMutableByteArray# mat (j *# bs) mat (i *# bs) bs s1 of + s2 -> copyMutableByteArray# vec 0# mat (j *# bs) bs s2 + where + bs = n *# SIZEOF_HSFLOAT# + +-- | Starting from i-th row and i+1-th column, substract a multiple of i-th column from i+1 .. m columns, +-- such that there are only zeroes in i-th row and i+1..m columns elements. +clearRowEnd# :: Int# -- n + -> Int# -- m + -> Int# -- ith column to remove from all others + -> MutableByteArray# s -- byte array of matrix + -> State# s -- previous state + -> (# State# s, Float# #) -- next state and a diagonal element +clearRowEnd# n m i mat s0 = (# loop' (i +# 1#) s1, y' #) + where + y0 = (n +# 1#) *# i +# 1# -- first element in source column + (# s1, y' #) = readFloatArray# mat ((n +# 1#) *# i) s0 -- diagonal element, must be non-zero + yrc = 1.0# `divideFloat#` y' + n' = n -# i -# 1# + loop' k s | isTrue# (k >=# m) = s + | otherwise = loop' (k +# 1#) + ( let x0 = k *# n +# i + (# s', a' #) = readFloatArray# mat x0 s + s'' = writeFloatArray# mat x0 0.0# s' + a = a' `timesFloat#` yrc + in multNRem# n' (x0 +# 1#) y0 a mat s'' + ) + +-- | Substract a multiple of i-th column from 0 .. i-1 and i+1 .. m columns, +-- such that there are only zeroes in i-th row and i+1..m columns elements. +-- Assuming that elements 0..i-1 in i-th row are zeroes, so they do not affect other columns. +-- After all columns updated, divide i-th row by its diagonal element +clearRowAll# :: Int# -- n + -> Int# -- m + -> Int# -- ith column to remove from all others + -> MutableByteArray# s -- byte array of matrix + -> State# s -- previous state + -> (# State# s, Float# #) -- next state and a diagonal element +clearRowAll# n m i mat s0 = (# divLoop (i +# 1#) (writeFloatArray# mat ((n +# 1#) *# i) 1.0# (loop' 0# i (loop' (i +# 1#) m s1))), y' #) + where + y0 = (n +# 1#) *# i +# 1# -- first element in source column + (# s1, y' #) = readFloatArray# mat ((n +# 1#) *# i) s0 -- diagonal element, must be non-zero + yrc = 1.0# `divideFloat#` y' + n' = n -# i -# 1# + loop' k km s | isTrue# (k >=# km) = s + | otherwise = loop' (k +# 1#) km + ( let x0 = k *# n +# i + (# s', a' #) = readFloatArray# mat x0 s + s'' = writeFloatArray# mat x0 0.0# s' + a = a' `timesFloat#` yrc + in multNRem# n' (x0 +# 1#) y0 a mat s'' + ) + divLoop k s | isTrue# (k >=# n) = s + | otherwise = divLoop (k +# 1#) + ( let x0 = n *# i +# k + (# s', x #) = readFloatArray# mat x0 s + in writeFloatArray# mat x0 (timesFloat# x yrc) s' + ) + + +-- | Remove a multiple of one row from another one. +-- do: xi = xi - yi*a +multNRem# :: Int# -- n - nr of elements to go through + -> Int# -- start idx of y + -> Int# -- start idx of x + -> Float# -- multiplier a + -> MutableByteArray# s -- byte array of matrix + -> State# s -- previous state + -> State# s -- next state +multNRem# 0# _ _ _ _ s = s +multNRem# n x0 y0 a mat s = multNRem# (n -# 1#) (x0 +# 1#) (y0 +# 1#) a mat + ( case readFloatArray# mat y0 s of + (# s1, y #) -> case readFloatArray# mat x0 s1 of + (# s2, x #) -> writeFloatArray# mat x0 (x `minusFloat#` timesFloat# y a) s2 + ) + + + +-- | Gives index of maximum (absolute) element in i-th row, starting from i-th element only. +-- If i >= m then returns i. +maxInRowRem# :: Int# -- n + -> Int# -- m + -> Int# -- ith column to start to search for and a row to look in + -> MutableByteArray# s -- byte array of matrix + -> State# s -- previous state + -> (# State# s, Int# #) -- next state +maxInRowRem# n m i mat s0 = loop' i (abs# v) i s1 + where + (# s1, v #) = readFloatArray# mat ((n +# 1#) *# i) s0 + abs# x = if isTrue# (x `geFloat#` 0.0#) then x else negateFloat# x + loop' ok ov k s | isTrue# (k >=# m) = (# s, ok #) + | otherwise = case readFloatArray# mat (n *# k +# i) s of + (# s', v' #) -> if isTrue# (abs# v' `gtFloat#` ov) + then loop' k (abs# v') (k +# 1#) s' + else loop' ok ov (k +# 1#) s' + + + + +--accumV2Float :: (KnownNat n, KnownNat m) => (Float# -> Float# -> Float# -> Float#) -> MFloatXNM n m -> MFloatXNM n m -> Float# -> Float# +--accumV2Float f x@(MFloatXNM a) (MFloatXNM b) = loop' 0# +-- where +-- loop' i acc | isTrue# (i ==# n) = acc +-- | otherwise = loop' (i +# 1#) (f (indexFloatArray# a i) (indexFloatArray# b i) acc) +-- n = dimN# x *# dimM# x +-- +-- +--accumVReverse :: (KnownNat n, KnownNat m) => (Float# -> a -> a) -> MFloatXNM n m -> a -> a +--accumVReverse f x@(MFloatXNM a) = loop' (n -# 1#) +-- where +-- loop' i acc | isTrue# (i ==# -1#) = acc +-- | otherwise = loop' (i -# 1#) (f (indexFloatArray# a i) acc) +-- n = dimN# x *# dimM# x 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,265 @@ +{-# LANGUAGE MagicHash, UnboxedTuples, DataKinds #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE MultiParamTypeClasses, DataKinds, KindSignatures #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Matrix.Class +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Matrix.Class + ( MatrixCalculus (..) + , SquareMatrixCalculus (..) + , Matrix2x2 (..) + , MatrixProduct (..) + , MatrixInverse (..) + , prodF, prodD + , prodI8, prodI16, prodI32, prodI64 + , prodW8, prodW16, prodW32, prodW64 + ) where + +import GHC.Base (runRW#) +import GHC.Prim +import GHC.Types + +#include "MachDeps.h" +#include "HsBaseConfig.h" + +import Numeric.Commons +import Numeric.Vector.Class +import Numeric.Vector.Family (Vector) +import Numeric.Matrix.Family (Matrix) + + +class MatrixCalculus t (n :: Nat) (m :: Nat) v | v -> t, v -> n, v -> m where + -- | Fill Mat with the same value + broadcastMat :: t -> v + -- | Get element by its index + indexMat :: Int -> Int -> v -> t + -- | Transpose Mat + transpose :: (MatrixCalculus t m n w, PrimBytes w) => v -> w + -- | First dimension size of a matrix + dimN :: v -> Int + -- | Second dimension size of a matrix + dimM :: v -> Int + -- | Get vector column by its index + indexCol :: (VectorCalculus t n w, PrimBytes w) => Int -> v -> w + -- | Get vector row by its index + indexRow :: (VectorCalculus t m w, PrimBytes w) => Int -> v -> w + +class SquareMatrixCalculus t (n :: Nat) v | v -> t, v -> n where + -- | Mat with 1 on diagonal and 0 elsewhere + eye :: v + -- | Put the same value on the Mat diagonal, 0 otherwise + diag :: t -> v + -- | Determinant of Mat + det :: v -> t + -- | Sum of diagonal elements + trace :: v -> t + -- | Get the diagonal elements from Mat into Vec + fromDiag :: (VectorCalculus t n w, PrimBytes w) => v -> w + -- | Set Vec values into the diagonal elements of Mat + toDiag :: (VectorCalculus t n w, PrimBytes w) => w -> v + + +class Matrix2x2 t where + -- | Compose a 2x2D matrix + mat22 :: Vector t 2 -> Vector t 2 -> Matrix t 2 2 + rowsOfM22 :: Matrix t 2 2 -> (Vector t 2, Vector t 2) + colsOfM22 :: Matrix t 2 2 -> (Vector t 2, Vector t 2) + + + +class MatrixProduct a b c where + -- | matrix-matrix or matrix-vector product + prod :: a -> b -> c + + +class MatrixInverse a where + inverse :: a -> a + + + +prodF :: (FloatBytes a, FloatBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodF n m k x y = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> + let loop' i j l r | isTrue# (l ==# m) = r + | otherwise = loop' i j (l +# 1#) (r `plusFloat#` timesFloat# (ixF (i +# n *# l) x) + (ixF (l +# m *# j) y)) + in case loop2# n k + (\i j s' -> writeFloatArray# marr (i +# n *# j) (loop' i j 0# 0.0#) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + bs = n *# k *# SIZEOF_HSFLOAT# +{-# INLINE prodF #-} + +prodD :: (DoubleBytes a, DoubleBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodD n m k x y= case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> + let loop' i j l r | isTrue# (l ==# m) = r + | otherwise = loop' i j (l +# 1#) (r +## (*##) (ixD (i +# n *# l) x) + (ixD (l +# m *# j) y)) + in case loop2# n k + (\i j s' -> writeDoubleArray# marr (i +# n *# j) (loop' i j 0# 0.0##) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + bs = n *# k *# SIZEOF_HSDOUBLE# +{-# INLINE prodD #-} + +prodI8 :: (IntBytes a, IntBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodI8 n m k x y= case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> + let loop' i j l r | isTrue# (l ==# m) = r + | otherwise = loop' i j (l +# 1#) (r +# (*#) (ixI (i +# n *# l) x) + (ixI (l +# m *# j) y)) + in case loop2# n k + (\i j s' -> writeInt8Array# marr (i +# n *# j) (loop' i j 0# 0#) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + bs = n *# k +{-# INLINE prodI8 #-} + + +prodI16 :: (IntBytes a, IntBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodI16 n m k x y= case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> + let loop' i j l r | isTrue# (l ==# m) = r + | otherwise = loop' i j (l +# 1#) (r +# (*#) (ixI (i +# n *# l) x) + (ixI (l +# m *# j) y)) + in case loop2# n k + (\i j s' -> writeInt16Array# marr (i +# n *# j) (loop' i j 0# 0#) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + bs = n *# k +{-# INLINE prodI16 #-} + + +prodI32 :: (IntBytes a, IntBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodI32 n m k x y= case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> + let loop' i j l r | isTrue# (l ==# m) = r + | otherwise = loop' i j (l +# 1#) (r +# (*#) (ixI (i +# n *# l) x) + (ixI (l +# m *# j) y)) + in case loop2# n k + (\i j s' -> writeInt32Array# marr (i +# n *# j) (loop' i j 0# 0#) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + bs = n *# k +{-# INLINE prodI32 #-} + + +prodI64 :: (IntBytes a, IntBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodI64 n m k x y= case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> + let loop' i j l r | isTrue# (l ==# m) = r + | otherwise = loop' i j (l +# 1#) (r +# (*#) (ixI (i +# n *# l) x) + (ixI (l +# m *# j) y)) + in case loop2# n k + (\i j s' -> writeInt64Array# marr (i +# n *# j) (loop' i j 0# 0#) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + bs = n *# k +{-# INLINE prodI64 #-} + + +prodW8 :: (WordBytes a, WordBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodW8 n m k x y = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> + let loop' i j l r | isTrue# (l ==# m) = r + | otherwise = loop' i j (l +# 1#) (r `plusWord#` timesWord# (ixW (i +# n *# l) x) + (ixW (l +# m *# j) y)) + in case loop2# n k + (\i j s' -> writeWord8Array# marr (i +# n *# j) (loop' i j 0# 0##) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + bs = n *# k +{-# INLINE prodW8 #-} + + +prodW16 :: (WordBytes a, WordBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodW16 n m k x y = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> + let loop' i j l r | isTrue# (l ==# m) = r + | otherwise = loop' i j (l +# 1#) (r `plusWord#` timesWord# (ixW (i +# n *# l) x) + (ixW (l +# m *# j) y)) + in case loop2# n k + (\i j s' -> writeWord16Array# marr (i +# n *# j) (loop' i j 0# 0##) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + bs = n *# k +{-# INLINE prodW16 #-} + +prodW32 :: (WordBytes a, WordBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodW32 n m k x y = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> + let loop' i j l r | isTrue# (l ==# m) = r + | otherwise = loop' i j (l +# 1#) (r `plusWord#` timesWord# (ixW (i +# n *# l) x) + (ixW (l +# m *# j) y)) + in case loop2# n k + (\i j s' -> writeWord32Array# marr (i +# n *# j) (loop' i j 0# 0##) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + bs = n *# k +{-# INLINE prodW32 #-} + +prodW64 :: (WordBytes a, WordBytes b, PrimBytes c) => Int# -> Int# -> Int# -> a -> b -> c +prodW64 n m k x y = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> + let loop' i j l r | isTrue# (l ==# m) = r + | otherwise = loop' i j (l +# 1#) (r `plusWord#` timesWord# (ixW (i +# n *# l) x) + (ixW (l +# m *# j) y)) + in case loop2# n k + (\i j s' -> writeWord64Array# marr (i +# n *# j) (loop' i j 0# 0##) s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> fromBytes r + where + bs = n *# k +{-# INLINE prodW64 #-} + + +-- | 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 +loop2# n m f = loop' 0# 0# + where + loop' i j s | isTrue# (j ==# m) = s + | 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/Numeric/Matrix/Family.hs b/src/Numeric/Matrix/Family.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Matrix/Family.hs @@ -0,0 +1,29 @@ +{-# LANGUAGE KindSignatures, DataKinds, TypeFamilyDependencies, MagicHash #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Matrix.Family +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Matrix.Family + ( Matrix + , MFloatXNM (..) + ) where + +import GHC.TypeLits +import GHC.Prim + + + +-- | Family of all matrix types, specialized on low-dimensional matrices +type family Matrix t (n :: Nat) (m :: Nat) = v | v -> t n m where + Matrix Float n m = MFloatXNM n m + + +-- | NxM D matrix +data MFloatXNM (n::Nat) (m::Nat) = MFloatXNM ByteArray# diff --git a/src/Numeric/Tensor.hs b/src/Numeric/Tensor.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Tensor.hs @@ -0,0 +1,290 @@ +{-# LANGUAGE UndecidableInstances #-} + +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE TypeFamilies, TypeFamilyDependencies #-} +{-# LANGUAGE MultiParamTypeClasses, MagicHash #-} +{-# LANGUAGE KindSignatures, DataKinds #-} +{-# LANGUAGE TypeOperators, FlexibleInstances, ScopedTypeVariables #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Tensor +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Tensor where + +import GHC.TypeLits +import GHC.Prim +import Data.Proxy + + + +data Dim (ds :: [Nat]) = Dim + +class Dimensions (ds :: [Nat]) where + dims :: Dim ds -> [Int] + + + +headDim :: Dim (d ': ds) -> Proxy d +headDim _ = Proxy + +tailDim :: Dim (d ': ds) -> Dim ds +tailDim _ = Dim + + +instance Dimensions '[] where + dims _ = [] + +instance (KnownNat d, Dimensions ds) => Dimensions (d ': ds) where + dims x = (fromIntegral . natVal $ headDim x) : dims (tailDim x) + + +printCrazy :: Dimensions d => Dim d -> String +printCrazy d = show $ dims d + +contraDimsType :: Tensor t n m -> Dim n +contraDimsType _ = Dim + +coDimsType :: Tensor t n m -> Dim m +coDimsType _ = Dim + +contraDims :: Dimensions n => Tensor t n m -> [Int] +contraDims = dims . contraDimsType + +coDims :: Dimensions m => Tensor t n m -> [Int] +coDims = dims . coDimsType + + +type Vec2 = Tensor Double '[2] '[] + + +vec2 :: Double -> Double -> Vec2 +vec2 x y = T10 $ Vector2 x y + + +class TensorCalculus t (ns :: [Nat]) (ms :: [Nat]) where + data Tensor t ns ms + type TensorStore t ns ms +-- -- | Add a contravariant rank +-- infixr 5 .<. +-- (.<.) :: Tensor t ns ms -> Tensor t ns ms -> Tensor t (2 ': ns) ms +-- -- | Add a covariant rank +-- infixr 5 .>. +-- (.>.) :: Tensor t ns ms -> Tensor t ns ms -> Tensor t ns (2 ': ms) +-- -- | Append dimension of the first contravariant rank +-- infixr 5 .< +-- (.<) :: AppendDim (TensorStore t ns ms) (TensorStore t (nb ': ns) ms) (TensorStore t ((nb + 1) ': ns) ms) +-- => Tensor t ns ms -> Tensor t (nb ': ns) ms -> Tensor t ((nb + 1) ': ns) ms +-- -- | Append dimension of the first covariant rank +-- infixr 5 .> +-- (.>) :: AppendDim (TensorStore t ns ms) (TensorStore t ns (mb ': ms)) (TensorStore t ns ((mb + 1) ': ms)) +-- => Tensor t ns ms -> Tensor t ns (mb ': ms) -> Tensor t ns ((mb + 1) ': ms) + +-- AppendDim (Tensor t ns ms) (Tensor t ns (mb ': ms)) (Tensor t ns ((mb + 1) ': ms)) + + +instance TensorCalculus t '[] '[] where + newtype Tensor t '[] '[] = T00 t deriving (Bounded, Enum, Eq, Integral, Num, Fractional, Floating, Ord, Read, Real, RealFrac, RealFloat, Show) + type TensorStore t '[] '[] = t +-- T00 a .<. T00 b = T10 $ Vector2 a b +-- T00 a .>. T00 b = T01 $ Vector2 a b +-- T00 a .< T10 b = T10 $ appendDim a b +-- T00 a .> T01 b = T01 $ appendDim a b +instance TensorCalculus t '[n] '[] where + newtype Tensor t '[n] '[] = T10 (SomeVector t n) + type TensorStore t '[n] '[] = SomeVector t n +-- contraV (T00 a) (T00 b) = T10 $ Vector2 a b +instance TensorCalculus t '[] '[m] where + newtype Tensor t '[] '[m] = T01 (SomeVector t m) + type TensorStore t '[] '[m] = SomeVector t m +instance TensorCalculus t '[n1, n2] '[] where + newtype Tensor t '[n1, n2] '[] = T20 (SomeMatrix t n1 n2) + type TensorStore t '[n1, n2] '[] = SomeMatrix t n1 n2 +instance TensorCalculus t '[n] '[m] where + newtype Tensor t '[n] '[m] = T11 (SomeMatrix t n m) + type TensorStore t '[n] '[m] = SomeMatrix t n m +instance TensorCalculus t '[] '[m1,m2] where + newtype Tensor t '[] '[m1,m2] = T02 (SomeMatrix t m1 m2) + type TensorStore t '[] '[m1,m2] = SomeMatrix t m1 m2 +instance TensorCalculus t (n1 ': n2 ': n3 ': ns) '[] where + newtype Tensor t (n1 ': n2 ': n3 ': ns) '[] = Tn0 (NDArray t) deriving Show + type TensorStore t (n1 ': n2 ': n3 ': ns) '[] = NDArray t +instance TensorCalculus t (n1 ': n2 ': n3 ': ns) '[m0] where + newtype Tensor t (n1 ': n2 ': n3 ': ns) '[m0] = Tn1 (NDArray t) deriving Show + type TensorStore t (n1 ': n2 ': n3 ': ns) '[m0] = NDArray t +instance TensorCalculus t (n1 ': n2 ': n3 ': ns) '[m0,m1] where + newtype Tensor t (n1 ': n2 ': n3 ': ns) '[m0, m1] = Tn2 (NDArray t) deriving Show + type TensorStore t (n1 ': n2 ': n3 ': ns) '[m0, m1] = NDArray t +instance TensorCalculus t '[] (m1 ': m2 ': m3 ': ms) where + newtype Tensor t '[] (m1 ': m2 ': m3 ': ms) = T0m (NDArray t) deriving Show + type TensorStore t '[] (m1 ': m2 ': m3 ': ms) = NDArray t +instance TensorCalculus t '[n1] (m1 ': m2 ': m3 ': ms) where + newtype Tensor t '[n1] (m1 ': m2 ': m3 ': ms) = T1m (NDArray t) deriving Show + type TensorStore t '[n1] (m1 ': m2 ': m3 ': ms) = NDArray t +instance TensorCalculus t '[n1, n2] (m1 ': m2 ': m3 ': ms) where + newtype Tensor t '[n1, n2] (m1 ': m2 ': m3 ': ms) = T2m (NDArray t) deriving Show + type TensorStore t '[n1, n2] (m1 ': m2 ': m3 ': ms) = NDArray t +instance TensorCalculus t (n1 ': n2 ': n3 ': ns) (m1 ': m2 ': m3 ': ms) where + newtype Tensor t (n1 ': n2 ': n3 ': ns) (m1 ': m2 ': m3 ': ms) = Tnm (NDArray t) deriving Show + type TensorStore t (n1 ': n2 ': n3 ': ns) (m1 ': m2 ': m3 ': ms) = NDArray t + + + +deriving instance Show (SomeVector t n) => Show (Tensor t '[n] '[]) +deriving instance Plus (SomeVector t n) => Plus (Tensor t '[n] '[]) +deriving instance Show (SomeVector t m) => Show (Tensor t '[] '[m]) +deriving instance Plus (SomeVector t m) => Plus (Tensor t '[] '[m]) +deriving instance Num (SomeVector t n) => Num (Tensor t '[n] '[]) +deriving instance Show (SomeMatrix t n0 n1) => Show (Tensor t '[n0, n1] '[]) +deriving instance Show (SomeMatrix t n0 m0) => Show (Tensor t '[n0] '[m0]) +deriving instance Show (SomeMatrix t m0 m1) => Show (Tensor t '[] '[m0,m1]) + + +--contraV :: Tensor t ns ms -> Tensor t ns ms -> Tensor t (2 ':. ns) ms +--contraV + + +class Plus a where + plus :: a -> a -> a + +instance Num t => Plus (Vector1 t) where + plus (Vector1 a) (Vector1 b) = Vector1 (a+b) + +instance Num t => Plus (Vector2 t) where + plus (Vector2 a1 a2) (Vector2 b1 b2) = Vector2 (a1+b1) (a2+b2) + +instance Num t => Plus (Vector3 t) where + plus (Vector3 a1 a2 a3) (Vector3 b1 b2 b3) = Vector3 (a1+b1) (a2+b2) (a3+b3) + +instance Num t => Plus (Vector4 t) where + plus (Vector4 a1 a2 a3 a4) (Vector4 b1 b2 b3 b4) = Vector4 (a1+b1) (a2+b2) (a3+b3) (a4+b4) + +instance Num t => Plus (VectorN t n) where + plus (VectorN as) (VectorN bs) = VectorN $ zipWith (+) as bs + + + + +newtype Vector1 t = Vector1 t + deriving Show +data Vector2 t = Vector2 t t + deriving Show +data Vector3 t = Vector3 t t t + deriving Show +data Vector4 t = Vector4 t t t t + deriving Show +newtype VectorN t (n::Nat) = VectorN [t] + deriving Show + +newtype Matrix1x1 t = Matrix1x1 t + deriving Show +data Matrix2x2 t = Matrix2x2 t t t t + deriving Show +newtype MatrixNxM t (n::Nat) (m::Nat) = MatrixNxM [[t]] + deriving Show + + +data NDArray t = NDArray ByteArray# + +instance Show (NDArray t) where + show _ = "Big array" + +--data TensorContraOnly t (ns :: [Nat]) = TContra Addr# +--data TensorCoOnly t (ns :: [Nat]) = TCo Addr# +--data TensorCo t (ns :: [Nat]) (ms :: [Nat]) = T Addr# + +type family SomeVector t (n :: Nat) = v | v -> t n where + SomeVector t 1 = Vector1 t + SomeVector t 2 = Vector2 t + SomeVector t 3 = Vector3 t + SomeVector t 4 = Vector4 t + SomeVector t n = VectorN t n + + +type family SomeMatrix t (n :: Nat) (m :: Nat) = v | v -> t n m where + SomeMatrix t 1 1 = Matrix1x1 t + SomeMatrix t 2 2 = Matrix2x2 t + SomeMatrix t n m = MatrixNxM t n m + + +class AppendDim a b c | a b -> c where + appendDim :: a -> b -> c + +instance AppendDim t (Vector2 t) (Vector3 t) where + appendDim a (Vector2 b1 b2) = Vector3 a b1 b2 +instance AppendDim t (Vector3 t) (Vector4 t) where + appendDim a (Vector3 b1 b2 b3) = Vector4 a b1 b2 b3 +instance AppendDim t (Vector4 t) (VectorN t 5) where + appendDim a (Vector4 b1 b2 b3 b4) = VectorN [a,b1,b2,b3,b4] +instance (m ~ (n+1), 5 <= m, 4 <= n) => AppendDim t (VectorN t n) (VectorN t m) where + appendDim a (VectorN bs) = VectorN $ a : bs + + +--class VectorOps a b c | a b -> c where +-- appendVecs :: a -> b -> c +-- +-- +--instance VectorOps (Vector1 t) (Vector1 t) (Vector2 t) where +-- appendVecs (Vector1 a) (Vector1 b) = Vector2 a b +--instance VectorOps (Vector2 t) (Vector1 t) (Vector3 t) where +-- appendVecs (Vector2 a1 a2) (Vector1 b) = Vector3 a1 a2 b +--instance VectorOps (Vector3 t) (Vector1 t) (Vector4 t) where +-- appendVecs (Vector3 a1 a2 a3) (Vector1 b) = Vector4 a1 a2 a3 b +--instance VectorOps (Vector4 t) (Vector1 t) (VectorN t 5) where +-- appendVecs (Vector4 a1 a2 a3 a4) (Vector1 b) = VectorN [a1,a2,a3,a4,b] +--instance m ~ (n+1) => VectorOps (VectorN t n) (Vector1 t) (VectorN t m) where +-- appendVecs (VectorN as) (Vector1 b) = VectorN $ as ++ [b] +-- +--instance VectorOps (Vector1 t) (Vector2 t) (Vector3 t) where +-- appendVecs (Vector1 a) (Vector2 b1 b2) = Vector3 a b1 b2 +--instance VectorOps (Vector2 t) (Vector2 t) (Vector4 t) where +-- appendVecs (Vector2 a1 a2) (Vector2 b1 b2) = Vector4 a1 a2 b1 b2 +--instance VectorOps (Vector3 t) (Vector2 t) (VectorN t 5) where +-- appendVecs (Vector3 a1 a2 a3) (Vector2 b1 b2) = VectorN [a1,a2,a3,b1,b2] +--instance VectorOps (Vector4 t) (Vector2 t) (VectorN t 6) where +-- appendVecs (Vector4 a1 a2 a3 a4) (Vector2 b1 b2) = VectorN [a1,a2,a3,a4,b1,b2] +--instance m ~ (n+2) => VectorOps (VectorN t n) (Vector2 t) (VectorN t m) where +-- appendVecs (VectorN as) (Vector2 b1 b2) = VectorN $ as ++ [b1,b2] +-- +--instance VectorOps (Vector1 t) (Vector3 t) (Vector4 t) where +-- appendVecs (Vector1 a) (Vector3 b1 b2 b3) = Vector4 a b1 b2 b3 +--instance VectorOps (Vector2 t) (Vector3 t) (VectorN t 5) where +-- appendVecs (Vector2 a1 a2) (Vector3 b1 b2 b3) = VectorN [a1,a2,b1,b2,b3] +--instance VectorOps (Vector3 t) (Vector3 t) (VectorN t 6) where +-- appendVecs (Vector3 a1 a2 a3) (Vector3 b1 b2 b3) = VectorN [a1,a2,a3,b1,b2,b3] +--instance VectorOps (Vector4 t) (Vector3 t) (VectorN t 7) where +-- appendVecs (Vector4 a1 a2 a3 a4) (Vector3 b1 b2 b3) = VectorN [a1,a2,a3,a4,b1,b2,b3] +--instance m ~ (n+3) => VectorOps (VectorN t n) (Vector3 t) (VectorN t m) where +-- appendVecs (VectorN as) (Vector3 b1 b2 b3) = VectorN $ as ++ [b1,b2,b3] +-- +--instance VectorOps (Vector1 t) (Vector4 t) (VectorN t 5) where +-- appendVecs (Vector1 a) (Vector4 b1 b2 b3 b4) = VectorN [a,b1,b2,b3,b4] +--instance VectorOps (Vector2 t) (Vector4 t) (VectorN t 6) where +-- appendVecs (Vector2 a1 a2) (Vector4 b1 b2 b3 b4) = VectorN [a1,a2,b1,b2,b3,b4] +--instance VectorOps (Vector3 t) (Vector4 t) (VectorN t 7) where +-- appendVecs (Vector3 a1 a2 a3) (Vector4 b1 b2 b3 b4) = VectorN [a1,a2,a3,b1,b2,b3,b4] +--instance VectorOps (Vector4 t) (Vector4 t) (VectorN t 8) where +-- appendVecs (Vector4 a1 a2 a3 a4) (Vector4 b1 b2 b3 b4) = VectorN [a1,a2,a3,a4,b1,b2,b3,b4] +--instance m ~ (n+4) => VectorOps (VectorN t n) (Vector4 t) (VectorN t m) where +-- appendVecs (VectorN as) (Vector4 b1 b2 b3 b4) = VectorN $ as ++ [b1,b2,b3,b4] +-- +--instance k ~ (m+1) => VectorOps (Vector1 t) (VectorN t m) (VectorN t k) where +-- appendVecs (Vector1 a) (VectorN bs) = VectorN $ a : bs +--instance k ~ (m+2) => VectorOps (Vector2 t) (VectorN t m) (VectorN t 6) where +-- appendVecs (Vector2 a1 a2) (VectorN bs) = VectorN $ a1 : a2 : bs +--instance k ~ (m+3) => VectorOps (Vector3 t) (VectorN t m) (VectorN t 7) where +-- appendVecs (Vector3 a1 a2 a3) (VectorN bs) = VectorN $ a1 : a2 : a3 : bs +--instance k ~ (m+4) => VectorOps (Vector4 t) (VectorN t m) (VectorN t 8) where +-- appendVecs (Vector4 a1 a2 a3 a4) (VectorN bs) = VectorN $ a1 : a2 : a3 : a4 : bs +--instance k ~ (m+n) => VectorOps (VectorN t n) (VectorN t m) (VectorN t m) where +-- appendVecs (VectorN as) (VectorN bs) = VectorN $ as ++ bs + diff --git a/src/Numeric/Vector.hs b/src/Numeric/Vector.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Vector.hs @@ -0,0 +1,59 @@ +{-# LANGUAGE MagicHash, UnboxedTuples, DataKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeOperators #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Vector +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Vector + ( Vector + -- * Type classes + , VectorCalculus (..), Vector2D (..), Vector3D (..), Vector4D (..) + -- * Type abbreviations + , Vec2f, Vec3f, Vec4f + -- * Common functions on vectors + , (<:>) + ) where + +import GHC.Base (runRW#) +import GHC.Prim +import GHC.TypeLits + +import Numeric.Commons +import Numeric.Vector.Class (VectorCalculus(..), Vector2D (..), Vector3D (..), Vector4D (..)) +import Numeric.Vector.Family (Vector) + + +-- Import instances +import Numeric.Vector.Base.FloatX2 () +import Numeric.Vector.Base.FloatXN () + + +-- Type abbreviations + +type Vec2f = Vector Float 2 +type Vec3f = Vector Float 3 +type Vec4f = Vector Float 4 + + +-- | Append one vector to another, adding up their dimensionality +(<:>) :: ( PrimBytes (Vector t n) + , PrimBytes (Vector t m) + , PrimBytes (Vector t (n+m)) + ) + => Vector t n -> Vector t m -> Vector t (n + m) +a <:> b = case (# toBytes a, toBytes b, byteSize a, byteSize b #) of + (# arr1, arr2, n, m #) -> case runRW# + ( \s0 -> case newByteArray# (n +# m) s0 of + (# s1, marr #) -> case copyByteArray# arr1 0# marr 0# n s1 of + s2 -> case copyByteArray# arr2 0# marr n m s2 of + s3 -> unsafeFreezeByteArray# marr s3 + ) of (# _, r #) -> fromBytes r +infixl 5 <:> diff --git a/src/Numeric/Vector/Base/FloatX2.hs b/src/Numeric/Vector/Base/FloatX2.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Vector/Base/FloatX2.hs @@ -0,0 +1,257 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MagicHash, UnboxedTuples, DataKinds #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Vector.Base.FloatX2 +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Vector.Base.FloatX2 () where + +#include "MachDeps.h" +#include "HsBaseConfig.h" + +import GHC.Base (runRW#) +import GHC.Prim +import GHC.Types + +import Numeric.Commons +import Numeric.Vector.Class +import Numeric.Vector.Family + + + + +instance Show VFloatX2 where + show (VFloatX2 a1 a2) = "{ " ++ show (F# a1) + ++ ", " ++ show (F# a2) + ++ " }" + + + +instance Eq VFloatX2 where + VFloatX2 a1 a2 == VFloatX2 b1 b2 = isTrue# ( (a1 `eqFloat#` b1) + `andI#` (a2 `eqFloat#` b2) + ) + {-# INLINE (==) #-} + VFloatX2 a1 a2 /= VFloatX2 b1 b2 = isTrue# ( (a1 `neFloat#` b1) + `orI#` (a2 `neFloat#` b2) + ) + {-# INLINE (/=) #-} + + + +-- | Implement partial ordering for `>`, `<`, `>=`, `<=` and lexicographical ordering for `compare` +instance Ord VFloatX2 where + VFloatX2 a1 a2 > VFloatX2 b1 b2 = isTrue# ( (a1 `gtFloat#` b1) + `andI#` (a2 `gtFloat#` b2) + ) + {-# INLINE (>) #-} + VFloatX2 a1 a2 < VFloatX2 b1 b2 = isTrue# ( (a1 `ltFloat#` b1) + `andI#` (a2 `ltFloat#` b2) + ) + {-# INLINE (<) #-} + VFloatX2 a1 a2 >= VFloatX2 b1 b2 = isTrue# ( (a1 `geFloat#` b1) + `andI#` (a2 `geFloat#` b2) + ) + {-# INLINE (>=) #-} + VFloatX2 a1 a2 <= VFloatX2 b1 b2 = isTrue# ( (a1 `leFloat#` b1) + `andI#` (a2 `leFloat#` b2) + ) + {-# INLINE (<=) #-} + -- | Compare lexicographically + compare (VFloatX2 a1 a2) (VFloatX2 b1 b2) + | isTrue# (a1 `gtFloat#` b1) = GT + | isTrue# (a1 `ltFloat#` b1) = LT + | isTrue# (a2 `gtFloat#` b2) = GT + | isTrue# (a2 `ltFloat#` b2) = LT + | otherwise = EQ + {-# INLINE compare #-} + + + +-- | element-wise operations for vectors +instance Num VFloatX2 where + VFloatX2 a1 a2 + VFloatX2 b1 b2 = VFloatX2 (plusFloat# a1 b1) (plusFloat# a2 b2) + {-# INLINE (+) #-} + VFloatX2 a1 a2 - VFloatX2 b1 b2 = VFloatX2 (minusFloat# a1 b1) (minusFloat# a2 b2) + {-# INLINE (-) #-} + VFloatX2 a1 a2 * VFloatX2 b1 b2 = VFloatX2 (timesFloat# a1 b1) (timesFloat# a2 b2) + {-# INLINE (*) #-} + negate (VFloatX2 a1 a2) = VFloatX2 (negateFloat# a1) (negateFloat# a2) + {-# INLINE negate #-} + abs (VFloatX2 a1 a2) = VFloatX2 (if isTrue# (a1 `geFloat#` 0.0#) then a1 else negateFloat# a1) + (if isTrue# (a2 `geFloat#` 0.0#) then a2 else negateFloat# a2) + {-# INLINE abs #-} + signum (VFloatX2 a1 a2) = VFloatX2 (if isTrue# (a1 `gtFloat#` 0.0#) then 1.0# else if isTrue# (a1 `ltFloat#` 0.0#) then -1.0# else 0.0# ) + (if isTrue# (a2 `gtFloat#` 0.0#) then 1.0# else if isTrue# (a1 `ltFloat#` 0.0#) then -1.0# else 0.0# ) + {-# INLINE signum #-} + fromInteger n = case fromInteger n of F# x -> VFloatX2 x x + {-# INLINE fromInteger #-} + + + +instance Fractional VFloatX2 where + VFloatX2 a1 a2 / VFloatX2 b1 b2 = VFloatX2 (divideFloat# a1 b1) + (divideFloat# a2 b2) + {-# INLINE (/) #-} + recip (VFloatX2 a1 a2) = VFloatX2 (divideFloat# 1.0# a1) + (divideFloat# 1.0# a2) + {-# INLINE recip #-} + fromRational r = case fromRational r of F# x -> VFloatX2 x x + {-# INLINE fromRational #-} + + + +instance Floating VFloatX2 where + pi = VFloatX2 3.141592653589793238# 3.141592653589793238# + {-# INLINE pi #-} + exp (VFloatX2 a1 a2) = VFloatX2 (expFloat# a1) + (expFloat# a2) + {-# INLINE exp #-} + log (VFloatX2 a1 a2) = VFloatX2 (logFloat# a1) + (logFloat# a2) + {-# INLINE log #-} + sqrt (VFloatX2 a1 a2) = VFloatX2 (sqrtFloat# a1) + (sqrtFloat# a2) + {-# INLINE sqrt #-} + sin (VFloatX2 a1 a2) = VFloatX2 (sinFloat# a1) + (sinFloat# a2) + {-# INLINE sin #-} + cos (VFloatX2 a1 a2) = VFloatX2 (cosFloat# a1) + (cosFloat# a2) + {-# INLINE cos #-} + tan (VFloatX2 a1 a2) = VFloatX2 (tanFloat# a1) + (tanFloat# a2) + {-# INLINE tan #-} + asin (VFloatX2 a1 a2) = VFloatX2 (asinFloat# a1) + (asinFloat# a2) + {-# INLINE asin #-} + acos (VFloatX2 a1 a2) = VFloatX2 (acosFloat# a1) + (acosFloat# a2) + {-# INLINE acos #-} + atan (VFloatX2 a1 a2) = VFloatX2 (atanFloat# a1) + (atanFloat# a2) + {-# INLINE atan #-} + sinh (VFloatX2 a1 a2) = VFloatX2 (sinFloat# a1) + (sinFloat# a2) + {-# INLINE sinh #-} + cosh (VFloatX2 a1 a2) = VFloatX2 (coshFloat# a1) + (coshFloat# a2) + {-# INLINE cosh #-} + tanh (VFloatX2 a1 a2) = VFloatX2 (tanhFloat# a1) + (tanhFloat# a2) + {-# INLINE tanh #-} + VFloatX2 a1 a2 ** VFloatX2 b1 b2 = VFloatX2 (powerFloat# a1 b1) + (powerFloat# 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 #-} + + + +-- log1p (VFloatX2 a1 a2) = case ( log1p (F# a1), log1p (F# a2) ) of +-- (F# x1, F# x2) -> VFloatX2 x1 x2 +-- expm1 (VFloatX2 a1 a2) = case ( expm1 (F# a1), expm1 (F# a2) ) of +-- (F# x1, F# x2) -> VFloatX2 x1 x2 +-- +-- log1mexp a +-- | a <= log 2 = log (negate (expm1Float a)) +-- | otherwise = log1p (negate (exp a)) +-- {-# INLINE log1mexp #-} +-- log1pexp a +-- | a <= 18 = log1p (exp a) +-- | a <= 100 = a + exp (negate a) +-- | otherwise = a +-- {-# INLINE log1pexp #-} + + + +instance VectorCalculus Float 2 VFloatX2 where + broadcastVec (F# x) = VFloatX2 x x + {-# INLINE broadcastVec #-} + VFloatX2 a1 a2 .*. VFloatX2 b1 b2 = case timesFloat# a1 b1 + `plusFloat#` timesFloat# a2 b2 of + x -> VFloatX2 x x + {-# INLINE (.*.) #-} + VFloatX2 a1 a2 `dot` VFloatX2 b1 b2 = F# ( timesFloat# a1 b1 + `plusFloat#` timesFloat# a2 b2 + ) + {-# INLINE dot #-} + indexVec 1 (VFloatX2 a1 _) = F# a1 + indexVec 2 (VFloatX2 _ a2) = F# a2 + indexVec i _ = error $ "Bad index " ++ show i ++ " for 2D vector" + {-# INLINE indexVec #-} + normL1 v = case abs v of + VFloatX2 a1 a2 -> F# (a1 `plusFloat#` a2) + {-# INLINE normL1 #-} + normL2 v = sqrt $ dot v v + {-# INLINE normL2 #-} + normLPInf (VFloatX2 a1 a2) = F# (if isTrue# (a1 `gtFloat#` a2) then a1 else a2) + {-# INLINE normLPInf #-} + normLNInf (VFloatX2 a1 a2) = F# (if isTrue# (a1 `gtFloat#` a2) then a2 else a1) + {-# INLINE normLNInf #-} + normLP n (VFloatX2 a1 a2) = case realToFrac n of + F# x -> F# ( powerFloat# (divideFloat# 1.0# x) + ( powerFloat# a1 x + `plusFloat#` powerFloat# a2 x + ) + ) + {-# INLINE normLP #-} + dim _ = 2 + {-# INLINE dim #-} + + + + +instance Vector2D Float where + vec2 (F# x) (F# y) = VFloatX2 x y + {-# INLINE vec2 #-} + det2 (VFloatX2 a1 a2) (VFloatX2 b1 b2) = F# (timesFloat# a1 b2 `minusFloat#` timesFloat# a2 b1) + {-# INLINE det2 #-} + + +instance PrimBytes VFloatX2 where + toBytes (VFloatX2 a1 a2) = case runRW# + ( \s0 -> case newByteArray# (SIZEOF_HSFLOAT# *# 2#) s0 of + (# s1, marr #) -> case writeFloatArray# marr 0# a1 s1 of + s2 -> case writeFloatArray# marr 1# a2 s2 of + s3 -> unsafeFreezeByteArray# marr s3 + ) of (# _, a #) -> a + {-# INLINE toBytes #-} + fromBytes arr = VFloatX2 + (indexFloatArray# arr 0#) + (indexFloatArray# arr 1#) + {-# INLINE fromBytes #-} + byteSize _ = SIZEOF_HSFLOAT# *# 2# + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_HSFLOAT# + {-# INLINE byteAlign #-} + +instance FloatBytes VFloatX2 where + ixF 0# (VFloatX2 a1 _) = a1 + ixF 1# (VFloatX2 _ a2) = a2 + ixF _ _ = undefined + {-# INLINE ixF #-} + + + + + + + diff --git a/src/Numeric/Vector/Base/FloatXN.hs b/src/Numeric/Vector/Base/FloatXN.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Vector/Base/FloatXN.hs @@ -0,0 +1,268 @@ +{-# LANGUAGE TypeOperators, TypeFamilies #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MagicHash, UnboxedTuples, DataKinds #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Vector.Base.FloatXN +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Vector.Base.FloatXN () where + +#include "MachDeps.h" +#include "HsBaseConfig.h" + +import GHC.Base (runRW#) +import GHC.Prim +import GHC.Types +import GHC.TypeLits + +import Numeric.Commons +import Numeric.Vector.Class +import Numeric.Vector.Family + + +instance KnownNat n => Show (VFloatXN n) where + show x = "{" ++ drop 1 + ( accumVReverse (\a s -> ", " ++ show (F# a) ++ s) x " }" + ) + +instance KnownNat n => Eq (VFloatXN n) where + a == b = accumV2 (\x y r -> r && isTrue# (x `eqFloat#` y)) a b True + {-# INLINE (==) #-} + a /= b = accumV2 (\x y r -> r || isTrue# (x `neFloat#` y)) a b False + {-# INLINE (/=) #-} + +-- | Implement partial ordering for `>`, `<`, `>=`, `<=` and lexicographical ordering for `compare` +instance KnownNat n => Ord (VFloatXN n) where + a > b = accumV2 (\x y r -> r && isTrue# (x `gtFloat#` y)) a b True + {-# INLINE (>) #-} + a < b = accumV2 (\x y r -> r && isTrue# (x `ltFloat#` y)) a b True + {-# INLINE (<) #-} + a >= b = accumV2 (\x y r -> r && isTrue# (x `geFloat#` y)) a b True + {-# INLINE (>=) #-} + a <= b = accumV2 (\x y r -> r && isTrue# (x `leFloat#` y)) a b True + {-# INLINE (<=) #-} + -- | Compare lexicographically + compare a b = accumV2 (\x y r -> r `mappend` + if isTrue# (x `gtFloat#` y) + then GT + else if isTrue# (x `ltFloat#` y) + then LT + else EQ + ) a b EQ + {-# INLINE compare #-} + + + +instance (KnownNat n, 3 <= n) => Num (VFloatXN n) where + (+) = zipV plusFloat# + {-# INLINE (+) #-} + (-) = zipV minusFloat# + {-# INLINE (-) #-} + (*) = zipV timesFloat# + {-# INLINE (*) #-} + negate = mapV negateFloat# + {-# INLINE negate #-} + abs = mapV (\x -> if isTrue# (x `geFloat#` 0.0#) then x else negateFloat# x) + {-# INLINE abs #-} + signum = mapV (\x -> if isTrue# (x `gtFloat#` 0.0#) then 1.0# else if isTrue# (x `ltFloat#` 0.0#) then -1.0# else 0.0#) + {-# INLINE signum #-} + fromInteger = broadcastVec . fromInteger + {-# INLINE fromInteger #-} + + + +instance (KnownNat n, 3 <= n) => Fractional (VFloatXN n) where + (/) = zipV divideFloat# + {-# INLINE (/) #-} + recip = mapV (divideFloat# 1.0#) + {-# INLINE recip #-} + fromRational = broadcastVec . fromRational + {-# INLINE fromRational #-} + + + +instance (KnownNat n, 3 <= n) => Floating (VFloatXN n) where + pi = broadcastVec pi + {-# INLINE pi #-} + exp = mapV expFloat# + {-# INLINE exp #-} + log = mapV logFloat# + {-# INLINE log #-} + sqrt = mapV sqrtFloat# + {-# INLINE sqrt #-} + sin = mapV sinFloat# + {-# INLINE sin #-} + cos = mapV cosFloat# + {-# INLINE cos #-} + tan = mapV tanFloat# + {-# INLINE tan #-} + asin = mapV asinFloat# + {-# INLINE asin #-} + acos = mapV acosFloat# + {-# INLINE acos #-} + atan = mapV atanFloat# + {-# INLINE atan #-} + sinh = mapV sinFloat# + {-# INLINE sinh #-} + cosh = mapV coshFloat# + {-# INLINE cosh #-} + tanh = mapV tanhFloat# + {-# INLINE tanh #-} + (**) = zipV powerFloat# + {-# INLINE (**) #-} + + logBase = zipV (\x y -> logFloat# y `divideFloat#` logFloat# x) + {-# INLINE logBase #-} + asinh = mapV (\x -> logFloat# (x `plusFloat#` sqrtFloat# (1.0# `plusFloat#` timesFloat# x x))) + {-# INLINE asinh #-} + acosh = mapV (\x -> case plusFloat# x 1.0# of + y -> logFloat# ( x `plusFloat#` timesFloat# y (sqrtFloat# (minusFloat# x 1.0# `divideFloat#` y))) + ) + {-# INLINE acosh #-} + atanh = mapV (\x -> 0.5# `timesFloat#` logFloat# (plusFloat# 1.0# x `divideFloat#` minusFloat# 1.0# x)) + {-# INLINE atanh #-} + + + +instance (KnownNat n, 3 <= n) => VectorCalculus Float n (VFloatXN n) where + broadcastVec (F# x) = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop# n + (\i s' -> writeFloatArray# marr i x s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> VFloatXN r + where + n = dim# (undefined :: VFloatXN n) + bs = n *# 4# + {-# INLINE broadcastVec #-} + a .*. b = broadcastVec $ dot a b + {-# INLINE (.*.) #-} + dot a b = F# (accumV2Float (\x y r -> r `plusFloat#` timesFloat# x y) a b 0.0#) + {-# INLINE dot #-} + indexVec (I# i) _x@(VFloatXN arr) +#ifndef UNSAFE_INDICES + | isTrue# ( (i ># dim# _x) + `orI#` (i <=# 0#) + ) = error $ "Bad index " ++ show (I# i) ++ " for " ++ show (dim _x) ++ "D vector" + | otherwise +#endif + = F# (indexFloatArray# arr (i -# 1#)) + {-# INLINE indexVec #-} + normL1 v = F# (accumVFloat (\x a -> a `plusFloat#` (if isTrue# (x `geFloat#` 0.0#) then x else negateFloat# x)) v 0.0#) + {-# INLINE normL1 #-} + normL2 v = sqrt $ F# (accumVFloat (\x a -> a `plusFloat#` timesFloat# x x) v 0.0#) + {-# INLINE normL2 #-} + normLPInf v@(VFloatXN arr) = F# (accumVFloat (\x a -> if isTrue# (x `geFloat#` a) then x else a) v (indexFloatArray# arr 0#)) + {-# INLINE normLPInf #-} + normLNInf v@(VFloatXN arr) = F# (accumVFloat (\x a -> if isTrue# (x `leFloat#` a) then x else a) v (indexFloatArray# arr 0#)) + {-# INLINE normLNInf #-} + normLP n v = case realToFrac n of + F# p -> F# (powerFloat# (divideFloat# 1.0# p) (accumVFloat (\x a -> a `plusFloat#` powerFloat# x p) v 0.0#)) + {-# INLINE normLP #-} + dim x = I# (dim# x) + {-# INLINE dim #-} + + + + +instance KnownNat n => PrimBytes (VFloatXN n) where + toBytes (VFloatXN a) = a + {-# INLINE toBytes #-} + fromBytes = VFloatXN + {-# INLINE fromBytes #-} + byteSize x = SIZEOF_HSFLOAT# *# dim# x + {-# INLINE byteSize #-} + byteAlign _ = ALIGNMENT_HSFLOAT# + {-# INLINE byteAlign #-} + + +instance FloatBytes (VFloatXN n) where + ixF i (VFloatXN a) = indexFloatArray# a i + {-# INLINE ixF #-} + +----------------------------------------------------------------------------- +-- Helpers +----------------------------------------------------------------------------- + + +dim# :: KnownNat n => VFloatXN n -> Int# +dim# x = case fromInteger (natVal x) of I# n -> n +{-# INLINE dim# #-} + +-- | Do something in a loop for int i from 0 to n-1 +loop# :: Int# -> (Int# -> State# s -> State# s) -> State# s -> State# s +loop# n f = loop' 0# + where + loop' i s | isTrue# (i ==# n) = s + | otherwise = case f i s of s1 -> loop' (i +# 1#) s1 +{-# INLINE loop# #-} + + +zipV :: KnownNat n => (Float# -> Float# -> Float#) -> VFloatXN n -> VFloatXN n -> VFloatXN n +zipV f x@(VFloatXN a) (VFloatXN b) = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop# n + (\i s' -> case f (indexFloatArray# a i) (indexFloatArray# b i) of + r -> writeFloatArray# marr i r s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> VFloatXN r + where + n = dim# x + bs = n *# 4# + +mapV :: KnownNat n => (Float# -> Float#) -> VFloatXN n -> VFloatXN n +mapV f x@(VFloatXN a) = case runRW# + ( \s0 -> case newByteArray# bs s0 of + (# s1, marr #) -> case loop# n + (\i s' -> case f (indexFloatArray# a i) of + r -> writeFloatArray# marr i r s' + ) s1 of + s2 -> unsafeFreezeByteArray# marr s2 + ) of (# _, r #) -> VFloatXN r + where + n = dim# x + bs = n *# 4# + + +accumVFloat :: KnownNat n => (Float# -> Float# -> Float#) -> VFloatXN n -> Float# -> Float# +accumVFloat f x@(VFloatXN a) = loop' 0# + where + loop' i acc | isTrue# (i ==# n) = acc + | otherwise = loop' (i +# 1#) (f (indexFloatArray# a i) acc) + n = dim# x + +accumV2 :: KnownNat n => (Float# -> Float# -> a -> a) -> VFloatXN n -> VFloatXN n -> a -> a +accumV2 f x@(VFloatXN a) (VFloatXN b) = loop' 0# + where + loop' i acc | isTrue# (i ==# n) = acc + | otherwise = loop' (i +# 1#) (f (indexFloatArray# a i) (indexFloatArray# b i) acc) + n = dim# x + + +accumV2Float :: KnownNat n => (Float# -> Float# -> Float# -> Float#) -> VFloatXN n -> VFloatXN n -> Float# -> Float# +accumV2Float f x@(VFloatXN a) (VFloatXN b) = loop' 0# + where + loop' i acc | isTrue# (i ==# n) = acc + | otherwise = loop' (i +# 1#) (f (indexFloatArray# a i) (indexFloatArray# b i) acc) + n = dim# x + + +accumVReverse :: KnownNat n => (Float# -> a -> a) -> VFloatXN n -> a -> a +accumVReverse f x@(VFloatXN a) = loop' (n -# 1#) + where + loop' i acc | isTrue# (i ==# -1#) = acc + | otherwise = loop' (i -# 1#) (f (indexFloatArray# a i) acc) + n = dim# x diff --git a/src/Numeric/Vector/Class.hs b/src/Numeric/Vector/Class.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Vector/Class.hs @@ -0,0 +1,70 @@ +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE MultiParamTypeClasses, DataKinds, KindSignatures #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Vector.Class +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Vector.Class + ( VectorCalculus (..) + , Vector2D (..) + , Vector3D (..) + , Vector4D (..) + ) where + + +import GHC.TypeLits + +import Numeric.Vector.Family (Vector) + +class VectorCalculus t (n :: Nat) v | v -> t, v -> n where + -- | Fill Vec with the same value + broadcastVec :: t -> v + -- | Scalar product -- sum of Vecs' components products, propagated into whole Vec + infixl 7 .*. + (.*.) :: v -> v -> v + -- | Scalar product -- sum of Vecs' components products -- a scalar + dot :: v -> v -> t + -- | Get element by its index + indexVec :: Int -> v -> t + -- | Sum of absolute values + normL1 :: v -> t + -- | hypot function (square root of squares) + normL2 :: v -> t + -- | Maximum of absolute values + normLPInf :: v -> t + -- | Minimum of absolute values + normLNInf :: v -> t + -- | Norm in Lp space + normLP :: Int -> v -> t + -- | Dimensionality of a vector + dim :: v -> Int + +{-# RULES +"normLP/L1" normLP 1 = normL1 +"normLP/L2" normLP 2 = normL2 + #-} + +class Vector2D t where + -- | Compose a 2D vector + vec2 :: t -> t -> Vector t 2 + -- | Take a determinant of a matrix composed from two 2D vectors. + -- Like a cross product in 2D. + det2 :: Vector t 2 -> Vector t 2 -> t + + +class Vector3D t where + -- | Compose a 3D vector + vec3 :: t -> t -> t -> Vector t 3 + -- | Cross product + cross :: Vector t 3 -> Vector t 3 -> Vector t 3 + +class Vector4D t where + -- | Compose a 4D vector + vec4 :: t -> t -> t -> t -> Vector t 4 diff --git a/src/Numeric/Vector/Family.hs b/src/Numeric/Vector/Family.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Vector/Family.hs @@ -0,0 +1,32 @@ +{-# LANGUAGE KindSignatures, DataKinds, TypeFamilyDependencies, MagicHash #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Vector.Family +-- Copyright : (c) Artem Chirkin +-- License : MIT +-- +-- Maintainer : chirkin@arch.ethz.ch +-- +-- +----------------------------------------------------------------------------- + +module Numeric.Vector.Family + ( Vector + , VFloatX2 (..), VFloatXN (..) + ) where + +import GHC.TypeLits +import GHC.Prim + + +-- | Family of all vector types, specialized on low-dimensional vectors +type family Vector t (n :: Nat) = v | v -> t n where + Vector Float 2 = VFloatX2 + Vector Float n = VFloatXN n + + +-- | 2D Float vector +data VFloatX2 = VFloatX2 Float# Float# + +-- | ND vector +data VFloatXN (n::Nat) = VFloatXN ByteArray# diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,2 @@ +main :: IO () +main = putStrLn "Test suite not yet implemented"