diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,4 @@ +# Changelog + +## [0.1.0.0] - 2020-07-07 +Initial release. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2020 Nils Alex + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,13 @@ +# safe-tensor +Dependently typed tensor algebra in Haskell. Useful for applications in field theory, e.g., carrying out calculations for https://doi.org/10.1103/PhysRevD.101.084025 + +## Rationale +Tensor calculus is reflected in the type system. Let us view a tensor as a multilinear map from a product of vector spaces and duals thereof to the common field. The type of each tensor is its *generalized rank*, describing the vector spaces it acts on and assigning a label to each vector space. There are a few rules for tensor operations: + +- Only tensors of the same type may be added. The result is a tensor of this type. +- Tensors may be multiplied if the resulting generalized rank does not contain repeated labels for the same (dual) vector space. +- The contraction of a tensors removes pairs of vector space and dual vector space with the same label from the generalized rank. + +It is thus impossible to perform inconsistent tensor operations. + +There is also an existentially typed variant of the tensor type useful for runtime computations. These computations take place in the Error monad, throwing errors if operand types are not consistent. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/safe-tensor.cabal b/safe-tensor.cabal new file mode 100644 --- /dev/null +++ b/safe-tensor.cabal @@ -0,0 +1,59 @@ +cabal-version: 1.12 + +-- This file has been generated from package.yaml by hpack version 0.33.0. +-- +-- see: https://github.com/sol/hpack +-- +-- hash: d8db1e3bc409f22390d910f02253e624ef929b8fb88ef85dd8dce87340fc3b65 + +name: safe-tensor +version: 0.1.0.0 +synopsis: Dependently typed tensor algebra +description: Please see the README on GitHub at +category: Math +homepage: https://github.com/nilsalex/safe-tensor#readme +bug-reports: https://github.com/nilsalex/safe-tensor/issues +author: Nils Alex +maintainer: nils.alex@fau.de +copyright: 2020 Nils Alex +license: MIT +license-file: LICENSE +build-type: Simple +extra-source-files: + README.md + CHANGELOG.md + +source-repository head + type: git + location: https://github.com/nilsalex/safe-tensor + +library + exposed-modules: + Math.Tensor + Math.Tensor.Basic + Math.Tensor.Basic.Area + Math.Tensor.Basic.Delta + Math.Tensor.Basic.Epsilon + Math.Tensor.Basic.Sym2 + Math.Tensor.Basic.TH + Math.Tensor.LinearAlgebra + Math.Tensor.LinearAlgebra.Equations + Math.Tensor.LinearAlgebra.Matrix + Math.Tensor.LinearAlgebra.Scalar + Math.Tensor.Safe + Math.Tensor.Safe.Proofs + Math.Tensor.Safe.TH + Math.Tensor.Safe.Vector + other-modules: + Paths_safe_tensor + hs-source-dirs: + src + ghc-options: -Wall + build-depends: + base >=4.7 && <5 + , constraints >=0.10 && <0.13 + , containers >=0.6 && <0.7 + , hmatrix >=0.20 && <0.21 + , mtl >=2.2 && <2.3 + , singletons >=2.5 && <2.8 + default-language: Haskell2010 diff --git a/src/Math/Tensor.hs b/src/Math/Tensor.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor.hs @@ -0,0 +1,326 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE StandaloneDeriving #-} + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor +Description : Existentially quantified wrapper for Math.Tensor.Safe. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Existentially quantified wrapper around the safe interface from 'Math.Tensor.Safe'. +In contrast to the safe interface, all tensor operations are fair game, but potentially +unsafe operations take place in the Error monad 'Control.Monad.Except' and may fail +with an error message. +-} +----------------------------------------------------------------------------- +module Math.Tensor + ( -- * Existentially quantified around 'Tensor' + T(..) + , Label + -- * Fundamental operations + , rankT + , scalar + , zero + -- * Conversion from and to lists + , toListT + , fromListT + , removeZerosT + -- * Tensor algebra + , (.*) + , (.+) + , (.-) + , (#.) + -- * Tensor operations + , contractT + , transposeT + , transposeMultT + , relabelT + , -- * Constructing ranks + conRank + , covRank + , conCovRank + ) where + +import Math.Tensor.Safe +import Math.Tensor.Safe.TH + +import Data.Kind (Type) + +import Data.Singletons + ( Sing, SingI (sing), Demote + , withSomeSing, withSingI + , fromSing + ) +import Data.Singletons.Prelude + ( SBool (STrue, SFalse) + , SMaybe (SNothing, SJust) + ) +import Data.Singletons.Prelude.Maybe + ( sIsJust + ) +import Data.Singletons.Decide + ( Decision(Proved, Disproved) + , (:~:) (Refl), (%~) + ) +import Data.Singletons.TypeLits + ( Nat + , Symbol + ) + +import Data.Bifunctor (first) + +import Data.List.NonEmpty (NonEmpty((:|)),sort) + +import Control.Monad.Except (MonadError, throwError) + +-- |'T' wraps around 'Tensor' and exposes only the value type @v@. +data T :: Type -> Type where + T :: forall (r :: Rank) v. SingI r => Tensor r v -> T v + +deriving instance Show v => Show (T v) + +instance Functor T where + fmap f (T t) = T $ fmap f t + +-- |Produces a 'Scalar' of given value. Result is pure +-- because there is only one possible rank: @'[]@ +scalar :: v -> T v +scalar v = T $ Scalar v + +-- |Produces a 'ZeroTensor' of given rank @r@. Throws an +-- error if @'Sane' r ~ ''False'@. +zero :: MonadError String m => Demote Rank -> m (T v) +zero dr = + withSomeSing dr $ \sr -> + case sSane sr %~ STrue of + Proved Refl -> + case sr of + (_ :: Sing r) -> withSingI sr $ return $ T (ZeroTensor :: Tensor r v) + Disproved _ -> throwError $ "Illegal index list for zero : " ++ show dr + +vecToList :: Vec n a -> [a] +vecToList VNil = [] +vecToList (x `VCons` xs) = x : vecToList xs + +vecFromList :: forall (n :: N) a m. + MonadError String m => Sing n -> [a] -> m (Vec n a) +vecFromList SZ [] = return VNil +vecFromList (SS _) [] = throwError "List provided for vector reconstruction is too short." +vecFromList SZ (_:_) = throwError "List provided for vector reconstruction is too long." +vecFromList (SS sn) (x:xs) = do + xs' <- vecFromList sn xs + return $ x `VCons` xs' + +-- |Pure function removing all zeros from a tensor. Wraps around 'removeZeros'. +removeZerosT :: (Eq v, Num v) => T v -> T v +removeZerosT o = + case o of + T t -> T $ removeZeros t + +-- |Tensor product. Returns an error if ranks overlap, i.e. +-- @MergeR r1 r2 ~ 'Nothing@. Wraps around '(&*)'. +(.*) :: (Num v, MonadError String m) => T v -> T v -> m (T v) +(.*) o1 o2 = + case o1 of + T (t1 :: Tensor r1 v) -> + case o2 of + T (t2 :: Tensor r2 v) -> + let sr1 = sing :: Sing r1 + sr2 = sing :: Sing r2 + in case sMergeR sr1 sr2 of + SNothing -> throwError "Tensors have overlapping indices. Cannot multiply." + SJust sr' -> withSingI sr' $ return $ T (t1 &* t2) +infixl 7 .* + +-- |Scalar multiplication of a tensor. +(#.) :: Num v => v -> T v -> T v +(#.) s = fmap (*s) + +infixl 7 #. + +-- |Tensor addition. Returns an error if summand ranks do not coincide. Wraps around '(&+)'. +(.+) :: (Eq v, Num v, MonadError String m) => T v -> T v -> m (T v) +(.+) o1 o2 = + case o1 of + T (t1 :: Tensor r1 v) -> + case o2 of + T (t2 :: Tensor r2 v) -> + let sr1 = sing :: Sing r1 + sr2 = sing :: Sing r2 + in case sr1 %~ sr2 of + Proved Refl -> case sSane sr1 %~ STrue of + Proved Refl -> return $ T (t1 &+ t2) + Disproved _ -> throwError "Rank of summands is not sane." + Disproved _ -> throwError "Generalized tensor ranks do not match. Cannot add." +infixl 6 .+ + +-- |Tensor subtraction. Returns an error if summand ranks do not coincide. Wraps around '(&-)'. +(.-) :: (Eq v, Num v, MonadError String m) => T v -> T v -> m (T v) +(.-) o1 o2 = + case o1 of + T (t1 :: Tensor r1 v) -> + case o2 of + T (t2 :: Tensor r2 v) -> + let sr1 = sing :: Sing r1 + sr2 = sing :: Sing r2 + in case sr1 %~ sr2 of + Proved Refl -> case sSane sr1 %~ STrue of + Proved Refl -> return $ T (t1 &- t2) + Disproved _ -> throwError "Rank of operands is not sane." + Disproved _ -> throwError "Generalized tensor ranks do not match. Cannot add." + +-- |Tensor contraction. Pure function, because a tensor of any rank can be contracted. +-- Wraps around 'contract'. +contractT :: (Num v, Eq v) => T v -> T v +contractT o = + case o of + T (t :: Tensor r v) -> + let sr = sing :: Sing r + sr' = sContractR sr + in withSingI sr' $ T $ contract t + +-- |Tensor transposition. Returns an error if given indices cannot be transposed. +-- Wraps around 'transpose'. +transposeT :: MonadError String m => + Demote (VSpace Symbol Nat) -> Demote (Ix Symbol) -> Demote (Ix Symbol) -> + T v -> m (T v) +transposeT v ia ib o = + case o of + T (t :: Tensor r v) -> + let sr = sing :: Sing r + in withSingI sr $ + withSomeSing v $ \sv -> + withSomeSing ia $ \sia -> + withSomeSing ib $ \sib -> + case sCanTranspose sv sia sib sr of + STrue -> return $ T $ transpose sv sia sib t + SFalse -> throwError $ "Cannot transpose indices " ++ show v ++ " " ++ show ia ++ " " ++ show ib ++ "!" + +-- |Transposition of multiple indices. Returns an error if given indices cannot be transposed. +-- Wraps around 'transposeMult'. +transposeMultT :: MonadError String m => + Demote (VSpace Symbol Nat) -> Demote [(Symbol,Symbol)] -> Demote [(Symbol,Symbol)] -> T v -> m (T v) +transposeMultT _ [] [] _ = throwError "Empty lists for transpositions!" +transposeMultT v (con:cons) [] o = + case o of + T (t :: Tensor r v) -> + let sr = sing :: Sing r + cons' = sort $ con :| cons + tr = TransCon (fmap fst cons') (fmap snd cons') + in withSingI sr $ + withSomeSing v $ \sv -> + withSomeSing tr $ \str -> + case sIsJust (sTranspositions sv str sr) %~ STrue of + Proved Refl -> return $ T $ transposeMult sv str t + Disproved _ -> throwError $ "Cannot transpose indices " ++ show v ++ " " ++ show tr ++ "!" +transposeMultT v [] (cov:covs) o = + case o of + T (t :: Tensor r v) -> + let sr = sing :: Sing r + covs' = sort $ cov :| covs + tr = TransCov (fmap fst covs') (fmap snd covs') + in withSingI sr $ + withSomeSing v $ \sv -> + withSomeSing tr $ \str -> + case sIsJust (sTranspositions sv str sr) %~ STrue of + Proved Refl -> return $ T $ transposeMult sv str t + Disproved _ -> throwError $ "Cannot transpose indices " ++ show v ++ " " ++ show tr ++ "!" +transposeMultT _ _ _ _ = throwError "Simultaneous transposition of contravariant and covariant indices not yet supported!" + +-- |Relabelling of tensor indices. Returns an error if given relabellings are not allowed. +-- Wraps around 'relabel'. +relabelT :: MonadError String m => + Demote (VSpace Symbol Nat) -> Demote [(Symbol,Symbol)] -> T v -> m (T v) +relabelT _ [] _ = throwError "Empty list for relabelling!" +relabelT v (r:rs) o = + case o of + T (t :: Tensor r v) -> + let sr = sing :: Sing r + rr = sort $ r :| rs + in withSingI sr $ + withSomeSing v $ \sv -> + withSomeSing rr $ \srr -> + case sRelabelR sv srr sr of + SJust sr' -> + withSingI sr' $ + case sSane sr' %~ STrue of + Proved Refl -> return $ T $ relabel sv srr t + Disproved _ -> throwError $ "Cannot relabel indices " ++ show v ++ " " ++ show rr ++ "!" + _ -> throwError $ "Cannot relabel indices " ++ show v ++ " " ++ show rr ++ "!" + +-- |Exposes the hidden rank over which 'T' quantifies. Possible because of the @'SingI' r@ constraint. +rankT :: T v -> Demote Rank +rankT o = + case o of + T (_ :: Tensor r v) -> + let sr = sing :: Sing r + in fromSing sr + +-- |Assocs list of the tensor. +toListT :: T v -> [([Int], v)] +toListT o = + case o of + T (t :: Tensor r v) -> let sr = sing :: Sing r + sn = sLengthR sr + in withSingI sn $ + first vecToList <$> toList t + +-- |Constructs a tensor from a rank and an assocs list. Returns an error for illegal ranks +-- or incompatible assocs lists. +fromListT :: MonadError String m => Demote Rank -> [([Int], v)] -> m (T v) +fromListT r xs = + withSomeSing r $ \sr -> + withSingI sr $ + let sn = sLengthR sr + in case sSane sr %~ STrue of + Proved Refl -> T . fromList' sr <$> + mapM (\(vec, val) -> do + vec' <- vecFromList sn vec + return (vec', val)) xs + Disproved _ -> throwError $ "Insane tensor rank : " <> show r + +-- |The unrefined type of labels. +-- +-- @ Demote Symbol ~ Text @ +type Label = Demote Symbol + +-- |Lifts sanity check of ranks into the error monad. +saneRank :: (Ord s, Ord n, MonadError String m) => GRank s n -> m (GRank s n) +saneRank r + | sane r = pure r + | otherwise = throwError "Index lists must be strictly ascending." + +-- |Creates contravariant rank from vector space labe, vector space dimension, +-- and list of index labels. Returns an error for illegal ranks. +conRank :: (MonadError String m, Integral a, Ord s, Ord n, Num n) => + s -> a -> [s] -> m (GRank s n) +conRank _ _ [] = throwError "Generalized rank must have non-vanishing index list!" +conRank v d (i:is) = saneRank [(VSpace v (fromIntegral d), Con (i :| is))] + +-- |Creates covariant rank from vector space labe, vector space dimension, +-- and list of index labels. Returns an error for illegal ranks. +covRank :: (MonadError String m, Integral a, Ord s, Ord n, Num n) => + s -> a -> [s] -> m (GRank s n) +covRank _ _ [] = throwError "Generalized rank must have non-vanishing index list!" +covRank v d (i:is) = saneRank [(VSpace v (fromIntegral d), Cov (i :| is))] + +-- |Creates mixed rank from vector space label, vector space dimension, +-- and lists of index labels. Returns an error for illegal ranks. +conCovRank :: (MonadError String m, Integral a, Ord s, Ord n, Num n) => + s -> a -> [s] -> [s] -> m (GRank s n) +conCovRank _ _ _ [] = throwError "Generalized rank must have non-vanishing index list!" +conCovRank _ _ [] _ = throwError "Generalized rank must have non-vanishing index list!" +conCovRank v d (i:is) (j:js) = saneRank [(VSpace v (fromIntegral d), ConCov (i :| is) (j :| js))] diff --git a/src/Math/Tensor/Basic.hs b/src/Math/Tensor/Basic.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/Basic.hs @@ -0,0 +1,28 @@ +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.Basic +Description : Definitions of basic tensors. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Definitions of basic tensors, re-exported for convenience. +-} +----------------------------------------------------------------------------- + +module Math.Tensor.Basic + ( -- * Kronecker delta + module Math.Tensor.Basic.Delta + , -- * Epsilon tensor densities + module Math.Tensor.Basic.Epsilon + , -- * Symmetric tensors + module Math.Tensor.Basic.Sym2 + , -- * Area-symmetric tensors + module Math.Tensor.Basic.Area + ) where + +import Math.Tensor.Basic.Delta +import Math.Tensor.Basic.Epsilon +import Math.Tensor.Basic.Sym2 +import Math.Tensor.Basic.Area diff --git a/src/Math/Tensor/Basic/Area.hs b/src/Math/Tensor/Basic/Area.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/Basic/Area.hs @@ -0,0 +1,366 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TypeApplications #-} + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.Basic.Area +Description : Definitions of area-symmetric tensors. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Definitions of area-symmetric tensors. +-} +----------------------------------------------------------------------------- +module Math.Tensor.Basic.Area + ( -- * Area metric tensor induced from flat Lorentzian metric + flatAreaCon + , someFlatAreaCon + , -- * Injections from \(AS(V)\) into \(V\times V\times V\times V\) + injAreaCon' + , injAreaCov' + , someInjAreaCon + , someInjAreaCov + , -- * Surjections from \(V\times V\times V\times V\) onto \(AS(V)\) + surjAreaCon' + , surjAreaCov' + , someSurjAreaCon + , someSurjAreaCov + , -- * Vertical coefficients for functions on \(AS(V)\) + someInterAreaCon + , someInterAreaCov + , -- * Kronecker delta on \(AS(V)\) + someDeltaArea + , -- * Internals + trianMapArea + , facMapArea + , areaSign + , sortArea + ) where + +import Math.Tensor +import Math.Tensor.Safe +import Math.Tensor.Safe.TH +import Math.Tensor.Basic.TH +import Math.Tensor.Basic.Delta + +import Data.Singletons + ( Sing + , SingI (sing) + , Demote + , withSomeSing + , withSingI + ) +import Data.Singletons.Prelude + ( SBool (STrue) + , PSemigroup ((<>)) + , SMaybe (SJust) + , (%<>) + ) +import Data.Singletons.Decide + ( (:~:) (Refl) + , Decision (Proved) + , (%~) + ) +import Data.Singletons.TypeLits + ( Symbol + , withKnownSymbol + ) + +import Data.Maybe (catMaybes) +import Data.Ratio (Ratio, numerator, denominator) +import qualified Data.Map.Strict as Map + ( Map + , lookup + , fromList + ) +import Data.List.NonEmpty (NonEmpty((:|))) +import Control.Monad.Except + ( MonadError + , runExcept + ) + +trianMapArea :: Map.Map (Vec ('S ('S ('S ('S 'Z)))) Int) Int +trianMapArea = Map.fromList $ zip indices4 indices1 + where + indices1 :: [Int] + indices1 = [0..] + indices4 :: [Vec ('S ('S ('S ('S 'Z)))) Int] + indices4 = [a `VCons` (b `VCons` (c `VCons` (d `VCons` VNil))) | + a <- [0..2], + b <- [a+1..3], + c <- [a..2], + d <- [c+1..3], + not (a == c && b > d) ] + +facMapArea :: forall b.Num b => Map.Map (Vec ('S ('S ('S ('S 'Z)))) Int) b +facMapArea = Map.fromList $ [(a `VCons` (b `VCons` (c `VCons` (d `VCons` VNil))), fac a b c d) | + a <- [0..2], + b <- [a+1..3], + c <- [a..2], + d <- [c+1..3], + not (a == c && b > d)] + where + fac :: Int -> Int -> Int -> Int -> b + fac a b c d + | a == c && b == d = 4 + | otherwise = 8 + +areaSign :: (Ord a, Num v) => a -> a -> a -> a -> Maybe v +areaSign a b c d + | a == b = Nothing + | c == d = Nothing + | a > b = ((-1) *) <$> areaSign b a c d + | c > d = ((-1) *) <$> areaSign a b d c + | otherwise = Just 1 + +sortArea :: Ord a => a -> a -> a -> a -> Vec ('S ('S ('S ('S 'Z)))) a +sortArea a b c d + | a > b = sortArea b a c d + | c > d = sortArea a b d c + | (a,b) > (c,d) = sortArea c d a b + | otherwise = a `VCons` (b `VCons` (c `VCons` (d `VCons` VNil))) + + +injAreaCon' :: forall (id :: Symbol) (a :: Symbol) (b :: Symbol) ( c :: Symbol) (d :: Symbol) + (i :: Symbol) (r :: Rank) v. + ( + InjAreaConRank id a b c d i ~ 'Just r, + SingI r, + Num v + ) => Sing id -> Sing a -> Sing b -> Sing c -> Sing d -> Sing i -> Tensor r v +injAreaCon' _ _ _ _ _ _ = + case sLengthR sr of + SS (SS (SS (SS (SS SZ)))) -> + case sSane sr %~ STrue of + Proved Refl -> fromList assocs + where + sr = sing :: Sing r + tm = trianMapArea + assocs = catMaybes $ + (\a b c d -> + do + s <- areaSign a b c d + let v = sortArea a b c d + i <- Map.lookup v tm + return (a `VCons` (b `VCons` (c `VCons` (d `VCons` (i `VCons` VNil)))), s :: v)) + <$> [0..3] <*> [0..3] <*> [0..3] <*> [0..3] + +injAreaCov' :: forall (id :: Symbol) (a :: Symbol) (b :: Symbol) ( c :: Symbol) (d :: Symbol) + (i :: Symbol) (r :: Rank) v. + ( + InjAreaCovRank id a b c d i ~ 'Just r, + SingI r, + Num v + ) => Sing id -> Sing a -> Sing b -> Sing c -> Sing d -> Sing i -> Tensor r v +injAreaCov' _ _ _ _ _ _ = + case sLengthR sr of + SS (SS (SS (SS (SS SZ)))) -> + case sSane sr %~ STrue of + Proved Refl -> fromList assocs + where + sr = sing :: Sing r + tm = trianMapArea + assocs = catMaybes $ + (\a b c d -> + do + s <- areaSign a b c d + let v = sortArea a b c d + i <- Map.lookup v tm + return (a `VCons` (b `VCons` (c `VCons` (d `VCons` (i `VCons` VNil)))), s :: v)) + <$> [0..3] <*> [0..3] <*> [0..3] <*> [0..3] + +surjAreaCon' :: forall (id :: Symbol) (a :: Symbol) (b :: Symbol) ( c :: Symbol) (d :: Symbol) + (i :: Symbol) (r :: Rank) v. + ( + SurjAreaConRank id a b c d i ~ 'Just r, + SingI r, + Fractional v + ) => Sing id -> Sing a -> Sing b -> Sing c -> Sing d -> Sing i -> Tensor r v +surjAreaCon' _ _ _ _ _ _ = + case sLengthR sr of + SS (SS (SS (SS (SS SZ)))) -> + case sSane sr %~ STrue of + Proved Refl -> fromList assocs + where + sr = sing :: Sing r + tm = trianMapArea + fm = facMapArea @v + assocs = catMaybes $ + (\a b c d -> + do + s <- areaSign a b c d + let v = sortArea a b c d + i <- Map.lookup v tm + f <- Map.lookup v fm + return (a `VCons` (b `VCons` (c `VCons` (d `VCons` (i `VCons` VNil)))), s/f :: v)) + <$> [0..3] <*> [0..3] <*> [0..3] <*> [0..3] + +surjAreaCov' :: forall (id :: Symbol) (a :: Symbol) (b :: Symbol) ( c :: Symbol) (d :: Symbol) + (i :: Symbol) (r :: Rank) v. + ( + SurjAreaCovRank id a b c d i ~ 'Just r, + SingI r, + Fractional v + ) => Sing id -> Sing a -> Sing b -> Sing c -> Sing d -> Sing i -> Tensor r v +surjAreaCov' _ _ _ _ _ _ = + case sLengthR sr of + SS (SS (SS (SS (SS SZ)))) -> + case sSane sr %~ STrue of + Proved Refl -> fromList assocs + where + sr = sing :: Sing r + tm = trianMapArea + fm = facMapArea @v + assocs = catMaybes $ + (\a b c d -> + do + s <- areaSign a b c d + let v = sortArea a b c d + i <- Map.lookup v tm + f <- Map.lookup v fm + return (a `VCons` (b `VCons` (c `VCons` (d `VCons` (i `VCons` VNil)))), s/f :: v)) + <$> [0..3] <*> [0..3] <*> [0..3] <*> [0..3] + +_injAreaCon :: Num v => Demote Symbol -> Demote Symbol -> T v +_injAreaCon vid i = + withSomeSing vid $ \sid -> + withSomeSing i $ \si -> + withSomeSing " 01" $ \s01 -> + withSomeSing " 02" $ \s02 -> + withSomeSing " 03" $ \s03 -> + withSomeSing " 04" $ \s04 -> + case sInjAreaConRank sid s01 s02 s03 s04 si of + SJust sr -> withSingI sr $ T $ injAreaCon' sid s01 s02 s03 s04 si + +_injAreaCov :: Num v => Demote Symbol -> Demote Symbol -> T v +_injAreaCov vid i = + withSomeSing vid $ \sid -> + withSomeSing i $ \si -> + withSomeSing " 01" $ \s01 -> + withSomeSing " 02" $ \s02 -> + withSomeSing " 03" $ \s03 -> + withSomeSing " 04" $ \s04 -> + case sInjAreaCovRank sid s01 s02 s03 s04 si of + SJust sr -> withSingI sr $ T $ injAreaCov' sid s01 s02 s03 s04 si + +_surjAreaCon :: Fractional v => Demote Symbol -> Demote Symbol -> T v +_surjAreaCon vid i = + withSomeSing vid $ \sid -> + withSomeSing i $ \si -> + withSomeSing " 01" $ \s01 -> + withSomeSing " 02" $ \s02 -> + withSomeSing " 03" $ \s03 -> + withSomeSing " 04" $ \s04 -> + case sSurjAreaConRank sid s01 s02 s03 s04 si of + SJust sr -> withSingI sr $ T $ surjAreaCon' sid s01 s02 s03 s04 si + +_surjAreaCov :: Fractional v => Demote Symbol -> Demote Symbol -> T v +_surjAreaCov vid i = + withSomeSing vid $ \sid -> + withSomeSing i $ \si -> + withSomeSing " 01" $ \s01 -> + withSomeSing " 02" $ \s02 -> + withSomeSing " 03" $ \s03 -> + withSomeSing " 04" $ \s04 -> + case sSurjAreaCovRank sid s01 s02 s03 s04 si of + SJust sr -> withSingI sr $ T $ surjAreaCov' sid s01 s02 s03 s04 si + +someInjAreaCon :: forall v m.(Num v, MonadError String m) => + Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + m (T v) +someInjAreaCon vid a b c d i = relabelT (VSpace vid 4) [(" 01",a), (" 02",b), (" 03",c), (" 04",d)] t + where + t = _injAreaCon vid i :: T v + +someInjAreaCov :: forall v m.(Num v, MonadError String m) => + Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + m (T v) +someInjAreaCov vid a b c d i = relabelT (VSpace vid 4) [(" 01",a), (" 02",b), (" 03",c), (" 04",d)] t + where + t = _injAreaCov vid i :: T v + +someSurjAreaCon :: forall v m.(Fractional v, MonadError String m) => + Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + m (T v) +someSurjAreaCon vid a b c d i = relabelT (VSpace vid 4) [(" 01",a), (" 02",b), (" 03",c), (" 04",d)] t + where + t = _surjAreaCon vid i :: T v + +someSurjAreaCov :: forall v m.(Fractional v, MonadError String m) => + Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + m (T v) +someSurjAreaCov vid a b c d i = relabelT (VSpace vid 4) [(" 01",a), (" 02",b), (" 03",c), (" 04",d)] t + where + t = _surjAreaCov vid i :: T v + +someInterAreaCon :: Num v => + Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + T v +someInterAreaCon vid m n a b = t + where + Right t = runExcept $ + do + j :: T (Ratio Int) <- someSurjAreaCon vid " 1" " 2" " 3" n a + i :: T (Ratio Int) <- someInjAreaCon vid " 1" " 2" " 3" m b + prod <- i .* j + let res = contractT $ fmap (*(-4)) prod + return $ fmap (\v -> if denominator v == 1 + then fromIntegral (numerator v) + else error "someInterAreaCon is not fraction-free, as it should be!") res + +someInterAreaCov :: Num v => + Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + T v +someInterAreaCov vid m n a b = t + where + Right t = runExcept $ + do + j :: T (Ratio Int) <- someSurjAreaCov vid " 1" " 2" " 3" m a + i :: T (Ratio Int) <- someInjAreaCov vid " 1" " 2" " 3" n b + prod <- i .* j + let res = contractT $ fmap (*4) prod + return $ fmap (\v -> if denominator v == 1 + then fromIntegral (numerator v) + else error "someInterAreaCov is not fraction-free, as it should be!") res + +someDeltaArea :: Num v => Demote Symbol -> Demote Symbol -> Demote Symbol -> T v +someDeltaArea vid = someDelta (vid <> "Area") 21 + +flatAreaCon :: forall (id :: Symbol) (a :: Symbol) (r :: Rank) v. + ( + '[ '( 'VSpace (id <> "Area") 21, 'Con (a ':| '[]))] ~ r, + Num v + ) => Sing id -> Sing a -> Tensor r v +flatAreaCon sid sa = + withKnownSymbol (sid %<> (sing :: Sing "Area")) $ + withKnownSymbol sa $ + fromList [(0 `VCons` VNil, -1), (5 `VCons` VNil, 1), -- 0 1 2 3 4 5 + (6 `VCons` VNil, -1), (9 `VCons` VNil, -1), -- 6 7 8 9 10 + (11 `VCons` VNil, -1), (12 `VCons` VNil, 1), -- 11 12 13 14 + (15 `VCons` VNil, 1), -- 15 16 17 + (18 `VCons` VNil, 1), -- 18 19 + (20 `VCons` VNil, 1)] -- 20 + + + +someFlatAreaCon :: Num v => Demote Symbol -> Demote Symbol -> T v +someFlatAreaCon vid a = + withSomeSing vid $ \sid -> + withSomeSing a $ \sa -> + withKnownSymbol (sid %<> (sing :: Sing "Area")) $ + withKnownSymbol sa $ + T $ flatAreaCon sid sa diff --git a/src/Math/Tensor/Basic/Delta.hs b/src/Math/Tensor/Basic/Delta.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/Basic/Delta.hs @@ -0,0 +1,118 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.Basic.Delta +Description : Definitions of Kronecker deltas. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Definitions of Kronecker deltas \(\delta^{a}_{\hphantom ab}\) +(identity automorphisms) for arbitrary vector spaces. +-} +----------------------------------------------------------------------------- +module Math.Tensor.Basic.Delta + ( -- * Kronecker delta + delta + , delta' + , someDelta + ) where + +import Math.Tensor +import Math.Tensor.Safe +import Math.Tensor.Safe.TH +import Math.Tensor.Basic.TH + +import Data.Singletons + ( Sing + , SingI (sing) + , Demote + , withSomeSing + ) +import Data.Singletons.Prelude + ( SList (SNil) + , SBool (STrue) + ) +import Data.Singletons.Decide + ( (:~:) (Refl) + , Decision (Proved) + , (%~) + ) +import Data.Singletons.TypeLits + ( Symbol + , Nat + , KnownNat + , withKnownNat + , natVal + , withKnownSymbol + ) + +import Data.List.NonEmpty (NonEmpty ((:|))) + +-- |The Kronecker delta \(\delta^a_{\hphantom ab} \) for a given +-- @'VSpace' id n@ with contravariant +-- index label @a@ and covariant index label @b@. Labels and dimension +-- are passed explicitly as singletons. +delta' :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) (r :: Rank) v. + ( + KnownNat n, + Num v, + '[ '( 'VSpace id n, 'ConCov (a ':| '[]) (b ':| '[])) ] ~ r, + Tail' (Tail' r) ~ '[], + Sane (Tail' r) ~ 'True + ) => + Sing id -> Sing n -> Sing a -> Sing b -> + Tensor r v +delta' _ _ _ _ = delta + +-- |The Kronecker delta \(\delta^a_{\hphantom ab} \) for a given +-- @'VSpace' id n@ with contravariant +-- index label @a@ and covariant index label @b@. +delta :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) (r :: Rank) v. + ( + '[ '( 'VSpace id n, 'ConCov (a ':| '[]) (b ':| '[]))] ~ r, + Tail' (Tail' r) ~ '[], + Sane (Tail' r) ~ 'True, + SingI n, + Num v + ) => Tensor r v +delta = case (sing :: Sing n) of + sn -> let x = fromIntegral $ withKnownNat sn $ natVal sn + in Tensor (f x) + where + f :: Int -> [(Int, Tensor (Tail' r) v)] + f x = map (\i -> (i, Tensor [(i, Scalar 1)])) [0..x - 1] + +-- |The Kronecker delta \(\delta^a_{\hphantom ab} \) for a given +-- @'VSpace' id n@ with contravariant +-- index label @a@ and covariant index label @b@. Labels and dimension +-- are passed as values. Result is existentially quantified. +someDelta :: Num v => + Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> + T v +someDelta vid vdim a b = + withSomeSing vid $ \svid -> + withSomeSing vdim $ \svdim -> + withSomeSing a $ \sa -> + withSomeSing b $ \sb -> + withKnownNat svdim $ + withKnownSymbol svid $ + withKnownSymbol sa $ + withKnownSymbol sb $ + let sl = sDeltaRank svid svdim sa sb + in case sTail' (sTail' sl) of + SNil -> case sSane (sTail' sl) %~ STrue of + Proved Refl -> T $ delta' svid svdim sa sb diff --git a/src/Math/Tensor/Basic/Epsilon.hs b/src/Math/Tensor/Basic/Epsilon.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/Basic/Epsilon.hs @@ -0,0 +1,184 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE OverloadedStrings #-} + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.Basic.Epsilon +Description : Definitions of epsilon tensor densities. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Definitions of covariant and contravariant epsilon tensor densities +like \(\epsilon_{abc}\). +-} +----------------------------------------------------------------------------- +module Math.Tensor.Basic.Epsilon + ( -- * Epsilon tensor densities + epsilon' + , someEpsilon + , epsilonInv' + , someEpsilonInv + , -- * Internals + permSign + ) where + +import Math.Tensor +import Math.Tensor.Safe +import Math.Tensor.Safe.TH +import Math.Tensor.Basic.TH + +import Data.Singletons + ( Sing + , SingI (sing) + , Demote + , withSomeSing + , withSingI + , fromSing + ) +import Data.Singletons.Prelude + ( SBool (STrue) + , SMaybe (SNothing, SJust) + ) +import Data.Singletons.Decide + ( (:~:) (Refl) + , Decision (Proved) + , (%~) + ) +import Data.Singletons.TypeLits + ( Symbol + , Nat + , KnownNat + , withKnownNat + ) + +import Data.List (sort,permutations) +import qualified Data.List.NonEmpty as NE (NonEmpty((:|)),sort) + +import Control.Monad.Except (MonadError, throwError) + +-- |Sign of a permutation: +-- +-- @ +-- permSign [1,2,3] = 1 +-- permSign [2,1,3] = -1 +-- @ +permSign :: (Num v, Ord a) => [a] -> v +permSign (x:xs) + | even (length defects) = permSign xs + | odd (length defects) = (-1) * permSign xs + where + defects = filter ( + Sing id -> Sing n -> Sing is -> Tensor r v +epsilon' _ sn _ = + case sLengthR sr %~ sn' of + Proved Refl -> + case sSane sr %~ STrue of + Proved Refl -> fromList xs + where + sr = sing :: Sing r + sn' = sFromNat sn + n = fromSing sn + perms = sort $ permutations $ take (fromIntegral n) [(0 :: Int)..] + xs = fmap (\p -> (vecFromListUnsafe sn' p, permSign p :: v)) perms + +-- |Totally antisymmetric contravariant tensor density of weight +1 +-- such that +-- +-- \[ \epsilon^{12\dots n} = 1. \] +-- +-- Vector space label, vector space dimension and index labels +-- are passed as singletons. +epsilonInv' :: forall (id :: Symbol) (n :: Nat) (is :: NE.NonEmpty Symbol) (r :: Rank) v. + ( + KnownNat n, + Num v, + EpsilonInvRank id n is ~ 'Just r, + SingI r + ) => + Sing id -> Sing n -> Sing is -> Tensor r v +epsilonInv' _ sn _ = + case sLengthR sr %~ sn' of + Proved Refl -> + case sSane sr %~ STrue of + Proved Refl -> fromList xs + where + sr = sing :: Sing r + sn' = sFromNat sn + n = fromSing sn + perms = sort $ permutations $ take (fromIntegral n) [(0 :: Int)..] + xs = fmap (\p -> (vecFromListUnsafe sn' p, permSign p :: v)) perms + +-- |Totally antisymmetric covariant tensor density of weight -1 +-- such that +-- +-- \[ \epsilon_{12\dots n} = 1. \] +-- +-- Vector space label, vector space dimension and index labels +-- are passed as values. Result is existentially quantified. +someEpsilon :: forall v m.(Num v, MonadError String m) => + Demote Symbol -> Demote Nat -> [Demote Symbol] -> + m (T v) +someEpsilon _ _ [] = throwError "Empty index list!" +someEpsilon vid vdim (i:is) = + let sign = permSign (i:is) :: v + in withSomeSing vid $ \svid -> + withSomeSing vdim $ \svdim -> + withSomeSing (NE.sort (i NE.:| is)) $ \sis -> + withKnownNat svdim $ + case sEpsilonRank svid svdim sis of + SJust sr -> + withSingI sr $ + return $ T $ (* sign) <$> epsilon' svid svdim sis + SNothing -> throwError $ "Illegal index list " ++ show (i:is) ++ "!" + +-- |Totally antisymmetric contravariant tensor density of weight +1 +-- such that +-- +-- \[ \epsilon^{12\dots n} = 1. \] +-- +-- Vector space label, vector space dimension and index labels +-- are passed as values. Result is existentially quantified. +someEpsilonInv :: forall v m.(Num v, MonadError String m) => + Demote Symbol -> Demote Nat -> [Demote Symbol] -> + m (T v) +someEpsilonInv _ _ [] = throwError "Empty index list!" +someEpsilonInv vid vdim (i:is) = + let sign = permSign (i:is) :: v + in withSomeSing vid $ \svid -> + withSomeSing vdim $ \svdim -> + withSomeSing (NE.sort (i NE.:| is)) $ \sis -> + withKnownNat svdim $ + case sEpsilonInvRank svid svdim sis of + SJust sr -> + withSingI sr $ + return $ T $ (* sign) <$> epsilonInv' svid svdim sis + SNothing -> throwError $ "Illegal index list " ++ show (i:is) ++ "!" diff --git a/src/Math/Tensor/Basic/Sym2.hs b/src/Math/Tensor/Basic/Sym2.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/Basic/Sym2.hs @@ -0,0 +1,501 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE OverloadedStrings #-} + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.Basic.Sym2 +Description : Definitions of symmetric tensors. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Definitions of symmetric tensors. +-} +----------------------------------------------------------------------------- +module Math.Tensor.Basic.Sym2 + ( -- * Flat positive-definite metric + gamma + , gamma' + , someGamma + , gammaInv + , gammaInv' + , someGammaInv + , -- * Flat Lorentzian metric + eta + , eta' + , someEta + , etaInv + , etaInv' + , someEtaInv + , -- * Injections from \(S^2V\) into \(V\times V\) + injSym2Con' + , injSym2Cov' + , someInjSym2Con + , someInjSym2Cov + , -- * Surjections from \(V\times V\) onto \(S^2V\) + surjSym2Con' + , surjSym2Cov' + , someSurjSym2Con + , someSurjSym2Cov + , -- * Vertical coefficients for functions on \(S^2V\) + someInterSym2Con + , someInterSym2Cov + , -- * Kronecker delta on \(S^2V\) + someDeltaSym2 + , -- * Internals + trianMapSym2 + , facMapSym2 + , sym2Assocs + , sym2AssocsFac + ) where + +import Math.Tensor +import Math.Tensor.Safe +import Math.Tensor.Safe.TH +import Math.Tensor.Basic.TH +import Math.Tensor.Basic.Delta + +import Data.Singletons + ( Sing + , SingI (sing) + , Demote + , withSomeSing + , withSingI + ) +import Data.Singletons.Prelude + ( SBool (STrue) + , POrd ((<)) + , SOrdering (SLT) + , SMaybe (SJust) + , sCompare + ) +import Data.Singletons.Decide + ( (:~:) (Refl) + , Decision (Proved) + , (%~) + ) +import Data.Singletons.TypeLits + ( Symbol + , Nat + , withKnownNat + , natVal + , withKnownSymbol + ) + +import Data.Ratio (Ratio, numerator, denominator) +import Data.List.NonEmpty (NonEmpty((:|))) +import qualified Data.Map.Strict as Map + ( Map + , fromList + , lookup + ) + +import Control.Monad.Except + ( MonadError + , throwError + , runExcept + ) + +trianMapSym2 :: forall a.Integral a => a -> Map.Map (Vec ('S ('S 'Z)) Int) Int +trianMapSym2 n = Map.fromList $ zip indices2 indices1 + where + maxInd :: Int + maxInd = fromIntegral n - 1 + indices1 :: [Int] + indices1 = [0..] + indices2 = [a `VCons` (b `VCons` VNil) | a <- [0..maxInd], b <- [a..maxInd] ] + +facMapSym2 :: forall a b.(Integral a, Num b) => a -> Map.Map (Vec ('S ('S 'Z)) Int) b +facMapSym2 n = Map.fromList $ [(a `VCons` (b `VCons` VNil), fac a b) | + a <- [0..maxInd], b <- [a..maxInd] ] + where + maxInd = fromIntegral n - 1 + fac :: Int -> Int -> b + fac a b + | a == b = 1 + | otherwise = 2 + +sym2Assocs :: forall (n :: Nat) v.Num v => Sing n -> [(Vec ('S ('S ('S 'Z))) Int, v)] +sym2Assocs sn = assocs + where + n = withKnownNat sn (natVal sn) + tm = trianMapSym2 n + maxInd = fromIntegral n - (1 :: Int) + assocs = (\a b -> let v = vec a b + in case Map.lookup v tm of + Just i -> (a `VCons` (b `VCons` (i `VCons` VNil)), 1 :: v) + _ -> error $ "indices " ++ show (min a b, max a b) ++ + " not present in triangle map " ++ show tm) + <$> [0..maxInd] <*> [0..maxInd] + + vec :: Int -> Int -> Vec ('S ('S 'Z)) Int + vec a b = min a b `VCons` (max a b `VCons` VNil) + +sym2AssocsFac :: forall (n :: Nat) v.Fractional v => Sing n -> [(Vec ('S ('S ('S 'Z))) Int, v)] +sym2AssocsFac sn = assocs + where + n = withKnownNat sn (natVal sn) + tm = trianMapSym2 n + fm = facMapSym2 n + maxInd = fromIntegral (withKnownNat sn (natVal sn)) - (1 :: Int) + assocs = (\a b -> case + (do + let v = vec a b + i <- Map.lookup v tm + f <- Map.lookup v fm + return (a `VCons` (b `VCons` (i `VCons` VNil)), 1 / f :: v)) of + Just x -> x + Nothing -> error "sym2AssocsFac are not fraction-free, as they should be!") + <$> [0..maxInd] <*> [0..maxInd] + + vec :: Int -> Int -> Vec ('S ('S 'Z)) Int + vec a b = min a b `VCons` (max a b `VCons` VNil) + +gamma' :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) (r :: Rank) v. + ( + '[ '( 'VSpace id n, 'Cov (a ':| '[b])) ] ~ r, + (a < b) ~ 'True, + SingI n, + Num v + ) => + Sing id -> Sing n -> Sing a -> Sing b -> + Tensor r v +gamma' _ _ _ _ = gamma + +gamma :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) (r :: Rank) v. + ( + '[ '( 'VSpace id n, 'Cov (a ':| '[b])) ] ~ r, + (a < b) ~ 'True, + SingI n, + Num v + ) => Tensor r v +gamma = case (sing :: Sing n) of + sn -> let x = fromIntegral $ withKnownNat sn $ natVal sn + in Tensor (f x) + where + f :: Int -> [(Int, Tensor (Tail' r) v)] + f x = map (\i -> (i, Tensor [(i, Scalar 1)])) [0..x - 1] + +eta' :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) (r :: Rank) v. + ( + '[ '( 'VSpace id n, 'Cov (a ':| '[b])) ] ~ r, + (a < b) ~ 'True, + SingI n, + Num v + ) => + Sing id -> Sing n -> Sing a -> Sing b -> + Tensor r v +eta' _ _ _ _ = eta + +eta :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) (r :: Rank) v. + ( + '[ '( 'VSpace id n, 'Cov (a ':| '[b])) ] ~ r, + (a < b) ~ 'True, + SingI n, + Num v + ) => Tensor r v +eta = case (sing :: Sing n) of + sn -> let x = fromIntegral $ withKnownNat sn $ natVal sn + in Tensor (f x) + where + f :: Int -> [(Int, Tensor (Tail' r) v)] + f x = map (\i -> (i, Tensor [(i, Scalar (if i == 0 then 1 else -1))])) [0..x - 1] + +gammaInv' :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) (r :: Rank) v. + ( + '[ '( 'VSpace id n, 'Con (a ':| '[b])) ] ~ r, + (a < b) ~ 'True, + SingI n, + Num v + ) => + Sing id -> Sing n -> Sing a -> Sing b -> + Tensor r v +gammaInv' _ _ _ _ = gammaInv + +gammaInv :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) (r :: Rank) v. + ( + '[ '( 'VSpace id n, 'Con (a ':| '[b])) ] ~ r, + (a < b) ~ 'True, + SingI n, + Num v + ) => Tensor r v +gammaInv = case (sing :: Sing n) of + sn -> let x = fromIntegral $ withKnownNat sn $ natVal sn + in Tensor (f x) + where + f :: Int -> [(Int, Tensor (Tail' r) v)] + f x = map (\i -> (i, Tensor [(i, Scalar 1)])) [0..x - 1] + +etaInv' :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) (r :: Rank) v. + ( + '[ '( 'VSpace id n, 'Con (a ':| '[b])) ] ~ r, + (a < b) ~ 'True, + SingI n, + Num v + ) => + Sing id -> Sing n -> Sing a -> Sing b -> + Tensor r v +etaInv' _ _ _ _ = etaInv + +etaInv :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) (r :: Rank) v. + ( + '[ '( 'VSpace id n, 'Con (a ':| '[b])) ] ~ r, + (a < b) ~ 'True, + SingI n, + Num v + ) => Tensor r v +etaInv = case (sing :: Sing n) of + sn -> let x = fromIntegral $ withKnownNat sn $ natVal sn + in Tensor (f x) + where + f :: Int -> [(Int, Tensor (Tail' r) v)] + f x = map (\i -> (i, Tensor [(i, Scalar (if i == 0 then 1 else -1))])) [0..x - 1] + +injSym2Con' :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) + (i :: Symbol) (r :: Rank) v. + ( + InjSym2ConRank id n a b i ~ 'Just r, + SingI r, + Num v + ) => Sing id -> Sing n -> Sing a -> Sing b -> Sing i -> Tensor r v +injSym2Con' _ svdim _ _ _ = + case sSane sr %~ STrue of + Proved Refl -> + case sLengthR sr of + SS (SS (SS SZ)) -> fromList $ sym2Assocs svdim + where + sr = sing :: Sing r + +injSym2Cov' :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) + (i :: Symbol) (r :: Rank) v. + ( + InjSym2CovRank id n a b i ~ 'Just r, + SingI r, + Num v + ) => Sing id -> Sing n -> Sing a -> Sing b -> Sing i -> Tensor r v +injSym2Cov' _ svdim _ _ _ = + case sSane sr %~ STrue of + Proved Refl -> + case sLengthR sr of + SS (SS (SS SZ)) -> fromList $ sym2Assocs svdim + where + sr = sing :: Sing r + +surjSym2Con' :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) + (i :: Symbol) (r :: Rank) v. + ( + SurjSym2ConRank id n a b i ~ 'Just r, + SingI r, + Fractional v + ) => Sing id -> Sing n -> Sing a -> Sing b -> Sing i -> Tensor r v +surjSym2Con' _ svdim _ _ _ = + case sSane sr %~ STrue of + Proved Refl -> + case sLengthR sr of + SS (SS (SS SZ)) -> fromList $ sym2AssocsFac svdim + where + sr = sing :: Sing r + +surjSym2Cov' :: forall (id :: Symbol) (n :: Nat) (a :: Symbol) (b :: Symbol) + (i :: Symbol) (r :: Rank) v. + ( + SurjSym2CovRank id n a b i ~ 'Just r, + SingI r, + Fractional v + ) => Sing id -> Sing n -> Sing a -> Sing b -> Sing i -> Tensor r v +surjSym2Cov' _ svdim _ _ _ = + case sSane sr %~ STrue of + Proved Refl -> + case sLengthR sr of + SS (SS (SS SZ)) -> fromList $ sym2AssocsFac svdim + where + sr = sing :: Sing r + +someGamma :: (Num v, MonadError String m) => + Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> + m (T v) +someGamma vid vdim a b + | a > b = someGamma vid vdim b a + | a == b = throwError $ "cannot construct gamma with indices " ++ show vid ++ " " ++ show vdim ++ " " ++ show a ++ " " ++ show b + | otherwise = + withSomeSing vid $ \svid -> + withSomeSing vdim $ \svdim -> + withSomeSing a $ \sa -> + withSomeSing b $ \sb -> + withKnownNat svdim $ + withKnownSymbol svid $ + withKnownSymbol sa $ + withKnownSymbol sb $ + case sCompare sa sb of + SLT -> return $ T $ gamma' svid svdim sa sb + +someGammaInv :: (Num v, MonadError String m) => + Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> + m (T v) +someGammaInv vid vdim a b + | a > b = someGammaInv vid vdim b a + | a == b = throwError $ "cannot construct gamma with indices " ++ show vid ++ " " ++ show vdim ++ " " ++ show a ++ " " ++ show b + | otherwise = + withSomeSing vid $ \svid -> + withSomeSing vdim $ \svdim -> + withSomeSing a $ \sa -> + withSomeSing b $ \sb -> + withKnownNat svdim $ + withKnownSymbol svid $ + withKnownSymbol sa $ + withKnownSymbol sb $ + case sCompare sa sb of + SLT -> return $ T $ gammaInv' svid svdim sa sb + +someEta :: (Num v, MonadError String m) => + Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> + m (T v) +someEta vid vdim a b + | a > b = someEta vid vdim b a + | a == b = throwError $ "cannot construct eta with indices " ++ show vid ++ " " ++ show vdim ++ " " ++ show a ++ " " ++ show b + | otherwise = + withSomeSing vid $ \svid -> + withSomeSing vdim $ \svdim -> + withSomeSing a $ \sa -> + withSomeSing b $ \sb -> + withKnownNat svdim $ + withKnownSymbol svid $ + withKnownSymbol sa $ + withKnownSymbol sb $ + case sCompare sa sb of + SLT -> return $ T $ eta' svid svdim sa sb + +someEtaInv :: (Num v, MonadError String m) => + Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> + m (T v) +someEtaInv vid vdim a b + | a > b = someEtaInv vid vdim b a + | a == b = throwError $ "cannot construct eta with indices " ++ show vid ++ " " ++ show vdim ++ " " ++ show a ++ " " ++ show b + | otherwise = + withSomeSing vid $ \svid -> + withSomeSing vdim $ \svdim -> + withSomeSing a $ \sa -> + withSomeSing b $ \sb -> + withKnownNat svdim $ + withKnownSymbol svid $ + withKnownSymbol sa $ + withKnownSymbol sb $ + case sCompare sa sb of + SLT -> return $ T $ etaInv' svid svdim sa sb + +someInjSym2Con :: (Num v, MonadError String m) => + Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + m (T v) +someInjSym2Con vid dim a b i + | a > b = someInjSym2Con vid dim b a i + | a == b = throwError $ "Invalid spacetime index for sym2 intertwiner: " ++ show a ++ " " ++ show b ++ "!" + | otherwise = + withSomeSing vid $ \svid -> + withSomeSing dim $ \sdim -> + withSomeSing a $ \sa -> + withSomeSing b $ \sb -> + withSomeSing i $ \si -> + case sInjSym2ConRank svid sdim sa sb si of + SJust sr -> + withSingI sr $ + case sSane sr %~ STrue of + Proved Refl -> return $ T $ injSym2Con' svid sdim sa sb si + +someInjSym2Cov :: (Num v, MonadError String m) => + Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + m (T v) +someInjSym2Cov vid dim a b i + | a > b = someInjSym2Cov vid dim b a i + | a == b = throwError $ "Invalid spacetime index for sym2 intertwiner: " ++ show a ++ " " ++ show b ++ "!" + | otherwise = + withSomeSing vid $ \svid -> + withSomeSing dim $ \sdim -> + withSomeSing a $ \sa -> + withSomeSing b $ \sb -> + withSomeSing i $ \si -> + case sInjSym2CovRank svid sdim sa sb si of + SJust sr -> + withSingI sr $ + case sSane sr %~ STrue of + Proved Refl -> return $ T $ injSym2Cov' svid sdim sa sb si + +someSurjSym2Con :: (Fractional v, MonadError String m) => + Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + m (T v) +someSurjSym2Con vid dim a b i + | a > b = someSurjSym2Con vid dim b a i + | a == b = throwError $ "Invalid spacetime index for sym2 intertwiner: " ++ show a ++ " " ++ show b ++ "!" + | otherwise = + withSomeSing vid $ \svid -> + withSomeSing dim $ \sdim -> + withSomeSing a $ \sa -> + withSomeSing b $ \sb -> + withSomeSing i $ \si -> + case sSurjSym2ConRank svid sdim sa sb si of + SJust sr -> + withSingI sr $ + case sSane sr %~ STrue of + Proved Refl -> return $ T $ surjSym2Con' svid sdim sa sb si + +someSurjSym2Cov :: (Fractional v, MonadError String m) => + Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + m (T v) +someSurjSym2Cov vid dim a b i + | a > b = someSurjSym2Cov vid dim b a i + | a == b = throwError $ "Invalid spacetime index for sym2 intertwiner: " ++ show a ++ " " ++ show b ++ "!" + | otherwise = + withSomeSing vid $ \svid -> + withSomeSing dim $ \sdim -> + withSomeSing a $ \sa -> + withSomeSing b $ \sb -> + withSomeSing i $ \si -> + case sSurjSym2CovRank svid sdim sa sb si of + SJust sr -> + withSingI sr $ + case sSane sr %~ STrue of + Proved Refl -> return $ T $ surjSym2Cov' svid sdim sa sb si + +someInterSym2Con :: Num v => + Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + T v +someInterSym2Con vid dim m n a b = t + where + Right t = runExcept $ + do + (j :: T (Ratio Int)) <- someSurjSym2Con vid dim " " n a + (i :: T (Ratio Int)) <- someInjSym2Con vid dim " " m b + prod <- i .* j + let res = contractT $ fmap ((-2) *) prod + return $ fmap (\v -> if denominator v == 1 + then fromIntegral (numerator v) + else error "someInterSym2Con is not fraction-free, as it should be!") res + +someInterSym2Cov :: Num v => + Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> Demote Symbol -> Demote Symbol -> + T v +someInterSym2Cov vid dim m n a b = t + where + Right t = runExcept $ + do + (j :: T (Ratio Int)) <- someSurjSym2Cov vid dim " " m a + (i :: T (Ratio Int)) <- someInjSym2Cov vid dim " " n b + prod <- i .* j + let res = contractT $ fmap (2*) prod + return $ fmap (\v -> if denominator v == 1 + then fromIntegral (numerator v) + else error "someInterSym2Cov is not fraction-free, as it should be!") res + +someDeltaSym2 :: Num v => Demote Symbol -> Demote Nat -> Demote Symbol -> Demote Symbol -> T v +someDeltaSym2 vid n = someDelta (vid <> "Sym2") ((n*(n+1)) `div` 2) diff --git a/src/Math/Tensor/Basic/TH.hs b/src/Math/Tensor/Basic/TH.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/Basic/TH.hs @@ -0,0 +1,141 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE DefaultSignatures #-} +{-# LANGUAGE EmptyCase #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE InstanceSigs #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE OverloadedStrings #-} + +{-# LANGUAGE CPP #-} +#if MIN_VERSION_base(4,14,0) +{-# LANGUAGE StandaloneKindSignatures #-} +#endif + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.Basic.TH +Description : Template Haskell for Math.Tensor.Basic +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Template Haskell for 'Math.Tensor.Basic'. +-} +----------------------------------------------------------------------------- +module Math.Tensor.Basic.TH where + +import Math.Tensor.Safe.TH + +import Data.Singletons.Prelude +import Data.Singletons.Prelude.Enum +import Data.Singletons.Prelude.List.NonEmpty hiding (sLength) +import Data.Singletons.TH +import Data.Singletons.TypeLits + +import Data.List.NonEmpty (NonEmpty((:|))) + +$(singletons [d| + + -- ############# + -- ### delta ### + -- ############# + deltaRank :: Symbol -> Nat -> Symbol -> Symbol -> Rank + deltaRank vid vdim a b = [(VSpace vid vdim, ConCov (a :| []) (b :| []))] + + -- ############### + -- ### epsilon ### + -- ############### + epsilonRank :: Symbol -> Nat -> NonEmpty Symbol -> Maybe Rank + epsilonRank vid vdim is = + case isLengthNE is vdim of + True -> + case isAscendingNE is of + True -> Just [(VSpace vid vdim, Cov is)] + False -> Nothing + False -> Nothing + + epsilonInvRank :: Symbol -> Nat -> NonEmpty Symbol -> Maybe Rank + epsilonInvRank vid vdim is = + case isLengthNE is vdim of + True -> + case isAscendingNE is of + True -> Just [(VSpace vid vdim, Con is)] + False -> Nothing + False -> Nothing + + -- ############ + -- ### sym2 ### + -- ############ + sym2Dim :: Nat -> Nat + sym2Dim = go 0 + where + go :: Nat -> Nat -> Nat + go acc n = case n == 0 of + True -> acc + False -> go (acc + n) (pred n) + + injSym2ConRank :: Symbol -> Nat -> Symbol -> Symbol -> Symbol -> Maybe Rank + injSym2ConRank vid vdim a b i = + let r = [(VSpace vid vdim, Con (a :| [b])), (VSpace (vid <> "Sym2") (sym2Dim vdim), Cov (i :| []))] + in case sane r of + True -> Just r + False -> Nothing + + injSym2CovRank :: Symbol -> Nat -> Symbol -> Symbol -> Symbol -> Maybe Rank + injSym2CovRank vid vdim a b i = + let r = [(VSpace vid vdim, Cov (a :| [b])), (VSpace (vid <> "Sym2") (sym2Dim vdim), Con (i :| []))] + in case sane r of + True -> Just r + False -> Nothing + + surjSym2ConRank :: Symbol -> Nat -> Symbol -> Symbol -> Symbol -> Maybe Rank + surjSym2ConRank = injSym2CovRank + + surjSym2CovRank :: Symbol -> Nat -> Symbol -> Symbol -> Symbol -> Maybe Rank + surjSym2CovRank = injSym2ConRank + + -- ############ + -- ### Area ### + -- ############ + injAreaConRank :: Symbol -> Symbol -> Symbol -> Symbol -> Symbol -> Symbol -> Maybe Rank + injAreaConRank vid a b c d i = + let r = [(VSpace vid (4 :: Nat), Con (a :| [b,c,d])), (VSpace (vid <> "Area") (21 :: Nat), Cov (i :| []))] + in case sane r of + True -> Just r + False -> Nothing + + injAreaCovRank :: Symbol -> Symbol -> Symbol -> Symbol -> Symbol -> Symbol -> Maybe Rank + injAreaCovRank vid a b c d i = + let r = [(VSpace vid (4 :: Nat), Cov (a :| [b,c,d])), (VSpace (vid <> "Area") (21 :: Nat), Con (i :| []))] + in case sane r of + True -> Just r + False -> Nothing + + surjAreaConRank :: Symbol -> Symbol -> Symbol -> Symbol -> Symbol -> Symbol -> Maybe Rank + surjAreaConRank vid a b c d i = + let r = [(VSpace vid (4 :: Nat), Cov (a :| [b,c,d])), (VSpace (vid <> "Area") (21 :: Nat), Con (i :| []))] + in case sane r of + True -> Just r + False -> Nothing + + surjAreaCovRank :: Symbol -> Symbol -> Symbol -> Symbol -> Symbol -> Symbol -> Maybe Rank + surjAreaCovRank vid a b c d i = + let r = [(VSpace vid (4 :: Nat), Con (a :| [b,c,d])), (VSpace (vid <> "Area") (21 :: Nat), Cov (i :| []))] + in case sane r of + True -> Just r + False -> Nothing + + |]) diff --git a/src/Math/Tensor/LinearAlgebra.hs b/src/Math/Tensor/LinearAlgebra.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/LinearAlgebra.hs @@ -0,0 +1,47 @@ +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.LinearAlgebra +Description : Linear algebra for tensor equations. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Linear algebra for tensor equations. +-} +----------------------------------------------------------------------------- +module Math.Tensor.LinearAlgebra + ( -- * Linear combinations and polynomials + -- ** Data types + Lin(..) + , Poly(..) + , -- ** Construction, inspection, modification + singletonPoly + , polyMap + , getVars + , shiftVars + , normalize + , -- * Tensor equations + -- ** Extracting tensor equations and matrix representations + Equation + , tensorToEquations + , tensorsToSparseMat + , tensorsToMat + , -- ** Rank of a linear tensor equation system + systemRank + , -- ** Solutions + Solution + , solveTensor + , solveSystem + , redefineIndets + , -- ** Internals + equationFromRational + , equationsToSparseMat + , equationsToMat + , fromRref + , fromRow + , applySolution + ) where + +import Math.Tensor.LinearAlgebra.Scalar +import Math.Tensor.LinearAlgebra.Equations diff --git a/src/Math/Tensor/LinearAlgebra/Equations.hs b/src/Math/Tensor/LinearAlgebra/Equations.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/LinearAlgebra/Equations.hs @@ -0,0 +1,227 @@ +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.LinearAlgebra.Equations +Description : Linear tensor equations. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Linear tensor equations. +-} +----------------------------------------------------------------------------- +module Math.Tensor.LinearAlgebra.Equations + ( Equation + , tensorToEquations + , equationFromRational + , equationsToSparseMat + , equationsToMat + , tensorsToSparseMat + , tensorsToMat + , systemRank + , Solution + , fromRref + , fromRow + , applySolution + , solveTensor + , solveSystem + , redefineIndets + ) where + +import Math.Tensor + ( T + , removeZerosT + , toListT + ) +import Math.Tensor.LinearAlgebra.Scalar + ( Poly (Const, Affine, NotSupported) + , Lin (Lin) + , polyMap + , normalize + ) +import Math.Tensor.LinearAlgebra.Matrix + ( rref + , isrref + , verify + ) + +import qualified Numeric.LinearAlgebra.Data as HM + ( Matrix + , R + , Z + , fromLists + , toLists + ) +import Numeric.LinearAlgebra (rank) + +import Data.Maybe (mapMaybe) +import qualified Data.IntMap.Strict as IM + ( IntMap + , foldl' + , map + , assocs + , null + , findWithDefault + , lookupMax + , keys + , fromList + , (!) + , difference + , intersectionWith + , mapKeys + ) +import Data.List (nub, sort) +import Data.Ratio (numerator, denominator) + +-- |A linear equation is a mapping from variable +-- indices to coefficients +type Equation a = IM.IntMap a + +-- |Extract linear equations from tensor components. +-- The equations are normalized, sorted, and made unique. +tensorToEquations :: Integral a => T (Poly Rational) -> [Equation a] +tensorToEquations = nub . sort . fmap (equationFromRational . normalize . snd) . toListT + +-- |Extract linear equation with integral coefficients from polynomial +-- tensor component with rational coefficients. +-- Made made integral by multiplying with the @lcm@ of all denominators. +equationFromRational :: forall a.Integral a => Poly Rational -> Equation a +equationFromRational (Affine x (Lin lin)) + | x == 0 = lin' + | otherwise = error "affine equation not supported for the moment!" + where + fac :: a + fac = IM.foldl' (\acc v -> lcm (fromIntegral (denominator v)) acc) 1 lin + lin' = IM.map (\v -> fromIntegral (numerator v) * (fac `div` fromIntegral (denominator v))) lin +equationFromRational _ = error "equation can only be extracted from linear scalar!" + +-- |Convert list of equations to sparse matrix representation of the +-- linear system. +equationsToSparseMat :: [Equation a] -> [((Int,Int), a)] +equationsToSparseMat xs = concat $ zipWith (\i m -> fmap (\(j,v) -> ((i,j),v)) (IM.assocs m)) [1..] xs + +-- |Convert list of equations to dense matrix representation of the +-- linear system. +equationsToMat :: Integral a => [Equation a] -> [[a]] +equationsToMat eqns = mapMaybe (\m -> if IM.null m + then Nothing + else Just $ fmap (\j -> IM.findWithDefault 0 j m) [1..maxVar]) eqns + where + maxVar = maximum $ mapMaybe (fmap fst . IM.lookupMax) eqns + +-- |Extract sparse matrix representation for the linear system given +-- by a list of existentially quantified tensors with polynomial values. +tensorsToSparseMat :: Integral a => [T (Poly Rational)] -> [((Int,Int), a)] +tensorsToSparseMat = equationsToSparseMat . concatMap tensorToEquations + +-- |Extract dense matrix representation for the linear system given +-- by a list of existentially quantified tensors with polynomial values. +tensorsToMat :: Integral a => [T (Poly Rational)] -> [[a]] +tensorsToMat = equationsToMat . concatMap tensorToEquations + +-- |Rank of a linear system. Uses dense svd provided by hmatrix. +matRank :: forall a.Integral a => [[a]] -> Int +matRank [] = 0 +matRank mat = rank (hmat :: HM.Matrix HM.R) + where + hmat = HM.fromLists $ fmap (fmap (fromIntegral @a @HM.R)) mat + +-- |Rank of the linear system given by a list of existentially +-- quantified tensors with polynomial values. +systemRank :: [T (Poly Rational)] -> Int +systemRank = matRank . tensorsToMat @Int + +-- |The solution to a linear system is represented as a list of +-- substitution rules, stored as @'IM.IntMap' ('Poly' 'Rational')@. +type Solution = IM.IntMap (Poly Rational) + +-- |Read substitution rules from reduced row echelon form +-- of a linear system. +fromRref :: HM.Matrix HM.Z -> Solution +fromRref ref = IM.fromList assocs + where + rows = HM.toLists ref + assocs = mapMaybe fromRow rows + +-- |Read single substitution rule from single +-- row of reduced row echelon form. +fromRow :: forall a.Integral a => [a] -> Maybe (Int, Poly Rational) +fromRow xs = case assocs of + [] -> Nothing + [(i,_)] -> Just (i, Const 0) + (i, v):assocs' -> let assocs'' = fmap (\(i',v') -> (i', - fromIntegral @a @Rational v' / fromIntegral @a @Rational v)) assocs' + in Just (i, Affine 0 (Lin (IM.fromList assocs''))) + where + assocs = filter ((/=0). snd) $ zip [(1::Int)..] xs + +-- |Apply substitution rules to tensor component. +applySolution :: Solution -> Poly Rational -> Poly Rational +applySolution s (Affine x (Lin lin)) + | x == 0 = case p of + Affine xFin (Lin linFin) -> if IM.null linFin + then Const xFin + else p + _ -> p + | otherwise = error "affine equations not yet supported" + where + s' = IM.intersectionWith (\row v -> if v == 0 + then error "value 0 encountered in linear scalar" + else polyMap (v*) row) s lin + lin' = IM.difference lin s + p0 = if IM.null lin' + then Const 0 + else Affine 0 (Lin lin') + + p = IM.foldl' (+) p0 s' + + {- + lin' = IM.foldlWithKey' + (\lin' i sub -> + case Affine 0 (Lin lin') + sub of + Affine 0 (Lin lin'') -> IM.delete i lin'' + Const 0 -> IM.empty + _ -> error "affine equations not yet supported") + lin s' + -} +applySolution _ _ = error "only linear equations supported" + +-- |Apply substitution rules to all components of a tensor. +solveTensor :: Solution -> T (Poly Rational) -> T (Poly Rational) +solveTensor sol = removeZerosT . fmap (applySolution sol) + +-- |Solve a linear system and apply solution to the tensorial +-- indeterminants. +solveSystem :: + [T (Poly Rational)] -- ^ Tensorial linear system + -> [T (Poly Rational)] -- ^ List of indeterminant tensors + -> [T (Poly Rational)] -- ^ Solved indeterminant tensors +solveSystem system indets + | wrongSolution = error "Wrong solution found. May be an Int64 overflow." + | otherwise = indets' + where + mat = HM.fromLists $ tensorsToMat @HM.Z system + ref = rref mat + wrongSolution = not (isrref ref && verify mat ref) + sol = fromRref ref + indets' = fmap (solveTensor sol) indets + +-- |Relabelling of the indeterminants present in a list of tensors. +-- Redefines the labels of @n@ indeterminants as @[1..n]@, preserving +-- the previous order. +redefineIndets :: [T (Poly v)] -> [T (Poly v)] +redefineIndets indets = fmap (fmap (\case + Const c -> Const c + NotSupported -> NotSupported + Affine a (Lin lin) -> + Affine a (Lin (IM.mapKeys (varMap IM.!) lin)))) indets + where + comps = snd <$> concatMap toListT indets + vars = nub $ concat $ mapMaybe (\case + Affine _ (Lin lin) -> Just $ IM.keys lin + _ -> Nothing) comps + varMap = IM.fromList $ zip vars [(1::Int)..] diff --git a/src/Math/Tensor/LinearAlgebra/Matrix.hs b/src/Math/Tensor/LinearAlgebra/Matrix.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/LinearAlgebra/Matrix.hs @@ -0,0 +1,440 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE FlexibleContexts #-} + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.LinearAlgebra.Matrix +Description : Gaussian elimination subroutines. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Gaussian elimination subroutines. +-} +----------------------------------------------------------------------------- +module Math.Tensor.LinearAlgebra.Matrix + ( gaussianST + , gaussianFFST + , gaussian + , gaussianFF + , rrefST + , rref + , independentColumns + , independentColumnsFF + , independentColumnsRREF + , independentColumnsVerifiedFF + , independentColumnsMat + , independentColumnsMatFF + , independentColumnsMatRREF + , pivotsU + , pivotsUFF + , findPivotMax + , findPivotMaxFF + , findRowPivot + , isref + , isrref + , isrref' + , verify + ) where + +import Numeric.LinearAlgebra + ( Matrix + , Vector + , Container + , Extractor (All, Take, Drop) + , Z + , toLists + , rows + , cols + , find + , (¿) + , (??) + , (><) + , (===) + , rank + , fromZ + ) +import Numeric.LinearAlgebra.Devel + ( STMatrix + , RowOper (AXPY, SCAL, SWAP) + , ColRange (FromCol) + , RowRange (Row) + , freezeMatrix + , thawMatrix + , modifyMatrix + , readMatrix + , rowOper + ) + +import Data.List (maximumBy) + +import Control.Monad (foldM) +import Control.Monad.ST + ( ST + , runST + ) + +-- | 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 e (i,j) of + Nothing -> [] + Just (i', j') -> j' : go (i'+1, j'+1) + maxAbs = maximum $ map (maximum . map abs) $ toLists mat + e = eps * maxAbs + +pivotsUFF :: Matrix Z -> [Int] +pivotsUFF mat = go (0,0) + where + go (i,j) + = case findPivotFF mat (i,j) of + Nothing -> [] + Just (i', j') -> j' : go (i'+1, j'+1) + +eps :: Double +eps = 1e-12 + +-- find next pivot in upper triangular matrix +findPivotFF :: Matrix Z -> (Int, Int) -> Maybe (Int, Int) +findPivotFF mat (i, j) + | n == j = Nothing + | m == i = Nothing + | otherwise = case nonZeros of + [] -> if n == j+1 + then Nothing + else findPivotFF 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 (/= 0) col + +findPivot :: Matrix Double -> Double -> (Int, Int) -> Maybe (Int, Int) +findPivot mat e (i, j) + | n == j = Nothing + | m == i = Nothing + | otherwise = case nonZeros of + [] -> if n == j+1 + then Nothing + else findPivot mat e (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 ((>= e) . abs) col + +-- | Find pivot element below position (i, j) with greatest absolute value in the ST monad. +findPivotMaxFF :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe (Int, Int)) +findPivotMaxFF m n i j mat + | n == j = return Nothing + | m == i = return Nothing + | otherwise = + do + col <- mapM (\i' -> do + x <- readMatrix mat i' j + return (i', x)) + [i..m-1] + let nonZeros = filter ((/= 0) . snd) col + case nonZeros of + [] -> if n == j+1 + then return Nothing + else findPivotMaxFF m n i (j+1) mat + (pi_,_):_ -> return $ Just (pi_, j) + +findPivotMax :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s (Maybe (Int, Int)) +findPivotMax m n i j mat + | n == j = return Nothing + | m == i = return Nothing + | otherwise = + do + col <- mapM (\i' -> do + x <- readMatrix mat i' j + return (i', abs x)) + [i..m-1] + let nonZeros = filter ((>= eps) . abs . snd) col + let (pi_, _) = maximumBy (\(_, x) (_, y) -> x `compare` y) nonZeros + case nonZeros of + [] -> if n == j+1 + then return Nothing + else findPivotMax m n i (j+1) mat + _ -> return $ Just (pi_, j) + +findRowPivot :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe Int) +findRowPivot m n i j mat + | j + 1 > n = error "out of bounds" -- return Nothing + | i + 1 > min m n = error "out of bounds" -- return Nothing + | otherwise = + do + row <- mapM (\j' -> do + x <- readMatrix mat i j' + return (j', x)) + [0 .. j] + let nonZeros = filter ((/=0) . snd) row + case nonZeros of + [] -> return Nothing + (pj, _):_ -> return $ Just pj + +backwardFF' :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s () +backwardFF' m n i j mat + | i == 0 = return () + | otherwise = do + iPivot' <- findRowPivot m n i j mat + case iPivot' of + Nothing -> backwardFF' m n (i-1) j mat + Just pj -> do + pv <- readMatrix mat i pj + mapM_ (reduce pv pj) [0 .. i-1] + backwardFF' m n (i-1) (pj-1) mat + where + reduce pv pj r = do + Just pr <- findRowPivot m n r pj mat + -- prv <- readMatrix mat r pr + pjv <- readMatrix mat r pj + if pjv == 0 + then return () + else + let op1 = SCAL pv (Row r) (FromCol pr) + op2 = AXPY (-pjv) i r (FromCol pj) + in do + rowOper op1 mat + rowOper op2 mat + g <- foldM (\acc c -> gcd acc <$> readMatrix mat r c) 0 [pr .. n-1] + if g == 0 + then return() + else mapM_ (\c -> modifyMatrix mat r c (`quot` g)) [pr .. n-1] + +gaussianFF' :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s () +gaussianFF' m n i j mat = do + iPivot' <- findPivotMaxFF m n i j mat + case iPivot' of + Nothing -> return () + Just (r, p) -> do + rowOper (SWAP i r (FromCol j)) mat + pv <- readMatrix mat i p + mapM_ (reduce pv p) [i+1 .. m-1] + gaussianFF' m n (i+1) (p+1) mat + where + reduce pv p r = do + rv <- readMatrix mat r p + if rv == 0 + then return () + else + let op1 = SCAL pv (Row r) (FromCol p) + op2 = AXPY (-rv) i r (FromCol p) + in do + rowOper op1 mat + rowOper op2 mat + g <- foldM (\acc c -> gcd acc <$> readMatrix mat r c) 0 [p .. n-1] + if g == 0 + then return() + else mapM_ (\c -> modifyMatrix mat r c (`quot` g)) [p .. n-1] + +-- gaussian elimination of sub matrix below position (i, j) +gaussian' :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s () +gaussian' m n i j mat = do + iPivot' <- findPivotMax m n i j mat + case iPivot' of + Nothing -> return () + Just (r, p) -> do + rowOper (SWAP i r (FromCol j)) mat + pv <- readMatrix mat i p + mapM_ (reduce pv p) [i+1 .. m-1] + gaussian' m n (i+1) (p+1) mat + where + reduce pv p r = do + rv <- readMatrix mat r p + if abs rv < eps + then return () + else + let frac = -rv / pv + op = AXPY frac i r (FromCol p) + in do + rowOper op mat + mapM_ (\j' -> modifyMatrix mat r j' (\x -> if abs x < eps then 0 else x)) [p..n-1] + +gaussianFFST :: Int -> Int -> STMatrix s Z -> ST s () +gaussianFFST m n = gaussianFF' m n 0 0 + +-- | Gaussian elimination perfomed in-place in the @'ST'@ monad. +gaussianST :: Int -> Int -> STMatrix s Double -> ST s () +gaussianST m n = gaussian' m n 0 0 + +rrefST :: Int -> Int -> STMatrix s Z -> ST s () +rrefST m n mat = do + gaussianFF' m n 0 0 mat + backwardFF' m n (r'-1) (n-1) mat + where + r' = min m n + +isref :: (Num a, Ord a, Container Vector a) => Matrix a -> Bool +isref mat = case pivot of + [] -> True + (r,p):_ -> (r <= 0) + && + (let leftMat = mat ?? (Drop 1, Take (p+1)) + rightMat = mat ?? (Drop 1, Drop (p+1)) + leftZero = null $ find (/=0) leftMat + rightRef = isref rightMat + in leftZero && rightRef) + where + pivot = find (/=0) mat + +isrref' :: (Num a, Ord a, Container Vector a) => Int -> Matrix a -> Bool +isrref' r mat = case pivot of + [] -> True + (r',p):_ -> (r' <= 0) + && (let leftMat = subMat ?? (Drop 1, Take (p+1)) + col = mat ¿ [p] + colSingleton = case find (/=0) col of + [_] -> True + _ -> False + leftZero = null $ find (/=0) leftMat + nextRref = isrref' (r+1) mat + in leftZero && colSingleton && nextRref) + where + subMat = mat ?? (Drop r, All) + pivot = find (/=0) subMat + +isrref :: (Num a, Ord a, Container Vector a) => Matrix a -> Bool +isrref = isrref' 0 + +rref :: Matrix Z -> Matrix Z +rref mat = runST $ do + matST <- thawMatrix mat + rrefST m n matST + freezeMatrix matST + where + m = rows mat + n = cols mat + +gaussianFF :: Matrix Z -> Matrix Z +gaussianFF mat = runST $ do + matST <- thawMatrix mat + gaussianFFST m n matST + freezeMatrix matST + where + m = rows mat + n = cols mat + +-- | 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 + matST <- thawMatrix mat + gaussianST m n matST + freezeMatrix matST + where + m = rows mat + n = cols mat + +independentColumnsRREF :: Matrix Z -> [Int] +independentColumnsRREF mat = pivotsUFF mat' + where + mat' = rref mat + +independentColumnsFF :: Matrix Z -> [Int] +independentColumnsFF mat = pivotsUFF mat' + where + mat' = gaussianFF mat + +independentColumnsVerifiedFF :: Matrix Z -> [Int] +independentColumnsVerifiedFF mat + | isref mat' && verify mat mat' + = pivotsUFF mat' + | otherwise = error "could not verify row echelon form" + where + mat' = gaussianFF mat + +-- | 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 + +verify :: Matrix Z -> Matrix Z -> Bool +verify mat ref = rank1 == rank2 && rank1 == rank3 + where + matD = fromZ mat :: Matrix Double + refD = fromZ ref :: Matrix Double + rank1 = rank matD + rank2 = rank refD + rank3 = rank $ matD === refD + +independentColumnsMatRREF :: Matrix Z -> Matrix Z +independentColumnsMatRREF mat = + case independentColumnsRREF mat of + [] -> (rows mat >< 1) $ repeat 0 + cs -> mat ¿ cs + +independentColumnsMatFF :: Matrix Z -> Matrix Z +independentColumnsMatFF mat = + case independentColumnsFF mat of + [] -> (rows mat >< 1) $ repeat 0 + cs -> mat ¿ cs + +-- | 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/LinearAlgebra/Scalar.hs b/src/Math/Tensor/LinearAlgebra/Scalar.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/LinearAlgebra/Scalar.hs @@ -0,0 +1,117 @@ +{-# LANGUAGE Safe #-} +{-# LANGUAGE MultiParamTypeClasses #-} + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.LinearAlgebra.Scalar +Description : Scalar types for usage as Tensor values. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Scalar types for usage as Tensor values. +-} +----------------------------------------------------------------------------- +module Math.Tensor.LinearAlgebra.Scalar + ( Lin(..) + , Poly(..) + , singletonPoly + , polyMap + , getVars + , shiftVars + , normalize + ) where + +import qualified Data.IntMap.Strict as IM + ( IntMap + , singleton + , null + , keys + , map + , filter + , mapKeysMonotonic + , unionWith + , findMin + ) + +-- |Linear combination represented as mapping from +-- variable number to prefactor. +newtype Lin a = Lin (IM.IntMap a) deriving (Show, Ord, Eq) + +-- |Polynomial: Can be constant, affine, or something of higher +-- rank which is not yet implemented. +data Poly a = Const !a -- ^ constant value + | Affine !a !(Lin a) -- ^ constant value plus linear term + | NotSupported -- ^ higher rank + deriving (Show, Ord, Eq) + +-- |Produces an affine value \(c + a\cdot x_i\) +singletonPoly :: a -- ^ constant + -> Int -- ^ variable number + -> a -- ^ prefactor + -> Poly a +singletonPoly a i v = Affine a $ Lin $ IM.singleton i v + +-- |Maps over 'Poly' +polyMap :: (a -> b) -> Poly a -> Poly b +polyMap f (Const a) = Const (f a) +polyMap f (Affine a (Lin lin)) = Affine (f a) $ Lin $ IM.map f lin +polyMap _ _ = NotSupported + +instance (Num a, Eq a) => Num (Poly a) where + Const a + Const b = Const $ a + b + Const a + Affine b lin = Affine (a+b) lin + Affine a lin + Const b = Affine (a+b) lin + Affine a (Lin m1) + Affine b (Lin m2) + | IM.null m' = Const $ a + b + | otherwise = Affine (a+b) (Lin m') + where + m' = IM.filter (/=0) $ IM.unionWith (+) m1 m2 + NotSupported + _ = NotSupported + _ + NotSupported = NotSupported + + negate = polyMap negate + + abs (Const a) = Const $ abs a + abs _ = NotSupported + + signum (Const a) = Const $ signum a + signum _ = NotSupported + + fromInteger = Const . fromInteger + + Const a * Const b = Const $ a * b + Const a * Affine b (Lin lin) + | a == 0 = Const 0 + | otherwise = Affine (a*b) $ Lin $ IM.map (a*) lin + Affine a (Lin lin) * Const b + | b == 0 = Const 0 + | otherwise = Affine (a*b) $ Lin $ IM.map (*b) lin + _ * _ = NotSupported + +-- |Returns list of variable numbers present in the polynomial. +getVars :: Poly a -> [Int] +getVars (Const _) = [] +getVars NotSupported = [] +getVars (Affine _ (Lin lm)) = IM.keys lm + +-- |Shifts variable numbers in the polynomial by a constant value. +shiftVars :: Int -> Poly a -> Poly a +shiftVars _ (Const a) = Const a +shiftVars _ NotSupported = NotSupported +shiftVars s (Affine a (Lin lin)) = + Affine a $ Lin $ IM.mapKeysMonotonic (+s) lin + +-- |Normalizes a polynomial: +-- \[ +-- \mathrm{normalize}(c) = 1 \\ +-- \mathrm{normalize}(c + a_1\cdot x_1 + a_2\cdot x_2 + \dots + a_n\cdot x_n) = \frac{c}{a_1} + 1\cdot x_1 + \frac{a_2}{a_1}\cdot x_2 + \dots + \frac{a_n}{a_1}\cdot x_n +-- \] +normalize :: Fractional a => Poly a -> Poly a +normalize (Const _) = Const 1 +normalize NotSupported = NotSupported +normalize (Affine a (Lin lin)) = Affine (a/v) $ Lin $ IM.map (/v) lin + where + (_,v) = IM.findMin lin + diff --git a/src/Math/Tensor/Safe.hs b/src/Math/Tensor/Safe.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/Safe.hs @@ -0,0 +1,574 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE StandaloneDeriving #-} + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.Safe +Description : Dependently typed tensor algebra. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Dependently typed tensor algebra. +-} +----------------------------------------------------------------------------- +module Math.Tensor.Safe + ( -- * The Tensor GADT + Tensor(..) + -- * Generic Rank of a Tensor + -- |A vector space is the product of a label and a dimension. + , VSpace(..) + -- |The generic tensor rank is a list of vector spaces with label, dimension and + -- associated index list. + , GRank + -- |The rank of a tensor is a generic rank specialized to 'Symbol' and 'Nat' + , Rank + -- * Length-typed assocs lists + -- |Type-level naturals used internally. + , N(..) + , -- |Length-typed vector used internally. + Vec(..) + , vecFromListUnsafe + , -- * Conversion from and to lists + fromList + , fromList' + , toList + , -- * Tensor algebra + (&+), (&-), (&*), removeZeros + , -- * Contraction + contract + , -- * Transpositions + transpose + , transposeMult + , -- * Relabelling + relabel + ) where + +import Math.Tensor.Safe.TH +import Math.Tensor.Safe.Proofs +import Math.Tensor.Safe.Vector + +import Data.Kind (Type) + +import Data.Constraint (Dict(Dict), (:-)(Sub)) + +import Data.Singletons + ( Sing + , SingI (sing) + , withSingI, fromSing + ) +import Data.Singletons.Prelude + ( SBool (STrue, SFalse) + , SList (SNil) + , SMaybe (SJust) + , SOrdering (SLT, SEQ, SGT) + , STuple2 (STuple2) + , Tail + , sFst, sSnd, sHead, sTail + , sCompare, (%==) + ) +import Data.Singletons.Prelude.Maybe + ( IsJust + , sIsJust + ) +import Data.Singletons.Decide + ( Decision (Proved, Disproved) + , (:~:) (Refl) + , (%~) + ) +import Data.Singletons.TypeLits (Nat, Symbol) + +import Data.Maybe (catMaybes) +import Data.Bifunctor (first,second) +import Data.List (foldl',groupBy,sortBy) + +-- |The 'Tensor' type is parameterized by its generalized 'Rank' @r@ and holds +-- arbitrary values @v@. +data Tensor :: Rank -> Type -> Type where + ZeroTensor :: forall (r :: Rank) v . Sane r ~ 'True => Tensor r v -- ^ + -- A tensor of any sane rank type can be zero. + Scalar :: forall v. !v -> Tensor '[] v -- ^ + -- A tensor of empty rank is a scalar holding some value. + Tensor :: forall (r :: Rank) (r' :: Rank) v. + (Sane r ~ 'True, Tail' r ~ r') => + [(Int, Tensor r' v)] -> Tensor r v -- ^ + -- A non-zero tensor of sane non-empty rank is represented as an assocs list of + -- component-value pairs. The keys must be unique and in ascending order. + -- The values are tensors of the next-lower rank. + +deriving instance Eq v => Eq (Tensor r v) +deriving instance Show v => Show (Tensor r v) + +instance Functor (Tensor r) where + fmap _ ZeroTensor = ZeroTensor + fmap f (Scalar s) = Scalar $ f s + fmap f (Tensor ms) = Tensor $ fmap (fmap (fmap f)) ms + +-- |Union of assocs lists with a merging function if a component is present in both lists +-- and two functions to treat components only present in either list. +unionWith :: (a -> b -> c) -> (a -> c) -> (b -> c) -> [(Int, a)] -> [(Int, b)] -> [(Int, c)] +unionWith _ _ f [] ys = fmap (fmap f) ys +unionWith _ f _ xs [] = fmap (fmap f) xs +unionWith f g h xs@((ix,vx):xs') ys@((iy,vy):ys') = + case ix `compare` iy of + LT -> (ix, g vx) : unionWith f g h xs' ys + EQ -> (ix, f vx vy) : unionWith f g h xs' ys' + GT -> (iy, h vy) : unionWith f g h xs ys' + +-- |Given a 'Num' and 'Eq' instance, remove all zero values from the tensor, +-- eventually replacing a zero @Scalar@ or an empty @Tensor@ with @ZeroTensor@. +removeZeros :: (Num v, Eq v) => Tensor r v -> Tensor r v +removeZeros ZeroTensor = ZeroTensor +removeZeros (Scalar s) = if s == 0 then ZeroTensor else Scalar s +removeZeros (Tensor ms) = + case ms' of + [] -> ZeroTensor + _ -> Tensor ms' + where + ms' = filter + (\(_, t) -> + case t of + ZeroTensor -> False + _ -> True) $ + fmap (fmap removeZeros) ms + +-- |Tensor addition. Ranks of summands and sum coincide. +-- Zero values are removed from the result. +(&+) :: forall (r :: Rank) (r' :: Rank) v. + ((r ~ r'), Num v, Eq v) => + Tensor r v -> Tensor r' v -> Tensor r v +(&+) ZeroTensor t = t +(&+) t ZeroTensor = t +(&+) (Scalar s) (Scalar s') = + if s'' == 0 then ZeroTensor else Scalar s'' + where + s'' = s + s' +(&+) (Tensor xs) (Tensor xs') = removeZeros $ Tensor xs'' + where + xs'' = unionWith (&+) id id xs xs' +(&+) _ _ = error "Cannot add scalar and tensor! Should have been caught by the type system!" + +infixl 6 &+ + +-- |Tensor subtraction. Ranks of operands and difference coincide. +-- Zero values are removed from the result. +(&-) :: forall (r :: Rank) (r' :: Rank) v. + ((r ~ r'), Num v, Eq v) => + Tensor r v -> Tensor r' v -> Tensor r v +(&-) t1 t2 = t1 &+ fmap negate t2 + +infixl 6 &- + +-- |Tensor multiplication, ranks of factors passed explicitly as singletons. +mult :: forall (r :: Rank) (r' :: Rank) (r'' :: Rank) v. + (Num v, 'Just r'' ~ MergeR r r') => + Sing r -> Sing r' -> Tensor r v -> Tensor r' v -> Tensor r'' v +mult _ _ (Scalar s) (Scalar t) = Scalar (s*t) +mult sr sr' (Scalar s) (Tensor ms) = + case saneTail'Proof sr' of + Sub Dict -> Tensor $ fmap (fmap (mult sr (sTail' sr') (Scalar s))) ms +mult sr sr' (Tensor ms) (Scalar s) = + case saneTail'Proof sr of + Sub Dict -> Tensor $ fmap (fmap (\t -> mult (sTail' sr) sr' t (Scalar s))) ms +mult sr sr' (Tensor ms) (Tensor ms') = + let sh = sHead' sr + sh' = sHead' sr' + st = sTail' sr + st' = sTail' sr' + in case saneMergeRProof sr sr' of + Sub Dict -> + case sh of + STuple2 sv si -> + case sh' of + STuple2 sv' si' -> + case sCompare sv sv' of + SLT -> case proofMergeLT sr sr' of + Sub Dict -> + case saneTail'Proof sr of + Sub Dict -> Tensor $ fmap (fmap (\t -> mult st sr' t (Tensor ms'))) ms + SGT -> case proofMergeGT sr sr' of + Sub Dict -> + case saneTail'Proof sr' of + Sub Dict -> Tensor $ fmap (fmap (mult sr st' (Tensor ms))) ms' + SEQ -> case proofMergeIxNotEQ sr sr' of + Sub Dict -> + case sIxCompare si si' of + SLT -> case proofMergeIxLT sr sr' of + Sub Dict -> + case saneTail'Proof sr of + Sub Dict -> Tensor $ fmap (fmap (\t -> mult st sr' t (Tensor ms'))) ms + SGT -> case proofMergeIxGT sr sr' of + Sub Dict -> + case saneTail'Proof sr' of + Sub Dict -> Tensor $ fmap (fmap (mult sr st' (Tensor ms))) ms' +mult sr sr' ZeroTensor ZeroTensor = + case saneMergeRProof sr sr' of + Sub Dict -> ZeroTensor +mult sr sr' ZeroTensor (Scalar _) = + case saneMergeRProof sr sr' of + Sub Dict -> ZeroTensor +mult sr sr' ZeroTensor (Tensor _) = + case saneMergeRProof sr sr' of + Sub Dict -> ZeroTensor +mult sr sr' (Scalar _) ZeroTensor = + case saneMergeRProof sr sr' of + Sub Dict -> ZeroTensor +mult sr sr' (Tensor _) ZeroTensor = + case saneMergeRProof sr sr' of + Sub Dict -> ZeroTensor + +-- |Tensor multiplication. Ranks of factors must not overlap. The Product +-- rank is the merged rank of the factors. +(&*) :: forall (r :: Rank) (r' :: Rank) (r'' :: Rank) v. + (Num v, 'Just r'' ~ MergeR r r', SingI r, SingI r') => + Tensor r v -> Tensor r' v -> Tensor r'' v +(&*) = mult (sing :: Sing r) (sing :: Sing r') + +infixl 7 &* + +contract' :: forall (r :: Rank) (r' :: Rank) v. + (r' ~ ContractR r, Num v, Eq v) + => Sing r -> Tensor r v -> Tensor r' v +contract' sr t = case sContractR sr %~ sr of + Proved Refl -> t + Disproved _ -> contract'' sr t + +contract'' :: forall (r :: Rank) (r' :: Rank) v. + (r' ~ ContractR r, Num v, Eq v) + => Sing r -> Tensor r v -> Tensor r' v +contract'' sr ZeroTensor = + case saneContractProof sr of + Sub Dict -> ZeroTensor +contract'' _ (Scalar v) = Scalar v +contract'' sr (Tensor ms) = + case sTail' sr of + SNil -> + case singletonContractProof sr of + Sub Dict -> Tensor ms + st -> + case saneContractProof sr of + Sub Dict -> + let st' = sTail' st + sh = sHead' sr + sv = sFst sh + si = sSnd sh + sh' = sHead' st + sv' = sFst sh' + si' = sSnd sh' + in case sv %== sv' of + SFalse -> + case contractTailDiffVProof sr of + Sub Dict -> removeZeros $ Tensor $ fmap (fmap (contract'' st)) ms + STrue -> case si of + SICon sa -> case si' of + SICov sb -> case sa %== sb of + STrue -> + let ms' = fmap (\(i, v) -> case v of + Tensor vs -> + case filter (\(i', _) -> i == i') vs of + [] -> Nothing + [(_, v')] -> Just v' + _ -> error "duplicate key in tensor assoc list") ms + ms'' = catMaybes ms' :: [Tensor (Tail' (Tail' r)) v] + in case saneTail'Proof sr of + Sub Dict -> + case saneTail'Proof st of + Sub Dict -> + case contractTailSameVSameIProof sr of + Sub Dict -> contract' st' $ foldl' (&+) ZeroTensor ms'' + SFalse -> + case contractTailSameVDiffIProof sr of + Sub Dict -> removeZeros $ Tensor $ fmap (fmap (contract'' st)) ms + SICon _ -> + case contractTailSameVNoCovProof sr of + Sub Dict -> removeZeros $ Tensor $ fmap (fmap (contract'' st)) ms + SICov _ -> + case contractTailSameVNoConProof sr of + Sub Dict -> removeZeros $ Tensor $ fmap (fmap (contract'' st)) ms + +-- |Tensor contraction. Contracting a tensor is the identity function on non-contractible tensors. +-- Otherwise, the result is the contracted tensor with the contracted labels removed from the rank. +contract :: forall (r :: Rank) (r' :: Rank) v. + (r' ~ ContractR r, SingI r, Num v, Eq v) + => Tensor r v -> Tensor r' v +contract = contract' (sing :: Sing r) + +-- |Tensor transposition. Given a vector space and two index labels, the result is a tensor with +-- the corresponding entries swapped. Only possible if the indices are part of the rank. The +-- rank remains untouched. +transpose :: forall (vs :: VSpace Symbol Nat) (a :: Ix Symbol) (b :: Ix Symbol) (r :: Rank) v. + (CanTranspose vs a b r ~ 'True, SingI r) => + Sing vs -> Sing a -> Sing b -> Tensor r v -> Tensor r v +transpose _ _ _ ZeroTensor = ZeroTensor +transpose _ _ _ (Scalar _) = error "This is not possible, might yet have to convince the type system." +transpose v a b t@(Tensor ms) = + case a `sCompare` b of + SEQ -> t + SGT -> case sCanTranspose v b a (sing :: Sing r) %~ STrue of + Proved Refl -> transpose v b a t + SLT -> + let sr = sing :: Sing r + sh = sHead' sr + sv = sFst sh + si = sSnd sh + st = sTail' sr + in withSingI st $ + case sv %~ v of + Proved Refl -> case si %~ a of + Proved Refl -> let sr' = sRemoveUntil b sr + in withSingI sr' $ + case sSane sr' %~ STrue of + Proved Refl -> + let tl = toTListUntil b t + tl' = fmap (\(i:is, val) -> (last is : (init is ++ [i]),val)) tl + tl'' = sortBy (\(i,_) (i',_) -> i `compare` i') tl' + in fromTList tl'' + Disproved _ -> case sCanTranspose v a b st of + STrue -> Tensor $ fmap (fmap (transpose v a b)) ms + Disproved _ -> case sCanTranspose v a b st of + STrue -> Tensor $ fmap (fmap (transpose v a b)) ms + +-- |Transposition of multiple labels. Given a vector space and a list of transpositions, the +-- result is a tensor with the corresponding entries swapped. Only possible if the indices are +-- part of the rank. The rank remains untouched. +transposeMult :: forall (vs :: VSpace Symbol Nat) (tl :: TransList Symbol) (r :: Rank) v. + (IsJust (Transpositions vs tl r) ~ 'True, SingI r) => + Sing vs -> Sing tl -> Tensor r v -> Tensor r v +transposeMult _ _ ZeroTensor = ZeroTensor +transposeMult sv stl tens@(Tensor ms) = + let sr = sing :: Sing r + sh = sHead' sr + st = sTail' sr + sr' = sTail sr + sts = sTranspositions sv stl sr + in case sv %~ sFst sh of + Proved Refl -> + case sSane sr' %~ STrue of + Proved Refl -> + case sts of + SJust sts' -> + withSingI (sFst (sHead sr)) $ + withSingI sr' $ + let sn = sLengthIL (sSnd (sHead sr)) + n = fromSing sn + ts = fromSing sts' + ts' = go ts $ take' n 0 + xs = toTListWhile tens + xs' = fmap (first (transposeIndices ts')) xs + xs'' = sortBy (\(i,_) (i',_) -> i `compare` i') xs' + in fromTList xs'' + Disproved _ -> + withSingI st $ + case sIsJust (sTranspositions sv stl st) %~ STrue of + Proved Refl -> Tensor $ fmap (fmap (transposeMult sv stl)) ms + where + take' :: N -> Int -> [Int] + take' Z i = [i] + take' (S n) i = i : take' n (i+1) + + transposeIndices :: [Int] -> [Int] -> [Int] + transposeIndices ts' is = fmap snd $ + sortBy (\(i,_) (i',_) -> i `compare` i') $ + zip ts' is + + go :: [(N,N)] -> [Int] -> [Int] + go [] is = is + go ((s,t):ts) (i:is) = + case s' `compare` i of + EQ -> t' : go ts is + GT -> i : go ((s,t):ts) is + LT -> error $ "illegal permutation" <> show ((s,t):ts) <> "\t" <> show (i:is) + where + s' = toInt s + t' = toInt t + go _ [] = error "cannot transpose elements of empty list" + +-- |Tensor relabelling. Given a vector space and a list of relabellings, the result is a tensor +-- with the resulting rank after relabelling. Only possible if labels to be renamed are part of +-- the rank and if uniqueness of labels after relabelling is preserved. +relabel :: forall (vs :: VSpace Symbol Nat) (rl :: RelabelList Symbol) (r1 :: Rank) (r2 :: Rank) v. + (RelabelR vs rl r1 ~ 'Just r2, Sane r2 ~ 'True, SingI r1, SingI r2) => + Sing vs -> Sing rl -> Tensor r1 v -> Tensor r2 v +relabel _ _ ZeroTensor = ZeroTensor +relabel sv srl tens@(Tensor ms) = + let sr1 = sing :: Sing r1 + sr2 = sing :: Sing r2 + sh = sHead' sr1 + sr1' = sTail' sr1 + sr2' = sTail' sr2 + sr1'' = sTail sr1 + sts = sRelabelTranspositions srl (sSnd (sHead sr1)) + in case sv %~ sFst sh of + Proved Refl -> + case sSane sr1'' %~ STrue of + Proved Refl -> + case sts of + SJust sts' -> + withSingI (sFst (sHead sr1)) $ + withSingI sr1'' $ + let sn = sLengthIL (sSnd (sHead sr1)) + n = fromSing sn + ts = fromSing sts' + ts' = go ts $ take' n 0 + xs = toTListWhile tens + xs' = fmap (first (transposeIndices ts')) xs + xs'' = sortBy (\(i,_) (i',_) -> i `compare` i') xs' + in fromTList xs'' + Disproved _ -> + case sRelabelR sv srl sr1' %~ SJust sr2' of + Proved Refl -> + case sSane sr2' %~ STrue of + Proved Refl -> withSingI sr1' $ withSingI sr2' $ Tensor $ fmap (fmap (relabel sv srl)) ms + where + take' :: N -> Int -> [Int] + take' Z i = [i] + take' (S n) i = i : take' n (i+1) + + transposeIndices :: [Int] -> [Int] -> [Int] + transposeIndices ts' is = fmap snd $ + sortBy (\(i,_) (i',_) -> i `compare` i') $ + zip ts' is + + go :: [(N,N)] -> [Int] -> [Int] + go [] is = is + go ((s,t):ts) (i:is) = + case s' `compare` i of + EQ -> t' : go ts is + GT -> i : go ((s,t):ts) is + LT -> error $ "illegal permutation" <> show ((s,t):ts) <> "\t" <> show (i:is) + where + s' = toInt s + t' = toInt t + go _ [] = error "cannot transpose elements of empty list" + +-- |Get assocs list from tensor. Keys are length-indexed vectors of indices. +toList :: forall r v n. + (SingI r, SingI n, LengthR r ~ n) => + Tensor r v -> [(Vec n Int, v)] +toList ZeroTensor = [] +toList (Scalar s) = [(VNil, s)] +toList (Tensor ms) = + let st = sTail' (sing :: Sing r) + sn = sing :: Sing n + sm = sLengthR st + in case st of + SNil -> + case sn of + SS SZ -> fmap (\(i, Scalar s) -> (VCons i VNil, s)) ms + _ -> + case sn of + SS sm' -> + withSingI sm' $ + case sm %~ sm' of + Proved Refl -> + concatMap (\(i, v) -> case v of Tensor _ -> fmap (first (VCons i)) (withSingI st $ toList v)) ms + +fromList' :: forall r v n. + (Sane r ~ 'True, LengthR r ~ n) => + Sing r -> [(Vec n Int, v)] -> Tensor r v +fromList' _ [] = ZeroTensor +fromList' sr xs = + let sn = sLengthR sr + st = sTail' sr + sm = sLengthR st + in case sn of + SZ -> + case sr %~ SNil of + Proved Refl -> Scalar $ snd (head xs) + SS sm' -> + withSingI sm' $ + case sm %~ sm' of + Proved Refl -> + withSingI st $ + case sSane st %~ STrue of + Proved Refl -> + case fmap (\(i `VCons` is,v) -> (i,(is ,v))) xs of + xs' -> Tensor $ fmap (fromList' st) <$> myGroup xs' + where + myGroup :: Eq k => [(k,a)] -> [(k, [a])] + myGroup ys = + let ys' = groupBy (\(i,_) (i',_) -> i == i') ys + in fmap (\x -> (fst $ head x, fmap snd x)) ys' + +-- |Construct 'Tensor' from assocs list. Keys are length-indexed vectors of indices. +fromList :: forall r v n. + (SingI r, Sane r ~ 'True, LengthR r ~ n) => + [(Vec n Int, v)] -> Tensor r v +fromList = + let sr = sing :: Sing r + in fromList' sr + +-- |Decompose tensor into assocs list with keys being lists of indices for the first vector space +-- and values being the tensors with lower rank for the remaining vector spaces. +toTListWhile :: forall r v. + (SingI r, Sane r ~ 'True) => + Tensor r v -> [([Int], Tensor (Tail r) v)] +toTListWhile (Tensor ms) = + let sr = sing :: Sing r + st = sTail' sr + in case st %~ sTail sr of + Proved Refl -> fmap (first pure) ms + Disproved _ -> + case sSane st %~ STrue of + Proved Refl -> + case sTail sr %~ sTail st of + Proved Refl -> + withSingI st $ + withSingI (sFst (sHead st)) $ + let ms' = fmap (second toTListWhile) ms + in concatMap (\(i, xs) -> fmap (first (i :)) xs) ms' + +-- |Decompose tensor into assocs list with keys being lists of indices up to and including the +-- desired label, and values being tensors of corresponding lower rank. +toTListUntil :: forall (a :: Ix Symbol) r r' v. + (SingI r, SingI r', RemoveUntil a r ~ r', Sane r ~ 'True, Sane r' ~ 'True) => + Sing a -> Tensor r v -> [([Int], Tensor r' v)] +toTListUntil sa (Tensor ms) = + let sr = sing :: Sing r + st = sTail' sr + sh = sHead' sr + in case sSnd sh %~ sa of + Proved Refl -> withSingI st $ + case st %~ (sing :: Sing r') of + Proved Refl -> fmap (first pure) ms + Disproved _ -> + withSingI st $ + case sSane st %~ STrue of + Proved Refl -> + case sRemoveUntil sa st %~ (sing :: Sing r') of + Proved Refl -> + let ms' = fmap (second (toTListUntil sa)) ms + in concatMap (\(i, xs) -> fmap (first (i :)) xs) ms' + +-- |Construct tensor from assocs list. Keys are lists of indices, values are +-- tensors of lower rank. Used internally for tensor algebra. +fromTList :: forall r r' v.(Sane r ~ 'True, Sane r' ~ 'True, SingI r, SingI r') => + [([Int], Tensor r v)] -> Tensor r' v +fromTList [] = ZeroTensor +fromTList xs@((i0,t0):ys) + | null i0 = if null ys + then case (sing :: Sing r) %~ (sing :: Sing r') of + Proved Refl -> t0 + else error $ "illegal assocs in fromTList : " ++ show (fmap fst xs) + | otherwise = + let sr' = sing :: Sing r' + st' = sTail' sr' + in withSingI st' $ + case sSane st' of + STrue -> Tensor $ fmap (fmap fromTList) xs''' + where + xs' = fmap (\(i:is,v) -> (i,(is,v))) xs + xs'' = groupBy (\(i,_) (i',_) -> i == i') xs' + xs''' = fmap (\x -> (fst $ head x, map snd x)) xs'' diff --git a/src/Math/Tensor/Safe/Proofs.hs b/src/Math/Tensor/Safe/Proofs.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/Safe/Proofs.hs @@ -0,0 +1,267 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.Safe.Proofs +Description : Identities for functions on generalized tensor ranks. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Identities for functions on generalized tensor ranks. +-} +----------------------------------------------------------------------------- +module Math.Tensor.Safe.Proofs + ( -- * Tails of sane ranks are sane + saneTail'Proof + , singITail'Proof + , -- * Properties of merged ranks + saneMergeRProof + , proofMergeLT + , proofMergeGT + , proofMergeIxNotEQ + , proofMergeIxLT + , proofMergeIxGT + , -- * Properties of contractions + saneContractProof + , singletonContractProof + , contractTailDiffVProof + , contractTailSameVNoConProof + , contractTailSameVNoCovProof + , contractTailSameVDiffIProof + , contractTailSameVSameIProof + ) where + +import Math.Tensor.Safe.TH + +import Data.Constraint + ( Dict (Dict) + , (:-) (Sub) + ) +import Unsafe.Coerce (unsafeCoerce) + +import Data.Singletons.Prelude + ( Sing, SingI + , PEq ((==)) + , Symbol + , Fst, Snd, Compare + ) + +-- |The 'Tail'' of a sane rank type is sane. +{-# INLINE saneTail'Proof #-} +saneTail'Proof :: forall (r :: Rank).Sing r -> (Sane r ~ 'True) :- (Sane (Tail' r) ~ 'True) +saneTail'Proof _ = Sub axiom + where + axiom :: Sane r ~ 'True => Dict (Sane (Tail' r) ~ 'True) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |If a rank type has a 'SingI' instance, the tail has a 'SingI' instance. +{-# INLINE singITail'Proof #-} +singITail'Proof :: forall (r :: Rank).Sing r -> SingI r :- SingI (Tail' r) +singITail'Proof _ = Sub axiom + where + axiom :: SingI r => Dict (SingI (Tail' r)) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |Successfully merging two sane rank types (result is not @Nothing@) yields a sane rank type. +{-# INLINE saneMergeRProof #-} +saneMergeRProof :: forall (r :: Rank) (r' :: Rank) (r'' :: Rank). + Sing r -> Sing r' -> + (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'') :- (Sane r'' ~ 'True) +saneMergeRProof _ _ = Sub axiom + where + axiom :: (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'') => + Dict (Sane r'' ~ 'True) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |If two rank types can be merged and the first 'VSpace' of the first rank type is less than +-- the first 'VSpace' of the second rank type, the 'Tail'' of the merged rank type is equal to +-- the tail of the first rank type merged with the second rank type. +{-# INLINE proofMergeLT #-} +proofMergeLT :: forall (r :: Rank) (r' :: Rank) (r'' :: Rank). + Sing r -> Sing r' -> + (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'', + Compare (Fst (Head' r)) (Fst (Head' r')) ~ 'LT) + :- (MergeR (Tail' r) r' ~ 'Just (Tail' r'')) +proofMergeLT _ _ = Sub axiom + where + axiom :: (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'', + Compare (Fst (Head' r)) (Fst (Head' r')) ~ 'LT) + => Dict (MergeR (Tail' r) r' ~ 'Just (Tail' r'')) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |If two rank types can be merged and the first 'VSpace' of the first rank type coincides with +-- the first 'VSpace' of the second rank type, the first index of the first rank type cannot +-- equal the first index of the second rank type. +{-# INLINE proofMergeIxNotEQ #-} +proofMergeIxNotEQ :: forall (r :: Rank) (r' :: Rank) (r'' :: Rank). + Sing r -> Sing r' -> + (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'', + Compare (Fst (Head' r)) (Fst (Head' r')) ~ 'EQ) + :- ((IxCompare (Snd (Head' r)) (Snd (Head' r')) == 'EQ) ~ 'False) +proofMergeIxNotEQ _ _ = Sub axiom + where + axiom :: (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'', + Compare (Fst (Head' r)) (Fst (Head' r')) ~ 'EQ) + => Dict ((IxCompare (Snd (Head' r)) (Snd (Head' r')) == 'EQ) ~ 'False) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |If two rank types can be merged and the first 'VSpace' of the first rank type coincides with +-- the first 'VSpace' of the second rank type and the first index of the first rank type compares +-- less than the first index of the second rank type, the 'Tail'' of the merged rank type is equal +-- to the tail of the first rank type merged with the second rank type. +{-# INLINE proofMergeIxLT #-} +proofMergeIxLT :: forall (r :: Rank) (r' :: Rank) (r'' :: Rank). + Sing r -> Sing r' -> + (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'', + Compare (Fst (Head' r)) (Fst (Head' r')) ~ 'EQ, + IxCompare (Snd (Head' r)) (Snd (Head' r')) ~ 'LT) + :- (MergeR (Tail' r) r' ~ 'Just (Tail' r'')) +proofMergeIxLT _ _ = Sub axiom + where + axiom :: (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'', + Compare (Fst (Head' r)) (Fst (Head' r')) ~ 'EQ, + IxCompare (Snd (Head' r)) (Snd (Head' r')) ~ 'LT) + => Dict (MergeR (Tail' r) r' ~ 'Just (Tail' r'')) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |If two rank types can be merged and the first 'VSpace' of the first rank type is greater than +-- the first 'VSpace' of the second rank type, the 'Tail'' of the merged rank type is equal to +-- the first rank type merged with the tail of the second rank type. +{-# INLINE proofMergeGT #-} +proofMergeGT :: forall (r :: Rank) (r' :: Rank) (r'' :: Rank). + Sing r -> Sing r' -> + (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'', + Compare (Fst (Head' r)) (Fst (Head' r')) ~ 'GT) + :- (MergeR r (Tail' r') ~ 'Just (Tail' r'')) +proofMergeGT _ _ = Sub axiom + where + axiom :: (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'', + Compare (Fst (Head' r)) (Fst (Head' r')) ~ 'GT) + => Dict (MergeR r (Tail' r') ~ 'Just (Tail' r'')) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |If two rank types can be merged and the first 'VSpace' of the first rank type coincides with +-- the first 'VSpace' of the second rank type and the first index of the first rank type compares +-- greater than the first index of the second rank type, the 'Tail'' of the merged rank type is equal +-- to the first rank type merged with the tail of the second rank type. +{-# INLINE proofMergeIxGT #-} +proofMergeIxGT :: forall (r :: Rank) (r' :: Rank) (r'' :: Rank). + Sing r -> Sing r' -> + (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'', + Compare (Fst (Head' r)) (Fst (Head' r')) ~ 'EQ, + IxCompare (Snd (Head' r)) (Snd (Head' r')) ~ 'GT) + :- (MergeR r (Tail' r') ~ 'Just (Tail' r'')) +proofMergeIxGT _ _ = Sub axiom + where + axiom :: (Sane r ~ 'True, Sane r' ~ 'True, MergeR r r' ~ 'Just r'', + Compare (Fst (Head' r)) (Fst (Head' r')) ~ 'EQ, + IxCompare (Snd (Head' r)) (Snd (Head' r')) ~ 'GT) + => Dict (MergeR r (Tail' r') ~ 'Just (Tail' r'')) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |If a rank type is sane, its contraction is also sane. +{-# INLINE saneContractProof #-} +saneContractProof :: forall (r :: Rank).Sing r -> (Sane r ~ 'True) :- (Sane (ContractR r) ~ 'True) +saneContractProof _ = Sub axiom + where + axiom :: Sane r ~ 'True => Dict (Sane (ContractR r) ~ 'True) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |The contraction of the empty rank type is the empty rank type. +{-# INLINE singletonContractProof #-} +singletonContractProof :: forall (r :: Rank). + Sing r -> (Tail' r ~ '[]) :- (ContractR r ~ r) +singletonContractProof _ = Sub axiom + where + axiom :: Tail' r ~ '[] => Dict (ContractR r ~ r) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |If the first two labels of a rank type cannot be contracted because they belong to +-- different 'VSpace's, the 'Tail'' of the contracted rank type is equal to the contraction +-- of the 'Tail'' of the rank type. +{-# INLINE contractTailDiffVProof #-} +contractTailDiffVProof :: forall (r :: Rank) (t :: Rank) (t' :: Rank). + Sing r -> + (t ~ Tail' r, t' ~ Tail' t, (Fst (Head' r) == Fst (Head' t)) ~ 'False) + :- (Tail' (ContractR r) ~ ContractR t) +contractTailDiffVProof _ = Sub axiom + where + axiom :: (t ~ Tail' r, t' ~ Tail' t, (Fst (Head' r) == Fst (Head' t)) ~ 'False) + => Dict (Tail' (ContractR r) ~ ContractR t) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + + +-- |If the first two labels of a rank type cannot be contracted because the first label is +-- covariant, the 'Tail'' of the contracted rank type is equal to the contraction +-- of the 'Tail'' of the rank type. +{-# INLINE contractTailSameVNoConProof #-} +contractTailSameVNoConProof :: forall (r :: Rank) (t :: Rank) (t' :: Rank) (i :: Symbol). + Sing r -> + (t ~ Tail' r, t' ~ Tail' t, (Fst (Head' r) == Fst (Head' t)) ~ 'True, + Snd (Head' r) ~ 'ICov i) + :- (Tail' (ContractR r) ~ ContractR t) +contractTailSameVNoConProof _ = Sub axiom + where + axiom :: (t ~ Tail' r, t' ~ Tail' t, (Fst (Head' r) == Fst (Head' t)) ~ 'True, + Snd (Head' r) ~ 'ICov i) + => Dict (Tail' (ContractR r) ~ ContractR t) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |If the first two labels of a rank type cannot be contracted because the second label is +-- covariant, the 'Tail'' of the contracted rank type is equal to the contraction +-- of the 'Tail'' of the rank type. +{-# INLINE contractTailSameVNoCovProof #-} +contractTailSameVNoCovProof :: forall (r :: Rank) (t :: Rank) (t' :: Rank) (i :: Symbol). + Sing r -> + (t ~ Tail' r, t' ~ Tail' t, (Fst (Head' r) == Fst (Head' t)) ~ 'True, + Snd (Head' t) ~ 'ICon i) + :- (Tail' (ContractR r) ~ ContractR t) +contractTailSameVNoCovProof _ = Sub axiom + where + axiom :: (t ~ Tail' r, t' ~ Tail' t, (Fst (Head' r) == Fst (Head' t)) ~ 'True, + Snd (Head' t) ~ 'ICon i) + => Dict (Tail' (ContractR r) ~ ContractR t) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |If the first two labels of a rank type cannot be contracted because they differ, +-- the 'Tail'' of the contracted rank type is equal to the contraction of the 'Tail'' of the rank type. +{-# INLINE contractTailSameVDiffIProof #-} +contractTailSameVDiffIProof :: forall (r :: Rank) (t :: Rank) (t' :: Rank) (i :: Symbol) (j :: Symbol). + Sing r -> + (t ~ Tail' r, t' ~ Tail' t, (Fst (Head' r) == Fst (Head' t)) ~ 'True, + Snd (Head' r) ~ 'ICon i, Snd (Head' t) ~ 'ICov j, (i == j) ~ 'False) + :- (Tail' (ContractR r) ~ ContractR t) +contractTailSameVDiffIProof _ = Sub axiom + where + axiom :: (t ~ Tail' r, t' ~ Tail' t, (Fst (Head' r) == Fst (Head' t)) ~ 'True, + Snd (Head' r) ~ 'ICon i, Snd (Head' t) ~ 'ICov j, (i == j) ~ 'False) + => Dict (Tail' (ContractR r) ~ ContractR t) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) + +-- |If the first two labels of a rank type can be contracted, the contracted rank type is equal +-- to the contraction of the tail. +{-# INLINE contractTailSameVSameIProof #-} +contractTailSameVSameIProof :: forall (r :: Rank) (t :: Rank) (t' :: Rank) (i :: Symbol) (j :: Symbol). + Sing r -> + (t ~ Tail' r, t' ~ Tail' t, (Fst (Head' r) == Fst (Head' t)) ~ 'True, + Snd (Head' r) ~ 'ICon i, Snd (Head' t) ~ 'ICov j, (i == j) ~ 'True) + :- (ContractR t' ~ ContractR r) +contractTailSameVSameIProof _ = Sub axiom + where + axiom :: (t ~ Tail' r, t' ~ Tail' t, (Fst (Head' r) == Fst (Head' t)) ~ 'True, + Snd (Head' r) ~ 'ICon i, Snd (Head' t) ~ 'ICov j, (i == j) ~ 'True) + => Dict (ContractR t' ~ ContractR r) + axiom = unsafeCoerce (Dict :: Dict (a ~ a)) diff --git a/src/Math/Tensor/Safe/TH.hs b/src/Math/Tensor/Safe/TH.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/Safe/TH.hs @@ -0,0 +1,568 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE DefaultSignatures #-} +{-# LANGUAGE EmptyCase #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE InstanceSigs #-} +{-# LANGUAGE NoStarIsType #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE OverloadedStrings #-} + +{-# LANGUAGE CPP #-} +#if MIN_VERSION_base(4,14,0) +{-# LANGUAGE StandaloneKindSignatures #-} +#endif + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.Safe.TH +Description : Generalized rank lifted to the type level. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Generalized rank lifted to the type level. For documentation see re-exports +in "Math.Tensor.Safe". +-} +----------------------------------------------------------------------------- +module Math.Tensor.Safe.TH where + +import Data.Singletons.Prelude +import Data.Singletons.Prelude.Enum +import Data.Singletons.Prelude.List.NonEmpty hiding (sLength) +import Data.Singletons.Prelude.Ord +import Data.Singletons.TH +import Data.Singletons.TypeLits + +import Data.List.NonEmpty (NonEmpty((:|)),sort,sortBy,(<|)) + +$(singletons [d| + data N where + Z :: N + S :: N -> N + + fromNat :: Nat -> N + fromNat n = case n == 0 of + True -> Z + False -> S $ fromNat (pred n) + + deriving instance Show N + deriving instance Eq N + instance Ord N where + Z <= _ = True + (S _) <= Z = False + (S n) <= (S m) = n <= m + + instance Num N where + Z + n = n + (S n) + m = S $ n + m + + n - Z = n + Z - S _ = error "cannot subtract (S n) from Z!" + (S n) - (S m) = n - m + + negate Z = Z + negate _ = error "cannot negate (S n)!" + + Z * _ = Z + (S n) * m = m + n * m + + abs n = n + signum n = n + + fromInteger n + | n == 0 = Z + | otherwise = S $ fromInteger (n-1) + + data VSpace a b = VSpace { vId :: a, + vDim :: b } + deriving (Show, Ord, Eq) + + data Ix a = ICon a | ICov a deriving (Show, Ord, Eq) + + ixCompare :: Ord a => Ix a -> Ix a -> Ordering + ixCompare (ICon a) (ICon b) = compare a b + ixCompare (ICon a) (ICov b) = case compare a b of + LT -> LT + EQ -> LT + GT -> GT + ixCompare (ICov a) (ICon b) = case compare a b of + LT -> LT + EQ -> GT + GT -> GT + ixCompare (ICov a) (ICov b) = compare a b + + data IList a = ConCov (NonEmpty a) (NonEmpty a) | + Cov (NonEmpty a) | + Con (NonEmpty a) + deriving (Show, Ord, Eq) + + type GRank s n = [(VSpace s n, IList s)] + type Rank = GRank Symbol Nat + + isAscending :: Ord a => [a] -> Bool + isAscending [] = True + isAscending [_] = True + isAscending (x:y:xs) = x < y && isAscending (y:xs) + + isAscendingNE :: Ord a => NonEmpty a -> Bool + isAscendingNE (x :| xs) = isAscending (x:xs) + + isAscendingI :: Ord a => IList a -> Bool + isAscendingI (ConCov x y) = isAscendingNE x && isAscendingNE y + isAscendingI (Con x) = isAscendingNE x + isAscendingI (Cov y) = isAscendingNE y + + isLengthNE :: NonEmpty a -> Nat -> Bool + isLengthNE (_ :| []) l = l == 1 + isLengthNE (_ :| (x:xs)) l = isLengthNE (x :| xs) (pred l) + + lengthNE :: NonEmpty a -> N + lengthNE (_ :| []) = S Z + lengthNE (_ :| (x:xs)) = S Z + lengthNE (x :| xs) + + lengthIL :: IList a -> N + lengthIL (ConCov xs ys) = lengthNE xs + lengthNE ys + lengthIL (Con xs) = lengthNE xs + lengthIL (Cov ys) = lengthNE ys + + lengthR :: GRank s n -> N + lengthR [] = Z + lengthR ((_,x):xs) = lengthIL x + lengthR xs + + sane :: (Ord a, Ord b) => [(VSpace a b, IList a)] -> Bool + sane [] = True + sane [(_, is)] = isAscendingI is + sane ((v, is):(v', is'):xs) = + v < v' && isAscendingI is && sane ((v',is'):xs) + + head' :: Ord s => GRank s n -> (VSpace s n, Ix s) + head' ((v, l):_) = (v, case l of + ConCov (a :| _) (b :| _) -> + case compare a b of + LT -> ICon a + EQ -> ICon a + GT -> ICov b + Con (a :| _) -> ICon a + Cov (a :| _) -> ICov a) + head' [] = error "head' of empty list" + + tail' :: Ord s => GRank s n -> GRank s n + tail' ((v, l):ls) = + let l' = case l of + ConCov (a :| []) (b :| []) -> + case compare a b of + LT -> Just $ Cov (b :| []) + EQ -> Just $ Cov (b :| []) + GT -> Just $ Con (a :| []) + + ConCov (a :| (a':as)) (b :| []) -> + case compare a b of + LT -> Just $ ConCov (a' :| as) (b :| []) + EQ -> Just $ ConCov (a' :| as) (b :| []) + GT -> Just $ Con (a :| (a':as)) + + ConCov (a :| []) (b :| (b':bs)) -> + case compare a b of + LT -> Just $ Cov (b :| (b':bs)) + EQ -> Just $ Cov (b :| (b':bs)) + GT -> Just $ ConCov (a :| []) (b' :| bs) + + ConCov (a :| (a':as)) (b :| (b':bs)) -> + case compare a b of + LT -> Just $ ConCov (a' :| as) (b :| (b':bs)) + EQ -> Just $ ConCov (a' :| as) (b :| (b':bs)) + GT -> Just $ ConCov (a :| (a':as)) (b' :| bs) + + Con (_ :| []) -> Nothing + Con (_ :| (a':as)) -> Just $ Con (a' :| as) + + Cov (_ :| []) -> Nothing + Cov (_ :| (a':as)) -> Just $ Cov (a' :| as) + in case l' of + Just l'' -> (v, l''):ls + Nothing -> ls + tail' [] = error "tail' of empty list" + + mergeR :: (Ord s, Ord n) => GRank s n -> GRank s n -> Maybe (GRank s n) + mergeR [] ys = Just ys + mergeR xs [] = Just xs + mergeR ((xv,xl):xs) ((yv,yl):ys) = + case compare xv yv of + LT -> ((xv,xl) :) <$> mergeR xs ((yv,yl):ys) + EQ -> do + xl' <- mergeIL xl yl + xs' <- mergeR xs ys + Just $ (xv, xl') : xs' + GT -> ((yv,yl) :) <$> mergeR ((xv,xl):xs) ys + + mergeIL :: Ord a => IList a -> IList a -> Maybe (IList a) + mergeIL (ConCov xs ys) (ConCov xs' ys') = do + xs'' <- mergeNE xs xs' + ys'' <- mergeNE ys ys' + Just $ ConCov xs'' ys'' + mergeIL (ConCov xs ys) (Con xs') = do + xs'' <- mergeNE xs xs' + Just $ ConCov xs'' ys + mergeIL (ConCov xs ys) (Cov ys') = do + ys'' <- mergeNE ys ys' + Just $ ConCov xs ys'' + mergeIL (Con xs) (ConCov xs' ys) = do + xs'' <- mergeNE xs xs' + Just $ ConCov xs'' ys + mergeIL (Con xs) (Con xs') = Con <$> mergeNE xs xs' + mergeIL (Con xs) (Cov ys) = Just $ ConCov xs ys + mergeIL (Cov ys) (ConCov xs ys') = do + ys'' <- mergeNE ys ys' + Just $ ConCov xs ys'' + mergeIL (Cov ys) (Con xs) = Just $ ConCov xs ys + mergeIL (Cov ys) (Cov ys') = Cov <$> mergeNE ys ys' + + merge :: Ord a => [a] -> [a] -> Maybe [a] + merge [] ys = Just ys + merge xs [] = Just xs + merge (x:xs) (y:ys) = + case compare x y of + LT -> (x :) <$> merge xs (y:ys) + EQ -> Nothing + GT -> (y :) <$> merge (x:xs) ys + + mergeNE :: Ord a => NonEmpty a -> NonEmpty a -> Maybe (NonEmpty a) + mergeNE (x :| xs) (y :| ys) = + case compare x y of + LT -> (x :|) <$> merge xs (y:ys) + EQ -> Nothing + GT -> (y :|) <$> merge (x:xs) ys + + contractR :: Ord s => GRank s n -> GRank s n + contractR [] = [] + contractR ((v, is):xs) = case contractI is of + Nothing -> contractR xs + Just is' -> (v, is') : contractR xs + + prepICon :: a -> IList a -> IList a + prepICon a (ConCov (x:|xs) y) = ConCov (a:|(x:xs)) y + prepICon a (Con (x:|xs)) = Con (a:|(x:xs)) + prepICon a (Cov y) = ConCov (a:|[]) y + + prepICov :: a -> IList a -> IList a + prepICov a (ConCov x (y:|ys)) = ConCov x (a:|(y:ys)) + prepICov a (Con x) = ConCov x (a:|[]) + prepICov a (Cov (y:|ys)) = Cov (a:|(y:ys)) + + contractI :: Ord a => IList a -> Maybe (IList a) + contractI (ConCov (x:|xs) (y:|ys)) = + case compare x y of + EQ -> case xs of + [] -> case ys of + [] -> Nothing + (y':ys') -> Just $ Cov (y' :| ys') + (x':xs') -> case ys of + [] -> Just $ Con (x' :| xs') + (y':ys') -> contractI $ ConCov (x':|xs') (y':|ys') + LT -> case xs of + [] -> Just $ ConCov (x:|xs) (y:|ys) + (x':xs') -> case contractI $ ConCov (x':|xs') (y:|ys) of + Nothing -> Just $ Con (x:|[]) + Just i -> Just $ prepICon x i + GT -> case ys of + [] -> Just $ ConCov (x:|xs) (y:|ys) + (y':ys') -> case contractI $ ConCov (x:|xs) (y':|ys') of + Nothing -> Just $ Cov (y:|[]) + Just i -> Just $ prepICov y i + contractI (Con x) = Just $ Con x + contractI (Cov x) = Just $ Cov x + + subsetNE :: Ord a => NonEmpty a -> NonEmpty a -> Bool + subsetNE (x :| []) ys = x `elemNE` ys + subsetNE (x :| (x':xs)) ys = (x `elemNE` ys) && ((x' :| xs) `subsetNE` ys) + + elemNE :: Ord a => a -> NonEmpty a -> Bool + elemNE a (x :| []) = a == x + elemNE a (x :| (x':xs)) = case compare a x of + LT -> False + EQ -> True + GT -> elemNE a (x' :| xs) + + canTransposeCon :: (Ord s, Ord n) => VSpace s n -> s -> s -> GRank s n -> Bool + canTransposeCon _ _ _ [] = False + canTransposeCon v a b ((v',il):r) = + case compare v v' of + LT -> False + GT -> canTransposeCon v a b r + EQ -> case il of + Cov _ -> canTransposeCon v a b r + Con cs -> case elemNE a cs of + True -> case elemNE b cs of + True -> True + False -> False + False -> case elemNE b cs of + True -> False + False -> canTransposeCon v a b r + ConCov cs _ -> case elemNE a cs of + True -> case elemNE b cs of + True -> True + False -> False + False -> case elemNE b cs of + True -> False + False -> canTransposeCon v a b r + + canTransposeCov :: (Ord s, Ord n) => VSpace s n -> s -> s -> GRank s n -> Bool + canTransposeCov _ _ _ [] = False + canTransposeCov v a b ((v',il):r) = + case compare v v' of + LT -> False + GT -> canTransposeCov v a b r + EQ -> case il of + Con _ -> canTransposeCov v a b r + Cov cs -> case elemNE a cs of + True -> case elemNE b cs of + True -> True + False -> False + False -> case elemNE b cs of + True -> False + False -> canTransposeCov v a b r + ConCov _ cs -> case elemNE a cs of + True -> case elemNE b cs of + True -> True + False -> False + False -> case elemNE b cs of + True -> False + False -> canTransposeCov v a b r + + canTranspose :: (Ord s, Ord n) => VSpace s n -> Ix s -> Ix s -> GRank s n -> Bool + canTranspose v (ICon a) (ICon b) r = case compare a b of + LT -> canTransposeCon v a b r + EQ -> True + GT -> canTransposeCov v b a r + canTranspose v (ICov a) (ICov b) r = case compare a b of + LT -> canTransposeCov v a b r + EQ -> True + GT -> canTransposeCov v b a r + canTranspose _ (ICov _) (ICon _) _ = False + canTranspose _ (ICon _) (ICov _) _ = False + + removeUntil :: Ord s => Ix s -> GRank s n -> GRank s n + removeUntil i r = go i r + where + go i' r' + | snd (head' r') == i' = tail' r' + | otherwise = go i $ tail' r' + + data TransList a = TransCon (NonEmpty a) (NonEmpty a) | + TransCov (NonEmpty a) (NonEmpty a) + deriving (Show, Eq) + + saneTransList :: Ord a => TransList a -> Bool + saneTransList tl = + case tl of + TransCon sources targets -> isAscendingNE sources && + sort targets == sources + TransCov sources targets -> isAscendingNE sources && + sort targets == sources + + canTransposeMult :: (Ord s, Ord n) => VSpace s n -> TransList s -> GRank s n -> Bool + canTransposeMult vs tl r = case transpositions vs tl r of + Just _ -> True + Nothing -> False + + transpositions :: (Ord s, Ord n) => VSpace s n -> TransList s -> GRank s n -> Maybe [(N,N)] + transpositions _ _ [] = Nothing + transpositions vs tl ((vs',il):r) = + case saneTransList tl of + False -> Nothing + True -> + case compare vs vs' of + LT -> Nothing + GT -> transpositions vs tl r + EQ -> + case il of + Con xs -> + case tl of + TransCon sources targets -> transpositions' sources targets (fmap Just xs) + TransCov _ _ -> Nothing + Cov xs -> + case tl of + TransCov sources targets -> transpositions' sources targets (fmap Just xs) + TransCon _ _ -> Nothing + ConCov xsCon xsCov -> + case tl of + TransCon sources targets -> transpositions' sources targets (xsCon `zipCon` xsCov) + TransCov sources targets -> transpositions' sources targets (xsCon `zipCov` xsCov) + + zipCon :: Ord a => NonEmpty a -> NonEmpty a -> NonEmpty (Maybe a) + zipCon (x :| xs) (y :| ys) = + case ICon x `ixCompare` ICov y of + LT -> case xs of + [] -> Just x :| fmap (const Nothing) (y:ys) + (x':xs') -> Just x <| zipCon (x' :| xs') (y :| ys) + GT -> case ys of + [] -> Nothing :| fmap Just (x : xs) + (y':ys') -> Nothing <| zipCon (x :| xs) (y' :| ys') + + zipCov :: Ord a => NonEmpty a -> NonEmpty a -> NonEmpty (Maybe a) + zipCov (x :| xs) (y :| ys) = + case ICon x `ixCompare` ICov y of + LT -> case xs of + [] -> Nothing :| fmap Just (y:ys) + (x':xs') -> Nothing <| zipCov (x' :| xs') (y :| ys) + GT -> case ys of + [] -> Just y :| fmap (const Nothing) (x : xs) + (y':ys') -> Just y <| zipCov (x :| xs) (y' :| ys') + + transpositions' :: Eq a => NonEmpty a -> NonEmpty a -> NonEmpty (Maybe a) -> Maybe [(N,N)] + transpositions' sources targets xs = + do + ss <- mapM (`find` xs') sources + ts <- mapM (`find` xs') targets + zip' ss ts + where + xs' = go' Z xs + + go' :: N -> NonEmpty a -> NonEmpty (N,a) + go' n (y :| []) = (n,y) :| [] + go' n (y :| (y':ys')) = (n,y) <| go' (S n) (y' :| ys') + + find :: Eq a => a -> NonEmpty (N, Maybe a) -> Maybe N + find _ ((_,Nothing) :| []) = Nothing + find a ((_,Nothing) :| (y':ys')) = find a (y' :| ys') + find a ((n,Just y) :| ys) = + case a == y of + True -> Just n + False -> + case ys of + [] -> Nothing + y':ys' -> find a (y' :| ys') + + zip' :: NonEmpty a -> NonEmpty b -> Maybe [(a,b)] + zip' (a :| []) (b :| []) = Just [(a,b)] + zip' (_ :| (_:_)) (_ :| []) = Nothing + zip' (_ :| []) (_ :| (_:_)) = Nothing + zip' (y:|(y':ys')) (z:|(z':zs')) = ((y,z):) <$> zip' (y':|ys') (z':|zs') + + type RelabelList s = NonEmpty (s,s) + + saneRelabelList :: Ord a => NonEmpty (a,a) -> Bool + saneRelabelList xs = isAscendingNE xs && + isAscendingNE xs' + where + xs' = sort $ fmap (\(a,b) -> (b,a)) xs + + relabelNE :: Ord a => NonEmpty (a,a) -> NonEmpty a -> Maybe (NonEmpty (a,a)) + relabelNE = go + where + go :: Ord a => NonEmpty (a,a) -> NonEmpty a -> Maybe (NonEmpty (a,a)) + go ((source,target) :| ms) (x :| xs) = + case source `compare` x of + LT -> + case ms of + [] -> Just $ (\a -> (a,a)) <$> (x :| xs) + m':ms' -> go (m' :| ms') (x :| xs) + EQ -> + case ms of + [] -> Just $ ((target,source) :|) $ fmap (\a -> (a,a)) xs + m':ms' -> + case xs of + [] -> Just $ (target,source) :| [] + x':xs' -> ((target,source) <|) <$> go (m' :| ms') (x' :| xs') + GT -> + case xs of + [] -> Just $ (x,x) :| [] + x':xs' -> ((x,x) <|) <$> go ((source,target) :| ms) (x' :| xs') + + relabelR :: (Ord s, Ord n) => VSpace s n -> RelabelList s -> GRank s n -> Maybe (GRank s n) + relabelR _ _ [] = Nothing + relabelR vs rls ((vs',il):r) = + case vs `compare` vs' of + LT -> Nothing + EQ -> (\il' -> (vs',il'):r) <$> relabelIL rls il + GT -> ((vs',il) :) <$> relabelR vs rls r + + relabelIL :: Ord a => NonEmpty (a,a) -> IList a -> Maybe (IList a) + relabelIL rl is = case relabelIL' rl is of + Nothing -> Nothing + Just (Con is') -> Just $ Con $ fst <$> is' + Just (Cov is') -> Just $ Cov $ fst <$> is' + Just (ConCov is' js') -> Just $ ConCov (fst <$> is') (fst <$> js') + + relabelIL' :: Ord a => NonEmpty (a,a) -> IList a -> Maybe (IList (a,a)) + relabelIL' rl (Con is) = + do + is' <- Con . sort <$> relabelNE rl is + case isAscendingI is' of + True -> return is' + False -> Nothing + relabelIL' rl (Cov is) = + do + is' <- Cov . sort <$> relabelNE rl is + case isAscendingI is' of + True -> return is' + False -> Nothing + relabelIL' rl (ConCov is js) = + do + is' <- sort <$> relabelNE rl is + js' <- sort <$> relabelNE rl js + let l' = ConCov is' js' + case isAscendingI l' of + True -> return l' + False -> Nothing + + relabelTranspositions :: Ord a => NonEmpty (a,a) -> IList a -> Maybe [(N,N)] + relabelTranspositions rl is = + case relabelIL' rl is of + Nothing -> Nothing + Just (Con is') -> Just $ relabelTranspositions' is' + Just (Cov is') -> Just $ relabelTranspositions' is' + Just (ConCov is' js') -> Just $ relabelTranspositions' $ is' `zipConCov` js' + + zipConCov :: Ord a => NonEmpty a -> NonEmpty a -> NonEmpty a + zipConCov = go + where + go :: Ord a => NonEmpty a -> NonEmpty a -> NonEmpty a + go (i :| is) (j :| js) = + case i `compare` j of + LT -> case is of + [] -> i <| (j:|js) + i':is' -> i <| go (i' :| is') (j :| js) + EQ -> case is of + [] -> i <| (j:|js) + i':is' -> i <| go (i' :| is') (j :| js) + GT -> case js of + [] -> j <| (i:|is) + j':js' -> j <| go (i :| is) (j' :| js') + + relabelTranspositions' :: Ord a => NonEmpty (a,a) -> [(N,N)] + relabelTranspositions' is = go'' is''' + where + is' = go Z is + is'' = sortBy (\a b -> snd a `compare` snd b) is' + is''' = go' Z is'' + --is'''' = sort is''' + + go :: N -> NonEmpty (a,b) -> NonEmpty (N,b) + go n ((_,y) :| []) = (n,y) :| [] + go n ((_,y) :| (j:js)) = (n,y) <| go (S n) (j :| js) + + go' :: N -> NonEmpty (a,b) -> NonEmpty (a,N) + go' n ((x,_) :| []) = (x,n) :| [] + go' n ((x,_) :| (j:js)) = (x,n) <| go' (S n) (j :| js) + + go'' :: NonEmpty (a,a) -> [(a,a)] + go'' ((x1,x2) :| []) = [(x2,x1)] + go'' ((x1,x2) :| (y:ys)) = (x2,x1) : go'' (y :| ys) + |]) + +toInt :: N -> Int +toInt Z = 0 +toInt (S n) = 1 + toInt n diff --git a/src/Math/Tensor/Safe/Vector.hs b/src/Math/Tensor/Safe/Vector.hs new file mode 100644 --- /dev/null +++ b/src/Math/Tensor/Safe/Vector.hs @@ -0,0 +1,53 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE StandaloneDeriving #-} + +----------------------------------------------------------------------------- +{-| +Module : Math.Tensor.Safe.Vector +Description : Length-typed vector. +Copyright : (c) Nils Alex, 2020 +License : MIT +Maintainer : nils.alex@fau.de +Stability : experimental + +Length-typed vector. +-} +----------------------------------------------------------------------------- +module Math.Tensor.Safe.Vector + ( Vec(..) + , vecFromListUnsafe + ) where + +import Math.Tensor.Safe.TH + +import Data.Kind (Type) +import Data.Singletons (Sing) + +data Vec :: N -> Type -> Type where + VNil :: Vec 'Z a + VCons :: a -> Vec n a -> Vec ('S n) a + +deriving instance Show a => Show (Vec n a) + +instance Eq a => Eq (Vec n a) where + VNil == VNil = True + (x `VCons` xs) == (y `VCons` ys) = x == y && xs == ys + +instance Ord a => Ord (Vec n a) where + VNil `compare` VNil = EQ + (x `VCons` xs) `compare` (y `VCons` ys) = + case x `compare` y of + LT -> LT + EQ -> xs `compare` ys + GT -> GT + +vecFromListUnsafe :: forall (n :: N) a. + Sing n -> [a] -> Vec n a +vecFromListUnsafe SZ [] = VNil +vecFromListUnsafe (SS sn) (x:xs) = + let xs' = vecFromListUnsafe sn xs + in x `VCons` xs' +vecFromListUnsafe _ _ = error "cannot reconstruct vector from list: incompatible lengths"