diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,7 @@ +# Revision history for accelerate-typelit + +## 2016-04-23 -- 0.1.0.0 -- initial version + +* Defining data types to handle type-safe matrices and vectors +* Providing vector/matrix/scalar functions for addition, multiplication + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,13 @@ +Copyright (c) 2016 Martin Heuschober + +Permission to use, copy, modify, and/or distribute this software for any purpose +with or without fee is hereby granted, provided that the above copyright notice +and this permission notice appear in all copies. + +THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH +REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, +INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS +OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER +TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF +THIS SOFTWARE. diff --git a/Readme.md b/Readme.md new file mode 100644 --- /dev/null +++ b/Readme.md @@ -0,0 +1,69 @@ +Accelerate TypeLits +=================== + +[![Build Status](https://travis-ci.org/epsilonhalbe/accelerate-typelits.svg?branch=master)](https://travis-ci.org/epsilonhalbe/accelerate-typelits) +Synopsis +-------- + +This library provides a high level interface to `accelerate` for matrix +computations. + +Installation +------------ + +The simplest way to install this library is using `cabal` or `cabal-sandbox` + +``` +> cabal install accelerate-typelits +``` + +If you want to have the most recent version, the project is on github so you can +checkout the project. + +``` +> git clone https://github.com/epsilonhalbe/accelerate-typelits.git +> cd accelerate-typelits +> cabal install +``` + +--- + +``` +> git clone https://github.com/epsilonhalbe/accelerate-typelits.git +> cd accelerate-typelits +> cabal sandbox init +> cabal install +``` + +There is also a `stack.yaml` file included, so one can also use [stack][1] in +order to compile this library. + +``` +> git clone https://github.com/epsilonhalbe/accelerate-typelits.git +> cd accelerate-typelits +> stack --stack-yaml stack-7.10.yaml build +``` + +--- + +The operators have been designed to give a visual hint of the respective +parameters. + +- `#` for matrices +- `^` for vectors +- `.` for scalars + +So for example `#*^` represents the multiplication of a matrix with a vector, +analogously `^*#` works the other way around. Other examples would be `#*#` for +matrix-matrix multiplication and `.*^` scalar multiplication of a vector. + +Operator precedence is usually the same as the numeric equivalence. + +Credits +------- + +The matrix-vector and matrix-matrix products have been inspired by Henning +Thielemann's [`accelerate-arithmetic`][2] library + +[1]: https://haskellstack.com +[2]: https://hackage.haskell.org/package/accelerate-arithmetic 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/accelerate-typelits.cabal b/accelerate-typelits.cabal new file mode 100644 --- /dev/null +++ b/accelerate-typelits.cabal @@ -0,0 +1,74 @@ +--------------------------------------------------------------------------------- +-- -- --------------------------- -- -- +-- -- --- Accelerate TypeLits --- -- -- +-- -- --------------------------- -- -- +--------------------------------------------------------------------------------- + +name: accelerate-typelits +version: 0.1.0.0 +synopsis: a typesafe way encode accelerate matrices and vectors +description: a small wrapper plus convenience functions on top of + accelerate to represent matrices with their dimensions +stability: experimental +license: ISC +license-file: LICENSE +author: Martin Heuschober +maintainer: Martin Heuschober +bug-reports: http://github.com/epsilonhalbe/accelerate-typelit/issues +copyright: (c) 2016 Martin Heuschober +category: Math +build-type: Simple +extra-source-files: ChangeLog.md + , Readme.md + , stack-7.10.yaml +cabal-version: >=1.22 +tested-with: GHC == 7.10.3 + +source-repository head + type: git + location: git://github.com/epsilonhalbe/accelerate-typelit + +--------------------------------------------------------------------------------- +-- Library -- +--------------------------------------------------------------------------------- + +library + exposed-modules: Data.Array.Accelerate.TypeLits + , Data.Array.Accelerate.TypeLits.System.Random.MWC + other-modules: Data.Array.Accelerate.TypeLits.Internal + build-depends: base >=4.8 && <4.9 + , accelerate + , accelerate-random + , mwc-random + , smallcheck + , QuickCheck + + hs-source-dirs: src + default-language: Haskell2010 + ghc-options: -Wall + +--------------------------------------------------------------------------------- +-- Tests -- +--------------------------------------------------------------------------------- + +test-suite tests + type: exitcode-stdio-1.0 + main-is: Test.hs + --exposed-modules: + other-modules: Test.Data.Array.Accelerate.TypeLits + build-depends: base >=4.8 && <4.9 + , accelerate + , accelerate-random + , accelerate-typelits + , HUnit-Plus + , mwc-random + , QuickCheck + , smallcheck + , tasty + , tasty-hunit + , tasty-quickcheck + , tasty-smallcheck + + hs-source-dirs: test + default-language: Haskell2010 + diff --git a/src/Data/Array/Accelerate/TypeLits.hs b/src/Data/Array/Accelerate/TypeLits.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/TypeLits.hs @@ -0,0 +1,356 @@ +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE PartialTypeSignatures #-} +module Data.Array.Accelerate.TypeLits + ( + -- * Types + AccScalar, + AccVector, + AccMatrix, + -- * Classes + AccFunctor(..), + -- * Constructors + mkMatrix, + mkVector, + mkScalar, + unsafeMkMatrix, + unsafeMkVector, + unMatrix, + unVector, + unScalar, + identityMatrix, + zeroV, + zeroM, + -- * Functions + -- ** Scalar & X + (.*^), + (./^), + (.*#), + (./#), + -- ** AccMatrix & Vector + (#*^), + (^*#), + -- ** AccVector & Vector + (^+^), + (^-^), + (^*^), + -- ** AccMatrix & Matrix + (#+#), + (#-#), + (#*#), + (#**.), + -- ** Utility functions + transpose, + zipWithV, + zipWithM, + ) + where + +import qualified Data.Array.Accelerate as A + +import Data.Proxy (Proxy(..)) +import GHC.TypeLits (KnownNat, natVal) +import Data.Array.Accelerate.TypeLits.Internal +import Data.Array.Accelerate ( (:.)((:.)) + , Exp + , DIM2, DIM3, Z(Z) + , IsFloating, IsNum, Elt + , All(All), Any(Any)) + +identityMatrix :: forall n a. (KnownNat n, IsNum a, Elt a) => AccMatrix n n a +-- | constructor for the nxn dimensional identity matrix, given by +-- +-- > ⎛ 1 0 … 0 0 ⎞ +-- > ⎜ 0 1 … 0 0 ⎟ +-- > ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ +-- > ⎜ 0 0 … 1 0 ⎟ +-- > ⎝ 0 0 … 0 1 ⎠ + +identityMatrix = AccMatrix $ A.use $ A.fromFunction (Z:.n':.n') aux + where aux :: DIM2 -> a + aux (Z:.i:.j) = if i == j then 1 else 0 + n' = fromIntegral $ natVal (Proxy :: Proxy n) + +zeroV :: forall n a. (KnownNat n, IsNum a, Elt a) => AccVector n a +-- | constructor for the n dimensional zero vector, given by +-- +-- > ⎛ 0 ⎞ +-- > ⎜ . ⎟ +-- > ⎜ . ⎟ +-- > ⎜ . ⎟ +-- > ⎜ . ⎟ +-- > ⎜ . ⎟ +-- > ⎝ 0 ⎠ + +zeroV = unsafeMkVector $ replicate n' 0 + where n' = fromIntegral $ natVal (Proxy :: Proxy n) + +zeroM :: forall m n a. (KnownNat m, KnownNat n, IsNum a, Elt a) => AccMatrix m n a +-- | constructor for the mxn dimensional zero matrix, given by +-- +-- > ⎛ 0 0 … 0 0 ⎞ +-- > ⎜ 0 0 … 0 0 ⎟ +-- > ⎜ . . . . ⎟ +-- > ⎜ 0 0 … 0 0 ⎟ +-- > ⎝ 0 0 … 0 0 ⎠ + +zeroM = unsafeMkMatrix $ replicate (m'*n') 0 + where n' = fromIntegral $ natVal (Proxy :: Proxy n) + m' = fromIntegral $ natVal (Proxy :: Proxy m) + + +(#*^) :: forall m n a. (KnownNat m, KnownNat n, IsNum a, Elt a) + => AccMatrix m n a -> AccVector n a -> AccVector n a +-- | the usual matrix-vector product +-- +-- > ⎛ w₁₁ w₁₂ … w₁ₙ ⎞ ⎛x₁⎞ ⎛ w₁₁*x₁ + w₁₂*x₂ + … w₁ₙ*xₙ ⎞ +-- > ⎜ w₂₁ w₂₂ … w₂ₙ ⎟ ⎜x₂⎟ ⎜ w₂₁*x₁ + w₂₂*x₂ + … w₂ₙ*xₙ ⎟ +-- > ⎜ . . . ⎟ ⎜. ⎟ ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ ✕ ⎜. ⎟ = ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ ⎜. ⎟ ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ ⎜. ⎟ ⎜ . . . ⎟ +-- > ⎝ wₘ₁ wₘ₂ … wₘₙ ⎠ ⎝xₙ⎠ ⎝ wₘ₁*x₁ + wₘ₂*x₂ + … wₘₙ*xₙ ⎠ + +ma #*^ va = let ma' = unMatrix ma + va' = unVector va + in AccVector $ A.fold1 (+) + $ A.zipWith (*) + ma' + (A.replicate (A.lift $ Z :. m' :. All) va') + where m' = fromIntegral $ natVal (Proxy :: Proxy m) :: Int + +infixl 7 #*^ + +(^*#) :: forall m n a. (KnownNat m, KnownNat n, IsNum a, Elt a) + => AccVector m a -> AccMatrix m n a -> AccVector n a +-- | the usual vector-matrix product +-- +-- > ⎛x₁⎞T ⎛w₁₁ w₁₂ … w₁ₙ ⎞ ⎛ x₁*w₁₁ + x₂*w₁₂ + … xₙ*w₁ₙ ⎞ +-- > ⎜x₂⎟ ⎜w₂₁ w₂₂ … w₂ₙ ⎟ ⎜ x₁*w₂₁ + x₂*w₂₂ + … xₙ*w₂ₙ ⎟ +-- > ⎜. ⎟ ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎜. ⎟ ✕ ⎜ . . . ⎟ = ⎜ . . . ⎟ +-- > ⎜. ⎟ ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎜. ⎟ ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎝xₘ⎠ ⎝wₘ₁ wₘ₂ … wₘₙ ⎠ ⎝ x₁*wₘ₁ + x₂*wₘ₂ + … xₙ*wₘₙ ⎠ + +va ^*# ma = let va' = unVector va + ma' = unMatrix ma + in AccVector $ A.fold1 (+) + $ A.zipWith (*) + (A.replicate (A.lift $ Z :. n' :. All) va') + ma' + where n' = fromIntegral $ natVal (Proxy :: Proxy n) :: Int + +infixr 7 ^*# + +(^+^) :: forall n a. (KnownNat n, IsNum a, Elt a) + => AccVector n a -> AccVector n a -> AccVector n a +-- | the usual vector addition +-- +-- > ⎛v₁⎞ ⎛w₁⎞ ⎛ v₁+w₁ ⎞ +-- > ⎜v₂⎟ ⎜w₂⎟ ⎜ v₂+w₁ ⎟ +-- > ⎜. ⎟ ⎜. ⎟ ⎜ . ⎟ +-- > ⎜. ⎟ + ⎜. ⎟ = ⎜ . ⎟ +-- > ⎜. ⎟ ⎜. ⎟ ⎜ . ⎟ +-- > ⎜. ⎟ ⎜. ⎟ ⎜ . ⎟ +-- > ⎝vₙ⎠ ⎝wₙ⎠ ⎝ vₙ+wₙ ⎠ + +v ^+^ w = AccVector $ A.zipWith (+) (unVector v) (unVector w) +-- | the usual vector subtraction +-- +-- > ⎛v₁⎞ ⎛w₁⎞ ⎛ v₁-w₁ ⎞ +-- > ⎜v₂⎟ ⎜w₂⎟ ⎜ v₂-w₁ ⎟ +-- > ⎜. ⎟ ⎜. ⎟ ⎜ . ⎟ +-- > ⎜. ⎟ - ⎜. ⎟ = ⎜ . ⎟ +-- > ⎜. ⎟ ⎜. ⎟ ⎜ . ⎟ +-- > ⎜. ⎟ ⎜. ⎟ ⎜ . ⎟ +-- > ⎝vₙ⎠ ⎝wₙ⎠ ⎝ vₙ-wₙ ⎠ + +(^-^) :: forall n a. (KnownNat n, IsNum a, Elt a) + => AccVector n a -> AccVector n a -> AccVector n a +v ^-^ w = AccVector $ A.zipWith (-) (unVector v) (unVector w) + +infixl 6 ^+^ +infixl 6 ^-^ + +(^*^) :: forall n a. (KnownNat n, IsNum a, Elt a) + => AccVector n a -> AccVector n a -> AccScalar a +-- | the usual inner product of two vectors +-- +-- > ⎛v₁⎞ ⎛w₁⎞ +-- > ⎜v₂⎟ ⎜w₂⎟ +-- > ⎜. ⎟ ⎜. ⎟ +-- > ⎜. ⎟ * ⎜. ⎟ = v₁*w₁ + v₂*w₁ + … + vₙ*wₙ +-- > ⎜. ⎟ ⎜. ⎟ +-- > ⎜. ⎟ ⎜. ⎟ +-- > ⎝vₙ⎠ ⎝wₙ⎠ + +v ^*^ w = AccScalar $ A.sum $ A.zipWith (*) (unVector v) (unVector w) + +infixl 7 ^*^ + +(#+#) :: forall m n a. (KnownNat m, KnownNat n, IsNum a, Elt a) + => AccMatrix m n a -> AccMatrix m n a -> AccMatrix m n a +-- | the usual matrix addition/subtraction +-- +-- > ⎛ v₁₁ v₁₂ … v₁ₙ ⎞ ⎛ w₁₁ w₁₂ … w₁ₙ ⎞ ⎛ v₁₁+w₁₁ v₁₂+w₁₂ … v₁ₙ+w₁ₙ ⎞ +-- > ⎜ v₂₁ v₂₂ … v₂ₙ ⎟ ⎜ w₂₁ w₂₂ … w₂ₙ ⎟ ⎜ v₂₁+w₂₁ v₂₂+w₂₂ … v₂ₙ+w₂ₙ ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ + ⎜ . . . ⎟ = ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎝ vₘ₁ vₘ₂ … vₘₙ ⎠ ⎝ wₘ₁ wₘ₂ … wₘₙ ⎠ ⎝ vₘ₁+wₘ₁ wₘ₂+vₘ₂ … vₘₙ+wₘₙ ⎠ + +v #+# w = AccMatrix $ A.zipWith (+) (unMatrix v) (unMatrix w) + +(#-#) :: forall m n a. (KnownNat m, KnownNat n, IsNum a, Elt a) + => AccMatrix m n a -> AccMatrix m n a -> AccMatrix m n a +-- | the usual matrix addition/subtraction +-- +-- > ⎛ v₁₁ v₁₂ … v₁ₙ ⎞ ⎛ w₁₁ w₁₂ … w₁ₙ ⎞ ⎛ v₁₁+w₁₁ v₁₂+w₁₂ … v₁ₙ+w₁ₙ ⎞ +-- > ⎜ v₂₁ v₂₂ … v₂ₙ ⎟ ⎜ w₂₁ w₂₂ … w₂ₙ ⎟ ⎜ v₂₁+w₂₁ v₂₂+w₂₂ … v₂ₙ+w₂ₙ ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ + ⎜ . . . ⎟ = ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎝ vₘ₁ vₘ₂ … vₘₙ ⎠ ⎝ wₘ₁ wₘ₂ … wₘₙ ⎠ ⎝ vₘ₁+wₘ₁ wₘ₂+vₘ₂ … vₘₙ+wₘₙ ⎠ + +v #-# w = AccMatrix $ A.zipWith (-) (unMatrix v) (unMatrix w) + +infixl 6 #+# +infixl 6 #-# + +(#*#) :: forall k m n a. (KnownNat k, KnownNat m, KnownNat n, IsNum a, Elt a) + => AccMatrix k m a -> AccMatrix m n a -> AccMatrix k n a +-- | the usual matrix multiplication +-- +-- > ⎛ v₁₁ v₁₂ … v₁ₘ ⎞ ⎛ w₁₁ w₁₂ … w₁ₙ ⎞ ⎛ (v₁₁*w₁₁+v₁₂*w₂₁+…+v₁ₘ*wₘ₁) . . . (v₁₁*w₁ₙ+v₁₂*w₂ₙ+…+v₁ₘ*wₘₙ) ⎞ +-- > ⎜ v₂₁ v₂₂ … v₂ₘ ⎟ ⎜ w₂₁ w₂₂ … w₂ₙ ⎟ ⎜ . . ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ ⎜ . . ⎟ +-- > ⎜ . . . ⎟ * ⎜ . . . ⎟ = ⎜ . . ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ ⎜ . . ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ ⎜ . . ⎟ +-- > ⎝ vₖ₁ vₖ₂ … vₖₘ ⎠ ⎝ wₘ₁ wₘ₂ … wₘₙ ⎠ ⎝ (vₖ₁*w₁₁+vₖ₂*w₂₁+…+vₖₘ*wₘ₁) . . . (vₖ₁*w₁ₙ+vₖ₂*w₂ₙ+…+vₖₘ*wₘₙ) ⎠ + +v #*# w = AccMatrix $ A.fold1 (+) + $ A.backpermute (A.lift $ Z:.ek:.en:.em ) reindex + $ A.zipWith (*) v' w' + where [k',m',n'] = map fromIntegral [ natVal (Proxy :: Proxy k) + , natVal (Proxy :: Proxy m) + , natVal (Proxy :: Proxy n)] :: [Int] + [ek,em,en] = map fromIntegral [k',m',n'] :: [Exp Int] + v' = A.replicate (A.lift $ Any:.All:.All:.k') (unMatrix v) + w' = A.replicate (A.lift $ Any:.n':.All:.All) (unMatrix w) + reindex :: Exp DIM3 -> Exp DIM3 + reindex ix = let (Z:.i:.t:.j) = A.unlift ix + in A.lift (Z:.i:.j:.t :: Z :. Exp Int :. Exp Int :. Exp Int) + +infixl 7 #*# + +(.*^) :: forall n a. (KnownNat n, IsNum a, Elt a) + => Exp a -> AccVector n a -> AccVector n a +-- | the usual multiplication of a scalar with a vector +-- +-- > ⎛x₁⎞ ⎛ a*x₁ ⎞ +-- > ⎜x₂⎟ ⎜ a*x₂ ⎟ +-- > ⎜. ⎟ ⎜ . ⎟ +-- > a • ⎜. ⎟ = ⎜ . ⎟ +-- > ⎜. ⎟ ⎜ . ⎟ +-- > ⎜. ⎟ ⎜ . ⎟ +-- > ⎝xₙ⎠ ⎝ a*xₙ ⎠ + +a .*^ v = let v' = unVector v + in AccVector $ A.map (* a) v' + +(./^) :: forall n a. (KnownNat n, IsFloating a, Elt a) + => Exp a -> AccVector n a -> AccVector n a +-- | a convenient helper deviding every element of a vector +-- +-- > ⎛x₁⎞ ⎛ x₁/a ⎞ +-- > ⎜x₂⎟ ⎜ x₂/a ⎟ +-- > ⎜. ⎟ ⎜ . ⎟ +-- > a / ⎜. ⎟ = ⎜ . ⎟ +-- > ⎜. ⎟ ⎜ . ⎟ +-- > ⎜. ⎟ ⎜ . ⎟ +-- > ⎝xₙ⎠ ⎝ xₙ/a ⎠ +a ./^ v = let v' = unVector v + in AccVector $ A.map (/ a) v' + +infixl 7 .*^ +infixl 7 ./^ + +(.*#) :: forall m n a. (KnownNat m, KnownNat n, IsNum a, Elt a) + => Exp a -> AccMatrix m n a -> AccMatrix m n a +-- | the usual multiplication of a scalar with a matrix +-- +-- > ⎛ w₁₁ w₁₂ … w₁ₙ ⎞ ⎛ a*w₁₁ a*w₁₂ … a*w₁ₙ ⎞ +-- > ⎜ w₂₁ w₂₂ … w₂ₙ ⎟ ⎜ a*w₂₁ a*w₂₂ … a*w₂ₙ ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > a • ⎜ . . . ⎟ = ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎝ wₘ₁ wₘ₂ … wₘₙ ⎠ ⎝ a*wₘ₁ a*wₘ₂ … a*wₘₙ ⎠ + +a .*# v = let v' = unMatrix v + in AccMatrix $ A.map (* a) v' + +(./#) :: forall m n a. (KnownNat m ,KnownNat n, IsFloating a, Elt a) + => Exp a -> AccMatrix m n a -> AccMatrix m n a +-- | a convenient helper deviding every element of a matrix +-- +-- > ⎛ w₁₁ w₁₂ … w₁ₙ ⎞ ⎛ w₁₁/a w₁₂/a … w₁ₙ/a ⎞ +-- > ⎜ w₂₁ w₂₂ … w₂ₙ ⎟ ⎜ w₂₁/a w₂₂/a … w₂ₙ/a ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > a / ⎜ . . . ⎟ = ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ ⎜ . . . ⎟ +-- > ⎝ wₘ₁ wₘ₂ … wₘₙ ⎠ ⎝ wₘ₁/a wₘ₂/a … wₘₙ/a ⎠ +a ./# v = let v' = unMatrix v + in AccMatrix $ A.map (/ a) v' + +infixl 7 .*# +infixl 7 ./# + +(#**.) :: forall n a. (KnownNat n, IsNum a, Elt a) + => AccMatrix n n a -> Int -> AccMatrix n n a +-- | the exponentiation of a square matrix with an `Int`. Negative exponents +-- raise an error - as inverse matrices are not yet implemented. +-- +-- > ⎛ v₁₁ v₁₂ … v₁ₙ ⎞ k +-- > ⎜ v₂₁ v₂₂ … v₂ₙ ⎟ +-- > ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ +-- > ⎜ . . . ⎟ +-- > ⎝ vₙ₁ vₙ₂ … vₙₙ ⎠ + +_ #**. 0 = identityMatrix +v #**. 1 = v +v #**. i | i < 0 = error $ "no negative exponents allowed in matrix exponetiation," + ++ "inverse function not yet implemented" + | otherwise = (v#**. (i-1)) #*# v + +infixr 8 #**. + +transpose :: forall m n a. (KnownNat m, KnownNat n, Elt a) + => AccMatrix m n a -> AccMatrix n m a +-- | transpose for matrices - note the dimension of the matrix change. +transpose = AccMatrix . A.transpose . unMatrix + + +zipWithM :: forall m n a b c. (KnownNat m, KnownNat n, Elt a, Elt b, Elt c) + => (Exp a -> Exp b -> Exp c) -> AccMatrix m n a -> AccMatrix m n b -> AccMatrix m n c +-- | the pendant of the usual zipWith function for matrices, but can only be +-- used with the same dimensions for both input +zipWithM f ma mb = AccMatrix $ A.zipWith f (unMatrix ma) (unMatrix mb) + +zipWithV :: forall n a b c. (KnownNat n, Elt a, Elt b, Elt c) + => (Exp a -> Exp b -> Exp c) -> AccVector n a -> AccVector n b -> AccVector n c +-- | the pendant of the usual zipWith function for vectors, but can only be +-- used with the same dimensions for both input +zipWithV f ma mb = AccVector $ A.zipWith f (unVector ma) (unVector mb) diff --git a/src/Data/Array/Accelerate/TypeLits/Internal.hs b/src/Data/Array/Accelerate/TypeLits/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/TypeLits/Internal.hs @@ -0,0 +1,124 @@ +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE PartialTypeSignatures #-} +module Data.Array.Accelerate.TypeLits.Internal where + +import GHC.TypeLits ( Nat, KnownNat, natVal) + +import Control.Monad (replicateM) + +import qualified Data.Array.Accelerate as A +import qualified Data.Array.Accelerate.Interpreter as I +import Data.Proxy (Proxy(..)) +import Data.Array.Accelerate ( (:.)((:.)), Array + , Exp + , DIM0, DIM1, DIM2, Z(Z) + , Elt, Acc + ) + +import Test.SmallCheck.Series +import Test.QuickCheck.Arbitrary + +newtype AccScalar a = AccScalar { unScalar :: Acc (Array DIM0 a)} + deriving (Show) + +instance forall a. (Eq a, Elt a) => Eq (AccScalar a) where + s == t = let s' = I.run $ unScalar s + t' = I.run $ unScalar t + in A.toList s' == A.toList t' + +-- | A typesafe way to represent an AccVector and its dimension +newtype AccVector (dim :: Nat) a = AccVector { unVector :: Acc (Array DIM1 a)} + deriving (Show) + +instance forall n a. (KnownNat n, Eq a, Elt a) => Eq (AccVector n a) where + v == w = let v' = I.run $ unVector v + w' = I.run $ unVector w + in A.toList v' == A.toList w' + +instance forall mm n a. (Serial mm a, KnownNat n, Eq a, Elt a) + => Serial mm (AccVector n a) where + series = AccVector . A.use . A.fromList (Z:.n') <$> cons1 (replicate n') + where n' = fromIntegral $ natVal (Proxy :: Proxy n) + +instance forall n a. (KnownNat n, Arbitrary a, Eq a, Elt a) + => Arbitrary (AccVector n a) where + arbitrary = AccVector . A.use . A.fromList (Z:.n') <$> replicateM n' arbitrary + where n' = fromIntegral $ natVal (Proxy :: Proxy n) + +-- | A typesafe way to represent an AccMatrix and its rows/colums +newtype AccMatrix (rows :: Nat) (cols :: Nat) a = AccMatrix {unMatrix :: Acc (Array DIM2 a)} + deriving (Show) + +instance forall m n a. (KnownNat m, KnownNat n, Eq a, Elt a) => Eq (AccMatrix m n a) where + v == w = let v' = I.run $ unMatrix v + w' = I.run $ unMatrix w + in A.toList v' == A.toList w' + +instance forall mm m n a. (Serial mm a, KnownNat m, KnownNat n, Eq a, Elt a) + => Serial mm (AccMatrix m n a) where + series = AccMatrix . A.use . A.fromList (Z:.m':.n') <$> cons1 (replicate $ m'*n') + where m' = fromIntegral $ natVal (Proxy :: Proxy m) + n' = fromIntegral $ natVal (Proxy :: Proxy n) + +instance forall m n a. (KnownNat m, KnownNat n, Arbitrary a, Eq a, Elt a) + => Arbitrary (AccMatrix m n a) where + arbitrary = AccMatrix . A.use . A.fromList (Z:.m':.n') <$> replicateM (m'*n') arbitrary + where m' = fromIntegral $ natVal (Proxy :: Proxy m) + n' = fromIntegral $ natVal (Proxy :: Proxy n) + +-- | a functor like instance for a functor like instance for Accelerate computations +-- instead of working with simple functions `(a -> b)` this uses (Exp a -> Exp b) +class AccFunctor f where + afmap :: forall a b. (Elt a, Elt b) => (Exp a -> Exp b) -> f a -> f b + +instance AccFunctor AccScalar where + afmap f (AccScalar a) = AccScalar (A.map f a) + +instance forall n. (KnownNat n) => AccFunctor (AccVector n) where + afmap f (AccVector a) = AccVector (A.map f a) + +instance forall m n. (KnownNat m, KnownNat n) => AccFunctor (AccMatrix m n) where + afmap f (AccMatrix a) = AccMatrix (A.map f a) + +mkVector :: forall n a. (KnownNat n, Elt a) => [a] -> Maybe (AccVector n a) +-- | a smart constructor to generate Vectors - returning Nothing +-- if the input list is not as long as the dimension of the Vector +mkVector as = if length as == n' + then Just $ unsafeMkVector as + else Nothing + where n' = fromIntegral $ natVal (Proxy :: Proxy n) + +unsafeMkVector :: forall n a. (KnownNat n, Elt a) => [a] -> AccVector n a +-- | unsafe smart constructor to generate Vectors +-- the length of the input list is not checked +unsafeMkVector as = AccVector (A.use $ A.fromList (Z:.n') as) + where n' = fromIntegral $ natVal (Proxy :: Proxy n) + +mkMatrix :: forall m n a. (KnownNat m, KnownNat n, Elt a) + => [a] -> Maybe (AccMatrix m n a) +-- | a smart constructor to generate Matrices - returning Nothing +-- if the input list is not as long as the "length" of the Matrix, i.e. rows*colums +mkMatrix as = if length as == m'*n' + then Just $ unsafeMkMatrix as + else Nothing + where m' = fromIntegral $ natVal (Proxy :: Proxy m) + n' = fromIntegral $ natVal (Proxy :: Proxy n) + +unsafeMkMatrix :: forall m n a. (KnownNat m, KnownNat n, Elt a) + => [a] -> AccMatrix m n a +-- | unsafe smart constructor to generate Matrices +-- the length of the input list is not checked +unsafeMkMatrix as = AccMatrix (A.use $ A.fromList (Z:. m':.n') as) + where m' = fromIntegral $ natVal (Proxy :: Proxy m) + n' = fromIntegral $ natVal (Proxy :: Proxy n) + +mkScalar :: forall a. Elt a => Exp a -> AccScalar a +-- | a smart constructor to generate scalars +mkScalar = AccScalar . A.unit diff --git a/src/Data/Array/Accelerate/TypeLits/System/Random/MWC.hs b/src/Data/Array/Accelerate/TypeLits/System/Random/MWC.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/TypeLits/System/Random/MWC.hs @@ -0,0 +1,36 @@ +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# OPTIONS_GHC -fno-warn-unused-imports #-} +module Data.Array.Accelerate.TypeLits.System.Random.MWC + ( rndMatrixWith + , rndVectorWith + , module Distributions + ) where + +import Data.Array.Accelerate.TypeLits.Internal + +import GHC.TypeLits +import Data.Proxy +import qualified Data.Array.Accelerate as A +import Data.Array.Accelerate (Elt, Z(..), (:.)(..)) +import Data.Array.Accelerate.System.Random.MWC +import System.Random.MWC.Distributions as Distributions + + +rndMatrixWith :: forall m n e. (KnownNat m, KnownNat n, Elt e) => (GenIO -> IO e) -> IO (AccMatrix m n e) +-- | mwc random provides a fast and "statistically-safe" random distribution to +-- work with this +rndMatrixWith cdf = do r <- randomArray (const cdf) sh + return $ AccMatrix $ A.use r + where m' = fromInteger $ natVal (Proxy :: Proxy m) + n' = fromInteger $ natVal (Proxy :: Proxy n) + sh = Z:.m':.n' + +rndVectorWith :: forall n e. (KnownNat n, Elt e) => (GenIO -> IO e) -> IO (AccVector n e) +-- | mwc random provides a fast and "statistically-safe" random distribution to +-- work with this +rndVectorWith cdf = do r <- randomArray (const cdf) sh + return $ AccVector $ A.use r + where n' = fromInteger $ natVal (Proxy :: Proxy n) + sh = Z:.n' diff --git a/stack-7.10.yaml b/stack-7.10.yaml new file mode 100644 --- /dev/null +++ b/stack-7.10.yaml @@ -0,0 +1,12 @@ +resolver: lts-5.12 + +packages: + - '.' + +extra-deps: + - accelerate-random-0.15.0.0 + - HUnit-Plus-1.1.0 + +flags: {} + +extra-package-dbs: [] diff --git a/test/Test.hs b/test/Test.hs new file mode 100644 --- /dev/null +++ b/test/Test.hs @@ -0,0 +1,18 @@ +module Main where + +import Test.Tasty + +import qualified Test.Data.Array.Accelerate.TypeLits +import qualified Test.Data.Array.Accelerate.TypeLits.Internal +import qualified Test.Data.Array.Accelerate.TypeLits.System.Random.MWC + + +main :: IO () +main = defaultMain tests + +tests :: TestTree +tests = testGroup "Tests" + [ Test.Data.Array.Accelerate.TypeLits.tests + , Test.Data.Array.Accelerate.TypeLits.Internal.tests + , Test.Data.Array.Accelerate.TypeLits.System.Random.MWC.tests + ] diff --git a/test/Test/Data/Array/Accelerate/TypeLits.hs b/test/Test/Data/Array/Accelerate/TypeLits.hs new file mode 100644 --- /dev/null +++ b/test/Test/Data/Array/Accelerate/TypeLits.hs @@ -0,0 +1,310 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# OPTIONS_GHC -fno-warn-unused-imports #-} +module Test.Data.Array.Accelerate.TypeLits where + +import Test.Tasty +import Test.Tasty.SmallCheck as SC +import Test.Tasty.QuickCheck as QC +import Test.Tasty.HUnit +import qualified Test.HUnitPlus.Base as H +import qualified Test.SmallCheck.Series as SC + +import Data.Array.Accelerate.TypeLits + +import Data.Maybe (isJust, isNothing) +import Control.Exception + +tests :: TestTree +tests = testGroup "Data.Array.Accelerate.TypeLits" [properties, unitTests] + +properties :: TestTree +properties = testGroup "Properties" [scProps, qcProps] + +scProps :: TestTree +scProps = testGroup "(checked by SmallCheck)" + [ testGroup "* Functions" + [ testGroup "** Classes" + [ testGroup "AccFunctor" + [ SC.testProperty "Vector: afmap id = id" $ + \v -> afmap id (v :: AccVector 3 Int) == v + , SC.testProperty "Matrix: afmap id = id" $ + \ma -> afmap id (ma :: AccMatrix 3 3 Int) == ma + --scprop + ] + ] + , testGroup "** Scalar & X" + [ testGroup ".*^" + [ SC.testProperty "universal property of zero" $ + \v -> 0 .*^ v == (zeroV :: AccVector 2 Int) + , SC.testProperty "universal property of 1" $ + \v -> 1 .*^ v == (v :: AccVector 2 Int) + , SC.testProperty "universal property of 2" $ + \v -> 2 .*^ v == v ^+^ (v :: AccVector 2 Int) + --scprop + ] + , testGroup "./^" + [ SC.testProperty "universal property of 1" $ + \v -> 1 ./^ v == (v :: AccVector 2 Float) + --scprop + ] + , testGroup ".*#" + [ SC.testProperty "universal property of zero" $ + \ma -> 0 .*# ma == (zeroM :: AccMatrix 3 2 Int) + , SC.testProperty "universal property of 1" $ + \ma -> 1 .*# ma == (ma :: AccMatrix 3 2 Int) + , SC.testProperty "universal property of 2" $ + \ma -> 2 .*# ma == ma #+# (ma :: AccMatrix 3 2 Int) + --scprop + ] + , testGroup "./#" + [ SC.testProperty "universal property of 1" $ + \ma -> 1 ./# ma == (ma :: AccMatrix 3 2 Float) + --scprop + ] + ] + , testGroup "** AccMatrix & AccVector" + [ testGroup "#*^" + [ SC.testProperty "id . v == v" $ + \v -> identityMatrix #*^ v == (v :: AccVector 2 Int) + , SC.testProperty "0 . v == 0" $ + \v -> (zeroM :: AccMatrix 2 2 Int) #*^ v == zeroV + , SC.testProperty "ma . 0 == 0" $ + \ma -> (ma :: AccMatrix 2 2 Int) #*^ zeroV == zeroV + --scprop + ] + , testGroup "^*#" + [ SC.testProperty "v . id == v" $ + \v -> v ^*# identityMatrix == (v :: AccVector 2 Int) + , SC.testProperty "v . 0 == 0" $ + \v -> v ^*# (zeroM :: AccMatrix 2 2 Int) == zeroV + , SC.testProperty "ma . 0 == 0" $ + \ma -> zeroV ^*# (ma :: AccMatrix 2 2 Int) == zeroV + --scprop + ] + ] + , testGroup "** AccVector & AccVector" + [ testGroup "^+^" + [ SC.testProperty "0 + v = v" $ + \v -> zeroV ^+^ v == (v :: AccVector 2 Int) + , SC.testProperty "v + 0 = v" $ + \v -> zeroV ^+^ v == (v :: AccVector 2 Int) + , SC.testProperty "v + (-v) = 0 with operator precedence" $ + \v -> v ^+^ (-1) .*^ v == (zeroV :: AccVector 2 Int) + , SC.testProperty "(-v) + v = 0 with operator precedence" $ + \v -> (-1) .*^ v ^+^ v == (zeroV :: AccVector 2 Int) + , SC.testProperty "commutativity: v + w = w + v" $ + \v w -> v ^+^ w == w ^+^ (v :: AccVector 2 Int) + , SC.testProperty "associativity: (u + v) + w = u + (v + w)" $ + \u v w -> (u ^+^ v) ^+^ w == (u ^+^ (v ^+^ w) :: AccVector 2 Int) + --scprop + ] + , testGroup "^-^" + [ SC.testProperty "v - v = 0" $ + \v -> v ^-^ v == (zeroV :: AccVector 2 Int) + ] + , testGroup "^*^" + [ SC.testProperty "zeroV is orthogonal for everything" $ + \v -> zeroV ^*^ (v :: AccVector 2 Int) == mkScalar 0 + , SC.testProperty "zeroV is orthogonal for everything" $ + \v -> (v :: AccVector 2 Int) ^*^ zeroV == mkScalar 0 + , SC.testProperty "inner product is commutative" $ + \v w -> v ^*^ w == w ^*^ (v :: AccVector 2 Int) + --scprop + ] + ] + , testGroup "** AccMatrix & AccMatrix" + [ testGroup "#+#" + [ SC.testProperty "0 + ma = ma" $ + \ma -> zeroM #+# ma == (ma :: AccMatrix 3 2 Int) + , SC.testProperty "ma + 0 = ma" $ + \ma -> zeroM #+# ma == (ma :: AccMatrix 3 2 Int) + , SC.testProperty "ma + (-ma) = 0 with operator precedence" $ + \ma -> ma #+# (-1) .*# ma == (zeroM :: AccMatrix 3 2 Int) + , SC.testProperty "(-ma) + ma = 0 with operator precedence" $ + \ma -> (-1) .*# ma #+# ma == (zeroM :: AccMatrix 3 2 Int) + , SC.testProperty "commutativity: ma + -mb = -mb + ma" $ + \ma mb -> ma #+# mb == mb #+# (ma :: AccMatrix 3 2 Int) + , SC.testProperty "associativity: (mc + ma) + -mb = mc + (ma + -mb)" $ + \ma mb mc -> (ma #+# mb) #+# mc == (ma #+# (mb #+# mc) :: AccMatrix 3 2 Int) + --scprop + ] + , testGroup "#-#" + [ SC.testProperty "ma - ma = 0" $ + \ma -> ma #-# ma == (zeroM :: AccMatrix 3 2 Int) + ] + , testGroup "#*#" + [ SC.testProperty "id * ma = ma" $ + \ma -> identityMatrix #*# ma == (ma :: AccMatrix 3 2 Int) + , SC.testProperty "ma * id = ma" $ + \ma -> ma #*# identityMatrix == (ma :: AccMatrix 3 2 Int) + --scprop + ] + , testGroup "#**." + [ SC.testProperty "id^n = id" $ + \n -> let n' = SC.getNonNegative n + in identityMatrix #**. n' == (identityMatrix :: AccMatrix 2 2 Int) + , SC.testProperty "exponential law: a^n * a^m = a^(n+m)" $ + \n m ma -> let m' = SC.getNonNegative m + n' = SC.getNonNegative n + in (ma #**. m') #*# (ma #**. n') == (ma #**. (m' +n') :: AccMatrix 2 2 Int) + , SC.testProperty "exponentiation of diagonal matrices" $ + \k l m n -> let n' = SC.getNonNegative n + in unsafeMkMatrix [k,0,0 + ,0,l,0 + ,0,0,m] #**. n' + == (unsafeMkMatrix [k^n', 0 , 0 + , 0 ,l^n', 0 + , 0 , 0 ,m^n'] :: AccMatrix 3 3 Int) + --scprop + ] + ] + , testGroup "** Utility functions" + [ testGroup "transpose" + [ SC.testProperty "transpose . transpose = id" $ + \ma -> transpose (transpose ma) == (ma :: AccMatrix 3 2 Int) + --scprop + ] + , testGroup "zipWithV" + [ SC.testProperty "zipWithV const v w = v" $ + \v w -> zipWithV const (v :: AccVector 3 Int) (w :: AccVector 3 Int) == v + --scprop + ] + , testGroup "zipWithM" + [ SC.testProperty "zipWithM const ma mb = ma" $ + \ma mb -> zipWithM const (ma :: AccMatrix 3 2 Int) (mb :: AccMatrix 3 2 Int) == ma + --scprop + ] + ] + ] + ] + +qcProps :: TestTree +qcProps = testGroup "(checked by QuickCheck)" + [--qcprop + ] + +unitTests :: TestTree +unitTests = testGroup "HUnit tests" + [ testGroup "* Constructors" + [ testGroup "identity" + [ testCase "dim 2" $ + unsafeMkMatrix [1,0 + ,0,1] + @=? (identityMatrix :: AccMatrix 2 2 Int) + , testCase "dim 3" $ + unsafeMkMatrix [1,0,0 + ,0,1,0 + ,0,0,1] + @=? (identityMatrix :: AccMatrix 3 3 Int) + , testCase "dim 4" $ + unsafeMkMatrix [1,0,0,0 + ,0,1,0,0 + ,0,0,1,0 + ,0,0,0,1] + @=? (identityMatrix :: AccMatrix 4 4 Int) + --hunit + ] + , testGroup "mkMatrix" + [ testCase "dim4 - ok" $ + isJust (mkMatrix [1..16] :: Maybe (AccMatrix 4 4 Int)) @=? True + , testCase "dim4 - fail - too short" $ + isNothing (mkMatrix [1..15] :: Maybe (AccMatrix 4 4 Int)) @=? True + , testCase "dim4 - fail - too long" $ + isNothing (mkMatrix [1..17] :: Maybe (AccMatrix 4 4 Int)) @=? True + --hunit + ] + , testGroup "mkVector" + [ testCase "dim4 - ok" $ + isJust (mkVector [1..4] :: Maybe (AccVector 4 Int)) @=? True + , testCase "dim4 - fail - too short" $ + isNothing (mkVector [1..3] :: Maybe (AccVector 4 Int)) @=? True + , testCase "dim4 - fail - too long" $ + isNothing (mkVector [1..5] :: Maybe (AccVector 4 Int)) @=? True + --hunit + ] + ] + , testGroup "* Functions" + [ testGroup "** Scalar & X" + [ testGroup ".*^" + [ testCase "2 .*^ <1,2>" $ + 2 .*^ (unsafeMkVector [1,2] :: AccVector 2 Int) @?= (unsafeMkVector [2,4] :: AccVector 2 Int) + , testCase "2 .*^ <1,2>" $ + (-2) .*^ (unsafeMkVector [1,2] :: AccVector 2 Int) @?= (unsafeMkVector [-2,-4] :: AccVector 2 Int) + --hunit + ] + , testGroup "./^" + [--hunit + ] + , testGroup ".*#" + [--hunit + ] + , testGroup "./#" + [--hunit + ] + ] + , testGroup "** AccMatrix & AccVector" + [ testGroup "#*^" + [ testCase "m #*^ v" $ + fmap (i4 #*^) v4 @?= v4 + --hunit + ] + , testGroup "^*#" + [--hunit + ] + ] + , testGroup "** AccVector & AccVector" + [ testGroup "^+^" + [--hunit + ] + , testGroup "^-^" + [--hunit + ] + , testGroup "^*^" + [--hunit + ] + ] + , testGroup "** AccMatrix & AccMatrix" + [ testGroup "#+#" + [--hunit + ] + , testGroup "#-#" + [--hunit + ] + , testGroup "#*#" + [ testCase "i4 #*# i4 == i4 " $ + i4 #*# i4 @?= i4 + , testCase "a #*# i4 == i4 " $ + i4 #*# i4 @?= i4 + --hunit + ] + , testGroup "#**." + [testCase "id ^ (-1) = fail" $ + H.assertThrowsExact (ErrorCall "") + (return (identityMatrix #**. (-1)) :: IO (AccMatrix 3 3 Int)) + ] + ] + , testGroup "** Utility functions" + [ testGroup "transpose" + [--hunit + ] + , testGroup "zipWithV" + [--hunit + ] + , testGroup "zipWithM" + [--hunit + ] + ] + ] + ] +i4 :: AccMatrix 4 4 Int +i4 = identityMatrix + +j4 :: Maybe (AccMatrix 4 4 Int) +j4 = Just identityMatrix + +v4 :: Maybe (AccVector 4 Int) +v4 = mkVector [1,1,1,1] + +m4 :: Maybe (AccMatrix 4 4 Int) +m4 = mkMatrix [1,0,0,0 ,0,1,0,0 ,0,0,1,0 ,0,0,0,1]