diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,3 +1,8 @@ +# 0.4.0.0, 2018-11-24 +- added multi-core parallelism, sample benchmark is available in ./benchmark/results.zip +- more robust QuickCheck tests +- some API changes; simplified Multilinear class + # 0.3.2.0, 2018-11-18 - added filter and zipWith functions diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,4 +1,5 @@ -# Multilinear [![Hackage](https://img.shields.io/hackage/v/multilinear.svg)](https://hackage.haskell.org/package/multilinear) [![Build Status](https://travis-ci.org/ArturB/multilinear.svg?branch=master)](https://travis-ci.org/ArturB/multilinear) +# Multilinear [![Hackage](https://img.shields.io/hackage/v/multilinear.svg)](https://hackage.haskell.org/package/multilinear) [![Build Status](https://travis-ci.org/ArturB/multilinear.svg?branch=master)](https://travis-ci.org/ArturB/multilinear) ![BuildStatus](https://ci.appveyor.com/api/projects/status/github/ArturB/multilinear) + ## Summary Multilinear is general - purpose linear algebra and multi-dimensional array library for Haskell. It provides generic and efficient implementation of linear algebra operations on vectors, linear functionals, matrices and its higher - rank analoges: tensors. It can also be used as simply a miltidimensional arrays for everyone. diff --git a/Setup.hs b/Setup.hs deleted file mode 100644 --- a/Setup.hs +++ /dev/null @@ -1,2 +0,0 @@ -import Distribution.Simple -main = defaultMain diff --git a/benchmark/memory/Bench.hs b/benchmark/memory/Bench.hs deleted file mode 100644 --- a/benchmark/memory/Bench.hs +++ /dev/null @@ -1,71 +0,0 @@ -{-| -Module : Bench -Description : Benchmark of Multilinear library -Copyright : (c) Artur M. Brodzki, 2018 -License : BSD3 -Maintainer : artur@brodzki.org -Stability : experimental -Portability : Windows/POSIX - --} - -module Main ( - main -) where - -import Weigh -import qualified Multilinear.Matrix as Matrix -import qualified Multilinear.Vector as Vector - --- | Simple generator function for benchmarked matrices -gen :: Int -> Int -> Double -gen j k = sin (fromIntegral j) + cos (fromIntegral k) - --- matrix sizes -s1 :: Int -s1 = 64 -s2 :: Int -s2 = 256 -s3 :: Int -s3 = 1024 - --- | ENTRY POINT -main :: IO () -main = mainWith (do - setColumns [Case, Allocated, GCs, Live, Max] - - -- Benchmarking small vectors - func "vector 1 elem generation" (Vector.fromIndices "i" 1) id - func "vector 2 elem generation" (Vector.fromIndices "i" 2) id - func "vector 3 elem generation" (Vector.fromIndices "i" 3) id - - -- Benchmarking matrix generators - func "matrix 64 x 64 generation" - (Matrix.fromIndices "ij" s1 s1) gen - func "matrix 256 x 256 generation" - (Matrix.fromIndices "ij" s2 s2) gen - func "matrix 1024 x 1024 generation" - (Matrix.fromIndices "ij" s3 s3) gen - - -- Benchmarking matrix addition - func "matrix 64 x 64 addition" - (+ Matrix.fromIndices "ab" s1 s1 gen) - (Matrix.fromIndices "ab" s1 s1 (\a b -> fromIntegral a + fromIntegral b)) - func "matrix 256 x 256 addition" - (+ Matrix.fromIndices "ab" s2 s2 gen) - (Matrix.fromIndices "ab" s2 s2 (\a b -> fromIntegral a + fromIntegral b)) - func "matrix 1024 x 1024 addition" - (+ Matrix.fromIndices "ab" s3 s3 gen) - (Matrix.fromIndices "ab" s3 s3 (\a b -> fromIntegral a + fromIntegral b)) - - -- Benchmarking matrix multiplication - func "matrix 40 x 4,000 multiplication" - (* Matrix.fromIndices "jk" 4000 40 gen) - (Matrix.fromIndices "ij" 40 4000 gen) - func "matrix 40 x 16,000 multiplication" - (* Matrix.fromIndices "jk" 16000 40 gen) - (Matrix.fromIndices "ij" 40 16000 gen) - func "matrix 40 x 64,000 multiplication" - (* Matrix.fromIndices "jk" 64000 40 gen) - (Matrix.fromIndices "ij" 40 64000 gen) - ) diff --git a/benchmark/multicore/memory/Bench.hs b/benchmark/multicore/memory/Bench.hs new file mode 100644 --- /dev/null +++ b/benchmark/multicore/memory/Bench.hs @@ -0,0 +1,73 @@ +{-| +Module : Bench +Description : Benchmark of Multilinear library +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Main ( + main +) where + +import Weigh +import Multilinear.Class +import Multilinear.Generic.MultiCore +import qualified Multilinear.Matrix as Matrix +import qualified Multilinear.Vector as Vector + +-- | Simple generator function for benchmarked matrices +gen :: Int -> Int -> Double +gen j k = sin (fromIntegral j) + cos (fromIntegral k) + +-- matrix sizes +s1 :: Int +s1 = 64 +s2 :: Int +s2 = 256 +s3 :: Int +s3 = 1024 + +-- | ENTRY POINT +main :: IO () +main = mainWith (do + setColumns [Case, Allocated, GCs, Live, Max] + + -- Benchmarking small vectors + --value "vector 1 elem generation" (Vector.fromIndices "i" 1 fromIntegral :: Tensor Double) + --value "vector 2 elem generation" (Vector.fromIndices "i" 2 fromIntegral :: Tensor Double) + --value "vector 3 elem generation" (Vector.fromIndices "i" 3 fromIntegral :: Tensor Double) + + -- Benchmarking matrix generators + --value "matrix 64 x 64 generation" + -- (Matrix.fromIndices "ij" s1 s1 gen :: Tensor Double) + --value "matrix 256 x 256 generation" + -- (Matrix.fromIndices "ij" s2 s2 gen :: Tensor Double) + --value "matrix 1024 x 1024 generation" + -- (Matrix.fromIndices "ij" s3 s3 gen :: Tensor Double) + + -- Benchmarking matrix addition + value "matrix 64 x 64 addition" $ + (+ Matrix.fromIndices "ab" s1 s1 gen) + (Matrix.fromIndices "ab" s1 s1 (\a b -> fromIntegral a + fromIntegral b) :: Tensor Double) + --value "matrix 256 x 256 addition" $ + -- (+ Matrix.fromIndices "ab" s2 s2 gen) + -- (Matrix.fromIndices "ab" s2 s2 (\a b -> fromIntegral a + fromIntegral b) :: Tensor Double) + --value "matrix 1024 x 1024 addition" $ + -- (+ Matrix.fromIndices "ab" s3 s3 gen) + -- (Matrix.fromIndices "ab" s3 s3 (\a b -> fromIntegral a + fromIntegral b) :: Tensor Double) + + -- Benchmarking matrix multiplication + --value "matrix 40 x 4,000 multiplication" $ + -- (* Matrix.fromIndices "jk" 4000 40 gen) + -- (Matrix.fromIndices "ij" 40 4000 gen :: Tensor Double) + --value "matrix 40 x 16,000 multiplication" $ + -- (* Matrix.fromIndices "jk" 16000 40 gen) + -- (Matrix.fromIndices "ij" 40 16000 gen :: Tensor Double) + --value "matrix 40 x 64,000 multiplication" $ + -- (* Matrix.fromIndices "jk" 64000 40 gen) + -- (Matrix.fromIndices "ij" 40 64000 gen :: Tensor Double) + ) diff --git a/benchmark/multicore/profile/Bench.hs b/benchmark/multicore/profile/Bench.hs new file mode 100644 --- /dev/null +++ b/benchmark/multicore/profile/Bench.hs @@ -0,0 +1,27 @@ +{-| +Module : Bench +Description : Benchmark of Multilinear library +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Main ( + main +) where + +import Control.DeepSeq +import Multilinear.Class +import Multilinear.Generic.MultiCore +import qualified Multilinear.Matrix as Matrix + +gen :: Int -> Int -> Double +gen j k = sin (fromIntegral j) + cos (fromIntegral k) + +main :: IO () +main = do + let m = (Matrix.fromIndices "ij" 1000 1000 gen :: Tensor Double) * (Matrix.fromIndices "jk" 1000 1000 gen :: Tensor Double) + m `deepseq` putStrLn $ "All done! Indices of m:" ++ show (indices m) diff --git a/benchmark/multicore/time/Bench.hs b/benchmark/multicore/time/Bench.hs new file mode 100644 --- /dev/null +++ b/benchmark/multicore/time/Bench.hs @@ -0,0 +1,46 @@ +{-| +Module : Bench +Description : Benchmark of Multilinear library +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Main ( + main +) where + +import Criterion.Main +import Multilinear.Class +import Multilinear.Generic.MultiCore +import qualified Multilinear.Matrix as Matrix + +-- | Simple generator function for bencharking matrices +gen :: Int -> Int -> Double +gen j k = sin (fromIntegral j) + cos (fromIntegral k) + +-- | Generate benchmark of matrix multiplication +sizedMatrixMultBench :: + Int -- ^ size of square matrix to multiplicate + -> Benchmark +sizedMatrixMultBench s = + bench (show s ++ "x" ++ show s) $ + nf ((Matrix.fromIndices "ij" s s gen :: Tensor Double) *) (Matrix.fromIndices "jk" s s gen :: Tensor Double) + +-- | Generate benchmark of matrix addition +sizedMatrixAddBench :: + Int -- ^ size of square matrix to add + -> Benchmark +sizedMatrixAddBench s = + bench (show s ++ "x" ++ show s) $ + nf ((Matrix.fromIndices "ij" s s gen :: Tensor Double) +) (Matrix.fromIndices "ij" s s gen :: Tensor Double) + +-- | ENTRY POINT +main :: IO () +main = defaultMain [ + bgroup "matrix addition" $ sizedMatrixAddBench <$> [64, 128, 256, 512, 1024], + bgroup "matrix multiplication" $ sizedMatrixMultBench <$> [64, 128, 256, 512, 1024] + ] diff --git a/benchmark/profile/Bench.hs b/benchmark/profile/Bench.hs deleted file mode 100644 --- a/benchmark/profile/Bench.hs +++ /dev/null @@ -1,27 +0,0 @@ -{-| -Module : Bench -Description : Benchmark of Multilinear library -Copyright : (c) Artur M. Brodzki, 2018 -License : BSD3 -Maintainer : artur@brodzki.org -Stability : experimental -Portability : Windows/POSIX - --} - -module Main ( - main -) where - -import Control.DeepSeq -import Multilinear.Class -import qualified Multilinear.Matrix as Matrix -import qualified Multilinear.Vector as Vector - -gen :: Int -> Int -> Double -gen j k = sin (fromIntegral j) + cos (fromIntegral k) - -main :: IO () -main = do - let m = (Matrix.fromIndices "ij" 40 6000 gen) * (Matrix.fromIndices "jk" 6000 40 gen) - m `deepseq` putStrLn $ "All done! Indices of m:" ++ show (indices m) diff --git a/benchmark/sequential/memory/Bench.hs b/benchmark/sequential/memory/Bench.hs new file mode 100644 --- /dev/null +++ b/benchmark/sequential/memory/Bench.hs @@ -0,0 +1,72 @@ +{-| +Module : Bench +Description : Benchmark of Multilinear library +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Main ( + main +) where + +import Weigh +import Multilinear +import qualified Multilinear.Matrix as Matrix +import qualified Multilinear.Vector as Vector + +-- | Simple generator function for benchmarked matrices +gen :: Int -> Int -> Double +gen j k = sin (fromIntegral j) + cos (fromIntegral k) + +-- matrix sizes +s1 :: Int +s1 = 64 +s2 :: Int +s2 = 256 +s3 :: Int +s3 = 1024 + +-- | ENTRY POINT +main :: IO () +main = mainWith (do + setColumns [Case, Allocated, GCs, Live, Max] + + -- Benchmarking small vectors + value "vector 1 elem generation" (Vector.fromIndices "i" 1 fromIntegral :: Tensor Double) + value "vector 2 elem generation" (Vector.fromIndices "i" 2 fromIntegral :: Tensor Double) + value "vector 3 elem generation" (Vector.fromIndices "i" 3 fromIntegral :: Tensor Double) + + -- Benchmarking matrix generators + value "matrix 64 x 64 generation" + (Matrix.fromIndices "ij" s1 s1 gen :: Tensor Double) + value "matrix 256 x 256 generation" + (Matrix.fromIndices "ij" s2 s2 gen :: Tensor Double) + value "matrix 1024 x 1024 generation" + (Matrix.fromIndices "ij" s3 s3 gen :: Tensor Double) + + -- Benchmarking matrix addition + value "matrix 64 x 64 addition" $ + (+ Matrix.fromIndices "ab" s1 s1 gen) + (Matrix.fromIndices "ab" s1 s1 (\a b -> fromIntegral a + fromIntegral b) :: Tensor Double) + value "matrix 256 x 256 addition" $ + (+ Matrix.fromIndices "ab" s2 s2 gen) + (Matrix.fromIndices "ab" s2 s2 (\a b -> fromIntegral a + fromIntegral b) :: Tensor Double) + value "matrix 1024 x 1024 addition" $ + (+ Matrix.fromIndices "ab" s3 s3 gen) + (Matrix.fromIndices "ab" s3 s3 (\a b -> fromIntegral a + fromIntegral b) :: Tensor Double) + + -- Benchmarking matrix multiplication + value "matrix 40 x 4,000 multiplication" $ + (* Matrix.fromIndices "jk" 4000 40 gen) + (Matrix.fromIndices "ij" 40 4000 gen :: Tensor Double) + value "matrix 40 x 16,000 multiplication" $ + (* Matrix.fromIndices "jk" 16000 40 gen) + (Matrix.fromIndices "ij" 40 16000 gen :: Tensor Double) + value "matrix 40 x 64,000 multiplication" $ + (* Matrix.fromIndices "jk" 64000 40 gen) + (Matrix.fromIndices "ij" 40 64000 gen :: Tensor Double) + ) diff --git a/benchmark/sequential/profile/Bench.hs b/benchmark/sequential/profile/Bench.hs new file mode 100644 --- /dev/null +++ b/benchmark/sequential/profile/Bench.hs @@ -0,0 +1,26 @@ +{-| +Module : Bench +Description : Benchmark of Multilinear library +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Main ( + main +) where + +import Control.DeepSeq +import Multilinear +import qualified Multilinear.Matrix as Matrix + +gen :: Int -> Int -> Double +gen j k = sin (fromIntegral j) + cos (fromIntegral k) + +main :: IO () +main = do + let m = (Matrix.fromIndices "ij" 1000 1000 gen :: Tensor Double) + (Matrix.fromIndices "jk" 1000 1000 gen :: Tensor Double) + m `deepseq` putStrLn $ "All done! Indices of m:" ++ show (indices m) diff --git a/benchmark/sequential/time/Bench.hs b/benchmark/sequential/time/Bench.hs new file mode 100644 --- /dev/null +++ b/benchmark/sequential/time/Bench.hs @@ -0,0 +1,45 @@ +{-| +Module : Bench +Description : Benchmark of Multilinear library +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Main ( + main +) where + +import Criterion.Main +import Multilinear +import qualified Multilinear.Matrix as Matrix + +-- | Simple generator function for bencharking matrices +gen :: Int -> Int -> Double +gen j k = sin (fromIntegral j) + cos (fromIntegral k) + +-- | Generate benchmark of matrix multiplication +sizedMatrixMultBench :: + Int -- ^ size of square matrix to multiplicate + -> Benchmark +sizedMatrixMultBench s = + bench (show s ++ "x" ++ show s) $ + nf ((Matrix.fromIndices "ij" s s gen :: Tensor Double) *) (Matrix.fromIndices "jk" s s gen :: Tensor Double) + +-- | Generate benchmark of matrix addition +sizedMatrixAddBench :: + Int -- ^ size of square matrix to add + -> Benchmark +sizedMatrixAddBench s = + bench (show s ++ "x" ++ show s) $ + nf ((Matrix.fromIndices "ij" s s gen :: Tensor Double) +) (Matrix.fromIndices "ij" s s gen :: Tensor Double) + +-- | ENTRY POINT +main :: IO () +main = defaultMain [ + bgroup "matrix addition" $ sizedMatrixAddBench <$> [64, 128, 256, 512, 1024], + bgroup "matrix multiplication" $ sizedMatrixMultBench <$> [64, 128, 256, 512, 1024] + ] diff --git a/benchmark/time/Bench.hs b/benchmark/time/Bench.hs deleted file mode 100644 --- a/benchmark/time/Bench.hs +++ /dev/null @@ -1,44 +0,0 @@ -{-| -Module : Bench -Description : Benchmark of Multilinear library -Copyright : (c) Artur M. Brodzki, 2018 -License : BSD3 -Maintainer : artur@brodzki.org -Stability : experimental -Portability : Windows/POSIX - --} - -module Main ( - main -) where - -import Criterion.Main -import qualified Multilinear.Matrix as Matrix - --- | Simple generator function for bencharking matrices -gen :: Int -> Int -> Double -gen j k = sin (fromIntegral j) + cos (fromIntegral k) - --- | Generate benchmark of matrix multiplication -sizedMatrixMultBench :: - Int -- ^ size of square matrix to multiplicate - -> Benchmark -sizedMatrixMultBench s = - bench ((show s) ++ "x" ++ (show s)) $ - nf ((Matrix.fromIndices "ij" s s gen) *) (Matrix.fromIndices "jk" s s gen) - --- | Generate benchmark of matrix addition -sizedMatrixAddBench :: - Int -- ^ size of square matrix to add - -> Benchmark -sizedMatrixAddBench s = - bench ((show s) ++ "x" ++ (show s)) $ - nf ((Matrix.fromIndices "ij" s s gen) +) (Matrix.fromIndices "ij" s s gen) - --- | ENTRY POINT -main :: IO () -main = defaultMain [ - bgroup "matrix multiplication" $ sizedMatrixMultBench <$> [64, 128, 256, 512], - bgroup "matrix addition" $ sizedMatrixAddBench <$> [64, 128, 256, 512] - ] diff --git a/multilinear.cabal b/multilinear.cabal --- a/multilinear.cabal +++ b/multilinear.cabal @@ -2,10 +2,10 @@ -- -- see: https://github.com/sol/hpack -- --- hash: d2d892ab946256d37232b1576c7d134de32865ac69b2c022d7d67f404bf52401 +-- hash: 033c8dbef9eab24df7b3396f684551d557968e2a8963865f404a7ca247732d7d name: multilinear -version: 0.3.2.0 +version: 0.4.0.0 synopsis: Comprehensive and efficient (multi)linear algebra implementation. description: Comprehensive and efficient (multi)linear algebra implementation, based on generic tensor formalism and concise Ricci-Curbastro index syntax. More information available on GitHub: category: Machine learning @@ -32,9 +32,10 @@ Multilinear.Class Multilinear.Form Multilinear.Generic + Multilinear.Generic.MultiCore + Multilinear.Generic.Sequential Multilinear.Index Multilinear.Index.Finite - Multilinear.Index.Infinite Multilinear.Matrix Multilinear.NForm Multilinear.NVector @@ -44,27 +45,25 @@ Paths_multilinear hs-source-dirs: src - default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving + default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs GeneralizedNewtypeDeriving LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving ghc-options: -O2 -W build-depends: base >=4.7 && <5 , containers , deepseq - , mwc-random - , primitive - , statistics + , parallel , vector default-language: Haskell2010 -test-suite all +test-suite multicore type: exitcode-stdio-1.0 main-is: Spec.hs other-modules: - Test.QuickCheck.Multilinear + Test.QuickCheck.Multilinear.Generic.MultiCore Paths_multilinear hs-source-dirs: - test - default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving + test/multicore + default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs GeneralizedNewtypeDeriving LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving ghc-options: -O2 -W -threaded -rtsopts -with-rtsopts=-N build-depends: QuickCheck @@ -76,14 +75,34 @@ , quickcheck-instances default-language: Haskell2010 -benchmark memory +test-suite sequential type: exitcode-stdio-1.0 + main-is: Spec.hs + other-modules: + Test.QuickCheck.Multilinear.Generic.Sequential + Paths_multilinear + hs-source-dirs: + test/sequential + default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs GeneralizedNewtypeDeriving LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving + ghc-options: -O2 -W -threaded -rtsopts -with-rtsopts=-N + build-depends: + QuickCheck + , base >=4.7 && <5 + , containers + , deepseq + , generic-random + , multilinear + , quickcheck-instances + default-language: Haskell2010 + +benchmark multicore-memory + type: exitcode-stdio-1.0 main-is: Bench.hs other-modules: Paths_multilinear hs-source-dirs: - benchmark/memory - default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving + benchmark/multicore/memory + default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs GeneralizedNewtypeDeriving LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving ghc-options: -O2 -W -threaded -rtsopts -with-rtsopts=-N build-depends: base >=4.7 && <5 @@ -91,14 +110,14 @@ , weigh default-language: Haskell2010 -benchmark profile +benchmark multicore-profile type: exitcode-stdio-1.0 main-is: Bench.hs other-modules: Paths_multilinear hs-source-dirs: - benchmark/profile - default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving + benchmark/multicore/profile + default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs GeneralizedNewtypeDeriving LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving ghc-options: -O2 -W -threaded -rtsopts -with-rtsopts=-N build-depends: base >=4.7 && <5 @@ -106,15 +125,60 @@ , multilinear default-language: Haskell2010 -benchmark time +benchmark multicore-time type: exitcode-stdio-1.0 main-is: Bench.hs other-modules: Paths_multilinear hs-source-dirs: - benchmark/time - default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving + benchmark/multicore/time + default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs GeneralizedNewtypeDeriving LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving ghc-options: -O2 -W -threaded -rtsopts -with-rtsopts=-N + build-depends: + base >=4.7 && <5 + , criterion + , multilinear + default-language: Haskell2010 + +benchmark sequential-memory + type: exitcode-stdio-1.0 + main-is: Bench.hs + other-modules: + Paths_multilinear + hs-source-dirs: + benchmark/sequential/memory + default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs GeneralizedNewtypeDeriving LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving + ghc-options: -O2 -W + build-depends: + base >=4.7 && <5 + , multilinear + , weigh + default-language: Haskell2010 + +benchmark sequential-profile + type: exitcode-stdio-1.0 + main-is: Bench.hs + other-modules: + Paths_multilinear + hs-source-dirs: + benchmark/sequential/profile + default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs GeneralizedNewtypeDeriving LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving + ghc-options: -O2 -W + build-depends: + base >=4.7 && <5 + , deepseq + , multilinear + default-language: Haskell2010 + +benchmark sequential-time + type: exitcode-stdio-1.0 + main-is: Bench.hs + other-modules: + Paths_multilinear + hs-source-dirs: + benchmark/sequential/time + default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs GeneralizedNewtypeDeriving LambdaCase MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving + ghc-options: -O2 -W build-depends: base >=4.7 && <5 , criterion diff --git a/src/Multilinear.hs b/src/Multilinear.hs --- a/src/Multilinear.hs +++ b/src/Multilinear.hs @@ -142,49 +142,12 @@ -} module Multilinear ( - module Multilinear.Class - {-module Form, - module Generic, - module Index, - module Index.Finite, - module Index.Infinite, - module Matrix, - module NForm, - module Tensor, - module Vector, - module X-} + module Multilinear.Class, + module Multilinear.Generic, + module Multilinear.Generic.Sequential ) where --- Re-export other library modules +-- Re-export basic library modules import Multilinear.Class -{-import qualified Multilinear.Form as Form -import Multilinear.Generic as Generic -import qualified Multilinear.Index as Index -import qualified Multilinear.Index.Finite as Index.Finite -import qualified Multilinear.Index.Infinite as Index.Infinite -import qualified Multilinear.Matrix as Matrix -import qualified Multilinear.NForm as NForm -import qualified Multilinear.NVector as NVector -import qualified Multilinear.Tensor as Tensor -import qualified Multilinear.Vector as Vector - -import Statistics.Distribution as X -import Statistics.Distribution.Beta as X -import Statistics.Distribution.Binomial as X -import Statistics.Distribution.CauchyLorentz as X -import Statistics.Distribution.ChiSquared as X -import Statistics.Distribution.Exponential as X -import Statistics.Distribution.FDistribution as X -import Statistics.Distribution.Gamma as X -import Statistics.Distribution.Geometric as X -import Statistics.Distribution.Hypergeometric as X -import Statistics.Distribution.Laplace as X -import Statistics.Distribution.Normal as X -import Statistics.Distribution.StudentT as X -import Statistics.Distribution.Uniform as X - -import System.IO.Unsafe as X -import Control.Monad.Trans.Class as X -import Control.Monad.Trans.Either as X -import Control.Exception.Base as X --} +import Multilinear.Generic +import Multilinear.Generic.Sequential diff --git a/src/Multilinear/Class.hs b/src/Multilinear/Class.hs --- a/src/Multilinear/Class.hs +++ b/src/Multilinear/Class.hs @@ -143,8 +143,7 @@ -} module Multilinear.Class ( - Multilinear(..), - Accessible(..) + Multilinear(..) ) where import Data.Maybe @@ -154,33 +153,41 @@ {-| Multidimensional array treated as multilinear map - tensor -} class ( - Unboxed.Unbox a, - Num (t a) -- Tensors may be added, subtracted and multiplicated + Unboxed.Unbox a ) => Multilinear t a where - {-| Add scalar @a@ to each element of tensor @t@ -} - infixl 4 +. - (+.) :: a -> t a -> t a - - {-| Subtract each element of tensor @t@ from scalar scalar left -} - infixl 4 -. - (-.) :: a -> t a -> t a - - {-| Multiply scalar @a@ by each element of tensor @t@ -} - infixl 5 *. - (*.) :: a -> t a -> t a + {-| Generic tensor constructor, using combinator function on its indices -} + fromIndices :: + String -- ^ Upper indices names (one character per index) + -> String -- ^ Lower indices names (one character per index) + -> [Int] -- ^ Upper indices sizes + -> [Int] -- ^ Lower indices sizes + -> ([Int] -> [Int] -> a) -- ^ Combinator function (f [u1,u2,...] [d1,d2,...] returns a tensor element at t [u1,u2,...] [d1,d2,...]) + -> t a -- ^ Generated tensor - {-| Add each element of tensor @t@ to scalar @a@ -} - infixl 4 .+ - (.+) :: t a -> a -> t a + {-| Tensor constructor, that returns tensor with all elments equal to v. -} + {-# INLINE const #-} + const :: + String -- ^ Upper indices names (one character per index) + -> String -- ^ Lower indices names (one character per index) + -> [Int] -- ^ Upper indices sizes + -> [Int] -- ^ Lower indices sizes + -> a -- Value of tensor elements + -> t a -- ^ Generated tensor + const u d usize dsize v = fromIndices u d usize dsize (\_ _ -> v) - {-| Subtract scalar @a@ from each element of tensor @t@ -} - infixl 4 .- - (.-) :: t a -> a -> t a + {-| Accessing tensor elements -} + {-| @el ["i","j"] t [4,5]@ returns all tensor elements which index @i@ is equal to 4 and index @j@ is equal to 5. + Values of other indices are insignificant -} + {-| If given index value is out of range, then modulo operation is performed: + el ["i","j"] t [40 50] = t[40 mod size i, 50 mod size j] -} + el :: t a -> (String,[Int]) -> t a - {-| Multiply each element of tensor @t@ by scalar @a@ -} - infixl 5 .* - (.*) :: t a -> a -> t a + {-| Infix equivalent for el -} + {-# INLINE ($$|) #-} + infixl 7 $$| + ($$|) :: t a -> (String,[Int]) -> t a + t $$| is = el t is {-| List of all tensor indices -} indices :: t a -> [TIndex] @@ -304,54 +311,3 @@ infixl 6 <<<| (<<<|) :: t a -> String -> t a t <<<| n = shiftLeftmost t n - - {-| Simple mapping -} - {-| @map f t@ returns tensor @t2@ in which @t2[i1,i2,...] = f t[i1,i2,...]@ -} - map :: Unboxed.Unbox b => (a -> b) -> t a -> t b - - {-| Filtering tensor. - Filtering multi-dimensional arrray may be dangerous, as we always assume, - that on each recursion level, all tensors have the same size (length). - To disable invalid filters, filtering is done over indices, not tensor elements. - Filter function takes and index name and index value and if it returns True, this index value remains in result tensor. - This allows to remove whole columns or rows of eg. a matrix: - filter (\i n -> i /= "a" || i > 10) filters all rows of "a" index (because if i /= "a", filter returns True) - and for "a" index filter elements with index value <= 10 - But this disallow to remove particular matrix element. - If for some index all elements are removed, the index itself is removed from tensor. -} - filter :: - (String -> Int -> Bool) -- ^ filter function - -> t a -- ^ tensor to filter - -> t a - - {-| Filtering one index of tensor. -} - filterIndex :: - String -- ^ Index name to filter - -> (Int -> Bool) -- ^ filter function - -> t a -- ^ tensor to filter - -> t a - filterIndex iname f = Multilinear.Class.filter (\i n -> i /= iname || f n) - - {-| Zip tensors with binary combinator -} - zipWith :: ( - Unboxed.Unbox b, Unboxed.Unbox c - ) => (a -> b -> c) - -> t a - -> t b - -> t c - -{-| If container on which tensor instance is built, allows for random access of its elements, then the tensor can be instanced as Accessible -} -class Multilinear t a => Accessible t a where - - {-| Accessing tensor elements -} - {-| @el ["i","j"] t [4,5]@ returns all tensor elements which index @i@ is equal to 4 and index @j@ is equal to 5. - Values of other indices are insignificant -} - {-| If given index value is out of range, then modulo operation is performed: - el ["i","j"] t [40 50] = t[40 mod size i, 50 mod size j] -} - el :: t a -> (String,[Int]) -> t a - - {-| Infix equivalent for el -} - {-# INLINE ($$|) #-} - infixl 7 $$| - ($$|) :: t a -> (String,[Int]) -> t a - t $$| is = el t is diff --git a/src/Multilinear/Form.hs b/src/Multilinear/Form.hs --- a/src/Multilinear/Form.hs +++ b/src/Multilinear/Form.hs @@ -14,133 +14,32 @@ module Multilinear.Form ( -- * Generators - -- ** Finite functionals Multilinear.Form.fromIndices, - Multilinear.Form.const, - Multilinear.Form.randomDouble, - Multilinear.Form.randomDoubleSeed, - Multilinear.Form.randomInt, - Multilinear.Form.randomIntSeed + Multilinear.Form.const ) where -import Control.Monad.Primitive import qualified Data.Vector.Unboxed as Unboxed -import Multilinear.Generic -import Multilinear.Tensor as Tensor -import Statistics.Distribution +import Multilinear invalidIndices :: String invalidIndices = "Indices and its sizes not compatible with structure of linear functional!" --- * Finite functional generators - {-| Generate linear functional as function of indices -} - fromIndices :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> (Int -> a) -- ^ Generator function - returns a linear functional component at index @i@ - -> Tensor a -- ^ Generated linear functional - -fromIndices [i] s f = Tensor.fromIndices ([],[]) ([i],[s]) $ \[] [x] -> f x + Num a, Unboxed.Unbox a, Multilinear t a + ) => String -- ^ Index name (one character) + -> Int -- ^ Number of elements + -> (Int -> a) -- ^ Generator function - returns a linear functional component at index @i@ + -> t a -- ^ Generated linear functional +fromIndices [i] s f = Multilinear.fromIndices [] [i] [] [s] $ \[] [x] -> f x fromIndices _ _ _ = error invalidIndices {-| Generate linear functional with all components equal to some @v@ -} - const :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> a -- ^ Value of each element - -> Tensor a -- ^ Generated linear functional - -const [i] s = Tensor.const ([],[]) ([i],[s]) + Num a, Unboxed.Unbox a, Multilinear t a + ) => String -- ^ Index name (one character) + -> Int -- ^ Number of elements + -> a -- ^ Value of each element + -> t a -- ^ Generated linear functional +const [i] s = Multilinear.const [] [i] [] [s] const _ _ = \_ -> error invalidIndices - -{-| Generate linear functional with random real components with given probability distribution. -The functional is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDouble :: ( - ContGen d - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Double) -- ^ Generated linear functional - -randomDouble [i] s = Tensor.randomDouble ([],[]) ([i],[s]) -randomDouble _ _ = \_ -> return $ error invalidIndices - -{-| Generate linear functional with random integer components with given probability distribution. -The functional is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomInt :: ( - DiscreteGen d - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Int) -- ^ Generated linear functional - -randomInt [i] s = Tensor.randomInt ([],[]) ([i],[s]) -randomInt _ _ = \_ -> return $ error invalidIndices - -{-| Generate linear functional with random real components with given probability distribution and given seed. -The functional is wrapped in a monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDoubleSeed :: ( - ContGen d, PrimMonad m - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Double) -- ^ Generated linear functional - -randomDoubleSeed [i] s = Tensor.randomDoubleSeed ([],[]) ([i],[s]) -randomDoubleSeed _ _ = \_ _ -> return $ error invalidIndices - -{-| Generate linear functional with random integer components with given probability distribution and given seed. -The functional is wrapped in a monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomIntSeed :: ( - DiscreteGen d, PrimMonad m - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Int) -- ^ Generated linear functional - -randomIntSeed [i] s = Tensor.randomIntSeed ([],[]) ([i],[s]) -randomIntSeed _ _ = \_ _ -> return $ error invalidIndices - diff --git a/src/Multilinear/Generic.hs b/src/Multilinear/Generic.hs --- a/src/Multilinear/Generic.hs +++ b/src/Multilinear/Generic.hs @@ -1,704 +1,16 @@ {-| -Module : Multilinear.Generic.AsArray -Description : Generic array tensor +Module : Multilinear.Generic +Description : Re-export default tensor implementation Copyright : (c) Artur M. Brodzki, 2018 License : BSD3 Maintainer : artur@brodzki.org Stability : experimental Portability : Windows/POSIX -- This module contains generic implementation of tensor defined as nested arrays - -} module Multilinear.Generic ( - Tensor(..), (!), - isScalar, isSimple, isFiniteTensor, - dot, _elemByElem, contractionErr, tensorIndex, _standardize + module DefaultTensorImplementation ) where -import Control.DeepSeq -import Data.Foldable -import Data.List -import Data.Maybe -import qualified Data.Vector as Boxed -import qualified Data.Vector.Unboxed as Unboxed -import GHC.Generics -import Multilinear.Class as Multilinear -import qualified Multilinear.Index as Index -import qualified Multilinear.Index.Finite as Finite - -{-| ERROR MESSAGE -} -incompatibleTypes :: String -incompatibleTypes = "Incompatible tensor types!" - -{-| ERROR MESSAGE -} -scalarIndices :: String -scalarIndices = "Scalar has no indices!" - -{-| ERROR MESSAGE -} -indexNotFound :: String -indexNotFound = "This tensor has not such index!" - -{-| ERROR MESSAGE -} -tensorOfScalars :: String -tensorOfScalars = "Tensor construction error! Vector of scalars" - -{-| Tensor defined recursively as scalar or list of other tensors -} -{-| @c@ is type of a container, @i@ is type of index size and @a@ is type of tensor elements -} -data Tensor a where - {-| Scalar -} - Scalar :: { - scalarVal :: a - } -> Tensor a - {-| Simple, one-dimensional finite tensor -} - SimpleFinite :: { - tensorFiniteIndex :: Finite.Index, - tensorScalars :: Unboxed.Vector a - } -> Tensor a - {-| Finite array of other tensors -} - FiniteTensor :: { - {-| Finite index "Mutltilinear.Index.Finite" of tensor -} - tensorFiniteIndex :: Finite.Index, - {-| Array of tensors on deeper recursion level -} - tensorsFinite :: Boxed.Vector (Tensor a) - } -> Tensor a - deriving (Eq, Generic) - - -{-| Return true if tensor is a scalar -} -{-# INLINE isScalar #-} -isScalar :: Unboxed.Unbox a => Tensor a -> Bool -isScalar x = case x of - Scalar _ -> True - _ -> False - -{-| Return true if tensor is a simple tensor -} -{-# INLINE isSimple #-} -isSimple :: Unboxed.Unbox a => Tensor a -> Bool -isSimple x = case x of - SimpleFinite _ _ -> True - _ -> False - -{-| Return True if tensor is a complex tensor -} -{-# INLINE isFiniteTensor #-} -isFiniteTensor :: Unboxed.Unbox a => Tensor a -> Bool -isFiniteTensor x = case x of - FiniteTensor _ _ -> True - _ -> False - -{-| Return generic tensor index -} -{-# INLINE tensorIndex #-} -tensorIndex :: Unboxed.Unbox a => Tensor a -> Index.TIndex -tensorIndex x = case x of - Scalar _ -> error scalarIndices - SimpleFinite i _ -> Index.toTIndex i - FiniteTensor i _ -> Index.toTIndex i - -{-| Returns sample tensor on deeper recursion level.Used to determine some features common for all tensors -} -{-# INLINE firstTensor #-} -firstTensor :: Unboxed.Unbox a => Tensor a -> Tensor a -firstTensor x = case x of - FiniteTensor _ ts -> Boxed.head ts - _ -> x - -{-| Recursive indexing on list tensor - @t ! i = t[i]@ -} -{-# INLINE (!) #-} -(!) :: Unboxed.Unbox a => Tensor a -- ^ tensor @t@ - -> Int -- ^ index @i@ - -> Tensor a -- ^ tensor @t[i]@ -t ! i = case t of - Scalar _ -> error scalarIndices - SimpleFinite ind ts -> - if i >= Finite.indexSize ind then - error ("Index + " ++ show ind ++ " out of bonds!") - else Scalar $ ts Unboxed.! i - FiniteTensor ind ts -> - if i >= Finite.indexSize ind then - error ("Index + " ++ show ind ++ " out of bonds!") - else ts Boxed.! i - --- | NFData instance -instance NFData a => NFData (Tensor a) - --- | move contravariant indices to lower recursion level -_standardize :: (Num a, Unboxed.Unbox a, NFData a) => Tensor a -> Tensor a -_standardize tens = foldr' f tens $ indices tens - where - f i t = if Index.isContravariant i then - t <<<| Index.indexName i - else t - --- | Print tensor -instance ( - Unboxed.Unbox a, Show a, Num a, NFData a - ) => Show (Tensor a) where - - -- merge errors first and then print whole tensor - show = show' . _standardize - where - show' x = case x of - -- Scalar is showed simply as its value - Scalar v -> show v - -- SimpleFinite is shown dependent on its index... - SimpleFinite index ts -> show index ++ "S: " ++ case index of - -- If index is contravariant, show tensor components vertically - Finite.Contravariant _ _ -> "\n" ++ tail (Unboxed.foldl' (\string e -> string ++ "\n |" ++ show e) "" ts) - -- If index is covariant or indifferent, show tensor compoments horizontally - _ -> "[" ++ tail (Unboxed.foldl' (\string e -> string ++ "," ++ show e) "" ts) ++ "]" - -- FiniteTensor is shown dependent on its index... - FiniteTensor index ts -> show index ++ "T: " ++ case index of - -- If index is contravariant, show tensor components vertically - Finite.Contravariant _ _ -> "\n" ++ tail (Boxed.foldl' (\string e -> string ++ "\n |" ++ show e) "" ts) - -- If index is covariant or indifferent, show tensor compoments horizontally - _ -> "[" ++ tail (Boxed.foldl' (\string e -> string ++ "," ++ show e) "" ts) ++ "]" - -{-| Merge FiniteTensor of Scalars to SimpleFinite tensor for performance improvement -} -{-# INLINE _mergeScalars #-} -_mergeScalars :: Unboxed.Unbox a => Tensor a -> Tensor a -_mergeScalars x = case x of - (FiniteTensor index1 ts1) -> case ts1 Boxed.! 0 of - Scalar _ -> SimpleFinite index1 $ Unboxed.generate (Boxed.length ts1) (\i -> scalarVal (ts1 Boxed.! i)) - _ -> FiniteTensor index1 $ _mergeScalars <$> ts1 - _ -> x - --- | Generic map function, which does not require a,b types to be Num -_map :: ( - Unboxed.Unbox a, Unboxed.Unbox b - ) => (a -> b) - -> Tensor a - -> Tensor b -_map f x = case x of - -- Mapping scalar simply maps its value - Scalar v -> Scalar $ f v - -- Mapping complex tensor does mapping element by element - SimpleFinite index ts -> SimpleFinite index (f `Unboxed.map` ts) - FiniteTensor index ts -> FiniteTensor index $ _map f <$> ts - -{-| Apply a tensor operator (here denoted by (+) ) elem by elem, trying to connect as many common indices as possible -} -{-# INLINE _elemByElem' #-} -_elemByElem' :: (Num a, Unboxed.Unbox a, NFData a) - => Tensor a -- ^ First argument of operator - -> Tensor a -- ^ Second argument of operator - -> (a -> a -> a) -- ^ Operator on tensor elements if indices are different - -> (Tensor a -> Tensor a -> Tensor a) -- ^ Tensor operator called if indices are the same - -> Tensor a -- ^ Result tensor - --- @Scalar x + Scalar y = Scalar x + y@ -_elemByElem' (Scalar x1) (Scalar x2) f _ = Scalar $ f x1 x2 --- @Scalar x + Tensor t[i] = Tensor r[i] | r[i] = x + t[i]@ -_elemByElem' (Scalar x) t f _ = (x `f`) `Multilinear.map` t --- @Tensor t[i] + Scalar x = Tensor r[i] | r[i] = t[i] + x@ -_elemByElem' t (Scalar x) f _ = (`f` x) `Multilinear.map` t - --- Two simple tensors case -_elemByElem' t1@(SimpleFinite index1 v1) t2@(SimpleFinite index2 _) f op - | Index.indexName index1 == Index.indexName index2 = op t1 t2 - | otherwise = FiniteTensor index1 $ Boxed.generate (Unboxed.length v1) - (\i -> (\x -> f x `Multilinear.map` t2) (v1 Unboxed.! i)) - --- Two finite tensors case -_elemByElem' t1@(FiniteTensor index1 v1) t2@(FiniteTensor index2 v2) f op - | Index.indexName index1 == Index.indexName index2 = op t1 t2 - | Index.indexName index1 `Data.List.elem` indicesNames t2 = - FiniteTensor index2 $ (\x -> _elemByElem' t1 x f op) <$> v2 - | otherwise = FiniteTensor index1 $ (\x -> _elemByElem' x t2 f op) <$> v1 - --- Simple and finite tensor case -_elemByElem' t1@(SimpleFinite index1 _) t2@(FiniteTensor index2 v2) f op - | Index.indexName index1 == Index.indexName index2 = op t1 t2 - | otherwise = FiniteTensor index2 $ (\x -> _elemByElem' t1 x f op) <$> v2 - --- Finite and simple tensor case -_elemByElem' t1@(FiniteTensor index1 v1) t2@(SimpleFinite index2 _) f op - | Index.indexName index1 == Index.indexName index2 = op t1 t2 - | otherwise = FiniteTensor index1 $ (\x -> _elemByElem' x t2 f op) <$> v1 - -{-| Apply a tensor operator elem by elem and merge scalars to simple tensor at the and -} -{-# INLINE _elemByElem #-} -_elemByElem :: (Num a, Unboxed.Unbox a, NFData a) - => Tensor a -- ^ First argument of operator - -> Tensor a -- ^ Second argument of operator - -> (a -> a -> a) -- ^ Operator on tensor elements if indices are different - -> (Tensor a -> Tensor a -> Tensor a) -- ^ Tensor operator called if indices are the same - -> Tensor a -- ^ Result tensor -_elemByElem t1 t2 f op = - let commonIndices = Data.List.filter (`Data.List.elem` indicesNames t2) $ indicesNames t1 - t1' = foldl' (|>>>) t1 commonIndices - t2' = foldl' (|>>>) t2 commonIndices - in _mergeScalars $ _elemByElem' t1' t2' f op - --- | Zipping two tensors with a combinator, assuming they have the same indices. -{-# INLINE zipT #-} -zipT :: ( - Num a, NFData a, Unboxed.Unbox a - ) => (a -> a -> a) -- ^ The zipping combinator - -> Tensor a -- ^ First tensor to zip - -> Tensor a -- ^ Second tensor to zip - -> Tensor a -- ^ Result tensor - --- Two simple tensors case -zipT f t1@(SimpleFinite index1 v1) t2@(SimpleFinite index2 v2) = - if index1 == index2 then - SimpleFinite index1 $ Unboxed.zipWith f v1 v2 - else dot t1 t2 - ---Two finite tensors case -zipT f t1@(FiniteTensor index1 v1) t2@(FiniteTensor index2 v2) = - if index1 == index2 then - FiniteTensor index1 $ Boxed.zipWith (zipT f) v1 v2 - else dot t1 t2 - --- Finite and simple tensor case -zipT f t1@(FiniteTensor index1 v1) t2@(SimpleFinite index2 v2) = - if index1 == index2 then - let f' t s = (`f` s) `_map` t - in FiniteTensor index1 $ Boxed.generate (Finite.indexSize index1) (\i -> f' (v1 Boxed.! i) (v2 Unboxed.! i)) - else dot t1 t2 - --- Simple and finite tensor case -zipT f t1@(SimpleFinite index1 v1) t2@(FiniteTensor index2 v2) = - if index1 == index2 then - let f' s t = (s `f`) `_map` t - in FiniteTensor index1 $ Boxed.generate (Finite.indexSize index1) (\i -> f' (v1 Unboxed.! i) (v2 Boxed.! i)) - else dot t1 t2 - --- Zipping something with scalar is impossible -zipT _ _ _ = error $ "zipT: " ++ scalarIndices - --- | zipT error -{-# INLINE zipErr #-} -zipErr :: String -- ^ zipT function variant where the error occured - -> Index.TIndex -- ^ Index of first dot product parameter - -> Index.TIndex -- ^ Index of second dot product parameter - -> Tensor a -- ^ Erorr message -zipErr variant i1' i2' = error $ - "zipT: " ++ variant ++ " - " ++ incompatibleTypes ++ - " - index1 is " ++ show i1' ++ - " and index2 is " ++ show i2' - - --- | dot product of two tensors -{-# INLINE dot #-} -dot :: (Num a, Unboxed.Unbox a, NFData a) - => Tensor a -- ^ First dot product argument - -> Tensor a -- ^ Second dot product argument - -> Tensor a -- ^ Resulting dot product - --- Two simple tensors product -dot (SimpleFinite i1@(Finite.Covariant count1 _) ts1') (SimpleFinite i2@(Finite.Contravariant count2 _) ts2') - | count1 == count2 = - Scalar $ Unboxed.sum $ Unboxed.zipWith (*) ts1' ts2' - | otherwise = contractionErr "simple-simple" (Index.toTIndex i1) (Index.toTIndex i2) -dot (SimpleFinite i1@(Finite.Contravariant count1 _) ts1') (SimpleFinite i2@(Finite.Covariant count2 _) ts2') - | count1 == count2 = - Scalar $ Unboxed.sum $ Unboxed.zipWith (*) ts1' ts2' - | otherwise = contractionErr "simple-simple" (Index.toTIndex i1) (Index.toTIndex i2) -dot t1@(SimpleFinite _ _) t2@(SimpleFinite _ _) = zipT (*) t1 t2 - --- Two finite tensors product -dot (FiniteTensor i1@(Finite.Covariant count1 _) ts1') (FiniteTensor i2@(Finite.Contravariant count2 _) ts2') - | count1 == count2 = Boxed.sum $ Boxed.zipWith (*) ts1' ts2' - | otherwise = contractionErr "finite-finite" (Index.toTIndex i1) (Index.toTIndex i2) -dot (FiniteTensor i1@(Finite.Contravariant count1 _) ts1') (FiniteTensor i2@(Finite.Covariant count2 _) ts2') - | count1 == count2 = Boxed.sum $ Boxed.zipWith (*) ts1' ts2' - | otherwise = contractionErr "finite-finite" (Index.toTIndex i1) (Index.toTIndex i2) -dot t1@(FiniteTensor _ _) t2@(FiniteTensor _ _) = zipT (*) t1 t2 - --- Simple tensor and finite tensor product -dot (SimpleFinite i1@(Finite.Covariant count1 _) ts1') (FiniteTensor i2@(Finite.Contravariant count2 _) ts2') - | count1 == count2 = Boxed.sum $ Boxed.generate count1 (\i -> (ts1' Unboxed.! i) *. (ts2' Boxed.! i)) - | otherwise = contractionErr "simple-finite" (Index.toTIndex i1) (Index.toTIndex i2) -dot (SimpleFinite i1@(Finite.Contravariant count1 _) ts1') (FiniteTensor i2@(Finite.Covariant count2 _) ts2') - | count1 == count2 = Boxed.sum $ Boxed.generate count1 (\i -> (ts1' Unboxed.! i) *. (ts2' Boxed.! i)) - | otherwise = contractionErr "simple-finite" (Index.toTIndex i1) (Index.toTIndex i2) -dot t1@(SimpleFinite _ _) t2@(FiniteTensor _ _) = zipT (*) t1 t2 - --- Finite tensor and simple tensor product -dot (FiniteTensor i1@(Finite.Covariant count1 _) ts1') (SimpleFinite i2@(Finite.Contravariant count2 _) ts2') - | count1 == count2 = Boxed.sum $ Boxed.generate count1 (\i -> (ts1' Boxed.! i) .* (ts2' Unboxed.! i)) - | otherwise = contractionErr "finite-simple" (Index.toTIndex i1) (Index.toTIndex i2) -dot (FiniteTensor i1@(Finite.Contravariant count1 _) ts1') (SimpleFinite i2@(Finite.Covariant count2 _) ts2') - | count1 == count2 = Boxed.sum $ Boxed.generate count1 (\i -> (ts1' Boxed.! i) .* (ts2' Unboxed.! i)) - | otherwise = contractionErr "finite-simple" (Index.toTIndex i1) (Index.toTIndex i2) -dot t1@(FiniteTensor _ _) t2@(SimpleFinite _ _) = zipT (*) t1 t2 - --- Other cases cannot happen! -dot t1' t2' = contractionErr "other" (tensorIndex t1') (tensorIndex t2') - --- | contraction error -{-# INLINE contractionErr #-} -contractionErr :: String -- ^ dot function variant where error occured - -> Index.TIndex -- ^ Index of first dot product parameter - -> Index.TIndex -- ^ Index of second dot product parameter - -> Tensor a -- ^ Erorr message - -contractionErr variant i1' i2' = error $ - "Tensor product: " ++ variant ++ " - " ++ incompatibleTypes ++ - " - index1 is " ++ show i1' ++ - " and index2 is " ++ show i2' - -{-| Transpose Vector of Vectors, analogous to Data.List.transpose function. It is assumed, that all vectors on deeper recursion level have the same length. -} -_transpose :: Boxed.Vector (Boxed.Vector a) -- ^ Vector of vectors to transpose - -> Boxed.Vector (Boxed.Vector a) -_transpose v = - let outerS = Boxed.length v - innerS = Boxed.length $ v Boxed.! 0 - in Boxed.generate innerS (\i -> Boxed.generate outerS $ \j -> v Boxed.! j Boxed.! i) - --- | Tensors can be added, subtracted and multiplicated -instance (Unboxed.Unbox a, Num a, NFData a) => Num (Tensor a) where - - -- Adding - element by element - {-# INLINE (+) #-} - t1 + t2 = _elemByElem t1 t2 (+) $ zipT (+) - - -- Subtracting - element by element - {-# INLINE (-) #-} - t1 - t2 = _elemByElem t1 t2 (-) $ zipT (-) - - -- Multiplicating is treated as tensor product - -- Tensor product applies Einstein summation convention - {-# INLINE (*) #-} - t1 * t2 = _elemByElem t1 t2 (*) dot - - -- Absolute value - element by element - {-# INLINE abs #-} - abs t = abs `Multilinear.map` t - - -- Signum operation - element by element - {-# INLINE signum #-} - signum t = signum `Multilinear.map` t - - -- Simple integer can be conveted to Scalar - {-# INLINE fromInteger #-} - fromInteger x = Scalar $ fromInteger x - --- | Tensors can be divided by each other -instance (Unboxed.Unbox a, Fractional a, NFData a) => Fractional (Tensor a) where - - {-# INLINE (/) #-} - -- Scalar division return result of division of its values - Scalar x1 / Scalar x2 = Scalar $ x1 / x2 - -- Tensor and scalar are divided value by value - Scalar x1 / t2 = (x1 /) `Multilinear.map` t2 - t1 / Scalar x2 = (/ x2) `Multilinear.map` t1 - -- Two complex tensors cannot be (for now) simply divided - -- // TODO - tensor division and inversion - _ / _ = error "TODO" - - -- A scalar can be generated from rational number - {-# INLINE fromRational #-} - fromRational x = Scalar $ fromRational x - --- Real-number functions on tensors. --- Function of tensor is tensor of function of its elements --- E.g. exp [1,2,3,4] = [exp 1, exp2, exp3, exp4] -instance (Unboxed.Unbox a, Floating a, NFData a) => Floating (Tensor a) where - - {-| PI number -} - {-# INLINE pi #-} - pi = Scalar pi - - {-| Exponential function. (exp t)[i] = exp( t[i] ) -} - {-# INLINE exp #-} - exp t = exp `Multilinear.map` t - - {-| Natural logarithm. (log t)[i] = log( t[i] ) -} - {-# INLINE log #-} - log t = log `Multilinear.map` t - - {-| Sinus. (sin t)[i] = sin( t[i] ) -} - {-# INLINE sin #-} - sin t = sin `Multilinear.map` t - - {-| Cosinus. (cos t)[i] = cos( t[i] ) -} - {-# INLINE cos #-} - cos t = cos `Multilinear.map` t - - {-| Inverse sinus. (asin t)[i] = asin( t[i] ) -} - {-# INLINE asin #-} - asin t = asin `Multilinear.map` t - - {-| Inverse cosinus. (acos t)[i] = acos( t[i] ) -} - {-# INLINE acos #-} - acos t = acos `Multilinear.map` t - - {-| Inverse tangent. (atan t)[i] = atan( t[i] ) -} - {-# INLINE atan #-} - atan t = atan `Multilinear.map` t - - {-| Hyperbolic sinus. (sinh t)[i] = sinh( t[i] ) -} - {-# INLINE sinh #-} - sinh t = sinh `Multilinear.map` t - - {-| Hyperbolic cosinus. (cosh t)[i] = cosh( t[i] ) -} - {-# INLINE cosh #-} - cosh t = cosh `Multilinear.map` t - - {-| Inverse hyperbolic sinus. (asinh t)[i] = asinh( t[i] ) -} - {-# INLINE asinh #-} - asinh t = acosh `Multilinear.map` t - - {-| Inverse hyperbolic cosinus. (acosh t)[i] = acosh (t[i] ) -} - {-# INLINE acosh #-} - acosh t = acosh `Multilinear.map` t - - {-| Inverse hyperbolic tangent. (atanh t)[i] = atanh( t[i] ) -} - {-# INLINE atanh #-} - atanh t = atanh `Multilinear.map` t - --- Multilinear operations -instance (Unboxed.Unbox a, Num a, NFData a) => Multilinear Tensor a where - - -- Add scalar right - {-# INLINE (.+) #-} - t .+ x = (+x) `Multilinear.map` t - - -- Subtract scalar right - {-# INLINE (.-) #-} - t .- x = (\p -> p - x) `Multilinear.map` t - - -- Multiplicate by scalar right - {-# INLINE (.*) #-} - t .* x = (*x) `Multilinear.map` t - - -- Add scalar left - {-# INLINE (+.) #-} - x +. t = (x+) `Multilinear.map` t - - -- Subtract scalar left - {-# INLINE (-.) #-} - x -. t = (x-) `Multilinear.map` t - - -- Multiplicate by scalar left - {-# INLINE (*.) #-} - x *. t = (x*) `Multilinear.map` t - - -- List of all tensor indices - {-# INLINE indices #-} - indices x = case x of - Scalar _ -> [] - FiniteTensor i ts -> Index.toTIndex i : indices (head $ toList ts) - SimpleFinite i _ -> [Index.toTIndex i] - - -- Get tensor order [ (contravariant,covariant)-type ] - {-# INLINE order #-} - order x = case x of - Scalar _ -> (0,0) - SimpleFinite index _ -> case index of - Finite.Contravariant _ _ -> (1,0) - Finite.Covariant _ _ -> (0,1) - Finite.Indifferent _ _ -> (0,0) - _ -> let (cnvr, covr) = order $ firstTensor x - in case tensorIndex x of - Index.Contravariant _ _ -> (cnvr+1,covr) - Index.Covariant _ _ -> (cnvr,covr+1) - Index.Indifferent _ _ -> (cnvr,covr) - - -- Get size of tensor index or Left if index is infinite or tensor has no such index - {-# INLINE size #-} - size t iname = case t of - Scalar _ -> error scalarIndices - SimpleFinite index _ -> - if Index.indexName index == iname - then Finite.indexSize index - else error indexNotFound - FiniteTensor index _ -> - if Index.indexName index == iname - then Finite.indexSize index - else size (firstTensor t) iname - - -- Rename tensor indices - {-# INLINE ($|) #-} - - Scalar x $| _ = Scalar x - SimpleFinite (Finite.Contravariant isize _) ts $| (u:_, _) = SimpleFinite (Finite.Contravariant isize [u]) ts - SimpleFinite (Finite.Covariant isize _) ts $| (_, d:_) = SimpleFinite (Finite.Covariant isize [d]) ts - FiniteTensor (Finite.Contravariant isize _) ts $| (u:us, ds) = FiniteTensor (Finite.Contravariant isize [u]) $ ($| (us,ds)) <$> ts - FiniteTensor (Finite.Covariant isize _) ts $| (us, d:ds) = FiniteTensor (Finite.Covariant isize [d]) $ ($| (us,ds)) <$> ts - t $| _ = t - - -- Raise an index - {-# INLINE (/\) #-} - Scalar x /\ _ = Scalar x - FiniteTensor index ts /\ n - | Index.indexName index == n = - FiniteTensor (Finite.Contravariant (Finite.indexSize index) n) $ (/\ n) <$> ts - | otherwise = - FiniteTensor index $ (/\ n) <$> ts - t1@(SimpleFinite index ts) /\ n - | Index.indexName index == n = - SimpleFinite (Finite.Contravariant (Finite.indexSize index) n) ts - | otherwise = t1 - - -- Lower an index - {-# INLINE (\/) #-} - Scalar x \/ _ = Scalar x - FiniteTensor index ts \/ n - | Index.indexName index == n = - FiniteTensor (Finite.Covariant (Finite.indexSize index) n) $ (\/ n) <$> ts - | otherwise = - FiniteTensor index $ (\/ n) <$> ts - t1@(SimpleFinite index ts) \/ n - | Index.indexName index == n = - SimpleFinite (Finite.Covariant (Finite.indexSize index) n) ts - | otherwise = t1 - - {-| Transpose a tensor (switch all indices types) -} - {-# INLINE transpose #-} - transpose (Scalar x) = Scalar x - - transpose (FiniteTensor (Finite.Covariant count name) ts) = - FiniteTensor (Finite.Contravariant count name) (Multilinear.transpose <$> ts) - transpose (FiniteTensor (Finite.Contravariant count name) ts) = - FiniteTensor (Finite.Covariant count name) (Multilinear.transpose <$> ts) - transpose (FiniteTensor (Finite.Indifferent count name) ts) = - FiniteTensor (Finite.Indifferent count name) (Multilinear.transpose <$> ts) - - transpose (SimpleFinite (Finite.Covariant count name) ts) = - SimpleFinite (Finite.Contravariant count name) ts - transpose (SimpleFinite (Finite.Contravariant count name) ts) = - SimpleFinite (Finite.Covariant count name) ts - transpose (SimpleFinite (Finite.Indifferent count name) ts) = - SimpleFinite (Finite.Indifferent count name) ts - - - {-| Shift tensor index right -} - {-| Moves given index one level deeper in recursion -} - -- Scalar has no indices to shift - Scalar x `shiftRight` _ = Scalar x - -- Simple tensor has only one index which cannot be shifted - t1@(SimpleFinite _ _) `shiftRight` _ = t1 - -- Finite tensor is shifted by converting to list and transposing it - t1@(FiniteTensor index1 ts1) `shiftRight` ind - -- We don't shift this index - | Data.List.length (indicesNames t1) > 1 && Index.indexName index1 /= ind = - FiniteTensor index1 $ (|>> ind) <$> ts1 - -- We found index to shift - | Data.List.length (indicesNames t1) > 1 && Index.indexName index1 == ind = - -- Next index - let index2 = tensorFiniteIndex (ts1 Boxed.! 0) - -- Elements to transpose - dane = if isSimple (ts1 Boxed.! 0) - then (\un -> Boxed.generate (Unboxed.length un) (\i -> Scalar $ un Unboxed.! i)) <$> - (tensorScalars <$> ts1) - else tensorsFinite <$> ts1 - -- reconstruct tensor with transposed elements - in _mergeScalars $ FiniteTensor index2 $ FiniteTensor index1 <$> (_transpose dane) - -- there is only one index and therefore it cannot be shifted - | otherwise = t1 - - {-| Map function, as in Functor typeclass -} - {-# INLINE map #-} - map f x = case x of - -- Mapping scalar simply maps its value - Scalar v -> Scalar $ f v - -- Mapping complex tensor does mapping element by element - SimpleFinite index ts -> SimpleFinite index (f `Unboxed.map` ts) - FiniteTensor index ts -> FiniteTensor index $ Multilinear.map f <$> ts - - {-| Filter functions -} - filter _ (Scalar x) = Scalar x - - filter f (SimpleFinite index ts) = - let iname = Finite.indexName' index - ts' = (\i _ -> f iname i) `Unboxed.ifilter` ts - in SimpleFinite index { Finite.indexSize = Unboxed.length ts' } ts' - - filter f (FiniteTensor index ts) = - let iname = Finite.indexName' index - ts' = Multilinear.filter f <$> ((\i _ -> f iname i) `Boxed.ifilter` ts) - ts'' = - (\case - (SimpleFinite _ ts) -> not $ Unboxed.null ts - (FiniteTensor _ ts) -> not $ Boxed.null ts - _ -> error $ "Filter: " ++ tensorOfScalars - ) `Boxed.filter` ts' - in FiniteTensor index { Finite.indexSize = Boxed.length ts'' } ts'' - - {-| Zip tensors with binary combinator -} - - -- Zippin two Scalars simply combines their values - zipWith f (Scalar x1) (Scalar x2) = Scalar $ f x1 x2 - - -- zipping complex tensor with scalar - zipWith f t (Scalar x) = (`f` x) `_map` t - - -- zipping scalar with complex tensor - zipWith f (Scalar x) t = (x `f`) `_map` t - - -- Two simple tensors case - zipWith f (SimpleFinite index1 v1) (SimpleFinite index2 v2) = - if index1 == index2 then - SimpleFinite index1 $ Unboxed.zipWith f v1 v2 - else zipErr "simple-simple" (Index.toTIndex index1) (Index.toTIndex index2) - - --Two finite tensors case - zipWith f (FiniteTensor index1 v1) (FiniteTensor index2 v2) = - if index1 == index2 then - FiniteTensor index1 $ Boxed.zipWith (Multilinear.zipWith f) v1 v2 - else zipErr "finite-finite" (Index.toTIndex index1) (Index.toTIndex index2) - - -- Finite and simple tensor case - zipWith f (FiniteTensor index1 v1) (SimpleFinite index2 v2) = - if index1 == index2 then - let f' t s = (`f` s) `_map` t - in FiniteTensor index1 $ Boxed.generate (Finite.indexSize index1) (\i -> f' (v1 Boxed.! i) (v2 Unboxed.! i)) - else zipErr "finite-simple" (Index.toTIndex index1) (Index.toTIndex index2) - - -- Simple and finite tensor case - zipWith f (SimpleFinite index1 v1) (FiniteTensor index2 v2) = - if index1 == index2 then - let f' s t = (s `f`) `_map` t - in FiniteTensor index1 $ Boxed.generate (Finite.indexSize index1) (\i -> f' (v1 Unboxed.! i) (v2 Boxed.! i)) - else zipErr "simple-finite" (Index.toTIndex index1) (Index.toTIndex index2) - - - - - - - - -{-| List allows for random access to tensor elements -} -instance (Unboxed.Unbox a, Num a, NFData a) => Accessible Tensor a where - - {-| Accessing tensor elements -} - {-# INLINE el #-} - - -- Scalar has only one element - el (Scalar x) _ = Scalar x - - -- simple tensor case - el t1@(SimpleFinite index1 _) (inds,vals) = - -- zip indices with their given values - let indvals = zip inds vals - -- find value for simple tensor index if given - val = Data.List.find (\(n,_) -> [n] == Index.indexName index1) indvals - -- if value for current index is given - in if isJust val - -- then get it from current tensor - then t1 ! snd (fromJust val) - -- otherwise return whole tensor - no filtering defined - else t1 - - -- finite tensor case - el t1@(FiniteTensor index1 v1) (inds,vals) = - -- zip indices with their given values - let indvals = zip inds vals - -- find value for current index if given - val = Data.List.find (\(n,_) -> [n] == Index.indexName index1) indvals - -- and remove used index from indices list - indvals1 = Data.List.filter (\(n,_) -> [n] /= Index.indexName index1) indvals - -- indices unused so far - inds1 = Data.List.map fst indvals1 - -- and its corresponding values - vals1 = Data.List.map snd indvals1 - -- if value for current index was given - in if isJust val - -- then get it from current tensor and recursively process other indices - then el (t1 ! snd (fromJust val)) (inds1,vals1) - -- otherwise recursively access elements of all child tensors - else FiniteTensor index1 $ (\t -> el t (inds,vals)) <$> v1 +import Multilinear.Generic.Sequential as DefaultTensorImplementation diff --git a/src/Multilinear/Generic/MultiCore.hs b/src/Multilinear/Generic/MultiCore.hs new file mode 100644 --- /dev/null +++ b/src/Multilinear/Generic/MultiCore.hs @@ -0,0 +1,765 @@ +{-| +Module : Multilinear.Generic.MultiCore +Description : Generic implementation of tensor as nested arrays, evaluated in sequential manner +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Multilinear.Generic.MultiCore ( + -- * Generic tensor datatype and its instances + Tensor(..), + -- * Auxiliary functions + (!), isScalar, isSimple, isFiniteTensor, + tensorIndex, _standardize, _mergeScalars, + _map, _contractedIndices, _elemByElem, zipT, + -- * Additional functions + (.+), (.-), (.*), (+.), (-.), (*.), + Multilinear.Generic.MultiCore.map, + Multilinear.Generic.MultiCore.filter, + Multilinear.Generic.MultiCore.filterIndex, + Multilinear.Generic.MultiCore.zipWith +) where + +import Control.DeepSeq +import qualified Control.Parallel.Strategies as Parallel +import Data.Foldable +import Data.List +import Data.Maybe +import qualified Data.Set as Set +import qualified Data.Vector as Boxed +import qualified Data.Vector.Unboxed as Unboxed +import GHC.Generics +import Multilinear.Class as Multilinear +import qualified Multilinear.Index as Index +import qualified Multilinear.Index.Finite as Finite + +{-| ERROR MESSAGE -} +incompatibleTypes :: String +incompatibleTypes = "Incompatible tensor types!" + +{-| ERROR MESSAGE -} +scalarIndices :: String +scalarIndices = "Scalar has no indices!" + +{-| ERROR MESSAGE -} +indexNotFound :: String +indexNotFound = "This tensor has not such index!" + +{-| ERROR MESSAGE -} +tensorOfScalars :: String +tensorOfScalars = "Tensor construction error! Vector of scalars" + +{-| Tensor defined recursively as scalar or list of other tensors -} +{-| @c@ is type of a container, @i@ is type of index size and @a@ is type of tensor elements -} +data Tensor a where + {-| Scalar -} + Scalar :: { + scalarVal :: a + } -> Tensor a + {-| Simple, one-dimensional finite tensor -} + SimpleFinite :: { + tensorFiniteIndex :: Finite.Index, + tensorScalars :: Unboxed.Vector a + } -> Tensor a + {-| Finite array of other tensors -} + FiniteTensor :: { + {-| Finite index "Mutltilinear.Index.Finite" of tensor -} + tensorFiniteIndex :: Finite.Index, + {-| Array of tensors on deeper recursion level -} + tensorsFinite :: Boxed.Vector (Tensor a) + } -> Tensor a + deriving (Eq, Generic) + +{-| Return true if tensor is a scalar -} +{-# INLINE isScalar #-} +isScalar :: Unboxed.Unbox a => Tensor a -> Bool +isScalar x = case x of + Scalar _ -> True + _ -> False + +{-| Return true if tensor is a simple tensor -} +{-# INLINE isSimple #-} +isSimple :: Unboxed.Unbox a => Tensor a -> Bool +isSimple x = case x of + SimpleFinite _ _ -> True + _ -> False + +{-| Return True if tensor is a complex tensor -} +{-# INLINE isFiniteTensor #-} +isFiniteTensor :: Unboxed.Unbox a => Tensor a -> Bool +isFiniteTensor x = case x of + FiniteTensor _ _ -> True + _ -> False + +{-| Return generic tensor index -} +{-# INLINE tensorIndex #-} +tensorIndex :: Unboxed.Unbox a => Tensor a -> Index.TIndex +tensorIndex x = case x of + Scalar _ -> error scalarIndices + SimpleFinite i _ -> Index.toTIndex i + FiniteTensor i _ -> Index.toTIndex i + +{-| Returns sample tensor on deeper recursion level.Used to determine some features common for all tensors -} +{-# INLINE firstTensor #-} +firstTensor :: Unboxed.Unbox a => Tensor a -> Tensor a +firstTensor x = case x of + FiniteTensor _ ts -> Boxed.head ts + _ -> x + +{-| Recursive indexing on list tensor. If index is greater than index size, performs modulo indexing + @t ! i = t[i]@ -} +{-# INLINE (!) #-} +(!) :: Unboxed.Unbox a => Tensor a -- ^ tensor @t@ + -> Int -- ^ index @i@ + -> Tensor a -- ^ tensor @t[i]@ +t ! i = case t of + Scalar _ -> error scalarIndices + SimpleFinite ind ts -> Scalar $ ts Unboxed.! (i `mod` Finite.indexSize ind) + FiniteTensor ind ts -> ts Boxed.! (i `mod` Finite.indexSize ind) + +-- | NFData instance +instance NFData a => NFData (Tensor a) + +-- | Move contravariant indices to lower recursion level +{-# INLINE _standardize #-} +_standardize :: (Unboxed.Unbox a) => Tensor a -> Tensor a +_standardize tens = foldl' (<<<|) tens $ Index.indexName <$> (Index.isContravariant `Prelude.filter` indices tens) + +-- | Print tensor +instance ( + Unboxed.Unbox a, Show a + ) => Show (Tensor a) where + + -- merge errors first and then print whole tensor + show = show' . _standardize + where + show' x = case x of + -- Scalar is showed simply as its value + Scalar v -> show v + -- SimpleFinite is shown dependent on its index... + SimpleFinite index ts -> show index ++ "S: " ++ case index of + -- If index is contravariant, show tensor components vertically + Finite.Contravariant _ _ -> "\n" ++ tail (Unboxed.foldl' (\string e -> string ++ "\n |" ++ show e) "" ts) + -- If index is covariant or indifferent, show tensor compoments horizontally + _ -> "[" ++ tail (Unboxed.foldl' (\string e -> string ++ "," ++ show e) "" ts) ++ "]" + -- FiniteTensor is shown dependent on its index... + FiniteTensor index ts -> show index ++ "T: " ++ case index of + -- If index is contravariant, show tensor components vertically + Finite.Contravariant _ _ -> "\n" ++ tail (Boxed.foldl' (\string e -> string ++ "\n |" ++ show e) "" ts) + -- If index is covariant or indifferent, show tensor compoments horizontally + _ -> "[" ++ tail (Boxed.foldl' (\string e -> string ++ "," ++ show e) "" ts) ++ "]" + +{-| Merge FiniteTensor of Scalars to SimpleFinite tensor for performance improvement -} +{-# INLINE _mergeScalars #-} +_mergeScalars :: Unboxed.Unbox a => Tensor a -> Tensor a +_mergeScalars x = case x of + (FiniteTensor index1 ts1) -> case ts1 Boxed.! 0 of + Scalar _ -> SimpleFinite index1 $ Unboxed.generate (Boxed.length ts1) (\i -> scalarVal (ts1 Boxed.! i)) + _ -> FiniteTensor index1 $ _mergeScalars <$> ts1 + _ -> x + +-- | Generic map function, which does not require a,b types to be Num +_map :: ( + Unboxed.Unbox a, Unboxed.Unbox b, NFData b + ) => (a -> b) + -> Tensor a + -> Tensor b +_map f x = case x of + -- Mapping scalar simply maps its value + Scalar v -> Scalar $ f v + -- Mapping complex tensor does mapping element by element + SimpleFinite index ts -> SimpleFinite index (f `Unboxed.map` ts) + FiniteTensor index ts -> + let len = Boxed.length ts + lts = Boxed.toList $ _map f <$> ts + ltsp = lts `Parallel.using` Parallel.parListChunk (len `div` 8) Parallel.rdeepseq + in FiniteTensor index $ Boxed.fromList ltsp + +{-| Transpose Vector of Vectors, analogous to Data.List.transpose function. It is assumed, that all vectors on deeper recursion level have the same length. -} +_transpose :: Boxed.Vector (Boxed.Vector a) -- ^ Vector of vectors to transpose + -> Boxed.Vector (Boxed.Vector a) +_transpose v = + let outerS = Boxed.length v + innerS = Boxed.length $ v Boxed.! 0 + in Boxed.generate innerS (\i -> Boxed.generate outerS $ \j -> v Boxed.! j Boxed.! i) + +-- | Contracted indices have to be consumed in result tensor. +_contractedIndices :: + Tensor Double -- ^ first tensor to contract + -> Tensor Double -- ^ second tensor to contract + -> Set.Set String +_contractedIndices t1 t2 = + let iContravariantNames1 = Set.fromList $ Index.indexName <$> (Index.isContravariant `Prelude.filter` indices t1) + iCovariantNames1 = Set.fromList $ Index.indexName <$> (Index.isCovariant `Prelude.filter` indices t1) + iContravariantNames2 = Set.fromList $ Index.indexName <$> (Index.isContravariant `Prelude.filter` indices t2) + iCovariantNames2 = Set.fromList $ Index.indexName <$> (Index.isCovariant `Prelude.filter` indices t2) + in -- contracted are indices covariant in the first tensor and contravariant in the second + Set.intersection iCovariantNames1 iContravariantNames2 `Set.union` + -- or contravariant in the first tensor and covariant in the second + Set.intersection iContravariantNames1 iCovariantNames2 + +{-| Apply a tensor operator (here denoted by (+) ) elem by elem, trying to connect as many common indices as possible -} +{-# INLINE _elemByElem' #-} +_elemByElem' :: (Num a, Unboxed.Unbox a, NFData a) + => Tensor a -- ^ First argument of operator + -> Tensor a -- ^ Second argument of operator + -> (a -> a -> a) -- ^ Operator on tensor elements if indices are different + -> (Tensor a -> Tensor a -> Tensor a) -- ^ Tensor operator called if indices are the same + -> Tensor a -- ^ Result tensor +-- @Scalar x + Scalar y = Scalar x + y@ +_elemByElem' (Scalar x1) (Scalar x2) f _ = Scalar $ f x1 x2 +-- @Scalar x + Tensor t[i] = Tensor r[i] | r[i] = x + t[i]@ +_elemByElem' (Scalar x) t f _ = (x `f`) `Multilinear.Generic.MultiCore.map` t +-- @Tensor t[i] + Scalar x = Tensor r[i] | r[i] = t[i] + x@ +_elemByElem' t (Scalar x) f _ = (`f` x) `Multilinear.Generic.MultiCore.map` t +-- Two simple tensors case +_elemByElem' t1@(SimpleFinite index1 v1) t2@(SimpleFinite index2 _) f op + | Index.indexName index1 == Index.indexName index2 = op t1 t2 + | otherwise = FiniteTensor index1 $ Boxed.generate (Unboxed.length v1) + (\i -> (\x -> f x `Multilinear.Generic.MultiCore.map` t2) (v1 Unboxed.! i)) +-- Two finite tensors case +_elemByElem' t1@(FiniteTensor index1 v1) t2@(FiniteTensor index2 v2) f op + | Index.indexName index1 == Index.indexName index2 = op t1 t2 + | Index.indexName index1 `Data.List.elem` indicesNames t2 = + let len2 = Finite.indexSize index2 + rl = Boxed.toList $ (\x -> _elemByElem' t1 x f op) <$> v2 + rlp = rl `Parallel.using` Parallel.parListChunk (len2 `div` 8) Parallel.rdeepseq + in FiniteTensor index2 $ Boxed.fromList rlp + | otherwise = + let len1 = Finite.indexSize index1 + rl = Boxed.toList $ (\x -> _elemByElem' x t2 f op) <$> v1 + rlp = rl `Parallel.using` Parallel.parListChunk (len1 `div` 8) Parallel.rdeepseq + in FiniteTensor index1 $ Boxed.fromList rlp +-- Simple and finite tensor case +_elemByElem' t1@(SimpleFinite index1 _) t2@(FiniteTensor index2 v2) f op + | Index.indexName index1 == Index.indexName index2 = op t1 t2 + | otherwise = + let len2 = Finite.indexSize index2 + rl = Boxed.toList $ (\x -> _elemByElem' t1 x f op) <$> v2 + rlp = rl `Parallel.using` Parallel.parListChunk (len2 `div` 8) Parallel.rdeepseq + in FiniteTensor index2 $ Boxed.fromList rlp +-- Finite and simple tensor case +_elemByElem' t1@(FiniteTensor index1 v1) t2@(SimpleFinite index2 _) f op + | Index.indexName index1 == Index.indexName index2 = op t1 t2 + | otherwise = + let len1 = Finite.indexSize index1 + rl = Boxed.toList $ (\x -> _elemByElem' x t2 f op) <$> v1 + rlp = rl `Parallel.using` Parallel.parListChunk (len1 `div` 8) Parallel.rdeepseq + in FiniteTensor index1 $ Boxed.fromList rlp + +{-| Apply a tensor operator elem by elem and merge scalars to simple tensor at the and -} +{-# INLINE _elemByElem #-} +_elemByElem :: (Num a, Unboxed.Unbox a, NFData a) + => Tensor a -- ^ First argument of operator + -> Tensor a -- ^ Second argument of operator + -> (a -> a -> a) -- ^ Operator on tensor elements if indices are different + -> (Tensor a -> Tensor a -> Tensor a) -- ^ Tensor operator called if indices are the same + -> Tensor a -- ^ Result tensor +_elemByElem t1 t2 f op = + let commonIndices = + if indices t1 /= indices t2 then + Data.List.filter (`Data.List.elem` indicesNames t2) $ indicesNames t1 + else [] + t1' = foldl' (|>>>) t1 commonIndices + t2' = foldl' (|>>>) t2 commonIndices + in _mergeScalars $ _elemByElem' t1' t2' f op + +-- | Zipping two tensors with a combinator, assuming they have the same indices. +{-# INLINE zipT #-} +zipT :: ( + Num a, NFData a, Unboxed.Unbox a + ) => (a -> a -> a) -- ^ The zipping combinator + -> Tensor a -- ^ First tensor to zip + -> Tensor a -- ^ Second tensor to zip + -> Tensor a -- ^ Result tensor +-- Two simple tensors case +zipT f t1@(SimpleFinite index1 v1) t2@(SimpleFinite index2 v2) = + if index1 == index2 then + SimpleFinite index1 $ Unboxed.zipWith f v1 v2 + else dot t1 t2 +--Two finite tensors case +zipT f t1@(FiniteTensor index1 v1) t2@(FiniteTensor index2 v2) = + if index1 == index2 then let + l1l = Boxed.length v1 + l3 = Boxed.toList (Boxed.zipWith (zipT f) v1 v2) `Parallel.using` Parallel.parListChunk (l1l `div` 8) Parallel.rdeepseq + in FiniteTensor index1 $ Boxed.fromList l3 + else dot t1 t2 +-- Zipping something with scalar is impossible +zipT _ _ _ = error $ "zipT: " ++ scalarIndices + +-- | dot product of two tensors +{-# INLINE dot #-} +dot :: (Num a, Unboxed.Unbox a, NFData a) + => Tensor a -- ^ First dot product argument + -> Tensor a -- ^ Second dot product argument + -> Tensor a -- ^ Resulting dot product +-- Two simple tensors product +dot (SimpleFinite i1@(Finite.Covariant count1 _) ts1') (SimpleFinite i2@(Finite.Contravariant count2 _) ts2') + | count1 == count2 = + Scalar $ Unboxed.sum $ Unboxed.zipWith (*) ts1' ts2' + | otherwise = contractionErr "simple-simple" (Index.toTIndex i1) (Index.toTIndex i2) +dot (SimpleFinite i1@(Finite.Contravariant count1 _) ts1') (SimpleFinite i2@(Finite.Covariant count2 _) ts2') + | count1 == count2 = + Scalar $ Unboxed.sum $ Unboxed.zipWith (*) ts1' ts2' + | otherwise = contractionErr "simple-simple" (Index.toTIndex i1) (Index.toTIndex i2) +dot t1@(SimpleFinite _ _) t2@(SimpleFinite _ _) = zipT (*) t1 t2 +-- Two finite tensors product +dot (FiniteTensor i1@(Finite.Covariant count1 _) ts1') (FiniteTensor i2@(Finite.Contravariant count2 _) ts2') + | count1 == count2 = + let zipList = Boxed.toList $ Boxed.zipWith (*) ts1' ts2' + zipListPar = zipList `Parallel.using` Parallel.parListChunk (count1 `div` 8) Parallel.rdeepseq + in Boxed.sum $ Boxed.fromList zipListPar + | otherwise = contractionErr "finite-finite" (Index.toTIndex i1) (Index.toTIndex i2) +dot (FiniteTensor i1@(Finite.Contravariant count1 _) ts1') (FiniteTensor i2@(Finite.Covariant count2 _) ts2') + | count1 == count2 = + let zipList = Boxed.toList $ Boxed.zipWith (*) ts1' ts2' + zipListPar = zipList `Parallel.using` Parallel.parListChunk (count1 `div` 8) Parallel.rdeepseq + in Boxed.sum $ Boxed.fromList zipListPar + | otherwise = contractionErr "finite-finite" (Index.toTIndex i1) (Index.toTIndex i2) +dot t1@(FiniteTensor _ _) t2@(FiniteTensor _ _) = zipT (*) t1 t2 +-- Other cases cannot happen! +dot t1' t2' = contractionErr "other" (tensorIndex t1') (tensorIndex t2') + +-- | contraction error +{-# INLINE contractionErr #-} +contractionErr :: String -- ^ dot function variant where error occured + -> Index.TIndex -- ^ Index of first dot product parameter + -> Index.TIndex -- ^ Index of second dot product parameter + -> Tensor a -- ^ Erorr message + +contractionErr variant i1' i2' = error $ + "Tensor product: " ++ variant ++ " - " ++ incompatibleTypes ++ + " - index1 is " ++ show i1' ++ + " and index2 is " ++ show i2' + +-- | zipping error +{-# INLINE zipErr #-} +zipErr :: String -- ^ zipT function variant where the error occured + -> Index.TIndex -- ^ Index of first dot product parameter + -> Index.TIndex -- ^ Index of second dot product parameter + -> Tensor a -- ^ Erorr message +zipErr variant i1' i2' = error $ + "zipT: " ++ variant ++ " - " ++ incompatibleTypes ++ + " - index1 is " ++ show i1' ++ + " and index2 is " ++ show i2' + + +-- | Tensors can be added, subtracted and multiplicated +instance (Unboxed.Unbox a, Num a, NFData a) => Num (Tensor a) where + + -- Adding - element by element + {-# INLINE (+) #-} + t1 + t2 = _elemByElem t1 t2 (+) $ zipT (+) + + -- Subtracting - element by element + {-# INLINE (-) #-} + t1 - t2 = _elemByElem t1 t2 (-) $ zipT (-) + + -- Multiplicating is treated as tensor product + -- Tensor product applies Einstein summation convention + {-# INLINE (*) #-} + t1 * t2 = _elemByElem t1 t2 (*) dot + + -- Absolute value - element by element + {-# INLINE abs #-} + abs t = abs `Multilinear.Generic.MultiCore.map` t + + -- Signum operation - element by element + {-# INLINE signum #-} + signum t = signum `Multilinear.Generic.MultiCore.map` t + + -- Simple integer can be conveted to Scalar + {-# INLINE fromInteger #-} + fromInteger x = Scalar $ fromInteger x + +-- | Tensors can be divided by each other +instance (Unboxed.Unbox a, Fractional a, NFData a) => Fractional (Tensor a) where + -- Tensor dividing: TODO + {-# INLINE (/) #-} + _ / _ = error "TODO" + + -- A scalar can be generated from rational number + {-# INLINE fromRational #-} + fromRational x = Scalar $ fromRational x + +-- Real-number functions on tensors. +-- Function of tensor is tensor of function of its elements +-- E.g. exp [1,2,3,4] = [exp 1, exp2, exp3, exp4] +instance (Unboxed.Unbox a, Floating a, NFData a) => Floating (Tensor a) where + + {-| PI number -} + {-# INLINE pi #-} + pi = Scalar pi + + {-| Exponential function. (exp t)[i] = exp( t[i] ) -} + {-# INLINE exp #-} + exp t = exp `Multilinear.Generic.MultiCore.map` t + + {-| Natural logarithm. (log t)[i] = log( t[i] ) -} + {-# INLINE log #-} + log t = log `Multilinear.Generic.MultiCore.map` t + + {-| Sinus. (sin t)[i] = sin( t[i] ) -} + {-# INLINE sin #-} + sin t = sin `Multilinear.Generic.MultiCore.map` t + + {-| Cosinus. (cos t)[i] = cos( t[i] ) -} + {-# INLINE cos #-} + cos t = cos `Multilinear.Generic.MultiCore.map` t + + {-| Inverse sinus. (asin t)[i] = asin( t[i] ) -} + {-# INLINE asin #-} + asin t = asin `Multilinear.Generic.MultiCore.map` t + + {-| Inverse cosinus. (acos t)[i] = acos( t[i] ) -} + {-# INLINE acos #-} + acos t = acos `Multilinear.Generic.MultiCore.map` t + + {-| Inverse tangent. (atan t)[i] = atan( t[i] ) -} + {-# INLINE atan #-} + atan t = atan `Multilinear.Generic.MultiCore.map` t + + {-| Hyperbolic sinus. (sinh t)[i] = sinh( t[i] ) -} + {-# INLINE sinh #-} + sinh t = sinh `Multilinear.Generic.MultiCore.map` t + + {-| Hyperbolic cosinus. (cosh t)[i] = cosh( t[i] ) -} + {-# INLINE cosh #-} + cosh t = cosh `Multilinear.Generic.MultiCore.map` t + + {-| Inverse hyperbolic sinus. (asinh t)[i] = asinh( t[i] ) -} + {-# INLINE asinh #-} + asinh t = acosh `Multilinear.Generic.MultiCore.map` t + + {-| Inverse hyperbolic cosinus. (acosh t)[i] = acosh (t[i] ) -} + {-# INLINE acosh #-} + acosh t = acosh `Multilinear.Generic.MultiCore.map` t + + {-| Inverse hyperbolic tangent. (atanh t)[i] = atanh( t[i] ) -} + {-# INLINE atanh #-} + atanh t = atanh `Multilinear.Generic.MultiCore.map` t + +-- Multilinear operations +instance (Unboxed.Unbox a) => Multilinear Tensor a where + -- Generic tensor constructor + -- If only one upper index is given, generate a SimpleFinite tensor with upper index + fromIndices [u] [] [s] [] f = + SimpleFinite (Finite.Contravariant s [u]) $ Unboxed.generate s $ \x -> f [x] [] + + -- If only one lower index is given, generate a SimpleFinite tensor with lower index + fromIndices [] [d] [] [s] f = + SimpleFinite (Finite.Covariant s [d]) $ Unboxed.generate s $ \x -> f [] [x] + + -- If many indices are given, first generate upper indices recursively from indices list + fromIndices (u:us) d (s:size) dsize f = + FiniteTensor (Finite.Contravariant s [u]) $ Boxed.generate s (\x -> fromIndices us d size dsize (\uss dss -> f (x:uss) dss) ) + + -- After upper indices, generate lower indices recursively from indices list + fromIndices u (d:ds) usize (s:size) f = + FiniteTensor (Finite.Covariant s [d]) $ Boxed.generate s (\x -> fromIndices u ds usize size (\uss dss -> f uss (x:dss)) ) + + -- If there are indices without size or sizes without names, throw an error + fromIndices _ _ _ _ _ = error "Indices and its sizes incompatible!" + + {-| Accessing tensor elements -} + {-# INLINE el #-} + -- Scalar has only one element + el (Scalar x) _ = Scalar x + -- simple tensor case + el t1@(SimpleFinite index1 _) (inds,vals) = + -- zip indices with their given values + let indvals = zip inds vals + -- find value for simple tensor index if given + val = Data.List.find (\(n,_) -> [n] == Index.indexName index1) indvals + -- if value for current index is given + in if isJust val + -- then get it from current tensor + then t1 ! snd (fromJust val) + -- otherwise return whole tensor - no filtering defined + else t1 + -- finite tensor case + el t1@(FiniteTensor index1 v1) (inds,vals) = + -- zip indices with their given values + let indvals = zip inds vals + -- find value for current index if given + val = Data.List.find (\(n,_) -> [n] == Index.indexName index1) indvals + -- and remove used index from indices list + indvals1 = Data.List.filter (\(n,_) -> [n] /= Index.indexName index1) indvals + -- indices unused so far + inds1 = Data.List.map fst indvals1 + -- and its corresponding values + vals1 = Data.List.map snd indvals1 + -- if value for current index was given + in if isJust val + -- then get it from current tensor and recursively process other indices + then el (t1 ! snd (fromJust val)) (inds1,vals1) + -- otherwise recursively access elements of all child tensors + else FiniteTensor index1 $ (\t -> el t (inds,vals)) <$> v1 + + -- List of all tensor indices + {-# INLINE indices #-} + indices x = case x of + Scalar _ -> [] + FiniteTensor i ts -> Index.toTIndex i : indices (head $ toList ts) + SimpleFinite i _ -> [Index.toTIndex i] + + -- Get tensor order [ (contravariant,covariant)-type ] + {-# INLINE order #-} + order x = case x of + Scalar _ -> (0,0) + SimpleFinite index _ -> case index of + Finite.Contravariant _ _ -> (1,0) + Finite.Covariant _ _ -> (0,1) + _ -> let (cnvr, covr) = order $ firstTensor x + in case tensorIndex x of + Index.Contravariant _ _ -> (cnvr+1,covr) + Index.Covariant _ _ -> (cnvr,covr+1) + + -- Get size of tensor index or Left if index is infinite or tensor has no such index + {-# INLINE size #-} + size t iname = case t of + Scalar _ -> error scalarIndices + SimpleFinite index _ -> + if Index.indexName index == iname + then Finite.indexSize index + else error indexNotFound + FiniteTensor index _ -> + if Index.indexName index == iname + then Finite.indexSize index + else size (firstTensor t) iname + + -- Rename tensor indices + {-# INLINE ($|) #-} + + Scalar x $| _ = Scalar x + SimpleFinite (Finite.Contravariant isize _) ts $| (u:_, _) = SimpleFinite (Finite.Contravariant isize [u]) ts + SimpleFinite (Finite.Covariant isize _) ts $| (_, d:_) = SimpleFinite (Finite.Covariant isize [d]) ts + FiniteTensor (Finite.Contravariant isize _) ts $| (u:us, ds) = FiniteTensor (Finite.Contravariant isize [u]) $ ($| (us,ds)) <$> ts + FiniteTensor (Finite.Covariant isize _) ts $| (us, d:ds) = FiniteTensor (Finite.Covariant isize [d]) $ ($| (us,ds)) <$> ts + t $| _ = t + + -- Raise an index + {-# INLINE (/\) #-} + Scalar x /\ _ = Scalar x + FiniteTensor index ts /\ n + | Index.indexName index == n = + FiniteTensor (Finite.Contravariant (Finite.indexSize index) n) $ (/\ n) <$> ts + | otherwise = + FiniteTensor index $ (/\ n) <$> ts + t1@(SimpleFinite index ts) /\ n + | Index.indexName index == n = + SimpleFinite (Finite.Contravariant (Finite.indexSize index) n) ts + | otherwise = t1 + + -- Lower an index + {-# INLINE (\/) #-} + Scalar x \/ _ = Scalar x + FiniteTensor index ts \/ n + | Index.indexName index == n = + FiniteTensor (Finite.Covariant (Finite.indexSize index) n) $ (\/ n) <$> ts + | otherwise = + FiniteTensor index $ (\/ n) <$> ts + t1@(SimpleFinite index ts) \/ n + | Index.indexName index == n = + SimpleFinite (Finite.Covariant (Finite.indexSize index) n) ts + | otherwise = t1 + + {-| Transpose a tensor (switch all indices types) -} + {-# INLINE transpose #-} + transpose (Scalar x) = Scalar x + + transpose (FiniteTensor (Finite.Covariant count name) ts) = + FiniteTensor (Finite.Contravariant count name) (Multilinear.transpose <$> ts) + transpose (FiniteTensor (Finite.Contravariant count name) ts) = + FiniteTensor (Finite.Covariant count name) (Multilinear.transpose <$> ts) + + transpose (SimpleFinite (Finite.Covariant count name) ts) = + SimpleFinite (Finite.Contravariant count name) ts + transpose (SimpleFinite (Finite.Contravariant count name) ts) = + SimpleFinite (Finite.Covariant count name) ts + + + {-| Shift tensor index right -} + {-| Moves given index one level deeper in recursion -} + -- Scalar has no indices to shift + Scalar x `shiftRight` _ = Scalar x + -- Simple tensor has only one index which cannot be shifted + t1@(SimpleFinite _ _) `shiftRight` _ = t1 + -- Finite tensor is shifted by converting to list and transposing it + t1@(FiniteTensor index1 ts1) `shiftRight` ind + -- We don't shift this index + | Data.List.length (indicesNames t1) > 1 && Index.indexName index1 /= ind = + FiniteTensor index1 $ (|>> ind) <$> ts1 + -- We found index to shift + | Data.List.length (indicesNames t1) > 1 && Index.indexName index1 == ind = + -- Next index + let index2 = tensorFiniteIndex (ts1 Boxed.! 0) + -- Elements to transpose + dane = if isSimple (ts1 Boxed.! 0) + then (\un -> Boxed.generate (Unboxed.length un) (\i -> Scalar $ un Unboxed.! i)) <$> + (tensorScalars <$> ts1) + else tensorsFinite <$> ts1 + result = FiniteTensor index2 $ FiniteTensor index1 <$> (_transpose dane) + -- reconstruct tensor with transposed elements + in _mergeScalars result + -- there is only one index and therefore it cannot be shifted + | otherwise = t1 + +-- Add scalar right +{-# INLINE (.+) #-} +(.+) :: ( + Unboxed.Unbox a, Num a, NFData a + ) => Tensor a + -> a + -> Tensor a +t .+ x = (+x) `Multilinear.Generic.MultiCore.map` t + +-- Subtract scalar right +{-# INLINE (.-) #-} +(.-) :: ( + Unboxed.Unbox a, Num a, NFData a + ) => Tensor a + -> a + -> Tensor a +t .- x = (\p -> p - x) `Multilinear.Generic.MultiCore.map` t + +-- Multiplicate by scalar right +{-# INLINE (.*) #-} +(.*) :: ( + Unboxed.Unbox a, Num a, NFData a + ) => Tensor a + -> a + -> Tensor a +t .* x = (*x) `Multilinear.Generic.MultiCore.map` t + +-- Add scalar left +{-# INLINE (+.) #-} +(+.) :: ( + Unboxed.Unbox a, Num a, NFData a + ) => a + -> Tensor a + -> Tensor a +x +. t = (x+) `Multilinear.Generic.MultiCore.map` t + +-- Subtract scalar left +{-# INLINE (-.) #-} +(-.) :: ( + Unboxed.Unbox a, Num a, NFData a + ) => a + -> Tensor a + -> Tensor a +x -. t = (x-) `Multilinear.Generic.MultiCore.map` t + +-- Multiplicate by scalar left +{-# INLINE (*.) #-} +(*.) :: ( + Unboxed.Unbox a, Num a, NFData a + ) => a + -> Tensor a + -> Tensor a +x *. t = (x*) `Multilinear.Generic.MultiCore.map` t + + +{-| Simple mapping -} +{-| @map f t@ returns tensor @t2@ in which @t2[i1,i2,...] = f t[i1,i2,...]@ -} +{-# INLINE map #-} +map :: ( + Unboxed.Unbox a, Unboxed.Unbox b, NFData b + ) => (a -> b) + -> Tensor a + -> Tensor b +map f x = case x of + -- Mapping scalar simply maps its value + Scalar v -> Scalar $ f v + -- Mapping complex tensor does mapping element by element + SimpleFinite index ts -> SimpleFinite index (f `Unboxed.map` ts) + FiniteTensor index ts -> + let len = Boxed.length ts + lts = Boxed.toList $ Multilinear.Generic.MultiCore.map f <$> ts + ltsp = lts `Parallel.using` Parallel.parListChunk (len `div` 8) Parallel.rdeepseq + in FiniteTensor index $ Boxed.fromList ltsp + +{-| Filtering tensor. + Filtering multi-dimensional arrray may be dangerous, as we always assume, + that on each recursion level, all tensors have the same size (length). + To disable invalid filters, filtering is done over indices, not tensor elements. + Filter function takes and index name and index value and if it returns True, this index value remains in result tensor. + This allows to remove whole columns or rows of eg. a matrix: + filter (\i n -> i /= "a" || i > 10) filters all rows of "a" index (because if i /= "a", filter returns True) + and for "a" index filter elements with index value <= 10 + But this disallow to remove particular matrix element. + If for some index all elements are removed, the index itself is removed from tensor. -} +{-# INLINE filter #-} +filter :: ( + Unboxed.Unbox a + ) => (String -> Int -> Bool) -- ^ filter function + -> Tensor a -- ^ tensor to filter + -> Tensor a +filter _ (Scalar x) = Scalar x +filter f (SimpleFinite index ts) = + let iname = Finite.indexName' index + ts' = (\i _ -> f iname i) `Unboxed.ifilter` ts + in SimpleFinite index { Finite.indexSize = Unboxed.length ts' } ts' +filter f (FiniteTensor index ts) = + let iname = Finite.indexName' index + ts' = Multilinear.Generic.MultiCore.filter f <$> ((\i _ -> f iname i) `Boxed.ifilter` ts) + ts'' = + (\case + (SimpleFinite _ ts) -> not $ Unboxed.null ts + (FiniteTensor _ ts) -> not $ Boxed.null ts + _ -> error $ "Filter: " ++ tensorOfScalars + ) `Boxed.filter` ts' + in FiniteTensor index { Finite.indexSize = Boxed.length ts'' } ts'' + +{-| Filtering one index of tensor. -} +{-# INLINE filterIndex #-} +filterIndex :: ( + Unboxed.Unbox a + ) => String -- ^ Index name to filter + -> (Int -> Bool) -- ^ filter function + -> Tensor a -- ^ tensor to filter + -> Tensor a +filterIndex iname f = Multilinear.Generic.MultiCore.filter (\i n -> i /= iname || f n) + +{-| Zip tensors with binary combinator, assuming they have all indices the same -} +{-# INLINE zipWith' #-} +zipWith' :: ( + Unboxed.Unbox a, Unboxed.Unbox b, Unboxed.Unbox c, NFData c + ) => (a -> b -> c) + -> Tensor a + -> Tensor b + -> Tensor c +-- Zipping two Scalars simply combines their values +zipWith' f (Scalar x1) (Scalar x2) = Scalar $ f x1 x2 +-- zipping complex tensor with scalar +zipWith' f t (Scalar x) = (`f` x) `_map` t +-- zipping scalar with complex tensor +zipWith' f (Scalar x) t = (x `f`) `_map` t +-- Two simple tensors case +zipWith' f (SimpleFinite index1 v1) (SimpleFinite index2 v2) = + if index1 == index2 then + SimpleFinite index1 $ Unboxed.zipWith f v1 v2 + else zipErr "simple-simple" (Index.toTIndex index1) (Index.toTIndex index2) +--Two finite tensors case +zipWith' f (FiniteTensor index1 v1) (FiniteTensor index2 v2) = + if index1 == index2 then + FiniteTensor index1 $ Boxed.zipWith (Multilinear.Generic.MultiCore.zipWith f) v1 v2 + else zipErr "finite-finite" (Index.toTIndex index1) (Index.toTIndex index2) +-- Other cases cannot happen! +zipWith' _ _ _ = error "Invalid indices to peroform zip!" + +{-# INLINE zipWith #-} +zipWith :: ( + Unboxed.Unbox a, Unboxed.Unbox b, Unboxed.Unbox c, NFData c + ) => (a -> b -> c) + -> Tensor a + -> Tensor b + -> Tensor c +zipWith f t1 t2 = + let t1' = _standardize t1 + t2' = _standardize t2 + in zipWith' f t1' t2' diff --git a/src/Multilinear/Generic/Sequential.hs b/src/Multilinear/Generic/Sequential.hs new file mode 100644 --- /dev/null +++ b/src/Multilinear/Generic/Sequential.hs @@ -0,0 +1,736 @@ +{-| +Module : Multilinear.Generic.Sequential +Description : Generic implementation of tensor as nested arrays, evaluated in sequential manner +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Multilinear.Generic.Sequential ( + -- * Generic tensor datatype and its instances + Tensor(..), + -- * Auxiliary functions + (!), isScalar, isSimple, isFiniteTensor, + tensorIndex, _standardize, _mergeScalars, + _map, _contractedIndices, _elemByElem, zipT, + -- * Additional functions + (.+), (.-), (.*), (+.), (-.), (*.), + Multilinear.Generic.Sequential.map, + Multilinear.Generic.Sequential.filter, + Multilinear.Generic.Sequential.filterIndex, + Multilinear.Generic.Sequential.zipWith +) where + +import Control.DeepSeq +import Data.Foldable +import Data.List +import Data.Maybe +import qualified Data.Set as Set +import qualified Data.Vector as Boxed +import qualified Data.Vector.Unboxed as Unboxed +import GHC.Generics +import Multilinear.Class as Multilinear +import qualified Multilinear.Index as Index +import qualified Multilinear.Index.Finite as Finite + +{-| ERROR MESSAGE -} +incompatibleTypes :: String +incompatibleTypes = "Incompatible tensor types!" + +{-| ERROR MESSAGE -} +scalarIndices :: String +scalarIndices = "Scalar has no indices!" + +{-| ERROR MESSAGE -} +indexNotFound :: String +indexNotFound = "This tensor has not such index!" + +{-| ERROR MESSAGE -} +tensorOfScalars :: String +tensorOfScalars = "Tensor construction error! Vector of scalars" + +{-| Tensor defined recursively as scalar or list of other tensors -} +{-| @c@ is type of a container, @i@ is type of index size and @a@ is type of tensor elements -} +data Tensor a where + {-| Scalar -} + Scalar :: { + scalarVal :: a + } -> Tensor a + {-| Simple, one-dimensional finite tensor -} + SimpleFinite :: { + tensorFiniteIndex :: Finite.Index, + tensorScalars :: Unboxed.Vector a + } -> Tensor a + {-| Finite array of other tensors -} + FiniteTensor :: { + {-| Finite index "Mutltilinear.Index.Finite" of tensor -} + tensorFiniteIndex :: Finite.Index, + {-| Array of tensors on deeper recursion level -} + tensorsFinite :: Boxed.Vector (Tensor a) + } -> Tensor a + deriving (Eq, Generic) + +{-| Return true if tensor is a scalar -} +{-# INLINE isScalar #-} +isScalar :: Unboxed.Unbox a => Tensor a -> Bool +isScalar x = case x of + Scalar _ -> True + _ -> False + +{-| Return true if tensor is a simple tensor -} +{-# INLINE isSimple #-} +isSimple :: Unboxed.Unbox a => Tensor a -> Bool +isSimple x = case x of + SimpleFinite _ _ -> True + _ -> False + +{-| Return True if tensor is a complex tensor -} +{-# INLINE isFiniteTensor #-} +isFiniteTensor :: Unboxed.Unbox a => Tensor a -> Bool +isFiniteTensor x = case x of + FiniteTensor _ _ -> True + _ -> False + +{-| Return generic tensor index -} +{-# INLINE tensorIndex #-} +tensorIndex :: Unboxed.Unbox a => Tensor a -> Index.TIndex +tensorIndex x = case x of + Scalar _ -> error scalarIndices + SimpleFinite i _ -> Index.toTIndex i + FiniteTensor i _ -> Index.toTIndex i + +{-| Returns sample tensor on deeper recursion level.Used to determine some features common for all tensors -} +{-# INLINE firstTensor #-} +firstTensor :: Unboxed.Unbox a => Tensor a -> Tensor a +firstTensor x = case x of + FiniteTensor _ ts -> Boxed.head ts + _ -> x + +{-| Recursive indexing on list tensor. If index is greater than index size, performs modulo indexing + @t ! i = t[i]@ -} +{-# INLINE (!) #-} +(!) :: Unboxed.Unbox a => Tensor a -- ^ tensor @t@ + -> Int -- ^ index @i@ + -> Tensor a -- ^ tensor @t[i]@ +t ! i = case t of + Scalar _ -> error scalarIndices + SimpleFinite ind ts -> Scalar $ ts Unboxed.! (i `mod` Finite.indexSize ind) + FiniteTensor ind ts -> ts Boxed.! (i `mod` Finite.indexSize ind) + +-- | NFData instance +instance NFData a => NFData (Tensor a) + +-- | Move contravariant indices to lower recursion level +{-# INLINE _standardize #-} +_standardize :: (Unboxed.Unbox a) => Tensor a -> Tensor a +_standardize tens = foldl' (<<<|) tens $ Index.indexName <$> (Index.isContravariant `Prelude.filter` indices tens) + +-- | Print tensor +instance ( + Unboxed.Unbox a, Show a + ) => Show (Tensor a) where + + -- merge errors first and then print whole tensor + show = show' . _standardize + where + show' x = case x of + -- Scalar is showed simply as its value + Scalar v -> show v + -- SimpleFinite is shown dependent on its index... + SimpleFinite index ts -> show index ++ "S: " ++ case index of + -- If index is contravariant, show tensor components vertically + Finite.Contravariant _ _ -> "\n" ++ tail (Unboxed.foldl' (\string e -> string ++ "\n |" ++ show e) "" ts) + -- If index is covariant or indifferent, show tensor compoments horizontally + _ -> "[" ++ tail (Unboxed.foldl' (\string e -> string ++ "," ++ show e) "" ts) ++ "]" + -- FiniteTensor is shown dependent on its index... + FiniteTensor index ts -> show index ++ "T: " ++ case index of + -- If index is contravariant, show tensor components vertically + Finite.Contravariant _ _ -> "\n" ++ tail (Boxed.foldl' (\string e -> string ++ "\n |" ++ show e) "" ts) + -- If index is covariant or indifferent, show tensor compoments horizontally + _ -> "[" ++ tail (Boxed.foldl' (\string e -> string ++ "," ++ show e) "" ts) ++ "]" + +{-| Merge FiniteTensor of Scalars to SimpleFinite tensor for performance improvement -} +{-# INLINE _mergeScalars #-} +_mergeScalars :: Unboxed.Unbox a => Tensor a -> Tensor a +_mergeScalars x = case x of + (FiniteTensor index1 ts1) -> case ts1 Boxed.! 0 of + Scalar _ -> SimpleFinite index1 $ Unboxed.generate (Boxed.length ts1) (\i -> scalarVal (ts1 Boxed.! i)) + _ -> FiniteTensor index1 $ _mergeScalars <$> ts1 + _ -> x + +-- | Generic map function, which does not require a,b types to be Num +_map :: ( + Unboxed.Unbox a, Unboxed.Unbox b + ) => (a -> b) + -> Tensor a + -> Tensor b +_map f x = case x of + -- Mapping scalar simply maps its value + Scalar v -> Scalar $ f v + -- Mapping complex tensor does mapping element by element + SimpleFinite index ts -> SimpleFinite index (f `Unboxed.map` ts) + FiniteTensor index ts -> FiniteTensor index $ _map f <$> ts + +{-| Transpose Vector of Vectors, analogous to Data.List.transpose function. It is assumed, that all vectors on deeper recursion level have the same length. -} +_transpose :: Boxed.Vector (Boxed.Vector a) -- ^ Vector of vectors to transpose + -> Boxed.Vector (Boxed.Vector a) +_transpose v = + let outerS = Boxed.length v + innerS = Boxed.length $ v Boxed.! 0 + in Boxed.generate innerS (\i -> Boxed.generate outerS $ \j -> v Boxed.! j Boxed.! i) + +-- | Contracted indices have to be consumed in result tensor. +_contractedIndices :: + Tensor Double -- ^ first tensor to contract + -> Tensor Double -- ^ second tensor to contract + -> Set.Set String +_contractedIndices t1 t2 = + let iContravariantNames1 = Set.fromList $ Index.indexName <$> (Index.isContravariant `Prelude.filter` indices t1) + iCovariantNames1 = Set.fromList $ Index.indexName <$> (Index.isCovariant `Prelude.filter` indices t1) + iContravariantNames2 = Set.fromList $ Index.indexName <$> (Index.isContravariant `Prelude.filter` indices t2) + iCovariantNames2 = Set.fromList $ Index.indexName <$> (Index.isCovariant `Prelude.filter` indices t2) + in -- contracted are indices covariant in the first tensor and contravariant in the second + Set.intersection iCovariantNames1 iContravariantNames2 `Set.union` + -- or contravariant in the first tensor and covariant in the second + Set.intersection iContravariantNames1 iCovariantNames2 + +{-| Apply a tensor operator (here denoted by (+) ) elem by elem, trying to connect as many common indices as possible -} +{-# INLINE _elemByElem' #-} +_elemByElem' :: (Num a, Unboxed.Unbox a, NFData a) + => Tensor a -- ^ First argument of operator + -> Tensor a -- ^ Second argument of operator + -> (a -> a -> a) -- ^ Operator on tensor elements if indices are different + -> (Tensor a -> Tensor a -> Tensor a) -- ^ Tensor operator called if indices are the same + -> Tensor a -- ^ Result tensor +-- @Scalar x + Scalar y = Scalar x + y@ +_elemByElem' (Scalar x1) (Scalar x2) f _ = Scalar $ f x1 x2 +-- @Scalar x + Tensor t[i] = Tensor r[i] | r[i] = x + t[i]@ +_elemByElem' (Scalar x) t f _ = (x `f`) `Multilinear.Generic.Sequential.map` t +-- @Tensor t[i] + Scalar x = Tensor r[i] | r[i] = t[i] + x@ +_elemByElem' t (Scalar x) f _ = (`f` x) `Multilinear.Generic.Sequential.map` t +-- Two simple tensors case +_elemByElem' t1@(SimpleFinite index1 v1) t2@(SimpleFinite index2 _) f op + | Index.indexName index1 == Index.indexName index2 = op t1 t2 + | otherwise = FiniteTensor index1 $ Boxed.generate (Unboxed.length v1) + (\i -> (\x -> f x `Multilinear.Generic.Sequential.map` t2) (v1 Unboxed.! i)) +-- Two finite tensors case +_elemByElem' t1@(FiniteTensor index1 v1) t2@(FiniteTensor index2 v2) f op + | Index.indexName index1 == Index.indexName index2 = op t1 t2 + | Index.indexName index1 `Data.List.elem` indicesNames t2 = + FiniteTensor index2 $ (\x -> _elemByElem' t1 x f op) <$> v2 + | otherwise = FiniteTensor index1 $ (\x -> _elemByElem' x t2 f op) <$> v1 +-- Simple and finite tensor case +_elemByElem' t1@(SimpleFinite index1 _) t2@(FiniteTensor index2 v2) f op + | Index.indexName index1 == Index.indexName index2 = op t1 t2 + | otherwise = FiniteTensor index2 $ (\x -> _elemByElem' t1 x f op) <$> v2 +-- Finite and simple tensor case +_elemByElem' t1@(FiniteTensor index1 v1) t2@(SimpleFinite index2 _) f op + | Index.indexName index1 == Index.indexName index2 = op t1 t2 + | otherwise = FiniteTensor index1 $ (\x -> _elemByElem' x t2 f op) <$> v1 + +{-| Apply a tensor operator elem by elem and merge scalars to simple tensor at the and -} +{-# INLINE _elemByElem #-} +_elemByElem :: (Num a, Unboxed.Unbox a, NFData a) + => Tensor a -- ^ First argument of operator + -> Tensor a -- ^ Second argument of operator + -> (a -> a -> a) -- ^ Operator on tensor elements if indices are different + -> (Tensor a -> Tensor a -> Tensor a) -- ^ Tensor operator called if indices are the same + -> Tensor a -- ^ Result tensor +_elemByElem t1 t2 f op = + let commonIndices = + if indices t1 /= indices t2 then + Data.List.filter (`Data.List.elem` indicesNames t2) $ indicesNames t1 + else [] + t1' = foldl' (|>>>) t1 commonIndices + t2' = foldl' (|>>>) t2 commonIndices + in _mergeScalars $ _elemByElem' t1' t2' f op + +-- | Zipping two tensors with a combinator, assuming they have the same indices. +{-# INLINE zipT #-} +zipT :: ( + Num a, NFData a, Unboxed.Unbox a + ) => (a -> a -> a) -- ^ The zipping combinator + -> Tensor a -- ^ First tensor to zip + -> Tensor a -- ^ Second tensor to zip + -> Tensor a -- ^ Result tensor +-- Two simple tensors case +zipT f t1@(SimpleFinite index1 v1) t2@(SimpleFinite index2 v2) = + if index1 == index2 then + SimpleFinite index1 $ Unboxed.zipWith f v1 v2 + else dot t1 t2 +--Two finite tensors case +zipT f t1@(FiniteTensor index1 v1) t2@(FiniteTensor index2 v2) = + if index1 == index2 then + FiniteTensor index1 $ Boxed.zipWith (zipT f) v1 v2 + else dot t1 t2 +-- Zipping something with scalar is impossible +zipT _ _ _ = error $ "zipT: " ++ scalarIndices + +-- | dot product of two tensors +{-# INLINE dot #-} +dot :: (Num a, Unboxed.Unbox a, NFData a) + => Tensor a -- ^ First dot product argument + -> Tensor a -- ^ Second dot product argument + -> Tensor a -- ^ Resulting dot product + +-- Two simple tensors product +dot (SimpleFinite i1@(Finite.Covariant count1 _) ts1') (SimpleFinite i2@(Finite.Contravariant count2 _) ts2') + | count1 == count2 = + Scalar $ Unboxed.sum $ Unboxed.zipWith (*) ts1' ts2' + | otherwise = contractionErr "simple-simple" (Index.toTIndex i1) (Index.toTIndex i2) +dot (SimpleFinite i1@(Finite.Contravariant count1 _) ts1') (SimpleFinite i2@(Finite.Covariant count2 _) ts2') + | count1 == count2 = + Scalar $ Unboxed.sum $ Unboxed.zipWith (*) ts1' ts2' + | otherwise = contractionErr "simple-simple" (Index.toTIndex i1) (Index.toTIndex i2) +dot t1@(SimpleFinite _ _) t2@(SimpleFinite _ _) = zipT (*) t1 t2 + +-- Two finite tensors product +dot (FiniteTensor i1@(Finite.Covariant count1 _) ts1') (FiniteTensor i2@(Finite.Contravariant count2 _) ts2') + | count1 == count2 = Boxed.sum $ Boxed.zipWith (*) ts1' ts2' + | otherwise = contractionErr "finite-finite" (Index.toTIndex i1) (Index.toTIndex i2) +dot (FiniteTensor i1@(Finite.Contravariant count1 _) ts1') (FiniteTensor i2@(Finite.Covariant count2 _) ts2') + | count1 == count2 = Boxed.sum $ Boxed.zipWith (*) ts1' ts2' + | otherwise = contractionErr "finite-finite" (Index.toTIndex i1) (Index.toTIndex i2) +dot t1@(FiniteTensor _ _) t2@(FiniteTensor _ _) = zipT (*) t1 t2 + +-- Other cases cannot happen! +dot t1' t2' = contractionErr "other" (tensorIndex t1') (tensorIndex t2') + +-- | contraction error +{-# INLINE contractionErr #-} +contractionErr :: String -- ^ dot function variant where error occured + -> Index.TIndex -- ^ Index of first dot product parameter + -> Index.TIndex -- ^ Index of second dot product parameter + -> Tensor a -- ^ Erorr message + +contractionErr variant i1' i2' = error $ + "Tensor product: " ++ variant ++ " - " ++ incompatibleTypes ++ + " - index1 is " ++ show i1' ++ + " and index2 is " ++ show i2' + +-- | zipping error +{-# INLINE zipErr #-} +zipErr :: String -- ^ zipT function variant where the error occured + -> Index.TIndex -- ^ Index of first dot product parameter + -> Index.TIndex -- ^ Index of second dot product parameter + -> Tensor a -- ^ Erorr message +zipErr variant i1' i2' = error $ + "zipT: " ++ variant ++ " - " ++ incompatibleTypes ++ + " - index1 is " ++ show i1' ++ + " and index2 is " ++ show i2' + + +-- | Tensors can be added, subtracted and multiplicated +instance (Unboxed.Unbox a, Num a, NFData a) => Num (Tensor a) where + + -- Adding - element by element + {-# INLINE (+) #-} + t1 + t2 = _elemByElem t1 t2 (+) $ zipT (+) + + -- Subtracting - element by element + {-# INLINE (-) #-} + t1 - t2 = _elemByElem t1 t2 (-) $ zipT (-) + + -- Multiplicating is treated as tensor product + -- Tensor product applies Einstein summation convention + {-# INLINE (*) #-} + t1 * t2 = _elemByElem t1 t2 (*) dot + + -- Absolute value - element by element + {-# INLINE abs #-} + abs t = abs `Multilinear.Generic.Sequential.map` t + + -- Signum operation - element by element + {-# INLINE signum #-} + signum t = signum `Multilinear.Generic.Sequential.map` t + + -- Simple integer can be conveted to Scalar + {-# INLINE fromInteger #-} + fromInteger x = Scalar $ fromInteger x + +-- | Tensors can be divided by each other +instance (Unboxed.Unbox a, Fractional a, NFData a) => Fractional (Tensor a) where + -- Tensor dividing: TODO + {-# INLINE (/) #-} + _ / _ = error "TODO" + + -- A scalar can be generated from rational number + {-# INLINE fromRational #-} + fromRational x = Scalar $ fromRational x + +-- Real-number functions on tensors. +-- Function of tensor is tensor of function of its elements +-- E.g. exp [1,2,3,4] = [exp 1, exp2, exp3, exp4] +instance (Unboxed.Unbox a, Floating a, NFData a) => Floating (Tensor a) where + + {-| PI number -} + {-# INLINE pi #-} + pi = Scalar pi + + {-| Exponential function. (exp t)[i] = exp( t[i] ) -} + {-# INLINE exp #-} + exp t = exp `Multilinear.Generic.Sequential.map` t + + {-| Natural logarithm. (log t)[i] = log( t[i] ) -} + {-# INLINE log #-} + log t = log `Multilinear.Generic.Sequential.map` t + + {-| Sinus. (sin t)[i] = sin( t[i] ) -} + {-# INLINE sin #-} + sin t = sin `Multilinear.Generic.Sequential.map` t + + {-| Cosinus. (cos t)[i] = cos( t[i] ) -} + {-# INLINE cos #-} + cos t = cos `Multilinear.Generic.Sequential.map` t + + {-| Inverse sinus. (asin t)[i] = asin( t[i] ) -} + {-# INLINE asin #-} + asin t = asin `Multilinear.Generic.Sequential.map` t + + {-| Inverse cosinus. (acos t)[i] = acos( t[i] ) -} + {-# INLINE acos #-} + acos t = acos `Multilinear.Generic.Sequential.map` t + + {-| Inverse tangent. (atan t)[i] = atan( t[i] ) -} + {-# INLINE atan #-} + atan t = atan `Multilinear.Generic.Sequential.map` t + + {-| Hyperbolic sinus. (sinh t)[i] = sinh( t[i] ) -} + {-# INLINE sinh #-} + sinh t = sinh `Multilinear.Generic.Sequential.map` t + + {-| Hyperbolic cosinus. (cosh t)[i] = cosh( t[i] ) -} + {-# INLINE cosh #-} + cosh t = cosh `Multilinear.Generic.Sequential.map` t + + {-| Inverse hyperbolic sinus. (asinh t)[i] = asinh( t[i] ) -} + {-# INLINE asinh #-} + asinh t = acosh `Multilinear.Generic.Sequential.map` t + + {-| Inverse hyperbolic cosinus. (acosh t)[i] = acosh (t[i] ) -} + {-# INLINE acosh #-} + acosh t = acosh `Multilinear.Generic.Sequential.map` t + + {-| Inverse hyperbolic tangent. (atanh t)[i] = atanh( t[i] ) -} + {-# INLINE atanh #-} + atanh t = atanh `Multilinear.Generic.Sequential.map` t + +-- Multilinear operations +instance (Unboxed.Unbox a) => Multilinear Tensor a where + -- Generic tensor constructor + -- If only one upper index is given, generate a SimpleFinite tensor with upper index + fromIndices [u] [] [s] [] f = + SimpleFinite (Finite.Contravariant s [u]) $ Unboxed.generate s $ \x -> f [x] [] + + -- If only one lower index is given, generate a SimpleFinite tensor with lower index + fromIndices [] [d] [] [s] f = + SimpleFinite (Finite.Covariant s [d]) $ Unboxed.generate s $ \x -> f [] [x] + + -- If many indices are given, first generate upper indices recursively from indices list + fromIndices (u:us) d (s:size) dsize f = + FiniteTensor (Finite.Contravariant s [u]) $ Boxed.generate s (\x -> fromIndices us d size dsize (\uss dss -> f (x:uss) dss) ) + + -- After upper indices, generate lower indices recursively from indices list + fromIndices u (d:ds) usize (s:size) f = + FiniteTensor (Finite.Covariant s [d]) $ Boxed.generate s (\x -> fromIndices u ds usize size (\uss dss -> f uss (x:dss)) ) + + -- If there are indices without size or sizes without names, throw an error + fromIndices _ _ _ _ _ = error "Indices and its sizes incompatible!" + + {-| Accessing tensor elements -} + {-# INLINE el #-} + -- Scalar has only one element + el (Scalar x) _ = Scalar x + -- simple tensor case + el t1@(SimpleFinite index1 _) (inds,vals) = + -- zip indices with their given values + let indvals = zip inds vals + -- find value for simple tensor index if given + val = Data.List.find (\(n,_) -> [n] == Index.indexName index1) indvals + -- if value for current index is given + in if isJust val + -- then get it from current tensor + then t1 ! snd (fromJust val) + -- otherwise return whole tensor - no filtering defined + else t1 + -- finite tensor case + el t1@(FiniteTensor index1 v1) (inds,vals) = + -- zip indices with their given values + let indvals = zip inds vals + -- find value for current index if given + val = Data.List.find (\(n,_) -> [n] == Index.indexName index1) indvals + -- and remove used index from indices list + indvals1 = Data.List.filter (\(n,_) -> [n] /= Index.indexName index1) indvals + -- indices unused so far + inds1 = Data.List.map fst indvals1 + -- and its corresponding values + vals1 = Data.List.map snd indvals1 + -- if value for current index was given + in if isJust val + -- then get it from current tensor and recursively process other indices + then el (t1 ! snd (fromJust val)) (inds1,vals1) + -- otherwise recursively access elements of all child tensors + else FiniteTensor index1 $ (\t -> el t (inds,vals)) <$> v1 + + -- List of all tensor indices + {-# INLINE indices #-} + indices x = case x of + Scalar _ -> [] + FiniteTensor i ts -> Index.toTIndex i : indices (head $ toList ts) + SimpleFinite i _ -> [Index.toTIndex i] + + -- Get tensor order [ (contravariant,covariant)-type ] + {-# INLINE order #-} + order x = case x of + Scalar _ -> (0,0) + SimpleFinite index _ -> case index of + Finite.Contravariant _ _ -> (1,0) + Finite.Covariant _ _ -> (0,1) + _ -> let (cnvr, covr) = order $ firstTensor x + in case tensorIndex x of + Index.Contravariant _ _ -> (cnvr+1,covr) + Index.Covariant _ _ -> (cnvr,covr+1) + + -- Get size of tensor index or Left if index is infinite or tensor has no such index + {-# INLINE size #-} + size t iname = case t of + Scalar _ -> error scalarIndices + SimpleFinite index _ -> + if Index.indexName index == iname + then Finite.indexSize index + else error indexNotFound + FiniteTensor index _ -> + if Index.indexName index == iname + then Finite.indexSize index + else size (firstTensor t) iname + + -- Rename tensor indices + {-# INLINE ($|) #-} + + Scalar x $| _ = Scalar x + SimpleFinite (Finite.Contravariant isize _) ts $| (u:_, _) = SimpleFinite (Finite.Contravariant isize [u]) ts + SimpleFinite (Finite.Covariant isize _) ts $| (_, d:_) = SimpleFinite (Finite.Covariant isize [d]) ts + FiniteTensor (Finite.Contravariant isize _) ts $| (u:us, ds) = FiniteTensor (Finite.Contravariant isize [u]) $ ($| (us,ds)) <$> ts + FiniteTensor (Finite.Covariant isize _) ts $| (us, d:ds) = FiniteTensor (Finite.Covariant isize [d]) $ ($| (us,ds)) <$> ts + t $| _ = t + + -- Raise an index + {-# INLINE (/\) #-} + Scalar x /\ _ = Scalar x + FiniteTensor index ts /\ n + | Index.indexName index == n = + FiniteTensor (Finite.Contravariant (Finite.indexSize index) n) $ (/\ n) <$> ts + | otherwise = + FiniteTensor index $ (/\ n) <$> ts + t1@(SimpleFinite index ts) /\ n + | Index.indexName index == n = + SimpleFinite (Finite.Contravariant (Finite.indexSize index) n) ts + | otherwise = t1 + + -- Lower an index + {-# INLINE (\/) #-} + Scalar x \/ _ = Scalar x + FiniteTensor index ts \/ n + | Index.indexName index == n = + FiniteTensor (Finite.Covariant (Finite.indexSize index) n) $ (\/ n) <$> ts + | otherwise = + FiniteTensor index $ (\/ n) <$> ts + t1@(SimpleFinite index ts) \/ n + | Index.indexName index == n = + SimpleFinite (Finite.Covariant (Finite.indexSize index) n) ts + | otherwise = t1 + + {-| Transpose a tensor (switch all indices types) -} + {-# INLINE transpose #-} + transpose (Scalar x) = Scalar x + + transpose (FiniteTensor (Finite.Covariant count name) ts) = + FiniteTensor (Finite.Contravariant count name) (Multilinear.transpose <$> ts) + transpose (FiniteTensor (Finite.Contravariant count name) ts) = + FiniteTensor (Finite.Covariant count name) (Multilinear.transpose <$> ts) + + transpose (SimpleFinite (Finite.Covariant count name) ts) = + SimpleFinite (Finite.Contravariant count name) ts + transpose (SimpleFinite (Finite.Contravariant count name) ts) = + SimpleFinite (Finite.Covariant count name) ts + + + {-| Shift tensor index right -} + {-| Moves given index one level deeper in recursion -} + -- Scalar has no indices to shift + Scalar x `shiftRight` _ = Scalar x + -- Simple tensor has only one index which cannot be shifted + t1@(SimpleFinite _ _) `shiftRight` _ = t1 + -- Finite tensor is shifted by converting to list and transposing it + t1@(FiniteTensor index1 ts1) `shiftRight` ind + -- We don't shift this index + | Data.List.length (indicesNames t1) > 1 && Index.indexName index1 /= ind = + FiniteTensor index1 $ (|>> ind) <$> ts1 + -- We found index to shift + | Data.List.length (indicesNames t1) > 1 && Index.indexName index1 == ind = + -- Next index + let index2 = tensorFiniteIndex (ts1 Boxed.! 0) + -- Elements to transpose + dane = if isSimple (ts1 Boxed.! 0) + then (\un -> Boxed.generate (Unboxed.length un) (\i -> Scalar $ un Unboxed.! i)) <$> + (tensorScalars <$> ts1) + else tensorsFinite <$> ts1 + result = FiniteTensor index2 $ FiniteTensor index1 <$> (_transpose dane) + -- reconstruct tensor with transposed elements + in _mergeScalars result + -- there is only one index and therefore it cannot be shifted + | otherwise = t1 + +-- Add scalar right +{-# INLINE (.+) #-} +(.+) :: ( + Unboxed.Unbox a, Num a + ) => Tensor a + -> a + -> Tensor a +t .+ x = (+x) `Multilinear.Generic.Sequential.map` t + +-- Subtract scalar right +{-# INLINE (.-) #-} +(.-) :: ( + Unboxed.Unbox a, Num a + ) => Tensor a + -> a + -> Tensor a +t .- x = (\p -> p - x) `Multilinear.Generic.Sequential.map` t + +-- Multiplicate by scalar right +{-# INLINE (.*) #-} +(.*) :: ( + Unboxed.Unbox a, Num a + ) => Tensor a + -> a + -> Tensor a +t .* x = (*x) `Multilinear.Generic.Sequential.map` t + +-- Add scalar left +{-# INLINE (+.) #-} +(+.) :: ( + Unboxed.Unbox a, Num a + ) => a + -> Tensor a + -> Tensor a +x +. t = (x+) `Multilinear.Generic.Sequential.map` t + +-- Subtract scalar left +{-# INLINE (-.) #-} +(-.) :: ( + Unboxed.Unbox a, Num a + ) => a + -> Tensor a + -> Tensor a +x -. t = (x-) `Multilinear.Generic.Sequential.map` t + +-- Multiplicate by scalar left +{-# INLINE (*.) #-} +(*.) :: ( + Unboxed.Unbox a, Num a + ) => a + -> Tensor a + -> Tensor a +x *. t = (x*) `Multilinear.Generic.Sequential.map` t + + +{-| Simple mapping -} +{-| @map f t@ returns tensor @t2@ in which @t2[i1,i2,...] = f t[i1,i2,...]@ -} +{-# INLINE map #-} +map :: ( + Unboxed.Unbox a, Unboxed.Unbox b + ) => (a -> b) + -> Tensor a + -> Tensor b +map f x = case x of + -- Mapping scalar simply maps its value + Scalar v -> Scalar $ f v + -- Mapping complex tensor does mapping element by element + SimpleFinite index ts -> SimpleFinite index (f `Unboxed.map` ts) + FiniteTensor index ts -> FiniteTensor index $ Multilinear.Generic.Sequential.map f <$> ts + +{-| Filtering tensor. + Filtering multi-dimensional arrray may be dangerous, as we always assume, + that on each recursion level, all tensors have the same size (length). + To disable invalid filters, filtering is done over indices, not tensor elements. + Filter function takes and index name and index value and if it returns True, this index value remains in result tensor. + This allows to remove whole columns or rows of eg. a matrix: + filter (\i n -> i /= "a" || i > 10) filters all rows of "a" index (because if i /= "a", filter returns True) + and for "a" index filter elements with index value <= 10 + But this disallow to remove particular matrix element. + If for some index all elements are removed, the index itself is removed from tensor. -} +{-# INLINE filter #-} +filter :: ( + Unboxed.Unbox a + ) => (String -> Int -> Bool) -- ^ filter function + -> Tensor a -- ^ tensor to filter + -> Tensor a +filter _ (Scalar x) = Scalar x +filter f (SimpleFinite index ts) = + let iname = Finite.indexName' index + ts' = (\i _ -> f iname i) `Unboxed.ifilter` ts + in SimpleFinite index { Finite.indexSize = Unboxed.length ts' } ts' +filter f (FiniteTensor index ts) = + let iname = Finite.indexName' index + ts' = Multilinear.Generic.Sequential.filter f <$> ((\i _ -> f iname i) `Boxed.ifilter` ts) + ts'' = + (\case + (SimpleFinite _ ts) -> not $ Unboxed.null ts + (FiniteTensor _ ts) -> not $ Boxed.null ts + _ -> error $ "Filter: " ++ tensorOfScalars + ) `Boxed.filter` ts' + in FiniteTensor index { Finite.indexSize = Boxed.length ts'' } ts'' + +{-| Filtering one index of tensor. -} +{-# INLINE filterIndex #-} +filterIndex :: ( + Unboxed.Unbox a + ) => String -- ^ Index name to filter + -> (Int -> Bool) -- ^ filter function + -> Tensor a -- ^ tensor to filter + -> Tensor a +filterIndex iname f = Multilinear.Generic.Sequential.filter (\i n -> i /= iname || f n) + +{-| Zip tensors with binary combinator, assuming they have all indices the same -} +{-# INLINE zipWith' #-} +zipWith' :: ( + Unboxed.Unbox a, Unboxed.Unbox b, Unboxed.Unbox c + ) => (a -> b -> c) + -> Tensor a + -> Tensor b + -> Tensor c +-- Zipping two Scalars simply combines their values +zipWith' f (Scalar x1) (Scalar x2) = Scalar $ f x1 x2 +-- zipping complex tensor with scalar +zipWith' f t (Scalar x) = (`f` x) `_map` t +-- zipping scalar with complex tensor +zipWith' f (Scalar x) t = (x `f`) `_map` t +-- Two simple tensors case +zipWith' f (SimpleFinite index1 v1) (SimpleFinite index2 v2) = + if index1 == index2 then + SimpleFinite index1 $ Unboxed.zipWith f v1 v2 + else zipErr "simple-simple" (Index.toTIndex index1) (Index.toTIndex index2) +--Two finite tensors case +zipWith' f (FiniteTensor index1 v1) (FiniteTensor index2 v2) = + if index1 == index2 then + FiniteTensor index1 $ Boxed.zipWith (Multilinear.Generic.Sequential.zipWith f) v1 v2 + else zipErr "finite-finite" (Index.toTIndex index1) (Index.toTIndex index2) +-- Other cases cannot happen! +zipWith' _ _ _ = error "Invalid indices to peroform zip!" + +{-# INLINE zipWith #-} +zipWith :: ( + Unboxed.Unbox a, Unboxed.Unbox b, Unboxed.Unbox c + ) => (a -> b -> c) + -> Tensor a + -> Tensor b + -> Tensor c +zipWith f t1 t2 = + let t1' = _standardize t1 + t2' = _standardize t2 + in zipWith' f t1' t2' diff --git a/src/Multilinear/Index.hs b/src/Multilinear/Index.hs --- a/src/Multilinear/Index.hs +++ b/src/Multilinear/Index.hs @@ -30,9 +30,6 @@ {-| Returns True if index is upper (contravariant), False otherwise. -} isContravariant :: i -> Bool - {-| Returns True if index if indifferent, False otherwise. -} - isIndifferent :: i -> Bool - {-| Returns True if two indices are equivalent, thus differs only by name, but share same size and type. -} equivI :: i -> i -> Bool @@ -53,18 +50,13 @@ Contravariant { indexSize :: Maybe Int, tIndexName :: String - } | - Indifferent { - indexSize :: Maybe Int, - tIndexName :: String - } + } deriving (Eq, Generic) {-| Show tensor index -} instance Show TIndex where show (Covariant c n) = "[" ++ n ++ ":" ++ show c ++ "]" show (Contravariant c n) = "<" ++ n ++ ":" ++ show c ++ ">" - show (Indifferent c n) = "(" ++ n ++ ":" ++ show c ++ ")" {-| Finite index is a Multilinear.Index instance -} instance Index TIndex where @@ -80,10 +72,6 @@ isContravariant (Contravariant _ _) = True isContravariant _ = False - {-| Return true if index is indifferent |-} - isIndifferent (Indifferent _ _) = True - isIndifferent _ = False - {-| Returns true if two indices are quivalent, i.e. differs only by name, but share same type and size. -} equivI (Covariant count1 _) (Covariant count2 _) | count1 == count2 = True @@ -91,9 +79,6 @@ equivI (Contravariant count1 _) (Contravariant count2 _) | count1 == count2 = True | otherwise = False - equivI (Indifferent count1 _) (Indifferent count2 _) - | count1 == count2 = True - | otherwise = False equivI _ _ = False {-| TIndex must not be converted to TIndex -} @@ -105,5 +90,3 @@ ind1 <= ind2 = tIndexName ind1 <= tIndexName ind2 || (tIndexName ind1 == tIndexName ind2 && indexSize ind1 <= indexSize ind2) - - diff --git a/src/Multilinear/Index/Finite.hs b/src/Multilinear/Index/Finite.hs --- a/src/Multilinear/Index/Finite.hs +++ b/src/Multilinear/Index/Finite.hs @@ -28,10 +28,6 @@ Contravariant { indexSize :: Int, indexName' :: String - } | - Indifferent { - indexSize :: Int, - indexName' :: String } deriving (Eq, Generic) @@ -39,7 +35,6 @@ instance Show Index where show (Covariant c n) = "[" ++ n ++ ":" ++ show c ++ "]" show (Contravariant c n) = "<" ++ n ++ ":" ++ show c ++ ">" - show (Indifferent c n) = "(" ++ n ++ ":" ++ show c ++ ")" {-| Finite index is a Multilinear.Index instance -} instance TIndex.Index Index where @@ -55,10 +50,6 @@ isContravariant (Contravariant _ _) = True isContravariant _ = False - {-| Return true if index is indifferent |-} - isIndifferent (Indifferent _ _) = True - isIndifferent _ = False - {-| Returns true if two indices are quivalent, i.e. differs only by name, but share same type and size. -} equivI (Covariant count1 _) (Covariant count2 _) | count1 == count2 = True @@ -66,15 +57,11 @@ equivI (Contravariant count1 _) (Contravariant count2 _) | count1 == count2 = True | otherwise = False - equivI (Indifferent count1 _) (Indifferent count2 _) - | count1 == count2 = True - | otherwise = False equivI _ _ = False {-| Convert to TIndex type -} toTIndex (Covariant size name) = TIndex.Covariant (Just size) name toTIndex (Contravariant size name) = TIndex.Contravariant (Just size) name - toTIndex (Indifferent size name) = TIndex.Indifferent (Just size) name {-| Indices can be compared by its name and size |-} {-| Used to allow to put tensors to typical ordered containers |-} diff --git a/src/Multilinear/Index/Infinite.hs b/src/Multilinear/Index/Infinite.hs deleted file mode 100644 --- a/src/Multilinear/Index/Infinite.hs +++ /dev/null @@ -1,71 +0,0 @@ -{-| -Module : Multilinear.Index.Infinite -Description : Infinite-dimensional tensor index. -Copyright : (c) Artur M. Brodzki, 2018 -License : BSD3 -Maintainer : artur@brodzki.org -Stability : experimental -Portability : Windows/POSIX - -Infinite-dimensional tensor index. - --} - -module Multilinear.Index.Infinite ( - Index(..) -) where - -import Control.DeepSeq -import GHC.Generics -import qualified Multilinear.Index as TIndex - -{-| Index of infinite-dimensional tensor -} -data Index = - Covariant { - indexName' :: String - } | - Contravariant { - indexName' :: String - } | - Indifferent { - indexName' :: String - } - deriving (Eq, Generic) - -{-| Show instance of Infinite index -} -instance Show Index where - show (Covariant n) = "[" ++ n ++ "]" - show (Contravariant n) = "<" ++ n ++ ">" - show (Indifferent n) = "(" ++ n ++ ")" - -{-| Infinite index is a Multilinear.Index instance -} -instance TIndex.Index Index where - - {-| Index name -} - indexName = indexName' - - {-| Return true if index is covariant |-} - isCovariant (Covariant _) = True - isCovariant _ = False - - {-| Return true if index is contravariant |-} - isContravariant (Contravariant _) = True - isContravariant _ = False - - {-| Return true if index is indifferent |-} - isIndifferent (Indifferent _) = True - isIndifferent _ = False - - {-| Returns true if two indices are quivalent, i.e. differs only by name, but share same type. -} - equivI (Covariant _) (Covariant _) = True - equivI (Contravariant _) (Contravariant _) = True - equivI (Indifferent _) (Indifferent _) = True - equivI _ _ = False - - {-| Convert to TIndex -} - toTIndex (Covariant name) = TIndex.Covariant Nothing name - toTIndex (Contravariant name) = TIndex.Contravariant Nothing name - toTIndex (Indifferent name) = TIndex.Indifferent Nothing name - --- NFData instance -instance NFData Index diff --git a/src/Multilinear/Matrix.hs b/src/Multilinear/Matrix.hs --- a/src/Multilinear/Matrix.hs +++ b/src/Multilinear/Matrix.hs @@ -15,134 +15,33 @@ module Multilinear.Matrix ( -- * Generators Multilinear.Matrix.fromIndices, - Multilinear.Matrix.const, - Multilinear.Matrix.randomDouble, - Multilinear.Matrix.randomDoubleSeed, - Multilinear.Matrix.randomInt, - Multilinear.Matrix.randomIntSeed, + Multilinear.Matrix.const ) where -import Control.Monad.Primitive import qualified Data.Vector.Unboxed as Unboxed -import Multilinear.Generic -import qualified Multilinear.Tensor as Tensor -import Statistics.Distribution +import Multilinear invalidIndices :: String invalidIndices = "Indices and its sizes not compatible with structure of matrix!" {-| Generate matrix as function of its indices -} - fromIndices :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Indices names (one character per index, first character: rows index, second character: columns index) - -> Int -- ^ Number of matrix rows - -> Int -- ^ Number of matrix columns - -> (Int -> Int -> a) -- ^ Generator function - returns a matrix component at @i,j@ - -> Tensor a -- ^ Generated matrix - -fromIndices [u,d] us ds f = Tensor.fromIndices ([u],[us]) ([d],[ds]) $ \[ui] [di] -> f ui di + Num a, Unboxed.Unbox a, Multilinear t a + ) => String -- ^ Indices names (one character per index, first character: rows index, second character: columns index) + -> Int -- ^ Number of matrix rows + -> Int -- ^ Number of matrix columns + -> (Int -> Int -> a) -- ^ Generator function - returns a matrix component at @i,j@ + -> t a -- ^ Generated matrix +fromIndices [u,d] us ds f = Multilinear.fromIndices [u] [d] [us] [ds] $ \[ui] [di] -> f ui di fromIndices _ _ _ _ = error invalidIndices {-| Generate matrix with all components equal to @v@ -} - const :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Indices names (one character per index, first character: rows index, second character: columns index) - -> Int -- ^ Number of matrix rows - -> Int -- ^ Number of matrix columns - -> a -- ^ Value of matrix components - -> Tensor a -- ^ Generated matrix - -const [u,d] us ds = Tensor.const ([u],[us]) ([d],[ds]) + Num a, Unboxed.Unbox a, Multilinear t a + ) => String -- ^ Indices names (one character per index, first character: rows index, second character: columns index) + -> Int -- ^ Number of matrix rows + -> Int -- ^ Number of matrix columns + -> a -- ^ Value of matrix components + -> t a -- ^ Generated matrix +const [u,d] us ds = Multilinear.const [u] [d] [us] [ds] const _ _ _ = \_ -> error invalidIndices - -{-| Generate matrix with random real components with given probability distribution. -The matrix is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDouble :: ( - ContGen d - ) => String -- ^ Indices names (one character per index, first character: rows index, second character: columns index) - -> Int -- ^ Number of matrix rows - -> Int -- ^ Number of matrix columns - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Double) -- ^ Generated matrix - -randomDouble [u,d] us ds = Tensor.randomDouble ([u],[us]) ([d],[ds]) -randomDouble _ _ _ = \_ -> return $ error invalidIndices - -{-| Generate matrix with random integer components with given probability distribution. -The matrix is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomInt :: ( - DiscreteGen d - ) => String -- ^ Indices names (one character per index, first character: rows index, second character: columns index) - -> Int -- ^ Number of matrix rows - -> Int -- ^ Number of matrix columns - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Int) -- ^ Generated matrix - -randomInt [u,d] us ds = Tensor.randomInt ([u],[us]) ([d],[ds]) -randomInt _ _ _ = \_ -> return $ error invalidIndices - -{-| Generate matrix with random real components with given probability distribution and given seed. -The matrix is wrapped in the a monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDoubleSeed :: ( - ContGen d, PrimMonad m - ) => String -- ^ Indices names (one character per index, first character: rows index, second character: columns index) - -> Int -- ^ Number of matrix rows - -> Int -- ^ Number of matrix columns - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Double) -- ^ Generated matrix - -randomDoubleSeed [u,d] us ds = Tensor.randomDoubleSeed ([u],[us]) ([d],[ds]) -randomDoubleSeed _ _ _ = \_ _ -> return $ error invalidIndices - -{-| Generate matrix with random integer components with given probability distribution. and given seed. -The matrix is wrapped in a monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomIntSeed :: ( - DiscreteGen d, PrimMonad m - ) => String -- ^ Indices names (one character per index, first character: rows index, second character: columns index) - -> Int -- ^ Number of matrix rows - -> Int -- ^ Number of matrix columns - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Int) -- ^ Generated matrix - -randomIntSeed [u,d] us ds = Tensor.randomIntSeed ([u],[us]) ([d],[ds]) -randomIntSeed _ _ _ = \_ _ -> return $ error invalidIndices diff --git a/src/Multilinear/NForm.hs b/src/Multilinear/NForm.hs --- a/src/Multilinear/NForm.hs +++ b/src/Multilinear/NForm.hs @@ -15,155 +15,26 @@ module Multilinear.NForm ( -- * Generators Multilinear.NForm.fromIndices, - Multilinear.NForm.const, - Multilinear.NForm.randomDouble, - Multilinear.NForm.randomDoubleSeed, - Multilinear.NForm.randomInt, - Multilinear.NForm.randomIntSeed, - -- * Common cases - Multilinear.NForm.dot, - Multilinear.NForm.cross + Multilinear.NForm.const ) where -import Control.Monad.Primitive import qualified Data.Vector.Unboxed as Unboxed -import Multilinear.Generic -import qualified Multilinear.Tensor as Tensor -import Statistics.Distribution - -invalidIndices :: String -invalidIndices = "Indices and its sizes incompatible with n-form structure!" - -invalidCrossProductIndices :: String -invalidCrossProductIndices = "Indices and its sizes incompatible with cross product structure!" +import Multilinear {-| Generate N-form as function of its indices -} - fromIndices :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Indices names (one characted per index) - -> [Int] -- ^ Indices sizes - -> ([Int] -> a) -- ^ Generator function - -> Tensor a -- ^ Generated N-form - -fromIndices d ds f = Tensor.fromIndices ([],[]) (d,ds) $ \[] -> f + Num a, Unboxed.Unbox a, Multilinear t a + ) => String -- ^ Indices names (one characted per index) + -> [Int] -- ^ Indices sizes + -> ([Int] -> a) -- ^ Generator function + -> t a -- ^ Generated N-form +fromIndices d ds f = Multilinear.fromIndices [] d [] ds $ \[] -> f {-| Generate N-form with all components equal to @v@ -} - const :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Indices names (one characted per index) - -> [Int] -- ^ Indices sizes - -> a -- ^ N-form elements value - -> Tensor a -- ^ Generated N-form - -const d ds = Tensor.const ([],[]) (d,ds) - -{-| Generate n-vector with random real components with given probability distribution. -The n-vector is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDouble :: ( - ContGen d - ) => String -- ^ Indices names (one character per index) - -> [Int] -- ^ Indices sizes - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Double) -- ^ Generated linear functional - -randomDouble d ds = Tensor.randomDouble ([],[]) (d,ds) - -{-| Generate n-vector with random integer components with given probability distribution. -The n-vector is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomInt :: ( - DiscreteGen d - ) => String -- ^ Indices names (one character per index) - -> [Int] -- ^ Indices sizes - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Int) -- ^ Generated n-vector - -randomInt d ds = Tensor.randomInt ([],[]) (d,ds) - -{-| Generate n-vector with random real components with given probability distribution and given seed. -The form is wrapped in a monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDoubleSeed :: ( - ContGen d, PrimMonad m - ) => String -- ^ Index name (one character) - -> [Int] -- ^ Number of elements - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Double) -- ^ Generated n-vector - -randomDoubleSeed d ds = Tensor.randomDoubleSeed ([],[]) (d,ds) - -{-| Generate n-vector with random integer components with given probability distribution and given seed. -The form is wrapped in a monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomIntSeed :: ( - DiscreteGen d, PrimMonad m - ) => String -- ^ Index name (one character) - -> [Int] -- ^ Number of elements - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Int) -- ^ Generated n-vector - -randomIntSeed d ds = Tensor.randomIntSeed ([],[]) (d,ds) - -{-| 2-form representing a dot product -} - -dot :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Indices names (one characted per index) - -> Int -- ^ Size of tensor (dot product is a square tensor) - -> Tensor a -- ^ Generated dot product - -dot [i1,i2] size = fromIndices [i1,i2] [size,size] (\[i,j] -> if i == j then 1 else 0) -dot _ _ = error invalidIndices - -{-| Tensor representing a cross product (Levi - Civita symbol). It also allows to compute a determinant of square matrix - determinant of matrix @M@ is a equal to length of cross product of all columns of @M@ -} --- // TODO - -cross :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Indices names (one characted per index) - -> Int -- ^ Size of tensor (dot product is a square tensor) - -> Tensor a -- ^ Generated dot product - -cross [i,j,k] size = - Tensor.fromIndices ([i],[size]) ([j,k],[size,size]) - (\[_] [_,_] -> 0) -cross _ _ = error invalidCrossProductIndices + Num a, Unboxed.Unbox a, Multilinear t a + ) => String -- ^ Indices names (one characted per index) + -> [Int] -- ^ Indices sizes + -> a -- ^ N-form elements value + -> t a -- ^ Generated N-form +const d = Multilinear.const [] d [] diff --git a/src/Multilinear/NVector.hs b/src/Multilinear/NVector.hs --- a/src/Multilinear/NVector.hs +++ b/src/Multilinear/NVector.hs @@ -15,121 +15,26 @@ module Multilinear.NVector ( -- * Generators Multilinear.NVector.fromIndices, - Multilinear.NVector.const, - Multilinear.NVector.randomDouble, - Multilinear.NVector.randomDoubleSeed, - Multilinear.NVector.randomInt, - Multilinear.NVector.randomIntSeed, + Multilinear.NVector.const ) where -import Control.Monad.Primitive import qualified Data.Vector.Unboxed as Unboxed -import Multilinear.Generic -import Multilinear.Tensor as Tensor -import Statistics.Distribution +import Multilinear {-| Generate n-vector as function of its indices -} - fromIndices :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Indices names (one characted per index) - -> [Int] -- ^ Indices sizes - -> ([Int] -> a) -- ^ Generator function - -> Tensor a -- ^ Generated n-vector - -fromIndices u us f = Tensor.fromIndices (u,us) ([],[]) $ \uis [] -> f uis + Num a, Unboxed.Unbox a, Multilinear t a + ) => String -- ^ Indices names (one characted per index) + -> [Int] -- ^ Indices sizes + -> ([Int] -> a) -- ^ Generator function + -> t a -- ^ Generated n-vector +fromIndices u usize f = Multilinear.fromIndices u [] usize [] $ \uis [] -> f uis {-| Generate n-vector with all components equal to @v@ -} - const :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Indices names (one characted per index) - -> [Int] -- ^ Indices sizes - -> a -- ^ n-vector elements value - -> Tensor a -- ^ Generated n-vector - -const u us = Tensor.const (u,us) ([],[]) - -{-| Generate n-vector with random real components with given probability distribution. -The n-vector is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDouble :: ( - ContGen d - ) => String -- ^ Indices names (one character per index) - -> [Int] -- ^ Indices sizes - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Double) -- ^ Generated linear functional - -randomDouble u us = Tensor.randomDouble (u,us) ([],[]) - -{-| Generate n-vector with random integer components with given probability distribution. -The n-vector is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomInt :: ( - DiscreteGen d - ) => String -- ^ Indices names (one character per index) - -> [Int] -- ^ Indices sizes - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Int) -- ^ Generated n-vector - -randomInt u us = Tensor.randomInt (u,us) ([],[]) - -{-| Generate n-vector with random real components with given probability distribution and given seed. -The form is wrapped in a monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDoubleSeed :: ( - ContGen d, PrimMonad m - ) => String -- ^ Index name (one character) - -> [Int] -- ^ Number of elements - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Double) -- ^ Generated n-vector - -randomDoubleSeed u us = Tensor.randomDoubleSeed (u,us) ([],[]) - -{-| Generate n-vector with random integer components with given probability distribution and given seed. -The form is wrapped in a monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomIntSeed :: ( - DiscreteGen d, PrimMonad m - ) => String -- ^ Index name (one character) - -> [Int] -- ^ Number of elements - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Int) -- ^ Generated n-vector - -randomIntSeed u us = Tensor.randomIntSeed (u,us) ([],[]) + Num a, Unboxed.Unbox a, Multilinear t a + ) => String -- ^ Indices names (one characted per index) + -> [Int] -- ^ Indices sizes + -> a -- ^ n-vector elements value + -> t a -- ^ Generated n-vector +const u usize = Multilinear.const u [] usize [] diff --git a/src/Multilinear/Tensor.hs b/src/Multilinear/Tensor.hs --- a/src/Multilinear/Tensor.hs +++ b/src/Multilinear/Tensor.hs @@ -14,62 +14,24 @@ module Multilinear.Tensor ( -- * Generators - Multilinear.Tensor.fromIndices, - Multilinear.Tensor.generate, - Multilinear.Tensor.const, - Multilinear.Tensor.randomDouble, - Multilinear.Tensor.randomDoubleSeed, - Multilinear.Tensor.randomInt, - Multilinear.Tensor.randomIntSeed + Multilinear.Tensor.generate ) where -import Control.Monad.Primitive import qualified Data.Vector as Boxed import qualified Data.Vector.Unboxed as Unboxed import Multilinear.Generic import Multilinear.Index.Finite as Finite -import Statistics.Distribution -import qualified System.Random.MWC as MWC invalidIndices :: (String, [Int]) -> (String, [Int]) -> String invalidIndices us ds = "Indices and its sizes incompatible, upper indices: " ++ show us ++", lower indices: " ++ show ds -{-| Generate tensor as functions of its indices -} - -fromIndices :: ( - Num a, Unboxed.Unbox a - ) => (String,[Int]) -- ^ Upper indices names (one character per index) and its sizes - -> (String,[Int]) -- ^ Lower indices names (one character per index) and its sizes - -> ([Int] -> [Int] -> a) -- ^ Generator function (f [u1,u2,...] [d1,d2,...] returns a tensor element at t [u1,u2,...] [d1,d2,...]) - -> Tensor a -- ^ Generated tensor - --- If only one upper index is given, generate a SimpleFinite tensor with upper index -fromIndices ([u],[s]) ([],[]) f = - SimpleFinite (Contravariant s [u]) $ Unboxed.generate s $ \x -> f [x] [] - --- If only one lower index is given, generate a SimpleFinite tensor with lower index -fromIndices ([],[]) ([d],[s]) f = - SimpleFinite (Covariant s [d]) $ Unboxed.generate s $ \x -> f [] [x] - --- If many indices are given, first generate upper indices recursively from indices list -fromIndices (u:us,s:size) d f = - FiniteTensor (Contravariant s [u]) $ Boxed.generate s (\x -> fromIndices (us,size) d (\uss dss -> f (x:uss) dss) ) - --- After upper indices, generate lower indices recursively from indices list -fromIndices u (d:ds,s:size) f = - FiniteTensor (Covariant s [d]) $ Boxed.generate s (\x -> fromIndices u (ds,size) (\uss dss -> f uss (x:dss)) ) - --- If there are indices without size or sizes without names, throw an error -fromIndices us ds _ = error $ invalidIndices us ds - {-| Generate tensor composed of other tensors -} - generate :: ( Num a, Unboxed.Unbox a - ) => (String,[Int]) -- ^ Upper indices names (one character per index) and its sizes - -> (String,[Int]) -- ^ Lower indices names (one character per index) and its sizes - -> ([Int] -> [Int] -> Tensor a) -- ^ Generator function (f [u1,u2,...] [d1,d2,...] returns a tensor element at t [u1,u2,...] [d1,d2,...]) - -> Tensor a -- ^ Generated tensor + ) => (String,[Int]) -- ^ Upper indices names (one character per index) and its sizes + -> (String,[Int]) -- ^ Lower indices names (one character per index) and its sizes + -> ([Int] -> [Int] -> Tensor a) -- ^ Generator function (f [u1,u2,...] [d1,d2,...] returns a tensor element at t [u1,u2,...] [d1,d2,...]) + -> Tensor a -- ^ Generated tensor -- If no indices are given, generate a tensor by using generator function generate ([],[]) ([],[]) f = f [] [] @@ -84,215 +46,3 @@ -- If there are indices without size or sizes without names, throw an error generate us ds _ = error $ invalidIndices us ds - -{-| Generate tensor with all components equal to @v@ -} - -const :: ( - Num a, Unboxed.Unbox a - ) => (String,[Int]) -- ^ Upper indices names (one character per index) and its sizes - -> (String,[Int]) -- ^ Lower indices names (one character per index) and its sizes - -> a -- ^ Tensor elements value - -> Tensor a -- ^ Generated tensor - --- If only one upper index is given, generate a SimpleFinite tensor with upper index -const ([u],[s]) ([],[]) v = - SimpleFinite (Contravariant s [u]) $ Unboxed.replicate s v - --- If only ine lower index is given, generate a SimpleFinite tensor with lower index -const ([],[]) ([d],[s]) v = - SimpleFinite (Covariant s [d]) $ Unboxed.replicate s v - --- If many indices are given, first generate upper indices recursively from indices list -const (u:us,s:size) d v = - FiniteTensor (Contravariant s [u]) $ Boxed.replicate (fromIntegral s) $ Multilinear.Tensor.const (us,size) d v - --- After upper indices, generate lower indices recursively from indices list -const u (d:ds,s:size) v = - FiniteTensor (Covariant s [d]) $ Boxed.replicate (fromIntegral s) $ Multilinear.Tensor.const u (ds,size) v - --- If there are indices without size or sizes without names, throw an error -const us ds _ = error $ invalidIndices us ds - -{-| Generate tensor with random real components with given probability distribution. -The tensor is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDouble :: ( - ContGen d - ) => (String,[Int]) -- ^ Upper indices names (one character per index) and its sizes - -> (String,[Int]) -- ^ Lower indices names (one character per index) and its sizes - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Double) -- ^ Generated tensor - --- If only one upper index is given, generate a SimpleFinite tensor with upper index -randomDouble ([u],[s]) ([],[]) distr = do - gen <- MWC.createSystemRandom - component' <- sequence $ Boxed.generate s $ \_ -> genContVar distr gen - let component = Unboxed.generate s $ \i -> component' Boxed.! i - return $ SimpleFinite (Contravariant s [u]) component - --- If only one lower index is given, generate a SimpleFinite tensor with lower index -randomDouble ([],[]) ([d],[s]) distr = do - gen <- MWC.createSystemRandom - component' <- sequence $ Boxed.generate s $ \_ -> genContVar distr gen - let component = Unboxed.generate s $ \i -> component' Boxed.! i - return $ SimpleFinite (Covariant s [d]) component - --- If many indices are given, first generate upper indices recursively from indices list -randomDouble (u:us,s:size) d distr = do - tensors <- sequence $ Boxed.generate s $ \_ -> randomDouble (us,size) d distr - return $ FiniteTensor (Contravariant s [u]) tensors - --- After upper indices, generate lower indices recursively from indices list -randomDouble u (d:ds,s:size) distr = do - tensors <- sequence $ Boxed.generate s $ \_ -> randomDouble u (ds,size) distr - return $ FiniteTensor (Covariant s [d]) tensors - --- If there are indices without size or sizes without names, throw an error -randomDouble us ds _ = return $ error $ invalidIndices us ds - -{-| Generate tensor with random integer components with given probability distribution. -The tensor is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomInt :: ( - DiscreteGen d - ) => (String,[Int]) -- ^ Upper indices names (one character per index) and its sizes - -> (String,[Int]) -- ^ Lower indices names (one character per index) and its sizes - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Int) -- ^ Generated tensor - --- If only one upper index is given, generate a SimpleFinite tensor with upper index -randomInt ([u],[s]) ([],[]) distr = do - gen <- MWC.createSystemRandom - component' <- sequence $ Boxed.generate s $ \_ -> genDiscreteVar distr gen - let component = Unboxed.generate s $ \i -> component' Boxed.! i - return $ SimpleFinite (Contravariant s [u]) component - --- If only one lower index is given, generate a SimpleFinite tensor with lower index -randomInt ([],[]) ([d],[s]) distr = do - gen <- MWC.createSystemRandom - component' <- sequence $ Boxed.generate s $ \_ -> genDiscreteVar distr gen - let component = Unboxed.generate s $ \i -> component' Boxed.! i - return $ SimpleFinite (Covariant s [d]) component - --- If many indices are given, first generate upper indices recursively from indices list -randomInt (u:us,s:size) d distr = do - tensors <- sequence $ Boxed.generate s $ \_ -> randomInt (us,size) d distr - return $ FiniteTensor (Contravariant s [u]) tensors - --- After upper indices, generate lower indices recursively from indices list -randomInt u (d:ds,s:size) distr = do - tensors <- sequence $ Boxed.generate s $ \_ -> randomInt u (ds,size) distr - return $ FiniteTensor (Covariant s [d]) tensors - --- If there are indices without size or sizes without names, throw an error -randomInt us ds _ = return $ error $ invalidIndices us ds - -{-| Generate tensor with random real components with given probability distribution and given seed. -The tensor is wrapped in a monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDoubleSeed :: ( - ContGen d, PrimMonad m - ) => (String,[Int]) -- ^ Upper indices names (one character per index) and its sizes - -> (String,[Int]) -- ^ Lower indices names (one character per index) and its sizes - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Double) -- ^ Generated tensor - --- If only one upper index is given, generate a SimpleFinite tensor with upper index -randomDoubleSeed ([u],[s]) ([],[]) distr seed = do - gen <- MWC.initialize (Boxed.singleton $ fromIntegral seed) - component' <- sequence $ Boxed.generate s $ \_ -> genContVar distr gen - let component = Unboxed.generate s $ \i -> component' Boxed.! i - return $ SimpleFinite (Contravariant s [u]) component - --- If only one lower index is given, generate a SimpleFinite tensor with lower index -randomDoubleSeed ([],[]) ([d],[s]) distr seed = do - gen <- MWC.initialize (Boxed.singleton $ fromIntegral seed) - component' <- sequence $ Boxed.generate s $ \_ -> genContVar distr gen - let component = Unboxed.generate s $ \i -> component' Boxed.! i - return $ SimpleFinite (Covariant s [d]) component - --- If many indices are given, first generate upper indices recursively from indices list -randomDoubleSeed (u:us,s:size) d distr seed = do - tensors <- sequence $ Boxed.generate s $ \_ -> randomDoubleSeed (us,size) d distr seed - return $ FiniteTensor (Contravariant s [u]) tensors - --- After upper indices, generate lower indices recursively from indices list -randomDoubleSeed u (d:ds,s:size) distr seed = do - tensors <- sequence $ Boxed.generate s $ \_ -> randomDoubleSeed u (ds,size) distr seed - return $ FiniteTensor (Covariant s [d]) tensors - --- If there are indices without size or sizes without names, throw an error -randomDoubleSeed us ds _ _ = return $ error $ invalidIndices us ds - -{-| Generate tensor with random integer components with given probability distribution and given seed. -The tensor is wrapped in a monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomIntSeed :: ( - DiscreteGen d, PrimMonad m - ) => (String,[Int]) -- ^ Index name (one character) - -> (String,[Int]) -- ^ Number of elements - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Int) -- ^ Generated tensor - --- If only one upper index is given, generate a SimpleFinite tensor with upper index -randomIntSeed ([u],[s]) ([],[]) distr seed = do - gen <- MWC.initialize (Boxed.singleton $ fromIntegral seed) - component' <- sequence $ Boxed.generate s $ \_ -> genDiscreteVar distr gen - let component = Unboxed.generate s $ \i -> component' Boxed.! i - return $ SimpleFinite (Contravariant s [u]) component - --- If only one lower index is given, generate a SimpleFinite tensor with lower index -randomIntSeed ([],[]) ([d],[s]) distr seed = do - gen <- MWC.initialize (Boxed.singleton $ fromIntegral seed) - component' <- sequence $ Boxed.generate s $ \_ -> genDiscreteVar distr gen - let component = Unboxed.generate s $ \i -> component' Boxed.! i - return $ SimpleFinite (Covariant s [d]) component - --- If many indices are given, first generate upper indices recursively from indices list -randomIntSeed (u:us,s:size) d distr seed = do - tensors <- sequence $ Boxed.generate s $ \_ -> randomIntSeed (us,size) d distr seed - return $ FiniteTensor (Contravariant s [u]) tensors - --- After upper indices, generate lower indices recursively from indices list -randomIntSeed u (d:ds,s:size) distr seed = do - tensors <- sequence $ Boxed.generate s $ \_ -> randomIntSeed u (ds,size) distr seed - return $ FiniteTensor (Covariant s [d]) tensors - --- If there are indices without size or sizes without names, throw an error -randomIntSeed us ds _ _ = return $ error $ invalidIndices us ds diff --git a/src/Multilinear/Vector.hs b/src/Multilinear/Vector.hs --- a/src/Multilinear/Vector.hs +++ b/src/Multilinear/Vector.hs @@ -15,129 +15,31 @@ module Multilinear.Vector ( -- * Generators Multilinear.Vector.fromIndices, - Multilinear.Vector.const, - Multilinear.Vector.randomDouble, - Multilinear.Vector.randomDoubleSeed, - Multilinear.Vector.randomInt, - Multilinear.Vector.randomIntSeed + Multilinear.Vector.const ) where -import Control.Monad.Primitive import qualified Data.Vector.Unboxed as Unboxed -import Multilinear.Generic -import Multilinear.Tensor as Tensor -import Statistics.Distribution +import Multilinear invalidIndices :: String invalidIndices = "Indices and its sizes not compatible with structure of vector!" {-| Generate vector as function of indices -} - fromIndices :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> (Int -> a) -- ^ Generator function - returns a vector component at index @i@ - -> Tensor a -- ^ Generated vector - -fromIndices [i] s f = Tensor.fromIndices ([i],[s]) ([],[]) $ \[x] [] -> f x + Num a, Unboxed.Unbox a, Multilinear t a + ) => String -- ^ Index name (one character) + -> Int -- ^ Number of elements + -> (Int -> a) -- ^ Generator function - returns a vector component at index @i@ + -> t a -- ^ Generated vector +fromIndices [i] s f = Multilinear.fromIndices [i] [] [s] [] $ \[x] [] -> f x fromIndices _ _ _ = error invalidIndices {-| Generate vector with all components equal to some @v@ -} - const :: ( - Num a, Unboxed.Unbox a - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> a -- ^ Value of each element - -> Tensor a -- ^ Generated vector - -const [i] s = Tensor.const ([i],[s]) ([],[]) + Num a, Unboxed.Unbox a, Multilinear t a + ) => String -- ^ Index name (one character) + -> Int -- ^ Number of elements + -> a -- ^ Value of each element + -> t a -- ^ Generated vector +const [i] s = Multilinear.const [i] [] [s] [] const _ _ = \_ -> error invalidIndices - -{-| Generate vector with random real components with given probability distribution. -The vector is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDouble :: ( - ContGen d - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Double) -- ^ Generated vector - -randomDouble [i] s = Tensor.randomDouble ([i],[s]) ([],[]) -randomDouble _ _ = \_ -> return $ error invalidIndices - -{-| Generate vector with random integer components with given probability distribution. -The vector is wrapped in the IO monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomInt :: ( - DiscreteGen d - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> IO (Tensor Int) -- ^ Generated vector - -randomInt [i] s = Tensor.randomInt ([i],[s]) ([],[]) -randomInt _ _ = \_ -> return $ error invalidIndices - -{-| Generate vector with random real components with given probability distribution and given seed. -The vector is wrapped in a monad. -} -{-| Available probability distributions: -} -{-| - Beta : "Statistics.Distribution.BetaDistribution" -} -{-| - Cauchy : "Statistics.Distribution.CauchyLorentz" -} -{-| - Chi-squared : "Statistics.Distribution.ChiSquared" -} -{-| - Exponential : "Statistics.Distribution.Exponential" -} -{-| - Gamma : "Statistics.Distribution.Gamma" -} -{-| - Normal : "Statistics.Distribution.Normal" -} -{-| - StudentT : "Statistics.Distribution.StudentT" -} -{-| - Uniform : "Statistics.Distribution.Uniform" -} -{-| - F : "Statistics.Distribution.FDistribution" -} -{-| - Laplace : "Statistics.Distribution.Laplace" -} - -randomDoubleSeed :: ( - ContGen d, PrimMonad m - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> d -- ^ Continuous probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Double) -- ^ Generated vector - -randomDoubleSeed [i] s = Tensor.randomDoubleSeed ([i],[s]) ([],[]) -randomDoubleSeed _ _ = \_ _ -> return $ error invalidIndices - -{-| Generate vector with random integer components with given probability distribution and given seed. -The vector is wrapped in a monad. -} -{-| Available probability distributions: -} -{-| - Binomial : "Statistics.Distribution.Binomial" -} -{-| - Poisson : "Statistics.Distribution.Poisson" -} -{-| - Geometric : "Statistics.Distribution.Geometric" -} -{-| - Hypergeometric: "Statistics.Distribution.Hypergeometric" -} - -randomIntSeed :: ( - DiscreteGen d, PrimMonad m - ) => String -- ^ Index name (one character) - -> Int -- ^ Number of elements - -> d -- ^ Discrete probability distribution (as from "Statistics.Distribution") - -> Int -- ^ Randomness seed - -> m (Tensor Int) -- ^ Generated vector - -randomIntSeed [i] s = Tensor.randomIntSeed ([i],[s]) ([],[]) -randomIntSeed _ _ = \_ _ -> return $ error invalidIndices - diff --git a/test/Spec.hs b/test/Spec.hs deleted file mode 100644 --- a/test/Spec.hs +++ /dev/null @@ -1,257 +0,0 @@ -{-| -Module : Main -Description : Test of Multilinear library -Copyright : (c) Artur M. Brodzki, 2018 -License : BSD3 -Maintainer : artur@brodzki.org -Stability : experimental -Portability : Windows/POSIX - --} - -module Main ( - main -) where - -import Data.Maybe -import qualified Data.Set as Set -import Multilinear -import Multilinear.Generic -import qualified Multilinear.Index as Index -import System.Exit -import System.IO -import Test.QuickCheck -import Test.QuickCheck.Multilinear() - --- | Default test number for property -defTestN :: Int -defTestN = 1000 - - ------------------------------- --- AUXILIARY TEST FUNCTIONS -- ------------------------------- - - --- quickCheck with parametrizable tests number -quickCheckN :: Testable prop => Int -> prop -> IO Result -quickCheckN n = quickCheckWithResult (Args - Nothing -- ^ Should we replay a previous test? No. - n -- ^ Maximum number of successful tests before succeeding set to N. - 1 -- ^ Maximum number of discarded tests per successful test before giving up - gave up after first failure. - n -- ^ Size to use for the biggest test cases. - True -- ^ Whether to print anything? yes. - 0) -- ^ Maximum number of shrinks to before giving up. Turn shrinking off. - --- | Execute property test and check result: --- | exit test suite with successs code if no errors occured --- | exit test suite with failure code if any error occured -executePropertyTest :: ( - Testable prop - ) => String -- ^ Tested property name - -> Int -- ^ Number of tests to do - -> prop -- ^ Property to test - -> IO () -executePropertyTest propName n f = do - putStr $ " Checking " ++ propName ++ " " - r <- quickCheckN n f - case r of - Success _ _ _ -> hFlush stdout - _ -> exitFailure - - ------------------------------- --- TESTED TENSOR PROPERTIES -- ------------------------------- - - --- | Unary operator applied on any tensor, --- | must preserve tensor indices in the result. -preserveIndicesUnary :: - (Tensor Double -> - Tensor Double) -- ^ Unary tensor operator to test - -> Tensor Double -- ^ Operator argument - -> Bool -preserveIndicesUnary f t = indices t == indices (f t) - --- | Binary operator applied on any two tensors which have all the same indices, --- | must preserve set union of these indices in the result. -preserveIndicesBinary :: - (Tensor Double -> - Tensor Double -> - Tensor Double) -- ^ Binary tensor operator to test - -> Tensor Double -- ^ First operator argument - -> Tensor Double -- ^ Second operator argument - -> Bool -preserveIndicesBinary f t1 t2 = - let i1 = Set.fromList $ indices t1 - i2 = Set.fromList $ indices t2 - in i1 /= i2 || i1 == Set.fromList (indices $ f t1 t2) - --- | Binary operator other than tensor product cannot contract (or consume) any index --- | it means, that in operators other than (*), the indices of result tensor are set union of arguments indices -mergeCommonIndices :: - (Tensor Double -> - Tensor Double -> - Tensor Double) -- ^ Binary tensor operator to test - -> Tensor Double -- ^ First operator argument - -> Tensor Double -- ^ Second operator argument - -> Bool -mergeCommonIndices f t1 t2 = - let indices1 = Set.fromList $ indices t1 - indices2 = Set.fromList $ indices t2 - inames1 = Set.fromList $ Index.indexName <$> indices t1 - inames2 = Set.fromList $ Index.indexName <$> indices t2 - - commonIndices = Set.intersection indices1 indices2 - commonIndicesNames = Set.intersection inames1 inames2 - - expectedIndices = Set.union inames1 inames2 - resultIndices = Set.fromList $ Index.indexName <$> indices (f t1 t2) - - -- if we have indices, which have the same name but different type, it is forbidden and test passed - in Set.size commonIndices /= Set.size commonIndicesNames || - -- otherwise, the result indices set must be union of arguments indices - expectedIndices == resultIndices - - --- | Contracted indices have to be consumed in result tensor. -consumeContractedIndices :: - Tensor Double -- ^ first tensor to contract - -> Tensor Double -- ^ second tensor to contract - -> Bool -consumeContractedIndices t1 t2 = - let inames1 = Set.fromList $ Index.indexName <$> indices t1 - inames2 = Set.fromList $ Index.indexName <$> indices t2 - - iContravariantNames1 = Set.fromList $ Index.indexName <$> (Index.isContravariant `Prelude.filter` indices t1) - iCovariantNames1 = Set.fromList $ Index.indexName <$> (Index.isCovariant `Prelude.filter` indices t1) - - iContravariantNames2 = Set.fromList $ Index.indexName <$> (Index.isContravariant `Prelude.filter` indices t2) - iCovariantNames2 = Set.fromList $ Index.indexName <$> (Index.isCovariant `Prelude.filter` indices t2) - - contractedIndices = - -- contracted are indices covariant in the first tensor and contravariant in the second - Set.intersection iCovariantNames1 iContravariantNames2 `Set.union` - -- or contravariant in the first tensor and covariant in the second - Set.intersection iContravariantNames1 iCovariantNames2 - - expectedIndices = Set.difference (Set.union inames1 inames2) contractedIndices - resultIndices = Set.fromList $ Index.indexName <$> indices (t1 * t2) - - in expectedIndices == resultIndices - --- | Order of the tensor must be equal to number of its covariant and contravariant indices -orderIndices :: Tensor Double -> Bool -orderIndices t = - let (conv, cov) = order t - iConv = Set.fromList $ Index.isContravariant `Prelude.filter` indices t - iCov = Set.fromList $ Index.isCovariant `Prelude.filter` indices t - in conv == Set.size iConv && cov == Set.size iCov - --- | Tensor must be equivalent in terms of its indices after any index shift -shiftEquiv :: Tensor Double -> Bool -shiftEquiv t = - let inames = indicesNames t - rShiftedTs = (\i -> t |>> i) <$> inames - lShiftedTs = (\i -> t <<| i) <$> inames - rtShiftedTs = (\i -> t |>>> i) <$> inames - ltShiftedTs = (\i -> t <<<| i) <$> inames - allShiftedTs = rShiftedTs ++ lShiftedTs ++ rtShiftedTs ++ ltShiftedTs ++ [t] - allPairs = pure (,) <*> allShiftedTs <*> allShiftedTs - in all (uncurry (|==|)) allPairs - -{-| After rename, index must hold a new name - This property assumes, tensor have max 5 indices of each type -} -renameTest :: Tensor Double -> Bool -renameTest t = - let (conv, cov) = order t - convNs = take conv ['m' .. ] - covNs = take cov ['s' .. ] - renamedT = t $| (convNs, covNs) - inamesAfter = concat $ indicesNames renamedT - in all (\i -> elem i convNs || elem i covNs) inamesAfter - --- | After any raising or lowering index, it must be a valid type. -raiseLowerTest :: Tensor Double -> Bool -raiseLowerTest t = - let inames = indicesNames t - lowered = inames `zip` ((t \/) <$> inames) - raised = inames `zip` ((t /\) <$> inames) - isLowered (i,tl) = i `elem` (Index.indexName <$> (Index.isCovariant `Prelude.filter` indices tl)) - isRaised (i,tr) = i `elem` (Index.indexName <$> (Index.isContravariant `Prelude.filter` indices tr)) - in all isLowered lowered && all isRaised raised - - --- | Filter second half of elements for each tensor index and check if they disappeared -filterIndexTest :: - Tensor Double -> Bool -filterIndexTest s@(Scalar _) = s == filterIndex "c" (const True) s -filterIndexTest t = - let indsT = indices t - -- filter second half of an index - filteredHalf i = filterIndex (Index.indexName i) (< (fromJust (Index.indexSize i) `div` 2)) t - fts = indsT `zip` (filteredHalf <$> indsT) -- tensors with filtered indices, paired with respective transformed indices - in all (\(i,ft) -> - size ft (Index.indexName i) == (fromJust (Index.indexSize i) `div` 2) - ) fts - --- | ENTRY POINT -main :: IO () -main = do - - -- PRINT PROBABILITY DISTRIBUTION OF TESTED TENSORS ORDER - executePropertyTest "probability distribution of order of tested tensors" 5000 $ - \(t :: Tensor Double) -> collect (order t) $ preserveIndicesUnary abs - - putStrLn "\nTesting multilinear library...\n" - - --------------------------- - -- CHECKING NUM INSTANCE -- - --------------------------- - - executePropertyTest "preserveIndicesBinary for (+)" defTestN $ preserveIndicesBinary (+) - executePropertyTest "preserveIndicesBinary for (-)" defTestN $ preserveIndicesBinary (-) - executePropertyTest "preserveIndicesBinary for (*)" defTestN $ preserveIndicesBinary (*) - executePropertyTest "preserveIndicesUnary for abs" defTestN $ preserveIndicesUnary abs - executePropertyTest "preserveIndicesUnary for signum" defTestN $ preserveIndicesUnary signum - - executePropertyTest "mergeCommonIndices for (+)" defTestN $ mergeCommonIndices (+) - executePropertyTest "mergeCommonIndices for (-)" defTestN $ mergeCommonIndices (-) - executePropertyTest "consumeContractedIndices" defTestN consumeContractedIndices - - -------------------------------- - -- CHECKING FLOATING INSTANCE -- - -------------------------------- - - executePropertyTest "preserveIndicesUnary for exp" defTestN $ preserveIndicesUnary exp - executePropertyTest "preserveIndicesUnary for log" defTestN $ preserveIndicesUnary log - executePropertyTest "preserveIndicesUnary for sin" defTestN $ preserveIndicesUnary sin - executePropertyTest "preserveIndicesUnary for cos" defTestN $ preserveIndicesUnary cos - executePropertyTest "preserveIndicesUnary for asin" defTestN $ preserveIndicesUnary asin - executePropertyTest "preserveIndicesUnary for acos" defTestN $ preserveIndicesUnary acos - executePropertyTest "preserveIndicesUnary for atan" defTestN $ preserveIndicesUnary atan - executePropertyTest "preserveIndicesUnary for sinh" defTestN $ preserveIndicesUnary sinh - executePropertyTest "preserveIndicesUnary for cosh" defTestN $ preserveIndicesUnary cosh - executePropertyTest "preserveIndicesUnary for asinh" defTestN $ preserveIndicesUnary asinh - executePropertyTest "preserveIndicesUnary for acosh" defTestN $ preserveIndicesUnary acosh - executePropertyTest "preserveIndicesUnary for atanh" defTestN $ preserveIndicesUnary atanh - - ----------------------------------- - -- CHECKING MULTILINEAR INSTANCE -- - ----------------------------------- - - executePropertyTest "preserveIndicesUnary for (+.)" defTestN $ preserveIndicesUnary (5 +.) - executePropertyTest "preserveIndicesUnary for (.+)" defTestN $ preserveIndicesUnary (.+ 5) - executePropertyTest "preserveIndicesUnary for (-.)" defTestN $ preserveIndicesUnary (5 -.) - executePropertyTest "preserveIndicesUnary for (.-)" defTestN $ preserveIndicesUnary (.- 5) - executePropertyTest "preserveIndicesUnary for (*.)" defTestN $ preserveIndicesUnary (5 *.) - executePropertyTest "preserveIndicesUnary for (.*)" defTestN $ preserveIndicesUnary (.* 5) - - executePropertyTest "orderIndices" defTestN orderIndices - executePropertyTest "shiftEquiv" defTestN shiftEquiv - executePropertyTest "renamedTest" defTestN renameTest - executePropertyTest "raiseLowerTest" defTestN raiseLowerTest - executePropertyTest "filterIndexTest" defTestN filterIndexTest - executePropertyTest "zipWithIndicesTest" defTestN $ preserveIndicesUnary (\t -> Multilinear.zipWith (+) t t) - diff --git a/test/Test/QuickCheck/Multilinear.hs b/test/Test/QuickCheck/Multilinear.hs deleted file mode 100644 --- a/test/Test/QuickCheck/Multilinear.hs +++ /dev/null @@ -1,81 +0,0 @@ -{-| -Module : Test.QuickCheck.Multilinear -Description : QucikCheck instances of Multilinear library -Copyright : (c) Artur M. Brodzki, 2018 -License : BSD3 -Maintainer : artur@brodzki.org -Stability : experimental -Portability : Windows/POSIX - --} - -module Test.QuickCheck.Multilinear ( - Arbitrary -) where - -import qualified Multilinear.Form as Form -import Multilinear.Generic -import qualified Multilinear.Vector as Vector -import Test.QuickCheck - --- | Sizes of indices used in Arbitrary Tensor instance -aS :: Int -aS = 12 -bS :: Int -bS = 12 -iS :: Int -iS = 10 -jS :: Int -jS = 15 -kS :: Int -kS = 10 - --- | Set of three Scalars, that can be used for building more complex tensors -scalars :: [Tensor Double] -scalars = [ - -- Scalars - Scalar (-1.0) - , Scalar 0.0 - , Scalar 1.0 - ] - --- | Set of 5 simple Scalars and 1D tensors for testing with only upper indices --- | We use set of a,b,i,j,k indices for further building more complex tensors sets -tensors1DUpper :: [Tensor Double] -tensors1DUpper = [ - -- Vectors with a,b,i,j,k indices - Vector.fromIndices "a" aS $ sin . fromIntegral - , Vector.fromIndices "b" bS $ cos . fromIntegral - , Vector.fromIndices "i" iS $ exp . fromIntegral - , Vector.fromIndices "j" jS $ cosh . fromIntegral - , Vector.fromIndices "k" kS $ tanh . fromIntegral - ] - --- | Set of 5 simple Scalars and 1D tensors for testing with only lower indices --- | We use set of a,b,i,j,k indices for further building more complex tensors sets -tensors1DLower :: [Tensor Double] -tensors1DLower = [ - -- Functional with a,b,i,j,k indices - can be contracted with vectors above or matrices below - Form.fromIndices "a" aS $ sin . fromIntegral - , Form.fromIndices "b" bS $ cos . fromIntegral - , Form.fromIndices "i" iS $ exp . fromIntegral - , Form.fromIndices "j" jS $ cosh . fromIntegral - , Form.fromIndices "k" kS $ tanh . fromIntegral - ] - --- | List sum of scalars and upper and lower indices simple tensors --- | List contains 18 tensors in total -tensors1D :: [Tensor Double] -tensors1D = scalars ++ tensors1DUpper ++ tensors1DLower - -{-| More complex (up to 3D) tensors, built as sums and differences of tensor products of all pairs from tensors list above - List contains 18^3 = 5832 tensors in total -} - -tensors3D :: [Tensor Double] -tensors3D = pure (*) <*> ts <*> tensors1D - where ts = pure (*) <*> tensors1D <*> tensors1D - --- | Arbitrary random generating instance of Tensor Double --- | Simply choose a tensot from tensors list above -instance Arbitrary (Tensor Double) where - arbitrary = elements tensors3D diff --git a/test/multicore/Spec.hs b/test/multicore/Spec.hs new file mode 100644 --- /dev/null +++ b/test/multicore/Spec.hs @@ -0,0 +1,446 @@ +{-| +Module : Main +Description : Test of multicore tensor +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Main ( + main +) where + +import Data.Maybe +import qualified Data.Set as Set +import Multilinear.Class +import Multilinear.Generic.MultiCore +import qualified Multilinear.Index as Index +import qualified Multilinear.Form as Form +import qualified Multilinear.Matrix as Matrix +import qualified Multilinear.Vector as Vector +import qualified Multilinear.NForm as NForm +import qualified Multilinear.NVector as NVector +import System.Exit +import System.IO +import Test.QuickCheck +import Test.QuickCheck.Multilinear.Generic.MultiCore() + +-- | Default test number for property +defTestN :: Int +defTestN = 1000 + + +------------------------------ +-- AUXILIARY TEST FUNCTIONS -- +------------------------------ + + +-- quickCheck with parametrizable tests number +quickCheckN :: Testable prop => Int -> prop -> IO Result +quickCheckN n = quickCheckWithResult (Args + Nothing -- ^ Should we replay a previous test? No. + n -- ^ Maximum number of successful tests before succeeding set to N. + 1 -- ^ Maximum number of discarded tests per successful test before giving up - gave up after first failure. + n -- ^ Size to use for the biggest test cases. + True -- ^ Whether to print anything? yes. + 0) -- ^ Maximum number of shrinks to before giving up. Turn shrinking off. + +-- | Execute property test and check result: +-- | exit test suite with successs code if no errors occured +-- | exit test suite with failure code if any error occured +executePropertyTest :: ( + Testable prop + ) => String -- ^ Tested property name + -> Int -- ^ Number of tests to do + -> prop -- ^ Property to test + -> IO () +executePropertyTest propName n f = do + putStr $ " Checking " ++ propName ++ " " + r <- quickCheckN n f + case r of + Success _ _ _ -> hFlush stdout + _ -> exitFailure + + +---------------------------------------------------- +-- TESTED TENSOR PROPERTIES FOR SEQUENTIAL TENSOR -- +---------------------------------------------------- + + +-- | Unary operator applied on any tensor, +-- | must preserve tensor indices in the result. +preserveIndicesUnary :: + (Tensor Double -> + Tensor Double) -- ^ Unary tensor operator to test + -> Tensor Double -- ^ Operator argument + -> Bool +preserveIndicesUnary f t = Set.fromList (indices t) == Set.fromList (indices (f t)) + +-- | Binary operator applied on any two tensors which have all the same indices, +-- | must preserve set union of these indices in the result. +preserveIndicesBinary :: + (Tensor Double -> + Tensor Double -> + Tensor Double) -- ^ Binary tensor operator to test + -> Tensor Double -- ^ First operator argument + -> Tensor Double -- ^ Second operator argument + -> Bool +preserveIndicesBinary f t1 t2 = + let i1 = Set.fromList $ indices t1 + i2 = Set.fromList $ indices t2 + in i1 /= i2 || i1 == Set.fromList (indices $ f t1 t2) + +-- | Binary operator other than tensor product cannot contract (or consume) any index +-- | it means, that in operators other than (*), the indices of result tensor are set union of arguments indices +mergeCommonIndices :: + (Tensor Double -> + Tensor Double -> + Tensor Double) -- ^ Binary tensor operator to test + -> Tensor Double -- ^ First operator argument + -> Tensor Double -- ^ Second operator argument + -> Bool +mergeCommonIndices f t1 t2 = + let indices1 = Set.fromList $ indices t1 + indices2 = Set.fromList $ indices t2 + inames1 = Set.fromList $ Index.indexName <$> indices t1 + inames2 = Set.fromList $ Index.indexName <$> indices t2 + + commonIndices = Set.intersection indices1 indices2 + commonIndicesNames = Set.intersection inames1 inames2 + + expectedIndices = Set.union inames1 inames2 + resultIndices = Set.fromList $ Index.indexName <$> indices (f t1 t2) + + -- if we have indices, which have the same name but different type, it is forbidden and test passed + in Set.size commonIndices /= Set.size commonIndicesNames || + -- otherwise, the result indices set must be union of arguments indices + expectedIndices == resultIndices + +-- | Contracted indices have to be consumed in result tensor. +consumeContractedIndices :: + Tensor Double -- ^ first tensor to contract + -> Tensor Double -- ^ second tensor to contract + -> Bool +consumeContractedIndices t1 t2 = + let inames1 = Set.fromList $ Index.indexName <$> indices t1 + inames2 = Set.fromList $ Index.indexName <$> indices t2 + contractedIndices = _contractedIndices t1 t2 + expectedIndices = Set.difference (Set.union inames1 inames2) contractedIndices + resultIndices = Set.fromList $ Index.indexName <$> indices (t1 * t2) + in expectedIndices == resultIndices + +-- | Test generic vector constructor indices +vectorConstructor :: Char -> Positive (Small Int) -> Bool +vectorConstructor c s = + let size = getSmall $ getPositive s + v :: Tensor Double = Vector.fromIndices [c] size fromIntegral + vConst :: Tensor Double = Vector.const [c] size (fromIntegral size) + in v `Multilinear.Class.size` [c] == size && vConst `Multilinear.Class.size` [c] == size + +-- | Test generic form constructor indices +formConstructor :: Char -> Positive (Small Int) -> Bool +formConstructor c s = + let size = getSmall $ getPositive s + f :: Tensor Double = Form.fromIndices [c] size fromIntegral + fConst :: Tensor Double = Form.const [c] size (fromIntegral size) + in f `Multilinear.Class.size` [c] == size && fConst `Multilinear.Class.size` [c] == size + +-- | Test generic matrix constructor indices +matrixConstructor :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +matrixConstructor c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = Matrix.fromIndices [c1,c2] size1 size2 (\x y -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = Matrix.const [c1,c2] size1 size2 (fromIntegral size1) + in c1 == c2 || ( v `Multilinear.Class.size` [c1] == size1 && v `Multilinear.Class.size` [c2] == size2 + && vConst `Multilinear.Class.size` [c1] == size1 && vConst `Multilinear.Class.size` [c2] == size2 ) + +-- | Test generic NVector constructor indices +nVectorConstructor :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +nVectorConstructor c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NVector.fromIndices [c1,c2] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = NVector.const [c1,c2] [size1,size2] (fromIntegral size1) + in c1 == c2 || ( v `Multilinear.Class.size` [c1] == size1 && v `Multilinear.Class.size` [c2] == size2 + && vConst `Multilinear.Class.size` [c1] == size1 && vConst `Multilinear.Class.size` [c2] == size2 ) + +-- | Test generic NForm constructor indices +nFormConstructor :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +nFormConstructor c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NForm.fromIndices [c1,c2] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = NForm.const [c1,c2] [size1,size2] (fromIntegral size1) + in c1 == c2 || ( v `Multilinear.Class.size` [c1] == size1 && v `Multilinear.Class.size` [c2] == size2 + && vConst `Multilinear.Class.size` [c1] == size1 && vConst `Multilinear.Class.size` [c2] == size2 ) + +-- | Test generic vector constructor indices error +vectorConstructorError :: Char -> Positive (Small Int) -> Property +vectorConstructorError c s = + let size = getSmall $ getPositive s + v :: Tensor Double = Vector.fromIndices [c,'a'] size fromIntegral + in expectFailure (total v) + +-- | Test generic form constructor indices error +formConstructorError :: Char -> Positive (Small Int) -> Property +formConstructorError c s = + let size = getSmall $ getPositive s + f :: Tensor Double = Form.fromIndices [c,'a'] size fromIntegral + in expectFailure (total f) + +-- | Test generic matrix constructor indices error +matrixConstructorError :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Property +matrixConstructorError c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = Matrix.fromIndices [c1,c2,'a'] size1 size2 (\x y -> fromIntegral x + fromIntegral y) + in expectFailure (total v) + +-- | Test generic NVector constructor indices error +nVectorConstructorError :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Property +nVectorConstructorError c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NVector.fromIndices [c1,c2,'a'] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + in expectFailure (total v) + +-- | Test generic NForm constructor indices error +nFormConstructorError :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Property +nFormConstructorError c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NForm.fromIndices [c1,c2,'a'] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + in expectFailure (total v) + + +-- | Test generic vector constructor values +vectorConstructorValues :: Char -> Positive (Small Int) -> Bool +vectorConstructorValues c s = + let size = getSmall $ getPositive s + v :: Tensor Double = Vector.fromIndices [c] size fromIntegral + vConst :: Tensor Double = Vector.const [c] size (fromIntegral size) + in all (\i -> v $$| ([c],[i]) == fromIntegral i) [0 .. size - 1] && + all (\i -> vConst $$| ([c],[i]) == fromIntegral size) [0 .. size - 1] + +-- | Test generic form constructor values +formConstructorValues :: Char -> Positive (Small Int) -> Bool +formConstructorValues c s = + let size = getSmall $ getPositive s + f :: Tensor Double = Form.fromIndices [c] size fromIntegral + fConst :: Tensor Double = Form.const [c] size (fromIntegral size) + in all (\i -> f $$| ([c],[i]) == fromIntegral i) [0.. size - 1] && + all (\i -> fConst $$| ([c],[i]) == fromIntegral size) [0 .. size - 1] + +-- | Test generic matrix constructor values +matrixConstructorValues :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +matrixConstructorValues c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = Matrix.fromIndices [c1,c2] size1 size2 (\x y -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = Matrix.const [c1,c2] size1 size2 (fromIntegral size1) + in c1 == c2 || ( + all (\(i1,i2) -> v $$| ([c1,c2],[i1,i2]) == fromIntegral i1 + fromIntegral i2) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) && + all (\(i1,i2) -> vConst $$| ([c1,c2],[i1,i2]) == fromIntegral size1) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) + ) + +-- | Test generic NVector constructor values +nVectorConstructorValues :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +nVectorConstructorValues c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NVector.fromIndices [c1,c2] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = NVector.const [c1,c2] [size1,size2] (fromIntegral size1) + in c1 == c2 || ( + all (\(i1,i2) -> v $$| ([c1,c2],[i1,i2]) == fromIntegral i1 + fromIntegral i2) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) && + all (\(i1,i2) -> vConst $$| ([c1,c2],[i1,i2]) == fromIntegral size1) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) + ) + +-- | Test generic NForm constructor values +nFormConstructorValues :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +nFormConstructorValues c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NForm.fromIndices [c1,c2] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = NForm.const [c1,c2] [size1,size2] (fromIntegral size1) + in c1 == c2 || ( + all (\(i1,i2) -> v $$| ([c1,c2],[i1,i2]) == fromIntegral i1 + fromIntegral i2) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) && + all (\(i1,i2) -> vConst $$| ([c1,c2],[i1,i2]) == fromIntegral size1) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) + ) + +-- | Check indices preservation if zipWith function +zipWithTest :: Tensor Double -> Tensor Double -> Bool +zipWithTest t1@(Scalar _) t2 = preserveIndicesUnary (\t -> Multilinear.Generic.MultiCore.zipWith (+) t1 t) t2 +zipWithTest t1 t2@(Scalar _) = preserveIndicesUnary (\t -> Multilinear.Generic.MultiCore.zipWith (+) t t2) t1 +zipWithTest t1 _ = preserveIndicesBinary (Multilinear.Generic.MultiCore.zipWith (+)) t1 t1 + +-- | Order of the tensor must be equal to number of its covariant and contravariant indices +orderIndices :: Tensor Double -> Bool +orderIndices t = + let (conv, cov) = order t + iConv = Set.fromList $ Index.isContravariant `Prelude.filter` indices t + iCov = Set.fromList $ Index.isCovariant `Prelude.filter` indices t + in conv == Set.size iConv && cov == Set.size iCov + +-- | Tensor must be equivalent in terms of its indices after any index shift +shiftEquiv :: Tensor Double -> Bool +shiftEquiv t = + let inames = indicesNames t + rShiftedTs = (\i -> t |>> i) <$> inames + lShiftedTs = (\i -> t <<| i) <$> inames + rtShiftedTs = (\i -> t |>>> i) <$> inames + ltShiftedTs = (\i -> t <<<| i) <$> inames + allShiftedTs = rShiftedTs ++ lShiftedTs ++ rtShiftedTs ++ ltShiftedTs ++ [t] + allPairs = pure (,) <*> allShiftedTs <*> allShiftedTs + in all (uncurry (|==|)) allPairs + +{-| After rename, index must hold a new name + This property assumes, tensor have max 5 indices of each type -} +renameTest :: Tensor Double -> Bool +renameTest t = + let (conv, cov) = order t + convNs = take conv ['m' .. ] + covNs = take cov ['s' .. ] + renamedT = t $| (convNs, covNs) + inamesAfter = concat $ indicesNames renamedT + in all (\i -> elem i convNs || elem i covNs) inamesAfter + +-- | After any raising or lowering index, it must be a valid type. +raiseLowerTest :: Tensor Double -> Bool +raiseLowerTest t = + let inames = indicesNames t + lowered = inames `zip` ((t `lower`) <$> inames) + raised = inames `zip` ((t `raise`) <$> inames) + isLowered (i,tl) = i `elem` (Index.indexName <$> (Index.isCovariant `Prelude.filter` indices tl)) + isRaised (i,tr) = i `elem` (Index.indexName <$> (Index.isContravariant `Prelude.filter` indices tr)) + in all isLowered lowered && all isRaised raised + +transposeTest :: Tensor Double -> Bool +transposeTest t = + let indices1 = indices t + inames1 = indicesNames t + t2 = transpose t + indices2 = indices t2 + inames2 = indicesNames t2 + in isScalar t || ( inames1 == inames2 && indices1 /= indices2 ) + +-- | Filter second half of elements for each tensor index and check if they disappeared +filterIndexTest :: + Tensor Double -> Bool +filterIndexTest s@(Scalar _) = s == filterIndex "c" (Prelude.const True) s +filterIndexTest t = + let indsT = indices t + -- filter second half of an index + filteredHalf i = filterIndex (Index.indexName i) (< (fromJust (Index.indexSize i) `div` 2)) t + fts = indsT `zip` (filteredHalf <$> indsT) -- tensors with filtered indices, paired with respective transformed indices + in all (\(i,ft) -> + size ft (Index.indexName i) == (fromJust (Index.indexSize i) `div` 2) + ) fts + +-- | Simple show test, just to check if function evaluates at all +showTest :: Tensor Double -> Bool +showTest t = length (show t) > 0 + +-- | ENTRY POINT +main :: IO () +main = do + + -- PRINT PROBABILITY DISTRIBUTION OF TESTED TENSORS ORDER + executePropertyTest "probability distribution of tensors order" 10000 $ + \(t :: Tensor Double) -> collect (order t) $ preserveIndicesUnary abs + executePropertyTest "probability distribution of contracted indices" 10000 $ + \(t1 :: Tensor Double, t2 :: Tensor Double) -> collect (length $ _contractedIndices t1 t2) $ preserveIndicesBinary (+) + + putStrLn "\nTesting multilinear library...\n" + + --------------------------- + -- CHECKING NUM INSTANCE -- + --------------------------- + + executePropertyTest "preserveIndicesBinary for (+)" defTestN $ preserveIndicesBinary (+) + executePropertyTest "preserveIndicesBinary for (-)" defTestN $ preserveIndicesBinary (-) + executePropertyTest "preserveIndicesBinary for (*)" defTestN $ preserveIndicesBinary (*) + executePropertyTest "preserveIndicesUnary for abs" defTestN $ preserveIndicesUnary abs + executePropertyTest "preserveIndicesUnary for signum" defTestN $ preserveIndicesUnary signum + + executePropertyTest "mergeCommonIndices for (+)" defTestN $ mergeCommonIndices (+) + executePropertyTest "mergeCommonIndices for (-)" defTestN $ mergeCommonIndices (-) + executePropertyTest "consumeContractedIndices" defTestN consumeContractedIndices + + ---------------------------------- + -- CHECKING FRACTIONAL INSTANCE -- + ---------------------------------- + + -- TODO + + -------------------------------- + -- CHECKING FLOATING INSTANCE -- + -------------------------------- + + executePropertyTest "preserveIndicesUnary for exp" defTestN $ preserveIndicesUnary exp + executePropertyTest "preserveIndicesUnary for log" defTestN $ preserveIndicesUnary log + executePropertyTest "preserveIndicesUnary for sin" defTestN $ preserveIndicesUnary sin + executePropertyTest "preserveIndicesUnary for cos" defTestN $ preserveIndicesUnary cos + executePropertyTest "preserveIndicesUnary for asin" defTestN $ preserveIndicesUnary asin + executePropertyTest "preserveIndicesUnary for acos" defTestN $ preserveIndicesUnary acos + executePropertyTest "preserveIndicesUnary for atan" defTestN $ preserveIndicesUnary atan + executePropertyTest "preserveIndicesUnary for sinh" defTestN $ preserveIndicesUnary sinh + executePropertyTest "preserveIndicesUnary for cosh" defTestN $ preserveIndicesUnary cosh + executePropertyTest "preserveIndicesUnary for asinh" defTestN $ preserveIndicesUnary asinh + executePropertyTest "preserveIndicesUnary for acosh" defTestN $ preserveIndicesUnary acosh + executePropertyTest "preserveIndicesUnary for atanh" defTestN $ preserveIndicesUnary atanh + + ----------------------------------- + -- CHECKING MULTILINEAR INSTANCE -- + ----------------------------------- + + executePropertyTest "orderIndices" defTestN orderIndices + executePropertyTest "shiftEquiv" defTestN shiftEquiv + executePropertyTest "renamedTest" defTestN renameTest + executePropertyTest "raiseLowerTest" defTestN raiseLowerTest + executePropertyTest "transposeTest" defTestN transposeTest + + ----------------------------------- + -- CHECKING AUXILIRARY FUNCTIONS -- + ----------------------------------- + + executePropertyTest "preserveIndicesUnary for (+.)" defTestN $ preserveIndicesUnary (5 +.) + executePropertyTest "preserveIndicesUnary for (.+)" defTestN $ preserveIndicesUnary (.+ 5) + executePropertyTest "preserveIndicesUnary for (-.)" defTestN $ preserveIndicesUnary (5 -.) + executePropertyTest "preserveIndicesUnary for (.-)" defTestN $ preserveIndicesUnary (.- 5) + executePropertyTest "preserveIndicesUnary for (*.)" defTestN $ preserveIndicesUnary (5 *.) + executePropertyTest "preserveIndicesUnary for (.*)" defTestN $ preserveIndicesUnary (.* 5) + executePropertyTest "filterIndexTest" defTestN filterIndexTest + executePropertyTest "zipWithTest" defTestN zipWithTest + executePropertyTest "showTest" 100 showTest + + ----------------------------------- + -- CHECKING GENERIC CONSTRUCTORS -- + ----------------------------------- + + executePropertyTest "vectorConstructor" defTestN vectorConstructor + executePropertyTest "formConstructor" defTestN formConstructor + executePropertyTest "matrixConstructor" defTestN matrixConstructor + executePropertyTest "nFormConstructor" defTestN nFormConstructor + executePropertyTest "nVectorConstructor" defTestN nVectorConstructor + + executePropertyTest "vectorConstructorError" defTestN vectorConstructorError + executePropertyTest "formConstructorError" defTestN formConstructorError + executePropertyTest "matrixConstructorError" defTestN matrixConstructorError + executePropertyTest "nFormConstructorError" defTestN nFormConstructorError + executePropertyTest "nVectorConstructorError" defTestN nVectorConstructorError + + executePropertyTest "vectorContructorValues" 100 vectorConstructorValues + executePropertyTest "formContructorValues" 100 formConstructorValues + executePropertyTest "matrixConstructorValues" 100 matrixConstructorValues + executePropertyTest "nFormConstructorValues" 100 nFormConstructorValues + executePropertyTest "nVectorConstructorValues" 100 nVectorConstructorValues diff --git a/test/multicore/Test/QuickCheck/Multilinear/Generic/MultiCore.hs b/test/multicore/Test/QuickCheck/Multilinear/Generic/MultiCore.hs new file mode 100644 --- /dev/null +++ b/test/multicore/Test/QuickCheck/Multilinear/Generic/MultiCore.hs @@ -0,0 +1,81 @@ +{-| +Module : Test.QuickCheck.Multilinear.Generic.MultiCore +Description : QucikCheck instances of MultiCore tensor +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Test.QuickCheck.Multilinear.Generic.MultiCore ( + Arbitrary +) where + +import qualified Multilinear.Form as Form +import Multilinear.Generic.MultiCore +import qualified Multilinear.Vector as Vector +import Test.QuickCheck + +-- | Sizes of indices used in Arbitrary Tensor instance +aS :: Int +aS = 12 +bS :: Int +bS = 12 +iS :: Int +iS = 10 +jS :: Int +jS = 15 +kS :: Int +kS = 10 + +-- | Set of three Scalars, that can be used for building more complex tensors +scalars :: [Tensor Double] +scalars = [ + -- Scalars + Scalar (-1.0) + , Scalar 0.0 + , Scalar 1.0 + ] + +-- | Set of 5 simple Scalars and 1D tensors for testing with only upper indices +-- | We use set of a,b,i,j,k indices for further building more complex tensors sets +tensors1DUpper :: [Tensor Double] +tensors1DUpper = [ + -- Vectors with a,b,i,j,k indices + Vector.fromIndices "a" aS $ sin . fromIntegral + , Vector.fromIndices "b" bS $ cos . fromIntegral + , Vector.fromIndices "i" iS $ exp . fromIntegral + , Vector.fromIndices "j" jS $ cosh . fromIntegral + , Vector.fromIndices "k" kS $ tanh . fromIntegral + ] + +-- | Set of 5 simple Scalars and 1D tensors for testing with only lower indices +-- | We use set of a,b,i,j,k indices for further building more complex tensors sets +tensors1DLower :: [Tensor Double] +tensors1DLower = [ + -- Functional with a,b,i,j,k indices - can be contracted with vectors above or matrices below + Form.fromIndices "a" aS $ sin . fromIntegral + , Form.fromIndices "b" bS $ cos . fromIntegral + , Form.fromIndices "i" iS $ exp . fromIntegral + , Form.fromIndices "j" jS $ cosh . fromIntegral + , Form.fromIndices "k" kS $ tanh . fromIntegral + ] + +-- | List sum of scalars and upper and lower indices simple tensors +-- | List contains 18 tensors in total +tensors1D :: [Tensor Double] +tensors1D = scalars ++ tensors1DUpper ++ tensors1DLower + +{-| More complex (up to 3D) tensors, built as sums and differences of tensor products of all pairs from tensors list above + List contains 18^3 = 5832 tensors in total -} + +tensors3D :: [Tensor Double] +tensors3D = pure (*) <*> ts <*> tensors1D + where ts = pure (*) <*> tensors1D <*> tensors1D + +-- | Arbitrary random generating instance of Tensor Double +-- | Simply choose a tensot from tensors list above +instance Arbitrary (Tensor Double) where + arbitrary = elements tensors3D diff --git a/test/sequential/Spec.hs b/test/sequential/Spec.hs new file mode 100644 --- /dev/null +++ b/test/sequential/Spec.hs @@ -0,0 +1,445 @@ +{-| +Module : Main +Description : Test of sequential tensor +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Main ( + main +) where + +import Data.Maybe +import qualified Data.Set as Set +import Multilinear +import qualified Multilinear.Index as Index +import qualified Multilinear.Form as Form +import qualified Multilinear.Matrix as Matrix +import qualified Multilinear.Vector as Vector +import qualified Multilinear.NForm as NForm +import qualified Multilinear.NVector as NVector +import System.Exit +import System.IO +import Test.QuickCheck +import Test.QuickCheck.Multilinear.Generic.Sequential() + +-- | Default test number for property +defTestN :: Int +defTestN = 1000 + + +------------------------------ +-- AUXILIARY TEST FUNCTIONS -- +------------------------------ + + +-- quickCheck with parametrizable tests number +quickCheckN :: Testable prop => Int -> prop -> IO Result +quickCheckN n = quickCheckWithResult (Args + Nothing -- ^ Should we replay a previous test? No. + n -- ^ Maximum number of successful tests before succeeding set to N. + 1 -- ^ Maximum number of discarded tests per successful test before giving up - gave up after first failure. + n -- ^ Size to use for the biggest test cases. + True -- ^ Whether to print anything? yes. + 0) -- ^ Maximum number of shrinks to before giving up. Turn shrinking off. + +-- | Execute property test and check result: +-- | exit test suite with successs code if no errors occured +-- | exit test suite with failure code if any error occured +executePropertyTest :: ( + Testable prop + ) => String -- ^ Tested property name + -> Int -- ^ Number of tests to do + -> prop -- ^ Property to test + -> IO () +executePropertyTest propName n f = do + putStr $ " Checking " ++ propName ++ " " + r <- quickCheckN n f + case r of + Success _ _ _ -> hFlush stdout + _ -> exitFailure + + +---------------------------------------------------- +-- TESTED TENSOR PROPERTIES FOR SEQUENTIAL TENSOR -- +---------------------------------------------------- + + +-- | Unary operator applied on any tensor, +-- | must preserve tensor indices in the result. +preserveIndicesUnary :: + (Tensor Double -> + Tensor Double) -- ^ Unary tensor operator to test + -> Tensor Double -- ^ Operator argument + -> Bool +preserveIndicesUnary f t = Set.fromList (indices t) == Set.fromList (indices (f t)) + +-- | Binary operator applied on any two tensors which have all the same indices, +-- | must preserve set union of these indices in the result. +preserveIndicesBinary :: + (Tensor Double -> + Tensor Double -> + Tensor Double) -- ^ Binary tensor operator to test + -> Tensor Double -- ^ First operator argument + -> Tensor Double -- ^ Second operator argument + -> Bool +preserveIndicesBinary f t1 t2 = + let i1 = Set.fromList $ indices t1 + i2 = Set.fromList $ indices t2 + in i1 /= i2 || i1 == Set.fromList (indices $ f t1 t2) + +-- | Binary operator other than tensor product cannot contract (or consume) any index +-- | it means, that in operators other than (*), the indices of result tensor are set union of arguments indices +mergeCommonIndices :: + (Tensor Double -> + Tensor Double -> + Tensor Double) -- ^ Binary tensor operator to test + -> Tensor Double -- ^ First operator argument + -> Tensor Double -- ^ Second operator argument + -> Bool +mergeCommonIndices f t1 t2 = + let indices1 = Set.fromList $ indices t1 + indices2 = Set.fromList $ indices t2 + inames1 = Set.fromList $ Index.indexName <$> indices t1 + inames2 = Set.fromList $ Index.indexName <$> indices t2 + + commonIndices = Set.intersection indices1 indices2 + commonIndicesNames = Set.intersection inames1 inames2 + + expectedIndices = Set.union inames1 inames2 + resultIndices = Set.fromList $ Index.indexName <$> indices (f t1 t2) + + -- if we have indices, which have the same name but different type, it is forbidden and test passed + in Set.size commonIndices /= Set.size commonIndicesNames || + -- otherwise, the result indices set must be union of arguments indices + expectedIndices == resultIndices + +-- | Contracted indices have to be consumed in result tensor. +consumeContractedIndices :: + Tensor Double -- ^ first tensor to contract + -> Tensor Double -- ^ second tensor to contract + -> Bool +consumeContractedIndices t1 t2 = + let inames1 = Set.fromList $ Index.indexName <$> indices t1 + inames2 = Set.fromList $ Index.indexName <$> indices t2 + contractedIndices = _contractedIndices t1 t2 + expectedIndices = Set.difference (Set.union inames1 inames2) contractedIndices + resultIndices = Set.fromList $ Index.indexName <$> indices (t1 * t2) + in expectedIndices == resultIndices + +-- | Test generic vector constructor indices +vectorConstructor :: Char -> Positive (Small Int) -> Bool +vectorConstructor c s = + let size = getSmall $ getPositive s + v :: Tensor Double = Vector.fromIndices [c] size fromIntegral + vConst :: Tensor Double = Vector.const [c] size (fromIntegral size) + in v `Multilinear.size` [c] == size && vConst `Multilinear.size` [c] == size + +-- | Test generic form constructor indices +formConstructor :: Char -> Positive (Small Int) -> Bool +formConstructor c s = + let size = getSmall $ getPositive s + f :: Tensor Double = Form.fromIndices [c] size fromIntegral + fConst :: Tensor Double = Form.const [c] size (fromIntegral size) + in f `Multilinear.size` [c] == size && fConst `Multilinear.size` [c] == size + +-- | Test generic matrix constructor indices +matrixConstructor :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +matrixConstructor c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = Matrix.fromIndices [c1,c2] size1 size2 (\x y -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = Matrix.const [c1,c2] size1 size2 (fromIntegral size1) + in c1 == c2 || ( v `Multilinear.size` [c1] == size1 && v `Multilinear.size` [c2] == size2 + && vConst `Multilinear.size` [c1] == size1 && vConst `Multilinear.size` [c2] == size2 ) + +-- | Test generic NVector constructor indices +nVectorConstructor :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +nVectorConstructor c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NVector.fromIndices [c1,c2] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = NVector.const [c1,c2] [size1,size2] (fromIntegral size1) + in c1 == c2 || ( v `Multilinear.size` [c1] == size1 && v `Multilinear.size` [c2] == size2 + && vConst `Multilinear.size` [c1] == size1 && vConst `Multilinear.size` [c2] == size2 ) + +-- | Test generic NForm constructor indices +nFormConstructor :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +nFormConstructor c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NForm.fromIndices [c1,c2] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = NForm.const [c1,c2] [size1,size2] (fromIntegral size1) + in c1 == c2 || ( v `Multilinear.size` [c1] == size1 && v `Multilinear.size` [c2] == size2 + && vConst `Multilinear.size` [c1] == size1 && vConst `Multilinear.size` [c2] == size2 ) + +-- | Test generic vector constructor indices error +vectorConstructorError :: Char -> Positive (Small Int) -> Property +vectorConstructorError c s = + let size = getSmall $ getPositive s + v :: Tensor Double = Vector.fromIndices [c,'a'] size fromIntegral + in expectFailure (total v) + +-- | Test generic form constructor indices error +formConstructorError :: Char -> Positive (Small Int) -> Property +formConstructorError c s = + let size = getSmall $ getPositive s + f :: Tensor Double = Form.fromIndices [c,'a'] size fromIntegral + in expectFailure (total f) + +-- | Test generic matrix constructor indices error +matrixConstructorError :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Property +matrixConstructorError c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = Matrix.fromIndices [c1,c2,'a'] size1 size2 (\x y -> fromIntegral x + fromIntegral y) + in expectFailure (total v) + +-- | Test generic NVector constructor indices error +nVectorConstructorError :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Property +nVectorConstructorError c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NVector.fromIndices [c1,c2,'a'] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + in expectFailure (total v) + +-- | Test generic NForm constructor indices error +nFormConstructorError :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Property +nFormConstructorError c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NForm.fromIndices [c1,c2,'a'] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + in expectFailure (total v) + + +-- | Test generic vector constructor values +vectorConstructorValues :: Char -> Positive (Small Int) -> Bool +vectorConstructorValues c s = + let size = getSmall $ getPositive s + v :: Tensor Double = Vector.fromIndices [c] size fromIntegral + vConst :: Tensor Double = Vector.const [c] size (fromIntegral size) + in all (\i -> v $$| ([c],[i]) == fromIntegral i) [0 .. size - 1] && + all (\i -> vConst $$| ([c],[i]) == fromIntegral size) [0 .. size - 1] + +-- | Test generic form constructor values +formConstructorValues :: Char -> Positive (Small Int) -> Bool +formConstructorValues c s = + let size = getSmall $ getPositive s + f :: Tensor Double = Form.fromIndices [c] size fromIntegral + fConst :: Tensor Double = Form.const [c] size (fromIntegral size) + in all (\i -> f $$| ([c],[i]) == fromIntegral i) [0.. size - 1] && + all (\i -> fConst $$| ([c],[i]) == fromIntegral size) [0 .. size - 1] + +-- | Test generic matrix constructor values +matrixConstructorValues :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +matrixConstructorValues c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = Matrix.fromIndices [c1,c2] size1 size2 (\x y -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = Matrix.const [c1,c2] size1 size2 (fromIntegral size1) + in c1 == c2 || ( + all (\(i1,i2) -> v $$| ([c1,c2],[i1,i2]) == fromIntegral i1 + fromIntegral i2) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) && + all (\(i1,i2) -> vConst $$| ([c1,c2],[i1,i2]) == fromIntegral size1) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) + ) + +-- | Test generic NVector constructor values +nVectorConstructorValues :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +nVectorConstructorValues c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NVector.fromIndices [c1,c2] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = NVector.const [c1,c2] [size1,size2] (fromIntegral size1) + in c1 == c2 || ( + all (\(i1,i2) -> v $$| ([c1,c2],[i1,i2]) == fromIntegral i1 + fromIntegral i2) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) && + all (\(i1,i2) -> vConst $$| ([c1,c2],[i1,i2]) == fromIntegral size1) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) + ) + +-- | Test generic NForm constructor values +nFormConstructorValues :: Char -> Char -> Positive (Small Int) -> Positive (Small Int) -> Bool +nFormConstructorValues c1 c2 s1 s2 = + let size1 = getSmall $ getPositive s1 + size2 = getSmall $ getPositive s2 + v :: Tensor Double = NForm.fromIndices [c1,c2] [size1,size2] (\[x,y] -> fromIntegral x + fromIntegral y) + vConst :: Tensor Double = NForm.const [c1,c2] [size1,size2] (fromIntegral size1) + in c1 == c2 || ( + all (\(i1,i2) -> v $$| ([c1,c2],[i1,i2]) == fromIntegral i1 + fromIntegral i2) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) && + all (\(i1,i2) -> vConst $$| ([c1,c2],[i1,i2]) == fromIntegral size1) + (pure (,) <*> [0 .. size1 - 1] <*> [0 .. size2 - 1]) + ) + +-- | Check indices preservation if zipWith function +zipWithTest :: Tensor Double -> Tensor Double -> Bool +zipWithTest t1@(Scalar _) t2 = preserveIndicesUnary (\t -> Multilinear.zipWith (+) t1 t) t2 +zipWithTest t1 t2@(Scalar _) = preserveIndicesUnary (\t -> Multilinear.zipWith (+) t t2) t1 +zipWithTest t1 _ = preserveIndicesBinary (Multilinear.zipWith (+)) t1 t1 + +-- | Order of the tensor must be equal to number of its covariant and contravariant indices +orderIndices :: Tensor Double -> Bool +orderIndices t = + let (conv, cov) = order t + iConv = Set.fromList $ Index.isContravariant `Prelude.filter` indices t + iCov = Set.fromList $ Index.isCovariant `Prelude.filter` indices t + in conv == Set.size iConv && cov == Set.size iCov + +-- | Tensor must be equivalent in terms of its indices after any index shift +shiftEquiv :: Tensor Double -> Bool +shiftEquiv t = + let inames = indicesNames t + rShiftedTs = (\i -> t |>> i) <$> inames + lShiftedTs = (\i -> t <<| i) <$> inames + rtShiftedTs = (\i -> t |>>> i) <$> inames + ltShiftedTs = (\i -> t <<<| i) <$> inames + allShiftedTs = rShiftedTs ++ lShiftedTs ++ rtShiftedTs ++ ltShiftedTs ++ [t] + allPairs = pure (,) <*> allShiftedTs <*> allShiftedTs + in all (uncurry (|==|)) allPairs + +{-| After rename, index must hold a new name + This property assumes, tensor have max 5 indices of each type -} +renameTest :: Tensor Double -> Bool +renameTest t = + let (conv, cov) = order t + convNs = take conv ['m' .. ] + covNs = take cov ['s' .. ] + renamedT = t $| (convNs, covNs) + inamesAfter = concat $ indicesNames renamedT + in all (\i -> elem i convNs || elem i covNs) inamesAfter + +-- | After any raising or lowering index, it must be a valid type. +raiseLowerTest :: Tensor Double -> Bool +raiseLowerTest t = + let inames = indicesNames t + lowered = inames `zip` ((t `lower`) <$> inames) + raised = inames `zip` ((t `raise`) <$> inames) + isLowered (i,tl) = i `elem` (Index.indexName <$> (Index.isCovariant `Prelude.filter` indices tl)) + isRaised (i,tr) = i `elem` (Index.indexName <$> (Index.isContravariant `Prelude.filter` indices tr)) + in all isLowered lowered && all isRaised raised + +transposeTest :: Tensor Double -> Bool +transposeTest t = + let indices1 = indices t + inames1 = indicesNames t + t2 = transpose t + indices2 = indices t2 + inames2 = indicesNames t2 + in isScalar t || ( inames1 == inames2 && indices1 /= indices2 ) + +-- | Filter second half of elements for each tensor index and check if they disappeared +filterIndexTest :: + Tensor Double -> Bool +filterIndexTest s@(Scalar _) = s == filterIndex "c" (Prelude.const True) s +filterIndexTest t = + let indsT = indices t + -- filter second half of an index + filteredHalf i = filterIndex (Index.indexName i) (< (fromJust (Index.indexSize i) `div` 2)) t + fts = indsT `zip` (filteredHalf <$> indsT) -- tensors with filtered indices, paired with respective transformed indices + in all (\(i,ft) -> + size ft (Index.indexName i) == (fromJust (Index.indexSize i) `div` 2) + ) fts + +-- | Simple show test, just to check if function evaluates at all +showTest :: Tensor Double -> Bool +showTest t = length (show t) > 0 + +-- | ENTRY POINT +main :: IO () +main = do + + -- PRINT PROBABILITY DISTRIBUTION OF TESTED TENSORS ORDER + executePropertyTest "probability distribution of tensors order" 10000 $ + \(t :: Tensor Double) -> collect (order t) $ preserveIndicesUnary abs + executePropertyTest "probability distribution of contracted indices" 10000 $ + \(t1 :: Tensor Double, t2 :: Tensor Double) -> collect (length $ _contractedIndices t1 t2) $ preserveIndicesBinary (+) + + putStrLn "\nTesting multilinear library...\n" + + --------------------------- + -- CHECKING NUM INSTANCE -- + --------------------------- + + executePropertyTest "preserveIndicesBinary for (+)" defTestN $ preserveIndicesBinary (+) + executePropertyTest "preserveIndicesBinary for (-)" defTestN $ preserveIndicesBinary (-) + executePropertyTest "preserveIndicesBinary for (*)" defTestN $ preserveIndicesBinary (*) + executePropertyTest "preserveIndicesUnary for abs" defTestN $ preserveIndicesUnary abs + executePropertyTest "preserveIndicesUnary for signum" defTestN $ preserveIndicesUnary signum + + executePropertyTest "mergeCommonIndices for (+)" defTestN $ mergeCommonIndices (+) + executePropertyTest "mergeCommonIndices for (-)" defTestN $ mergeCommonIndices (-) + executePropertyTest "consumeContractedIndices" defTestN consumeContractedIndices + + ---------------------------------- + -- CHECKING FRACTIONAL INSTANCE -- + ---------------------------------- + + -- TODO + + -------------------------------- + -- CHECKING FLOATING INSTANCE -- + -------------------------------- + + executePropertyTest "preserveIndicesUnary for exp" defTestN $ preserveIndicesUnary exp + executePropertyTest "preserveIndicesUnary for log" defTestN $ preserveIndicesUnary log + executePropertyTest "preserveIndicesUnary for sin" defTestN $ preserveIndicesUnary sin + executePropertyTest "preserveIndicesUnary for cos" defTestN $ preserveIndicesUnary cos + executePropertyTest "preserveIndicesUnary for asin" defTestN $ preserveIndicesUnary asin + executePropertyTest "preserveIndicesUnary for acos" defTestN $ preserveIndicesUnary acos + executePropertyTest "preserveIndicesUnary for atan" defTestN $ preserveIndicesUnary atan + executePropertyTest "preserveIndicesUnary for sinh" defTestN $ preserveIndicesUnary sinh + executePropertyTest "preserveIndicesUnary for cosh" defTestN $ preserveIndicesUnary cosh + executePropertyTest "preserveIndicesUnary for asinh" defTestN $ preserveIndicesUnary asinh + executePropertyTest "preserveIndicesUnary for acosh" defTestN $ preserveIndicesUnary acosh + executePropertyTest "preserveIndicesUnary for atanh" defTestN $ preserveIndicesUnary atanh + + ----------------------------------- + -- CHECKING MULTILINEAR INSTANCE -- + ----------------------------------- + + executePropertyTest "orderIndices" defTestN orderIndices + executePropertyTest "shiftEquiv" defTestN shiftEquiv + executePropertyTest "renamedTest" defTestN renameTest + executePropertyTest "raiseLowerTest" defTestN raiseLowerTest + executePropertyTest "transposeTest" defTestN transposeTest + + ----------------------------------- + -- CHECKING AUXILIRARY FUNCTIONS -- + ----------------------------------- + + executePropertyTest "preserveIndicesUnary for (+.)" defTestN $ preserveIndicesUnary (5 +.) + executePropertyTest "preserveIndicesUnary for (.+)" defTestN $ preserveIndicesUnary (.+ 5) + executePropertyTest "preserveIndicesUnary for (-.)" defTestN $ preserveIndicesUnary (5 -.) + executePropertyTest "preserveIndicesUnary for (.-)" defTestN $ preserveIndicesUnary (.- 5) + executePropertyTest "preserveIndicesUnary for (*.)" defTestN $ preserveIndicesUnary (5 *.) + executePropertyTest "preserveIndicesUnary for (.*)" defTestN $ preserveIndicesUnary (.* 5) + executePropertyTest "filterIndexTest" defTestN filterIndexTest + executePropertyTest "zipWithTest" defTestN zipWithTest + executePropertyTest "showTest" 100 showTest + + ----------------------------------- + -- CHECKING GENERIC CONSTRUCTORS -- + ----------------------------------- + + executePropertyTest "vectorConstructor" defTestN vectorConstructor + executePropertyTest "formConstructor" defTestN formConstructor + executePropertyTest "matrixConstructor" defTestN matrixConstructor + executePropertyTest "nFormConstructor" defTestN nFormConstructor + executePropertyTest "nVectorConstructor" defTestN nVectorConstructor + + executePropertyTest "vectorConstructorError" defTestN vectorConstructorError + executePropertyTest "formConstructorError" defTestN formConstructorError + executePropertyTest "matrixConstructorError" defTestN matrixConstructorError + executePropertyTest "nFormConstructorError" defTestN nFormConstructorError + executePropertyTest "nVectorConstructorError" defTestN nVectorConstructorError + + executePropertyTest "vectorContructorValues" 100 vectorConstructorValues + executePropertyTest "formContructorValues" 100 formConstructorValues + executePropertyTest "matrixConstructorValues" 100 matrixConstructorValues + executePropertyTest "nFormConstructorValues" 100 nFormConstructorValues + executePropertyTest "nVectorConstructorValues" 100 nVectorConstructorValues diff --git a/test/sequential/Test/QuickCheck/Multilinear/Generic/Sequential.hs b/test/sequential/Test/QuickCheck/Multilinear/Generic/Sequential.hs new file mode 100644 --- /dev/null +++ b/test/sequential/Test/QuickCheck/Multilinear/Generic/Sequential.hs @@ -0,0 +1,81 @@ +{-| +Module : Test.QuickCheck.Multilinear.Generic.Sequential +Description : QucikCheck instances of sequential Tensor +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Test.QuickCheck.Multilinear.Generic.Sequential ( + Arbitrary +) where + +import qualified Multilinear.Form as Form +import Multilinear.Generic.Sequential +import qualified Multilinear.Vector as Vector +import Test.QuickCheck + +-- | Sizes of indices used in Arbitrary Tensor instance +aS :: Int +aS = 12 +bS :: Int +bS = 12 +iS :: Int +iS = 10 +jS :: Int +jS = 15 +kS :: Int +kS = 10 + +-- | Set of three Scalars, that can be used for building more complex tensors +scalars :: [Tensor Double] +scalars = [ + -- Scalars + Scalar (-1.0) + , Scalar 0.0 + , Scalar 1.0 + ] + +-- | Set of 5 simple Scalars and 1D tensors for testing with only upper indices +-- | We use set of a,b,i,j,k indices for further building more complex tensors sets +tensors1DUpper :: [Tensor Double] +tensors1DUpper = [ + -- Vectors with a,b,i,j,k indices + Vector.fromIndices "a" aS $ sin . fromIntegral + , Vector.fromIndices "b" bS $ cos . fromIntegral + , Vector.fromIndices "i" iS $ exp . fromIntegral + , Vector.fromIndices "j" jS $ cosh . fromIntegral + , Vector.fromIndices "k" kS $ tanh . fromIntegral + ] + +-- | Set of 5 simple Scalars and 1D tensors for testing with only lower indices +-- | We use set of a,b,i,j,k indices for further building more complex tensors sets +tensors1DLower :: [Tensor Double] +tensors1DLower = [ + -- Functional with a,b,i,j,k indices - can be contracted with vectors above or matrices below + Form.fromIndices "a" aS $ sin . fromIntegral + , Form.fromIndices "b" bS $ cos . fromIntegral + , Form.fromIndices "i" iS $ exp . fromIntegral + , Form.fromIndices "j" jS $ cosh . fromIntegral + , Form.fromIndices "k" kS $ tanh . fromIntegral + ] + +-- | List sum of scalars and upper and lower indices simple tensors +-- | List contains 18 tensors in total +tensors1D :: [Tensor Double] +tensors1D = scalars ++ tensors1DUpper ++ tensors1DLower + +{-| More complex (up to 3D) tensors, built as sums and differences of tensor products of all pairs from tensors list above + List contains 18^3 = 5832 tensors in total -} + +tensors3D :: [Tensor Double] +tensors3D = pure (*) <*> ts <*> tensors1D + where ts = pure (*) <*> tensors1D <*> tensors1D + +-- | Arbitrary random generating instance of Tensor Double +-- | Simply choose a tensot from tensors list above +instance Arbitrary (Tensor Double) where + arbitrary = elements tensors3D