diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,3 +2,18 @@ ## [0.1.0.0] - 2019-08-16 Initial release. + +## [0.2] - 2019-08-21 +- removed singletons from dependencies +- removed eigen from dependencies + +- replaced eigen with hmatrix +- added implementation of gaussian elimination based on hmatrix + +- added test suite +- integration with travis ci + +- stack.yaml and cabal.project contain workaround for type-checker + error when using haddock with ghc 8.6 + +- ensured backwards compatibility down to ghc 8.0.2 (stackage lts-9.21) diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,5 +1,7 @@ # The `sparse-tensor` Package +![Hackage-Deps](https://img.shields.io/hackage-deps/v/sparse-tensor) [![Hackage](https://img.shields.io/hackage/v/sparse-tensor)](https://hackage.haskell.org/package/sparse-tensor) [![Build Status](https://travis-ci.org/tobireinhart/sparse-tensor.svg?branch=master)](https://travis-ci.org/tobireinhart/sparse-tensor) + `sparse-tensor` is a Haskell tensor algebra library. It defines the usual tensor algebra functions such as - addition diff --git a/Setup.hs b/Setup.hs deleted file mode 100644 --- a/Setup.hs +++ /dev/null @@ -1,2 +0,0 @@ -import Distribution.Simple -main = defaultMain diff --git a/cabal.project b/cabal.project new file mode 100644 --- /dev/null +++ b/cabal.project @@ -0,0 +1,4 @@ +packages: ./ + +package sparse-tensor + haddock-options: --optghc=-fdefer-type-errors diff --git a/sparse-tensor.cabal b/sparse-tensor.cabal --- a/sparse-tensor.cabal +++ b/sparse-tensor.cabal @@ -1,6 +1,6 @@ name: sparse-tensor -version: 0.1.0.0 -synopsis: typesafe tensor algebra library +version: 0.2 +synopsis: typesafe tensor algebra library description: . This package is intended to be used as a general purpose tensor algebra library. @@ -10,17 +10,6 @@ The implemented tensor data type is capable of being used with an arbitrary number of general abstract indices and can incorporate values of any type that allow for a meaningful addition, scaling, and multiplication. The package is thus very flexible and can easily be customised at wish. - . - This package performs best when the values a given tensor provides are explicitly needed. In particular all incorporated functions are - implemented such that they explicitly manipulate the values of the given tensors. This allows for fast evaluation of the individual values from given tensors. - . - In order to improve memory usage, the tensors in this package employ a sparse storage paradigm. In other words, one may define a tensor by only - providing its non-zero values explicitly. All remaining values are taken to be zero. - . - Furthermore this package also provides the functionality of fully automatically computing all linear independent expressions with given rank and symmetries - that can be composed by using the Minkowski metric and the totally antisymmetric tensor. In order words, the package includes functions that - construct a basis of the space of Lorentz invariant tensors of given symmetry and rank. With slight modifications the therefore implemented algorithms can - also be used for the computation of invariant tensor bases of any other special orthogonal group @SO(p,q)@. homepage: https://github.com/TobiReinhart/sparse-tensor#readme license: MIT license-file: LICENSE @@ -33,6 +22,8 @@ extra-source-files: README.md CHANGELOG.md + stack.yaml + cabal.project source-repository head type: git @@ -41,22 +32,36 @@ library hs-source-dirs: src default-language: Haskell2010 - build-depends: base >= 4.7 && < 5, - containers >= 0.5.7 && < 0.7, + build-depends: base >= 4.9 && < 5, + containers >= 0.5 && < 0.7, tf-random >= 0.5 && < 0.6, - ghc-typelits-natnormalise >= 0.6.2 && < 0.7, - ghc-typelits-knownnat >= 0.5.1 && < 0.7, - eigen == 3.3.4.1, - parallel >= 3.2.2 && < 3.3, - deepseq >= 1.4.2 && < 1.5, - cereal >= 0.5.7 && < 0.6, - singletons >= 2.4.1 && < 2.6, - bytestring >= 0.10.8 && < 0.11, - zlib >= 0.6.2 && < 0.7, - ad >= 4.3 && < 4.4 + ghc-typelits-natnormalise >= 0.5 && < 0.7, + ghc-typelits-knownnat >= 0.2 && < 0.7, + parallel >= 3.2 && < 3.3, + deepseq >= 1.1 && < 1.5, + cereal >= 0.4 && < 0.6, + bytestring >= 0.10 && < 0.11, + zlib >= 0.6 && < 0.7, + ad >= 4.2 && < 4.4, + hmatrix >= 0.16.1 && < 0.21 exposed-modules: Math.Tensor Math.Tensor.LorentzGenerator Math.Tensor.Examples.Gravity Math.Tensor.Examples.Gravity.Schwarzschild Math.Tensor.Examples.Gravity.SchwarzschildSymbolic Math.Tensor.Examples.Gravity.DiffeoSymEqns + Math.Tensor.Internal.LinearAlgebra + +test-suite test-sparse-tensor + type: exitcode-stdio-1.0 + hs-source-dirs: test + default-language: Haskell2010 + main-is: TestMain.hs + other-modules: LinearAlgebra, IndList, Ansatz, Serialization + build-depends: base >= 4.9 && < 5, + hmatrix >= 0.16.1 && < 0.21, + QuickCheck >= 2.8.2 && < 2.14, + tasty >= 0.11 && < 1.3, + tasty-hunit >= 0.9 && < 0.11, + tasty-quickcheck >= 0.8 && < 0.11, + sparse-tensor diff --git a/src/Math/Tensor.hs b/src/Math/Tensor.hs --- a/src/Math/Tensor.hs +++ b/src/Math/Tensor.hs @@ -31,14 +31,12 @@ {-# LANGUAGE GADTs #-} {-# LANGUAGE DataKinds #-} -{-# LANGUAGE KindSignatures #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE TypeSynonymInstances #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE AllowAmbiguousTypes #-} @@ -61,7 +59,7 @@ -- converting them to such it is advantageous when also these lists contain the additional length information in their type. IndList(..), -- -singletonInd, (+>), fromList, fromListMaybe, +singletonInd, (+>), fromList, fromListUnsafe, headInd, tailInd, sortInd, updateInd, -- ** The Tensor Type -- @@ -85,7 +83,7 @@ -- -- ** Index Type Class -- -TIndex(..), +TIndex, -- Ind3(..), Ind9(..), Ind20(..), -- @@ -299,14 +297,14 @@ -- __ Convert a @'Tensor'@ to a sparse matrix assocs list.__ toMatListT1', toMatListT2', toMatListT3', toMatListT4', toMatListT5', toMatListT6', toMatListT7', toMatListT8', -- --- __ Convert a @'Tensor'@ to a sparse "Eigen.SparseMatrix".__ -toEMatrixT1', toEMatrixT2', toEMatrixT3', toEMatrixT4', toEMatrixT5', toEMatrixT6', toEMatrixT7', toEMatrixT8', +-- __ Convert a @'Tensor'@ to a "HM.Matrix".__ +toMatrixT1', toMatrixT2', toMatrixT3', toMatrixT4', toMatrixT5', toMatrixT6', toMatrixT7', toMatrixT8', -- -- | __ Convert all @'Tensors'@ of a heterogeneous tensor List to a combined sparse matrix assocs list.__ toMatListT1, toMatListT2, toMatListT3, toMatListT4, toMatListT5, toMatListT6, toMatListT7, toMatListT8, -- --- | __ Convert all @'Tensors'@ of a heterogeneous tensor List to a combined sparse "Eigen.SparseMatrix".__ -toEMatrixT1, toEMatrixT2, toEMatrixT3, toEMatrixT4, toEMatrixT5, toEMatrixT6, toEMatrixT7, toEMatrixT8, +-- | __ Convert all @'Tensors'@ of a heterogeneous tensor List to a combined "HM.Matrix".__ +toMatrixT1, toMatrixT2, toMatrixT3, toMatrixT4, toMatrixT5, toMatrixT6, toMatrixT7, toMatrixT8, -- -- ** Tensor utility functions -- @@ -324,7 +322,7 @@ -- -- *** Rank Computations -- | Compute the rank of a tensor of @'AnsVarR'@ values. The tensor is converted to a matrix with columns labeling the individual variables that occur in 'ansVarR' --- and rows labeling the independent tensor components. The rank is then computed using 'Eigen' subroutines. These functions is for instance useful +-- and rows labeling the independent tensor components. The rank is then computed using 'hmatrix' subroutines. These functions is for instance useful -- when determining the rank of tensorial equations. -- -- __ Compute the rank of a single @'Tensor'@.__ @@ -343,30 +341,24 @@ import Data.Foldable (toList) import Control.Applicative (liftA2) -import Data.Ratio ((%), numerator, denominator) -import Data.List (nubBy, sortOn, intersect) +import Data.List (sortOn, intersect) import Data.Maybe (mapMaybe, fromMaybe) +import Data.Proxy +import Data.Type.Equality import qualified Data.IntMap.Strict as I -import qualified Data.Map.Strict as M -import Numeric.Natural (Natural(..)) +import Numeric.Natural (Natural) import qualified Numeric.AD.Rank1.Forward as AD import GHC.TypeLits import GHC.Generics (Generic(..)) import Control.DeepSeq (rnf, NFData(..)) import Data.Serialize (encodeLazy, decodeLazy, Serialize(..)) -import Data.Singletons (SingI(..), Proxy(..)) -import Data.Singletons.Decide -import Data.Singletons.Prelude.Enum -import Data.Singletons.TypeLits import Unsafe.Coerce (unsafeCoerce) -import qualified Data.ByteString.Lazy as BS (ByteString(..)) +import qualified Data.ByteString.Lazy as BS (ByteString) import Codec.Compression.GZip (compress, decompress) -import Data.Either (either) --for Linear Algebra subroutines - -import qualified Data.Eigen.SparseMatrix as Sparse -import qualified Data.Eigen.LA as Sol +import qualified Numeric.LinearAlgebra.Data as HM +import Numeric.LinearAlgebra (rank) --Length typed lists for the tensor indices. @@ -378,7 +370,7 @@ -- | Infix synonym for @'Append'@. -(+>) :: (Enum a) => Int -> IndList (n-1) a -> IndList n a +(+>) :: Enum a => Int -> IndList (n-1) a -> IndList n a (+>) i = Append (toEnum i) -- | Construct an @'IndList'@ that only contains a single value. @@ -393,41 +385,42 @@ NonZero :: (1 <= n) => IsZero n deriving instance Show (IsZero n) -isZero :: forall (n :: Nat). SNat n -> IsZero n -isZero n = case n %~ SNat @0 - of Proved Refl -> Zero - Disproved _ -> unsafeCoerce (NonZero @1) +isZero :: forall n. KnownNat n => Proxy n -> IsZero n +isZero _ = case sameNat (Proxy @n) (Proxy @0) + of Nothing -> unsafeCoerce (NonZero @1) + Just Refl -> Zero -fromList' :: forall (n :: Nat). SNat n -> forall (a :: *). [a] -> Maybe (IndList n a) -fromList' n xs = case isZero n - of Zero -> case xs - of [] -> Just Empty - _ -> Nothing - NonZero -> case xs - of [] -> Nothing - x:xs' -> case fromList' (sPred n) xs' - of Just v -> Just (x `Append` v) - Nothing -> Nothing +fromList' :: forall n. KnownNat n => Proxy n -> forall (a :: *). [a] -> Maybe (IndList n a) +fromList' (n :: Proxy n') xs = case isZero n + of Zero -> case xs + of [] -> Just Empty + _ -> Nothing + NonZero -> case xs + of [] -> Nothing + x:xs' -> case fromList' (undefined :: (Proxy (n' - 1))) xs' + of Just v -> Just (x `Append` v) + Nothing -> Nothing --- | Construction of a legth typed @'IndList'@ from an untyped list. +-- | Construction of a length typed @'IndList'@ from an untyped list. -fromList :: forall (n :: Nat). SingI n => forall (a :: *). [a] -> IndList n a -fromList = \case - Just v -> v - Nothing -> undefined - . fromList' sing +fromList :: forall n.KnownNat n => forall (a :: *). [a] -> Maybe (IndList n a) +fromList = fromList' (undefined :: Proxy n) --- | Construction from untyped lists returning a @'Maybe'@ value. +-- | Construction of a length typed @'IndList'@ (partial function). -fromListMaybe :: forall (n :: Nat). SingI n => forall (a :: *). [a] -> Maybe (IndList n a) -fromListMaybe = fromList' sing +fromListUnsafe :: forall n.KnownNat n => forall (a :: *). [a] -> IndList n a +fromListUnsafe xs = fromMaybe (error errorString) (fromList xs) + where + nVal = natVal (undefined :: Proxy n) + len = length xs + errorString = "Could not construct IndList " ++ show nVal ++ " from " ++ show len ++ " elements." --Instances of the length typed lists. instance (KnownNat n, Generic a) => Generic (IndList n a) where type Rep (IndList n a) = Rep [a] - to r = fromList $ to r + to r = fromListUnsafe $ to r from = from . toList deriving instance (KnownNat n, Generic a, Serialize a) => Serialize (IndList n a) @@ -443,11 +436,11 @@ deriving instance (Show a) => Show (IndList n a) instance Functor (IndList n) where - fmap f Empty = Empty + fmap _ Empty = Empty fmap f (Append x xs) = Append (f x) (fmap f xs) instance Foldable (IndList n) where - foldr f y Empty = y + foldr _ y Empty = y foldr f y (Append x xs) = f x (foldr f y xs) --list operations on length typed lists @@ -464,35 +457,30 @@ sortInd Empty = Empty sortInd (Append x xs) = insertSorted x $ sortInd xs -toListInd :: IndList n a -> [a] -toListInd = toList - -combineInds :: IndList n a -> IndList m a -> IndList (n+m) a -combineInds Empty l = l -combineInds (Append x xs) l = Append x $ combineInds xs l - -- | An implementation of the usual head function from "Data.List" for @'IndList'@s. The function returns the first element of the @'IndList'@. headInd :: IndList n a -> a -headInd (Append x xs) = x +headInd Empty = error "head of empty IndList" +headInd (Append x _) = x -- | An implementation of the usual tail function from "Data.List" for @'IndList'@s. The function removes the first element of the @'IndList'@. tailInd :: IndList n a -> IndList (n-1) a -tailInd (Append x xs) = xs +tailInd Empty = error "tail of empty IndList" +tailInd (Append _ xs) = xs -- | An implementation of the usual @(!!)@ function from "Data.List" for @'IndList'@s. The function returns the list element that is specified by the provided integer value. -- Indices start from zero. indexInd :: Int -> IndList n a -> a -indexInd 0 (Append x xs) = x -indexInd i (Append x xs) = indexInd (i-1) xs +indexInd 0 (Append x _) = x +indexInd i (Append _ xs) = indexInd (i-1) xs indexInd _ _ = error "Index is too large!" -- | The function replaces the element at index/position specified by its first argument with the element that is specified by its second argument. updateInd :: Int -> a -> IndList n a -> IndList n a -updateInd 0 s (Append x xs) = Append s xs +updateInd 0 s (Append _ xs) = Append s xs updateInd i s (Append x xs) = Append x $ updateInd (i-1) s xs updateInd _ _ _ = error "Index is too large!" @@ -504,12 +492,14 @@ where val = indexInd (i-1) xs newL = updateInd (i-1) x xs +swapHead _ Empty = error "cannot swap head of empty list" removeContractionInd :: (Eq a) => Int -> a -> (IndList n a, c) -> Maybe (IndList (n-1) a,c) removeContractionInd 0 ind1 (Append x xs, t) | ind1 == x = Just (xs,t) | otherwise = Nothing removeContractionInd i ind1 (Append x xs,t) = (\(m,n) -> (Append x m, n)) <$> removeContractionInd (i-1) ind1 (xs,t) +removeContractionInd _ _ (Empty, _) = error "cannot remove index from empty index list" {-- resort Inds in IndList according to the permutation given by [Int], length of [Int] must be n @@ -520,13 +510,13 @@ => [C, A, B] --} -resortInd :: (SingI n, Ord a) => [Int] -> IndList n a -> IndList n a +resortInd :: (KnownNat n, Ord a) => [Int] -> IndList n a -> IndList n a resortInd perm indList = newindList where l' = toList indList l'' = if length l' == length perm then zip perm l' else error "permutation has wrong length" lReSorted = sortOn fst l'' - newindList = fromList $ map snd lReSorted + newindList = fromListUnsafe $ map snd lReSorted --Index type class @@ -592,7 +582,7 @@ subS a b = a `addS` negateS b -- | Newtype wrapper that is used for representing scalar values. -newtype SField a = SField a deriving (Show, Eq, Ord) +newtype SField a = SField a deriving (Show, Eq, Ord, Generic, Serialize) instance Functor SField where fmap f (SField a) = SField $ f a @@ -612,7 +602,7 @@ -- | Newtype wrapper for symbolic values. -- Note that this version does not yet support simplification of symbolic values. -newtype SSymbolic = SSymbolic String deriving (Show, Eq, Ord) +newtype SSymbolic = SSymbolic String deriving (Show, Eq, Ord, Generic, Serialize) instance Num SSymbolic where SSymbolic s1 + SSymbolic s2 = SSymbolic $ "(" ++ s1 ++ ")+(" ++ s2 ++ ")" @@ -705,7 +695,7 @@ -- @'AnsVar'@ values describes a linear equation system. Using the functions @'toMatListT1'@, ... this equation system can -- then be transformed into a matrix with columns labeling the individual @'AnsVar'@s. -newtype AnsVar a = AnsVar (I.IntMap a) deriving Show +newtype AnsVar a = AnsVar (I.IntMap a) deriving (Show, Eq, Generic, Serialize) type AnsVarR = AnsVar (SField Rational) @@ -810,13 +800,9 @@ -- All future functions maintain this order when acting on a valid, i.e. ordered @'TMap'@s. type TMap k v = [(k,v)] ---check if the list is sorted -isValidTMap :: (Ord k, Eq v) => TMap k v -> Bool -isValidTMap l = l == sortOn fst l - --insert a new key value pair in the list, if already present the function (v -> v -> v) determines the new value at the ind insertWithTMap :: (Ord k) => (v -> v -> v) -> k -> v -> TMap k v -> TMap k v -insertWithTMap f key val [] = [(key,val)] +insertWithTMap _ key val [] = [(key,val)] insertWithTMap f key val ((k1,v1):xs) | key < k1 = (key,val) : ((k1,v1):xs) | key == k1 = (k1,f val v1) : xs @@ -824,8 +810,8 @@ --combine two sorted lists with combiner function addTMaps :: (Ord k) => (v -> v -> v) -> TMap k v -> TMap k v -> TMap k v -addTMaps f m1 [] = m1 -addTMaps f [] m2 = m2 +addTMaps _ m1 [] = m1 +addTMaps _ [] m2 = m2 addTMaps f ((k1,v1):xs) ((k2,v2):ys) | k1 < k2 = (k1,v1) : addTMaps f xs ((k2,v2):ys) | k2 < k1 = (k2,v2) : addTMaps f ((k1,v1):xs) ys @@ -879,36 +865,36 @@ --convert between typesafe and non typesafe tensors -lemma :: forall n m. (n-1 :~: m) -> (m+1 :~: n) -lemma _ = unsafeCoerce (Refl @n) - toRep :: forall n k v. KnownNat n => Tensor n k v -> TensorRep k v toRep (Scalar v) = ScalarR v -toRep (Tensor m) = -- case isZero (SNat @n) n == 0 is not possible, because types - -- of Zero -> undefined - -- NonZero -> - case lemma @n Refl - of Refl -> - let r = fromIntegral $ GHC.TypeLits.natVal (Proxy @n) - in TensorR r $ mapTMap (\(t :: Tensor (n-1) k v) -> toRep t) m +toRep (Tensor m) = case isZero (Proxy @n) of + NonZero -> let r = fromIntegral $ GHC.TypeLits.natVal (Proxy @n) + in TensorR r $ mapTMap (\(t :: Tensor (n-1) k v) -> toRep t) m toRep ZeroTensor = let r = fromIntegral $ GHC.TypeLits.natVal (Proxy @n) in ZeroR r -fromRep :: forall n k v. KnownNat n => TensorRep k v -> Tensor n k v -fromRep (ScalarR v) = case isZero (SNat @n) - of Zero -> Scalar v - NonZero -> undefined -fromRep (TensorR r m) = case someNatVal (fromIntegral r) - of Just l -> case l - of SomeNat (_ :: Proxy x) -> case isZero (SNat @x) - of NonZero -> case sameNat (Proxy @x) (Proxy @n) - of Nothing -> undefined - Just Refl -> Tensor (mapTMap (\t -> fromRep t :: Tensor (x-1) k v) m) - Zero -> undefined - Nothing -> undefined +fromRep :: forall n k v. KnownNat n => TensorRep k v -> Maybe (Tensor n k v) +fromRep (ScalarR v) = case isZero (Proxy @n) + of Zero -> Just (Scalar v) + NonZero -> Nothing +fromRep (TensorR r m) = case someNatVal (fromIntegral r) of + Nothing -> Nothing + Just l -> case l of + SomeNat (x' :: Proxy x) -> + case isZero x' of + Zero -> Nothing + NonZero -> case sameNat x' (Proxy @n) of + Nothing -> Nothing + Just Refl -> let tMap' = mapTMap (\t -> fromRep t :: Maybe (Tensor (x-1) k v)) m + tMap'' = filterTMap (\case Nothing -> False + _ -> True) tMap' + tMap''' = mapTMap (\(Just t) -> t) tMap'' + in case tMap''' of + [] -> Nothing + _ -> Just $ Tensor tMap''' fromRep (ZeroR r) = case someNatVal (fromIntegral r) - of Just l -> ZeroTensor - Nothing -> undefined + of Just _ -> Just ZeroTensor + Nothing -> Nothing --instances for the tensor data type @@ -920,22 +906,25 @@ instance KnownNat n => Generic (Tensor n k v) where type Rep (Tensor n k v) = Rep (TensorRep k v) from = from . toRep - to = fromRep . to + to = \case + Nothing -> error "Could not reconstruct tensor from representation." + Just t -> t + . fromRep . to deriving instance (KnownNat n, Ord k, Serialize k, Serialize v) => Serialize (Tensor n k v) instance Functor (Tensor n k) where fmap f (Scalar x) = Scalar (f x) fmap f (Tensor m) = Tensor (mapTMap (fmap f) m) - fmap f ZeroTensor = ZeroTensor + fmap _ ZeroTensor = ZeroTensor deriving instance (Show a, Show k) => Show (Tensor n k a) deriving instance (Eq a, Eq k) => Eq (Tensor n k a) -getTensorMap :: Tensor (n+1) k v -> TMap k (Tensor n k v) +getTensorMap :: forall n k v.(1 <= n) => Tensor n k v -> TMap k (Tensor (n-1) k v) getTensorMap (Tensor m) = m - +getTensorMap ZeroTensor = [] toMatListT1' :: (TIndex k1, TAdd a) => AbsTensor1 n1 k1 (AnsVar a) -> [((Int,Int),a)] toMatListT1' = collectMatList . toMatList1' @@ -961,29 +950,29 @@ toMatListT8' :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, TAdd a) => AbsTensor8 n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 (AnsVar a) -> [((Int,Int),a)] toMatListT8' = collectMatList . toMatList8' -toEMatrixT1' :: (TIndex k1, Real a) => AbsTensor1 n1 k1 (AnsVar (SField a)) -> Sparse.SparseMatrixXd -toEMatrixT1' = assocsToSparse . toMatListT1' +toMatrixT1' :: (TIndex k1, Real a) => AbsTensor1 n1 k1 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT1' = assocsToMatrix . toMatListT1' -toEMatrixT2' :: (TIndex k1, Real a) => AbsTensor2 n1 n2 k1 (AnsVar (SField a)) -> Sparse.SparseMatrixXd -toEMatrixT2' = assocsToSparse . toMatListT2' +toMatrixT2' :: (TIndex k1, Real a) => AbsTensor2 n1 n2 k1 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT2' = assocsToMatrix . toMatListT2' -toEMatrixT3' :: (TIndex k1, TIndex k2, Real a) => AbsTensor3 n1 n2 n3 k1 k2 (AnsVar (SField a)) -> Sparse.SparseMatrixXd -toEMatrixT3' = assocsToSparse . toMatListT3' +toMatrixT3' :: (TIndex k1, TIndex k2, Real a) => AbsTensor3 n1 n2 n3 k1 k2 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT3' = assocsToMatrix . toMatListT3' -toEMatrixT4' :: (TIndex k1, TIndex k2, Real a) => AbsTensor4 n1 n2 n3 n4 k1 k2 (AnsVar (SField a)) -> Sparse.SparseMatrixXd -toEMatrixT4' = assocsToSparse . toMatListT4' +toMatrixT4' :: (TIndex k1, TIndex k2, Real a) => AbsTensor4 n1 n2 n3 n4 k1 k2 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT4' = assocsToMatrix . toMatListT4' -toEMatrixT5' :: (TIndex k1, TIndex k2, TIndex k3, Real a) => AbsTensor5 n1 n2 n3 n4 n5 k1 k2 k3 (AnsVar (SField a)) -> Sparse.SparseMatrixXd -toEMatrixT5' = assocsToSparse . toMatListT5' +toMatrixT5' :: (TIndex k1, TIndex k2, TIndex k3, Real a) => AbsTensor5 n1 n2 n3 n4 n5 k1 k2 k3 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT5' = assocsToMatrix . toMatListT5' -toEMatrixT6' :: (TIndex k1, TIndex k2, TIndex k3, Real a) => AbsTensor6 n1 n2 n3 n4 n5 n6 k1 k2 k3 (AnsVar (SField a)) -> Sparse.SparseMatrixXd -toEMatrixT6' = assocsToSparse . toMatListT6' +toMatrixT6' :: (TIndex k1, TIndex k2, TIndex k3, Real a) => AbsTensor6 n1 n2 n3 n4 n5 n6 k1 k2 k3 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT6' = assocsToMatrix . toMatListT6' -toEMatrixT7' :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, Real a) => AbsTensor7 n1 n2 n3 n4 n5 n6 n7 k1 k2 k3 k4 (AnsVar (SField a)) -> Sparse.SparseMatrixXd -toEMatrixT7' = assocsToSparse . toMatListT7' +toMatrixT7' :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, Real a) => AbsTensor7 n1 n2 n3 n4 n5 n6 n7 k1 k2 k3 k4 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT7' = assocsToMatrix . toMatListT7' -toEMatrixT8' :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, Real a) => AbsTensor8 n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 (AnsVar (SField a)) -> Sparse.SparseMatrixXd -toEMatrixT8' = assocsToSparse . toMatListT8' +toMatrixT8' :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, Real a) => AbsTensor8 n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT8' = assocsToMatrix . toMatListT8' @@ -1002,15 +991,15 @@ -- | Same functionality as @'fromListT'@ but the indices of the various values are specified by usual list. This is certainly more flexible -- however all at the cost of reduced type safety compared to @'fromListT'@. -fromListT' :: forall n k v b. (TIndex k, TAdd v, SingI n) => [([k],v)] -> Tensor n k v +fromListT' :: forall n k v. (TIndex k, TAdd v, KnownNat n) => [([k],v)] -> Tensor n k v fromListT' l = fromListT indList where - indList = map (\(x,y) -> (fromList x, y)) l + indList = map (\(x,y) -> (fromListUnsafe x, y)) l fromListT1 :: (TIndex k1, TAdd v) => [(IndTuple1 n1 k1, v)] -> AbsTensor1 n1 k1 v fromListT1 = fromListT -fromListT1' :: forall n1 k1 v b. (SingI n1, TIndex k1, TAdd v) => [([k1],v)] -> AbsTensor1 n1 k1 v +fromListT1' :: forall n1 k1 v. (KnownNat n1, TIndex k1, TAdd v) => [([k1],v)] -> AbsTensor1 n1 k1 v fromListT1' = fromListT' fromListT2 :: (TIndex k1, TAdd v) => [(IndTuple2 n1 n2 k1, v)] -> AbsTensor2 n1 n2 k1 v @@ -1018,76 +1007,77 @@ where tensList = map mkTens2 l -fromListT2' :: forall n1 n2 k1 v b. (SingI n1, SingI n2, TIndex k1, TAdd v) => [(([k1],[k1]),v)] -> AbsTensor2 n1 n2 k1 v +fromListT2' :: forall n1 n2 k1 v. (KnownNat n1, KnownNat n2, TIndex k1, TAdd v) => [(([k1],[k1]),v)] -> AbsTensor2 n1 n2 k1 v fromListT2' l = fromListT2 indList where - indList = map (\((x1,x2),y) -> ((fromList x1, fromList x2),y)) l + indList = map (\((x1,x2),y) -> ((fromListUnsafe x1, fromListUnsafe x2),y)) l fromListT3 :: (TIndex k1, TIndex k2, TAdd v) => [(IndTuple3 n1 n2 n3 k1 k2, v)] -> AbsTensor3 n1 n2 n3 k1 k2 v fromListT3 l = foldr (&+) ZeroTensor tensList where tensList = map mkTens3 l -fromListT3' :: forall n1 n2 n3 k1 k2 v b. (SingI n1, SingI n2, SingI n3, TIndex k1, TIndex k2, TAdd v) => [(([k1],[k1],[k2]),v)] -> AbsTensor3 n1 n2 n3 k1 k2 v +fromListT3' :: forall n1 n2 n3 k1 k2 v. (KnownNat n1, KnownNat n2, KnownNat n3, TIndex k1, TIndex k2, TAdd v) => [(([k1],[k1],[k2]),v)] -> AbsTensor3 n1 n2 n3 k1 k2 v fromListT3' l = fromListT3 indList where - indList = map (\((x1,x2,x3),y) -> ((fromList x1, fromList x2, fromList x3),y)) l + indList = map (\((x1,x2,x3),y) -> ((fromListUnsafe x1, fromListUnsafe x2, fromListUnsafe x3),y)) l fromListT4 :: (TIndex k1, TIndex k2, TAdd v) => [(IndTuple4 n1 n2 n3 n4 k1 k2, v)] -> AbsTensor4 n1 n2 n3 n4 k1 k2 v fromListT4 l = foldr (&+) ZeroTensor tensList where tensList = map mkTens4 l -fromListT4' :: forall n1 n2 n3 n4 k1 k2 v b. (SingI n1, SingI n2, SingI n3, SingI n4, TIndex k1, TIndex k2, TAdd v) => [(([k1],[k1],[k2],[k2]),v)] -> AbsTensor4 n1 n2 n3 n4 k1 k2 v +fromListT4' :: forall n1 n2 n3 n4 k1 k2 v. (KnownNat n1, KnownNat n2, KnownNat n3, KnownNat n4, TIndex k1, TIndex k2, TAdd v) => [(([k1],[k1],[k2],[k2]),v)] -> AbsTensor4 n1 n2 n3 n4 k1 k2 v fromListT4' l = fromListT4 indList where - indList = map (\((x1,x2,x3,x4),y) -> ((fromList x1, fromList x2, fromList x3, fromList x4),y)) l + indList = map (\((x1,x2,x3,x4),y) -> ((fromListUnsafe x1, fromListUnsafe x2, fromListUnsafe x3, fromListUnsafe x4),y)) l fromListT5 :: (TIndex k1, TIndex k2, TIndex k3, TAdd v) => [(IndTuple5 n1 n2 n3 n4 n5 k1 k2 k3, v)] -> AbsTensor5 n1 n2 n3 n4 n5 k1 k2 k3 v fromListT5 l = foldr (&+) ZeroTensor tensList where tensList = map mkTens5 l -fromListT5' :: forall n1 n2 n3 n4 n5 k1 k2 k3 v b. (SingI n1, SingI n2, SingI n3, SingI n4, SingI n5, TIndex k1, TIndex k2, TIndex k3, TAdd v) => [(([k1],[k1],[k2],[k2],[k3]),v)] -> AbsTensor5 n1 n2 n3 n4 n5 k1 k2 k3 v +fromListT5' :: forall n1 n2 n3 n4 n5 k1 k2 k3 v. (KnownNat n1, KnownNat n2, KnownNat n3, KnownNat n4, KnownNat n5, TIndex k1, TIndex k2, TIndex k3, TAdd v) => [(([k1],[k1],[k2],[k2],[k3]),v)] -> AbsTensor5 n1 n2 n3 n4 n5 k1 k2 k3 v fromListT5' l = fromListT5 indList where - indList = map (\((x1,x2,x3,x4,x5),y) -> ((fromList x1, fromList x2, fromList x3, fromList x4, fromList x5),y)) l + indList = map (\((x1,x2,x3,x4,x5),y) -> ((fromListUnsafe x1, fromListUnsafe x2, fromListUnsafe x3, fromListUnsafe x4, fromListUnsafe x5),y)) l fromListT6 :: (TIndex k1, TIndex k2, TIndex k3, TAdd v) => [(IndTuple6 n1 n2 n3 n4 n5 n6 k1 k2 k3, v)] -> AbsTensor6 n1 n2 n3 n4 n5 n6 k1 k2 k3 v fromListT6 l = foldr (&+) ZeroTensor tensList where tensList = map mkTens6 l -fromListT6' :: forall n1 n2 n3 n4 n5 n6 k1 k2 k3 v b. (SingI n1, SingI n2, SingI n3, SingI n4, SingI n5, SingI n6, TIndex k1, TIndex k2, TIndex k3, TAdd v) => [(([k1],[k1],[k2],[k2],[k3],[k3]),v)] -> AbsTensor6 n1 n2 n3 n4 n5 n6 k1 k2 k3 v +fromListT6' :: forall n1 n2 n3 n4 n5 n6 k1 k2 k3 v. (KnownNat n1, KnownNat n2, KnownNat n3, KnownNat n4, KnownNat n5, KnownNat n6, TIndex k1, TIndex k2, TIndex k3, TAdd v) => [(([k1],[k1],[k2],[k2],[k3],[k3]),v)] -> AbsTensor6 n1 n2 n3 n4 n5 n6 k1 k2 k3 v fromListT6' l = fromListT6 indList where - indList = map (\((x1,x2,x3,x4,x5,x6),y) -> ((fromList x1, fromList x2, fromList x3, fromList x4, fromList x5, fromList x6),y)) l + indList = map (\((x1,x2,x3,x4,x5,x6),y) -> ((fromListUnsafe x1, fromListUnsafe x2, fromListUnsafe x3, fromListUnsafe x4, fromListUnsafe x5, fromListUnsafe x6),y)) l fromListT7 :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, TAdd v) => [(IndTuple7 n1 n2 n3 n4 n5 n6 n7 k1 k2 k3 k4, v)] -> AbsTensor7 n1 n2 n3 n4 n5 n6 n7 k1 k2 k3 k4 v fromListT7 l = foldr (&+) ZeroTensor tensList where tensList = map mkTens7 l -fromListT7' :: forall n1 n2 n3 n4 n5 n6 n7 k1 k2 k3 k4 v b. (SingI n1, SingI n2, SingI n3, SingI n4, SingI n5, SingI n6, SingI n7, TIndex k1, TIndex k2, TIndex k3, TIndex k4, TAdd v) => [(([k1],[k1],[k2],[k2],[k3],[k3],[k4]),v)] -> AbsTensor7 n1 n2 n3 n4 n5 n6 n7 k1 k2 k3 k4 v +fromListT7' :: forall n1 n2 n3 n4 n5 n6 n7 k1 k2 k3 k4 v. (KnownNat n1, KnownNat n2, KnownNat n3, KnownNat n4, KnownNat n5, KnownNat n6, KnownNat n7, TIndex k1, TIndex k2, TIndex k3, TIndex k4, TAdd v) => [(([k1],[k1],[k2],[k2],[k3],[k3],[k4]),v)] -> AbsTensor7 n1 n2 n3 n4 n5 n6 n7 k1 k2 k3 k4 v fromListT7' l = fromListT7 indList where - indList = map (\((x1,x2,x3,x4,x5,x6,x7),y) -> ((fromList x1, fromList x2, fromList x3, fromList x4, fromList x5, fromList x6, fromList x7),y)) l + indList = map (\((x1,x2,x3,x4,x5,x6,x7),y) -> ((fromListUnsafe x1, fromListUnsafe x2, fromListUnsafe x3, fromListUnsafe x4, fromListUnsafe x5, fromListUnsafe x6, fromListUnsafe x7),y)) l fromListT8 :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, TAdd v) => [(IndTuple8 n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4, v)] -> AbsTensor8 n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 v fromListT8 l = foldr (&+) ZeroTensor tensList where tensList = map mkTens8 l -fromListT8' :: forall n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 v b. (SingI n1, SingI n2, SingI n3, SingI n4, SingI n5, SingI n6, SingI n7, SingI n8, TIndex k1, TIndex k2, TIndex k3, TIndex k4, TAdd v) => [(([k1],[k1],[k2],[k2],[k3],[k3],[k4],[k4]),v)] -> AbsTensor8 n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 v +fromListT8' :: forall n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 v. (KnownNat n1, KnownNat n2, KnownNat n3, KnownNat n4, KnownNat n5, KnownNat n6, KnownNat n7, KnownNat n8, TIndex k1, TIndex k2, TIndex k3, TIndex k4, TAdd v) => [(([k1],[k1],[k2],[k2],[k3],[k3],[k4],[k4]),v)] -> AbsTensor8 n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 v fromListT8' l = fromListT8 indList where - indList = map (\((x1,x2,x3,x4,x5,x6,x7,x8),y) -> ((fromList x1, fromList x2, fromList x3, fromList x4, fromList x5, fromList x6, fromList x7, fromList x8),y)) l + indList = map (\((x1,x2,x3,x4,x5,x6,x7,x8),y) -> ((fromListUnsafe x1, fromListUnsafe x2, fromListUnsafe x3, fromListUnsafe x4, fromListUnsafe x5, fromListUnsafe x6, fromListUnsafe x7, fromListUnsafe x8),y)) l --helper function for tensor addition: adding one indices value pair to the tree structure of a given tensor insertOrAdd :: (TIndex k, TAdd v) => (IndList n k, v) -> Tensor n k v -> Tensor n k v insertOrAdd (Empty, a) (Scalar b) = Scalar $ a `addS` b +insertOrAdd (_, _) (Scalar _) = error "cannot add key value pair to scalar" insertOrAdd (Append x xs, a) (Tensor m) = Tensor $ insertWithTMap (\_ o -> insertOrAdd (xs, a) o) x indTens m where indTens = mkTens (xs, a) @@ -1105,6 +1095,7 @@ (&+) (Tensor m1) (Tensor m2) = Tensor $ addTMaps (&+) m1 m2 (&+) t1 ZeroTensor = t1 (&+) ZeroTensor t2 = t2 +(&+) _ _ = error "incompatible combination of summands" -- | Scalar multiplication of an arbitrary @'Tensor'@. Only requirement is that the corresponding values and the scalar type satisfy the -- @'Prod'@ constraint. @@ -1130,6 +1121,7 @@ (&-) (Tensor m1) (Tensor m2) = Tensor $ addTMaps (&-) m1 m2 (&-) t1 ZeroTensor = t1 (&-) ZeroTensor t2 = negateS t2 +(&-) _ _ = error "incompatible combination of summands" -- | Tensor product of two @'Tensor'@s. In the result for each index type the indices of the first @'Tensor'@ are arranged left of those in the second @'Tensor'@ argument. -- The values of the two @'Tensor'@s must satisfy the @'Prod'@ constraint. @@ -1140,8 +1132,8 @@ (&*) (Scalar x) (Scalar y) = Scalar $ prod x y (&*) (Scalar x) t2 = fmap (prod x) t2 (&*) (Tensor m) t2 = Tensor $ mapTMap (&* t2) m -(&*) t1 ZeroTensor = ZeroTensor -(&*) ZeroTensor t2 = ZeroTensor +(&*) _ ZeroTensor = ZeroTensor +(&*) ZeroTensor _ = ZeroTensor --encode and decode tensors as bytestrings for efficient storing @@ -1150,8 +1142,8 @@ encodeTensor = compress . encodeLazy -- | Utility function to decompress and de-serialize a @'Data.ByteString.Lazy.ByteString'@ into a @'Tensor'@, making use of the @'Serialize'@ instance. -decodeTensor :: (KnownNat n, Ord k, Serialize k, Serialize v) => BS.ByteString -> Tensor n k v -decodeTensor bs = either error id $ decodeLazy $ decompress bs +decodeTensor :: (KnownNat n, Ord k, Serialize k, Serialize v) => BS.ByteString -> Either String (Tensor n k v) +decodeTensor bs = decodeLazy $ decompress bs -- | The function computes partial derivatives of spacetime tensors of type @'STTens'@ with @'CFun'@ values. @@ -1281,7 +1273,8 @@ where l = map (\(x,y) -> (swapHead j x, y)) $ toListT t tensorTrans (i, j) (Tensor m) = Tensor $ mapTMap (tensorTrans (i-1, j-1)) m -tensorTrans (i ,j) ZeroTensor = ZeroTensor +tensorTrans _ ZeroTensor = ZeroTensor +tensorTrans _ (Scalar _) = error "attempt to transpose using a non-existing index" -- | > tensorTrans1 = tensorTrans tensorTrans1 :: (TIndex k1, TAdd v) => (Int,Int) -> AbsTensor1 n1 k1 v -> AbsTensor1 n1 k1 v @@ -1400,7 +1393,7 @@ -- | Evaluate a @'Tensor'@ for a specific value of its first contravariant index type returning the corresponding sub @'Tensor'@. -- The additional functions specified below apply the evaluation function @'evalTens'@ to the deeper @'Tensor'@ levels. -evalTens :: (KnownNat (n+1), TIndex k, TAdd v) => Int -> k -> Tensor (n+1) k v -> Tensor n k v +evalTens :: (KnownNat n, 1<=n, TIndex k, TAdd v) => Int -> k -> Tensor n k v -> Tensor (n-1) k v evalTens ind indVal (Tensor m) | ind > size -1 || ind < 0 = error "wrong index to evaluate" | ind == 0 = fromMaybe ZeroTensor $ lookup indVal m @@ -1409,6 +1402,7 @@ size = length $ fst $ head $ toListT' (Tensor m) l = [1..ind] ++ 0 : [ind+1..size -1] newTens = resortTens l (Tensor m) +evalTens _ _ ZeroTensor = ZeroTensor -- | > evalTens1 = evalTens evalTens1 :: (KnownNat n1, TIndex k1, TAdd v) => Int -> k1 -> AbsTensor1 (n1+1) k1 v -> AbsTensor1 n1 k1 v @@ -1510,9 +1504,9 @@ -- | Same functionality as @'symTens'@ but including the \( \frac{1}{2} \) in the result and thus defining a projection. -- The following functions apply @'symTensFac'@ to the index types of the deeper leaf levels, i.e. covariant indices of the 1st index type, contravariant indices of the 2nd index type, etc. -- --- > symTensFac inds t = (SField $ 1%2) &. symTens inds t +-- > symTensFac inds t = SField (1/2 :: Rational) &. symTens inds t symTensFac :: (TIndex k, TAdd v, Prod (SField Rational) v) => (Int,Int) -> Tensor n k v -> Tensor n k (TProd (SField Rational) v) -symTensFac inds t = (SField $ 1%2) &. symTens inds t +symTensFac inds t = SField (1/2 :: Rational) &. symTens inds t -- | > symATensFac1 = symTensFac symATensFac1 :: (TIndex k1, TAdd v, Prod (SField Rational) v) => @@ -1638,9 +1632,9 @@ -- | Same functionality as @'aSymTens'@ but including the \( \frac{1}{2} \) factors in the result and thus defining a projection. -- The following functions apply @'aSymTensFac'@ to the index types of the deeper leaf levels, i.e. covariant indices of the 1st index type, contravariant indices of the 2nd index type, etc. -- --- > aSymTensFac inds t = (SField $ 1%2)&. aSymTens inds t +-- > aSymTensFac inds t = SField (1/2 :: Rational) &. aSymTens inds t aSymTensFac :: (TIndex k, TAdd v, Prod (SField Rational) v) => (Int,Int) -> Tensor n k v -> Tensor n k (TProd (SField Rational) v) -aSymTensFac inds t = (SField $ 1%2)&. aSymTens inds t +aSymTensFac inds t = SField (1/2 :: Rational) &. aSymTens inds t -- | > aSymATensFac1 = aSymTensFac aSymATensFac1 :: (TIndex k1, TAdd v, Prod (SField Rational) v) => @@ -1766,9 +1760,9 @@ -- | Same functionality as @'symBlockTens'@ but including the usual factors of \( \frac{1}{2} \) in the result and thus defining a projection. -- The following functions apply @'symBlockTensFac'@ to the index types of the deeper leaf levels, i.e. covariant indices of the 1st index type, contravariant indices of the 2nd index type, etc. -- --- > symBlockTensFac inds t = (SField $ 1%2)&. symBlockTens inds t +-- > symBlockTensFac inds t = SField (1/2 :: Rational) &. symBlockTens inds t symBlockTensFac :: (TIndex k, TAdd v, Prod (SField Rational) v) => ([Int],[Int]) -> Tensor n k v -> Tensor n k (TProd (SField Rational) v) -symBlockTensFac inds t = (SField $ 1%2)&. symBlockTens inds t +symBlockTensFac inds t = SField (1/2 :: Rational) &. symBlockTens inds t -- | > symBlockATensFac1 = symBlockTensFac symBlockATensFac1 :: (TIndex k1, TAdd v, Prod (SField Rational) v) => @@ -1894,9 +1888,9 @@ -- | Same functionality as @'aSymBlockTens'@ but including the usual factors of \( \frac{1}{2} \) in the result and thus defining a projection. -- The following functions apply @'aSymBlockTensFac'@ to the index types of the deeper leaf levels, i.e. covariant indices of the 1st index type, contravariant indices of the 2nd index type, etc. -- --- > aSymBlockTensFac inds t = (SField $ 1%2)&. aSymBlockTens inds t +-- > aSymBlockTensFac inds t = SField (1/2 :: Rational) &. aSymBlockTens inds t aSymBlockTensFac :: (TIndex k, TAdd v, Prod (SField Rational) v) => ([Int],[Int]) -> Tensor n k v -> Tensor n k (TProd (SField Rational) v) -aSymBlockTensFac inds t = (SField $ 1%2)&. aSymBlockTens inds t +aSymBlockTensFac inds t = SField (1/2 :: Rational) &. aSymBlockTens inds t -- | > aSymBlockATensFac1 = aSymBlockTensFac aSymBlockATensFac1 :: (TIndex k1, TAdd v, Prod (SField Rational) v) => @@ -1963,6 +1957,7 @@ l = zip (repeat x) xs lNew = map pure l l' = getAllSwaps xs +getAllSwaps _ = error "cannot get swaps from empty permutations" getAllBlockSwaps :: [[Int]] -> [[([Int],[Int])]] getAllBlockSwaps [x,y] = [[(x,y)]] @@ -1971,6 +1966,7 @@ l = zip (repeat x) xs lNew = map pure l l' = getAllBlockSwaps xs +getAllBlockSwaps _ = error "cannot get swaps from empty permutations" factorial :: Int -> Int factorial 1 = 1 @@ -2053,7 +2049,7 @@ cyclicSymTensFac :: (TIndex k, TAdd v, Prod (SField Rational) v) => [Int] -> Tensor n k v -> Tensor n k (TProd (SField Rational) v) cyclicSymTensFac inds t = fac &. cyclicSymTens inds t where - fac = SField $ 1 % fromIntegral (factorial $ length inds) + fac = SField (1 / fromIntegral (factorial $ length inds) :: Rational) -- | > cyclicSymATensFac1 = cyclicSymTensFac cyclicSymATensFac1 :: (TIndex k1, TAdd v, Prod (SField Rational) v) => @@ -2120,7 +2116,7 @@ cyclicASymTens inds t = newTens where swapList = getAllSwaps inds - signList = map (\x -> (-1) ^ length x) swapList + signList = map (\x -> (-1 :: Integer) ^ length x) swapList tensList' = map (foldr tensorTrans t) swapList tensList = zipWith (\s v -> if s == -1 then negateS v else v) signList tensList' newTens = foldr (&+) t tensList @@ -2190,7 +2186,7 @@ cyclicASymTensFac :: (TIndex k, TAdd v, Prod (SField Rational) v) => [Int] -> Tensor n k v -> Tensor n k (TProd (SField Rational) v) cyclicASymTensFac inds t = fac &. cyclicASymTens inds t where - fac = SField $ 1 % fromIntegral (factorial $ length inds) + fac = SField (1 / fromIntegral (factorial $ length inds) :: Rational) -- | > cyclicASymATensFac1 = cyclicASymTensFac cyclicASymATensFac1 :: (TIndex k1, TAdd v, Prod (SField Rational) v) => @@ -2325,7 +2321,7 @@ cyclicBlockSymTensFac :: (TIndex k, TAdd v, Prod (SField Rational) v) => [[Int]] -> Tensor n k v -> Tensor n k (TProd (SField Rational) v) cyclicBlockSymTensFac inds t = fac &. cyclicBlockSymTens inds t where - fac = SField $ 1 % fromIntegral (factorial $ length inds) + fac = SField (1 / fromIntegral (factorial $ length inds) :: Rational) -- | > cyclicBlockSymATensFac1 = cyclicBlockSymTensFac cyclicBlockSymATensFac1 :: (TIndex k1, TAdd v, Prod (SField Rational) v) => @@ -2394,8 +2390,8 @@ l3 = filter (\(_,y) -> not (null y)) l2 tensList = map (\(x,y) -> (x, fromListT y)) l3 tensorContr (i,j) (Tensor m) = Tensor $ mapTMap (tensorContr (i-1,j)) m -tensorContr inds ZeroTensor = ZeroTensor -tensorContr inds (Scalar s) = error "cannot contract scalar!" +tensorContr _ ZeroTensor = ZeroTensor +tensorContr _ (Scalar _) = error "cannot contract scalar!" -- | @'contrATens1'@ is a synonym for @'tensorContr'@. It applies the contraction to the 1st index type of a tensor. -- @@ -2441,8 +2437,8 @@ --construct a tensor with a single value -mkTens1 :: (IndTuple1 n1 k1, v) -> AbsTensor1 n1 k1 v -mkTens1 = mkTens +_mkTens1 :: (IndTuple1 n1 k1, v) -> AbsTensor1 n1 k1 v +_mkTens1 = mkTens mkTens2 :: (IndTuple2 n1 n2 k1, v) -> AbsTensor2 n1 n2 k1 v mkTens2 ((i1,i2),s) = mkTens (i1,mkTens (i2,s)) @@ -2521,12 +2517,25 @@ concatMap (\(x,y) -> appendT2 x $ toListT y ) $ concatMap (\(x,y) -> appendT1 x $ toListT y ) $ toListT t +appendT1 :: a1 -> [(b1, b2)] -> [((a1, b1), b2)] appendT1 i = map (\(x,y) -> ((i,x),y)) + +appendT2 :: (a1, a2) -> [(b1, b2)] -> [((a1, a2, b1), b2)] appendT2 (i1,i2) = map (\(x,y) -> ((i1,i2,x),y)) + +appendT3 :: (a1, a2, a3) -> [(b1, b2)] -> [((a1, a2, a3, b1), b2)] appendT3 (i1,i2,i3) = map (\(x,y) -> ((i1,i2,i3,x),y)) + +appendT4 :: (a1, a2, a3, a4) -> [(b1, b2)] -> [((a1, a2, a3, a4, b1), b2)] appendT4 (i1,i2,i3,i4) = map (\(x,y) -> ((i1,i2,i3,i4,x),y)) + +appendT5 :: (a1, a2, a3, a4, a5) -> [(b1, b2)] -> [((a1, a2, a3, a4, a5, b1), b2)] appendT5 (i1,i2,i3,i4,i5) = map (\(x,y) -> ((i1,i2,i3,i4,i5,x),y)) + +appendT6 :: (a1, a2, a3, a4, a5, a6) -> [(b1, b2)] -> [((a1, a2, a3, a4, a5, a6, b1), b2)] appendT6 (i1,i2,i3,i4,i5,i6) = map (\(x,y) -> ((i1,i2,i3,i4,i5,i6,x),y)) + +appendT7 :: (a1, a2, a3, a4, a5, a6, a7) -> [(b1, b2)] -> [((a1, a2, a3, a4, a5, a6, a7, b1), b2)] appendT7 (i1,i2,i3,i4,i5,i6,i7) = map (\(x,y) -> ((i1,i2,i3,i4,i5,i6,i7,x),y)) -- | This function converts a given @'Tensor'@ to a non-typed index tuple list. @@ -2723,9 +2732,9 @@ toMatList8' :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, TAdd a) => AbsTensor8 n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 (AnsVar a) -> [[(Int, a)]] toMatList8' t = map snd $ toListShowVar8 t -normalize :: [(Int,Rational)] -> ([(Int,Rational)],Rational) -normalize [] = ([],1) -normalize ((a,b) : xs) = ((a,1) : map (\(x,y) -> (x,y / b)) xs,b) +_normalize :: [(Int,Rational)] -> ([(Int,Rational)],Rational) +_normalize [] = ([],1) +_normalize ((a,b) : xs) = ((a,1) : map (\(x,y) -> (x,y / b)) xs,b) data TensList1 k1 v where EmptyTList1 :: TensList1 k1 v @@ -2762,35 +2771,35 @@ -- | Usual map function for heterogeneous tensor lists. mapTensList1 :: (forall n1. AbsTensor1 n1 k1 v -> b ) -> TensList1 k1 v -> [b] -mapTensList1 f EmptyTList1 = [] +mapTensList1 _ EmptyTList1 = [] mapTensList1 f (AppendTList1 t l) = f t : mapTensList1 f l mapTensList2 :: (forall n1 n2. AbsTensor2 n1 n2 k1 v -> b ) -> TensList2 k1 v -> [b] -mapTensList2 f EmptyTList2 = [] +mapTensList2 _ EmptyTList2 = [] mapTensList2 f (AppendTList2 t l) = f t : mapTensList2 f l mapTensList3 :: (forall n1 n2 n3. AbsTensor3 n1 n2 n3 k1 k2 v -> b ) -> TensList3 k1 k2 v -> [b] -mapTensList3 f EmptyTList3 = [] +mapTensList3 _ EmptyTList3 = [] mapTensList3 f (AppendTList3 t l) = f t : mapTensList3 f l mapTensList4 :: (forall n1 n2 n3 n4. AbsTensor4 n1 n2 n3 n4 k1 k2 v -> b ) -> TensList4 k1 k2 v -> [b] -mapTensList4 f EmptyTList4 = [] +mapTensList4 _ EmptyTList4 = [] mapTensList4 f (AppendTList4 t l) = f t : mapTensList4 f l mapTensList5 :: (forall n1 n2 n3 n4 n5. AbsTensor5 n1 n2 n3 n4 n5 k1 k2 k3 v -> b ) -> TensList5 k1 k2 k3 v -> [b] -mapTensList5 f EmptyTList5 = [] +mapTensList5 _ EmptyTList5 = [] mapTensList5 f (AppendTList5 t l) = f t : mapTensList5 f l mapTensList6 :: (forall n1 n2 n3 n4 n5 n6. AbsTensor6 n1 n2 n3 n4 n5 n6 k1 k2 k3 v -> b ) -> TensList6 k1 k2 k3 v -> [b] -mapTensList6 f EmptyTList6 = [] +mapTensList6 _ EmptyTList6 = [] mapTensList6 f (AppendTList6 t l) = f t : mapTensList6 f l mapTensList7 :: (forall n1 n2 n3 n4 n5 n6 n7. AbsTensor7 n1 n2 n3 n4 n5 n6 n7 k1 k2 k3 k4 v -> b ) -> TensList7 k1 k2 k3 k4 v -> [b] -mapTensList7 f EmptyTList7 = [] +mapTensList7 _ EmptyTList7 = [] mapTensList7 f (AppendTList7 t l) = f t : mapTensList7 f l mapTensList8 :: (forall n1 n2 n3 n4 n5 n6 n7 n8. AbsTensor8 n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 v -> b ) -> TensList8 k1 k2 k3 k4 v -> [b] -mapTensList8 f EmptyTList8 = [] +mapTensList8 _ EmptyTList8 = [] mapTensList8 f (AppendTList8 t l) = f t : mapTensList8 f l infixr 5 ...> @@ -2981,7 +2990,7 @@ matList = concat $ mapTensList8 toMatList8' t ---convert to Eigen format for using LA subroutines +--convert to hmatrix format for using LA subroutines dims :: [((Int, Int), a)] -> (Int, Int) dims xs = (rows, cols) @@ -2989,92 +2998,90 @@ rows = maximum $ map (fst.fst) xs cols = maximum $ map (snd.fst) xs -assocsToSparse :: Real a => [((Int, Int), SField a)] -> Sparse.SparseMatrixXd -assocsToSparse [] = Sparse.fromList 1 1 [(0,0,0)] -assocsToSparse assocs = Sparse.fromList rows cols els +assocsToMatrix :: Real a => [((Int, Int), SField a)] -> HM.Matrix Double +assocsToMatrix [] = (1 HM.>< 1) [0] +assocsToMatrix assocs = HM.assoc (rows, cols) 0 els where (rows, cols) = dims assocs - els = map (\((x, y), SField z) -> (x-1, y-1, fromRational $ toRational z)) assocs + els = map (\((x, y), SField z) -> ((x-1, y-1), fromRational $ toRational z)) assocs -toEMatrixT1 :: (TIndex k1) => TensList1 k1 (AnsVar (SField Rational)) -> Sparse.SparseMatrixXd -toEMatrixT1 = assocsToSparse . toMatListT1 +toMatrixT1 :: (TIndex k1, Real a) => TensList1 k1 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT1 = assocsToMatrix . toMatListT1 -toEMatrixT2 :: (TIndex k1) => TensList2 k1 (AnsVar (SField Rational)) -> Sparse.SparseMatrixXd -toEMatrixT2 = assocsToSparse . toMatListT2 +toMatrixT2 :: (TIndex k1, Real a) => TensList2 k1 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT2 = assocsToMatrix . toMatListT2 -toEMatrixT3 :: (TIndex k1, TIndex k2) => TensList3 k1 k2 (AnsVar (SField Rational)) -> Sparse.SparseMatrixXd -toEMatrixT3 = assocsToSparse . toMatListT3 +toMatrixT3 :: (TIndex k1, TIndex k2, Real a) => TensList3 k1 k2 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT3 = assocsToMatrix . toMatListT3 -toEMatrixT4 :: (TIndex k1, TIndex k2) => TensList4 k1 k2 (AnsVar (SField Rational)) -> Sparse.SparseMatrixXd -toEMatrixT4 = assocsToSparse . toMatListT4 +toMatrixT4 :: (TIndex k1, TIndex k2, Real a) => TensList4 k1 k2 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT4 = assocsToMatrix . toMatListT4 -toEMatrixT5 :: (TIndex k1, TIndex k2, TIndex k3) => TensList5 k1 k2 k3 (AnsVar (SField Rational)) -> Sparse.SparseMatrixXd -toEMatrixT5 = assocsToSparse . toMatListT5 +toMatrixT5 :: (TIndex k1, TIndex k2, TIndex k3, Real a) => TensList5 k1 k2 k3 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT5 = assocsToMatrix . toMatListT5 -toEMatrixT6 :: (TIndex k1, TIndex k2, TIndex k3) => TensList6 k1 k2 k3 (AnsVar (SField Rational)) -> Sparse.SparseMatrixXd -toEMatrixT6 = assocsToSparse . toMatListT6 +toMatrixT6 :: (TIndex k1, TIndex k2, TIndex k3, Real a) => TensList6 k1 k2 k3 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT6 = assocsToMatrix . toMatListT6 -toEMatrixT7 :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4) => TensList7 k1 k2 k3 k4 (AnsVar (SField Rational)) -> Sparse.SparseMatrixXd -toEMatrixT7 = assocsToSparse . toMatListT7 +toMatrixT7 :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, Real a) => TensList7 k1 k2 k3 k4 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT7 = assocsToMatrix . toMatListT7 -toEMatrixT8 :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4) => TensList8 k1 k2 k3 k4 (AnsVar (SField Rational)) -> Sparse.SparseMatrixXd -toEMatrixT8 = assocsToSparse . toMatListT8 +toMatrixT8 :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, Real a) => TensList8 k1 k2 k3 k4 (AnsVar (SField a)) -> HM.Matrix Double +toMatrixT8 = assocsToMatrix . toMatListT8 --rank of the tensor can be computed with rank Sol.FullPivLU or Sol.JakobiSVD -tensorRank1' :: (TIndex k1) => AbsTensor1 n1 k1 (AnsVar (SField Rational)) -> Int -tensorRank1' t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT1 (singletonTList1 t) +tensorRank1' :: (TIndex k1, Real a, Real a) => AbsTensor1 n1 k1 (AnsVar (SField a)) -> Int +tensorRank1' t = rank $ toMatrixT1 (singletonTList1 t) -tensorRank1 :: (TIndex k1) => TensList1 k1 (AnsVar (SField Rational)) -> Int -tensorRank1 t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT1 t +tensorRank1 :: (TIndex k1, Real a) => TensList1 k1 (AnsVar (SField a)) -> Int +tensorRank1 t = rank $ toMatrixT1 t -tensorRank2' :: (TIndex k1) => AbsTensor2 n1 n2 k1 (AnsVar (SField Rational)) -> Int -tensorRank2' t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT2 (singletonTList2 t) +tensorRank2' :: (TIndex k1, Real a) => AbsTensor2 n1 n2 k1 (AnsVar (SField a)) -> Int +tensorRank2' t = rank $ toMatrixT2 (singletonTList2 t) -tensorRank2 :: (TIndex k1) => TensList2 k1 (AnsVar (SField Rational)) -> Int -tensorRank2 t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT2 t +tensorRank2 :: (TIndex k1, Real a) => TensList2 k1 (AnsVar (SField a)) -> Int +tensorRank2 t = rank $ toMatrixT2 t -tensorRank3' :: (TIndex k1, TIndex k2) => AbsTensor3 n1 n2 n3 k1 k2 (AnsVar (SField Rational)) -> Int -tensorRank3' t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT3 (singletonTList3 t) +tensorRank3' :: (TIndex k1, TIndex k2, Real a) => AbsTensor3 n1 n2 n3 k1 k2 (AnsVar (SField a)) -> Int +tensorRank3' t = rank $ toMatrixT3 (singletonTList3 t) -tensorRank3 :: (TIndex k1, TIndex k2) => TensList3 k1 k2 (AnsVar (SField Rational)) -> Int -tensorRank3 t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT3 t +tensorRank3 :: (TIndex k1, TIndex k2, Real a) => TensList3 k1 k2 (AnsVar (SField a)) -> Int +tensorRank3 t = rank $ toMatrixT3 t -tensorRank4' :: (TIndex k1, TIndex k2) => AbsTensor4 n1 n2 n3 n4 k1 k2 (AnsVar (SField Rational)) -> Int -tensorRank4' t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT4 (singletonTList4 t) +tensorRank4' :: (TIndex k1, TIndex k2, Real a) => AbsTensor4 n1 n2 n3 n4 k1 k2 (AnsVar (SField a)) -> Int +tensorRank4' t = rank $ toMatrixT4 (singletonTList4 t) -tensorRank4 :: (TIndex k1, TIndex k2) => TensList4 k1 k2 (AnsVar (SField Rational)) -> Int -tensorRank4 t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT4 t +tensorRank4 :: (TIndex k1, TIndex k2, Real a) => TensList4 k1 k2 (AnsVar (SField a)) -> Int +tensorRank4 t = rank $ toMatrixT4 t -tensorRank5' :: (TIndex k1, TIndex k2, TIndex k3) => AbsTensor5 n1 n2 n3 n4 n5 k1 k2 k3 (AnsVar (SField Rational)) -> Int -tensorRank5' t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT5 (singletonTList5 t) +tensorRank5' :: (TIndex k1, TIndex k2, TIndex k3, Real a) => AbsTensor5 n1 n2 n3 n4 n5 k1 k2 k3 (AnsVar (SField a)) -> Int +tensorRank5' t = rank $ toMatrixT5 (singletonTList5 t) -tensorRank5 :: (TIndex k1, TIndex k2, TIndex k3) => TensList5 k1 k2 k3 (AnsVar (SField Rational)) -> Int -tensorRank5 t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT5 t +tensorRank5 :: (TIndex k1, TIndex k2, TIndex k3, Real a) => TensList5 k1 k2 k3 (AnsVar (SField a)) -> Int +tensorRank5 t = rank $ toMatrixT5 t -tensorRank6' :: (TIndex k1, TIndex k2, TIndex k3) => AbsTensor6 n1 n2 n3 n4 n5 n6 k1 k2 k3 (AnsVar (SField Rational)) -> Int -tensorRank6' t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT6 (singletonTList6 t) +tensorRank6' :: (TIndex k1, TIndex k2, TIndex k3, Real a) => AbsTensor6 n1 n2 n3 n4 n5 n6 k1 k2 k3 (AnsVar (SField a)) -> Int +tensorRank6' t = rank $ toMatrixT6 (singletonTList6 t) -tensorRank6 :: (TIndex k1, TIndex k2, TIndex k3) => TensList6 k1 k2 k3 (AnsVar (SField Rational)) -> Int -tensorRank6 t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT6 t +tensorRank6 :: (TIndex k1, TIndex k2, TIndex k3, Real a) => TensList6 k1 k2 k3 (AnsVar (SField a)) -> Int +tensorRank6 t = rank $ toMatrixT6 t -tensorRank7' :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4) => AbsTensor7 n1 n2 n3 n4 n5 n6 n7 k1 k2 k3 k4 (AnsVar (SField Rational)) -> Int -tensorRank7' t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT7 (singletonTList7 t) +tensorRank7' :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, Real a) => AbsTensor7 n1 n2 n3 n4 n5 n6 n7 k1 k2 k3 k4 (AnsVar (SField a)) -> Int +tensorRank7' t = rank $ toMatrixT7 (singletonTList7 t) -tensorRank7 :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4) => TensList7 k1 k2 k3 k4 (AnsVar (SField Rational)) -> Int -tensorRank7 t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT7 t +tensorRank7 :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, Real a) => TensList7 k1 k2 k3 k4 (AnsVar (SField a)) -> Int +tensorRank7 t = rank $ toMatrixT7 t -tensorRank8' :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4) => AbsTensor8 n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 (AnsVar (SField Rational)) -> Int -tensorRank8' t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT8 (singletonTList8 t) - -tensorRank8 :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4) => TensList8 k1 k2 k3 k4 (AnsVar (SField Rational)) -> Int -tensorRank8 t = Sol.rank Sol.FullPivLU $ Sparse.toMatrix $ toEMatrixT8 t - +tensorRank8' :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, Real a) => AbsTensor8 n1 n2 n3 n4 n5 n6 n7 n8 k1 k2 k3 k4 (AnsVar (SField a)) -> Int +tensorRank8' t = rank $ toMatrixT8 (singletonTList8 t) +tensorRank8 :: (TIndex k1, TIndex k2, TIndex k3, TIndex k4, Real a) => TensList8 k1 k2 k3 k4 (AnsVar (SField a)) -> Int +tensorRank8 t = rank $ toMatrixT8 t diff --git a/src/Math/Tensor/Examples/Gravity.hs b/src/Math/Tensor/Examples/Gravity.hs --- a/src/Math/Tensor/Examples/Gravity.hs +++ b/src/Math/Tensor/Examples/Gravity.hs @@ -6,7 +6,7 @@ -- Maintainer : tobi.reinhart@fau.de, nils.alex@fau.de -- -- --- This module provides a variety of @'Tensor'@s that are currently predefined in the sparse.tensor package. +-- This module provides a variety of @'Tensor'@s that are currently predefined in the sparse-tensor package. -- -- Amongst many standard tensor from differential geometry and classical field theories such as Kronecker deltas \(\delta^a_b \) in multiple different dimensions, -- the Levi-Civita symbol \(\epsilon^{abcd} \) and the Minkowski metric \(\eta_{ab}\) and its inverse \(\eta^{ab}\), most included tensors were implemented during @@ -27,12 +27,12 @@ -- * Standard Tensors -- ** Kronecker Delta delta3, delta9, delta20, -delta3A, +delta3A, -- ** Minkowski Metric eta, invEta, etaA, invEtaA, etaAbs, -- ** Levi-Civita Symbol epsilon, epsilonInv, epsilonA, epsilonInvA, --- ** Generators of the Lorentz Group \( \mathrm{SO}(3,1)\) +-- ** Generators of the Lorentz Group -- | The following six tensors are a choice of generators of the Lorentz group \( \mathrm{SO}(3,1)\), i.e. they constitute a basis of the -- corresponding Lie algebra \( \mathrm{so}(3,1)\). -- @@ -80,7 +80,7 @@ import Math.Tensor liftSTtoATens :: STTens n1 n2 v -> ATens 0 0 0 0 n1 n2 v -liftSTtoATens = Scalar . Scalar . Scalar . Scalar +liftSTtoATens = Scalar . Scalar . Scalar . Scalar --start with deltas @@ -90,9 +90,9 @@ delta3 :: STTens 1 1 (SField Rational) delta3 = fromListT2 $ zip [(singletonInd (Ind3 i),singletonInd (Ind3 i)) | i <- [0..3]] (repeat $ SField 1) --- | Spacetime Kronecker delta as @'ATens'@. +-- | Spacetime Kronecker delta as @'ATens'@. delta3A :: ATens 0 0 0 0 1 1 (SField Rational) -delta3A = liftSTtoATens delta3 +delta3A = liftSTtoATens delta3 -- | Standard Kronecker delta for the @'Ind9'@ index type \(\delta^I_J\) as @'ATens' 0 0 1 1 0 0 ('SField' 'Rational')@. -- @@ -118,7 +118,7 @@ -- | Minkowski metric lifted to @'ATens'@. etaA :: ATens 0 0 0 0 0 2 (SField Rational) -etaA = liftSTtoATens eta +etaA = liftSTtoATens eta -- | Inverse spacetime Minkowski metric \(\eta^{ab}\) as @'ATens' 0 0 0 0 2 0 ('SField' 'Rational')@. The inverse Minkowski metric could -- also be defined as @'STTens' 2 0 ('SField' 'Rational')@ in similar fashion. @@ -141,13 +141,13 @@ l = map (\(x,y) -> ((Empty, Empty, Empty, singletonInd $ Ind9 x, Empty, Empty),SField y)) [(0,-1),(4,1),(7,1),(9,1)] --- | Covariant spacetime Levi-Civita symbol \(\epsilon_{abcd}\) as type @'ATTens' 0 4 ('SField' 'Rational')@. +-- | Covariant spacetime Levi-Civita symbol \(\epsilon_{abcd}\) as type @'ATTens' 0 4 ('SField' 'Rational')@. epsilon :: STTens 0 4 (SField Rational) -epsilon = fromListT2 l +epsilon = fromListT2 ls where - l = map (\([i,j,k,l],v) -> ((Empty, Append (Ind3 i) $ Append (Ind3 j) $ Append (Ind3 k) $ singletonInd (Ind3 l)),SField v)) epsL - epsSign [i,j,k,l] = (-1) ^ length (filter (==True) [j>i,k>i,l>i,k>j,l>j,l>k]) - epsL = map (\x -> (x, epsSign x)) $ permutations [0,1,2,3] + ls = map (\([i,j,k,l],v) -> ((Empty, Append (Ind3 i) $ Append (Ind3 j) $ Append (Ind3 k) $ singletonInd (Ind3 l)),SField v)) epsL + epsSign i j k l = (-1) ^ length (filter (==True) [j>i,k>i,l>i,k>j,l>j,l>k]) + epsL = map (\x@[i,j,k,l] -> (x, epsSign i j k l)) $ permutations [0,1,2,3] -- | Covariant Levi-Civita symbol lifted to @'ATens'@. epsilonA :: ATens 0 0 0 0 0 4 (SField Rational) @@ -155,11 +155,11 @@ -- | Contravariant spacetime Levi-Civita symbol \(\epsilon^{abcd}\) as type @'STTens'4 0 ('SField' 'Rational')@. T epsilonInv :: STTens 4 0 (SField Rational) -epsilonInv = fromListT2 l +epsilonInv = fromListT2 ls where - l = map (\([i,j,k,l],v) -> ((Append (Ind3 i) $ Append (Ind3 j) $ Append (Ind3 k) $ singletonInd (Ind3 l), Empty),SField v)) epsL - epsSign [i,j,k,l] = (-1) ^ length (filter (==True) [j>i,k>i,l>i,k>j,l>j,l>k]) - epsL = map (\x -> (x, epsSign x)) $ permutations [0,1,2,3] + ls = map (\([i,j,k,l],v) -> ((Append (Ind3 i) $ Append (Ind3 j) $ Append (Ind3 k) $ singletonInd (Ind3 l), Empty),SField v)) epsL + epsSign i j k l = (-1) ^ length (filter (==True) [j>i,k>i,l>i,k>j,l>j,l>k]) + epsL = map (\x@[i,j,k,l] -> (x, epsSign i j k l)) $ permutations [0,1,2,3] -- | Contravariant Levi-Civita symbol lifted to @'ATens'@. epsilonInvA :: ATens 0 0 0 0 4 0 (SField Rational) @@ -218,8 +218,8 @@ | a == b = 1 | otherwise = 1/2 -jMult3 :: (Eq a) => IndList 3 a -> Rational -jMult3 (Append a (Append b (Append c Empty))) +_jMult3 :: (Eq a) => IndList 3 a -> Rational +_jMult3 (Append a (Append b (Append c Empty))) | i == 1 = 1 | i == 2 = 1/3 | otherwise = 1/6 @@ -231,8 +231,8 @@ | a == c && b == d = 1/4 | otherwise = 1/8 -isZeroArea :: (Eq a) => IndList 4 a -> Bool -isZeroArea (Append a (Append b (Append c (Append d Empty)))) = a == b || c == d +_isZeroArea :: (Eq a) => IndList 4 a -> Bool +_isZeroArea (Append a (Append b (Append c (Append d Empty)))) = a == b || c == d areaSign :: (Eq a, Ord a) => IndList 4 a -> Rational areaSign (Append a (Append b (Append c (Append d Empty)))) = s1 * s2 @@ -249,7 +249,7 @@ -- | The tensor \(I^I_{ab} \) maps between covariant @'Ind9'@ indices and symmetric pairs of covariant @'Ind3'@ indices. interI2 :: ATens 0 0 1 0 0 2 (SField Rational) -interI2 = fromListT6 $ fmap (fmap SField) $ filter (\(i,k) -> k /= 0) $ map (\x -> (x,f x)) inds +interI2 = fromListT6 $ fmap (fmap SField) $ filter (\(_,k) -> k /= 0) $ map (\x -> (x,f x)) inds where trian2 = trianMap2 inds = [ (Empty, Empty, singletonInd $ Ind9 a, Empty, Empty, Append (Ind3 b) $ singletonInd $ Ind3 c) | a <- [0..9], b <- [0..3], c <- [0..3]] @@ -259,7 +259,7 @@ -- | The tensor \(J_I^{ab} \) maps between covariant @'Ind9'@ indices and pairs of covariant @'Ind3'@ indices. interJ2 :: ATens 0 0 0 1 2 0 (SField Rational) -interJ2 = fromListT6 $ fmap (fmap SField) $ filter (\(i,k) -> k /= 0) $ map (\x -> (x,f x)) inds +interJ2 = fromListT6 $ fmap (fmap SField) $ filter (\(_,k) -> k /= 0) $ map (\x -> (x,f x)) inds where trian2 = trianMap2 inds = [ (Empty, Empty, Empty, singletonInd $ Ind9 a, Append (Ind3 b) $ singletonInd $ Ind3 c, Empty) | a <- [0..9], b <- [0..3], c <- [0..3]] @@ -269,7 +269,7 @@ -- | The tensor \( I^A_{abcd}\) maps between covariant @'Ind20'@ indices and blocks of @4@ of covariant @'Ind3'@ indices. interIArea :: ATens 1 0 0 0 0 4 (SField Rational) -interIArea = fromListT6 $ fmap (fmap SField) $ filter (\(i,k) -> k /= 0) $ map (\x -> (x,f x)) inds +interIArea = fromListT6 $ fmap (fmap SField) $ filter (\(_,k) -> k /= 0) $ map (\x -> (x,f x)) inds where trianArea = trianMapArea inds = [ (singletonInd (Ind20 a), Empty, Empty, Empty, Empty, Append (Ind3 b) $ Append (Ind3 c) $ Append (Ind3 d) $ singletonInd $ Ind3 e) | a <- [0..20], b <- [0..3], c <- [0..3], d <- [0..3], e <- [0..3], not (b == c || d == e)] @@ -281,7 +281,7 @@ -- | The tensor \( J_A^{abcd}\) maps between contravariant @'Ind20'@ indices and blocks of @4@ of contravariant @'Ind3'@ indices. interJArea :: ATens 0 1 0 0 4 0 (SField Rational) -interJArea = fromListT6 $ fmap (fmap SField) $ filter (\(i,k) -> k /= 0) $ map (\x -> (x,f x)) inds +interJArea = fromListT6 $ fmap (fmap SField) $ filter (\(_,k) -> k /= 0) $ map (\x -> (x,f x)) inds where trianArea = trianMapArea inds = [ (Empty, singletonInd $ Ind20 a, Empty, Empty, Append (Ind3 b) $ Append (Ind3 c) $ Append (Ind3 d) $ singletonInd $ Ind3 e, Empty) | a <- [0..20], b <- [0..3], c <- [0..3], d <- [0..3], e <- [0..3], not (b == c || d == e)] @@ -516,13 +516,13 @@ return randList randArea :: IO (ATens 0 1 0 0 0 0 (SField Rational)) -randArea = do gen <- newTFGen +randArea = do randList <- randRats let inds = map (\i -> (Empty, singletonInd $ Ind20 i, Empty, Empty, Empty, Empty)) [0..20] return $ fromListT6 $ zip inds $ fmap SField randList randAxon :: IO (ATens 0 1 0 0 0 0 (SField Rational)) -randAxon = do gen <- newTFGen +randAxon = do randList <- randRats let inds = map (\i -> (Empty, singletonInd $ Ind20 i, Empty, Empty, Empty, Empty)) [5,9,12] let randInd = SField $ head randList @@ -531,7 +531,7 @@ return tens randFlatArea :: IO (ATens 0 1 0 0 0 0 (SField Rational)) -randFlatArea = do gen <- newTFGen +randFlatArea = do randList <- randRats let assocs = map (\(i,v) -> ( (Empty, singletonInd $ Ind20 i, Empty, Empty, Empty, Empty), SField v)) [(0, -1 * head randList),(5, randList !! 3),(6, -1 * randList !! 1),(9, -1 * randList !! 4),(11, -1 * randList !! 2),(12, randList !! 5),(15, head randList),(18, randList !! 1),(20, randList !! 2)] @@ -540,7 +540,7 @@ randAreaDerivative1 :: IO (ATens 0 1 0 0 0 1 (SField Rational)) -randAreaDerivative1 = do gen <- newTFGen +randAreaDerivative1 = do randList <- randRats let inds = [ f a p | a <- [0..20], p <- [0..3]] return $ fromListT6 $ zip inds $ fmap SField randList @@ -548,7 +548,7 @@ randAreaDerivative2 :: IO (ATens 0 1 0 1 0 0 (SField Rational)) -randAreaDerivative2 = do gen <- newTFGen +randAreaDerivative2 = do randList <- randRats let inds = [ f a i | a <- [0..20], i <- [0..9]] return $ fromListT6 $ zip inds $ fmap SField randList @@ -556,7 +556,7 @@ randMetric :: IO (ATens 0 0 0 1 0 0 (SField Rational)) -randMetric = do gen <- newTFGen +randMetric = do randList <- randRats let inds = map (\i -> (Empty, Empty, Empty, singletonInd $ Ind9 i, Empty, Empty)) [0..20] return $ fromListT6 $ zip inds $ fmap SField randList diff --git a/src/Math/Tensor/Examples/Gravity/DiffeoSymEqns.hs b/src/Math/Tensor/Examples/Gravity/DiffeoSymEqns.hs --- a/src/Math/Tensor/Examples/Gravity/DiffeoSymEqns.hs +++ b/src/Math/Tensor/Examples/Gravity/DiffeoSymEqns.hs @@ -10,7 +10,7 @@ -- These equations arise in the research topic of perturbative constructive gravity and are used there to encode perturbative diffeomorphism invariance. -- -- Providing further usage examples of the sparse-tensor package the equations included here nicely illustrate the syntax that is used when entering tensors. --- They also show how the sparse tensor package can be used to manipulate not only tenors but linear tensorial equations. The sparse-tensor package can for instance be used to extract the information that is contained in such a tensorial equation in matrix form. This then obviously allows one +-- They also show how the sparse-tensor package can be used to manipulate not only tenors but linear tensorial equations. The sparse-tensor package can for instance be used to extract the information that is contained in such a tensorial equation in matrix form. This then obviously allows one -- to computate the rank of the linear tensorial equation or even explicitly solve it. -- -- diff --git a/src/Math/Tensor/Examples/Gravity/Schwarzschild.hs b/src/Math/Tensor/Examples/Gravity/Schwarzschild.hs --- a/src/Math/Tensor/Examples/Gravity/Schwarzschild.hs +++ b/src/Math/Tensor/Examples/Gravity/Schwarzschild.hs @@ -34,8 +34,8 @@ [ ((Empty, Ind3 0 `Append` singletonInd (Ind3 0)), CFun $ \(_:r:_) -> r' r ), ((Empty, Ind3 1 `Append` singletonInd (Ind3 1)), CFun $ \(_:r:_) -> -1/r' r), - ((Empty, Ind3 2 `Append` singletonInd (Ind3 2)), CFun $ \(_:r:_) -> -(r^2)), - ((Empty, Ind3 3 `Append` singletonInd (Ind3 3)), CFun $ \(_:r:theta:_) -> -(r*sin theta)^2) + ((Empty, Ind3 2 `Append` singletonInd (Ind3 2)), CFun $ \(_:r:_) -> -(r**2)), + ((Empty, Ind3 3 `Append` singletonInd (Ind3 3)), CFun $ \(_:r:theta:_) -> -(r*sin theta)**2) ] where r' r = 1 - rs / r @@ -48,8 +48,8 @@ [ ((Ind3 0 `Append` singletonInd (Ind3 0), Empty), CFun $ \(_:r:_) -> 1/r' r), ((Ind3 1 `Append` singletonInd (Ind3 1), Empty), CFun $ \(_:r:_) -> - r' r), - ((Ind3 2 `Append` singletonInd (Ind3 2), Empty), CFun $ \(_:r:_) -> -1/(r^2)), - ((Ind3 3 `Append` singletonInd (Ind3 3), Empty), CFun $ \(_:r:theta:_) -> -1/(r*sin theta)^2) + ((Ind3 2 `Append` singletonInd (Ind3 2), Empty), CFun $ \(_:r:_) -> -1/(r**2)), + ((Ind3 3 `Append` singletonInd (Ind3 3), Empty), CFun $ \(_:r:theta:_) -> -1/(r*sin theta)**2) ] where r' r = 1 - rs / r diff --git a/src/Math/Tensor/Internal/LinearAlgebra.hs b/src/Math/Tensor/Internal/LinearAlgebra.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/Internal/LinearAlgebra.hs @@ -0,0 +1,193 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Math.Tensor.Internal.LinearAlgebra +-- Copyright : (c) 2019 Tobias Reinhart and Nils Alex +-- License : MIT +-- Maintainer : tobi.reinhart@fau.de, nils.alex@fau.de +-- +-- Gaussian elimination algorithm based on hmatrix. +----------------------------------------------------------------------------- +module Math.Tensor.Internal.LinearAlgebra ( +-- * Gaussian Elimination +gaussianST, +gaussian, +-- * Linearly Independent Columns +independentColumns, +independentColumnsMat, +-- * Pivots +pivotsU, +findPivotMax) + +where + +import Numeric.LinearAlgebra +import Numeric.LinearAlgebra.Data +import Numeric.LinearAlgebra.Devel + +import Data.List (maximumBy) + +import Control.Monad +import Control.Monad.ST + +-- | Returns the pivot columns of an upper triangular matrix. +-- +-- @ +-- λ let mat = (3 >< 4) [1, 0, 2, -3, 0, 0, 1, 0, 0, 0, 0, 0] +-- λ mat +-- (3><4) +-- [ 1.0, 0.0, 2.0, -3.0 +-- , 0.0, 0.0, 1.0, 0.0 +-- , 0.0, 0.0, 0.0, 0.0 ] +-- λ pivotsU mat +-- [0,2] +-- @ +-- + +pivotsU :: Matrix Double -> [Int] +pivotsU mat = go (0,0) + where + go (i,j) + = case findPivot mat (i,j) of + Nothing -> [] + Just (i', j') -> j' : go (i'+1, j'+1) + + +-- naive check for numerical zero + +eps :: Double -> Bool +eps = (< 1e-12) . abs + +-- find next pivot in upper triangular matrix + +findPivot :: Matrix Double -> (Int, Int) -> Maybe (Int, Int) +findPivot mat (i, j) + | n == j = Nothing + | m == i = Nothing + | otherwise = case nonZeros of + [] -> if n == j+1 + then Nothing + else findPivot mat (i, j+1) + (pi, pj):_ -> Just (pi, pj+j) + where + m = rows mat + n = cols mat + col = mat ¿ [j] + nonZeros = filter (\(i', _) -> i' >= i) $ find (not . eps) col + +-- | Find pivot element below position (i, j) with greatest absolute value. + +findPivotMax :: Matrix Double -> Int -> Int -> Maybe (Int, Int) +findPivotMax mat i j + | n == j = Nothing + | m == i = Nothing + | otherwise = case nonZeros of + [] -> if n == j+1 + then Nothing + else findPivotMax mat i (j+1) + _ -> Just (pi, pj+j) + where + m = rows mat + n = cols mat + col = mat ¿ [j] + nonZeros = filter (\(i', _) -> i' >= i) $ find (not . eps) col + (pi, pj) = maximumBy (\ix jx -> abs (col `atIndex` ix) + `compare` + abs (col `atIndex` jx)) + nonZeros + +-- gaussian elimination of sub matrix below position (i, j) + +gaussian' :: Int -> Int -> STMatrix s Double -> ST s () +gaussian' i j mat = do + m <- liftSTMatrix rows mat + n <- liftSTMatrix cols mat + iPivot' <- liftSTMatrix (\x -> findPivotMax x i j) mat + case iPivot' of + Nothing -> return () + Just (r, p) -> do + rowOper (SWAP i r (FromCol j)) mat + pv <- liftSTMatrix (`atIndex` (i, p)) mat + mapM_ (reduce pv p) [i+1 .. m-1] + gaussian' (i+1) (p+1) mat + where + reduce pv p r = do + rv <- liftSTMatrix (`atIndex` (r, p)) mat + if eps rv + then return () + else + let frac = -rv / pv + op = AXPY frac i r (FromCol p) + in rowOper op mat + +-- | Gaussian elimination perfomed in-place in the @'ST'@ monad. + +gaussianST :: STMatrix s Double -> ST s () +gaussianST = gaussian' 0 0 + + +-- | Gaussian elimination as pure function. Involves a copy of the input matrix. +-- +-- @ +-- λ let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1] +-- λ mat +-- (3><4) +-- [ 1.0, 1.0, -2.0, 0.0 +-- , 0.0, 2.0, -6.0, -4.0 +-- , 3.0, 0.0, 3.0, 1.0 ] +-- λ gaussian mat +-- (3><4) +-- [ 3.0, 0.0, 3.0, 1.0 +-- , 0.0, 2.0, -6.0, -4.0 +-- , 0.0, 0.0, 0.0, 1.6666666666666667 ] +-- @ +-- + +gaussian :: Matrix Double -> Matrix Double +gaussian mat = runST $ do + m <- thawMatrix mat + gaussianST m + freezeMatrix m + +-- | Returns the indices of a maximal linearly independent subset of the columns +-- in the matrix. +-- +-- @ +-- λ let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1] +-- λ mat +-- (3><4) +-- [ 1.0, 1.0, -2.0, 0.0 +-- , 0.0, 2.0, -6.0, -4.0 +-- , 3.0, 0.0, 3.0, 1.0 ] +-- λ independentColumns mat +-- [0,1,3] +-- @ +-- + +independentColumns :: Matrix Double -> [Int] +independentColumns mat = pivotsU mat' + where + mat' = gaussian mat + +-- | Returns a sub matrix containing a maximal linearly independent subset of +-- the columns in the matrix. +-- +-- @ +-- λ let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1] +-- λ mat +-- (3><4) +-- [ 1.0, 1.0, -2.0, 0.0 +-- , 0.0, 2.0, -6.0, -4.0 +-- , 3.0, 0.0, 3.0, 1.0 ] +-- λ independentColumnsMat mat +-- (3><3) +-- [ 1.0, 1.0, 0.0 +-- , 0.0, 2.0, -4.0 +-- , 3.0, 0.0, 1.0 ] +-- @ +-- + +independentColumnsMat :: Matrix Double -> Matrix Double +independentColumnsMat mat = + case independentColumns mat of + [] -> (rows mat >< 1) $ repeat 0 + cs -> mat ¿ cs diff --git a/src/Math/Tensor/LorentzGenerator.hs b/src/Math/Tensor/LorentzGenerator.hs --- a/src/Math/Tensor/LorentzGenerator.hs +++ b/src/Math/Tensor/LorentzGenerator.hs @@ -47,7 +47,7 @@ -- ** Node Types Eta(..), Epsilon(..), Var(..), -- ** Forest types -AnsatzForestEpsilon(..), AnsatzForestEta(..), +AnsatzForestEpsilon, AnsatzForestEta(..), -- ** Conversions of AnsatzForests -- *** List of Branches flattenForest, flattenForestEpsilon, forestEtaList, forestEpsList, forestEtaListLatex, forestEpsListLatex, @@ -77,11 +77,11 @@ -- ** The Memory Optimized Way -- The following functions essentially compute the same results as their __Fast__ counterparts, with the only distinction being that they employ a slightly different -- algorithm that avoids the problem of intermediate memory swelling and thus yields improved memory usage. All this is achieved at the cost of slightly higher computation times compared to the __Fast__ functions. -mkAnsatzTensorEigSym, mkAnsatzTensorEig, mkAnsatzTensorEigAbs, -mkAnsatzTensorEigSym', mkAnsatzTensorEig', +mkAnsatzTensorIncrementalSym, mkAnsatzTensorIncremental, mkAnsatzTensorIncrementalAbs, +mkAnsatzTensorIncrementalSym', mkAnsatzTensorIncremental', -- * Specifying Additional Data -- ** Symmetry Type -Symmetry(..), +Symmetry, -- ** Evaluation Lists -- *** Area Metric -- | The following provides an example of evaluation lists. @@ -95,29 +95,30 @@ symList4, symList6, symList8, symList10_1, symList10_2, symList12, symList14_1, symList14_2, -- *** Metric -- | The following are examples of symmetry lists. -metricsymList2, metricsymList4_1, metricsymList4_2, metricsymList6_1, metricsymList6_2, metricsymList6_3, metricsymList8_1, metricsymList8_2 +metricsymList2, metricsymList4_1, metricsymList4_2, metricsymList6_1, metricsymList6_2, metricsymList6_3, metricsymList8_1, metricsymList8_2, +symList16_1, areaList16_1 ) where import qualified Data.IntMap.Strict as I import qualified Data.Map.Strict as M -import Data.List (nub, permutations, foldl', (\\), elemIndex, nubBy, sortBy, insert, intersect, union, partition, delete) +import Data.List (nub, permutations, foldl', (\\), elemIndex, nubBy, sortBy, insert, intersect, union, partition, delete, maximumBy, splitAt) import Data.Maybe (fromJust, isNothing, fromMaybe, isJust, mapMaybe) -import Control.Parallel.Strategies (parListChunk, rdeepseq, runEval, NFData(..)) +import Control.Parallel.Strategies (parListChunk, rdeepseq, runEval, NFData) import Data.Serialize (encodeLazy, decodeLazy, Serialize(..)) import GHC.Generics -import qualified Data.ByteString.Lazy as BS (ByteString(..)) +import qualified Data.ByteString.Lazy as BS (ByteString) import Codec.Compression.GZip (compress, decompress) import Data.Either (either) import Data.Tuple (swap) import GHC.TypeLits -import Data.Singletons (SingI(..)) --LinearAlgebra subroutines -import qualified Data.Eigen.Matrix as Mat -import qualified Data.Eigen.SparseMatrix as Sparse -import qualified Data.Eigen.LA as Sol +import qualified Numeric.LinearAlgebra.Data as HM +import qualified Numeric.LinearAlgebra as Matrix +import Math.Tensor.Internal.LinearAlgebra (independentColumns) + import Math.Tensor {-- @@ -141,8 +142,6 @@ sortPair (a,b) = if a < b then (a,b) else (b,a) f1 = pairs ++ aPairs f2 = map (\(a,b) -> (head a, head b)) blocks - getPairs [a,b] = [(a,b)] - getPairs (x:xs) = (x, head xs) : getPairs xs f3 = concatMap getPairs cycles f4 = concatMap (getPairs . map head) blockCycles @@ -169,13 +168,14 @@ --} getExtraSyms1 :: [Int] -> Symmetry -> Symmetry -getExtraSyms1 [] syms = ([],[],[],[],[]) +getExtraSyms1 [] _ = ([],[],[],[],[]) +getExtraSyms1 [_] _ = error "cannot get extra syms from singleton index list" getExtraSyms1 (a:b:xs) (pairs,aPairs,blocks,cycles,blockCycles) = addSym (newPairs, [], newBlocks, [], []) (getExtraSyms1 xs newSyms) where allBlocks = blocks ++ concatMap mkBlocksFromBlockCycle blockCycles newBlocks' = map (\(x,y) -> unzip $ filter (\(c,d) -> (c,d) /= (a,b)) $ zip x y) allBlocks - (newBlocks, newPairs') = partition (\(a,b) -> length a > 1) newBlocks' - newPairs = map (\([a],[b]) -> (a,b)) newPairs' + (newBlocks, newPairs') = partition (\(a',_) -> length a' > 1) newBlocks' + newPairs = map (\([a'],[b']) -> (a',b')) newPairs' newSyms = addSym (pairs,aPairs,blocks,cycles,blockCycles) (newPairs, [], newBlocks, [], []) mkBlocksFromBlockCycle :: [[Int]] -> [([Int],[Int])] @@ -183,6 +183,7 @@ mkBlocksFromBlockCycle (x:xs) = l ++ mkBlocksFromBlockCycle xs where l = map (x,) xs +mkBlocksFromBlockCycle _ = error "cannot make block symmetries from empty block list" {-- Furthermore distributing a symmetric or antisymmetric pair of indices over 2 etas yields an additional symmetry or anti-symmetry @@ -194,7 +195,7 @@ --given one eta, if the eta contains an index from a symmetric or antisymmetric pair return the corresponding second index and the other index of the eta get2nd :: [Int] -> Symmetry -> (Maybe [(Int,Int)], Maybe [(Int,Int)]) -get2nd [a,b] (pairs,aPairs,blocks,cycles,blockCycles) = (sndPairs, sndAPairs) +get2nd [a,b] (pairs,aPairs,_,cycles,_) = (sndPairs, sndAPairs) where allPairs = pairs ++ concatMap mkSymsFromCycle cycles aPair = lookup a allPairs @@ -211,15 +212,16 @@ (Just x, Nothing) -> Just [(b,x)] (Nothing, Just y) -> Just [(a,y)] (Just x, Just y) -> if x == b then Nothing else Just [(b,x),(a,y)] +get2nd _ _ = error "given index list contains more or less than two indices" --find the eta that contains the computed second pair index and return the other indices of this eta get2ndSyms :: Maybe [(Int,Int)] -> Symmetry -> [[Int]] -> Symmetry -get2ndSyms Nothing syms etas = syms -get2ndSyms (Just i) (pairs,aPairs,blocks,cycles,blockCycles) etas = (newPairs,[],[],[],[]) +get2ndSyms Nothing syms _ = syms +get2ndSyms (Just i) _ etas = (newPairs,[],[],[],[]) where - get2ndInd l (i,j) = mapMaybe (\[a,b] -> if j == a then Just (i,b) else if j == b then Just (i,a) else Nothing) l + get2ndInd l (i',j) = mapMaybe (\[a,b] -> if j == a then Just (i',b) else if j == b then Just (i',a) else Nothing) l newPairs = concatMap (get2ndInd etas) i mkSymsFromCycle :: [Int] -> [(Int,Int)] @@ -227,15 +229,21 @@ mkSymsFromCycle (x:xs) = l ++ mkSymsFromCycle xs where l = map (x,) xs +mkSymsFromCycle _ = error "cannot make syms from empty cycle list" get2ndASyms :: Maybe [(Int,Int)] -> Symmetry -> [[Int]] -> Symmetry -get2ndASyms Nothing syms etas = syms -get2ndASyms (Just i) (pairs,aPairs,blocks,cycles,blockCycles) etas = ([], newAPairs,[],[],[]) +get2ndASyms Nothing syms _ = syms +get2ndASyms (Just i) _ etas = ([], newAPairs,[],[],[]) where - get2ndInd l (i,j) = mapMaybe (\[a,b] -> if j == a then Just (i,b) else if j == b then Just (i,a) else Nothing) l + get2ndInd l (i',j) = mapMaybe (\[a,b] -> if j == a then Just (i',b) else if j == b then Just (i',a) else Nothing) l newAPairs = concatMap (get2ndInd etas) i +mkEtas :: [Int] -> [[Int]] +mkEtas [] = [] +mkEtas [l,k] = [[l,k]] +mkEtas (l:k:ls) = [l,k] : mkEtas ls +mkEtas _ = error "cannot make etas from singleton index list" --apply to whole ind list @@ -243,14 +251,12 @@ getExtraSyms2 [] syms = syms getExtraSyms2 (a':b':xs) syms = addSym (getExtraSyms2 xs newSyms) newSyms where - mkEtas [] = [] - mkEtas [l,k] = [[l,k]] - mkEtas (l:k:ls) = [l,k] : mkEtas ls x = [a',b'] (i,j) = get2nd x syms (p,_,_,_,_) = get2ndSyms i syms (mkEtas xs) (_,a,_,_,_) = get2ndASyms j syms (mkEtas xs) newSyms = addSym (p,a,[],[],[]) syms +getExtraSyms2 _ _ = error "cannot get extra syms from singleton index list" --compute all extra symmetries @@ -262,11 +268,12 @@ getAllIndsEta :: [Int] -> [(Int,Int)] -> [[Int]] -getAllIndsEta [a,b] aSyms = [[a,b]] +getAllIndsEta [a,b] _ = [[a,b]] getAllIndsEta (x:xs) aSyms = concatMap res firstEta where firstEta = mapMaybe (\y -> if (x,y) `notElem` aSyms then Just ([x,y],delete y xs) else Nothing) xs res (a,b) = (++) a <$> getAllIndsEta b aSyms +getAllIndsEta _ _ = error "empty index list" filterEta :: [Int] -> Symmetry -> [(Int,Int)] -> Bool filterEta inds (p1,ap1,b1,c1,cb1) filters = filterSym inds totFilters && isNonZero @@ -274,16 +281,13 @@ (p2,ap2,b2,c2,cb2) = getAllExtraSyms inds (p1,ap1,b1,c1,cb1) extrafilters = mkFilters (p2,ap2,b2,c2,cb2) totFilters = filters `union` extrafilters - mkEtas [] = [] - mkEtas [l,k] = [(l,k)] - mkEtas (l:k:ls) = (l,k) : mkEtas ls - etas = mkEtas inds + etas = map (\[a,b] -> (a,b)) $ mkEtas inds isNonZero = null $ etas `intersect` union ap1 ap2 --construct a pre-reduced list of eta indices getEtaInds :: [Int] -> Symmetry -> [[Int]] -getEtaInds [] sym = [[]] +getEtaInds [] _ = [[]] getEtaInds inds (p,ap,b,c,bc) = filter (\x -> filterEta x (p,ap,b,c,bc) filters1) allInds where filters1 = mkFilters (p,ap,b,c,bc) @@ -305,29 +309,32 @@ --get all possible epsilon inds that are allowed under the above considerations getAllIndsEpsilon :: [Int] -> Symmetry -> [[Int]] -getAllIndsEpsilon inds (p,ap,b,cyc,cb) = [ [a,b,c,d] | a <- [1..i-3], b <- [a+1..i-2], c <- [b+1..i-1], d <- [c+1..i], +getAllIndsEpsilon inds (p,ap,bs,cyc,_) = [ [a,b,c,d] | a <- [1..numInds-3], b <- [a+1..numInds-2], c <- [b+1..numInds-1], d <- [c+1..numInds], not (isSym p [a,b,c,d]) && not (is3Area areaBlocks [a,b,c,d]) && isValid2Area areaBlocks [a,b,c,d] && not (is1Area areaBlocks [a,b,c,d]) && not (isSymCyc cyc [a,b,c,d]) ] where - i = length inds - blocks2 = filter (\x -> length (fst x) == 2) b + numInds = length inds + blocks2 = filter (\x -> length (fst x) == 2) bs areaBlocks = map (uncurry (++)) $ filter (\([a,b],[c,d]) -> (a,b) `elem` ap && (c,d) `elem` ap) blocks2 - isSym [] x = False + isSym [] _ = False isSym [(a,b)] [i,j,k,l] = length ([a,b] `intersect` [i,j,k,l]) == 2 isSym (x:xs) [i,j,k,l] | isSym [x] [i,j,k,l] = True | otherwise = isSym xs [i,j,k,l] - isSymCyc [] x = False + isSym _ _ = error "expected four indices" + isSymCyc [] _ = False isSymCyc [l'] [i,j,k,l] = length (l' `intersect` [i,j,k,l]) >= 2 isSymCyc (x:xs) [i,j,k,l] | isSymCyc [x] [i,j,k,l] = True | otherwise = isSymCyc xs [i,j,k,l] - is3Area [] i = False + isSymCyc _ _ = error "expected four indices" + is3Area [] _ = False is3Area [[a,b,c,d]] [i,j,k,l] = length ([a,b,c,d] `intersect` [i,j,k,l]) == 3 is3Area (x:xs) [i,j,k,l] | is3Area [x] [i,j,k,l] = True | otherwise = is3Area xs [i,j,k,l] - isValid2Area [] i = True + is3Area _ _ = error "expected four indices" + isValid2Area [] _ = True isValid2Area [[a,b,c,d]] [i,j,k,l] | length inter == 2 = inter == [a,b] | otherwise = True @@ -336,16 +343,18 @@ isValid2Area (x:xs) [i,j,k,l] | isValid2Area [x] [i,j,k,l] = isValid2Area xs [i,j,k,l] | otherwise = False - is1Area [] i = False + isValid2Area _ _ = error "expected four indices" + is1Area [] _ = False is1Area list [i,j,k,l] = maximum (map (length . ([i,j,k,l] `intersect`)) list) == 1 + is1Area _ _ = error "expected four indices" --a 2-block symmetry with the respectively first indices at an epsilon yields an additional anti-symmetry (note that we did not include higher block anti-symmetries) getExtraASymsEps :: [Int] -> Symmetry -> Symmetry -getExtraASymsEps eps (p,ap,blo,cyc,cb) = ([],newASyms, [], [], []) +getExtraASymsEps eps (_,_,blo,_,cb) = ([],newASyms, [], [], []) where allBlocks = blo ++ concatMap mkBlocksFromBlockCycle cb - blocks2 = filter (\(a,b) -> length a == 2) allBlocks + blocks2 = filter (\(a,_) -> length a == 2) allBlocks newASyms = mapMaybe (\([i,j],[k,l]) -> if length ([i,k] `intersect` eps) == 2 then Just (j,l) else if length ([j,l] `intersect` eps) == 2 then Just (i,k) else Nothing) blocks2 getEpsilonInds :: [Int] -> Symmetry -> [[Int]] @@ -353,11 +362,6 @@ where epsInds = getAllIndsEpsilon inds sym allInds = concat $ filter (not . null) $ map (\x -> map (x ++) $ getEtaInds (inds \\ x) (addSym sym (getExtraASymsEps x sym)) )epsInds - isSymP [] x = False - isSymP [(a,b)] [i,j,k,l] = length ([a,b] `intersect` [i,j,k,l]) == 2 - isSymP (x:xs) [i,j,k,l] - | isSymP [x] [i,j,k,l] = True - | otherwise = isSymP xs [i,j,k,l] filters = mkFilters sym allIndsRed = filter (\x -> let symEps = addSym (getExtraASymsEps (take 4 x) sym) sym symEta = addSym symEps (getAllExtraSyms (drop 4 x) symEps) @@ -432,16 +436,12 @@ decodeAnsatzForestEpsilon :: BS.ByteString -> AnsatzForestEpsilon decodeAnsatzForestEpsilon bs = either error id $ decodeLazy $ decompress bs -forestMap :: AnsatzForestEta -> M.Map Eta AnsatzForestEta -forestMap (ForestEta m) = m -{-# INLINEABLE forestMap #-} - --map a function over the nodes of the AnsatzTree (map over the tensors eta and epsilon) mapNodes :: (Eta -> Eta) -> AnsatzForestEta -> AnsatzForestEta -mapNodes f EmptyForest = EmptyForest +mapNodes _ EmptyForest = EmptyForest mapNodes f (ForestEta m) = ForestEta $ M.mapKeys f . M.map (mapNodes f) $ m -mapNodes f (Leaf x) = Leaf x +mapNodes _ (Leaf x) = Leaf x mapNodesEpsilon :: (Epsilon -> Epsilon) -> AnsatzForestEpsilon -> AnsatzForestEpsilon mapNodesEpsilon = M.mapKeys @@ -450,7 +450,7 @@ -- | Map a general function over all variables that are contained in the @'AnsatzForestEta'@. mapVars :: (Var -> Var) -> AnsatzForestEta -> AnsatzForestEta -mapVars f EmptyForest = EmptyForest +mapVars _ EmptyForest = EmptyForest mapVars f (Leaf var) = Leaf (f var) mapVars f (ForestEta m) = ForestEta $ M.map (mapVars f) m @@ -458,14 +458,6 @@ mapVarsEpsilon :: (Var -> Var) -> AnsatzForestEpsilon -> AnsatzForestEpsilon mapVarsEpsilon f = M.map (mapVars f) ---multiplying the vars with a fixed Int - -multVars :: Int -> AnsatzForestEta -> AnsatzForestEta -multVars i = mapVars (multVar i) - -multVarsEpsilon :: Int -> AnsatzForestEpsilon -> AnsatzForestEpsilon -multVarsEpsilon i = mapVarsEpsilon (multVar i) - --relabel and remove Vars in the Forest getLeafVals :: AnsatzForestEta -> [Var] @@ -473,12 +465,13 @@ getLeafVals (ForestEta m) = rest where rest = concatMap getLeafVals $ M.elems m +getLeafVals EmptyForest = [] getLeafValsEpsilon :: AnsatzForestEpsilon -> [Var] getLeafValsEpsilon m = concatMap getLeafVals $ M.elems m getVarLabels :: Var -> Int -getVarLabels (Var i j) = j +getVarLabels (Var _ j) = j -- | Return a list of the labels of all variables that are contained in the @'AnsatzForestEta'@. getForestLabels :: AnsatzForestEta -> [Int] @@ -514,7 +507,7 @@ | j `elem` vars = EmptyForest | otherwise = Leaf (Var i j) removeVarsEta vars (ForestEta m) = ForestEta $ M.filter (/= EmptyForest) $ M.map (removeVarsEta vars) m -removeVarsEta vars EmptyForest = EmptyForest +removeVarsEta _ EmptyForest = EmptyForest -- | Shift the variable labels of all variables that are contained in the @'AnsatzForestEpsilon'@ by the amount specified. relabelAnsatzForestEpsilon :: Int -> AnsatzForestEpsilon -> AnsatzForestEpsilon @@ -543,6 +536,7 @@ | otherwise = ForestEta newMap where newMap = M.filter (/= EmptyForest) $ M.unionWith addForests m1 m2 +addForests _ _ = error "cannot add Leaf and Forest" addForestsEpsilon :: AnsatzForestEpsilon -> AnsatzForestEpsilon -> AnsatzForestEpsilon addForestsEpsilon m1 m2 = M.filter (/= EmptyForest) $ M.unionWith addForests m1 m2 @@ -554,12 +548,14 @@ | otherwise = Leaf newLeafVal where newLeafVal = addVars var1 var2 -addList2Forest (ForestEta m1) (x:xs, var) = ForestEta $ M.insertWith (\a1 a2 -> addList2Forest a2 (xs, var)) x newVal m1 +addList2Forest (ForestEta m1) (x:xs, var) = ForestEta $ M.insertWith (\_ a2 -> addList2Forest a2 (xs, var)) x newVal m1 where newVal = mkForestFromAscList (xs,var) +addList2Forest (ForestEta f) ([], _) = ForestEta f +addList2Forest (Leaf _) _ = error "cannot add something to Leaf" addList2ForestEpsilon :: AnsatzForestEpsilon -> (Epsilon,[Eta],Var) -> AnsatzForestEpsilon -addList2ForestEpsilon m (eps,eta,var) = M.insertWith (\a1 a2 -> addList2Forest a2 (eta, var)) eps newVal m +addList2ForestEpsilon m (eps,eta,var) = M.insertWith (\_ a2 -> addList2Forest a2 (eta, var)) eps newVal m where newVal = mkForestFromAscList (eta,var) @@ -592,7 +588,7 @@ | M.size x == 1 = let [(a,b)] = M.assocs x in "|" : shift "`---- " " " (drawEtaTree a b) | otherwise = let (a,b) = head $ M.assocs x in "|" : shift "+---- " "| " (drawEtaTree a b) ++ drawSubTrees (M.delete a x) shift first other = zipWith (++) (first : repeat other) -drawEtaTree eta EmptyForest = [] +drawEtaTree _ EmptyForest = [] drawEpsilonTree :: Epsilon -> AnsatzForestEta -> [String] drawEpsilonTree (Epsilon i j k l) (Leaf (Var a b)) = ["(" ++ show i ++ "," ++ show j ++ "," ++ show k ++ "," ++ show l ++ ") * (" ++ show a ++ ") * x[" ++ show b ++ "]"] @@ -603,7 +599,7 @@ | M.size x == 1 = let [(a,b)] = M.assocs x in "|" : shift "`---- " " " (drawEtaTree a b) | otherwise = let (a,b) = head $ M.assocs x in "|" : shift "+---- " "| " (drawEtaTree a b) ++ drawSubTrees (M.delete a x) shift first other = zipWith (++) (first : repeat other) -drawEpsilonTree eps EmptyForest = [] +drawEpsilonTree _ EmptyForest = [] -- | Returns an ASCII drawing of the @'AnsatzForestEta'@ in the fashion explained in "Data.Tree". -- The ansatz \( x_1 \cdot 8 \{ \eta^{ac}\eta^{bd}\eta^{pq} - \eta^{ad}\eta^{bc}\eta^{pq} \} + x_2 \cdot 2 \{\eta^{ac}\eta^{bp}\eta^{dq} + \eta^{ac}\eta^{bq}\eta^{dp} - \eta^{bc}\eta^{ap}\eta^{dq} - \eta^{bc}\eta^{aq}\eta^{dp} - \eta^{ad}\eta^{bp}\eta^{cq} - \eta^{ad}\eta^{bq}\eta^{cp} + \eta^{bd}\eta^{ap}\eta^{cq} + \eta^{bd}\eta^{aq}\eta^{cp} \} \) is drawn to @@ -679,17 +675,17 @@ forestEtaList f = map fst fList'' where fList = flattenForest f - fList' = sortBy (\(e1, Var x1 y1 ) (e2, Var x2 y2) -> compare y1 y2) fList - fList'' = nubBy (\(e1, Var x1 y1 ) (e2, Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) fList' + fList' = sortBy (\(_, Var _ y1 ) (_, Var _ y2) -> compare y1 y2) fList + fList'' = nubBy (\(_, Var x1 y1 ) (_, Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) fList' -- | Return one representative, i.e. one individual product for each of the basis ansätze in an @'AnsatzForestEpsilon'@. The function thus returns the contained individual ansätze without -- their explicit symmetrization. forestEpsList :: AnsatzForestEpsilon -> [(Epsilon,[Eta])] -forestEpsList f = map (\(a,b,c) -> (a,b)) fList'' +forestEpsList f = map (\(a,b,_) -> (a,b)) fList'' where fList = flattenForestEpsilon f - fList' = sortBy (\(e1, e', Var x1 y1 ) (e2, e2', Var x2 y2) -> compare y1 y2) fList - fList'' = nubBy (\(e1, e1', Var x1 y1 ) (e2, e2', Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) fList' + fList' = sortBy (\(_, _, Var _ y1 ) (_, _, Var _ y2) -> compare y1 y2) fList + fList'' = nubBy (\(_, _, Var x1 y1 ) (_, _, Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) fList' --output in latex format @@ -702,9 +698,9 @@ forestEtaListLatex :: AnsatzForestEta -> String -> Char -> String forestEtaListLatex f inds var = tail $ concat etaL'' where - etaL = sortBy (\(e1, Var x1 y1 ) (e2, Var x2 y2) -> compare y1 y2) $ flattenForest f - etaL' = nubBy (\(e1, Var x1 y1 ) (e2, Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) etaL - etaL'' = map (\(a,Var x y) -> "+" ++ var : "_{" ++ show y ++ "}\\cdot" ++ concatMap (mkEtasLatex inds) a) etaL' + etaL = sortBy (\(_, Var _ y1 ) (_, Var _ y2) -> compare y1 y2) $ flattenForest f + etaL' = nubBy (\(_, Var x1 y1 ) (_, Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) etaL + etaL'' = map (\(a,Var _ y) -> "+" ++ var : "_{" ++ show y ++ "}\\cdot" ++ concatMap (mkEtasLatex inds) a) etaL' mkEpsLatex :: String -> Epsilon -> String mkEpsLatex inds (Epsilon i j k l) = "\\epsilon^{" ++ epsi : epsj : epsk : epsl : "}" @@ -715,9 +711,9 @@ forestEpsListLatex :: AnsatzForestEpsilon -> String -> Char -> String forestEpsListLatex f inds var = tail $ concat epsL'' where - epsL = sortBy (\(e1, e1', Var x1 y1 ) (e2, e2', Var x2 y2) -> compare y1 y2) $ flattenForestEpsilon f - epsL' = nubBy (\(e1, e1', Var x1 y1 ) (e2, e2', Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) epsL - epsL'' = map (\(a,b,Var x y) -> "+" ++ var : "_{" ++ show y ++ "}\\cdot" ++ mkEpsLatex inds a ++ concatMap (mkEtasLatex inds) b) epsL' + epsL = sortBy (\(_, _, Var _ y1 ) (_, _, Var _ y2) -> compare y1 y2) $ flattenForestEpsilon f + epsL' = nubBy (\(_, _, Var x1 y1 ) (_, _, Var x2 y2) -> if x1 == 0 || x2 == 0 then error "zeros!!" else y1 == y2) epsL + epsL'' = map (\(a,b,Var _ y) -> "+" ++ var : "_{" ++ show y ++ "}\\cdot" ++ mkEpsLatex inds a ++ concatMap (mkEtasLatex inds) b) epsL' --construct a forest of a given asclist @@ -913,13 +909,15 @@ --look up a 1d Forest (obtained from the index list) in the given Forest isElem :: [Eta] -> AnsatzForestEta -> Bool -isElem [] (Leaf x) = True -isElem x EmptyForest = False -isElem (x:xs) (ForestEta m) = case mForest of +isElem [] (Leaf _) = True +isElem _ EmptyForest = False +isElem (x:xs) (ForestEta m) = case mForest of Just forest -> xs `isElem` forest _ -> False where mForest = M.lookup x m +isElem _ (Leaf _) = error "cannot lookup eta in Leaf" +isElem [] (ForestEta _) = error "cannot look for eta in forest when no eta is specified" isElemEpsilon :: (Epsilon, [Eta]) -> AnsatzForestEpsilon -> Bool isElemEpsilon (eps,l) m = case mForest of @@ -981,20 +979,20 @@ l = [(I.!) iMap w, (I.!) iMap x, (I.!) iMap y, (I.!) iMap z] epsMap :: M.Map [Int] Int -epsMap = M.fromList $ map (\x -> (x, epsSign x)) $ permutations [0,1,2,3] +epsMap = M.fromList $ map (\x@[i,j,k,l] -> (x, epsSign i j k l)) $ permutations [0,1,2,3] where - epsSign [i,j,k,l] = (-1) ^ length (filter (==True) [j>i,k>i,l>i,k>j,l>j,l>k]) + epsSign i j k l = (-1) ^ length (filter (==True) [j>i,k>i,l>i,k>j,l>j,l>k]) --basic tree eval function evalAnsatzForestEta :: I.IntMap Int -> AnsatzForestEta -> I.IntMap Int -evalAnsatzForestEta evalM (Leaf (Var x y)) = I.singleton y x +evalAnsatzForestEta _ (Leaf (Var x y)) = I.singleton y x evalAnsatzForestEta evalM (ForestEta m) = M.foldlWithKey' foldF I.empty m where foldF b k a = let nodeVal = evalNodeEta evalM k in if isNothing nodeVal then b else I.unionWith (+) (I.map (fromJust nodeVal *) (evalAnsatzForestEta evalM a)) b -evalAnsatzForestEta evalM EmptyForest = I.empty +evalAnsatzForestEta _ EmptyForest = I.empty evalAnsatzForestEpsilon :: I.IntMap Int -> AnsatzForestEpsilon -> I.IntMap Int evalAnsatzForestEpsilon evalM = M.foldlWithKey' foldF I.empty @@ -1006,13 +1004,13 @@ --for a single Ansatz we do not need the IntMap to keep track of the VarLabels -> eval to a number eval1AnsatzForestEta :: I.IntMap Int -> AnsatzForestEta -> Int -eval1AnsatzForestEta evalM (Leaf (Var x _)) = x +eval1AnsatzForestEta _ (Leaf (Var x _)) = x eval1AnsatzForestEta evalM (ForestEta m) = M.foldlWithKey' foldF 0 m where foldF b k a = let nodeVal = evalNodeEta evalM k in if isNothing nodeVal then b else b + (fromJust nodeVal * eval1AnsatzForestEta evalM a) -eval1AnsatzForestEta evalM EmptyForest = 0 +eval1AnsatzForestEta _ EmptyForest = 0 eval1AnsatzForestEpsilon :: I.IntMap Int -> AnsatzForestEpsilon -> Int eval1AnsatzForestEpsilon evalM = M.foldlWithKey' foldF 0 @@ -1021,33 +1019,34 @@ in if isNothing nodeVal then b else b + (fromJust nodeVal * eval1AnsatzForestEta evalM a) ---eval a given 1Var ansatz to a sparse Matrix (a row vector) -> Eigen Indices start at 0 !! +--eval a given 1Var ansatz to a row vector -> HMatrix Indices start at 0 !! -mkVecList mkAns dofList evalM = vecList +mkVecList :: (Foldable t, NFData a1, Real a1) => + (a2 -> Maybe ((Int, Int), a1)) -> [a2] -> t a3 -> Maybe (HM.Matrix Double) +mkVecList mkAns dofList evalM = if null l + then Nothing + else Just $ HM.assoc (1,n) 0 $ map (fmap (\x -> fromRational $ toRational x / toRational maxVal)) l where l' = mapMaybe mkAns dofList l = runEval $ parListChunk 500 rdeepseq l' - lVals = map (\(x,y,z) -> z) l - max = maximum lVals + lVals = map (\((_,_),z) -> z) l + maxVal = maximum lVals n = length evalM - vecList = let vec = Sparse.fromList 1 n l in - if null l then Nothing else Just $ Sparse.scale (1/max) vec - -evalAnsatzEtaVecListEig :: [I.IntMap Int] -> AnsatzForestEta -> Maybe Sparse.SparseMatrixXd -evalAnsatzEtaVecListEig evalM EmptyForest = Nothing -evalAnsatzEtaVecListEig evalM f = mkVecList mkAns dofList evalM +evalAnsatzEtaVecListIncremental :: [I.IntMap Int] -> AnsatzForestEta -> Maybe (HM.Matrix Double) +evalAnsatzEtaVecListIncremental _ EmptyForest = Nothing +evalAnsatzEtaVecListIncremental evalM f = mkVecList mkAns dofList evalM where dofList = zip [0..] evalM mkAns (i,j) = let ansVal = eval1AnsatzForestEta j f - in if ansVal == 0 then Nothing else Just (0,i, fromIntegral ansVal) + in if ansVal == 0 then Nothing else Just ((0,i), ansVal) -evalAnsatzEpsilonVecListEig :: [I.IntMap Int] -> AnsatzForestEpsilon -> Maybe Sparse.SparseMatrixXd -evalAnsatzEpsilonVecListEig evalM f = if f == M.empty then Nothing else mkVecList mkAns dofList evalM +evalAnsatzEpsilonVecListIncremental :: [I.IntMap Int] -> AnsatzForestEpsilon -> Maybe (HM.Matrix Double) +evalAnsatzEpsilonVecListIncremental evalM f = if f == M.empty then Nothing else mkVecList mkAns dofList evalM where dofList = zip [0..] evalM mkAns (i,j) = let ansVal = eval1AnsatzForestEpsilon j f - in if ansVal == 0 then Nothing else Just (0,i, fromIntegral ansVal) + in if ansVal == 0 then Nothing else Just ((0,i), ansVal) --eval a given Forest for all inds @@ -1057,48 +1056,48 @@ evalAllEta :: [I.IntMap Int] -> AnsatzForestEta -> [[(Int,Int)]] -evalAllEta [] f = [] -evalAllEta evalMs EmptyForest = [] +evalAllEta [] _ = [] +evalAllEta _ EmptyForest = [] evalAllEta evalMs f = l' where - l = map (\x -> filter (\(a,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEta x f) evalMs + l = map (\x -> filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEta x f) evalMs l' = runEval $ parListChunk 500 rdeepseq l evalAllTensorEta :: (NFData a) => [(I.IntMap Int, a)] -> AnsatzForestEta -> AssocsList a -evalAllTensorEta [] f = [] -evalAllTensorEta evalMs EmptyForest = [] +evalAllTensorEta [] _ = [] +evalAllTensorEta _ EmptyForest = [] evalAllTensorEta evalMs f = l' where - l = map (\(x,z) -> (filter (\(a,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEta x f,z)) evalMs + l = map (\(x,z) -> (filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEta x f,z)) evalMs l' = runEval $ parListChunk 500 rdeepseq l evalAllEpsilon :: [I.IntMap Int] -> AnsatzForestEpsilon -> [[(Int,Int)]] -evalAllEpsilon [] f = [] +evalAllEpsilon [] _ = [] evalAllEpsilon evalMs f = if f == M.empty then [] else l' where - l = map (\x -> filter (\(a,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEpsilon x f) evalMs + l = map (\x -> filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEpsilon x f) evalMs l' = runEval $ parListChunk 500 rdeepseq l evalAllTensorEpsilon :: (NFData a) => [(I.IntMap Int, a)] -> AnsatzForestEpsilon -> AssocsList a -evalAllTensorEpsilon [] f = [] +evalAllTensorEpsilon [] _ = [] evalAllTensorEpsilon evalMs f = if f == M.empty then [] else l' where - l = map (\(x,z) -> ( filter (\(a,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEpsilon x f,z)) evalMs + l = map (\(x,z) -> ( filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEpsilon x f,z)) evalMs l' = runEval $ parListChunk 500 rdeepseq l evalAllTensorEtaAbs :: (NFData a) => [(I.IntMap Int, Int, a)] -> AnsatzForestEta -> AssocsListAbs a -evalAllTensorEtaAbs [] f = [] -evalAllTensorEtaAbs evalMs EmptyForest = [] +evalAllTensorEtaAbs [] _ = [] +evalAllTensorEtaAbs _ EmptyForest = [] evalAllTensorEtaAbs evalMs f = l' where - l = map (\(x,y,z) -> (filter (\(a,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEta x f, y,z)) evalMs + l = map (\(x,y,z) -> (filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEta x f, y,z)) evalMs l' = runEval $ parListChunk 500 rdeepseq l evalAllTensorEpsilonAbs :: (NFData a) => [(I.IntMap Int, Int, a)] -> AnsatzForestEpsilon -> AssocsListAbs a -evalAllTensorEpsilonAbs [] f = [] +evalAllTensorEpsilonAbs [] _ = [] evalAllTensorEpsilonAbs evalMs f = if f == M.empty then [] else l' where - l = map (\(x,y,z) -> ( filter (\(a,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEpsilon x f, y,z)) evalMs + l = map (\(x,y,z) -> ( filter (\(_,b) -> b /= 0) $ I.assocs $ evalAnsatzForestEpsilon x f, y,z)) evalMs l' = runEval $ parListChunk 500 rdeepseq l @@ -1114,66 +1113,63 @@ We start with the first way. --} -type RankDataEig = (Mat.MatrixXd, Sparse.SparseMatrixXd) +type RankDataIncremental = (HM.Matrix Double, HM.Matrix Double) -getVarNrEig :: RankDataEig -> Int -getVarNrEig = Sparse.rows . snd +getVarNrIncremental :: RankDataIncremental -> Int +getVarNrIncremental = HM.rows . snd --check in each step if the new ansatz vector is linear dependant w.r.t. the ansatz vectors obtained previously -checkNumericLinDepEig :: RankDataEig -> Maybe Sparse.SparseMatrixXd -> Maybe RankDataEig -checkNumericLinDepEig (lastMat, lastFullMat) (Just newVec) - | eigenRank < maxRank = Nothing +checkNumericLinDepIncremental :: RankDataIncremental -> Maybe (HM.Matrix Double) -> Maybe RankDataIncremental +checkNumericLinDepIncremental (lastMat, lastFullMat) (Just newVec) + | rk < maxRank = Nothing | otherwise = Just (newMat, newAnsatzMat) where - newVecTrans = Sparse.transpose newVec - scalar = Sparse.toMatrix $ Sparse.mul newVec newVecTrans - prodBlock = Sparse.toMatrix $ Sparse.mul lastFullMat newVecTrans - prodBlockTrans = Mat.transpose prodBlock - newMat = concatBlockMat lastMat prodBlock prodBlockTrans scalar - eigenRank = Sol.rank Sol.FullPivLU newMat - maxRank = min (Mat.cols newMat) (Mat.rows newMat) - newAnsatzMat = Sparse.fromRows $ Sparse.getRows lastFullMat ++ [newVec] -checkNumericLinDepEig (lastMat, lastFullMat) Nothing = Nothing - ---concat Matrices to a block Matrix - -concatBlockMat :: Mat.MatrixXd -> Mat.MatrixXd -> Mat.MatrixXd -> Mat.MatrixXd -> Mat.MatrixXd -concatBlockMat a b c d = newMat - where - newUpper = zipWith (++) (Mat.toList a) (Mat.toList b) - newLower = zipWith (++) (Mat.toList c) (Mat.toList d) - newMat = Mat.fromList $ newUpper ++ newLower + newVecTrans = HM.tr newVec + scalar = newVec Matrix.<> newVecTrans + prodBlock = lastFullMat Matrix.<> newVecTrans + prodBlockTrans = HM.tr prodBlock + newMat = HM.fromBlocks [[lastMat, prodBlock], + [prodBlockTrans, scalar ]] + rk = Matrix.rank newMat + maxRank = min (HM.cols newMat) (HM.rows newMat) + newAnsatzMat = lastFullMat HM.=== newVec +checkNumericLinDepIncremental _ Nothing = Nothing --in each step add the new AnsatzVector to the forest iff it is lin indep of the previous vectors {- alreadyPresentIO n total rDat - = putStrLn $ progress n total ++ " : " ++ "already present, not added, ansatz rank is " ++ show (getVarNrEig rDat) + = putStrLn $ progress n total ++ " : " ++ "already present, not added, ansatz rank is " ++ show (getVarNrIncremental rDat) notPresentNotAddedIO n total rDat - = putStrLn $ progress n total ++ " : " ++ "not present, linearly dependent, not added, ansatz rank is " ++ show (getVarNrEig rDat) + = putStrLn $ progress n total ++ " : " ++ "not present, linearly dependent, not added, ansatz rank is " ++ show (getVarNrIncremental rDat) notPresentAddedIO n total rDat - = putStrLn $ progress n total ++ " : " ++ "not present, linearly independent, added, ansatz rank is " ++ show (getVarNrEig rDat) + = putStrLn $ progress n total ++ " : " ++ "not present, linearly independent, added, ansatz rank is " ++ show (getVarNrIncremental rDat) progress n total = show n ++ " of " ++ show total -} +getNewRDat :: [I.IntMap Int] -> AnsatzForestEta -> RankDataIncremental -> Maybe RankDataIncremental getNewRDat evalM newAns rDat = newRDat where - newVec = evalAnsatzEtaVecListEig evalM newAns - newRDat = checkNumericLinDepEig rDat newVec + newVec = evalAnsatzEtaVecListIncremental evalM newAns + newRDat = checkNumericLinDepIncremental rDat newVec +getNewRDatEps :: [I.IntMap Int] -> AnsatzForestEpsilon -> RankDataIncremental -> Maybe RankDataIncremental getNewRDatEps evalM newAns rDat = newRDat where - newVec = evalAnsatzEpsilonVecListEig evalM newAns - newRDat = checkNumericLinDepEig rDat newVec + newVec = evalAnsatzEpsilonVecListIncremental evalM newAns + newRDat = checkNumericLinDepIncremental rDat newVec -getNewAns symList etaList rDat = symAnsatzForestEta symList $ mkForestFromAscList (etaList,Var 1 (getVarNrEig rDat + 1)) -getNewAnsEps symList epsList etaList rDat = symAnsatzForestEps symList $ mkForestFromAscListEpsilon (epsList,etaList,Var 1 (getVarNrEig rDat + 1)) +getNewAns :: Symmetry -> [Eta] -> RankDataIncremental -> AnsatzForestEta +getNewAns symList etaList rDat = symAnsatzForestEta symList $ mkForestFromAscList (etaList,Var 1 (getVarNrIncremental rDat + 1)) +getNewAnsEps :: Symmetry -> Epsilon -> [Eta] -> RankDataIncremental -> AnsatzForestEpsilon +getNewAnsEps symList epsList etaList rDat = symAnsatzForestEps symList $ mkForestFromAscListEpsilon (epsList,etaList,Var 1 (getVarNrIncremental rDat + 1)) + {- -addOrDiscardEtaEigIO :: Symmetry -> Int -> [I.IntMap Int] -> (AnsatzForestEta, RankDataEig) -> (Int, [Eta]) -> IO (AnsatzForestEta, RankDataEig) -addOrDiscardEtaEigIO symList len evalM (ans,rDat) (num, etaL) +addOrDiscardEtaIncrementalIO :: Symmetry -> Int -> [I.IntMap Int] -> (AnsatzForestEta, RankDataIncremental) -> (Int, [Eta]) -> IO (AnsatzForestEta, RankDataIncremental) +addOrDiscardEtaIncrementalIO symList len evalM (ans,rDat) (num, etaL) | isElem etaL ans = do alreadyPresentIO num len rDat return (ans,rDat) @@ -1190,8 +1186,8 @@ sumAns = addForests ans newAns -} -addOrDiscardEtaEig :: Symmetry -> [I.IntMap Int] -> (AnsatzForestEta, RankDataEig) -> [Eta] -> (AnsatzForestEta, RankDataEig) -addOrDiscardEtaEig symList evalM (ans,rDat) etaL +addOrDiscardEtaIncremental :: Symmetry -> [I.IntMap Int] -> (AnsatzForestEta, RankDataIncremental) -> [Eta] -> (AnsatzForestEta, RankDataIncremental) +addOrDiscardEtaIncremental symList evalM (ans,rDat) etaL | isElem etaL ans = (ans,rDat) | otherwise = case newRDat of Nothing -> (ans,rDat) @@ -1203,15 +1199,15 @@ {- -addOrDiscardEpsilonEigIO :: Symmetry -> Int -> [I.IntMap Int] -> (AnsatzForestEpsilon, RankDataEig) -> (Int,(Epsilon,[Eta])) -> IO (AnsatzForestEpsilon, RankDataEig) -addOrDiscardEpsilonEigIO symList len evalM (ans,rDat) (num,(epsL,etaL)) +addOrDiscardEpsilonIncrementalIO :: Symmetry -> Int -> [I.IntMap Int] -> (AnsatzForestEpsilon, RankDataIncremental) -> (Int,(Epsilon,[Eta])) -> IO (AnsatzForestEpsilon, RankDataIncremental) +addOrDiscardEpsilonIncrementalIO symList len evalM (ans,rDat) (num,(epsL,etaL)) | isElemEpsilon (epsL,etaL) ans = do alreadyPresentIO num len rDat return (ans,rDat) | otherwise = case newRDat of Nothing -> do notPresentNotAddedIO num len rDat - let r = getVarNrEig rDat + let r = getVarNrIncremental rDat return (ans,rDat) Just newRDat' -> do @@ -1223,8 +1219,8 @@ sumAns = addForestsEpsilon ans newAns -} -addOrDiscardEpsilonEig :: Symmetry -> [I.IntMap Int] -> (AnsatzForestEpsilon, RankDataEig) -> (Epsilon,[Eta]) -> (AnsatzForestEpsilon, RankDataEig) -addOrDiscardEpsilonEig symList evalM (ans,rDat) (epsL,etaL) +addOrDiscardEpsilonIncremental :: Symmetry -> [I.IntMap Int] -> (AnsatzForestEpsilon, RankDataIncremental) -> (Epsilon,[Eta]) -> (AnsatzForestEpsilon, RankDataIncremental) +addOrDiscardEpsilonIncremental symList evalM (ans,rDat) (epsL,etaL) | isElemEpsilon (epsL,etaL) ans = (ans,rDat) | otherwise = case newRDat of Nothing -> (ans,rDat) @@ -1238,97 +1234,97 @@ --construct the RankData from the first nonzero Ansatz {- -mk1stRankDataEtaEigIO :: Symmetry -> Int -> [(Int,[Eta])] -> [I.IntMap Int] -> IO (AnsatzForestEta,RankDataEig,[(Int,[Eta])]) -mk1stRankDataEtaEigIO symL numEta etaL evalM = +mk1stRankDataEtaIncrementalIO :: Symmetry -> Int -> [(Int,[Eta])] -> [I.IntMap Int] -> IO (AnsatzForestEta,RankDataIncremental,[(Int,[Eta])]) +mk1stRankDataEtaIncrementalIO symL numEta etaL evalM = do putStrLn $ show (fst $ head etaL) ++ " of " ++ show numEta let newAns = symAnsatzForestEta symL $ mkForestFromAscList (snd $ head etaL,Var 1 1) - let newVec = evalAnsatzEtaVecListEig evalM newAns + let newVec = evalAnsatzEtaVecListIncremental evalM newAns let restList = tail etaL case newVec of - Nothing -> if null restList then return (EmptyForest ,(Mat.fromList [], Sparse.fromList 0 0 []),[]) else mk1stRankDataEtaEigIO symL numEta restList evalM + Nothing -> if null restList then return (EmptyForest ,(HM.matrix 0 [], HM.matrix 0 []),[]) else mk1stRankDataEtaIncrementalIO symL numEta restList evalM Just newVec' -> return (newAns, (newMat, newVec'), restList) where - newVecTrans = Sparse.transpose newVec' - newMat = Sparse.toMatrix $ Sparse.mul newVec' newVecTrans + newVecTrans = HM.tr newVec' + newMat = newVec' Matrix.<> newVecTrans -} -mk1stRankDataEtaEig :: Symmetry -> [[Eta]] -> [I.IntMap Int] -> (AnsatzForestEta,RankDataEig,[[Eta]]) -mk1stRankDataEtaEig symL etaL evalM = output +mk1stRankDataEtaIncremental :: Symmetry -> [[Eta]] -> [I.IntMap Int] -> (AnsatzForestEta,RankDataIncremental,[[Eta]]) +mk1stRankDataEtaIncremental symL etaL evalM = output where newAns = symAnsatzForestEta symL $ mkForestFromAscList (head etaL,Var 1 1) - newVec = evalAnsatzEtaVecListEig evalM newAns + newVec = evalAnsatzEtaVecListIncremental evalM newAns restList = tail etaL output = case newVec of - Nothing -> if null restList then (EmptyForest,(Mat.fromList [], Sparse.fromList 0 0 []),[]) else mk1stRankDataEtaEig symL restList evalM + Nothing -> if null restList then (EmptyForest,(HM.matrix 0 [], HM.matrix 0 []),[]) else mk1stRankDataEtaIncremental symL restList evalM Just newVec' -> (newAns, (newMat, newVec'), restList) where - newVecTrans = Sparse.transpose newVec' - newMat = Sparse.toMatrix $ Sparse.mul newVec' newVecTrans + newVecTrans = HM.tr newVec' + newMat = newVec' Matrix.<> newVecTrans -mk1stRankDataEpsilonEig :: Symmetry -> [(Epsilon,[Eta])] -> [I.IntMap Int] -> (AnsatzForestEpsilon,RankDataEig,[(Epsilon,[Eta])]) -mk1stRankDataEpsilonEig symL epsL evalM = output +mk1stRankDataEpsilonIncremental :: Symmetry -> [(Epsilon,[Eta])] -> [I.IntMap Int] -> (AnsatzForestEpsilon,RankDataIncremental,[(Epsilon,[Eta])]) +mk1stRankDataEpsilonIncremental symL epsL evalM = output where newAns = symAnsatzForestEps symL $ mkForestFromAscListEpsilon (fst $ head epsL, snd $ head epsL,Var 1 1) - newVec = evalAnsatzEpsilonVecListEig evalM newAns + newVec = evalAnsatzEpsilonVecListIncremental evalM newAns restList = tail epsL output = case newVec of - Nothing -> if null restList then (M.empty,(Mat.fromList [], Sparse.fromList 0 0 []),[]) else mk1stRankDataEpsilonEig symL restList evalM + Nothing -> if null restList then (M.empty,(HM.matrix 0 [], HM.matrix 0 []),[]) else mk1stRankDataEpsilonIncremental symL restList evalM Just newVec' -> (newAns,(newMat, newVec'), restList) where - newVecTrans = Sparse.transpose newVec' - newMat = Sparse.toMatrix $ Sparse.mul newVec' newVecTrans + newVecTrans = HM.tr newVec' + newMat = newVec' Matrix.<> newVecTrans --finally reduce the ansatzList (IO versions print the current status for longer computations will follow with the next versions) -reduceAnsatzEtaEig :: Symmetry -> [[Eta]] -> [I.IntMap Int] -> (AnsatzForestEta,Sparse.SparseMatrixXd) -reduceAnsatzEtaEig symL etaL evalM - | null evalM = (EmptyForest, Sparse.fromList 0 0 []) - | null etaL = (EmptyForest, Sparse.fromList 0 0 []) +reduceAnsatzEtaIncremental :: Symmetry -> [[Eta]] -> [I.IntMap Int] -> (AnsatzForestEta, HM.Matrix Double) +reduceAnsatzEtaIncremental symL etaL evalM + | null evalM = (EmptyForest, HM.matrix 0 []) + | null etaL = (EmptyForest, HM.matrix 0 []) | otherwise = (finalForest, finalMat) where - (ans1,rDat1,restEtaL) = mk1stRankDataEtaEig symL etaL evalM - (finalForest, (_,finalMat)) = foldl' (addOrDiscardEtaEig symL evalM) (ans1,rDat1) restEtaL + (ans1,rDat1,restEtaL) = mk1stRankDataEtaIncremental symL etaL evalM + (finalForest, (_,finalMat)) = foldl' (addOrDiscardEtaIncremental symL evalM) (ans1,rDat1) restEtaL -reduceAnsatzEpsilonEig :: Symmetry -> [(Epsilon,[Eta])] -> [I.IntMap Int] -> (AnsatzForestEpsilon,Sparse.SparseMatrixXd) -reduceAnsatzEpsilonEig symL epsL evalM - | null evalM = (M.empty, Sparse.fromList 0 0 []) - | null epsL = (M.empty, Sparse.fromList 0 0 []) +reduceAnsatzEpsilonIncremental :: Symmetry -> [(Epsilon,[Eta])] -> [I.IntMap Int] -> (AnsatzForestEpsilon, HM.Matrix Double) +reduceAnsatzEpsilonIncremental symL epsL evalM + | null evalM = (M.empty, HM.matrix 0 []) + | null epsL = (M.empty, HM.matrix 0 []) | otherwise = (finalForest, finalMat) where - (ans1,rDat1,restEpsL) = mk1stRankDataEpsilonEig symL epsL evalM - (finalForest, (_,finalMat)) = foldl' (addOrDiscardEpsilonEig symL evalM) (ans1,rDat1) restEpsL + (ans1,rDat1,restEpsL) = mk1stRankDataEpsilonIncremental symL epsL evalM + (finalForest, (_,finalMat)) = foldl' (addOrDiscardEpsilonIncremental symL evalM) (ans1,rDat1) restEpsL --construct a basis ansatz forest -getEtaForestEig :: Int -> Symmetry -> [I.IntMap Int] -> (AnsatzForestEta,Sparse.SparseMatrixXd) -getEtaForestEig ord sym [] = (EmptyForest, Sparse.fromList 0 0 []) -getEtaForestEig ord sym evalMs - | null allEtaLists = (EmptyForest, Sparse.fromList 0 0 []) - | otherwise = reduceAnsatzEtaEig sym allEtaLists evalMs +getEtaForestIncremental :: Int -> Symmetry -> [I.IntMap Int] -> (AnsatzForestEta, HM.Matrix Double) +getEtaForestIncremental _ _ [] = (EmptyForest, HM.matrix 0 []) +getEtaForestIncremental ord sym evalMs + | null allEtaLists = (EmptyForest, HM.matrix 0 []) + | otherwise = reduceAnsatzEtaIncremental sym allEtaLists evalMs where allInds = getEtaInds [1..ord] sym allEtaLists = map mkEtaList allInds -getEpsForestEig :: Int -> Symmetry -> [I.IntMap Int] -> (AnsatzForestEpsilon,Sparse.SparseMatrixXd) -getEpsForestEig ord sym [] = (M.empty, Sparse.fromList 0 0 []) -getEpsForestEig ord sym evalMs - | null allEpsLists = (M.empty, Sparse.fromList 0 0 []) - | otherwise = reduceAnsatzEpsilonEig sym allEpsLists evalMs +getEpsForestIncremental :: Int -> Symmetry -> [I.IntMap Int] -> (AnsatzForestEpsilon, HM.Matrix Double) +getEpsForestIncremental _ _ [] = (M.empty, HM.matrix 0 []) +getEpsForestIncremental ord sym evalMs + | null allEpsLists = (M.empty, HM.matrix 0 []) + | otherwise = reduceAnsatzEpsilonIncremental sym allEpsLists evalMs where allInds = getEpsilonInds [1..ord] sym allEpsLists = map mkEpsilonList allInds --eta and eps forest combined -getFullForestEig :: Int -> Symmetry -> [I.IntMap Int] -> [I.IntMap Int] -> (AnsatzForestEta, AnsatzForestEpsilon, Sparse.SparseMatrixXd, Sparse.SparseMatrixXd) -getFullForestEig ord sym evalMEta evalMEps = (etaAns, epsAns, etaMat, epsMat) +getFullForestIncremental :: Int -> Symmetry -> [I.IntMap Int] -> [I.IntMap Int] -> (AnsatzForestEta, AnsatzForestEpsilon, HM.Matrix Double, HM.Matrix Double) +getFullForestIncremental ord sym evalMEta evalMEps = (etaAns, epsAns, etaMat, epsMat) where - (etaAns,etaMat) = getEtaForestEig ord sym evalMEta - (epsAns',epsMat) = getEpsForestEig ord sym evalMEps + (etaAns,etaMat) = getEtaForestIncremental ord sym evalMEta + (epsAns',epsMat) = getEpsForestIncremental ord sym evalMEps epsAns = relabelAnsatzForestEpsilon (1 + length (getForestLabels etaAns)) epsAns' {-- @@ -1340,7 +1336,7 @@ --} evalToTensSym :: Symmetry -> [(I.IntMap Int, IndTupleST n1 0)] -> [(I.IntMap Int, IndTupleST n1 0)] -> AnsatzForestEta -> AnsatzForestEpsilon -> STTens n1 0 AnsVarR -evalToTensSym (p,ap,b,c,bc) evalEta evalEps ansEta ansEps = symTens +evalToTensSym (p,ap,b,c,bc) evalEta evalEps ansEta ansEps = symT where p' = map (\(x,y) -> (x-1,y-1)) p ap' = map (\(x,y) -> (x-1,y-1)) ap @@ -1348,7 +1344,7 @@ c' = map (map (subtract 1)) c bc' = map (map (map (subtract 1))) bc tens = evalToTens evalEta evalEps ansEta ansEps - symTens = foldr cyclicBlockSymATens1 ( + symT = foldr cyclicBlockSymATens1 ( foldr cyclicSymATens1 ( foldr symBlockATens1 ( foldr aSymATens1 ( @@ -1389,10 +1385,10 @@ mkEvalMaps :: [[Int]] -> [I.IntMap Int] mkEvalMaps l = let s = length (head l) in map (mkEvalMap s) l -mkEvalMapsInds :: forall (n :: Nat). SingI n => [[Int]] -> [(I.IntMap Int, IndTupleST n 0)] -mkEvalMapsInds l = let s = length (head l) in map (\x -> (mkEvalMap s x, (fromList $ map toEnum x, Empty))) l +mkEvalMapsInds :: forall (n :: Nat). KnownNat n => [[Int]] -> [(I.IntMap Int, IndTupleST n 0)] +mkEvalMapsInds l = let s = length (head l) in map (\x -> (mkEvalMap s x, (fromListUnsafe $ map toEnum x, Empty))) l -mkAllEvalMaps :: forall (n :: Nat). SingI n => Symmetry -> [[Int]] -> ([I.IntMap Int], [I.IntMap Int], [(I.IntMap Int, IndTupleST n 0)], [(I.IntMap Int, IndTupleST n 0)]) +mkAllEvalMaps :: forall (n :: Nat). KnownNat n => Symmetry -> [[Int]] -> ([I.IntMap Int], [I.IntMap Int], [(I.IntMap Int, IndTupleST n 0)], [(I.IntMap Int, IndTupleST n 0)]) mkAllEvalMaps sym l = (evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) where evalLEta = filter isEtaList l @@ -1420,64 +1416,52 @@ evalMEpsInds = map (\(x,y,z) -> (mkEvalMap ord x, y, z)) evalLEps -- | The function is similar to @'mkAnsatzTensorFastSym'@ yet it uses an algorithm that prioritizes memory usage over fast computation times. -mkAnsatzTensorEigSym :: forall (n :: Nat). SingI n => Int -> Symmetry -> [[Int]] -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) -mkAnsatzTensorEigSym ord symmetries evalL = (ansEta, ansEps, tens) +mkAnsatzTensorIncrementalSym :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> [[Int]] -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) +mkAnsatzTensorIncrementalSym ord symmetries evalL = (ansEta, ansEps, tens) where (evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) = mkAllEvalMaps symmetries evalL - (ansEta, ansEps, _, _) = getFullForestEig ord symmetries evalMEtaRed evalMEpsRed + (ansEta, ansEps, _, _) = getFullForestIncremental ord symmetries evalMEtaRed evalMEpsRed tens = evalToTensSym symmetries evalMEtaInds evalMEpsInds ansEta ansEps -- | The function is similar to @'mkAnsatzTensorFast'@ yet it uses an algorithm that prioritizes memory usage over fast computation times. -mkAnsatzTensorEig :: forall (n :: Nat). SingI n => Int -> Symmetry -> [[Int]] -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) -mkAnsatzTensorEig ord symmetries evalL = (ansEta, ansEps, tens) +mkAnsatzTensorIncremental :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> [[Int]] -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) +mkAnsatzTensorIncremental ord symmetries evalL = (ansEta, ansEps, tens) where (evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) = mkAllEvalMaps symmetries evalL - (ansEta, ansEps, _, _) = getFullForestEig ord symmetries evalMEtaRed evalMEpsRed + (ansEta, ansEps, _, _) = getFullForestIncremental ord symmetries evalMEtaRed evalMEpsRed tens = evalToTens evalMEtaInds evalMEpsInds ansEta ansEps -- | The function is similar to @'mkAnsatzTensorFastAbs'@ yet it uses an algorithm that prioritizes memory usage over fast computation times. -mkAnsatzTensorEigAbs :: Int -> Symmetry -> [([Int], Int, [IndTupleAbs n1 0 n2 0 n3 0])] -> (AnsatzForestEta, AnsatzForestEpsilon, ATens n1 0 n2 0 n3 0 AnsVarR) -mkAnsatzTensorEigAbs ord symmetries evalL = (ansEta, ansEps, tens) +mkAnsatzTensorIncrementalAbs :: Int -> Symmetry -> [([Int], Int, [IndTupleAbs n1 0 n2 0 n3 0])] -> (AnsatzForestEta, AnsatzForestEpsilon, ATens n1 0 n2 0 n3 0 AnsVarR) +mkAnsatzTensorIncrementalAbs ord symmetries evalL = (ansEta, ansEps, tens) where (evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) = mkAllEvalMapsAbs symmetries evalL - (ansEta, ansEps, _, _) = getFullForestEig ord symmetries evalMEtaRed evalMEpsRed + (ansEta, ansEps, _, _) = getFullForestIncremental ord symmetries evalMEtaRed evalMEpsRed tens = evalToTensAbs evalMEtaInds evalMEpsInds ansEta ansEps --now we start with the second way -assocsToEig :: [[(Int,Int)]] -> Mat.MatrixXd -assocsToEig l = Sparse.toMatrix $ Sparse.fromList n m l' +assocsToMat :: [[(Int,Int)]] -> HM.Matrix Double +assocsToMat l = HM.assoc (m,n) 0 l' where - l' = concat $ zipWith (\r z -> map (\(x,y) -> (z-1, x-1, fromIntegral y)) r) l [1..] - n = maximum (map (\(x,_,_) -> x) l') + 1 - m = maximum (map (\(_,x,_) -> x) l') + 1 + l' = concat $ zipWith (\r z -> map (\(x,y) -> ((z-1, x-1), fromIntegral y)) r) l [1..] + sparse = M.fromList l' + m = maximum (map (\((x,_),_) -> x) l') + 1 + n = maximum (map (\((_,x),_) -> x) l') + 1 --filter the lin. dependant vars from the Assocs List -{- -optimized version, requires custom eigen build -getPivots' :: [[(Int,Int)]] -> [Int] -getPivots' l = map (1+) p - where - mat = assocsToEig l - p = Sol.pivots Sol.FullPivLU mat --} - getPivots :: [[(Int,Int)]] -> [Int] -getPivots l = map (1+) p +getPivots matList = map (1+) pivots where - mat = assocsToEig l - pMatTr = Mat.toList $ Mat.transpose $ Sol.image Sol.FullPivLU mat - matTr = Mat.toList $ Mat.transpose mat - p = mapMaybe (`elemIndex` matTr) pMatTr - - + mat = assocsToMat matList + pivots = independentColumns mat --reduce linear deps in the ansätze -reduceLinDepsFastEta :: [I.IntMap Int] -> Symmetry -> AnsatzForestEta -> AnsatzForestEta -reduceLinDepsFastEta evalM symL ansEta = newEtaAns +reduceLinDepsFastEta :: [I.IntMap Int] -> AnsatzForestEta -> AnsatzForestEta +reduceLinDepsFastEta evalM ansEta = newEtaAns where etaL = evalAllEta evalM ansEta etaVars = getPivots etaL @@ -1485,8 +1469,8 @@ remVarsEta = allEtaVars \\ etaVars newEtaAns = relabelAnsatzForest 1 $ removeVarsEta remVarsEta ansEta -reduceLinDepsFastEps :: [I.IntMap Int] -> Symmetry -> AnsatzForestEpsilon -> AnsatzForestEpsilon -reduceLinDepsFastEps evalM symL ansEps = newEpsAns +reduceLinDepsFastEps :: [I.IntMap Int] -> AnsatzForestEpsilon -> AnsatzForestEpsilon +reduceLinDepsFastEps evalM ansEps = newEpsAns where epsL = evalAllEpsilon evalM ansEps epsVars = getPivots epsL @@ -1501,15 +1485,15 @@ where ansEta = getEtaForestFast ord symmetries ansEpsilon = getEpsForestFast ord symmetries - ansEtaRed = reduceLinDepsFastEta evalMEtaRed symmetries ansEta - ansEpsRed' = reduceLinDepsFastEps evalMEpsRed symmetries ansEpsilon + ansEtaRed = reduceLinDepsFastEta evalMEtaRed ansEta + ansEpsRed' = reduceLinDepsFastEps evalMEpsRed ansEpsilon ansEpsRed = relabelAnsatzForestEpsilon (1 + length (getForestLabels ansEtaRed)) ansEpsRed' -- | The function computes all linear independent ansätze that have rank specified by the first integer argument and further satisfy the symmetry specified by the @'Symmetry'@ value. -- The additional argument of type @[['Int']]@ is used to provide the information of all (by means of the symmetry at hand) independent components of the ansätze. -- Explicit examples how this information can be computed are provided by the functions for @'areaList4'@, ... and also by @'metricList2'@, ... . -- The output is given as spacetime tensor @'STTens'@ and is explicitly symmetrized. -mkAnsatzTensorFastSym :: forall (n :: Nat). SingI n => Int -> Symmetry -> [[Int]]-> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) +mkAnsatzTensorFastSym :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> [[Int]]-> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) mkAnsatzTensorFastSym ord symmetries evalL = (ansEta, ansEps, tens) where (evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) = mkAllEvalMaps symmetries evalL @@ -1521,7 +1505,7 @@ -- | This function provides the same functionality as @'mkAnsatzTensorFast'@ but without explicit symmetrization of the result. In other words from each symmetrization sum only the first -- summand is returned. This is advantageous as for large expressions explicit symmetrization might be expensive and further is sometime simply not needed as the result might for instance be contracted against -- a symmetric object, which thus enforces the symmetry, in further steps of the computation. -mkAnsatzTensorFast :: forall (n :: Nat). SingI n => Int -> Symmetry -> [[Int]]-> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) +mkAnsatzTensorFast :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> [[Int]]-> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) mkAnsatzTensorFast ord symmetries evalL = (ansEta, ansEps, tens) where (evalMEtaRed, evalMEpsRed, evalMEtaInds, evalMEpsInds) = mkAllEvalMaps symmetries evalL @@ -1577,16 +1561,19 @@ filterASym :: [Int] -> (Int,Int) -> Bool filterASym inds (i,j) = (inds !! (i-1)) < (inds !! (j-1)) +getPairs :: [a] -> [(a, a)] +getPairs [a,b] = [(a,b)] +getPairs (x:xs) = (x, head xs) : getPairs xs +getPairs _ = error "invalid index combination" + filterCSym :: [Int] -> [Int] -> Bool filterCSym inds i = and boolL where - getPairs [a,b] = [(a,b)] - getPairs (x:xs) = (x, head xs) : getPairs xs pairL = getPairs i boolL = map (filterPSym inds) pairL filterBSym :: [Int] -> ([Int],[Int]) -> Bool -filterBSym inds ([],[]) = True +filterBSym _ ([],[]) = True filterBSym inds (x:xs,y:ys) | xVal < yVal = True | xVal == yVal = filterBSym inds (xs,ys) @@ -1594,12 +1581,11 @@ where xVal = inds !! (x-1) yVal = inds !! (y-1) +filterBSym _ _ = error "cannot non-empty list w.r.t. empty symmetries" filterBCSym :: [Int] -> [[Int]] -> Bool filterBCSym inds i = and boolL where - getPairs [a,b] = [(a,b)] - getPairs (x:xs) = (x, head xs) : getPairs xs pairL = getPairs i boolL = map (filterBSym inds) pairL @@ -1669,6 +1655,7 @@ swapBlocks (m1,m2) x = let m = I.fromList $ zip m1 m2 ++ zip m2 m1 in fromMaybe x $ I.lookup x m newMap = canonicalizeBlockPair (is,js) iMap +canonicalizeBlockPair _ _ = error "invalid index combination" canonicalizeIntMap :: Symmetry -> I.IntMap Int -> I.IntMap Int @@ -1683,11 +1670,11 @@ canonicalizeList sym inds = I.elems $ canonicalizeIntMap sym $ I.fromList $ zip [1..] inds allList' :: Int -> [(Int,Int)] -> [(Int,Int)] -> [(Int,Int)] -> [(Int,Int)] -> [[Int]] -allList' 1 syms aSyms symBounds aSymBounds = case (symB, aSymB) of +allList' 1 _ _ symBounds aSymBounds = case (symB, aSymB) of (Just j, Nothing) -> [[k] | k <- [j..3]] (Nothing, Just j) -> [[k] | k <- [j+1..3]] (Nothing, Nothing) -> [[0], [1], [2], [3]] - (Just j, Just k) -> [[k] | k <- [max j (k+1) .. 3]] + (Just j, Just k) -> [[k'] | k' <- [max j (k+1) .. 3]] where (symB,aSymB) = (lookup 1 symBounds, lookup 1 aSymBounds) allList' i syms aSyms symBounds aSymBounds = concatMap (\x -> (:) <$> [x] <*> allList' (i-1) newSyms newASyms (newSymBounds x) (newASymBounds x)) l @@ -1720,14 +1707,14 @@ --use the above functions to construct ansätze without providing eval lists by hand -- | The function is similar to @'mkAnsatzTensorFastSym''@ yet it uses an algorithm that prioritizes memory usage over fast computation times. -mkAnsatzTensorEigSym' :: forall (n :: Nat). SingI n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) -mkAnsatzTensorEigSym' ord symmetries = mkAnsatzTensorEigSym ord symmetries evalL +mkAnsatzTensorIncrementalSym' :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) +mkAnsatzTensorIncrementalSym' ord symmetries = mkAnsatzTensorIncrementalSym ord symmetries evalL where evalL = filter (`filterAllSym` symmetries) $ allList ord symmetries -- | Provides the same functionality as @'mkAnsatzTensorFastSym'@ with the difference that the list of independent index combinations is automatically computed form the present symmetry. -- Note that this yields slightly higher computation costs. -mkAnsatzTensorFastSym' :: forall (n :: Nat). SingI n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) +mkAnsatzTensorFastSym' :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) mkAnsatzTensorFastSym' ord symmetries = mkAnsatzTensorFastSym ord symmetries evalL where evalL = filter (`filterAllSym` symmetries) $ allList ord symmetries @@ -1735,14 +1722,14 @@ --and without explicit symmetrization -- | The function is similar to @'mkAnsatzTensorFast''@ yet it uses an algorithm that prioritizes memory usage over fast computation times. -mkAnsatzTensorEig' :: forall (n :: Nat). SingI n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) -mkAnsatzTensorEig' ord symmetries = mkAnsatzTensorEig ord symmetries evalL +mkAnsatzTensorIncremental' :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) +mkAnsatzTensorIncremental' ord symmetries = mkAnsatzTensorIncremental ord symmetries evalL where evalL = filter (`filterAllSym` symmetries) $ allList ord symmetries -- | Provides the same functionality as @'mkAnsatzTensorFast'@ with the difference that the list of independent index combinations is automatically computed form the present symmetry. -- Note that this yields slightly higher computation costs. -mkAnsatzTensorFast' :: forall (n :: Nat). SingI n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) +mkAnsatzTensorFast' :: forall (n :: Nat). KnownNat n => Int -> Symmetry -> (AnsatzForestEta, AnsatzForestEpsilon, STTens n 0 AnsVarR) mkAnsatzTensorFast' ord symmetries = mkAnsatzTensorFast ord symmetries evalL where evalL = filter (`filterAllSym` symmetries) $ allList ord symmetries @@ -1774,9 +1761,11 @@ areaMult [a,b,c,d] | a == c && b == d = 4 | otherwise = 8 +areaMult _ = error "expected four indices" iMult2 :: [Int] -> Int iMult2 [p,q] = if p == q then 1 else 2 +iMult2 _ = error "expected two indices" --Area metric eval lists @@ -1799,17 +1788,15 @@ areaList8 :: [([Int], Int, [IndTupleAbs 2 0 0 0 0 0])] areaList8 = list where - trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b') = ((I.!) trianArea a, (I.!) trianArea b) in (a' ++ b', areaMult a' * areaMult b', map (\[a,b] -> (Append (Ind20 $ a-1) $ singletonInd (Ind20 $ b-1), Empty, Empty, Empty, Empty, Empty)) $ nub $ permutations [a,b] ) | a <- [1..21], b <- [a..21]] + list = [ let (a',b') = ((I.!) trianArea a, (I.!) trianArea b) in (a' ++ b', areaMult a' * areaMult b', map (\[_a,_b] -> (Append (Ind20 $ _a-1) $ singletonInd (Ind20 $ _b-1), Empty, Empty, Empty, Empty, Empty)) $ nub $ permutations [a,b] ) | a <- [1..21], b <- [a..21]] -- | Evaluation list for \(a^{Ap Bq}\). Note that also when using the abstract indices this ansatz still features the \( (Ap) \leftrightarrow (Bq) \) symmetry. areaList10_1 :: [([Int], Int, [IndTupleAbs 2 0 0 0 2 0])] areaList10_1 = list where - trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b') = ((I.!) trianArea a, (I.!) trianArea b) in (a' ++ p : b' ++ [q], areaMult a' * areaMult b', map (\[[a,p],[b,q]] -> (Append (Ind20 $ a-1) $ singletonInd (Ind20 $ b-1), Empty, Empty, Empty, Append (Ind3 p) $ singletonInd (Ind3 q), Empty)) $ nub $ permutations [[a,p],[b,q]]) | a <- [1..21], b <- [a..21], p <- [0..3], q <- [0..3], not (a==b && p>q)] + list = [ let (a',b') = ((I.!) trianArea a, (I.!) trianArea b) in (a' ++ p : b' ++ [q], areaMult a' * areaMult b', map (\[[_a,_p],[_b,_q]] -> (Append (Ind20 $ _a-1) $ singletonInd (Ind20 $ _b-1), Empty, Empty, Empty, Append (Ind3 _p) $ singletonInd (Ind3 _q), Empty)) $ nub $ permutations [[a,p],[b,q]]) | a <- [1..21], b <- [a..21], p <- [0..3], q <- [0..3], not (a==b && p>q)] -- | Evaluation list for \(a^{ABI} \). areaList10_2 :: [([Int], Int, [IndTupleAbs 2 0 1 0 0 0])] @@ -1823,9 +1810,8 @@ areaList12 :: [([Int], Int, [IndTupleAbs 3 0 0 0 0 0])] areaList12 = list where - trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b',c') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c) in (a' ++ b' ++ c', areaMult a' * areaMult b' * areaMult c', map (\[a,b,c] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ b-1) $ singletonInd (Ind20 $ c-1), Empty, Empty, Empty, Empty, Empty)) $ nub $ permutations [a,b,c] )| a <- [1..21], b <- [a..21], c <- [b..21] ] + list = [ let (a',b',c') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c) in (a' ++ b' ++ c', areaMult a' * areaMult b' * areaMult c', map (\[_a,_b,_c] -> (Append (Ind20 $ _a-1) $ Append (Ind20 $ _b-1) $ singletonInd (Ind20 $ _c-1), Empty, Empty, Empty, Empty, Empty)) $ nub $ permutations [a,b,c] )| a <- [1..21], b <- [a..21], c <- [b..21] ] --AI:BJ areaList12_1 :: [([Int], Int, [IndTupleAbs 2 0 2 0 0 0])] @@ -1833,15 +1819,14 @@ where trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',i',b',j') = ((I.!) trianArea a, (I.!) trian2 i, (I.!) trianArea b, (I.!) trian2 j) in (a' ++ i' ++ b' ++ j' , areaMult a' * areaMult b' * iMult2 i' * iMult2 j', map (\[[a,i],[b,j]] -> (Append (Ind20 $ a-1) $ singletonInd (Ind20 $ b-1), Empty, Append (Ind9 $ i-1) $ singletonInd (Ind9 $ j-1), Empty, Empty, Empty)) $ nub $ permutations [[a,i],[b,j]] ) | a <- [1..21], b <- [a..21], i <- [1..10], j <- [1..10], not (a==b && i>j) ] + list = [ let (a',i',b',j') = ((I.!) trianArea a, (I.!) trian2 i, (I.!) trianArea b, (I.!) trian2 j) in (a' ++ i' ++ b' ++ j' , areaMult a' * areaMult b' * iMult2 i' * iMult2 j', map (\[[_a,_i],[_b,_j]] -> (Append (Ind20 $ _a-1) $ singletonInd (Ind20 $ _b-1), Empty, Append (Ind9 $ _i-1) $ singletonInd (Ind9 $ _j-1), Empty, Empty, Empty)) $ nub $ permutations [[a,i],[b,j]] ) | a <- [1..21], b <- [a..21], i <- [1..10], j <- [1..10], not (a==b && i>j) ] -- | Evaluation list for \(a^{ABp Cq}\). Note that also when using the abstract indices this ansatz still features the \( (Bp) \leftrightarrow (Cq) \) symmetry. areaList14_1 :: [([Int], Int, [IndTupleAbs 3 0 0 0 2 0])] areaList14_1 = list where - trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b',c') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c) in (a' ++ b' ++ p : c' ++ [q], areaMult a' * areaMult b' * areaMult c', map (\[[b,p],[c,q]] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ b-1) $ singletonInd (Ind20 $ c-1), Empty, Empty, Empty, Append (Ind3 p) $ singletonInd (Ind3 q), Empty)) $ nub $ permutations [[b,p],[c,q]]) | a <- [1..21], b <- [1..21], c <- [b..21], p <- [0..3], q <- [0..3], not (b==c && p>q) ] + list = [ let (a',b',c') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c) in (a' ++ b' ++ p : c' ++ [q], areaMult a' * areaMult b' * areaMult c', map (\[[_b,_p],[_c,_q]] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ _b-1) $ singletonInd (Ind20 $ _c-1), Empty, Empty, Empty, Append (Ind3 _p) $ singletonInd (Ind3 _q), Empty)) $ nub $ permutations [[b,p],[c,q]]) | a <- [1..21], b <- [1..21], c <- [b..21], p <- [0..3], q <- [0..3], not (b==c && p>q) ] -- | Evaluation list for \(a^{A B C I}\). Note that also when using the abstract indices this ansatz still features the \( (A) \leftrightarrow (B) \) symmetry. areaList14_2 :: [([Int], Int, [IndTupleAbs 3 0 1 0 0 0])] @@ -1849,7 +1834,7 @@ where trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b',c',i') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trian2 i) in ( a' ++ b' ++ c' ++ i', areaMult a' * areaMult b' * areaMult c' * iMult2 i', map (\[a,b] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ b-1) $ singletonInd (Ind20 $ c-1), Empty, singletonInd (Ind9 $ i-1), Empty, Empty, Empty)) $ nub $ permutations [a,b] ) | a <- [1..21], b <- [a..21], c <- [1..21], i <- [1..10] ] + list = [ let (a',b',c',i') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trian2 i) in ( a' ++ b' ++ c' ++ i', areaMult a' * areaMult b' * areaMult c' * iMult2 i', map (\[_a,_b] -> (Append (Ind20 $ _a-1) $ Append (Ind20 $ _b-1) $ singletonInd (Ind20 $ c-1), Empty, singletonInd (Ind9 $ i-1), Empty, Empty, Empty)) $ nub $ permutations [a,b] ) | a <- [1..21], b <- [a..21], c <- [1..21], i <- [1..10] ] --Ap:Bq:CI areaList16_1 :: [([Int], Int, [IndTupleAbs 3 0 1 0 2 0])] @@ -1857,7 +1842,7 @@ where trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b',c',i') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trian2 i) in (a' ++ p : b' ++ q : c' ++ i' , areaMult a' * areaMult b' * areaMult c' * iMult2 i', map (\[[a,p],[b,q]] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ b-1) $ singletonInd (Ind20 $ c-1), Empty, singletonInd (Ind9 $ i-1), Empty, Append (Ind3 p) $ singletonInd (Ind3 q), Empty)) $ nub $ permutations [[a,p],[b,q]]) | a <- [1..21], b <- [a..21], c <- [1..21], i <- [1..10], p <- [0..3], q <- [0..3], not (a==b && p>q) ] + list = [ let (a',b',c',i') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trian2 i) in (a' ++ p : b' ++ q : c' ++ i' , areaMult a' * areaMult b' * areaMult c' * iMult2 i', map (\[[_a,_p],[_b,_q]] -> (Append (Ind20 $ _a-1) $ Append (Ind20 $ _b-1) $ singletonInd (Ind20 $ c-1), Empty, singletonInd (Ind9 $ i-1), Empty, Append (Ind3 _p) $ singletonInd (Ind3 _q), Empty)) $ nub $ permutations [[a,p],[b,q]]) | a <- [1..21], b <- [a..21], c <- [1..21], i <- [1..10], p <- [0..3], q <- [0..3], not (a==b && p>q) ] --A:BI:CJ areaList16_2 :: [([Int], Int, [IndTupleAbs 3 0 2 0 0 0])] @@ -1865,7 +1850,7 @@ where trian2 = trianMap2 trianArea = trianMapArea - list = [let (a',b',c',i', j') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trian2 i, (I.!) trian2 j) in (a' ++ b' ++ i' ++ c' ++ j', areaMult a' * areaMult b' * areaMult c' * iMult2 i' * iMult2 j', map (\[[b,i],[c,j]] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ b-1) $ singletonInd (Ind20 $ c-1), Empty, Append (Ind9 $ i-1) $ singletonInd (Ind9 $ j-1), Empty, Empty, Empty) ) $ nub $ permutations [[b,i],[c,j]])| a <- [1..21], b <- [1..21], c <- [b..21], i <- [1..10], j <- [1..10], not (b==c && i>j)] + list = [let (a',b',c',i', j') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trian2 i, (I.!) trian2 j) in (a' ++ b' ++ i' ++ c' ++ j', areaMult a' * areaMult b' * areaMult c' * iMult2 i' * iMult2 j', map (\[[_b,_i],[_c,_j]] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ _b-1) $ singletonInd (Ind20 $ _c-1), Empty, Append (Ind9 $ _i-1) $ singletonInd (Ind9 $ _j-1), Empty, Empty, Empty) ) $ nub $ permutations [[b,i],[c,j]])| a <- [1..21], b <- [1..21], c <- [b..21], i <- [1..10], j <- [1..10], not (b==c && i>j)] --AI:BJ:CK areaList18 :: [([Int], Int, [IndTupleAbs 3 0 3 0 0 0])] @@ -1873,7 +1858,7 @@ where trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b',c',i', j', k') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trian2 i, (I.!) trian2 j, (I.!) trian2 k) in (a' ++ i' ++ b' ++ j' ++ c' ++ k', areaMult a' * areaMult b' * areaMult c' * iMult2 i' * iMult2 j' * iMult2 k', map (\[[a,i],[b,j],[c,k]] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ b-1) $ singletonInd (Ind20 $ c-1), Empty, Append (Ind9 $ i-1) $ Append (Ind9 $ j-1) $ singletonInd (Ind9 $ k-1), Empty, Empty, Empty) ) $ nub $ permutations [[a,i],[b,j],[c,k]]) | a <- [1..21], b <- [a..21], c <- [b..21], i <- [1..10], j <- [1..10], not (a==b && i>j), k <- [1..10], not (b==c && j>k) ] + list = [ let (a',b',c',i', j', k') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trian2 i, (I.!) trian2 j, (I.!) trian2 k) in (a' ++ i' ++ b' ++ j' ++ c' ++ k', areaMult a' * areaMult b' * areaMult c' * iMult2 i' * iMult2 j' * iMult2 k', map (\[[_a,_i],[_b,_j],[_c,_k]] -> (Append (Ind20 $ _a-1) $ Append (Ind20 $ _b-1) $ singletonInd (Ind20 $ _c-1), Empty, Append (Ind9 $ _i-1) $ Append (Ind9 $ _j-1) $ singletonInd (Ind9 $ _k-1), Empty, Empty, Empty) ) $ nub $ permutations [[a,i],[b,j],[c,k]]) | a <- [1..21], b <- [a..21], c <- [b..21], i <- [1..10], j <- [1..10], not (a==b && i>j), k <- [1..10], not (b==c && j>k) ] --order 4 @@ -1881,9 +1866,8 @@ areaList16 :: [([Int], Int, [IndTupleAbs 4 0 0 0 0 0])] areaList16 = list where - trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b',c', d') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trianArea d) in (a' ++ b' ++ c' ++ d', areaMult a' * areaMult b' * areaMult c' * areaMult d', map (\[a,b,c,d] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ b-1) $ Append (Ind20 $ c-1) $ singletonInd (Ind20 $ d-1), Empty, Empty, Empty, Empty, Empty)) $ nub $ permutations [a,b,c,d] )| a <- [1..21], b <- [a..21], c <- [b..21], d <- [c..21] ] + list = [ let (a',b',c', d') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trianArea d) in (a' ++ b' ++ c' ++ d', areaMult a' * areaMult b' * areaMult c' * areaMult d', map (\[_a,_b,_c,_d] -> (Append (Ind20 $ _a-1) $ Append (Ind20 $ _b-1) $ Append (Ind20 $ _c-1) $ singletonInd (Ind20 $ _d-1), Empty, Empty, Empty, Empty, Empty)) $ nub $ permutations [a,b,c,d] )| a <- [1..21], b <- [a..21], c <- [b..21], d <- [c..21] ] --A:B:C:DI areaList18_2 :: [( [Int], Int, [IndTupleAbs 4 0 1 0 0 0])] @@ -1891,24 +1875,22 @@ where trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b',c',d',i') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trianArea d, (I.!) trian2 i) in (a' ++ b' ++ c'++d'++i', areaMult a' * areaMult b' * areaMult c' * areaMult d' * iMult2 i', map (\[a,b,c] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ b-1) $ Append (Ind20 $ c-1) (singletonInd (Ind20 $ d-1)), Empty, singletonInd (Ind9 $ i-1), Empty, Empty, Empty) ) $ nub $ permutations [a,b,c] ) | a <- [1..21], b <- [a..21], c <- [b..21], d <- [1..21], i <- [1..10] ] + list = [ let (a',b',c',d',i') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trianArea d, (I.!) trian2 i) in (a' ++ b' ++ c'++d'++i', areaMult a' * areaMult b' * areaMult c' * areaMult d' * iMult2 i', map (\[_a,_b,_c] -> (Append (Ind20 $ _a-1) $ Append (Ind20 $ _b-1) $ Append (Ind20 $ _c-1) (singletonInd (Ind20 $ d-1)), Empty, singletonInd (Ind9 $ i-1), Empty, Empty, Empty) ) $ nub $ permutations [a,b,c] ) | a <- [1..21], b <- [a..21], c <- [b..21], d <- [1..21], i <- [1..10] ] --A:B:Cp:Dq areaList18_3 :: [([Int], Int, [IndTupleAbs 4 0 0 0 2 0])] areaList18_3 = list where - trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b',c',d') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trianArea d) in (a' ++ b' ++ c'++ p : d'++[q], areaMult a' * areaMult b' * areaMult c' * areaMult d', map ( \(a,b,c,p,d,q) -> (Append (Ind20 $ a-1) $ Append (Ind20 $ b-1) $ Append (Ind20 $ c-1) (singletonInd (Ind20 $ d-1)), Empty, Empty, Empty, Append (Ind3 p) (singletonInd (Ind3 q)), Empty) ) $ nub [(a,b,c,p,d,q),(b,a,c,p,d,q),(a,b,d,q,c,p),(b,a,d,q,c,p)] ) | a <- [1..21], b <- [a..21], c <- [1..21], d <- [c..21], p <- [0..3], q <- [0..3] , not (c == d && p > q) ] + list = [ let (a',b',c',d') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trianArea d) in (a' ++ b' ++ c'++ p : d'++[q], areaMult a' * areaMult b' * areaMult c' * areaMult d', map ( \(_a,_b,_c,_p,_d,_q) -> (Append (Ind20 $ _a-1) $ Append (Ind20 $ _b-1) $ Append (Ind20 $ _c-1) (singletonInd (Ind20 $ _d-1)), Empty, Empty, Empty, Append (Ind3 _p) (singletonInd (Ind3 _q)), Empty) ) $ nub [(a,b,c,p,d,q),(b,a,c,p,d,q),(a,b,d,q,c,p),(b,a,d,q,c,p)] ) | a <- [1..21], b <- [a..21], c <- [1..21], d <- [c..21], p <- [0..3], q <- [0..3] , not (c == d && p > q) ] --order 5 areaList20 :: [( [Int], Int, [IndTupleAbs 5 0 0 0 0 0])] areaList20 = list where - trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b',c', d', e') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trianArea d, (I.!) trianArea e) in (a' ++ b' ++ c' ++ d' ++ e', areaMult a' * areaMult b' * areaMult c' * areaMult d' * areaMult e', map (\[a,b,c,d,e] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ b-1) $ Append (Ind20 $ c-1) $ Append (Ind20 $ d-1) $ singletonInd (Ind20 $ e-1), Empty, Empty, Empty, Empty, Empty)) $ nub $ permutations [a,b,c,d,e] )| a <- [1..21], b <- [a..21], c <- [b..21], d <- [c..21], e <- [d..21] ] + list = [ let (a',b',c', d', e') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c, (I.!) trianArea d, (I.!) trianArea e) in (a' ++ b' ++ c' ++ d' ++ e', areaMult a' * areaMult b' * areaMult c' * areaMult d' * areaMult e', map (\[_a,_b,_c,_d,_e] -> (Append (Ind20 $ _a-1) $ Append (Ind20 $ _b-1) $ Append (Ind20 $ _c-1) $ Append (Ind20 $ _d-1) $ singletonInd (Ind20 $ _e-1), Empty, Empty, Empty, Empty, Empty)) $ nub $ permutations [a,b,c,d,e] )| a <- [1..21], b <- [a..21], c <- [b..21], d <- [c..21], e <- [d..21] ] --for the kinetic Ansätze for the Rom calculations -> extra symmetry @@ -1916,18 +1898,16 @@ areaList10Rom :: [( [Int], Int, [IndTupleAbs 2 0 0 0 2 0])] areaList10Rom = list where - trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b') = ((I.!) trianArea a, (I.!) trianArea b) in (a' ++ p : b' ++ [q], areaMult a' * areaMult b', map (\[a,p,b,q] -> (Append (Ind20 $ a-1) $ singletonInd (Ind20 $ b-1), Empty, Empty, Empty, Append (Ind3 p) $ singletonInd (Ind3 q), Empty)) $ nub [[a,p,b,q], [a,q,b,p], [b,p,a,q], [b,q,a,p]]) | a <- [1..21], b <- [a..21], p <- [0..3], q <- [p..3]] + list = [ let (a',b') = ((I.!) trianArea a, (I.!) trianArea b) in (a' ++ p : b' ++ [q], areaMult a' * areaMult b', map (\[_a,_p,_b,_q] -> (Append (Ind20 $ _a-1) $ singletonInd (Ind20 $ _b-1), Empty, Empty, Empty, Append (Ind3 _p) $ singletonInd (Ind3 _q), Empty)) $ nub [[a,p,b,q], [a,q,b,p], [b,p,a,q], [b,q,a,p]]) | a <- [1..21], b <- [a..21], p <- [0..3], q <- [p..3]] --Ap:Bq:C areaList14Rom :: [( [Int], Int, [IndTupleAbs 3 0 0 0 2 0])] areaList14Rom = list where - trian2 = trianMap2 trianArea = trianMapArea - list = [ let (a',b',c') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c) in (a' ++ p : b' ++ q : c' , areaMult a' * areaMult b' * areaMult c', map (\[[a,p],[b,q]] -> (Append (Ind20 $ a-1) $ Append (Ind20 $ b-1) $ singletonInd (Ind20 $ c-1), Empty, Empty, Empty, Append (Ind3 p) $ singletonInd (Ind3 q), Empty)) $ nub $ permutations [[a,p],[b,q]]) | a <- [1..21], b <- [a..21], c <- [1..21], p <- [0..3], q <- [0..3], not (a==b && p>q) ] + list = [ let (a',b',c') = ((I.!) trianArea a, (I.!) trianArea b, (I.!) trianArea c) in (a' ++ p : b' ++ q : c' , areaMult a' * areaMult b' * areaMult c', map (\[[_a,_p],[_b,_q]] -> (Append (Ind20 $ _a-1) $ Append (Ind20 $ _b-1) $ singletonInd (Ind20 $ c-1), Empty, Empty, Empty, Append (Ind3 _p) $ singletonInd (Ind3 _q), Empty)) $ nub $ permutations [[a,p],[b,q]]) | a <- [1..21], b <- [a..21], c <- [1..21], p <- [0..3], q <- [0..3], not (a==b && p>q) ] --now the same for the metric ansätze @@ -1955,7 +1935,7 @@ metricList4_2 = list where trianMetric = trianMap2 - list = [ let (a',b') = ((I.!) trianMetric a, (I.!) trianMetric b) in (a' ++ b', iMult2 a' * iMult2 b', map (\[a,b] -> (Empty, Empty, Append (Ind9 $ a-1) $ singletonInd (Ind9 $ b-1), Empty, Empty, Empty)) $ nub $ permutations [a,b] ) | a <- [1..10], b <- [a..10]] + list = [ let (a',b') = ((I.!) trianMetric a, (I.!) trianMetric b) in (a' ++ b', iMult2 a' * iMult2 b', map (\[_a,_b] -> (Empty, Empty, Append (Ind9 $ _a-1) $ singletonInd (Ind9 $ _b-1), Empty, Empty, Empty)) $ nub $ permutations [a,b] ) | a <- [1..10], b <- [a..10]] -- | Evaluation list for \(a^{Ap Bq}\). Note that also when using the abstract indices this ansatz still features the \( (Ap) \leftrightarrow (Bq) \) symmetry. @@ -1963,7 +1943,7 @@ metricList6_1 = list where trianMetric = trianMap2 - list = [ let (a',b') = ((I.!) trianMetric a, (I.!) trianMetric b) in (a' ++ p : b' ++ [q], iMult2 a' * iMult2 b', map (\[[a,p],[b,q]] -> (Empty, Empty, Append (Ind9 $ a-1) $ singletonInd (Ind9 $ b-1), Empty, Append (Ind3 p) $ singletonInd (Ind3 q), Empty)) $ nub $ permutations [[a,p],[b,q]]) | a <- [1..10], b <- [a..10], p <- [0..3], q <- [0..3], not (a==b && p>q)] + list = [ let (a',b') = ((I.!) trianMetric a, (I.!) trianMetric b) in (a' ++ p : b' ++ [q], iMult2 a' * iMult2 b', map (\[[_a,_p],[_b,_q]] -> (Empty, Empty, Append (Ind9 $ _a-1) $ singletonInd (Ind9 $ _b-1), Empty, Append (Ind3 _p) $ singletonInd (Ind3 _q), Empty)) $ nub $ permutations [[a,p],[b,q]]) | a <- [1..10], b <- [a..10], p <- [0..3], q <- [0..3], not (a==b && p>q)] -- | Evaluation list for \(a^{ABI} \). @@ -1979,7 +1959,7 @@ metricList6_3 = list where trianMetric = trianMap2 - list = [ let (a',b',c') = ((I.!) trianMetric a, (I.!) trianMetric b, (I.!) trianMetric c) in (a' ++ b' ++ c', iMult2 a' * iMult2 b' * iMult2 c', map (\[a,b,c] -> (Empty, Empty, Append (Ind9 $ a-1) $ Append (Ind9 $ b-1) $ singletonInd (Ind9 $ c-1), Empty, Empty, Empty)) $ nub $ permutations [a,b,c] )| a <- [1..10], b <- [a..10], c <- [b..10] ] + list = [ let (a',b',c') = ((I.!) trianMetric a, (I.!) trianMetric b, (I.!) trianMetric c) in (a' ++ b' ++ c', iMult2 a' * iMult2 b' * iMult2 c', map (\[_a,_b,_c] -> (Empty, Empty, Append (Ind9 $ _a-1) $ Append (Ind9 $ _b-1) $ singletonInd (Ind9 $ _c-1), Empty, Empty, Empty)) $ nub $ permutations [a,b,c] )| a <- [1..10], b <- [a..10], c <- [b..10] ] -- | Evaluation list for \(a^{ABp Cq}\). Note that also when using the abstract indices this ansatz still features the \( (Bp) \leftrightarrow (Cq) \) symmetry. @@ -1987,7 +1967,7 @@ metricList8_1 = list where trianMetric = trianMap2 - list = [ let (a',b',c') = ((I.!) trianMetric a, (I.!) trianMetric b, (I.!) trianMetric c) in (a' ++ b' ++ p : c' ++ [q], iMult2 a' * iMult2 b' * iMult2 c', map (\[[b,p],[c,q]] -> (Empty, Empty, Append (Ind9 $ a-1) $ Append (Ind9 $ b-1) $ singletonInd (Ind9 $ c-1), Empty, Append (Ind3 p) $ singletonInd (Ind3 q), Empty)) $ nub $ permutations [[b,p],[c,q]]) | a <- [1..10], b <- [1..10], c <- [b..10], p <- [0..3], q <- [0..3], not (b==c && p>q) ] + list = [ let (a',b',c') = ((I.!) trianMetric a, (I.!) trianMetric b, (I.!) trianMetric c) in (a' ++ b' ++ p : c' ++ [q], iMult2 a' * iMult2 b' * iMult2 c', map (\[[_b,_p],[_c,_q]] -> (Empty, Empty, Append (Ind9 $ a-1) $ Append (Ind9 $ _b-1) $ singletonInd (Ind9 $ _c-1), Empty, Append (Ind3 _p) $ singletonInd (Ind3 _q), Empty)) $ nub $ permutations [[b,p],[c,q]]) | a <- [1..10], b <- [1..10], c <- [b..10], p <- [0..3], q <- [0..3], not (b==c && p>q) ] -- | Evaluation list for \(a^{A B C I}\). Note that also when using the abstract indices this ansatz still features the \( (A) \leftrightarrow (B) \) symmetry. @@ -1995,7 +1975,7 @@ metricList8_2 = list where trianMetric = trianMap2 - list = [ let (a',b',c',i') = ((I.!) trianMetric a, (I.!) trianMetric b, (I.!) trianMetric c, (I.!) trianMetric i) in ( a' ++ b' ++ c' ++ i', iMult2 a' * iMult2 b' * iMult2 c' * iMult2 i', map (\[a,b] -> (Empty, Empty, Append (Ind9 $ a-1) $ Append (Ind9 $ b-1) $ Append (Ind9 $ c-1) $ singletonInd (Ind9 $ i-1), Empty, Empty, Empty)) $ nub $ permutations [a,b] ) | a <- [1..10], b <- [a..10], c <- [1..10], i <- [1..10] ] + list = [ let (a',b',c',i') = ((I.!) trianMetric a, (I.!) trianMetric b, (I.!) trianMetric c, (I.!) trianMetric i) in ( a' ++ b' ++ c' ++ i', iMult2 a' * iMult2 b' * iMult2 c' * iMult2 i', map (\[_a,_b] -> (Empty, Empty, Append (Ind9 $ _a-1) $ Append (Ind9 $ _b-1) $ Append (Ind9 $ c-1) $ singletonInd (Ind9 $ i-1), Empty, Empty, Empty)) $ nub $ permutations [a,b] ) | a <- [1..10], b <- [a..10], c <- [1..10], i <- [1..10] ] --symLists for the ansätze diff --git a/stack.yaml b/stack.yaml new file mode 100644 --- /dev/null +++ b/stack.yaml @@ -0,0 +1,6 @@ +resolver: lts-14.1 + +build: + haddock-arguments: + haddock-args: + - --optghc=-fdefer-type-errors diff --git a/test/Ansatz.hs b/test/Ansatz.hs new file mode 100644 --- /dev/null +++ b/test/Ansatz.hs @@ -0,0 +1,61 @@ +{-# LANGUAGE DataKinds #-} + +module Ansatz (ansatzTest) where + +import Math.Tensor +import Math.Tensor.LorentzGenerator + +import Test.Tasty +import Test.Tasty.HUnit + +ansatzForestAIEpsilon :: Bool +ansatzForestAIEpsilon = actual == expected + where + (_,eps6,_) = mkAnsatzTensorFastAbs 6 symList6 areaList6 :: (AnsatzForestEta, AnsatzForestEpsilon, ATens 1 0 1 0 0 0 AnsVarR) + actual = drawAnsatzEpsilon eps6 + expected = "(1,2,3,4)\n|\n`---- (5,6) * (16) * x[3]\n\n" + +testCase1 = testCase "AI epsilon ansatz string" $ + assertBool "epsilon ansatz forest AI does not match expected forest" ansatzForestAIEpsilon + +ansatz14_1Fast :: Bool +ansatz14_1Fast = rank == 110 + where + evalL = map (\(x,_,_) -> x) areaList14_1 + (_,_,ans) = mkAnsatzTensorFast 14 symList14_1 evalL :: (AnsatzForestEta, AnsatzForestEpsilon, STTens 14 0 AnsVarR) + rank = tensorRank2' ans + +testCase2 = testCase "ABpCq ansatz rank" $ + assertBool "rank of ABpCq ansatz does not equal 110" ansatz14_1Fast + +ansatz14_1Incr :: Bool +ansatz14_1Incr = rank == 110 + where + evalL = map (\(x,_,_) -> x) areaList14_1 + (_,_,ans) = mkAnsatzTensorIncremental 14 symList14_1 evalL :: (AnsatzForestEta, AnsatzForestEpsilon, STTens 14 0 AnsVarR) + rank = tensorRank2' ans + +testCase3 = testCase "ABpCq ansatz (incremental) rank" $ + assertBool "rank of ABpCq ansatz (incremental) does not equal 110" ansatz14_1Incr + +ansatz14_2Fast :: Bool +ansatz14_2Fast = rank == 72 + where + (_,_,ans) = mkAnsatzTensorFastAbs 14 symList14_2 areaList14_2 :: (AnsatzForestEta, AnsatzForestEpsilon, ATens 3 0 1 0 0 0 AnsVarR) + rank = tensorRank6' ans + +testCase4 = testCase "ABCI ansatz rank" $ + assertBool "rank of ABCI ansatz does not equal 72" ansatz14_2Fast + +ansatz14_2Incr :: Bool +ansatz14_2Incr = rank == 72 + where + (_,_,ans) = mkAnsatzTensorIncrementalAbs 14 symList14_2 areaList14_2 :: (AnsatzForestEta, AnsatzForestEpsilon, ATens 3 0 1 0 0 0 AnsVarR) + rank = tensorRank6' ans + +testCase5 = testCase "ABCI ansatz (incremental) rank" $ + assertBool "rank of ABCI ansatz (incremental) does not equal 72" ansatz14_2Incr + +testCase6 = testCase "" $ assertBool "" False + +ansatzTest = testGroup "AnsatzTest" [testCase1, testCase2, testCase3, testCase4, testCase5] diff --git a/test/IndList.hs b/test/IndList.hs new file mode 100644 --- /dev/null +++ b/test/IndList.hs @@ -0,0 +1,39 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE GADTs #-} + +module IndList (indListTest) where + +import Math.Tensor (IndList(..), fromList) + +import Test.Tasty +import Test.Tasty.HUnit + +il1 :: Maybe (IndList 5 Int) +il1 = fromList [1,2,3,4,5] + +il2 :: Maybe (IndList 7 Int) +il2 = fromList [1,2,3,4,5] + +il3 :: Maybe (IndList 1 Int) +il3 = fromList [1,2,3,4,5] + +testCase1 = testCase "IndListTest1" $ + assertBool + "fromList [1,2,3,4,5] :: Maybe (IndList 5 Int) failed" $ + case il1 of + Nothing -> False + Just (1 `Append` (2 `Append` (3 `Append` (4 `Append` (5 `Append` Empty))))) -> True + +testCase2 = testCase "IndListTest2" $ + assertBool "fromList [1,2,3,4,5] :: Maybe (IndList 5 Int) failed" $ + case il2 of + Nothing -> True + _ -> False + +testCase3 = testCase "IndListTest3" $ + assertBool "fromList [1,2,3,4,5] :: Maybe (IndList 5 Int) failed" $ + case il2 of + Nothing -> True + _ -> False + +indListTest = testGroup "IndListTest" [testCase1, testCase2, testCase3] diff --git a/test/LinearAlgebra.hs b/test/LinearAlgebra.hs new file mode 100644 --- /dev/null +++ b/test/LinearAlgebra.hs @@ -0,0 +1,66 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE LambdaCase #-} + +module LinearAlgebra (linearAlgebraTest) where + +import Test.QuickCheck +import Test.Tasty.QuickCheck + +import Test.Tasty + +import System.Exit + +import Numeric.LinearAlgebra (rank) +import qualified Numeric.LinearAlgebra.Data as Matrix + +import Math.Tensor.Internal.LinearAlgebra (independentColumnsMat) + +data SmallInt = S0 | S1 deriving (Show, Ord, Eq, Enum, Bounded) + +toSmall :: Int -> SmallInt +toSmall 0 = S0 +toSmall 1 = S1 +toSmall i = error $ "cannot convert " ++ show i ++ " to SmallInt" + +fromSmall :: Num a => SmallInt -> a +fromSmall S0 = 0 +fromSmall S1 = 1 + +instance Arbitrary SmallInt where + arbitrary = arbitraryBoundedEnum + +data MatrixData a = MatrixData (Positive Int) (Positive Int) [a] deriving Show + +instance Arbitrary a => Arbitrary (MatrixData a) where + arbitrary = do + m@(Positive m') <- arbitrary + n@(Positive n') <- arbitrary + xs <- vector (m'*n') + return $ MatrixData m n xs + +prop_smallValues :: MatrixData SmallInt -> Bool +prop_smallValues (MatrixData (Positive rows) (Positive cols) xs) = + rank mat' == rank mat + where + mat = (rows Matrix.>< cols) $ map fromSmall xs + mat' = independentColumnsMat mat + +prop_ints :: MatrixData Int -> Bool +prop_ints (MatrixData (Positive rows) (Positive cols) xs) = + rank mat' == rank mat + where + mat = (rows Matrix.>< cols) $ map fromIntegral xs + mat' = independentColumnsMat mat + +prop_doubles :: MatrixData Double -> Bool +prop_doubles (MatrixData (Positive rows) (Positive cols) xs) = + rank mat' == rank mat + where + mat = (rows Matrix.>< cols) xs + mat' = independentColumnsMat mat + +testCase1 = testProperty "prop_smallValues" prop_smallValues +testCase2 = testProperty "prop_ints" prop_ints +testCase3 = testProperty "prop_doubles" prop_doubles + +linearAlgebraTest = testGroup "LinearAlgebraTest" [testCase1, testCase2, testCase3] diff --git a/test/Serialization.hs b/test/Serialization.hs new file mode 100644 --- /dev/null +++ b/test/Serialization.hs @@ -0,0 +1,55 @@ +{-# LANGUAGE DataKinds #-} + +module Serialization (serializationTest) where + +import Math.Tensor +import Math.Tensor.LorentzGenerator + +import Test.Tasty +import Test.Tasty.HUnit + +ans4 :: ATens 1 0 0 0 0 0 AnsVarR +ans4 = ans4' + where + (_, _, ans4') = mkAnsatzTensorFastAbs 4 symList4 areaList4 :: (AnsatzForestEta, AnsatzForestEpsilon, ATens 1 0 0 0 0 0 AnsVarR) + +ans4Enc = encodeTensor ans4 + +ans4Success :: Bool +ans4Success = case ans4Dec of + Left _ -> False + Right ans4' -> ans4 == ans4' + where + ans4Dec = decodeTensor ans4Enc :: Either String (ATens 1 0 0 0 0 0 AnsVarR) + +testCase1 = testCase "encode and decode tensor (equal type)" $ + assertBool "decoded tensor does not match encoded tensor or failure while decoding" ans4Success + +ans4Fail1 :: Bool +ans4Fail1 = case ans4Dec of + Left _ -> True + Right _ -> False + where + ans4Dec = decodeTensor ans4Enc :: Either String (ATens 1 0 0 0 0 0 (SField Double)) + +testCase2 = testCase "encode and decode tensor (wrong type)" $ + assertBool "no error while decoding tensor of wrong type" ans4Fail1 + +ans4Fail2 :: Bool +ans4Fail2 = case ans4Dec of + Left _ -> True + Right _ -> False + where + ans4Dec = decodeTensor ans4Enc :: Either String (STTens 0 2 AnsVarR) + +testCase3 = testCase "encode and decode tensor (wrong type)" $ + assertBool "no error while decoding tensor of wrong type" ans4Fail2 + +ans4Fail3 :: Bool +ans4Fail3 = case ans4Dec of + Left _ -> True + Right _ -> False + where + ans4Dec = decodeTensor ans4Enc :: Either String (ATens 0 1 0 0 0 0 AnsVarR) + +serializationTest = testGroup "SerializationTest" [testCase1, testCase2, testCase3] diff --git a/test/TestMain.hs b/test/TestMain.hs new file mode 100644 --- /dev/null +++ b/test/TestMain.hs @@ -0,0 +1,19 @@ +module Main (main) where + +import Ansatz +import IndList +import LinearAlgebra +import Serialization + +import Test.Tasty + +test = testGroup "sparse-tensor tests" + [ + indListTest, + linearAlgebraTest, + ansatzTest, + serializationTest + ] + +main :: IO () +main = defaultMain test