multilinear 0.3.2.0 → 0.4.0.0

raw patch · 33 files changed

+3053/−2266 lines, 33 filesdep +paralleldep −mwc-randomdep −primitivedep −statisticssetup-changedPVP ok

version bump matches the API change (PVP)

Dependencies added: parallel

Dependencies removed: mwc-random, primitive, statistics

API changes (from Hackage documentation)

Files

• ChangeLog.md +5/−0
• README.md +2/−1
• Setup.hs +0/−2
• benchmark/memory/Bench.hs +0/−71
• benchmark/multicore/memory/Bench.hs +73/−0
• benchmark/multicore/profile/Bench.hs +27/−0
• benchmark/multicore/time/Bench.hs +46/−0
• benchmark/profile/Bench.hs +0/−27
• benchmark/sequential/memory/Bench.hs +72/−0
• benchmark/sequential/profile/Bench.hs +26/−0
• benchmark/sequential/time/Bench.hs +45/−0
• benchmark/time/Bench.hs +0/−44
• multilinear.cabal +84/−20
• src/Multilinear.hs +6/−43
• src/Multilinear/Class.hs +31/−75
• src/Multilinear/Form.hs +14/−115
• src/Multilinear/Generic.hs +4/−692
• src/Multilinear/Generic/MultiCore.hs +765/−0
• src/Multilinear/Generic/Sequential.hs +736/−0
• src/Multilinear/Index.hs +1/−18
• src/Multilinear/Index/Finite.hs +0/−13
• src/Multilinear/Index/Infinite.hs +0/−71
• src/Multilinear/Matrix.hs +16/−117
• src/Multilinear/NForm.hs +14/−143
• src/Multilinear/NVector.hs +14/−109
• src/Multilinear/Tensor.hs +5/−255
• src/Multilinear/Vector.hs +14/−112
• test/Spec.hs +0/−257
• test/Test/QuickCheck/Multilinear.hs +0/−81
• test/multicore/Spec.hs +446/−0
• test/multicore/Test/QuickCheck/Multilinear/Generic/MultiCore.hs +81/−0
• test/sequential/Spec.hs +445/−0
• test/sequential/Test/QuickCheck/Multilinear/Generic/Sequential.hs +81/−0

@@ -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
 

@@ -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. 

@@ -1,2 +0,0 @@-import Distribution.Simple-main = defaultMain

@@ -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)
-    )

@@ -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)
+    )

@@ -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)

@@ -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]
+    ]

@@ -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)

@@ -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)
+    )

@@ -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)

@@ -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]
+    ]

@@ -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]
-    ]

@@ -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: <https://github.com/ArturB/multilinear#readme> 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

@@ -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

@@ -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

@@ -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
-

@@ -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

@@ -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'

@@ -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'

@@ -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)
-
-

@@ -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 |-}

@@ -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

@@ -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

@@ -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 []

@@ -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 []

@@ -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

@@ -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
-

@@ -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)
-

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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