diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,3 +1,7 @@ +# 0.3.1.0, 2018-11-17 +- added more complex QuickCheck tests +- and therefore, some error fixes + # 0.3.0.0, 2018-11-14 - moved to Vector.Unboxed, to improve performance - simplified error handling - removed separate Err Tensor value diff --git a/benchmark/memory/Bench.hs b/benchmark/memory/Bench.hs --- a/benchmark/memory/Bench.hs +++ b/benchmark/memory/Bench.hs @@ -17,21 +17,55 @@ 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" 1024 1024) gen - func "matrix 10214 x 1024 addition" - (+ Matrix.fromIndices "ab" 1024 1024 gen) - (Matrix.fromIndices "ab" 1024 1024 (\a b -> fromIntegral a + fromIntegral b)) - func "matrix 40 x 60,000 multiplication" - (* Matrix.fromIndices "jk" 60000 40 gen) - (Matrix.fromIndices "ij" 40 60000 gen) + (Matrix.fromIndices "ij" s3 s3) gen + + -- Benchmarking matrix addition + func "matrix 64 x 64 addition" + (+ Matrix.fromIndices "ab" s1 s1 gen) + (Matrix.fromIndices "ab" s1 s1 (\a b -> fromIntegral a + fromIntegral b)) + func "matrix 256 x 256 addition" + (+ Matrix.fromIndices "ab" s2 s2 gen) + (Matrix.fromIndices "ab" s2 s2 (\a b -> fromIntegral a + fromIntegral b)) + func "matrix 1024 x 1024 addition" + (+ Matrix.fromIndices "ab" s3 s3 gen) + (Matrix.fromIndices "ab" s3 s3 (\a b -> fromIntegral a + fromIntegral b)) + + -- Benchmarking matrix multiplication + func "matrix 40 x 4,000 multiplication" + (* Matrix.fromIndices "jk" 4000 40 gen) + (Matrix.fromIndices "ij" 40 4000 gen) + func "matrix 40 x 16,000 multiplication" + (* Matrix.fromIndices "jk" 16000 40 gen) + (Matrix.fromIndices "ij" 40 16000 gen) + func "matrix 40 x 64,000 multiplication" + (* Matrix.fromIndices "jk" 64000 40 gen) + (Matrix.fromIndices "ij" 40 64000 gen) ) diff --git a/benchmark/profile/Bench.hs b/benchmark/profile/Bench.hs new file mode 100644 --- /dev/null +++ b/benchmark/profile/Bench.hs @@ -0,0 +1,27 @@ +{-| +Module : Bench +Description : Benchmark of Multilinear library +Copyright : (c) Artur M. Brodzki, 2018 +License : BSD3 +Maintainer : artur@brodzki.org +Stability : experimental +Portability : Windows/POSIX + +-} + +module Main ( + main +) where + +import Control.DeepSeq +import Multilinear.Class +import qualified Multilinear.Matrix as Matrix +import qualified Multilinear.Vector as Vector + +gen :: Int -> Int -> Double +gen j k = sin (fromIntegral j) + cos (fromIntegral k) + +main :: IO () +main = do + let m = (Matrix.fromIndices "ij" 40 6000 gen) * (Matrix.fromIndices "jk" 6000 40 gen) + m `deepseq` putStrLn $ "All done! Indices of m:" ++ show (indices m) diff --git a/benchmark/time/Bench.hs b/benchmark/time/Bench.hs --- a/benchmark/time/Bench.hs +++ b/benchmark/time/Bench.hs @@ -16,19 +16,27 @@ 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) -sizedMatrixMultBench :: Int -> Benchmark +-- | 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) -sizedMatrixAddBench :: Int -> Benchmark +-- | 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], diff --git a/multilinear.cabal b/multilinear.cabal --- a/multilinear.cabal +++ b/multilinear.cabal @@ -2,10 +2,10 @@ -- -- see: https://github.com/sol/hpack -- --- hash: 9a78adbc87a7e540f962949e44f5f5e34d8cbdfd6a8e7b8b77c7c70a029e6df0 +-- hash: 024729415c87b786fd9fc6cc27a0236cc6844c4e6fe16fa3277095336c73be40 name: multilinear -version: 0.3.0.0 +version: 0.3.1.0 synopsis: Comprehensive and efficient (multi)linear algebra implementation. description: Comprehensive and efficient (multi)linear algebra implementation, based on generic tensor formalism and concise Ricci-Curbastro index syntax. More information available on GitHub: category: Machine learning @@ -60,15 +60,20 @@ type: exitcode-stdio-1.0 main-is: Spec.hs other-modules: + Test.QuickCheck.Multilinear Paths_multilinear hs-source-dirs: test default-extensions: DeriveGeneric FlexibleContexts FlexibleInstances GADTs MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving ghc-options: -O2 -W -threaded -rtsopts -with-rtsopts=-N build-depends: - base >=4.7 && <5 + QuickCheck + , base >=4.7 && <5 + , containers , deepseq + , generic-random , multilinear + , quickcheck-instances default-language: Haskell2010 benchmark memory @@ -84,6 +89,21 @@ base >=4.7 && <5 , multilinear , weigh + default-language: Haskell2010 + +benchmark 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 MultiParamTypeClasses ScopedTypeVariables StandaloneDeriving + ghc-options: -O2 -W -threaded -rtsopts -with-rtsopts=-N + build-depends: + base >=4.7 && <5 + , deepseq + , multilinear default-language: Haskell2010 benchmark time diff --git a/src/Multilinear/Class.hs b/src/Multilinear/Class.hs --- a/src/Multilinear/Class.hs +++ b/src/Multilinear/Class.hs @@ -182,18 +182,6 @@ infixl 5 .* (.*) :: t a -> a -> t a - {-| Tensor adding - functionally equal to Num (+) but more efficient -} - infixl 4 .+. - (.+.) :: t a -> t a -> t a - - {-| Tensor subtracting - functionally equal to Num (-) but more efficient -} - infixl 4 .-. - (.-.) :: t a -> t a -> t a - - {-| Tensor multiplication - functionally equal to Num (*) but more efficient -} - infixl 5 .*. - (.*.) :: t a -> t a -> t a - {-| List of all tensor indices -} indices :: t a -> [TIndex] diff --git a/src/Multilinear/Generic.hs b/src/Multilinear/Generic.hs --- a/src/Multilinear/Generic.hs +++ b/src/Multilinear/Generic.hs @@ -28,13 +28,15 @@ import qualified Multilinear.Index as Index import qualified Multilinear.Index.Finite as Finite -{-| ERROR MESSAGES -} +{-| 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!" @@ -113,10 +115,10 @@ error ("Index + " ++ show ind ++ " out of bonds!") else ts Boxed.! i --- NFData instance +-- | NFData instance instance NFData a => NFData (Tensor a) --- move contravariant indices to lower recursion level +-- | 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 @@ -124,7 +126,7 @@ t <<<| Index.indexName i else t --- Print tensor +-- | Print tensor instance ( Unboxed.Unbox a, Show a, Num a, NFData a ) => Show (Tensor a) where @@ -148,24 +150,6 @@ -- If index is covariant or indifferent, show tensor compoments horizontally _ -> "[" ++ tail (Boxed.foldl' (\string e -> string ++ "," ++ show e) "" ts) ++ "]" --- Tensors can be compared lexigographically --- Allowes to put tensors in typical ordered containers -instance ( - Ord a, Unboxed.Unbox a - ) => Ord (Tensor a) where - - {-# INLINE (<=) #-} - -- Scalar is smaller than any complex tensor - -- Two scalars are compared by they values - Scalar x1 <= Scalar x2 = x1 <= x2 - Scalar _ <= _ = True - _ <= Scalar _ = False - -- Complex tensors are compared lexigographically - SimpleFinite _ ts1 <= SimpleFinite _ ts2 = ts1 <= ts2 - FiniteTensor _ ts1 <= FiniteTensor _ ts2 = ts1 <= ts2 - FiniteTensor _ _ <= SimpleFinite _ _ = False - SimpleFinite _ _ <= FiniteTensor _ _ = True - {-| Merge FiniteTensor of Scalars to SimpleFinite tensor for performance improvement -} {-# INLINE _mergeScalars #-} _mergeScalars :: Unboxed.Unbox a => Tensor a -> Tensor a @@ -180,7 +164,7 @@ _elemByElem' :: (Num a, Unboxed.Unbox a, NFData a) => Tensor a -- ^ First argument of operator -> Tensor a -- ^ Second argument of operator - -> (a -> a -> a) -- ^ Function on tensor elements if indices are different + -> (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 @@ -194,9 +178,8 @@ -- 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)) + | 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 @@ -215,7 +198,7 @@ | 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 -} +{-| 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 @@ -227,11 +210,11 @@ let commonIndices = filter (`Data.List.elem` indicesNames t2) $ indicesNames t1 t1' = foldl' (|>>>) t1 commonIndices t2' = foldl' (|>>>) t2 commonIndices - in t1' `deepseq` t2' `deepseq` _mergeScalars $ _elemByElem' t1' t2' f op + in _mergeScalars $ _elemByElem' t1' t2' f op --- Zipping two tensors with a combinator, assuming they have the same indices +-- | Zipping two tensors with a combinator, assuming they have the same indices. {-# INLINE zipT #-} -zipT :: (Num a, Unboxed.Unbox a, NFData a) +zipT :: (Num a, Unboxed.Unbox a) => (Tensor a -> Tensor a -> Tensor a) -- ^ Two tensors combinator -> (Tensor a -> a -> Tensor a) -- ^ Tensor and scalar combinator -> (a -> Tensor a -> Tensor a) -- ^ Scalar and tensor combinator @@ -244,32 +227,41 @@ zipT _ _ _ f (SimpleFinite index1 v1) (SimpleFinite index2 v2) = if index1 == index2 then SimpleFinite index1 $ Unboxed.zipWith f v1 v2 - else error incompatibleTypes + else zipErr (Index.toTIndex index1) (Index.toTIndex index2) --Two finite tensors case -zipT f _ _ _ (FiniteTensor index1 v1) (FiniteTensor index2 v2) = +zipT f _ _ _ (FiniteTensor index1 v1) (FiniteTensor index2 v2) = if index1 == index2 then FiniteTensor index1 $ Boxed.zipWith f v1 v2 - else error incompatibleTypes + else zipErr (Index.toTIndex index1) (Index.toTIndex index2) -- Finite and simple tensor case -zipT _ f _ _ (FiniteTensor index1 v1) (SimpleFinite index2 v2) = +zipT _ f _ _ (FiniteTensor index1 v1) (SimpleFinite index2 v2) = if index1 == index2 then - FiniteTensor index1 $ - Boxed.generate (Boxed.length v1) (\i -> f (v1 Boxed.! i) (v2 Unboxed.! i)) - else error incompatibleTypes + FiniteTensor index1 $ Boxed.generate (Finite.indexSize index1) (\i -> f (v1 Boxed.! i) (v2 Unboxed.! i)) + else zipErr (Index.toTIndex index1) (Index.toTIndex index2) -- Simple and finite tensor case -zipT _ _ f _ (SimpleFinite index1 v1) (FiniteTensor index2 v2) = +zipT _ _ f _ (SimpleFinite index1 v1) (FiniteTensor index2 v2) = if index1 == index2 then - FiniteTensor index1 $ - Boxed.generate (Unboxed.length v1) (\i -> f (v1 Unboxed.! i) (v2 Boxed.! i)) - else error incompatibleTypes + FiniteTensor index1 $ Boxed.generate (Finite.indexSize index1) (\i -> f (v1 Unboxed.! i) (v2 Boxed.! i)) + else zipErr (Index.toTIndex index1) (Index.toTIndex index2) -- Zipping something with scalar is impossible -zipT _ _ _ _ _ _ = error scalarIndices +zipT _ _ _ _ _ _ = error $ "zipT: " ++ scalarIndices --- dot product of two tensors +-- | zipT error +{-# INLINE zipErr #-} +zipErr :: Index.TIndex -- ^ Index of first dot product parameter + -> Index.TIndex -- ^ Index of second dot product parameter + -> Tensor a -- ^ Erorr message +zipErr i1' i2' = error $ + "zipT: " ++ 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 @@ -281,17 +273,44 @@ | count1 == count2 = Scalar $ Unboxed.sum $ Unboxed.zipWith (*) ts1' ts2' | otherwise = contractionErr (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 (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 (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 (Index.toTIndex i1) (Index.toTIndex i2) +dot t1@(FiniteTensor _ _) t2@(FiniteTensor _ _) = zipT (*) (.*) (*.) (*) t1 t2 --- Other cases cannot happen! +-- 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 (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 (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 (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 (Index.toTIndex i1) (Index.toTIndex i2) +dot t1@(FiniteTensor _ _) t2@(SimpleFinite _ _) = zipT (*) (.*) (*.) (*) t1 t2 + +-- Other cases cannot happen! dot t1' t2' = contractionErr (tensorIndex t1') (tensorIndex t2') --- contraction error ---{-# INLINE contractionErr #-} +-- | contraction error +{-# INLINE contractionErr #-} contractionErr :: Index.TIndex -- ^ Index of first dot product parameter -> Index.TIndex -- ^ Index of second dot product parameter -> Tensor a -- ^ Erorr message @@ -301,7 +320,15 @@ " - index1 is " ++ show i1' ++ " and index2 is " ++ show i2' --- Tensors can be added, subtracted and multiplicated +{-| 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 @@ -329,7 +356,7 @@ {-# INLINE fromInteger #-} fromInteger x = Scalar $ fromInteger x --- Tensors can be divided by each other +-- | Tensors can be divided by each other instance (Unboxed.Unbox a, Fractional a, NFData a) => Fractional (Tensor a) where {-# INLINE (/) #-} @@ -430,18 +457,6 @@ {-# INLINE (*.) #-} x *. t = (x*) `Multilinear.map` t - -- Two tensors sum - {-# INLINE (.+.) #-} - t1 .+. t2 = _elemByElem t1 t2 (+) $ zipT (+) (.+) (+.) (+) - - -- Two tensors difference - {-# INLINE (.-.) #-} - t1 .-. t2 = _elemByElem t1 t2 (-) $ zipT (+) (.+) (+.) (+) - - -- Tensor product - {-# INLINE (.*.) #-} - t1 .*. t2 = _elemByElem t1 t2 (+) dot - -- List of all tensor indices {-# INLINE indices #-} indices x = case x of @@ -551,14 +566,8 @@ then (\un -> Boxed.generate (Unboxed.length un) (\i -> Scalar $ un Unboxed.! i)) <$> (tensorScalars <$> ts1) else tensorsFinite <$> ts1 - -- Convert to list - daneList = Boxed.toList <$> Boxed.toList dane - -- and transpose tensor data (standard function available only for list) - transposedList = Data.List.transpose daneList - -- then reconvert to vector again - transposed = Boxed.fromList <$> Boxed.fromList transposedList - -- and reconstruct tensor with transposed elements - in _mergeScalars $ FiniteTensor index2 $ FiniteTensor index1 <$> transposed + -- 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 diff --git a/src/Multilinear/Index.hs b/src/Multilinear/Index.hs --- a/src/Multilinear/Index.hs +++ b/src/Multilinear/Index.hs @@ -16,6 +16,8 @@ TIndex(..) ) where +import GHC.Generics + {-| Tensor index class which may be lower (covariant), upper (contravariant) or indifferent. -} class Index i where @@ -56,7 +58,7 @@ indexSize :: Maybe Int, tIndexName :: String } - deriving Eq + deriving (Eq, Generic) {-| Show tensor index -} instance Show TIndex where @@ -97,9 +99,11 @@ {-| TIndex must not be converted to TIndex -} toTIndex = id -{-| Indices can be compared by its size |-} +{-| Indices can be compared by its name and size |-} {-| Used to allow to put tensors to typical ordered containers |-} instance Ord TIndex where - ind1 <= ind2 = indexSize ind1 <= indexSize ind2 + ind1 <= ind2 = + tIndexName ind1 <= tIndexName ind2 || + (tIndexName ind1 == tIndexName ind2 && indexSize ind1 <= indexSize ind2) diff --git a/src/Multilinear/Index/Finite.hs b/src/Multilinear/Index/Finite.hs --- a/src/Multilinear/Index/Finite.hs +++ b/src/Multilinear/Index/Finite.hs @@ -76,10 +76,12 @@ toTIndex (Contravariant size name) = TIndex.Contravariant (Just size) name toTIndex (Indifferent size name) = TIndex.Indifferent (Just size) name -{-| Indices can be compared by its size |-} +{-| Indices can be compared by its name and size |-} {-| Used to allow to put tensors to typical ordered containers |-} instance Ord Index where - ind1 <= ind2 = indexSize ind1 <= indexSize ind2 + ind1 <= ind2 = + indexName' ind1 <= indexName' ind2 || + (indexName' ind1 == indexName' ind2 && indexSize ind1 <= indexSize ind2) -- NFData instance instance NFData Index diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -13,44 +13,214 @@ main ) where -import Control.DeepSeq -import Multilinear.Generic +import qualified Data.Set as Set import Multilinear.Class -import qualified Multilinear.Matrix as Matrix -import qualified Multilinear.Vector as Vector +import Multilinear.Generic +import qualified Multilinear.Index as Index +import System.IO +import Test.QuickCheck +import Test.QuickCheck.Multilinear() -v_i :: Tensor Double -v_i = Vector.fromIndices "i" 10 fromIntegral +-- | Default test number for property +defTestN :: Int +defTestN = 1000 -v_j :: Tensor Double -v_j = Vector.fromIndices "j" 5 (\x -> fromIntegral x + 5.0) +-- quickCheck with parametrizable tests number +quickCheckN :: Testable prop => Int -> prop -> IO () +quickCheckN n = quickCheckWith (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. -v_k :: Tensor Double -v_k = Vector.fromIndices "k" 10 fromIntegral +-- Print property test result +printPropertyTest :: ( + Testable prop + ) => String -- ^ Tested property name + -> Int -- ^ Number of tests to do + -> prop -- ^ Property to test + -> IO () +printPropertyTest propName n f = do + putStr $ " Checking " ++ propName ++ " " + quickCheckN n f + hFlush stdout -m_ik :: Tensor Double -m_ik = Matrix.fromIndices "ik" 10 10 (\i j -> fromIntegral i + fromIntegral j) +-- | 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 must merge common indices in result tensor +-- | it means, that in operators other than (*), the result tensor indices 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 `filter` indices t1) + iCovariantNames1 = Set.fromList $ Index.indexName <$> (Index.isCovariant `filter` indices t1) + + iContravariantNames2 = Set.fromList $ Index.indexName <$> (Index.isContravariant `filter` indices t2) + iCovariantNames2 = Set.fromList $ Index.indexName <$> (Index.isCovariant `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 `filter` indices t + iCov = Set.fromList $ Index.isCovariant `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 + allEquivs = uncurry (|==|) <$> allPairs + in False `notElem` allEquivs + +-- | 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 + inamesValid = (\i -> elem i convNs || elem i covNs) <$> inamesAfter + in False `notElem` inamesValid + +-- | 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 `filter` indices tl))) <$> lowered + isRaised = (\(i,tr) -> i `elem` (Index.indexName <$> (Index.isContravariant `filter` indices tr))) <$> raised + in False `notElem` isLowered ++ isRaised + +-- | ENTRY POINT main :: IO () main = do - putStr "v^i = " - print v_i - putStr "v^j = " - print v_j - putStr "v^k = " - print v_k - putStr "Matrix m_ji = v^j + v_i = " - let m = v_j + (v_i \/ "i") - print m - putStr "m_ji * v^i = " - print $ m * v_i - putStr "m_ji * v^k = " - print $ m * v_k - putStr "Matrix m_ik = " - print m_ik - putStr "m_ik * v^k = " - print $ m_ik * v_k - putStr "m_ik |>>> i = " - print $ m_ik |>>> "i" - putStr "m_ji * m_ik" - print $ m * m_ik + --------------------------- + -- CHECKING NUM INSTANCE -- + --------------------------- + + printPropertyTest "preserveIndicesBinary for (+)" defTestN $ preserveIndicesBinary (+) + printPropertyTest "preserveIndicesBinary for (-)" defTestN $ preserveIndicesBinary (-) + printPropertyTest "preserveIndicesBinary for (*)" defTestN $ preserveIndicesBinary (*) + printPropertyTest "preserveIndicesUnary for abs" defTestN $ preserveIndicesUnary abs + printPropertyTest "preserveIndicesUnary for signum" defTestN $ preserveIndicesUnary signum + + printPropertyTest "mergeCommonIndices for (+)" defTestN $ mergeCommonIndices (+) + printPropertyTest "mergeCommonIndices for (-)" defTestN $ mergeCommonIndices (-) + printPropertyTest "consumeContractedIndices" defTestN consumeContractedIndices + + -------------------------------- + -- CHECKING FLOATING INSTANCE -- + -------------------------------- + + printPropertyTest "preserveIndicesUnary for exp" defTestN $ preserveIndicesUnary exp + printPropertyTest "preserveIndicesUnary for log" defTestN $ preserveIndicesUnary log + printPropertyTest "preserveIndicesUnary for sin" defTestN $ preserveIndicesUnary sin + printPropertyTest "preserveIndicesUnary for cos" defTestN $ preserveIndicesUnary cos + printPropertyTest "preserveIndicesUnary for asin" defTestN $ preserveIndicesUnary asin + printPropertyTest "preserveIndicesUnary for acos" defTestN $ preserveIndicesUnary acos + printPropertyTest "preserveIndicesUnary for atan" defTestN $ preserveIndicesUnary atan + printPropertyTest "preserveIndicesUnary for sinh" defTestN $ preserveIndicesUnary sinh + printPropertyTest "preserveIndicesUnary for cosh" defTestN $ preserveIndicesUnary cosh + printPropertyTest "preserveIndicesUnary for asinh" defTestN $ preserveIndicesUnary asinh + printPropertyTest "preserveIndicesUnary for acosh" defTestN $ preserveIndicesUnary acosh + printPropertyTest "preserveIndicesUnary for atanh" defTestN $ preserveIndicesUnary atanh + + ----------------------------------- + -- CHECKING MULTILINEAR INSTANCE -- + ----------------------------------- + + printPropertyTest "preserveIndicesUnary for (+.)" defTestN $ preserveIndicesUnary (5 +.) + printPropertyTest "preserveIndicesUnary for (.+)" defTestN $ preserveIndicesUnary (.+ 5) + printPropertyTest "preserveIndicesUnary for (-.)" defTestN $ preserveIndicesUnary (5 -.) + printPropertyTest "preserveIndicesUnary for (.-)" defTestN $ preserveIndicesUnary (.- 5) + printPropertyTest "preserveIndicesUnary for (*.)" defTestN $ preserveIndicesUnary (5 *.) + printPropertyTest "preserveIndicesUnary for (.*)" defTestN $ preserveIndicesUnary (.* 5) + + printPropertyTest "orderIndices" defTestN orderIndices + printPropertyTest "shiftEquiv" defTestN shiftEquiv + printPropertyTest "renamedTest" defTestN renameTest + printPropertyTest "raiseLowerTest" defTestN raiseLowerTest + diff --git a/test/Test/QuickCheck/Multilinear.hs b/test/Test/QuickCheck/Multilinear.hs new file mode 100644 --- /dev/null +++ b/test/Test/QuickCheck/Multilinear.hs @@ -0,0 +1,104 @@ +{-| +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 Multilinear.Class +import qualified Multilinear.Form as Form +import Multilinear.Generic +import qualified Multilinear.Matrix as Matrix +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 sample tensors for testing. +-- | All vectors, forms and matrices has indices from set [i,j,k] +-- | and sizes compatible with each other, suitable for all (+,-,*) operator. +-- | We have 33 tensors here, which allows to 33^2 circ. 1000 possible test cases for binary operators. +tensors :: [Tensor Double] +tensors = [ + -- Scalars + Scalar (-1.0) + , Scalar 0.0 + , Scalar 1.0 + -- Vectors with i,j,k indices + , Vector.fromIndices "i" iS fromIntegral + , Vector.fromIndices "j" jS (\x -> fromIntegral x - 5.0) + , Vector.fromIndices "k" kS fromIntegral + -- Functional with i,j,k indices - can be contracted with vectors above or matrices below + , Form.fromIndices "i" iS fromIntegral + , Form.fromIndices "j" jS (\x -> fromIntegral x - 5.0) + , Form.fromIndices "k" kS fromIntegral + -- Matrices with a,b indices + , Matrix.fromIndices "ab" aS bS (\i j -> fromIntegral i + fromIntegral j) + , Matrix.fromIndices "ab" aS bS (\i j -> 5 * fromIntegral i - fromIntegral j) + , Matrix.fromIndices "ab" aS bS (\_ _ -> 0.0) + -- The same matrices as above, but with changed indices order + , Matrix.fromIndices "ab" aS bS (\i j -> fromIntegral i + fromIntegral j) |>>> "a" + , Matrix.fromIndices "ab" aS bS (\i j -> 5 * fromIntegral i - fromIntegral j) |>>> "a" + , Matrix.fromIndices "ab" aS bS (\_ _ -> 0.0) |>>> "a" + -- Matrices with i,j,k indices and the same indices sizes as for vectors above + , Matrix.fromIndices "ij" iS jS (\i j -> fromIntegral i + fromIntegral j) + , Matrix.fromIndices "kj" kS jS (\i j -> fromIntegral i + fromIntegral j) + , Matrix.fromIndices "ik" iS kS (\i j -> fromIntegral i + fromIntegral j) + -- The same matrices as above, but with changed indices order + , Matrix.fromIndices "ij" iS jS (\i j -> fromIntegral i + fromIntegral j) |>>> "i" + , Matrix.fromIndices "kj" kS jS (\i j -> fromIntegral i + fromIntegral j) |>>> "k" + , Matrix.fromIndices "ik" iS kS (\i j -> fromIntegral i + fromIntegral j) |>>> "i" + -- Matrices with one i,j,k index and one a,b index + , Matrix.fromIndices "ja" jS aS (\i j -> fromIntegral i + fromIntegral j) + , Matrix.fromIndices "ak" aS kS (\i j -> fromIntegral i + fromIntegral j) + , Matrix.fromIndices "ai" aS iS (\i j -> fromIntegral i + fromIntegral j) + , Matrix.fromIndices "jb" jS bS (\i j -> fromIntegral i + fromIntegral j) + , Matrix.fromIndices "bk" bS kS (\i j -> fromIntegral i + fromIntegral j) + , Matrix.fromIndices "bi" bS iS (\i j -> fromIntegral i + fromIntegral j) + -- The same matrices as above, but with changed indices order + , Matrix.fromIndices "ja" jS aS (\i j -> fromIntegral i + fromIntegral j) |>>> "j" + , Matrix.fromIndices "ak" aS kS (\i j -> fromIntegral i + fromIntegral j) |>>> "a" + , Matrix.fromIndices "ai" aS iS (\i j -> fromIntegral i + fromIntegral j) |>>> "a" + , Matrix.fromIndices "jb" jS bS (\i j -> fromIntegral i + fromIntegral j) |>>> "j" + , Matrix.fromIndices "bk" bS kS (\i j -> fromIntegral i + fromIntegral j) |>>> "b" + , Matrix.fromIndices "bi" bS iS (\i j -> fromIntegral i + fromIntegral j) |>>> "b" + ] + +{- +-- | Second set of tensors; its indices have sizes incompatible with indices of tensors above. +-- | Multiplicating, adding and so on of tensors2 with tensors assumes an error. +tensors2 :: [Tensor Double] +tensors2 = [ + -- Matrices with i,j,k indices and different indices sizes as for vectors above + Matrix.fromIndices "ij" 11 16 (\i j -> fromIntegral i + fromIntegral j) + , Matrix.fromIndices "kj" 11 16 (\i j -> fromIntegral i + fromIntegral j) + , Matrix.fromIndices "ik" 11 11 (\i j -> fromIntegral i + fromIntegral j) + -- The same matrices as above but with changed indices order + , Matrix.fromIndices "ij" 11 16 (\i j -> fromIntegral i + fromIntegral j) |>>> "i" + , Matrix.fromIndices "kj" 11 16 (\i j -> fromIntegral i + fromIntegral j) |>>> "k" + , Matrix.fromIndices "ik" 11 11 (\i j -> fromIntegral i + fromIntegral j) |>>> "i" + ] +-} + +-- | Arbitrary random generating instance of Tensor Double +-- | Simply choose a tensot from tensors list above +instance Arbitrary (Tensor Double) where + arbitrary = elements tensors