diff --git a/Math/LinearMap/Asserted.hs b/Math/LinearMap/Asserted.hs --- a/Math/LinearMap/Asserted.hs +++ b/Math/LinearMap/Asserted.hs @@ -95,10 +95,12 @@ μ *^ LinearFunction f = LinearFunction $ (μ*^) . f instance VectorSpace w => Semimanifold (LinearFunction s v w) where type Needle (LinearFunction s v w) = LinearFunction s v w +#if !MIN_VERSION_manifolds_core(0,6,0) toInterior = pure fromInterior = id - (.+~^) = (^+^) translateP = Tagged (^+^) +#endif + (.+~^) = (^+^) instance VectorSpace w => PseudoAffine (LinearFunction s v w) where f.-~.g = return $ f^-^g diff --git a/Math/LinearMap/Category.hs b/Math/LinearMap/Category.hs --- a/Math/LinearMap/Category.hs +++ b/Math/LinearMap/Category.hs @@ -28,7 +28,7 @@ -- ** Function implementation LinearFunction (..), type (-+>)(), Bilinear - , lfun + , lfun, (-+$>) -- ** Tensor implementation , LinearMap (..), type (+>)() , (⊕), (>+<) @@ -69,6 +69,11 @@ , finishEigenSystem , Eigenvector(..) -- * The classes of suitable vector spaces + -- $classesExplanation + -- ** Deriving from basis + , module Math.LinearMap.Category.Instances.Deriving + -- ** The classes + , module Data.VectorSpace , LSpace , TensorSpace (..) , LinearSpace (..) @@ -83,6 +88,9 @@ , (.⊗) -- ** Hilbert space operations , (·), DualSpace, riesz, coRiesz, showsPrecAsRiesz, (.<) + -- ** Standard decompositions + , TensorDecomposable(..), RieszDecomposable(..) + , tensorDecomposeShowsPrec, rieszDecomposeShowsPrec -- ** Constraint synonyms , HilbertSpace, SimpleSpace, RealSpace , Num'(..) @@ -91,6 +99,7 @@ -- ** Double-dual, scalar-scalar etc. identity , ClosedScalarWitness(..), TrivialTensorWitness(..) , ScalarSpaceWitness(..), DualSpaceWitness(..), LinearManifoldWitness(..) + , DualFinitenessWitness(..) -- ** Misc , relaxNorm, transformNorm, transformVariance , findNormalLength, normalLength @@ -103,6 +112,7 @@ import Math.LinearMap.Category.Class import Math.LinearMap.Category.Instances +import Math.LinearMap.Category.Instances.Deriving import Math.LinearMap.Asserted import Math.VectorSpace.Docile import Math.LinearMap.Category.TensorQuot @@ -826,3 +836,13 @@ σm = mconcat [ Norm . arr $ m . (fmap ny $ m') | ((_,(_,Norm ny)), (m',m)) <- zip dataxy modelGens ] +-- $classesExplanation +-- To facilitate this library's abstract, base-less approach, the +-- spaces/types need to be instances of multiple typeclasses which +-- are rather tedious to write instances for. Due to the use of type +-- families, GHC's various deriving mechanisms cannot help much with this +-- either. +-- +-- In the common case that your spaces /do/ have a canonical basis, the +-- instances can be generated automatically with one of the Template Haskell +-- macros 'makeLinearSpaceFromBasis' or 'makeFiniteDimensionalFromBasis'. diff --git a/Math/LinearMap/Category/Class.hs b/Math/LinearMap/Category/Class.hs --- a/Math/LinearMap/Category/Class.hs +++ b/Math/LinearMap/Category/Class.hs @@ -14,6 +14,7 @@ {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE PatternSynonyms #-} @@ -24,6 +25,7 @@ {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE DefaultSignatures #-} +{-# LANGUAGE CPP #-} module Math.LinearMap.Category.Class where @@ -66,7 +68,11 @@ => ScalarSpaceWitness v data LinearManifoldWitness v where LinearManifoldWitness :: (Needle v ~ v, AffineSpace v, Diff v ~ v) - => BoundarylessWitness v -> LinearManifoldWitness v + => +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness v -> +#endif + LinearManifoldWitness v class (VectorSpace v, PseudoAffine v) => TensorSpace v where -- | The internal representation of a 'Tensor' product. @@ -225,6 +231,9 @@ , Scalar u ~ Scalar v, Scalar w ~ Scalar v ) => Bilinear ((v⊗u)+>w) (v⊗u) w + useTupleLinearSpaceComponents :: (v ~ (x,y)) + => ((LinearSpace x, LinearSpace y, Scalar x ~ Scalar y) => φ) -> φ + fmapLinearMap :: ∀ s v w x . ( LinearSpace v, TensorSpace w, TensorSpace x , Scalar v ~ s, Scalar w ~ s, Scalar x ~ s ) @@ -237,7 +246,10 @@ type TensorProduct (ZeroDim s) v = ZeroDim s scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of ClosedScalarWitness -> ScalarSpaceWitness - linearManifoldWitness = LinearManifoldWitness BoundarylessWitness + linearManifoldWitness = LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif zeroTensor = Tensor Origin toFlatTensor = LinearFunction $ \Origin -> Tensor Origin fromFlatTensor = LinearFunction $ \(Tensor Origin) -> Origin @@ -375,10 +387,12 @@ instance ∀ v w s . (LinearSpace v, TensorSpace w, Scalar v~s, Scalar w~s) => Semimanifold (LinearMap s v w) where type Needle (LinearMap s v w) = LinearMap s v w +#if !MIN_VERSION_manifolds_core(0,6,0) toInterior = pure fromInterior = id - (.+~^) = (^+^) translateP = Tagged (^+^) +#endif + (.+~^) = (^+^) instance ∀ v w s . (LinearSpace v, TensorSpace w, Scalar v~s, Scalar w~s) => PseudoAffine (LinearMap s v w) where f.-~.g = return $ f^-^g @@ -397,10 +411,12 @@ instance (TensorSpace v, TensorSpace w, Scalar v~s, Scalar w~s) => Semimanifold (Tensor s v w) where type Needle (Tensor s v w) = Tensor s v w +#if !MIN_VERSION_manifolds_core(0,6,0) toInterior = pure fromInterior = id - (.+~^) = (^+^) translateP = Tagged (^+^) +#endif + (.+~^) = (^+^) instance (TensorSpace v, TensorSpace w, Scalar v~s, Scalar w~s) => PseudoAffine (Tensor s v w) where f.-~.g = return $ f^-^g @@ -475,9 +491,19 @@ (ScalarSpaceWitness, ScalarSpaceWitness) -> ScalarSpaceWitness linearManifoldWitness = case ( linearManifoldWitness :: LinearManifoldWitness u , linearManifoldWitness :: LinearManifoldWitness v ) of - ( LinearManifoldWitness BoundarylessWitness - ,LinearManifoldWitness BoundarylessWitness ) - -> LinearManifoldWitness BoundarylessWitness + ( LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif + ,LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif + ) + -> LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif zeroTensor = zeroTensor <⊕ zeroTensor scaleTensor = bilinearFunction $ \μ (Tensor (v,w)) -> Tensor ( (scaleTensor-+$>μ)-+$>v, (scaleTensor-+$>μ)-+$>w ) @@ -567,6 +593,7 @@ \f (Tensor (tu,tv)) -> let LinearMap (fu,fv) = curryLinearMap $ f in ( (applyTensorLinMap-+$>uncurryLinearMap.asLinearMap $ fu)-+$>tu ) ^+^ ( (applyTensorLinMap-+$>uncurryLinearMap.asLinearMap $ fv)-+$>tv ) + useTupleLinearSpaceComponents r = r lfstBlock :: ( LSpace u, LSpace v, LSpace w , Scalar u ~ Scalar v, Scalar v ~ Scalar w ) @@ -616,9 +643,19 @@ , linearManifoldWitness :: LinearManifoldWitness u , linearManifoldWitness :: LinearManifoldWitness v ) of ( ScalarSpaceWitness - ,LinearManifoldWitness BoundarylessWitness - ,LinearManifoldWitness BoundarylessWitness ) - -> LinearManifoldWitness BoundarylessWitness + ,LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif + ,LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif + ) + -> LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif zeroTensor = deferLinearMap $ zeroV toFlatTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of ScalarSpaceWitness -> arr deferLinearMap . fmap toFlatTensor @@ -745,9 +782,19 @@ (ScalarSpaceWitness, ScalarSpaceWitness) -> ScalarSpaceWitness linearManifoldWitness = case ( linearManifoldWitness :: LinearManifoldWitness u , linearManifoldWitness :: LinearManifoldWitness v ) of - ( LinearManifoldWitness BoundarylessWitness - ,LinearManifoldWitness BoundarylessWitness ) - -> LinearManifoldWitness BoundarylessWitness + ( LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif + ,LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif + ) + -> LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif zeroTensor = lassocTensor $ zeroTensor toFlatTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of ScalarSpaceWitness -> arr lassocTensor . fmap toFlatTensor @@ -913,9 +960,19 @@ (ScalarSpaceWitness, ScalarSpaceWitness) -> ScalarSpaceWitness linearManifoldWitness = case ( linearManifoldWitness :: LinearManifoldWitness u , linearManifoldWitness :: LinearManifoldWitness v ) of - ( LinearManifoldWitness BoundarylessWitness - ,LinearManifoldWitness BoundarylessWitness ) - -> LinearManifoldWitness BoundarylessWitness + ( LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif + ,LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif + ) + -> LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif zeroTensor = fromLinearFn $ const0 toFlatTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of ScalarSpaceWitness -> fmap fromLinearFn $ applyDualVector @@ -1179,8 +1236,14 @@ ScalarSpaceWitness -> ScalarSpaceWitness linearManifoldWitness = case linearManifoldWitness :: LinearManifoldWitness (Needle (VRep m)) of - LinearManifoldWitness BoundarylessWitness - -> LinearManifoldWitness BoundarylessWitness + LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif + -> LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif zeroTensor = pseudoFmapTensorLHS GenericNeedle $ zeroTensor toFlatTensor = LinearFunction $ arr (pseudoFmapTensorLHS GenericNeedle) . getLinearFunction toFlatTensor @@ -1271,9 +1334,11 @@ => Semimanifold (GenericTupleDual f g p) where type Needle (GenericTupleDual f g p) = GenericTupleDual f g p (.+~^) = (^+^) +#if !MIN_VERSION_manifolds_core(0,6,0) fromInterior = id toInterior = pure translateP = Tagged (^+^) +#endif instance (AdditiveGroup (DualVector (f p)), AdditiveGroup (DualVector (g p))) => PseudoAffine (GenericTupleDual f g p) where p.-~.q = Just $ p.-.q @@ -1299,11 +1364,17 @@ (fmap fromTensor $ wellDefinedTensor (fromLinearMap $ gt)) scalarSpaceWitness = case scalarSpaceWitness :: ScalarSpaceWitness (f p) of ScalarSpaceWitness -> ScalarSpaceWitness - linearManifoldWitness = LinearManifoldWitness BoundarylessWitness + linearManifoldWitness = LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif zeroTensor = case ( linearManifoldWitness :: LinearManifoldWitness (f p) , dualSpaceWitness :: DualSpaceWitness (f p) , dualSpaceWitness :: DualSpaceWitness (g p) ) of - ( LinearManifoldWitness BoundarylessWitness + ( LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif ,DualSpaceWitness, DualSpaceWitness ) -> Tensor (fromTensor $ zeroTensor, fromTensor $ zeroTensor) toFlatTensor = case ( scalarSpaceWitness :: ScalarSpaceWitness (f p) @@ -1439,11 +1510,13 @@ (.+^) = (^+^) instance AdditiveGroup (DualVector (Needle (VRep m))) => Semimanifold (GenericNeedle' m) where - type Interior (GenericNeedle' m) = GenericNeedle' m type Needle (GenericNeedle' m) = GenericNeedle' m +#if !MIN_VERSION_manifolds_core(0,6,0) + type Interior (GenericNeedle' m) = GenericNeedle' m toInterior = pure fromInterior = id translateP = Tagged (^+^) +#endif (.+~^) = (^+^) instance AdditiveGroup (DualVector (Needle (VRep m))) => PseudoAffine (GenericNeedle' m) where @@ -1462,8 +1535,14 @@ ScalarSpaceWitness -> ScalarSpaceWitness linearManifoldWitness = case linearManifoldWitness :: LinearManifoldWitness (DualVector (Needle (VRep m))) of - LinearManifoldWitness BoundarylessWitness - -> LinearManifoldWitness BoundarylessWitness + LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif + -> LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif zeroTensor = pseudoFmapTensorLHS GenericNeedle' $ zeroTensor toFlatTensor = LinearFunction $ arr (pseudoFmapTensorLHS GenericNeedle') . getLinearFunction toFlatTensor diff --git a/Math/LinearMap/Category/Instances.hs b/Math/LinearMap/Category/Instances.hs --- a/Math/LinearMap/Category/Instances.hs +++ b/Math/LinearMap/Category/Instances.hs @@ -18,6 +18,8 @@ {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE CPP #-} {-# LANGUAGE TupleSections #-} +{-# LANGUAGE EmptyCase #-} +{-# LANGUAGE TypeApplications #-} module Math.LinearMap.Category.Instances where @@ -26,6 +28,9 @@ import Data.VectorSpace import Data.Basis +#if MIN_VERSION_manifolds_core(0,6,0) +import Math.Manifold.Core.Types (EmptyMfd) +#endif import Math.Manifold.Core.PseudoAffine import Prelude () @@ -61,19 +66,42 @@ import qualified GHC.Exts as GHC import qualified GHC.Generics as GHC + +#if MIN_VERSION_manifolds_core(0,6,0) +instance LinearSpace v => Semimanifold (EmptyMfd v) where + type Needle (EmptyMfd v) = v + p .+~^ _ = case p of {} + p .-~^ _ = case p of {} + semimanifoldWitness = case linearManifoldWitness @v of + LinearManifoldWitness -> SemimanifoldWitness +instance LinearSpace v => PseudoAffine (EmptyMfd v) where + p .-~. _ = case p of {} +#endif + + infixr 7 <.>^ (<.>^) :: LinearSpace v => DualVector v -> v -> Scalar v f<.>^v = (applyDualVector-+$>f)-+$>v - type ℝ = Double +autoLinearManifoldWitness :: (Semimanifold v, AffineSpace v, v ~ Needle v, v ~ Diff v +#if !MIN_VERSION_manifolds_core(0,6,0) + , v ~ Interior v +#endif + ) + => LinearManifoldWitness v +autoLinearManifoldWitness = LinearManifoldWitness +#if !MIN_VERSION_manifolds_core(0,6,0) + BoundarylessWitness +#endif + #define LinearScalarSpace(S) \ instance Num' (S) where {closedScalarWitness = ClosedScalarWitness}; \ instance TensorSpace (S) where { \ type TensorProduct (S) w = w; \ scalarSpaceWitness = ScalarSpaceWitness; \ - linearManifoldWitness = LinearManifoldWitness BoundarylessWitness; \ + linearManifoldWitness = autoLinearManifoldWitness; \ zeroTensor = Tensor zeroV; \ scaleTensor = bilinearFunction $ \μ (Tensor t) -> Tensor $ μ*^t; \ addTensors (Tensor v) (Tensor w) = Tensor $ v ^+^ w; \ @@ -111,19 +139,26 @@ LinearScalarSpace(Float) LinearScalarSpace(Rational) + +#if MIN_VERSION_manifolds_core(0,6,0) +#define FreeLinSpaceInteriorDecls +#else +#define FreeLinSpaceInteriorDecls \ + toInterior = pure; fromInterior = id; translateP = Tagged (^+^); +#endif + #define FreeLinearSpace(V, LV, tp, tenspl, tenid, dspan, contraction, contraaction) \ instance Num s => Semimanifold (V s) where { \ type Needle (V s) = V s; \ - toInterior = pure; fromInterior = id; \ - (.+~^) = (^+^); \ - translateP = Tagged (^+^) }; \ + FreeLinSpaceInteriorDecls \ + (.+~^) = (^+^) }; \ instance Num s => PseudoAffine (V s) where { \ v.-~.w = pure (v^-^w); (.-~!) = (^-^) }; \ instance ∀ s . (Num' s, Eq s) => TensorSpace (V s) where { \ type TensorProduct (V s) w = V w; \ scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of{ \ ClosedScalarWitness -> ScalarSpaceWitness}; \ - linearManifoldWitness = LinearManifoldWitness BoundarylessWitness; \ + linearManifoldWitness = autoLinearManifoldWitness; \ zeroTensor = Tensor $ pure zeroV; \ addTensors (Tensor m) (Tensor n) = Tensor $ liftA2 (^+^) m n; \ subtractTensors (Tensor m) (Tensor n) = Tensor $ liftA2 (^-^) m n; \ @@ -261,11 +296,12 @@ - instance (Num' n, UArr.Unbox n) => Semimanifold (FinSuppSeq n) where type Needle (FinSuppSeq n) = FinSuppSeq n - (.+~^) = (.+^); translateP = Tagged (.+^) - toInterior = pure; fromInterior = id + (.+~^) = (.+^) +#if !MIN_VERSION_manifolds_core(0,6,0) + translateP = Tagged (.+^); toInterior = pure; fromInterior = id +#endif instance (Num' n, UArr.Unbox n) => PseudoAffine (FinSuppSeq n) where v.-~.w = Just $ v.-.w; (.-~!) = (.-.) @@ -275,7 +311,7 @@ wellDefinedVector (FinSuppSeq v) = FinSuppSeq <$> UArr.mapM wellDefinedVector v scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness n of ClosedScalarWitness -> ScalarSpaceWitness - linearManifoldWitness = LinearManifoldWitness BoundarylessWitness + linearManifoldWitness = autoLinearManifoldWitness zeroTensor = Tensor [] toFlatTensor = LinearFunction $ Tensor . UArr.toList . getFiniteSeq fromFlatTensor = LinearFunction $ FinSuppSeq . UArr.fromList . getTensorProduct @@ -299,8 +335,10 @@ instance (Num' n, UArr.Unbox n) => Semimanifold (Sequence n) where type Needle (Sequence n) = Sequence n - (.+~^) = (.+^); translateP = Tagged (.+^) - toInterior = pure; fromInterior = id + (.+~^) = (.+^) +#if !MIN_VERSION_manifolds_core(0,6,0) + translateP = Tagged (.+^); toInterior = pure; fromInterior = id +#endif instance (Num' n, UArr.Unbox n) => PseudoAffine (Sequence n) where v.-~.w = Just $ v.-.w; (.-~!) = (.-.) @@ -313,7 +351,7 @@ wellDefinedTensor (Tensor a) = Tensor <$> Hask.traverse wellDefinedVector a scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness n of ClosedScalarWitness -> ScalarSpaceWitness - linearManifoldWitness = LinearManifoldWitness BoundarylessWitness + linearManifoldWitness = autoLinearManifoldWitness zeroTensor = Tensor [] toFlatTensor = LinearFunction $ Tensor . GHC.toList fromFlatTensor = LinearFunction $ GHC.fromList . getTensorProduct @@ -419,9 +457,11 @@ instance (TensorSpace v, Scalar v ~ s) => Semimanifold (SymmetricTensor s v) where type Needle (SymmetricTensor s v) = SymmetricTensor s v (.+~^) = (^+^) +#if !MIN_VERSION_manifolds_core(0,6,0) fromInterior = id toInterior = pure translateP = Tagged (^+^) +#endif instance (TensorSpace v, Scalar v ~ s) => PseudoAffine (SymmetricTensor s v) where (.-~!) = (^-^) instance (Num' s, TensorSpace v, Scalar v ~ s) => TensorSpace (SymmetricTensor s v) where @@ -429,7 +469,7 @@ wellDefinedVector (SymTensor t) = SymTensor <$> wellDefinedVector t scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of ClosedScalarWitness -> ScalarSpaceWitness - linearManifoldWitness = LinearManifoldWitness BoundarylessWitness + linearManifoldWitness = autoLinearManifoldWitness zeroTensor = Tensor zeroV toFlatTensor = case closedScalarWitness :: ClosedScalarWitness s of ClosedScalarWitness -> LinearFunction $ \(SymTensor t) @@ -548,9 +588,12 @@ , Monoidal f (LinearFunction (Scalar y)) (LinearFunction (Scalar y)) ) => Semimanifold (LinearApplicativeSpace f y) where type Needle (LinearApplicativeSpace f y) = LinearApplicativeSpace f y +#if !MIN_VERSION_manifolds_core(0,6,0) type Interior (LinearApplicativeSpace f y) = LinearApplicativeSpace f y toInterior = Just; fromInterior = id translateP = Tagged (^+^) +#endif + (.+~^) = (^+^) instance ( GHC.Generic1 f, TensorSpace y , TensorSpace (f y), Scalar (f y) ~ Scalar y @@ -563,6 +606,9 @@ instance (InnerSpace v, Scalar v ~ ℝ, TensorSpace v) => InnerSpace (Tensor ℝ ℝ v) where Tensor t <.> Tensor u = t <.> u +instance (InnerSpace v, TensorSpace v, Scalar v ~ ℝ) + => InnerSpace (LinearMap ℝ ℝ v) where + LinearMap f <.> LinearMap g = f<.>g instance (Show v) => Show (Tensor ℝ ℝ v) where showsPrec p (Tensor t) = showParen (p>9) $ ("Tensor "++) . showsPrec 10 t @@ -571,3 +617,18 @@ arbitrary = Tensor <$> QC.arbitrary shrink (Tensor t) = Tensor <$> QC.shrink t +instance (QC.Arbitrary v, Scalar v ~ ℝ) => QC.Arbitrary (LinearMap ℝ ℝ v) where + arbitrary = LinearMap <$> QC.arbitrary + shrink (LinearMap t) = LinearMap <$> QC.shrink t + +#define FreeArbitrarySpace(S) \ +instance (QC.Arbitrary v, Scalar v ~ ℝ) => QC.Arbitrary (Tensor ℝ (S ℝ) v) where { \ + arbitrary = Tensor <$> Hask.traverse (const QC.arbitrary) zeroV }; \ +instance (QC.Arbitrary v, Scalar v ~ ℝ) => QC.Arbitrary (LinearMap ℝ (S ℝ) v) where { \ + arbitrary = LinearMap <$> Hask.traverse (const QC.arbitrary) zeroV } + +FreeArbitrarySpace(V0) +FreeArbitrarySpace(V1) +FreeArbitrarySpace(V2) +FreeArbitrarySpace(V3) +FreeArbitrarySpace(V4) diff --git a/Math/LinearMap/Category/Instances/Deriving.hs b/Math/LinearMap/Category/Instances/Deriving.hs new file mode 100644 --- /dev/null +++ b/Math/LinearMap/Category/Instances/Deriving.hs @@ -0,0 +1,595 @@ +-- | +-- Module : Math.LinearMap.Instances.Deriving +-- Copyright : (c) Justus Sagemüller 2021 +-- License : GPL v3 +-- +-- Maintainer : (@) jsag $ hvl.no +-- Stability : experimental +-- Portability : portable +-- +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE DerivingStrategies #-} +{-# LANGUAGE GeneralisedNewtypeDeriving #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE TupleSections #-} + +module Math.LinearMap.Category.Instances.Deriving + ( makeLinearSpaceFromBasis, makeFiniteDimensionalFromBasis + -- * The instantiated classes + , AffineSpace(..), Semimanifold(..), PseudoAffine(..) + , TensorSpace(..), LinearSpace(..), FiniteDimensional(..), SemiInner(..) + -- * Internals + , BasisGeneratedSpace(..), LinearSpaceFromBasisDerivationConfig, def ) where + +import Math.LinearMap.Category.Class +import Math.VectorSpace.Docile + +import Data.VectorSpace +import Data.AffineSpace +import Data.Basis +import qualified Data.Map as Map +import Data.MemoTrie +import Data.Hashable + +import Prelude () +import qualified Prelude as Hask + +import Control.Category.Constrained.Prelude +import Control.Arrow.Constrained + +import Data.Coerce +import Data.Type.Coercion +import Data.Tagged +import Data.Traversable (traverse) +import Data.Default.Class + +import Math.Manifold.Core.PseudoAffine +import Math.LinearMap.Asserted +import Math.VectorSpace.ZeroDimensional +import Data.VectorSpace.Free + +import Language.Haskell.TH + +-- | Given a type @V@ that is already a 'VectorSpace' and 'HasBasis', generate +-- the other class instances that are needed to use the type with this +-- library. +-- +-- Prerequisites: (these can often be derived automatically, +-- using either the @newtype@ \/ @via@ strategy or generics \/ anyclass) +-- +-- @ +-- instance 'AdditiveGroup' V +-- +-- instance 'VectorSpace' V where +-- type Scalar V = -- a simple number type, usually 'Double' +-- +-- instance 'HasBasis' V where +-- type Basis V = -- a type with an instance of 'HasTrie' +-- @ +-- +-- Note that the 'Basis' does /not/ need to be orthonormal – in fact it +-- is not necessary to have a scalar product (i.e. an 'InnerSpace' instance) +-- at all. +-- +-- This macro, invoked like +-- @ +-- makeLinearSpaceFromBasis [t| V |] +-- @ +-- +-- will then generate @V@-instances for the classes 'Semimanifold', +-- 'PseudoAffine', 'AffineSpace', 'TensorSpace' and 'LinearSpace'. +makeLinearSpaceFromBasis :: Q Type -> DecsQ +makeLinearSpaceFromBasis v + = makeLinearSpaceFromBasis' def $ deQuantifyType v + +data LinearSpaceFromBasisDerivationConfig = LinearSpaceFromBasisDerivationConfig +instance Default LinearSpaceFromBasisDerivationConfig where + def = LinearSpaceFromBasisDerivationConfig + +-- | More general version of 'makeLinearSpaceFromBasis', that can be used with +-- parameterised types. +makeLinearSpaceFromBasis' :: LinearSpaceFromBasisDerivationConfig + -> Q (Cxt, Type) -> DecsQ +makeLinearSpaceFromBasis' _ cxtv = do + (cxt,v) <- do + (cxt', v') <- cxtv + return (pure cxt', pure v') + + exts <- extsEnabled + if not $ TypeFamilies`elem`exts && ScopedTypeVariables`elem`exts + then reportError "This macro requires -XTypeFamilies and -XScopedTypeVariables." + else pure () + + sequence + [ InstanceD Nothing <$> cxt <*> [t|Semimanifold $v|] <*> do + tySyns <- sequence [ +#if MIN_VERSION_template_haskell(2,15,0) + error "The TH type of TySynInstD has changed" +#else + TySynInstD ''Needle <$> do + TySynEqn . (:[]) <$> v <*> v + , TySynInstD ''Interior <$> do + TySynEqn . (:[]) <$> v <*> v +#endif + ] + methods <- [d| + $(varP 'toInterior) = pure + $(varP 'fromInterior) = id + $(varP 'translateP) = Tagged (^+^) + $(varP '(.+~^)) = (^+^) + $(varP 'semimanifoldWitness) = SemimanifoldWitness BoundarylessWitness + |] + return $ tySyns ++ methods + , InstanceD Nothing <$> cxt <*> [t|PseudoAffine $v|] <*> do + [d| + $(varP '(.-~!)) = (^-^) + |] + , InstanceD Nothing <$> cxt <*> [t|AffineSpace $v|] <*> do + tySyns <- sequence [ +#if MIN_VERSION_template_haskell(2,15,0) + error "The TH type of TySynInstD has changed" +#else + TySynInstD ''Diff <$> do + TySynEqn . (:[]) <$> v <*> v +#endif + ] + methods <- [d| + $(varP '(.+^)) = (^+^) + $(varP '(.-.)) = (^-^) + |] + return $ tySyns ++ methods + , InstanceD Nothing <$> cxt <*> [t|TensorSpace $v|] <*> do + tySyns <- sequence [ +#if MIN_VERSION_template_haskell(2,15,0) + error "The TH type of TySynInstD has changed" +#else + TySynInstD ''TensorProduct <$> do + wType <- VarT <$> newName "w" :: Q Type + TySynEqn . (:[wType]) <$> v + <*> [t| Basis $v :->: $(pure wType) |] +#endif + ] + methods <- [d| + $(varP 'wellDefinedVector) = \v + -> if v==v then Just v else Nothing + $(varP 'wellDefinedTensor) = \(Tensor v) + -> fmap (const $ Tensor v) . traverse (wellDefinedVector . snd) $ enumerate v + $(varP 'zeroTensor) = Tensor . trie $ const zeroV + $(varP 'toFlatTensor) = LinearFunction $ Tensor . trie . decompose' + $(varP 'fromFlatTensor) = LinearFunction $ \(Tensor t) + -> recompose $ enumerate t + $(varP 'scalarSpaceWitness) = ScalarSpaceWitness + $(varP 'linearManifoldWitness) = LinearManifoldWitness BoundarylessWitness + $(varP 'addTensors) = \(Tensor v) (Tensor w) + -> Tensor $ (^+^) <$> v <*> w + $(varP 'subtractTensors) = \(Tensor v) (Tensor w) + -> Tensor $ (^-^) <$> v <*> w + $(varP 'tensorProduct) = bilinearFunction + $ \v w -> Tensor . trie $ \bv -> decompose' v bv *^ w + $(varP 'transposeTensor) = LinearFunction $ \(Tensor t) + -> sumV [ (tensorProduct-+$>w)-+$>basisValue b + | (b,w) <- enumerate t ] + $(varP 'fmapTensor) = bilinearFunction + $ \(LinearFunction f) (Tensor t) + -> Tensor $ fmap f t + $(varP 'fzipTensorWith) = bilinearFunction + $ \(LinearFunction f) (Tensor tv, Tensor tw) + -> Tensor $ liftA2 (curry f) tv tw + $(varP 'coerceFmapTensorProduct) = \_ Coercion + -> error "Cannot yet coerce tensors defined from a `HasBasis` instance. This would require `RoleAnnotations` on `:->:`. Cf. https://gitlab.haskell.org/ghc/ghc/-/issues/8177" + |] + return $ tySyns ++ methods + , InstanceD Nothing <$> cxt <*> [t|BasisGeneratedSpace $v|] <*> do + [d| + $(varP 'proveTensorProductIsTrie) = \φ -> φ + |] + , InstanceD Nothing <$> cxt <*> [t|LinearSpace $v|] <*> do + tySyns <- sequence [ +#if MIN_VERSION_template_haskell(2,15,0) + error "The TH type of TySynInstD has changed" +#else + TySynInstD ''DualVector <$> do + TySynEqn . (:[]) <$> v + <*> [t| DualVectorFromBasis $v |] +#endif + ] + methods <- [d| + + + $(varP 'dualSpaceWitness) = case closedScalarWitness @(Scalar $v) of + ClosedScalarWitness -> DualSpaceWitness + $(varP 'linearId) = LinearMap . trie $ basisValue + $(varP 'tensorId) = tid + where tid :: ∀ w . (LinearSpace w, Scalar w ~ Scalar $v) + => ($v⊗w) +> ($v⊗w) + tid = case dualSpaceWitness @w of + DualSpaceWitness -> LinearMap . trie $ Tensor . \i + -> getTensorProduct $ + (fmapTensor @(DualVector w) + -+$>(LinearFunction $ \w -> Tensor . trie + $ (\j -> if i==j then w else zeroV) + :: $v⊗w)) + -+$> case linearId @w of + LinearMap lw -> Tensor lw :: DualVector w⊗w + $(varP 'applyDualVector) = bilinearFunction + $ \(DualVectorFromBasis f) v + -> sum [decompose' f i * vi | (i,vi) <- decompose v] + $(varP 'applyLinear) = bilinearFunction + $ \(LinearMap f) v + -> sumV [vi *^ untrie f i | (i,vi) <- decompose v] + $(varP 'applyTensorFunctional) = atf + where atf :: ∀ u . (LinearSpace u, Scalar u ~ Scalar $v) + => Bilinear (DualVector ($v ⊗ u)) + ($v ⊗ u) (Scalar $v) + atf = case dualSpaceWitness @u of + DualSpaceWitness -> bilinearFunction + $ \(LinearMap f) (Tensor t) + -> sum [ (applyDualVector-+$>fi)-+$>untrie t i + | (i, fi) <- enumerate f ] + $(varP 'applyTensorLinMap) = atlm + where atlm :: ∀ u w . ( LinearSpace u, TensorSpace w + , Scalar u ~ Scalar $v, Scalar w ~ Scalar $v ) + => Bilinear (($v ⊗ u) +> w) ($v ⊗ u) w + atlm = case dualSpaceWitness @u of + DualSpaceWitness -> bilinearFunction + $ \(LinearMap f) (Tensor t) + -> sumV [ (applyLinear-+$>(LinearMap fi :: u+>w)) + -+$> untrie t i + | (i, Tensor fi) <- enumerate f ] + $(varP 'useTupleLinearSpaceComponents) = error "Not a tuple type" + + |] + return $ tySyns ++ methods + ] + +data FiniteDimensionalFromBasisDerivationConfig + = FiniteDimensionalFromBasisDerivationConfig +instance Default FiniteDimensionalFromBasisDerivationConfig where + def = FiniteDimensionalFromBasisDerivationConfig + +-- | Like 'makeLinearSpaceFromBasis', but additionally generate instances for +-- 'FiniteDimensional' and 'SemiInner'. +makeFiniteDimensionalFromBasis :: Q Type -> DecsQ +makeFiniteDimensionalFromBasis v + = makeFiniteDimensionalFromBasis' def $ deQuantifyType v + +makeFiniteDimensionalFromBasis' :: FiniteDimensionalFromBasisDerivationConfig + -> Q (Cxt, Type) -> DecsQ +makeFiniteDimensionalFromBasis' _ cxtv = do + generalInsts <- makeLinearSpaceFromBasis' def cxtv + (cxt,v) <- do + (cxt', v') <- cxtv + return (pure cxt', pure v') + vtnameHash <- abs . hash . show <$> v + + fdInsts <- sequence + [ InstanceD Nothing <$> cxt <*> [t|FiniteDimensional $v|] <*> do + + -- This is a hack. Ideally, @newName@ should generate globally unique names, + -- but it doesn't, so we append a hash of the vector space type. + -- Cf. https://gitlab.haskell.org/ghc/ghc/-/issues/13054 + subBasisCstr <- newName $ "CompleteBasis"++show vtnameHash + + tySyns <- sequence [ +#if MIN_VERSION_template_haskell(2,15,0) + error "The TH type of TySynInstD has changed" +#else + DataInstD [] ''SubBasis + <$> ((:[]) <$> v) + <*> pure Nothing + <*> pure [NormalC subBasisCstr []] + <*> pure [] +#endif + ] + methods <- [d| + $(varP 'entireBasis) = $(conE subBasisCstr) + $(varP 'enumerateSubBasis) = + \ $(conP subBasisCstr []) -> basisValue . fst <$> enumerate (trie $ const ()) + $(varP 'tensorEquality) + = \(Tensor t) (Tensor t') -> and [ti == untrie t' i | (i,ti) <- enumerate t] + $(varP 'decomposeLinMap) = dlm + where dlm :: ∀ w . ($v+>w) + -> (SubBasis $v, [w]->[w]) + dlm (LinearMap f) = + ( $(conE subBasisCstr) + , (map snd (enumerate f) ++) ) + $(varP 'decomposeLinMapWithin) = dlm + where dlm :: ∀ w . SubBasis $v + -> ($v+>w) + -> Either (SubBasis $v, [w]->[w]) + ([w]->[w]) + dlm $(conP subBasisCstr []) (LinearMap f) = + (Right (map snd (enumerate f) ++) ) + $(varP 'recomposeSB) = rsb + where rsb :: SubBasis $v + -> [Scalar $v] + -> ($v, [Scalar $v]) + rsb $(conP subBasisCstr []) cs = first recompose + $ zipWith' (,) (fst <$> enumerate (trie $ const ())) cs + $(varP 'recomposeSBTensor) = rsbt + where rsbt :: ∀ w . (FiniteDimensional w, Scalar w ~ Scalar $v) + => SubBasis $v -> SubBasis w + -> [Scalar $v] + -> ($v⊗w, [Scalar $v]) + rsbt $(conP subBasisCstr []) sbw ws = + (first (\iws -> Tensor $ trie (Map.fromList iws Map.!)) + $ zipConsumeWith' (\i cs' -> first (\c->(i,c)) + $ recomposeSB sbw cs') + (fst <$> enumerate (trie $ const ())) ws) + $(varP 'recomposeLinMap) = rlm + where rlm :: ∀ w . SubBasis $v + -> [w] + -> ($v+>w, [w]) + rlm $(conP subBasisCstr []) ws = + (first (\iws -> LinearMap $ trie (Map.fromList iws Map.!)) + $ zipWith' (,) (fst <$> enumerate (trie $ const ())) ws) + $(varP 'recomposeContraLinMap) = rclm + where rclm :: ∀ w f . (LinearSpace w, Scalar w ~ Scalar $v, Hask.Functor f) + => (f (Scalar w) -> w) -> f (DualVectorFromBasis $v) + -> ($v+>w) + rclm f vs = + (LinearMap $ trie (\i -> f $ fmap (`decompose'`i) vs)) + $(varP 'recomposeContraLinMapTensor) = rclm + where rclm :: ∀ u w f + . ( FiniteDimensional u, LinearSpace w + , Scalar u ~ Scalar $v, Scalar w ~ Scalar $v, Hask.Functor f + ) + => (f (Scalar w) -> w) -> f ($v+>DualVector u) + -> (($v⊗u)+>w) + rclm f vus = case dualSpaceWitness @u of + DualSpaceWitness -> + ( + (LinearMap $ trie + (\i -> case recomposeContraLinMap @u @w @f f + $ fmap (\(LinearMap vu) -> untrie vu (i :: Basis $v)) vus of + LinearMap wuff -> Tensor wuff :: DualVector u⊗w ))) + $(varP 'uncanonicallyFromDual) = LinearFunction getDualVectorFromBasis + $(varP 'uncanonicallyToDual) = LinearFunction DualVectorFromBasis + + |] + return $ tySyns ++ methods + , InstanceD Nothing <$> cxt <*> [t|SemiInner $v|] <*> do + [d| + $(varP 'dualBasisCandidates) + = cartesianDualBasisCandidates + (enumerateSubBasis CompleteDualVBasis) + (\v -> map (abs . realToFrac . decompose' v . fst) + $ enumerate (trie $ const ()) ) + |] + ] + return $ generalInsts ++ fdInsts + + +deQuantifyType :: Q Type -> Q (Cxt, Type) +deQuantifyType t = do + t' <- t + return $ case t' of + ForallT _ cxt instT -> (cxt, instT) + _ -> ([], t') + + +newtype DualVectorFromBasis v = DualVectorFromBasis { getDualVectorFromBasis :: v } + deriving newtype (Eq, AdditiveGroup, VectorSpace, HasBasis) + +instance AdditiveGroup v => Semimanifold (DualVectorFromBasis v) where + type Needle (DualVectorFromBasis v) = DualVectorFromBasis v + type Interior (DualVectorFromBasis v) = DualVectorFromBasis v + toInterior = pure + fromInterior = id + translateP = Tagged (^+^) + (.+~^) = (^+^) + semimanifoldWitness = SemimanifoldWitness BoundarylessWitness + +instance AdditiveGroup v => AffineSpace (DualVectorFromBasis v) where + type Diff (DualVectorFromBasis v) = DualVectorFromBasis v + (.+^) = (^+^) + (.-.) = (^-^) + +instance AdditiveGroup v => PseudoAffine (DualVectorFromBasis v) where + (.-~!) = (^-^) + +instance ∀ v . ( HasBasis v, Num' (Scalar v) + , Scalar (Scalar v) ~ Scalar v + , HasTrie (Basis v) + , Eq v ) + => TensorSpace (DualVectorFromBasis v) where + type TensorProduct (DualVectorFromBasis v) w = Basis v :->: w + wellDefinedVector v + | v==v = Just v + | otherwise = Nothing + wellDefinedTensor (Tensor v) + = fmap (const $ Tensor v) . traverse (wellDefinedVector . snd) $ enumerate v + zeroTensor = Tensor . trie $ const zeroV + toFlatTensor = LinearFunction $ Tensor . trie . decompose' + fromFlatTensor = LinearFunction $ \(Tensor t) + -> recompose $ enumerate t + scalarSpaceWitness = ScalarSpaceWitness + linearManifoldWitness = LinearManifoldWitness BoundarylessWitness + addTensors (Tensor v) (Tensor w) = Tensor $ (^+^) <$> v <*> w + subtractTensors (Tensor v) (Tensor w) = Tensor $ (^-^) <$> v <*> w + tensorProduct = bilinearFunction + $ \v w -> Tensor . trie $ \bv -> decompose' v bv *^ w + transposeTensor = LinearFunction $ \(Tensor t) + -> sumV [ (tensorProduct-+$>w)-+$>basisValue b + | (b,w) <- enumerate t ] + fmapTensor = bilinearFunction + $ \(LinearFunction f) (Tensor t) + -> Tensor $ fmap f t + fzipTensorWith = bilinearFunction + $ \(LinearFunction f) (Tensor tv, Tensor tw) + -> Tensor $ liftA2 (curry f) tv tw + coerceFmapTensorProduct _ Coercion + = error "Cannot yet coerce tensors defined from a `HasBasis` instance. This would require `RoleAnnotations` on `:->:`. Cf. https://gitlab.haskell.org/ghc/ghc/-/issues/8177" + + +-- | Do not manually instantiate this class. It is used internally +-- by 'makeLinearSpaceFromBasis'. +class ( HasBasis v, Num' (Scalar v) + , LinearSpace v, DualVector v ~ DualVectorFromBasis v) + => BasisGeneratedSpace v where + proveTensorProductIsTrie + :: ∀ w φ . (TensorProduct v w ~ (Basis v :->: w) => φ) -> φ + +instance ∀ v . ( BasisGeneratedSpace v + , Scalar (Scalar v) ~ Scalar v + , HasTrie (Basis v) + , Eq v, Eq (Basis v) ) + => LinearSpace (DualVectorFromBasis v) where + type DualVector (DualVectorFromBasis v) = v + dualSpaceWitness = case closedScalarWitness @(Scalar v) of + ClosedScalarWitness -> DualSpaceWitness + linearId = proveTensorProductIsTrie @v @(DualVectorFromBasis v) + (LinearMap . trie $ DualVectorFromBasis . basisValue) + tensorId = tid + where tid :: ∀ w . (LinearSpace w, Scalar w ~ Scalar v) + => (DualVectorFromBasis v⊗w) +> (DualVectorFromBasis v⊗w) + tid = proveTensorProductIsTrie @v @(DualVector w⊗(DualVectorFromBasis v⊗w)) + ( case dualSpaceWitness @w of + DualSpaceWitness -> LinearMap . trie $ Tensor . \i + -> getTensorProduct $ + (fmapTensor @(DualVector w) + -+$>(LinearFunction $ \w -> Tensor . trie + $ (\j -> if i==j then w else zeroV) + :: DualVectorFromBasis v⊗w)) + -+$> case linearId @w of + LinearMap lw -> Tensor lw :: DualVector w⊗w ) + applyDualVector = proveTensorProductIsTrie @v @(DualVectorFromBasis v) + ( bilinearFunction $ \f (DualVectorFromBasis v) + -> sum [decompose' f i * vi | (i,vi) <- decompose v] ) + applyLinear = ali + where ali :: ∀ w . (TensorSpace w, Scalar w~Scalar v) + => Bilinear (DualVectorFromBasis v +> w) (DualVectorFromBasis v) w + ali = proveTensorProductIsTrie @v @w ( bilinearFunction + $ \(LinearMap f) (DualVectorFromBasis v) + -> sumV [vi *^ untrie f i | (i,vi) <- decompose v] ) + applyTensorFunctional = atf + where atf :: ∀ u . (LinearSpace u, Scalar u ~ Scalar v) + => Bilinear (DualVector (DualVectorFromBasis v ⊗ u)) + (DualVectorFromBasis v ⊗ u) (Scalar v) + atf = proveTensorProductIsTrie @v @(DualVector u) (case dualSpaceWitness @u of + DualSpaceWitness -> bilinearFunction + $ \(LinearMap f) (Tensor t) + -> sum [ (applyDualVector-+$>fi)-+$>untrie t i + | (i, fi) <- enumerate f ] + ) + applyTensorLinMap = atlm + where atlm :: ∀ u w . ( LinearSpace u, TensorSpace w + , Scalar u ~ Scalar v, Scalar w ~ Scalar v ) + => Bilinear ((DualVectorFromBasis v ⊗ u) +> w) + (DualVectorFromBasis v ⊗ u) w + atlm = proveTensorProductIsTrie @v @(DualVector u⊗w) ( + case dualSpaceWitness @u of + DualSpaceWitness -> bilinearFunction + $ \(LinearMap f) (Tensor t) + -> sumV [ (applyLinear-+$>(LinearMap fi :: u+>w)) + -+$> untrie t i + | (i, Tensor fi) <- enumerate f ] + ) + useTupleLinearSpaceComponents = error "Not a tuple type" + + +zipWith' :: (a -> b -> c) -> [a] -> [b] -> ([c], [b]) +zipWith' _ _ [] = ([], []) +zipWith' _ [] ys = ([], ys) +zipWith' f (x:xs) (y:ys) = first (f x y :) $ zipWith' f xs ys + +zipConsumeWith' :: (a -> [b] -> (c,[b])) -> [a] -> [b] -> ([c], [b]) +zipConsumeWith' _ _ [] = ([], []) +zipConsumeWith' _ [] ys = ([], ys) +zipConsumeWith' f (x:xs) ys + = case f x ys of + (z, ys') -> first (z :) $ zipConsumeWith' f xs ys' + +instance ∀ v . ( BasisGeneratedSpace v, FiniteDimensional v + , Scalar (Scalar v) ~ Scalar v + , HasTrie (Basis v), Ord (Basis v) + , Eq v, Eq (Basis v) ) + => FiniteDimensional (DualVectorFromBasis v) where + data SubBasis (DualVectorFromBasis v) = CompleteDualVBasis + entireBasis = CompleteDualVBasis + enumerateSubBasis CompleteDualVBasis + = basisValue . fst <$> enumerate (trie $ const ()) + tensorEquality (Tensor t) (Tensor t') + = and [ti == untrie t' i | (i,ti) <- enumerate t] + decomposeLinMap = dlm + where dlm :: ∀ w . (DualVectorFromBasis v+>w) + -> (SubBasis (DualVectorFromBasis v), [w]->[w]) + dlm (LinearMap f) = proveTensorProductIsTrie @v @w + ( CompleteDualVBasis + , (map snd (enumerate f) ++) ) + decomposeLinMapWithin = dlm + where dlm :: ∀ w . SubBasis (DualVectorFromBasis v) + -> (DualVectorFromBasis v+>w) + -> Either (SubBasis (DualVectorFromBasis v), [w]->[w]) + ([w]->[w]) + dlm CompleteDualVBasis (LinearMap f) = proveTensorProductIsTrie @v @w + (Right (map snd (enumerate f) ++) ) + recomposeSB = rsb + where rsb :: SubBasis (DualVectorFromBasis v) + -> [Scalar v] + -> (DualVectorFromBasis v, [Scalar v]) + rsb CompleteDualVBasis cs = first recompose + $ zipWith' (,) (fst <$> enumerate (trie $ const ())) cs + recomposeSBTensor = rsbt + where rsbt :: ∀ w . (FiniteDimensional w, Scalar w ~ Scalar v) + => SubBasis (DualVectorFromBasis v) -> SubBasis w + -> [Scalar v] + -> (DualVectorFromBasis v⊗w, [Scalar v]) + rsbt CompleteDualVBasis sbw ws = proveTensorProductIsTrie @v @w + (first (\iws -> Tensor $ trie (Map.fromList iws Map.!)) + $ zipConsumeWith' (\i cs' -> first (i,) $ recomposeSB sbw cs') + (fst <$> enumerate (trie $ const ())) ws) + recomposeLinMap = rlm + where rlm :: ∀ w . SubBasis (DualVectorFromBasis v) + -> [w] + -> (DualVectorFromBasis v+>w, [w]) + rlm CompleteDualVBasis ws = proveTensorProductIsTrie @v @w + (first (\iws -> LinearMap $ trie (Map.fromList iws Map.!)) + $ zipWith' (,) (fst <$> enumerate (trie $ const ())) ws) + recomposeContraLinMap = rclm + where rclm :: ∀ w f . (LinearSpace w, Scalar w ~ Scalar v, Hask.Functor f) + => (f (Scalar w) -> w) -> f v + -> (DualVectorFromBasis v+>w) + rclm f vs = proveTensorProductIsTrie @v @w + (LinearMap $ trie (\i -> f $ fmap (`decompose'`i) vs)) + recomposeContraLinMapTensor = rclm + where rclm :: ∀ u w f + . ( FiniteDimensional u, LinearSpace w + , Scalar u ~ Scalar v, Scalar w ~ Scalar v, Hask.Functor f + ) + => (f (Scalar w) -> w) -> f (DualVectorFromBasis v+>DualVector u) + -> ((DualVectorFromBasis v⊗u)+>w) + rclm f vus = case dualSpaceWitness @u of + DualSpaceWitness -> proveTensorProductIsTrie @v @(DualVector u) + (proveTensorProductIsTrie @v @(DualVector u⊗w) + (LinearMap $ trie + (\i -> case recomposeContraLinMap @u @w @f f + $ fmap (\(LinearMap vu) -> untrie vu (i :: Basis v)) vus of + LinearMap wuff -> Tensor wuff :: DualVector u⊗w ))) + uncanonicallyFromDual = LinearFunction DualVectorFromBasis + uncanonicallyToDual = LinearFunction getDualVectorFromBasis + + +instance ∀ v . ( BasisGeneratedSpace v, FiniteDimensional v + , Real (Scalar v), Scalar (Scalar v) ~ Scalar v + , HasTrie (Basis v), Ord (Basis v) + , Eq v, Eq (Basis v) ) + => SemiInner (DualVectorFromBasis v) where + dualBasisCandidates = cartesianDualBasisCandidates + (enumerateSubBasis entireBasis) + (\v -> map (abs . realToFrac . decompose' v . fst) + $ enumerate (trie $ const ()) ) + diff --git a/Math/VectorSpace/Docile.hs b/Math/VectorSpace/Docile.hs --- a/Math/VectorSpace/Docile.hs +++ b/Math/VectorSpace/Docile.hs @@ -13,6 +13,7 @@ {-# LANGUAGE TypeOperators #-} {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE GADTs #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE UndecidableInstances #-} @@ -23,6 +24,9 @@ {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE EmptyCase #-} +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE DefaultSignatures #-} module Math.VectorSpace.Docile where @@ -478,7 +482,13 @@ type DList x = [x]->[x] -class (LSpace v) => FiniteDimensional v where + +data DualFinitenessWitness v where + DualFinitenessWitness + :: FiniteDimensional (DualVector v) + => DualSpaceWitness v -> DualFinitenessWitness v + +class (LSpace v, Eq v) => FiniteDimensional v where -- | Whereas 'Basis'-values refer to a single basis vector, a single -- 'SubBasis' value represents a collection of such basis vectors, -- which can be used to associate a vector with a list of coefficients. @@ -528,8 +538,39 @@ -- library). uncanonicallyFromDual :: DualVector v -+> v uncanonicallyToDual :: v -+> DualVector v + + tensorEquality :: (TensorSpace w, Eq w, Scalar w ~ Scalar v) => v⊗w -> v⊗w -> Bool + dualFinitenessWitness :: DualFinitenessWitness v + default dualFinitenessWitness :: FiniteDimensional (DualVector v) + => DualFinitenessWitness v + dualFinitenessWitness = DualFinitenessWitness (dualSpaceWitness @v) + + +instance ( FiniteDimensional u, TensorSpace v + , Scalar u~s, Scalar v~s + , Eq u, Eq v ) => Eq (Tensor s u v) where + (==) = tensorEquality +instance ∀ s u v . ( FiniteDimensional u + , TensorSpace v + , Scalar u~s, Scalar v~s + , Eq v ) + => Eq (LinearMap s u v) where + LinearMap f == LinearMap g = case dualFinitenessWitness @u of + DualFinitenessWitness DualSpaceWitness + -> (Tensor f :: Tensor s (DualVector u) v) == Tensor g + +instance ∀ s u v . ( FiniteDimensional u + , TensorSpace v + , Scalar u~s, Scalar v~s + , Eq v ) + => Eq (LinearFunction s u v) where + f == g = (sampleLinearFunction-+$>f) == (sampleLinearFunction-+$>g) + + + + instance (Num' s) => FiniteDimensional (ZeroDim s) where data SubBasis (ZeroDim s) = ZeroBasis entireBasis = ZeroBasis @@ -544,6 +585,7 @@ recomposeContraLinMapTensor _ _ = LinearMap Origin uncanonicallyFromDual = id uncanonicallyToDual = id + tensorEquality (Tensor Origin) (Tensor Origin) = True instance (Num' s, Eq s, LinearSpace s) => FiniteDimensional (V0 s) where data SubBasis (V0 s) = V0Basis @@ -559,6 +601,7 @@ recomposeContraLinMapTensor _ _ = LinearMap V0 uncanonicallyFromDual = id uncanonicallyToDual = id + tensorEquality (Tensor V0) (Tensor V0) = True instance FiniteDimensional ℝ where data SubBasis ℝ = RealsBasis @@ -576,6 +619,7 @@ . recomposeContraLinMap fw . fmap getLinearMap uncanonicallyFromDual = id uncanonicallyToDual = id + tensorEquality (Tensor v) (Tensor w) = v==w #define FreeFiniteDimensional(V, VB, dimens, take, give) \ instance (Num' s, Eq s, LSpace s) \ @@ -602,7 +646,8 @@ ; rclmt DualSpaceWitness fw mv = LinearMap $ \ (\v -> fromLinearMap $ recomposeContraLinMap fw \ $ fmap (\(LinearMap q) -> foldl' (^+^) zeroV $ liftA2 (*^) v q) mv) \ - <$> Mat.identity } } + <$> Mat.identity }; \ + tensorEquality (Tensor s) (Tensor t) = s==t } FreeFiniteDimensional(V1, V1Basis, 1, c₀ , V1 c₀ ) FreeFiniteDimensional(V2, V2Basis, 2, c₀:c₁ , V2 c₀ c₁ ) FreeFiniteDimensional(V3, V3Basis, 3, c₀:c₁:c₂ , V3 c₀ c₁ c₂ ) @@ -674,6 +719,14 @@ , dualSpaceWitness :: DualSpaceWitness v ) of (DualSpaceWitness,DualSpaceWitness) -> uncanonicallyToDual *** uncanonicallyToDual + tensorEquality (Tensor (s₀,s₁)) (Tensor (t₀,t₁)) + = tensorEquality s₀ t₀ && tensorEquality s₁ t₁ + dualFinitenessWitness = case ( dualFinitenessWitness @u + , dualFinitenessWitness @v ) of + (DualFinitenessWitness DualSpaceWitness + , DualFinitenessWitness DualSpaceWitness) + -> DualFinitenessWitness DualSpaceWitness + deriving instance (Show (SubBasis u), Show (SubBasis v)) => Show (SubBasis (u,v)) @@ -748,7 +801,21 @@ >>> fmap uncanonicallyFromDual >>> transposeTensor >>> fmap uncanonicallyFromDual >>> transposeTensor + tensorEquality = tensTensorEquality + dualFinitenessWitness = case ( dualFinitenessWitness @u + , dualFinitenessWitness @v ) of + (DualFinitenessWitness DualSpaceWitness + , DualFinitenessWitness DualSpaceWitness) + -> DualFinitenessWitness DualSpaceWitness + +tensTensorEquality :: ∀ s u v w . ( FiniteDimensional u, FiniteDimensional v, TensorSpace w + , Scalar u ~ s, Scalar v ~ s, Scalar w ~ s + , Eq w ) + => Tensor s (Tensor s u v) w -> Tensor s (Tensor s u v) w -> Bool +tensTensorEquality (Tensor s) (Tensor t) + = tensorEquality (Tensor s :: Tensor s u (v⊗w)) (Tensor t) + tensorLinmapDecompositionhelpers :: ( FiniteDimensional v, LSpace w , Scalar v~s, Scalar w~s ) => ( [v+>w] -> (SubBasis v, DList w) @@ -776,6 +843,10 @@ deriving instance (Show (SubBasis u), Show (SubBasis v)) => Show (SubBasis (Tensor s u v)) +instance ∀ s v . (FiniteDimensional v, Scalar v ~ s) + => Eq (SymmetricTensor s v) where + SymTensor t == SymTensor u = t==u + instance ∀ s v . ( FiniteDimensional v, Scalar v~s, Scalar (DualVector v)~s , RealFloat' s ) @@ -791,8 +862,9 @@ -- dim Sym(𝑘,𝑉) = nCr (dim 𝑉 + 𝑘 - 1, 𝑘) -- dim Sym(2,𝑉) = nCr (𝑛 + 1, 2) = 𝑛⋅(𝑛+1)/2 where n = subbasisDimension b - decomposeLinMap = dclm dualSpaceWitness - where dclm (DualSpaceWitness :: DualSpaceWitness v) (LinearMap f) + decomposeLinMap = dclm dualFinitenessWitness + where dclm (DualFinitenessWitness DualSpaceWitness :: DualFinitenessWitness v) + (LinearMap f) = (SymTensBasis bf, rmRedundant 0 . symmetrise $ dlw []) where rmRedundant _ [] = id rmRedundant k (row:rest) @@ -819,8 +891,9 @@ oscld = (sqrt 0.5*)<$>o in sd₀ [] ++ [d] ++ oscld ++ mkSym (n-1) (zipWith (.) sds $ (:)<$>oscld) rest - recomposeLinMap = rclm dualSpaceWitness - where rclm (DualSpaceWitness :: DualSpaceWitness v) (SymTensBasis b) ws + recomposeLinMap = rclm dualFinitenessWitness + where rclm (DualFinitenessWitness DualSpaceWitness :: DualFinitenessWitness v) + (SymTensBasis b) ws = case recomposeLinMap (LinMapBasis b b) $ mkSym (subbasisDimension b) (repeat id) ws of (f, remws) -> (LinearMap $ rassocTensor . asTensor $ f, remws) @@ -847,13 +920,14 @@ in concat (sd₀ []) ++ d ++ concat oscld ++ mkSym nw (n-1) (zipWith (.) sds $ (:)<$>oscld) rest recomposeContraLinMap f tenss - = LinearMap . arr (rassocTensor . asTensor) . rcCLM dualSpaceWitness f + = LinearMap . arr (rassocTensor . asTensor) . rcCLM dualFinitenessWitness f $ fmap getSymmetricTensor tenss where rcCLM :: (Hask.Functor f, LinearSpace w, s~Scalar w) - => DualSpaceWitness v + => DualFinitenessWitness v -> (f s->w) -> f (Tensor s (DualVector v) (DualVector v)) -> LinearMap s (LinearMap s (DualVector v) v) w - rcCLM DualSpaceWitness f = recomposeContraLinMap f + rcCLM (DualFinitenessWitness DualSpaceWitness) f + = recomposeContraLinMap f recomposeContraLinMapTensor = rcCLMT' where rcCLMT' :: ∀ f u w . (Hask.Functor f, LinearSpace w, s~Scalar w , FiniteDimensional u, s~Scalar u) @@ -861,9 +935,9 @@ -> (SymmetricTensor s v ⊗ u) +> w rcCLMT' f tenss = LinearMap . arr (fmap rassocTensor . rassocTensor . asTensor) - . rcCLMT (dualSpaceWitness, dualSpaceWitness) f + . rcCLMT (dualFinitenessWitness, dualFinitenessWitness) f $ fmap getLinearMap tenss - where rcCLMT :: (DualSpaceWitness v, DualSpaceWitness u) + where rcCLMT :: (DualFinitenessWitness v, DualFinitenessWitness u) -> (f s->w) -> f (Tensor s (DualVector v) (Tensor s (DualVector v) (DualVector u))) -- -> LinearMap s (Tensor s (SymmetricTensor s v) u) w @@ -876,86 +950,110 @@ -- (DualVector v ⊗ DualVector u)) w -- ∼ Tensor s (DualVector v) -- (Tensor s (DualVector v ⊗ DualVector u) w) - rcCLMT (DualSpaceWitness, DualSpaceWitness) f = recomposeContraLinMap f - uncanonicallyFromDual = case dualSpaceWitness :: DualSpaceWitness v of - DualSpaceWitness -> LinearFunction + rcCLMT ( DualFinitenessWitness DualSpaceWitness + , DualFinitenessWitness DualSpaceWitness ) f + = recomposeContraLinMap f + uncanonicallyFromDual = case dualFinitenessWitness :: DualFinitenessWitness v of + DualFinitenessWitness DualSpaceWitness -> LinearFunction $ \(SymTensor t) -> SymTensor $ arr fromLinearMap . uncanonicallyFromDual $ t - uncanonicallyToDual = case dualSpaceWitness :: DualSpaceWitness v of - DualSpaceWitness -> LinearFunction + uncanonicallyToDual = case dualFinitenessWitness :: DualFinitenessWitness v of + DualFinitenessWitness DualSpaceWitness -> LinearFunction $ \(SymTensor t) -> SymTensor $ uncanonicallyToDual . arr asLinearMap $ t + dualFinitenessWitness = case dualFinitenessWitness @v of + DualFinitenessWitness DualSpaceWitness + -> DualFinitenessWitness DualSpaceWitness deriving instance (Show (SubBasis v)) => Show (SubBasis (SymmetricTensor s v)) instance ∀ s u v . - ( LSpace u, FiniteDimensional (DualVector u), FiniteDimensional v + ( LSpace u, FiniteDimensional u, FiniteDimensional v , Scalar u~s, Scalar v~s, Scalar (DualVector v)~s, Fractional' (Scalar v) ) => FiniteDimensional (LinearMap s u v) where data SubBasis (LinearMap s u v) = LinMapBasis !(SubBasis (DualVector u)) !(SubBasis v) - entireBasis = case ( dualSpaceWitness :: DualSpaceWitness u + entireBasis = case ( dualFinitenessWitness :: DualFinitenessWitness u , dualSpaceWitness :: DualSpaceWitness v ) of - (DualSpaceWitness, DualSpaceWitness) + (DualFinitenessWitness DualSpaceWitness, DualSpaceWitness) -> case entireBasis of TensorBasis bu bv -> LinMapBasis bu bv enumerateSubBasis - = case ( dualSpaceWitness :: DualSpaceWitness u + = case ( dualFinitenessWitness :: DualFinitenessWitness u , dualSpaceWitness :: DualSpaceWitness v ) of - (DualSpaceWitness, DualSpaceWitness) -> \(LinMapBasis bu bv) + (DualFinitenessWitness DualSpaceWitness, DualSpaceWitness) -> \(LinMapBasis bu bv) -> arr (fmap asLinearMap) . enumerateSubBasis $ TensorBasis bu bv - subbasisDimension (LinMapBasis bu bv) = subbasisDimension bu * subbasisDimension bv - decomposeLinMap = case ( dualSpaceWitness :: DualSpaceWitness u + subbasisDimension (LinMapBasis bu bv) + = case ( dualFinitenessWitness :: DualFinitenessWitness u ) of + (DualFinitenessWitness _) -> subbasisDimension bu * subbasisDimension bv + decomposeLinMap = case ( dualFinitenessWitness :: DualFinitenessWitness u , dualSpaceWitness :: DualSpaceWitness v ) of - (DualSpaceWitness, DualSpaceWitness) + (DualFinitenessWitness DualSpaceWitness, DualSpaceWitness) -> first (\(TensorBasis bu bv)->LinMapBasis bu bv) . decomposeLinMap . coerce - decomposeLinMapWithin = case ( dualSpaceWitness :: DualSpaceWitness u + decomposeLinMapWithin = case ( dualFinitenessWitness :: DualFinitenessWitness u , dualSpaceWitness :: DualSpaceWitness v ) of - (DualSpaceWitness, DualSpaceWitness) + (DualFinitenessWitness DualSpaceWitness, DualSpaceWitness) -> \(LinMapBasis bu bv) m -> case decomposeLinMapWithin (TensorBasis bu bv) (coerce m) of Right ws -> Right ws Left (TensorBasis bu' bv', ws) -> Left (LinMapBasis bu' bv', ws) - recomposeSB = case ( dualSpaceWitness :: DualSpaceWitness u + recomposeSB = case ( dualFinitenessWitness :: DualFinitenessWitness u , dualSpaceWitness :: DualSpaceWitness v ) of - (DualSpaceWitness, DualSpaceWitness) -> \(LinMapBasis bu bv) + (DualFinitenessWitness DualSpaceWitness, DualSpaceWitness) -> \(LinMapBasis bu bv) -> recomposeSB (TensorBasis bu bv) >>> first (arr fromTensor) - recomposeSBTensor = case ( dualSpaceWitness :: DualSpaceWitness u + recomposeSBTensor = case ( dualFinitenessWitness :: DualFinitenessWitness u , dualSpaceWitness :: DualSpaceWitness v ) of - (DualSpaceWitness, DualSpaceWitness) -> \(LinMapBasis bu bv) bw + (DualFinitenessWitness DualSpaceWitness, DualSpaceWitness) -> \(LinMapBasis bu bv) bw -> recomposeSBTensor (TensorBasis bu bv) bw >>> first coerce - recomposeLinMap = rlm dualSpaceWitness dualSpaceWitness + recomposeLinMap = rlm dualFinitenessWitness dualSpaceWitness where rlm :: ∀ w . (LSpace w, Scalar w ~ Scalar v) - => DualSpaceWitness u -> DualSpaceWitness w -> SubBasis (u+>v) -> [w] + => DualFinitenessWitness u -> DualSpaceWitness w -> SubBasis (u+>v) -> [w] -> ((u+>v)+>w, [w]) - rlm DualSpaceWitness DualSpaceWitness (LinMapBasis bu bv) ws + rlm (DualFinitenessWitness DualSpaceWitness) DualSpaceWitness (LinMapBasis bu bv) ws = ( coUncurryLinearMap . fromLinearMap $ fst . recomposeLinMap bu $ unfoldr (pure . recomposeLinMap bv) ws , drop (subbasisDimension bu * subbasisDimension bv) ws ) - recomposeContraLinMap = case ( dualSpaceWitness :: DualSpaceWitness u + recomposeContraLinMap = case ( dualFinitenessWitness :: DualFinitenessWitness u , dualSpaceWitness :: DualSpaceWitness v ) of - (DualSpaceWitness, DualSpaceWitness) -> \fw dds + (DualFinitenessWitness DualSpaceWitness, DualSpaceWitness) -> \fw dds -> argFromTensor $ recomposeContraLinMapTensor fw $ fmap (arr asLinearMap) dds - recomposeContraLinMapTensor = rclmt dualSpaceWitness dualSpaceWitness dualSpaceWitness + recomposeContraLinMapTensor = rclmt dualFinitenessWitness dualSpaceWitness dualSpaceWitness where rclmt :: ∀ f u' w . ( LinearSpace w, FiniteDimensional u' , Scalar w ~ s, Scalar u' ~ s , Hask.Functor f ) - => DualSpaceWitness u -> DualSpaceWitness v -> DualSpaceWitness u' + => DualFinitenessWitness u -> DualSpaceWitness v -> DualSpaceWitness u' -> (f (Scalar w) -> w) -> f ((u+>v)+>DualVector u') -> ((u+>v)⊗u')+>w - rclmt DualSpaceWitness DualSpaceWitness DualSpaceWitness fw dds + rclmt (DualFinitenessWitness DualSpaceWitness) + DualSpaceWitness DualSpaceWitness fw dds = uncurryLinearMap . coUncurryLinearMap . fmap curryLinearMap . coCurryLinearMap . argFromTensor $ recomposeContraLinMapTensor fw $ fmap (arr $ asLinearMap . coCurryLinearMap) dds - uncanonicallyToDual = case ( dualSpaceWitness :: DualSpaceWitness u + uncanonicallyToDual = case ( dualFinitenessWitness :: DualFinitenessWitness u , dualSpaceWitness :: DualSpaceWitness v ) of - (DualSpaceWitness, DualSpaceWitness) + (DualFinitenessWitness DualSpaceWitness, DualSpaceWitness) -> arr asTensor >>> fmap uncanonicallyToDual >>> transposeTensor >>> fmap uncanonicallyToDual >>> transposeTensor - uncanonicallyFromDual = case ( dualSpaceWitness :: DualSpaceWitness u + uncanonicallyFromDual = case ( dualFinitenessWitness :: DualFinitenessWitness u , dualSpaceWitness :: DualSpaceWitness v ) of - (DualSpaceWitness, DualSpaceWitness) + (DualFinitenessWitness DualSpaceWitness, DualSpaceWitness) -> arr fromTensor <<< fmap uncanonicallyFromDual <<< transposeTensor <<< fmap uncanonicallyFromDual <<< transposeTensor - + + tensorEquality = lmTensorEquality + dualFinitenessWitness = case ( dualFinitenessWitness @u + , dualFinitenessWitness @v ) of + (DualFinitenessWitness DualSpaceWitness + , DualFinitenessWitness DualSpaceWitness) + -> DualFinitenessWitness DualSpaceWitness + +lmTensorEquality :: ∀ s u v w . ( FiniteDimensional v, TensorSpace w + , FiniteDimensional u + , Scalar u ~ s, Scalar v ~ s, Scalar w ~ s + , Eq w ) + => Tensor s (LinearMap s u v) w -> Tensor s (LinearMap s u v) w -> Bool +lmTensorEquality (Tensor s) (Tensor t) = case dualFinitenessWitness @u of + DualFinitenessWitness DualSpaceWitness + -> tensorEquality (Tensor s :: Tensor s (DualVector u) (v⊗w)) (Tensor t) + deriving instance (Show (SubBasis (DualVector u)), Show (SubBasis v)) => Show (SubBasis (LinearMap s u v)) @@ -1038,7 +1136,7 @@ where m' = adjoint $ m :: DualVector v +> DualVector u -- | Invert an isomorphism. For other linear maps, the result is undefined. -unsafeInverse :: ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v ) +unsafeInverse :: ( FiniteDimensional u, SimpleSpace v, Scalar u ~ Scalar v ) => (u+>v) -> v+>u unsafeInverse m = recomposeContraLinMap (fst . recomposeSB mbas) $ [maybe zeroV id v' | v'<-v's] @@ -1048,12 +1146,12 @@ -- | The -- provides an isomorphism between a Hilbert space and its (continuous) dual space. -riesz :: ∀ v . (FiniteDimensional v, InnerSpace v) => DualVector v -+> v -riesz = case ( scalarSpaceWitness :: ScalarSpaceWitness v - , dualSpaceWitness :: DualSpaceWitness v ) of - (ScalarSpaceWitness,DualSpaceWitness) -> LinearFunction $ \dv -> - let (bas, compos) = decomposeLinMap $ sampleLinearFunction $ applyDualVector $ dv - in fst . recomposeSB bas $ compos [] +riesz :: ∀ v . ( FiniteDimensional v, InnerSpace v + , SimpleSpace v ) + => DualVector v -+> v +riesz = case dualFinitenessWitness @v of + DualFinitenessWitness DualSpaceWitness + -> arr . unsafeInverse $ arr coRiesz sRiesz :: ∀ v . FiniteDimensional v => DualSpace v -+> v sRiesz = case ( scalarSpaceWitness :: ScalarSpaceWitness v @@ -1138,6 +1236,22 @@ bw .< v = sampleLinearFunction $ LinearFunction $ \v' -> recompose [(bw, v<.>v')] +-- | This is the preferred method for showing linear maps, resulting in a +-- matrix view involving the '.<' operator. +-- We don't provide a generic `Show` instance; to make linear maps with +-- your own finite-dimensional type @V@ (with scalar @S@) showable, +-- this is the recommended way: +-- +-- @ +-- instance RieszDecomposable V where +-- rieszDecomposition = ... +-- instance (FiniteDimensional w, w ~ DualVector w, Scalar w ~ S, Show w) +-- => Show (LinearMap S w V) where +-- showsPrec = rieszDecomposeShowsPrec +-- @ +-- +-- Note that the custom type should always be the /codomain/ type, whereas +-- the domain should be kept parametric. rieszDecomposeShowsPrec :: ∀ u v s . ( RieszDecomposable u , FiniteDimensional v, v ~ DualVector v, Show v , Scalar u ~ s, Scalar v ~ s ) @@ -1145,10 +1259,10 @@ rieszDecomposeShowsPrec p m = case rieszDecomposition m of [] -> ("zeroV"++) ((b₀,dv₀):dvs) -> showParen (p>6) - $ \s -> showsPrecBasis ([]::[u]) 7 b₀ + $ \s -> showsPrecBasis @u 7 b₀ . (".<"++) . showsPrec 7 dv₀ $ foldr (\(b,dv) - -> (" ^+^ "++) . showsPrecBasis ([]::[u]) 7 b + -> (" ^+^ "++) . showsPrecBasis @u 7 b . (".<"++) . showsPrec 7 dv) s dvs instance Show (LinearMap s v (ZeroDim s)) where @@ -1187,40 +1301,40 @@ class (FiniteDimensional v, HasBasis v) => TensorDecomposable v where tensorDecomposition :: v⊗w -> [(Basis v, w)] - showsPrecBasis :: Hask.Functor p => p v -> Int -> Basis v -> ShowS + showsPrecBasis :: Int -> Basis v -> ShowS instance TensorDecomposable ℝ where tensorDecomposition (Tensor r) = [((), r)] - showsPrecBasis _ _ = shows -instance ( TensorDecomposable x, TensorDecomposable y - , Scalar x ~ Scalar y, Scalar (DualVector x) ~ Scalar (DualVector y) ) + showsPrecBasis _ = shows +instance ∀ x y . ( TensorDecomposable x, TensorDecomposable y + , Scalar x ~ Scalar y, Scalar (DualVector x) ~ Scalar (DualVector y) ) => TensorDecomposable (x,y) where tensorDecomposition (Tensor (tx,ty)) = map (first Left) (tensorDecomposition tx) ++ map (first Right) (tensorDecomposition ty) - showsPrecBasis proxy p (Left bx) - = showParen (p>9) $ ("Left "++) . showsPrecBasis (fst<$>proxy) 10 bx - showsPrecBasis proxy p (Right by) - = showParen (p>9) $ ("Right "++) . showsPrecBasis (snd<$>proxy) 10 by + showsPrecBasis p (Left bx) + = showParen (p>9) $ ("Left "++) . showsPrecBasis @x 10 bx + showsPrecBasis p (Right by) + = showParen (p>9) $ ("Right "++) . showsPrecBasis @y 10 by instance TensorDecomposable (ZeroDim ℝ) where tensorDecomposition _ = [] - showsPrecBasis _ _ = absurd + showsPrecBasis _ = absurd instance TensorDecomposable (V0 ℝ) where tensorDecomposition _ = [] - showsPrecBasis _ _ (Mat.E q) = (V0^.q ++) + showsPrecBasis _ (Mat.E q) = (V0^.q ++) instance TensorDecomposable (V1 ℝ) where tensorDecomposition (Tensor (V1 w)) = [(ex, w)] - showsPrecBasis _ _ (Mat.E q) = (V1"ex"^.q ++) + showsPrecBasis _ (Mat.E q) = (V1"ex"^.q ++) instance TensorDecomposable (V2 ℝ) where tensorDecomposition (Tensor (V2 x y)) = [ (ex, x), (ey, y) ] - showsPrecBasis _ _ (Mat.E q) = (V2"ex""ey"^.q ++) + showsPrecBasis _ (Mat.E q) = (V2"ex""ey"^.q ++) instance TensorDecomposable (V3 ℝ) where tensorDecomposition (Tensor (V3 x y z)) = [ (ex, x), (ey, y), (ez, z) ] - showsPrecBasis _ _ (Mat.E q) = (V3"ex""ey""ez"^.q ++) + showsPrecBasis _ (Mat.E q) = (V3"ex""ey""ez"^.q ++) instance TensorDecomposable (V4 ℝ) where tensorDecomposition (Tensor (V4 x y z w)) = [ (ex, x), (ey, y), (ez, z), (ew, w) ] - showsPrecBasis _ _ (Mat.E q) = (V4"ex""ey""ez""ew"^.q ++) + showsPrecBasis _ (Mat.E q) = (V4"ex""ey""ez""ew"^.q ++) tensorDecomposeShowsPrec :: ∀ u v s . ( TensorDecomposable u, FiniteDimensional v, Show v, Scalar u ~ s, Scalar v ~ s ) @@ -1228,10 +1342,10 @@ tensorDecomposeShowsPrec p t = case tensorDecomposition t of [] -> ("zeroV"++) ((b₀,dv₀):dvs) -> showParen (p>6) - $ \s -> showsPrecBasis ([]::[u]) 7 b₀ + $ \s -> showsPrecBasis @u 7 b₀ . (".⊗"++) . showsPrec 7 dv₀ $ foldr (\(b,dv) - -> (" ^+^ "++) . showsPrecBasis ([]::[u]) 7 b + -> (" ^+^ "++) . showsPrecBasis @u 7 b . (".⊗"++) . showsPrec 7 dv) s dvs instance Show (Tensor s (V0 s) v) where diff --git a/Math/VectorSpace/Dual.hs b/Math/VectorSpace/Dual.hs new file mode 100644 --- /dev/null +++ b/Math/VectorSpace/Dual.hs @@ -0,0 +1,64 @@ +-- | +-- Module : Math.VectorSpace.Dual +-- Copyright : (c) Justus Sagemüller 2020 +-- License : GPL v3 +-- +-- Maintainer : (@) jsag $ hvl.no +-- Stability : experimental +-- Portability : portable +-- + +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE KindSignatures #-} + +module Math.VectorSpace.Dual + ( Dualness(..), Dual, DualityWitness(..), ValidDualness(..) + , usingAnyDualness + ) where + +import Math.LinearMap.Category.Class +import Data.Kind (Type) + + +data Dualness = Vector | Functional + +type family Dual (dn :: Dualness) where + Dual Vector = Functional + Dual Functional = Vector + + +data DualityWitness (dn :: Dualness) where + DualityWitness :: (ValidDualness (Dual dn), Dual (Dual dn) ~ dn) + => DualityWitness dn + +data DualnessSingletons (dn :: Dualness) where + VectorWitness :: DualnessSingletons Vector + FunctionalWitness :: DualnessSingletons Functional + +class ValidDualness (dn :: Dualness) where + type Space dn v :: Type + dualityWitness :: DualityWitness dn + decideDualness :: DualnessSingletons dn +instance ValidDualness 'Vector where + type 'Vector`Space` v = v + dualityWitness = DualityWitness + decideDualness = VectorWitness +instance ValidDualness 'Functional where + type 'Functional`Space` v = DualVector v + dualityWitness = DualityWitness + decideDualness = FunctionalWitness + +usingAnyDualness :: ∀ rc dn . ValidDualness dn + => (rc 'Vector) + -> (rc 'Functional) + -> rc dn +usingAnyDualness vecCase ftnCase = case decideDualness @dn of + VectorWitness -> vecCase + FunctionalWitness -> ftnCase diff --git a/Math/VectorSpace/MiscUtil/MultiConstraints.hs b/Math/VectorSpace/MiscUtil/MultiConstraints.hs new file mode 100644 --- /dev/null +++ b/Math/VectorSpace/MiscUtil/MultiConstraints.hs @@ -0,0 +1,35 @@ +-- | +-- Module : Math.VectorSpace.MiscUtil.MultiConstraints +-- Copyright : (c) Justus Sagemüller 2020 +-- License : GPL v3 +-- +-- Maintainer : (@) jsag $ hvl.no +-- Stability : experimental +-- Portability : portable +-- + +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE UndecidableInstances #-} + +module Math.VectorSpace.MiscUtil.MultiConstraints + ( SameScalar + ) where + +import Data.VectorSpace +import Data.Kind (Type) +import GHC.Exts (Constraint) + + +type family AllWithScalar (s :: Type) (c :: Type -> Constraint) (vs :: [Type]) + :: Constraint where + AllWithScalar s c '[] = () + AllWithScalar s c (v ': vs) = (c v, Scalar v ~ s, AllWithScalar s c vs) + +type family SameScalar (c :: Type -> Constraint) (vs :: [Type]) :: Constraint where + SameScalar c '[] = () + SameScalar c (v ': vs) + = (c v, AllWithScalar (Scalar v) c vs) diff --git a/linearmap-category.cabal b/linearmap-category.cabal --- a/linearmap-category.cabal +++ b/linearmap-category.cabal @@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: linearmap-category -version: 0.4.0.1 +version: 0.4.1.0 synopsis: Native, complete, matrix-free linear algebra. description: The term /numerical linear algebra/ is often used almost synonymous with /matrix modifications/. However, what's interesting @@ -39,7 +39,10 @@ library exposed-modules: Math.LinearMap.Category + Math.LinearMap.Category.Instances.Deriving Math.VectorSpace.ZeroDimensional + Math.VectorSpace.Dual + Math.VectorSpace.MiscUtil.MultiConstraints Math.LinearMap.Category.Derivatives other-modules: Math.LinearMap.Category.Class Math.LinearMap.Asserted @@ -48,16 +51,31 @@ Math.VectorSpace.Docile other-extensions: FlexibleInstances, UndecidableInstances, FunctionalDependencies, TypeOperators, TypeFamilies build-depends: base >=4.8 && <5, - vector-space >=0.11 && <0.18, + vector-space >=0.11 && <0.18, MemoTrie, constrained-categories >=0.3 && <0.5, containers, vector, tagged, free-vector-spaces >= 0.1.4 && < 0.2, linear, lens, transformers, - manifolds-core >= 0.5.1.0 && < 0.6, - semigroups, + manifolds-core >= 0.5.1.0 && < 0.7, + semigroups, hashable, + data-default-class, ieee754 >= 0.7 && < 0.9, - call-stack, + call-stack, template-haskell >=2.13 && <2.15, QuickCheck >=2.11 && <2.15 -- hs-source-dirs: default-language: Haskell2010 + +test-suite tasty + main-is: test.hs + type: exitcode-stdio-1.0 + hs-source-dirs: test/tasty + build-depends: base, linearmap-category, vector-space + , QuickCheck + , manifolds-core + , linear + , constrained-categories + , tasty, tasty-quickcheck + ghc-options: -threaded "-with-rtsopts -N8 -M2G" + default-language: Haskell2010 + diff --git a/test/tasty/test.hs b/test/tasty/test.hs new file mode 100644 --- /dev/null +++ b/test/tasty/test.hs @@ -0,0 +1,123 @@ +-- | +-- Module : test +-- Copyright : (c) Justus Sagemüller 2021 +-- License : GPL v3 +-- +-- Maintainer : (@) jsag $ hvl.no +-- Stability : experimental +-- Portability : portable +-- + +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE GeneralisedNewtypeDeriving #-} +{-# LANGUAGE DerivingStrategies #-} + +import qualified Prelude as Hask +import Control.Category.Constrained.Prelude +import Control.Arrow.Constrained + +import Data.AffineSpace +import Linear.V4 +import Data.Basis +import Math.LinearMap.Category +import Math.Manifold.Core.Types +import Math.Manifold.Core.PseudoAffine + +import Test.Tasty +import Test.Tasty.QuickCheck +import qualified Test.QuickCheck as QC + + + +newtype ℝ⁵ a = ℝ⁵ { getℝ⁵ :: [ℝ] } + deriving (Eq, Show) + +instance AdditiveGroup (ℝ⁵ a) where + zeroV = ℝ⁵ $ replicate 5 0 + ℝ⁵ v ^+^ ℝ⁵ w = ℝ⁵ $ zipWith (+) v w + negateV (ℝ⁵ v) = ℝ⁵ $ map negate v + +instance VectorSpace (ℝ⁵ a) where + type Scalar (ℝ⁵ a) = ℝ + μ*^ℝ⁵ v = ℝ⁵ $ map (μ*) v + +instance InnerSpace (ℝ⁵ a) where + ℝ⁵ v <.> ℝ⁵ w = sum $ zipWith (*) v w + +type Z5 = Z3+Z2 +type Z3 = Z2+() +type Z2 = ()+() + +instance Num a => HasBasis (ℝ⁵ a) where + type Basis (ℝ⁵ a) = Z5 + basisValue (Left (Left (Left ()))) = ℝ⁵ [1,0,0,0,0] + basisValue (Left (Left (Right ()))) = ℝ⁵ [0,1,0,0,0] + basisValue (Left (Right () )) = ℝ⁵ [0,0,1,0,0] + basisValue (Right (Left () )) = ℝ⁵ [0,0,0,1,0] + basisValue (Right (Right () )) = ℝ⁵ [0,0,0,0,1] + decompose (ℝ⁵ [a,b,c,d,e]) = + [ (Left (Left (Left ())), a) + , (Left (Left (Right ())), b) + , (Left (Right () ), c) + , (Right (Left () ), d) + , (Right (Right () ), e) ] + decompose' (ℝ⁵ [a,b,c,d,e]) n = case n of + Left (Left (Left ())) -> a + Left (Left (Right ())) -> b + Left (Right () ) -> c + Right (Left () ) -> d + Right (Right () ) -> e + +instance Arbitrary (ℝ⁵ a) where + arbitrary = ℝ⁵ <$> QC.vectorOf 5 arbitrary + +makeFiniteDimensionalFromBasis [t| ∀ a . Num a => ℝ⁵ a |] + +newtype H¹ℝ⁵ = H¹ℝ⁵ { getH¹ℝ⁵ :: ℝ⁵ Int } + deriving newtype (Eq, Show, AdditiveGroup, VectorSpace, HasBasis, Arbitrary) + +makeFiniteDimensionalFromBasis [t| H¹ℝ⁵ |] + +derivative :: H¹ℝ⁵ -> ℝ⁵ Int +derivative (H¹ℝ⁵ (ℝ⁵ (x₀:xs))) = ℝ⁵ (x₀:xs) ^-^ ℝ⁵ (xs++[x₀]) + +instance InnerSpace H¹ℝ⁵ where + H¹ℝ⁵ v <.> H¹ℝ⁵ w = v<.>w + derivative (H¹ℝ⁵ v)<.>derivative (H¹ℝ⁵ w) + +instance Arbitrary (V4 ℝ) where + arbitrary = V4<$>arbitrary<*>arbitrary<*>arbitrary<*>arbitrary + + + +main :: IO () +main = do + defaultMain $ testGroup "Tests" + [ testGroup "Euclidean space" + [ testProperty "co-Riesz inversion" + $ \v -> (arr coRiesz\$coRiesz-+$>v) === (v :: V4 ℝ) + , testProperty "Random operator inversion" -- This isn't really expected to work + $ \f v -> (f \$ (f :: V4 ℝ+>V4 ℝ) $ v) ≈≈≈ v -- /always/, but singular matrices are + ] -- very seldom in the @Arbitrary@ instance. + , testGroup "Basis-derived space" + [ testProperty "Semimanifold addition" + $ \v w -> v.+~^w === (v^+^w :: ℝ⁵ Int) + , testProperty "Riesz representation, orthonormal basis" + $ \v -> (riesz-+$>coRiesz-+$>v) === (v :: ℝ⁵ Int) + , testProperty "Riesz representation, non-orthonormal basis" + $ \v -> (riesz-+$>coRiesz-+$>v) ≈≈≈ (v :: H¹ℝ⁵) + ] + ] + + +(≈≈≈) :: (InnerSpace v, Show v, Eq v, RealFrac (Scalar v)) + => v -> v -> QC.Property +v≈≈≈w + | magnitudeSq (v^-^w) < (magnitudeSq v + magnitudeSq w)*1e-8 = QC.property True + | otherwise = v===w