diff --git a/Math/LinearMap/Asserted.hs b/Math/LinearMap/Asserted.hs --- a/Math/LinearMap/Asserted.hs +++ b/Math/LinearMap/Asserted.hs @@ -94,7 +94,7 @@ instance Functor (LinearFunction s v) Coercion Coercion where fmap Coercion = Coercion -fmapScale :: ( VectorSpace w, Scalar w ~ s, VectorSpace s, Scalar s ~ s +fmapScale :: ( VectorSpace w, Scalar w ~ s, VectorSpace s , Functor f (LinearFunction s) (LinearFunction s) , Object (LinearFunction s) s , Object (LinearFunction s) w @@ -145,3 +145,7 @@ lApply :: Bilinear (v-+>w) v w lApply = bilinearFunction $ \(LinearFunction f) v -> f v + +infixr 0 -+$> +(-+$>) :: LinearFunction s v w -> v -> w +LinearFunction f -+$> v = f v diff --git a/Math/LinearMap/Category.hs b/Math/LinearMap/Category.hs --- a/Math/LinearMap/Category.hs +++ b/Math/LinearMap/Category.hs @@ -20,22 +20,23 @@ {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE TupleSections #-} {-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE ExplicitNamespaces #-} module Math.LinearMap.Category ( -- * Linear maps -- $linmapIntro -- ** Function implementation - LinearFunction (..), (-+>)(), Bilinear + LinearFunction (..), type (-+>)(), Bilinear -- ** Tensor implementation - , LinearMap (..), (+>)() + , LinearMap (..), type (+>)() , (⊕), (>+<) , adjoint -- ** Dual vectors -- $dualVectorIntro - , (<.>^) + , (<.>^), (-+|>) -- * Tensor spaces - , Tensor (..), (⊗)(), (⊗) + , Tensor (..), type (⊗)(), (⊗) -- * Norms -- $metricIntro , Norm(..), Seminorm @@ -71,17 +72,21 @@ -- * Utility -- ** Linear primitives , addV, scale, inner, flipBilin, bilinearFunction + -- ** Tensors with basis decomposition + , (.⊗) -- ** Hilbert space operations , DualSpace, riesz, coRiesz, showsPrecAsRiesz, (.<) -- ** Constraint synonyms , HilbertSpace, SimpleSpace - , Num', Num'', Num''' - , Fractional', Fractional'' + , Num' + , Fractional' , RealFrac', RealFloat' -- ** Misc , relaxNorm, transformNorm, transformVariance , findNormalLength, normalLength - , summandSpaceNorms, sumSubspaceNorms, sharedNormSpanningSystem + , summandSpaceNorms, sumSubspaceNorms + , sharedNormSpanningSystem, sharedSeminormSpanningSystem + , sharedSeminormSpanningSystem' ) where import Math.LinearMap.Category.Class @@ -95,6 +100,7 @@ import Data.Set (Set) import Data.Ord (comparing) import Data.List (maximumBy) +import Data.Maybe (catMaybes) import Data.Foldable (toList) import Data.Semigroup @@ -194,6 +200,17 @@ +-- | A linear map that simply projects from a dual vector in @u@ to a vector in @v@. +-- +-- @ +-- (du-+|>v) u ≡ v ^* (du<.>^u) +-- @ +infixr 7 -+|> +(-+|>) :: ( EnhancedCat f (LinearFunction s) + , LSpace u, LSpace v, Scalar u ~ s, Scalar v ~ s + , Object f u, Object f v ) + => DualVector u -> v -> f u v +du-+|>v = arr . LinearFunction $ (v^*) . (du<.>^) @@ -221,17 +238,22 @@ -- -- If the @dᵢ@ are a complete orthonormal system, you get the 'euclideanNorm' -- (in an inefficient form). -spanNorm :: LSpace v => [DualVector v] -> Seminorm v -spanNorm dvs = Norm . LinearFunction $ \v -> sumV [dv ^* (dv<.>^v) | dv <- dvs] +spanNorm :: ∀ v . LSpace v => [DualVector v] -> Seminorm v +spanNorm = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness + -> \dvs -> Norm . LinearFunction $ \v -> sumV [dv ^* (dv<.>^v) | dv <- dvs] -spanVariance :: LSpace v => [v] -> Variance v -spanVariance = spanNorm +spanVariance :: ∀ v . LSpace v => [v] -> Variance v +spanVariance = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> spanNorm -- | Modify a norm in such a way that the given vectors lie within its unit ball. -- (Not /optimally/ – the unit ball may be bigger than necessary.) -relaxNorm :: SimpleSpace v => Norm v -> [v] -> Norm v -relaxNorm me = \vs -> dualNorm . spanVariance $ vs' ++ vs - where vs' = normSpanningSystem' me +relaxNorm :: ∀ v . SimpleSpace v => Norm v -> [v] -> Norm v +relaxNorm = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness + -> \me vs -> let vs' = normSpanningSystem' me + in dualNorm . spanVariance $ vs' ++ vs -- | Scale the result of a norm with the absolute of the given number. -- @@ -240,8 +262,9 @@ -- @ -- -- Equivalently, this scales the norm's unit ball by the reciprocal of that factor. -scaleNorm :: LSpace v => Scalar v -> Norm v -> Norm v -scaleNorm μ (Norm n) = Norm $ μ^2 *^ n +scaleNorm :: ∀ v . LSpace v => Scalar v -> Norm v -> Norm v +scaleNorm = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> \μ (Norm n) -> Norm $ μ^2 *^ n -- | A positive (semi)definite symmetric bilinear form. This gives rise -- to a thus: @@ -291,32 +314,41 @@ -- (The orthonormal systems of the norm and its dual are mutually conjugate.) -- The dual norm of a seminorm is undefined. dualNorm :: SimpleSpace v => Norm v -> Variance v -dualNorm (Norm m) = Norm . arr . pseudoInverse $ arr m +dualNorm = spanVariance . normSpanningSystem' -transformNorm :: (LSpace v, LSpace w, Scalar v~Scalar w) => (v+>w) -> Norm w -> Norm v -transformNorm f (Norm m) = Norm . arr $ (adjoint $ f) . (fmap m $ f) +transformNorm :: ∀ v w . (LSpace v, LSpace w, Scalar v~Scalar w) + => (v+>w) -> Norm w -> Norm v +transformNorm = case ( dualSpaceWitness :: DualSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness w ) of + (DualSpaceWitness, DualSpaceWitness) + -> \f (Norm m) -> Norm . arr $ (adjoint $ f) . (fmap m $ f) -transformVariance :: (LSpace v, LSpace w, Scalar v~Scalar w) +transformVariance :: ∀ v w . (LSpace v, LSpace w, Scalar v~Scalar w) => (v+>w) -> Variance v -> Variance w -transformVariance f (Norm m) = Norm . arr $ f . (fmap m $ adjoint $ f) +transformVariance = case ( dualSpaceWitness :: DualSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness w ) of + (DualSpaceWitness, DualSpaceWitness) + -> \f (Norm m) -> Norm . arr $ f . (fmap m $ adjoint $ f) infixl 6 ^% (^%) :: (LSpace v, Floating (Scalar v)) => v -> Norm v -> v -v ^% Norm m = v ^/ sqrt ((m$v)<.>^v) +v ^% Norm m = v ^/ sqrt ((m-+$>v)<.>^v) -- | The unique positive number whose norm is 1 (if the norm is not constant zero). -findNormalLength :: RealFrac' s => Norm s -> Maybe s -findNormalLength (Norm m) = case m $ 1 of - o | o > 0 -> Just . sqrt $ recip o - | otherwise -> Nothing +findNormalLength :: ∀ s . RealFrac' s => Norm s -> Maybe s +findNormalLength (Norm m) = case ( closedScalarWitness :: ClosedScalarWitness s + , m-+$>1 ) of + (ClosedScalarWitness, o) | o > 0 -> Just . sqrt $ recip o + | otherwise -> Nothing -- | Unsafe version of 'findNormalLength', only works reliable if the norm -- is actually positive definite. -normalLength :: RealFrac' s => Norm s -> s -normalLength (Norm m) = case m $ 1 of - o | o >= 0 -> sqrt $ recip o - | o < 0 -> error "Norm fails to be positive semidefinite." - | otherwise -> error "Norm yields NaN." +normalLength :: ∀ s . RealFrac' s => Norm s -> s +normalLength (Norm m) = case ( closedScalarWitness :: ClosedScalarWitness s + , m-+$>1 ) of + (ClosedScalarWitness, o) | o >= 0 -> sqrt $ recip o + | o < 0 -> error "Norm fails to be positive semidefinite." + | otherwise -> error "Norm yields NaN." infixr 0 <$|, |$| -- | “Partially apply” a norm, yielding a dual vector @@ -326,12 +358,12 @@ -- ('euclideanNorm' '<$|' v) '<.>^' w ≡ v '<.>' w -- @ (<$|) :: LSpace v => Norm v -> v -> DualVector v -Norm m <$| v = m $ v +Norm m <$| v = m-+$>v -- | The squared norm. More efficient than '|$|' because that needs to take -- the square root. normSq :: LSpace v => Seminorm v -> v -> Scalar v -normSq (Norm m) v = (m$v)<.>^v +normSq (Norm m) v = (m-+$>v)<.>^v -- | Use a 'Norm' to measure the length / norm of a vector. -- @@ -345,26 +377,30 @@ -- is similar to the dimension of the space, or even larger than it. -- Use this function to optimise the underlying operator to a dense -- matrix representation. -densifyNorm :: LSpace v => Norm v -> Norm v -densifyNorm (Norm m) = Norm . arr $ sampleLinearFunction $ m +densifyNorm :: ∀ v . LSpace v => Norm v -> Norm v +densifyNorm = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness + -> \(Norm m) -> Norm . arr $ sampleLinearFunction $ m data OrthonormalSystem v = OrthonormalSystem { orthonormalityNorm :: Norm v , orthonormalVectors :: [v] } -orthonormaliseFussily :: (LSpace v, RealFloat (Scalar v)) +orthonormaliseFussily :: ∀ v . (LSpace v, RealFloat (Scalar v)) => Scalar v -> Norm v -> [v] -> [v] -orthonormaliseFussily fuss me = go [] - where go _ [] = [] - go ws (v₀:vs) - | mvd > fuss = let μ = 1/sqrt mvd - v = vd^*μ - in v : go ((v,dvd^*μ):ws) vs - | otherwise = go ws vs - where vd = orthogonalComplementProj' ws $ v₀ - dvd = applyNorm me $ vd - mvd = dvd<.>^vd +orthonormaliseFussily = onf dualSpaceWitness + where onf :: DualSpaceWitness v -> Scalar v -> Norm v -> [v] -> [v] + onf DualSpaceWitness fuss me = go [] + where go _ [] = [] + go ws (v₀:vs) + | mvd > fuss = let μ = 1/sqrt mvd + v = vd^*μ + in v : go ((v,dvd^*μ):ws) vs + | otherwise = go ws vs + where vd = orthogonalComplementProj' ws $ v₀ + dvd = applyNorm me $ vd + mvd = dvd<.>^vd orthogonalComplementProj' :: LSpace v => [(v, DualVector v)] -> (v-+>v) orthogonalComplementProj' ws = LinearFunction $ \v₀ @@ -372,7 +408,7 @@ orthogonalComplementProj :: LSpace v => Norm v -> [v] -> (v-+>v) orthogonalComplementProj (Norm m) - = orthogonalComplementProj' . map (id &&& (m$)) + = orthogonalComplementProj' . map (id &&& (m-+$>)) @@ -380,8 +416,8 @@ ev_Eigenvalue :: Scalar v -- ^ The estimated eigenvalue @λ@. , ev_Eigenvector :: v -- ^ Normalised vector @v@ that gets mapped to a multiple, namely: , ev_FunctionApplied :: v -- ^ @f $ v ≡ λ *^ v @. - , ev_Deviation :: v -- ^ Deviation of these two supposedly equivalent expressions. - , ev_Badness :: Scalar v -- ^ Squared norm of the deviation, normalised by the eigenvalue. + , ev_Deviation :: v -- ^ Deviation of @v@ to @(f$v)^/λ@. Ideally, this would of course be equal. + , ev_Badness :: Scalar v -- ^ Squared norm of the deviation. } deriving instance (Show v, Show (Scalar v)) => Show (Eigenvector v) @@ -401,29 +437,31 @@ -> [v] -- ^ Starting vector(s) for the power method. -> [[Eigenvector v]] -- ^ Infinite sequence of ever more accurate approximations -- to the eigensystem of the operator. -constructEigenSystem me@(Norm m) ε₀ f = iterate ( +constructEigenSystem me ε₀ f = iterate ( sortBy (comparing $ negate . abs . ev_Eigenvalue) . map asEV - . orthonormaliseFussily (1/4) (Norm m) + . orthonormaliseFussily (1/4) me . newSys) . map (asEV . (^%me)) where newSys [] = [] newSys (Eigenvector λ v fv dv ε : evs) | ε>ε₀ = case newSys evs of - [] -> [fv^/λ, dv^*(sqrt $ λ^2/ε)] - vn:vns -> fv^/λ : vn : dv^*(sqrt $ λ^2/ε) : vns + [] -> [fv^/λ, dv^/sqrt(ε+ε₀)] + vn:vns -> fv^/λ : vn : dv^/sqrt(ε+ε₀) : vns | ε>=0 = v : newSys evs | otherwise = newSys evs asEV v = Eigenvector λ v fv dv ε - where λ = v'<.>^fv - ε = normSq me dv / (λ^2 + ε₀) + where λ² = fv'<.>^fv + λ = fv'<.>^v + ε = normSq me dv fv = f $ v - dv = v^*λ ^-^ fv - v' = m $ v + fv' = me<$|fv + dv | λ²>0 = v ^-^ fv^*(λ/λ²) -- for stability reasons + | otherwise = zeroV -finishEigenSystem :: (LSpace v, RealFloat (Scalar v)) +finishEigenSystem :: ∀ v . (LSpace v, RealFloat (Scalar v)) => Norm v -> [Eigenvector v] -> [Eigenvector v] finishEigenSystem me = go . sortBy (comparing $ negate . ev_Eigenvalue) where go [] = [] @@ -440,7 +478,7 @@ fShift₁v₀ = fv₀ ^-^ λ₁*^v₀ - (μ₀₀,μ₀₁) = normalized ( λ₀ - λ₁ + (μ₀₀,μ₀₁) = normalised ( λ₀ - λ₁ , (me <$| fShift₁v₀)<.>^v₁ ) (μ₁₀,μ₁₁) = (-μ₀₁, μ₀₀) @@ -457,7 +495,9 @@ where λ = (me<$|v)<.>^fv dv = v^*λ ^-^ fv ε = normSq me dv / λ^2 - + + normalised (x,y) = (x/r, y/r) + where r = sqrt $ x^2 + y^2 -- | Find a system of vectors that approximate the eigensytem, in the sense that: -- each true eigenvalue is represented by an approximate one, and that is closer @@ -476,7 +516,7 @@ | normSq me vPerp > fpε = case evss of evs':_ | length evs' > oldDim -> go (v:vs) (length evs) evss - _ -> let evss' = constructEigenSystem me fpε (arr f) + _ -> let evss' = tail . constructEigenSystem me fpε (arr f) $ map ev_Eigenvector (head $ evss++[evs]) ++ [vPerp] in go vs (length evs) evss' | otherwise = go vs oldDim (evs:evss) @@ -510,11 +550,13 @@ normSpanningSystem :: SimpleSpace v - => Norm v -> [DualVector v] -normSpanningSystem = dualBasis . normSpanningSystem' + => Seminorm v -> [DualVector v] +normSpanningSystem me@(Norm m) + = catMaybes . map snd . orthonormaliseDuals 0 + . map (id&&&(m-+$>)) $ normSpanningSystem' me normSpanningSystem' :: (FiniteDimensional v, IEEE (Scalar v)) - => Norm v -> [v] + => Seminorm v -> [v] normSpanningSystem' me = orthonormaliseFussily 0 me $ enumerateSubBasis entireBasis @@ -531,35 +573,92 @@ -- @ -- n₁ = 'spanNorm' [dv^*η | (dv,η)<-shSys] -- @ +-- +-- A rather crude approximation ('roughEigenSystem') is used in this function, so do +-- not expect the above equations to hold with great accuracy. sharedNormSpanningSystem :: SimpleSpace v - => Norm v -> Norm v -> [(DualVector v, Scalar v)] -sharedNormSpanningSystem (Norm n) (Norm m) - = sep =<< roughEigenSystem (Norm n) (pseudoInverse (arr n) . arr m) - where sep (Eigenvector λ _ λv _ _) - | λ>0 = [(n$v, sqrt λ)] + => Norm v -> Seminorm v -> [(DualVector v, Scalar v)] +sharedNormSpanningSystem nn@(Norm n) nm + = first (n-+$>) <$> sharedNormSpanningSystem' 0 (nn, dualNorm nn) nm + +sharedNormSpanningSystem' :: ∀ v . SimpleSpace v + => Int -> (Norm v, Variance v) -> Seminorm v -> [(v, Scalar v)] +sharedNormSpanningSystem' = snss dualSpaceWitness + where snss :: DualSpaceWitness v -> Int -> (Norm v, Variance v) + -> Seminorm v -> [(v, Scalar v)] + snss DualSpaceWitness nRefine (nn@(Norm n), Norm n') (Norm m) + = sep =<< iterate (finishEigenSystem nn) + (roughEigenSystem nn $ arr n' . arr m) !! nRefine + sep (Eigenvector λ v λv _ _) + | λ>=0 = [(v, sqrt λ)] | otherwise = [] - where v = λv ^/ λ +-- | Like 'sharedNormSpanningSystem n₀ n₁', but allows /either/ of the norms +-- to be singular. +-- +-- @ +-- n₀ = 'spanNorm' [dv | (dv, Just _)<-shSys] +-- @ +-- +-- and +-- +-- @ +-- n₁ = 'spanNorm' $ [dv^*η | (dv, Just η)<-shSys] +-- ++ [ dv | (dv, Nothing)<-shSys] +-- @ +-- +-- You may also interpret a @Nothing@ here as an “infinite eigenvalue”, i.e. +-- it is so small as an spanning vector of @n₀@ that you would need to scale it +-- by ∞ to use it for spanning @n₁@. +sharedSeminormSpanningSystem :: ∀ v . SimpleSpace v + => Seminorm v -> Seminorm v -> [(DualVector v, Maybe (Scalar v))] +sharedSeminormSpanningSystem nn nm + = finalise dualSpaceWitness + <$> sharedNormSpanningSystem' 1 (combined, dualNorm combined) nn + where combined = densifyNorm $ nn<>nm + finalise :: DualSpaceWitness v -> (v, Scalar v) -> (DualVector v, Maybe (Scalar v)) + finalise DualSpaceWitness (v, μn) + | μn^2 > epsilon = (v'^*μn, Just $ sqrt (1 - μn^2)/μn) + | otherwise = (v', Nothing) + where v' = combined<$|v +-- | A system of vectors which are orthogonal with respect to both of the given +-- seminorms. (In general they are not /orthonormal/ to either of them.) +sharedSeminormSpanningSystem' :: ∀ v . SimpleSpace v + => Seminorm v -> Seminorm v -> [v] +sharedSeminormSpanningSystem' nn nm + = fst <$> sharedNormSpanningSystem' 1 (combined, dualNorm combined) nn + where combined = densifyNorm $ nn<>nm + + -- | Interpret a variance as a covariance between two subspaces, and -- normalise it by the variance on @u@. The result is effectively -- the linear regression coefficient of a simple regression of the vectors -- spanning the variance. -dependence :: (SimpleSpace u, SimpleSpace v, Scalar u~Scalar v) +dependence :: ∀ u v . (SimpleSpace u, SimpleSpace v, Scalar u~Scalar v) => Variance (u,v) -> (u+>v) -dependence (Norm m) = fmap ( snd . m . (id&&&zeroV) ) - $ pseudoInverse (arr $ fst . m . (id&&&zeroV)) +dependence = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness,DualSpaceWitness) + -> \(Norm m) -> fmap ( snd . m . (id&&&zeroV) ) + $ pseudoInverse (arr $ fst . m . (id&&&zeroV)) -summandSpaceNorms :: (SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v) +summandSpaceNorms :: ∀ u v . (SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v) => Norm (u,v) -> (Norm u, Norm v) -summandSpaceNorms nuv = ( densifyNorm $ spanNorm (fst<$>spanSys) - , densifyNorm $ spanNorm (snd<$>spanSys) ) - where spanSys = normSpanningSystem nuv +summandSpaceNorms = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness,DualSpaceWitness) + -> \nuv -> let spanSys = normSpanningSystem nuv + in ( densifyNorm $ spanNorm (fst<$>spanSys) + , densifyNorm $ spanNorm (snd<$>spanSys) ) -sumSubspaceNorms :: (LSpace u, LSpace v, Scalar u~Scalar v) +sumSubspaceNorms :: ∀ u v . (LSpace u, LSpace v, Scalar u~Scalar v) => Norm u -> Norm v -> Norm (u,v) -sumSubspaceNorms (Norm nu) (Norm nv) = Norm $ nu *** nv +sumSubspaceNorms = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness,DualSpaceWitness) + -> \(Norm nu) (Norm nv) -> Norm $ nu *** nv 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 @@ -21,6 +21,7 @@ {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE TupleSections #-} {-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE GADTs #-} module Math.LinearMap.Category.Class where @@ -38,9 +39,15 @@ import Math.LinearMap.Asserted import Math.VectorSpace.ZeroDimensional -type Num' s = (Num s, VectorSpace s, Scalar s ~ s) -type Num'' s = (Num' s, LinearSpace s) -type Num''' s = (Num s, InnerSpace s, Scalar s ~ s, LSpace' s, DualVector s ~ s) +data ClosedScalarWitness s where + ClosedScalarWitness :: (Scalar s ~ s, DualVector s ~ s) => ClosedScalarWitness s + +class (Num s, LinearSpace s) => Num' s where + closedScalarWitness :: ClosedScalarWitness s + +data ScalarSpaceWitness v where + ScalarSpaceWitness :: (Num' (Scalar v), Scalar (Scalar v) ~ Scalar v) + => ScalarSpaceWitness v class (VectorSpace v) => TensorSpace v where -- | The internal representation of a 'Tensor' product. @@ -49,26 +56,27 @@ -- scalar field in the @v@ vector with an entire @w@ vector. I.e., you have -- then a “nested vector” or, if @v@ is a @DualVector@ / “row vector”, a matrix. type TensorProduct v w :: * - zeroTensor :: (LSpace w, Scalar w ~ Scalar v) + scalarSpaceWitness :: ScalarSpaceWitness v + zeroTensor :: (TensorSpace w, Scalar w ~ Scalar v) => v ⊗ w toFlatTensor :: v -+> (v ⊗ Scalar v) fromFlatTensor :: (v ⊗ Scalar v) -+> v - addTensors :: (LSpace w, Scalar w ~ Scalar v) + addTensors :: (TensorSpace w, Scalar w ~ Scalar v) => (v ⊗ w) -> (v ⊗ w) -> v ⊗ w - subtractTensors :: (LSpace v, LSpace w, Num''' (Scalar v), Scalar w ~ Scalar v) + subtractTensors :: (TensorSpace v, TensorSpace w, Scalar w ~ Scalar v) => (v ⊗ w) -> (v ⊗ w) -> v ⊗ w - subtractTensors m n = addTensors m (negateTensor $ n) - scaleTensor :: (LSpace w, Scalar w ~ Scalar v) + subtractTensors m n = addTensors m (getLinearFunction negateTensor n) + scaleTensor :: (TensorSpace w, Scalar w ~ Scalar v) => Bilinear (Scalar v) (v ⊗ w) (v ⊗ w) - negateTensor :: (LSpace w, Scalar w ~ Scalar v) + negateTensor :: (TensorSpace w, Scalar w ~ Scalar v) => (v ⊗ w) -+> (v ⊗ w) - tensorProduct :: (LSpace w, Scalar w ~ Scalar v) + tensorProduct :: (TensorSpace w, Scalar w ~ Scalar v) => Bilinear v w (v ⊗ w) - transposeTensor :: (LSpace w, Scalar w ~ Scalar v) + transposeTensor :: (TensorSpace w, Scalar w ~ Scalar v) => (v ⊗ w) -+> (w ⊗ v) - fmapTensor :: (LSpace w, LSpace x, Scalar w ~ Scalar v, Scalar x ~ Scalar v) + fmapTensor :: (TensorSpace w, TensorSpace x, Scalar w ~ Scalar v, Scalar x ~ Scalar v) => Bilinear (w -+> x) (v⊗w) (v⊗x) - fzipTensorWith :: ( LSpace u, LSpace w, LSpace x + fzipTensorWith :: ( TensorSpace u, TensorSpace w, TensorSpace x , Scalar u ~ Scalar v, Scalar w ~ Scalar v, Scalar x ~ Scalar v ) => Bilinear ((w,x) -+> u) (v⊗w, v⊗x) (v⊗u) coerceFmapTensorProduct :: Hask.Functor p @@ -77,14 +85,18 @@ infixl 7 ⊗ -- | Infix version of 'tensorProduct'. -(⊗) :: (LSpace v, LSpace w, Scalar w ~ Scalar v) +(⊗) :: ∀ v w . (TensorSpace v, TensorSpace w, Scalar w ~ Scalar v, Num' (Scalar v)) => v -> w -> v ⊗ w -v⊗w = (tensorProduct $ v) $ w +v⊗w = (tensorProduct-+$>v)-+$>w +data DualSpaceWitness v where + DualSpaceWitness :: ( LinearSpace (Scalar v), DualVector (Scalar v) ~ Scalar v + , LinearSpace (DualVector v), Scalar (DualVector v) ~ Scalar v + , DualVector (DualVector v) ~ v ) + => DualSpaceWitness v + -- | The class of vector spaces @v@ for which @'LinearMap' s v w@ is well-implemented. -class ( TensorSpace v, TensorSpace (DualVector v) - , Num' (Scalar v), Scalar (DualVector v) ~ Scalar v ) - => LinearSpace v where +class (TensorSpace v, Num (Scalar v)) => LinearSpace v where -- | Suitable representation of a linear map from the space @v@ to its field. -- -- For the usual euclidean spaces, you can just define @'DualVector' v = v@. @@ -92,66 +104,101 @@ -- @v@-vectors as “column vectors”. 'LinearMap' will then effectively have -- a matrix layout.) type DualVector v :: * + + dualSpaceWitness :: DualSpaceWitness v linearId :: v +> v - idTensor :: LSpace v => v ⊗ DualVector v - idTensor = transposeTensor $ asTensor $ linearId + idTensor :: v ⊗ DualVector v + idTensor = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> transposeTensor-+$>asTensor $ linearId - sampleLinearFunction :: (LSpace v, LSpace w, Scalar v ~ Scalar w) + sampleLinearFunction :: (TensorSpace w, Scalar v ~ Scalar w) => (v-+>w) -+> (v+>w) - sampleLinearFunction = LinearFunction $ \f -> fmap f $ id + sampleLinearFunction = case ( scalarSpaceWitness :: ScalarSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness, DualSpaceWitness) -> LinearFunction + $ \f -> getLinearFunction (fmap f) id - toLinearForm :: LSpace v => DualVector v -+> (v+>Scalar v) - toLinearForm = toFlatTensor >>> arr fromTensor + toLinearForm :: DualVector v -+> (v+>Scalar v) + toLinearForm = case ( scalarSpaceWitness :: ScalarSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness,DualSpaceWitness) -> toFlatTensor >>> arr fromTensor - fromLinearForm :: LSpace v => (v+>Scalar v) -+> DualVector v - fromLinearForm = arr asTensor >>> fromFlatTensor + fromLinearForm :: (v+>Scalar v) -+> DualVector v + fromLinearForm = case ( scalarSpaceWitness :: ScalarSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness,DualSpaceWitness) -> arr asTensor >>> fromFlatTensor coerceDoubleDual :: Coercion v (DualVector (DualVector v)) - - blockVectSpan :: (LSpace w, Scalar w ~ Scalar v) - => w -+> (v⊗(v+>w)) - blockVectSpan' :: (LSpace v, LSpace w, Num''' (Scalar v), Scalar v ~ Scalar w) - => w -+> (v+>(v⊗w)) - blockVectSpan' = LinearFunction $ \w -> fmap (flipBilin tensorProduct $ w) $ id + coerceDoubleDual = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> Coercion - trace :: LSpace v => (v+>v) -+> Scalar v - trace = flipBilin contractLinearMapAgainst $ id + trace :: (v+>v) -+> Scalar v + trace = case scalarSpaceWitness :: ScalarSpaceWitness v of + ScalarSpaceWitness -> flipBilin contractLinearMapAgainst-+$>id - contractTensorMap :: (LSpace w, Scalar w ~ Scalar v) + contractTensorMap :: (TensorSpace w, Scalar w ~ Scalar v) => (v+>(v⊗w)) -+> w - contractMapTensor :: (LSpace w, Scalar w ~ Scalar v) + contractTensorMap = case scalarSpaceWitness :: ScalarSpaceWitness v of + ScalarSpaceWitness -> arr deferLinearMap >>> transposeTensor + >>> fmap trace >>> fromFlatTensor + contractMapTensor :: (TensorSpace w, Scalar w ~ Scalar v) => (v⊗(v+>w)) -+> w - contractFnTensor :: (LSpace v, LSpace w, Scalar w ~ Scalar v) - => (v⊗(v-+>w)) -+> w - contractFnTensor = fmap sampleLinearFunction >>> contractMapTensor - contractTensorFn :: (LSpace v, LSpace w, Scalar w ~ Scalar v) + contractMapTensor = case ( scalarSpaceWitness :: ScalarSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness,DualSpaceWitness) + -> arr (coUncurryLinearMap>>>asTensor) + >>> transposeTensor >>> fmap (arr asLinearMap >>> trace) + >>> fromFlatTensor + contractTensorFn :: ∀ w . (TensorSpace w, Scalar w ~ Scalar v) => (v-+>(v⊗w)) -+> w - contractTensorFn = sampleLinearFunction >>> contractTensorMap - contractTensorWith :: (LSpace v, LSpace w, Scalar w ~ Scalar v) - => Bilinear (v⊗w) (DualVector w) v - contractTensorWith = flipBilin $ LinearFunction - (\dw -> fromFlatTensor . fmap (flipBilin applyDualVector$dw)) - contractLinearMapAgainst :: (LSpace w, Scalar w ~ Scalar v) + contractTensorFn = LinearFunction $ getLinearFunction sampleLinearFunction + >>> getLinearFunction contractTensorMap + contractLinearMapAgainst :: (LinearSpace w, Scalar w ~ Scalar v) => Bilinear (v+>w) (w-+>v) (Scalar v) + contractLinearMapAgainst = case ( scalarSpaceWitness :: ScalarSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness,DualSpaceWitness) -> arr asTensor >>> transposeTensor + >>> applyDualVector >>> LinearFunction (. sampleLinearFunction) - applyDualVector :: LSpace v + applyDualVector :: LinearSpace v => Bilinear (DualVector v) v (Scalar v) - applyLinear :: (LSpace w, Scalar w ~ Scalar v) + applyLinear :: (TensorSpace w, Scalar w ~ Scalar v) => Bilinear (v+>w) v w - composeLinear :: ( LSpace w, LSpace x + composeLinear :: ( LinearSpace w, TensorSpace x , Scalar w ~ Scalar v, Scalar x ~ Scalar v ) => Bilinear (w+>x) (v+>w) (v+>x) + composeLinear = case scalarSpaceWitness :: ScalarSpaceWitness v of + ScalarSpaceWitness -> LinearFunction $ \f -> fmap (applyLinear-+$>f) + + tensorId :: (LinearSpace w, Scalar w ~ Scalar v) + => (v⊗w)+>(v⊗w) + + applyTensorFunctional :: ( LinearSpace u, Scalar u ~ Scalar v ) + => Bilinear (DualVector (v⊗u)) (v⊗u) (Scalar v) + + applyTensorLinMap :: ( LinearSpace u, TensorSpace w + , Scalar u ~ Scalar v, Scalar w ~ Scalar v ) + => Bilinear ((v⊗u)+>w) (v⊗u) w + +fmapLinearMap :: ∀ s v w x . ( LinearSpace v, TensorSpace w, TensorSpace x + , Scalar v ~ s, Scalar w ~ s, Scalar x ~ s ) + => Bilinear (LinearFunction s w x) (v+>w) (v+>x) +fmapLinearMap = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> bilinearFunction + $ \f -> arr asTensor >>> getLinearFunction (fmapTensor-+$>f) >>> arr fromTensor -instance Num''' s => TensorSpace (ZeroDim s) where +instance Num' s => TensorSpace (ZeroDim s) where type TensorProduct (ZeroDim s) v = ZeroDim s + scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of + ClosedScalarWitness -> ScalarSpaceWitness zeroTensor = Tensor Origin toFlatTensor = LinearFunction $ \Origin -> Tensor Origin fromFlatTensor = LinearFunction $ \(Tensor Origin) -> Origin - negateTensor = const0 + negateTensor = LinearFunction id scaleTensor = biConst0 addTensors (Tensor Origin) (Tensor Origin) = Tensor Origin subtractTensors (Tensor Origin) (Tensor Origin) = Tensor Origin @@ -160,19 +207,23 @@ fmapTensor = biConst0 fzipTensorWith = biConst0 coerceFmapTensorProduct _ Coercion = Coercion -instance Num''' s => LinearSpace (ZeroDim s) where +instance Num' s => LinearSpace (ZeroDim s) where type DualVector (ZeroDim s) = ZeroDim s + dualSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of + ClosedScalarWitness -> DualSpaceWitness linearId = LinearMap Origin idTensor = Tensor Origin + tensorId = LinearMap Origin + toLinearForm = LinearFunction . const $ LinearMap Origin fromLinearForm = const0 coerceDoubleDual = Coercion contractTensorMap = const0 contractMapTensor = const0 - contractTensorWith = biConst0 contractLinearMapAgainst = biConst0 - blockVectSpan = const0 applyDualVector = biConst0 applyLinear = biConst0 + applyTensorFunctional = biConst0 + applyTensorLinMap = biConst0 composeLinear = biConst0 @@ -206,12 +257,14 @@ fromTensor :: Coercion (Tensor s (DualVector v) w) (LinearMap s v w) fromTensor = Coercion -asLinearMap :: ∀ s v w . (LSpace v, Scalar v ~ s) +asLinearMap :: ∀ s v w . (LinearSpace v, Scalar v ~ s) => Coercion (Tensor s v w) (LinearMap s (DualVector v) w) -asLinearMap = Coercion -fromLinearMap :: ∀ s v w . (LSpace v, Scalar v ~ s) +asLinearMap = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> Coercion +fromLinearMap :: ∀ s v w . (LinearSpace v, Scalar v ~ s) => Coercion (LinearMap s (DualVector v) w) (Tensor s v w) -fromLinearMap = Coercion +fromLinearMap = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> Coercion -- | Infix synonym for 'LinearMap', without explicit mention of the scalar type. type v +> w = LinearMap (Scalar v) v w @@ -219,42 +272,50 @@ -- | Infix synonym for 'Tensor', without explicit mention of the scalar type. type v ⊗ w = Tensor (Scalar v) v w -type LSpace' v = ( LinearSpace v, LinearSpace (Scalar v) - , LinearSpace (DualVector v), DualVector (DualVector v) ~ v ) - -- | The workhorse of this package: most functions here work on vector --- spaces that fulfill the @'LSpace' v@ constraint. In summary, this is: +-- spaces that fulfill the @'LSpace' v@ constraint. -- --- * A 'VectorSpace' whose 'Scalar' is a 'Num'''' (i.e. a number type that --- has itself all the vector-space instances). --- * You have an implementation for @'TensorProduct' v w@, for any other space @w@. --- * You have a 'DualVector' space that fulfills @'DualVector' ('DualVector' v) ~ v@. +-- In summary, this is a 'VectorSpace' with an implementation for @'TensorProduct' v w@, +-- for any other space @w@, and with a 'DualVector' space. This fulfills +-- @'DualVector' ('DualVector' v) ~ v@ (this constraint is encapsulated in +-- 'DualSpaceWitness'). -- --- To make a new space of yours an 'LSpace', you must define instances of --- 'TensorSpace' and 'LinearSpace'. -type LSpace v = (LSpace' v, Num''' (Scalar v)) +-- To make a new space of yours an 'LSpace', you must define instances of +-- 'TensorSpace' and 'LinearSpace'. In fact, 'LSpace' is equivalent to +-- 'LinearSpace', but makes the condition explicit that the scalar and dual vectors +-- also form a linear space. 'LinearSpace' only stores that constraint in +-- 'dualSpaceWitness' (to avoid UndecidableSuperclasses). +type LSpace v = ( LinearSpace v, LinearSpace (Scalar v), LinearSpace (DualVector v) + , Num' (Scalar v) ) -instance (LSpace v, LSpace w, Scalar v~s, Scalar w~s) +instance (LinearSpace v, TensorSpace w, Scalar v~s, Scalar w~s) => AdditiveGroup (LinearMap s v w) where - zeroV = fromTensor $ zeroTensor - m^+^n = fromTensor $ (asTensor$m) ^+^ (asTensor$n) - m^-^n = fromTensor $ (asTensor$m) ^-^ (asTensor$n) - negateV = (fromTensor$) . negateV . (asTensor$) -instance (LSpace v, LSpace w, Scalar v~s, Scalar w~s) + zeroV = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> fromTensor $ zeroTensor + m^+^n = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> fromTensor $ (asTensor$m) ^+^ (asTensor$n) + m^-^n = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> fromTensor $ (asTensor$m) ^-^ (asTensor$n) + negateV = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> (fromTensor$) . negateV . (asTensor$) +instance ∀ v w s . (LinearSpace v, TensorSpace w, Scalar v~s, Scalar w~s) => VectorSpace (LinearMap s v w) where type Scalar (LinearMap s v w) = s - μ*^v = arr fromTensor . (scaleTensor$μ) . arr asTensor $ v + μ*^v = case ( dualSpaceWitness :: DualSpaceWitness v + , scalarSpaceWitness :: ScalarSpaceWitness w ) of + (DualSpaceWitness, ScalarSpaceWitness) + -> fromTensor $ (scaleTensor-+$>μ) -+$> asTensor $ v -instance (LSpace v, LSpace w, Scalar v~s, Scalar w~s) +instance (TensorSpace v, TensorSpace w, Scalar v~s, Scalar w~s) => AdditiveGroup (Tensor s v w) where zeroV = zeroTensor (^+^) = addTensors (^-^) = subtractTensors - negateV = arr negateTensor -instance (LSpace v, LSpace w, Scalar v~s, Scalar w~s) + negateV = getLinearFunction negateTensor +instance (TensorSpace v, TensorSpace w, Scalar v~s, Scalar w~s) => VectorSpace (Tensor s v w) where type Scalar (Tensor s v w) = s - μ*^t = (scaleTensor $ μ) $ t + μ*^t = (scaleTensor-+$>μ)-+$>t infixr 6 ⊕, >+<, <⊕ @@ -274,28 +335,42 @@ instance Category (LinearMap s) where - type Object (LinearMap s) v = (LSpace v, Scalar v ~ s) + type Object (LinearMap s) v = (LinearSpace v, Scalar v ~ s) id = linearId - (.) = arr . arr composeLinear -instance Num''' s => Cartesian (LinearMap s) where + (.) = lmc dualSpaceWitness + where lmc :: ∀ v w x . ( LinearSpace v, Scalar v ~ s + , LinearSpace w, Scalar w ~ s + , TensorSpace x, Scalar x ~ s ) + => DualSpaceWitness v + -> LinearMap s w x -> LinearMap s v w -> LinearMap s v x + lmc DualSpaceWitness = getLinearFunction . getLinearFunction composeLinear +instance Num' s => Cartesian (LinearMap s) where type UnitObject (LinearMap s) = ZeroDim s swap = (fmap (const0&&&id) $ id) ⊕ (fmap (id&&&const0) $ id) attachUnit = fmap (id&&&const0) $ id detachUnit = fst regroup = sampleLinearFunction $ LinearFunction regroup regroup' = sampleLinearFunction $ LinearFunction regroup' -instance Num''' s => Morphism (LinearMap s) where +instance Num' s => Morphism (LinearMap s) where f *** g = (fmap (id&&&const0) $ f) ⊕ (fmap (const0&&&id) $ g) -instance Num''' s => PreArrow (LinearMap s) where - f &&& g = fromTensor $ (fzipTensorWith$id) $ (asTensor $ f, asTensor $ g) +instance ∀ s . Num' s => PreArrow (LinearMap s) where + (&&&) = lmFanout + where lmFanout :: ∀ u v w . ( LinearSpace u, LinearSpace v, LinearSpace w + , Scalar u~s, Scalar v~s, Scalar w~s ) + => LinearMap s u v -> LinearMap s u w -> LinearMap s u (v,w) + lmFanout f g = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness w ) of + (DualSpaceWitness, DualSpaceWitness, DualSpaceWitness) + -> fromTensor $ (fzipTensorWith$id) $ (asTensor $ f, asTensor $ g) terminal = zeroV fst = sampleLinearFunction $ fst snd = sampleLinearFunction $ snd -instance Num''' s => EnhancedCat (->) (LinearMap s) where +instance Num' s => EnhancedCat (->) (LinearMap s) where arr m = arr $ applyLinear $ m -instance Num''' s => EnhancedCat (LinearFunction s) (LinearMap s) where +instance Num' s => EnhancedCat (LinearFunction s) (LinearMap s) where arr m = applyLinear $ m -instance Num''' s => EnhancedCat (LinearMap s) (LinearFunction s) where +instance Num' s => EnhancedCat (LinearMap s) (LinearFunction s) where arr m = sampleLinearFunction $ m @@ -303,63 +378,116 @@ -instance ∀ u v . ( Num''' (Scalar v), LSpace u, LSpace v, Scalar u ~ Scalar v ) +instance ∀ u v . ( TensorSpace u, TensorSpace v, Scalar u ~ Scalar v ) => TensorSpace (u,v) where type TensorProduct (u,v) w = (u⊗w, v⊗w) + scalarSpaceWitness = case ( scalarSpaceWitness :: ScalarSpaceWitness u + , scalarSpaceWitness :: ScalarSpaceWitness v ) of + (ScalarSpaceWitness, ScalarSpaceWitness) -> ScalarSpaceWitness zeroTensor = zeroTensor <⊕ zeroTensor - scaleTensor = scaleTensor&&&scaleTensor >>> LinearFunction ( - uncurry (***) >>> pretendLike Tensor ) - negateTensor = pretendLike Tensor $ negateTensor *** negateTensor + scaleTensor = bilinearFunction $ \μ (Tensor (v,w)) -> + Tensor ( (scaleTensor-+$>μ)-+$>v, (scaleTensor-+$>μ)-+$>w ) + negateTensor = LinearFunction $ \(Tensor (v,w)) + -> Tensor (negateTensor-+$>v, negateTensor-+$>w) addTensors (Tensor (fu, fv)) (Tensor (fu', fv')) = (fu ^+^ fu') <⊕ (fv ^+^ fv') subtractTensors (Tensor (fu, fv)) (Tensor (fu', fv')) = (fu ^-^ fu') <⊕ (fv ^-^ fv') - toFlatTensor = follow Tensor <<< toFlatTensor *** toFlatTensor - fromFlatTensor = flout Tensor >>> fromFlatTensor *** fromFlatTensor - tensorProduct = LinearFunction $ \(u,v) -> - (tensorProduct$u) &&& (tensorProduct$v) >>> follow Tensor - transposeTensor = flout Tensor >>> transposeTensor *** transposeTensor >>> fzip - fmapTensor = LinearFunction $ - \f -> pretendLike Tensor $ (fmapTensor$f) *** (fmapTensor$f) + toFlatTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> follow Tensor <<< toFlatTensor *** toFlatTensor + fromFlatTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> flout Tensor >>> fromFlatTensor *** fromFlatTensor + tensorProduct = bilinearFunction $ \(u,v) w -> + Tensor ((tensorProduct-+$>u)-+$>w, (tensorProduct-+$>v)-+$>w) + transposeTensor = LinearFunction $ \(Tensor (uw,vw)) + -> (fzipTensorWith-+$>id)-+$>(transposeTensor-+$>uw,transposeTensor-+$>vw) + fmapTensor = bilinearFunction $ + \f (Tensor (uw,vw)) -> Tensor ((fmapTensor-+$>f)-+$>uw, (fmapTensor-+$>f)-+$>vw) fzipTensorWith = bilinearFunction $ \f (Tensor (uw, vw), Tensor (ux, vx)) - -> Tensor ( (fzipTensorWith$f)$(uw,ux) - , (fzipTensorWith$f)$(vw,vx) ) + -> Tensor ( (fzipTensorWith-+$>f)-+$>(uw,ux) + , (fzipTensorWith-+$>f)-+$>(vw,vx) ) coerceFmapTensorProduct p cab = case ( coerceFmapTensorProduct (fst<$>p) cab , coerceFmapTensorProduct (snd<$>p) cab ) of (Coercion, Coercion) -> Coercion -instance ∀ u v . ( LinearSpace u, LinearSpace (DualVector u), DualVector (DualVector u) ~ u - , LinearSpace v, LinearSpace (DualVector v), DualVector (DualVector v) ~ v - , Scalar u ~ Scalar v, Num''' (Scalar u) ) +instance ∀ u v . ( LinearSpace u, LinearSpace v, Scalar u ~ Scalar v ) => LinearSpace (u,v) where type DualVector (u,v) = (DualVector u, DualVector v) - linearId = (fmap (id&&&const0) $ id) ⊕ (fmap (const0&&&id) $ id) - -- idTensor = fmapTensor linearCoFst idTensor <⊕ fmapTensor linearCoSnd idTensor - sampleLinearFunction = LinearFunction $ \f -> (sampleLinearFunction $ f . lCoFst) - ⊕ (sampleLinearFunction $ f . lCoSnd) - coerceDoubleDual = Coercion - blockVectSpan = (blockVectSpan >>> fmap lfstBlock) &&& (blockVectSpan >>> fmap lsndBlock) - >>> follow Tensor - contractTensorMap = flout LinearMap - >>> contractTensorMap . fmap (fst . flout Tensor) . arr fromTensor - ***contractTensorMap . fmap (snd . flout Tensor) . arr fromTensor - >>> addV - contractMapTensor = flout Tensor - >>> contractMapTensor . fmap (arr fromTensor . fst . flout LinearMap) - ***contractMapTensor . fmap (arr fromTensor . snd . flout LinearMap) - >>> addV - contractTensorWith = LinearFunction $ \(Tensor (fu, fv)) - -> (contractTensorWith$fu) &&& (contractTensorWith$fv) - contractLinearMapAgainst = flout LinearMap >>> bilinearFunction - (\(mu,mv) f -> ((contractLinearMapAgainst$fromTensor$mu)$(fst.f)) - + ((contractLinearMapAgainst$fromTensor$mv)$(snd.f)) ) - applyDualVector = LinearFunction $ \(du,dv) + + dualSpaceWitness = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> DualSpaceWitness + linearId = case ( scalarSpaceWitness :: ScalarSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness, DualSpaceWitness, DualSpaceWitness) + -> (fmap (id&&&const0)-+$>id) ⊕ (fmap (const0&&&id)-+$>id) + tensorId = tI scalarSpaceWitness dualSpaceWitness dualSpaceWitness dualSpaceWitness + where tI :: ∀ w . (LinearSpace w, Scalar w ~ Scalar v) + => ScalarSpaceWitness u -> DualSpaceWitness u + -> DualSpaceWitness v -> DualSpaceWitness w + -> ((u,v)⊗w)+>((u,v)⊗w) + tI ScalarSpaceWitness DualSpaceWitness DualSpaceWitness DualSpaceWitness + = LinearMap + ( rassocTensor . fromLinearMap . argFromTensor + $ fmap (LinearFunction $ \t -> Tensor (t,zeroV)) -+$> tensorId + , rassocTensor . fromLinearMap . argFromTensor + $ fmap (LinearFunction $ \t -> Tensor (zeroV,t)) -+$> tensorId ) + sampleLinearFunction = case ( scalarSpaceWitness :: ScalarSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness, DualSpaceWitness, DualSpaceWitness) + -> LinearFunction $ \f -> (sampleLinearFunction -+$> f . lCoFst) + ⊕ (sampleLinearFunction -+$> f . lCoSnd) +--blockVectSpan = case ( dualSpaceWitness :: DualSpaceWitness u +-- , dualSpaceWitness :: DualSpaceWitness v ) of +-- (DualSpaceWitness, DualSpaceWitness) +-- -> (blockVectSpan >>> fmap lfstBlock) &&& (blockVectSpan >>> fmap lsndBlock) +-- >>> follow Tensor +--contractTensorMap = flout LinearMap +-- >>> contractTensorMap . fmap (fst . flout Tensor) . arr fromTensor +-- ***contractTensorMap . fmap (snd . flout Tensor) . arr fromTensor +-- >>> addV +--contractMapTensor = flout Tensor +-- >>> contractMapTensor . fmap (arr fromTensor . fst . flout LinearMap) +-- ***contractMapTensor . fmap (arr fromTensor . snd . flout LinearMap) +-- >>> addV +--contractTensorWith = LinearFunction $ \(Tensor (fu, fv)) +-- -> (contractTensorWith$fu) &&& (contractTensorWith$fv) +--contractLinearMapAgainst = flout LinearMap >>> bilinearFunction +-- (\(mu,mv) f -> ((contractLinearMapAgainst$fromTensor$mu)$(fst.f)) +-- + ((contractLinearMapAgainst$fromTensor$mv)$(snd.f)) ) + applyDualVector = case ( scalarSpaceWitness :: ScalarSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness, DualSpaceWitness, DualSpaceWitness) + -> LinearFunction $ \(du,dv) -> (applyDualVector$du) *** (applyDualVector$dv) >>> addV - applyLinear = LinearFunction $ \(LinearMap (fu, fv)) -> - (applyLinear $ (asLinearMap $ fu)) *** (applyLinear $ (asLinearMap $ fv)) + applyLinear = case ( scalarSpaceWitness :: ScalarSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness, DualSpaceWitness, DualSpaceWitness) + -> LinearFunction $ \(LinearMap (fu, fv)) -> + (applyLinear -+$> (asLinearMap $ fu)) *** (applyLinear -+$> (asLinearMap $ fv)) >>> addV - composeLinear = bilinearFunction $ \f (LinearMap (fu, fv)) - -> f . (asLinearMap $ fu) ⊕ f . (asLinearMap $ fv) + composeLinear = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) + -> bilinearFunction $ \f (LinearMap (fu, fv)) + -> ((composeLinear-+$>f)-+$>asLinearMap $ fu) + ⊕ ((composeLinear-+$>f)-+$>asLinearMap $ fv) + applyTensorFunctional = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> bilinearFunction $ + \(LinearMap (fu,fv)) (Tensor (tu,tv)) + -> ((applyTensorFunctional-+$>asLinearMap$fu)-+$>tu) + + ((applyTensorFunctional-+$>asLinearMap$fv)-+$>tv) + applyTensorLinMap = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> bilinearFunction`id` + \f (Tensor (tu,tv)) -> let LinearMap (fu,fv) = curryLinearMap $ f + in ( (applyTensorLinMap-+$>uncurryLinearMap.asLinearMap $ fu)-+$>tu ) + ^+^ ( (applyTensorLinMap-+$>uncurryLinearMap.asLinearMap $ fv)-+$>tv ) lfstBlock :: ( LSpace u, LSpace v, LSpace w , Scalar u ~ Scalar v, Scalar v ~ Scalar w ) @@ -371,6 +499,20 @@ lsndBlock = LinearFunction (zeroV⊕) +-- | @((v'⊗w)+>x) -> ((v+>w)+>x) +argFromTensor :: ∀ s v w x . (LinearSpace v, LinearSpace w, Scalar v ~ s, Scalar w ~ s) + => Coercion (LinearMap s (Tensor s (DualVector v) w) x) + (LinearMap s (LinearMap s v w) x) +argFromTensor = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> curryLinearMap >>> fromLinearMap >>> coUncurryLinearMap + +-- | @((v+>w)+>x) -> ((v'⊗w)+>x)@ +argAsTensor :: ∀ s v w x . (LinearSpace v, LinearSpace w, Scalar v ~ s, Scalar w ~ s) + => Coercion (LinearMap s (LinearMap s v w) x) + (LinearMap s (Tensor s (DualVector v) w) x) +argAsTensor = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> uncurryLinearMap <<< asLinearMap <<< coCurryLinearMap + -- | @(u+>(v⊗w)) -> (u+>v)⊗w@ deferLinearMap :: Coercion (LinearMap s u (Tensor s v w)) (Tensor s (LinearMap s u v) w) deferLinearMap = Coercion @@ -385,99 +527,163 @@ rassocTensor :: Coercion (Tensor s (Tensor s u v) w) (Tensor s u (Tensor s v w)) rassocTensor = Coercion -instance ∀ s u v . ( Num''' s, LSpace u, LSpace v, Scalar u ~ s, Scalar v ~ s ) +instance ∀ s u v . ( LinearSpace u, TensorSpace v, Scalar u ~ s, Scalar v ~ s ) => TensorSpace (LinearMap s u v) where type TensorProduct (LinearMap s u v) w = TensorProduct (DualVector u) (Tensor s v w) + scalarSpaceWitness = case ( scalarSpaceWitness :: ScalarSpaceWitness u + , scalarSpaceWitness :: ScalarSpaceWitness v ) of + (ScalarSpaceWitness, ScalarSpaceWitness) -> ScalarSpaceWitness zeroTensor = deferLinearMap $ zeroV - toFlatTensor = arr deferLinearMap . fmap toFlatTensor - fromFlatTensor = fmap fromFlatTensor . arr hasteLinearMap + toFlatTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> arr deferLinearMap . fmap toFlatTensor + fromFlatTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> fmap fromFlatTensor . arr hasteLinearMap addTensors t₁ t₂ = deferLinearMap $ (hasteLinearMap$t₁) ^+^ (hasteLinearMap$t₂) subtractTensors t₁ t₂ = deferLinearMap $ (hasteLinearMap$t₁) ^-^ (hasteLinearMap$t₂) - scaleTensor = LinearFunction $ \μ -> arr deferLinearMap . scaleWith μ . arr hasteLinearMap + scaleTensor = bilinearFunction $ \μ t + -> deferLinearMap $ scaleWith μ -+$> hasteLinearMap $ t negateTensor = arr deferLinearMap . lNegateV . arr hasteLinearMap - transposeTensor -- t :: (u +> v) ⊗ w - = arr hasteLinearMap -- u +> (v ⊗ w) + transposeTensor = case ( scalarSpaceWitness :: ScalarSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness u ) of + (ScalarSpaceWitness,DualSpaceWitness)-> --(u +> v) ⊗ w + arr hasteLinearMap -- u +> (v ⊗ w) >>> fmap transposeTensor -- u +> (w ⊗ v) >>> arr asTensor -- u' ⊗ (w ⊗ v) >>> transposeTensor -- (w ⊗ v) ⊗ u' >>> arr rassocTensor -- w ⊗ (v ⊗ u') >>> fmap transposeTensor -- w ⊗ (u' ⊗ v) >>> arr (fmap fromTensor) -- w ⊗ (u +> v) - tensorProduct = LinearFunction $ \t -> arr deferLinearMap - . (flipBilin composeLinear $ t) . blockVectSpan' - fmapTensor = LinearFunction $ \f + tensorProduct = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> bilinearFunction $ \f s + -> deferLinearMap $ fmap (flipBilin tensorProduct-+$>s)-+$>f + fmapTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> LinearFunction $ \f -> arr deferLinearMap <<< fmap (fmap f) <<< arr hasteLinearMap - fzipTensorWith = LinearFunction $ \f + fzipTensorWith = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> LinearFunction $ \f -> arr deferLinearMap <<< fzipWith (fzipWith f) <<< arr hasteLinearMap *** arr hasteLinearMap - coerceFmapTensorProduct = cftlp - where cftlp :: ∀ a b p . p (LinearMap s u v) -> Coercion a b + coerceFmapTensorProduct = cftlp dualSpaceWitness + where cftlp :: ∀ a b p . DualSpaceWitness u -> p (LinearMap s u v) -> Coercion a b -> Coercion (TensorProduct (DualVector u) (Tensor s v a)) (TensorProduct (DualVector u) (Tensor s v b)) - cftlp _ c = coerceFmapTensorProduct ([]::[DualVector u]) + cftlp DualSpaceWitness _ c + = coerceFmapTensorProduct ([]::[DualVector u]) (fmap c :: Coercion (v⊗a) (v⊗b)) --- | @((u+>v)+>w) -> v+>(u⊗w)@ -coCurryLinearMap :: Coercion (LinearMap s (LinearMap s u v) w) (LinearMap s v (Tensor s u w)) -coCurryLinearMap = Coercion +-- | @((u+>v)+>w) -> u⊗(v+>w)@ +coCurryLinearMap :: ∀ s u v w . ( LinearSpace u, Scalar u ~ s + , LinearSpace v, Scalar v ~ s ) => + Coercion (LinearMap s (LinearMap s u v) w) (Tensor s u (LinearMap s v w)) +coCurryLinearMap = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) + -> asTensor >>> rassocTensor >>> fmap asLinearMap --- | @(u+>(v⊗w)) -> (v+>u)+>w@ -coUncurryLinearMap :: Coercion (LinearMap s u (Tensor s v w)) (LinearMap s (LinearMap s v u) w) -coUncurryLinearMap = Coercion +-- | @(u⊗(v+>w)) -> (u+>v)+>w@ +coUncurryLinearMap :: ∀ s u v w . ( LinearSpace u, Scalar u ~ s + , LinearSpace v, Scalar v ~ s ) => + Coercion (Tensor s u (LinearMap s v w)) (LinearMap s (LinearMap s u v) w) +coUncurryLinearMap = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) + -> fromTensor <<< lassocTensor <<< fmap fromLinearMap -curryLinearMap :: (Num''' s, LSpace u, Scalar u ~ s) +-- | @((u⊗v)+>w) -> (u+>(v+>w))@ +curryLinearMap :: ∀ u v w s . ( LinearSpace u, Scalar u ~ s ) => Coercion (LinearMap s (Tensor s u v) w) (LinearMap s u (LinearMap s v w)) -curryLinearMap = fmap fromTensor . fromTensor . rassocTensor . asTensor +curryLinearMap = case dualSpaceWitness :: DualSpaceWitness u of + DualSpaceWitness -> (Coercion :: Coercion ((u⊗v)+>w) + ((DualVector u)⊗(Tensor s (DualVector v) w)) ) + >>> fmap fromTensor >>> fromTensor -uncurryLinearMap :: (Num''' s, LSpace u, Scalar u ~ s) +-- | @(u+>(v+>w)) -> ((u⊗v)+>w)@ +uncurryLinearMap :: ∀ u v w s . ( LinearSpace u, Scalar u ~ s ) => Coercion (LinearMap s u (LinearMap s v w)) (LinearMap s (Tensor s u v) w) -uncurryLinearMap = fromTensor . lassocTensor . asTensor . fmap asTensor +uncurryLinearMap = case dualSpaceWitness :: DualSpaceWitness u of + DualSpaceWitness -> (Coercion :: Coercion + ((DualVector u)⊗(Tensor s (DualVector v) w)) + ((u⊗v)+>w) ) + <<< fmap asTensor <<< asTensor -uncurryLinearFn :: ( Num''' s, LSpace u, LSpace v, LSpace w +uncurryLinearFn :: ( Num' s, LSpace u, LSpace v, LSpace w , Scalar u ~ s, Scalar v ~ s, Scalar w ~ s ) => LinearFunction s u (LinearMap s v w) -+> LinearFunction s (Tensor s u v) w uncurryLinearFn = bilinearFunction $ \f t -> contractMapTensor . fmap f . transposeTensor $ t -instance ∀ s u v . (Num''' s, LSpace u, LSpace v, Scalar u ~ s, Scalar v ~ s) +instance ∀ s u v . (LinearSpace u, LinearSpace v, Scalar u ~ s, Scalar v ~ s) => LinearSpace (LinearMap s u v) where - type DualVector (LinearMap s u v) = LinearMap s v u - linearId = coUncurryLinearMap $ fmap blockVectSpan $ id - coerceDoubleDual = Coercion - blockVectSpan = arr deferLinearMap - . fmap (arr (fmap coUncurryLinearMap) . blockVectSpan) - . blockVectSpan' - applyLinear = bilinearFunction $ \f g -> contractTensorMap $ (coCurryLinearMap$f) . g - applyDualVector = contractLinearMapAgainst >>> LinearFunction (. applyLinear) - composeLinear = bilinearFunction $ \f g - -> coUncurryLinearMap $ fmap (fmap $ applyLinear $ f) $ (coCurryLinearMap$g) - contractTensorMap = contractTensorMap . fmap (contractMapTensor . arr (fmap hasteLinearMap)) - . arr coCurryLinearMap - contractMapTensor = contractTensorMap . fmap (contractMapTensor . arr (fmap coCurryLinearMap)) - . arr hasteLinearMap - contractTensorWith = arr hasteLinearMap >>> bilinearFunction (\l dw - -> fmap (flipBilin contractTensorWith $ dw) $ l ) - contractLinearMapAgainst = arr coCurryLinearMap >>> bilinearFunction (\l f - -> (contractLinearMapAgainst . fmap transposeTensor $ l) - . uncurryLinearFn $f ) + type DualVector (LinearMap s u v) = Tensor s u (DualVector v) + dualSpaceWitness = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> DualSpaceWitness + linearId = case dualSpaceWitness :: DualSpaceWitness u of + DualSpaceWitness -> fromTensor . lassocTensor . fromLinearMap . fmap asTensor + . curryLinearMap . fmap fromTensor $ tensorId + tensorId = uncurryLinearMap . coUncurryLinearMap . fmap curryLinearMap + . coCurryLinearMap . fmap deferLinearMap $ id + coerceDoubleDual = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> Coercion +--blockVectSpan = arr deferLinearMap +-- . fmap (arr (fmap coUncurryLinearMap) . blockVectSpan) +-- . blockVectSpan' + applyLinear = case dualSpaceWitness :: DualSpaceWitness u of + DualSpaceWitness -> bilinearFunction $ \f g + -> let tf = argAsTensor $ f + in (applyTensorLinMap-+$>tf)-+$>fromLinearMap $ g + applyDualVector = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> flipBilin applyTensorFunctional + applyTensorFunctional = atf scalarSpaceWitness dualSpaceWitness dualSpaceWitness + where atf :: ∀ w . (LinearSpace w, Scalar w ~ s) + => ScalarSpaceWitness u -> DualSpaceWitness u -> DualSpaceWitness w + -> Bilinear ((u+>v)+>DualVector w) ((u+>v)⊗w) s + atf ScalarSpaceWitness DualSpaceWitness DualSpaceWitness + = arr (coCurryLinearMap >>> asLinearMap) + >>> applyTensorFunctional >>> bilinearFunction`id`\f t + -> f . arr (asTensor . hasteLinearMap) -+$> t + applyTensorLinMap = case dualSpaceWitness :: DualSpaceWitness u of + DualSpaceWitness -> LinearFunction $ + arr (curryLinearMap>>>coCurryLinearMap + >>>fmap uncurryLinearMap>>>coUncurryLinearMap>>>argAsTensor) + >>> \f -> LinearFunction $ \g + -> (applyTensorLinMap-+$>f) + . arr (asTensor . hasteLinearMap) -+$> g +-- -> coUncurryLinearMap $ fmap (fmap $ applyLinear $ f) $ (coCurryLinearMap$g) +--contractTensorWith = arr hasteLinearMap >>> bilinearFunction (\l dw +-- -> fmap (flipBilin contractTensorWith $ dw) $ l ) +--contractLinearMapAgainst = arr coCurryLinearMap >>> bilinearFunction (\l f +-- -> (contractLinearMapAgainst . fmap transposeTensor $ l) +-- . uncurryLinearFn $f ) -instance ∀ s u v . (Num''' s, LSpace u, LSpace v, Scalar u ~ s, Scalar v ~ s) +instance ∀ s u v . (TensorSpace u, TensorSpace v, Scalar u ~ s, Scalar v ~ s) => TensorSpace (Tensor s u v) where type TensorProduct (Tensor s u v) w = TensorProduct u (Tensor s v w) + scalarSpaceWitness = case ( scalarSpaceWitness :: ScalarSpaceWitness u + , scalarSpaceWitness :: ScalarSpaceWitness v ) of + (ScalarSpaceWitness, ScalarSpaceWitness) -> ScalarSpaceWitness zeroTensor = lassocTensor $ zeroTensor - toFlatTensor = arr lassocTensor . fmap toFlatTensor - fromFlatTensor = fmap fromFlatTensor . arr rassocTensor + toFlatTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> arr lassocTensor . fmap toFlatTensor + fromFlatTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> fmap fromFlatTensor . arr rassocTensor addTensors t₁ t₂ = lassocTensor $ (rassocTensor$t₁) ^+^ (rassocTensor$t₂) subtractTensors t₁ t₂ = lassocTensor $ (rassocTensor$t₁) ^-^ (rassocTensor$t₂) - scaleTensor = LinearFunction $ \μ -> arr lassocTensor . scaleWith μ . arr rassocTensor + scaleTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> + LinearFunction $ \μ -> arr lassocTensor . scaleWith μ . arr rassocTensor negateTensor = arr lassocTensor . lNegateV . arr rassocTensor - tensorProduct = flipBilin $ LinearFunction $ \w - -> arr lassocTensor . fmap (flipBilin tensorProduct $ w) - transposeTensor = fmap transposeTensor . arr rassocTensor + tensorProduct = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> flipBilin $ LinearFunction $ \w + -> arr lassocTensor . fmap (flipBilin tensorProduct-+$>w) + transposeTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> fmap transposeTensor . arr rassocTensor . transposeTensor . fmap transposeTensor . arr rassocTensor - fmapTensor = LinearFunction $ \f + fmapTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> LinearFunction $ \f -> arr lassocTensor <<< fmap (fmap f) <<< arr rassocTensor - fzipTensorWith = LinearFunction $ \f + fzipTensorWith = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> LinearFunction $ \f -> arr lassocTensor <<< fzipWith (fzipWith f) <<< arr rassocTensor *** arr rassocTensor coerceFmapTensorProduct = cftlp @@ -486,60 +692,80 @@ (TensorProduct u (Tensor s v b)) cftlp _ c = coerceFmapTensorProduct ([]::[u]) (fmap c :: Coercion (v⊗a) (v⊗b)) -instance ∀ s u v . (Num''' s, LSpace u, LSpace v, Scalar u ~ s, Scalar v ~ s) +instance ∀ s u v . (LinearSpace u, LinearSpace v, Scalar u ~ s, Scalar v ~ s) => LinearSpace (Tensor s u v) where - type DualVector (Tensor s u v) = Tensor s (DualVector u) (DualVector v) - linearId = uncurryLinearMap $ fmap (fmap transposeTensor . blockVectSpan') $ id - coerceDoubleDual = Coercion - blockVectSpan = arr lassocTensor . arr (fmap $ fmap uncurryLinearMap) - . fmap (transposeTensor . arr deferLinearMap) . blockVectSpan - . arr deferLinearMap . fmap transposeTensor . blockVectSpan' - applyLinear = LinearFunction $ \f -> contractMapTensor - . fmap (applyLinear$curryLinearMap$f) . transposeTensor - applyDualVector = bilinearFunction $ \f t - -> (contractLinearMapAgainst $ (fromTensor$f)) - . contractTensorWith $ t - composeLinear = bilinearFunction $ \f g - -> uncurryLinearMap $ fmap (fmap $ applyLinear $ f) $ (curryLinearMap$g) - contractTensorMap = contractTensorMap + type DualVector (Tensor s u v) = LinearMap s u (DualVector v) + linearId = tensorId + tensorId = fmap lassocTensor . uncurryLinearMap . uncurryLinearMap + . fmap curryLinearMap . curryLinearMap $ tensorId + coerceDoubleDual = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> Coercion + dualSpaceWitness = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> DualSpaceWitness +--blockVectSpan = arr lassocTensor . arr (fmap $ fmap uncurryLinearMap) +-- . fmap (transposeTensor . arr deferLinearMap) . blockVectSpan +-- . arr deferLinearMap . fmap transposeTensor . blockVectSpan' + applyLinear = applyTensorLinMap + applyDualVector = applyTensorFunctional + applyTensorFunctional = atf scalarSpaceWitness dualSpaceWitness + where atf :: ∀ w . (LinearSpace w, Scalar w ~ s) + => ScalarSpaceWitness u -> DualSpaceWitness w + -> Bilinear (LinearMap s (Tensor s u v) (DualVector w)) + (Tensor s (Tensor s u v) w) + s + atf ScalarSpaceWitness DualSpaceWitness + = arr curryLinearMap >>> applyTensorFunctional + >>> LinearFunction`id`\f -> f . arr rassocTensor + applyTensorLinMap = LinearFunction $ arr (curryLinearMap>>>curryLinearMap + >>>fmap uncurryLinearMap>>>uncurryLinearMap) + >>> \f -> (applyTensorLinMap-+$>f) . arr rassocTensor + composeLinear = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> bilinearFunction $ \f g + -> uncurryLinearMap $ fmap (fmap $ applyLinear-+$>f) $ (curryLinearMap$g) + contractTensorMap = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> contractTensorMap . fmap (transposeTensor . contractTensorMap . fmap (arr rassocTensor . transposeTensor . arr rassocTensor)) . arr curryLinearMap - contractMapTensor = contractTensorMap . fmap transposeTensor . contractMapTensor + contractMapTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> contractTensorMap . fmap transposeTensor . contractMapTensor . fmap (arr (curryLinearMap . hasteLinearMap) . transposeTensor) . arr rassocTensor - contractTensorWith = arr rassocTensor >>> bilinearFunction (\l dw - -> fmap (flipBilin contractTensorWith $ dw) $ l ) - contractLinearMapAgainst = arr curryLinearMap >>> bilinearFunction (\l f - -> (contractLinearMapAgainst $ l) - $ contractTensorMap . fmap (transposeTensor . f) ) +--contractTensorWith = arr rassocTensor >>> bilinearFunction (\l dw +-- -> fmap (flipBilin contractTensorWith $ dw) $ l ) +--contractLinearMapAgainst = arr curryLinearMap >>> bilinearFunction (\l f +-- -> (contractLinearMapAgainst $ l) +-- $ contractTensorMap . fmap (transposeTensor . f) ) type DualSpace v = v+>Scalar v -type Fractional' s = (Fractional s, Eq s, VectorSpace s, Scalar s ~ s) -type Fractional'' s = (Fractional' s, LSpace s) +type Fractional' s = (Num' s, Fractional s, Eq s, VectorSpace s) -instance (Num''' s, LSpace v, Scalar v ~ s) +instance (TensorSpace v, Num' s, Scalar v ~ s) => Functor (Tensor s v) (LinearFunction s) (LinearFunction s) where - fmap f = fmapTensor $ f -instance (Num''' s, LSpace v, Scalar v ~ s) + fmap f = getLinearFunction fmapTensor f +instance (Num' s, TensorSpace v, Scalar v ~ s) => Monoidal (Tensor s v) (LinearFunction s) (LinearFunction s) where pureUnit = const0 - fzipWith f = fzipTensorWith $ f + fzipWith f = getLinearFunction fzipTensorWith f -instance (Num''' s, LSpace v, Scalar v ~ s) +instance (LinearSpace v, Num' s, Scalar v ~ s) => Functor (LinearMap s v) (LinearFunction s) (LinearFunction s) where - fmap f = arr fromTensor . fmap f . arr asTensor -instance (Num''' s, LSpace v, Scalar v ~ s) + fmap = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> \f -> arr fromTensor . fmap f . arr asTensor +instance (Num' s, LinearSpace v, Scalar v ~ s) => Monoidal (LinearMap s v) (LinearFunction s) (LinearFunction s) where pureUnit = const0 - fzipWith f = arr asTensor *** arr asTensor >>> fzipWith f >>> arr fromTensor + fzipWith = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> \f -> arr asTensor *** arr asTensor >>> fzipWith f >>> arr fromTensor -instance (Num''' s, TensorSpace v, Scalar v ~ s) +instance (TensorSpace v, Scalar v ~ s) => Functor (Tensor s v) Coercion Coercion where fmap = crcFmap where crcFmap :: ∀ s v a b . (TensorSpace v, Scalar v ~ s) @@ -547,28 +773,30 @@ crcFmap f = case coerceFmapTensorProduct ([]::[v]) f of Coercion -> Coercion -instance (LSpace v, Num''' s, Scalar v ~ s) +instance (LinearSpace v, Scalar v ~ s) => Functor (LinearMap s v) Coercion Coercion where - fmap = crcFmap - where crcFmap :: ∀ s v a b . (LSpace v, Num''' s, Scalar v ~ s) - => Coercion a b -> Coercion (LinearMap s v a) (LinearMap s v b) - crcFmap f = case coerceFmapTensorProduct ([]::[DualVector v]) f of + fmap = crcFmap dualSpaceWitness + where crcFmap :: ∀ s v a b . (LinearSpace v, Scalar v ~ s) + => DualSpaceWitness v -> Coercion a b + -> Coercion (LinearMap s v a) (LinearMap s v b) + crcFmap DualSpaceWitness f + = case coerceFmapTensorProduct ([]::[DualVector v]) f of Coercion -> Coercion instance Category (LinearFunction s) where - type Object (LinearFunction s) v = (LSpace v, Scalar v ~ s) + type Object (LinearFunction s) v = (TensorSpace v, Scalar v ~ s) id = LinearFunction id LinearFunction f . LinearFunction g = LinearFunction $ f.g -instance Num''' s => Cartesian (LinearFunction s) where +instance Num' s => Cartesian (LinearFunction s) where type UnitObject (LinearFunction s) = ZeroDim s swap = LinearFunction swap attachUnit = LinearFunction (, Origin) detachUnit = LinearFunction fst regroup = LinearFunction regroup regroup' = LinearFunction regroup' -instance Num''' s => Morphism (LinearFunction s) where +instance Num' s => Morphism (LinearFunction s) where LinearFunction f***LinearFunction g = LinearFunction $ f***g -instance Num''' s => PreArrow (LinearFunction s) where +instance Num' s => PreArrow (LinearFunction s) where LinearFunction f&&&LinearFunction g = LinearFunction $ f&&&g fst = LinearFunction fst; snd = LinearFunction snd terminal = const0 @@ -577,84 +805,119 @@ instance EnhancedCat (LinearFunction s) Coercion where arr = LinearFunction . coerceWith -instance (LSpace w, Scalar w ~ s) +instance (LinearSpace w, Num' s, Scalar w ~ s) => Functor (LinearFunction s w) (LinearFunction s) (LinearFunction s) where fmap f = LinearFunction (f.) -deferLinearFn :: Coercion (LinearFunction s u (Tensor s v w)) - (Tensor s (LinearFunction s u v) w) -deferLinearFn = Coercion +sampleLinearFunctionFn :: ( LinearSpace u, LinearSpace v, TensorSpace w + , Scalar u ~ Scalar v, Scalar v ~ Scalar w) + => ((u-+>v)-+>w) -+> ((u+>v)+>w) +sampleLinearFunctionFn = LinearFunction $ + \f -> sampleLinearFunction -+$> f . applyLinear -hasteLinearFn :: Coercion (Tensor s (LinearFunction s u v) w) - (LinearFunction s u (Tensor s v w)) -hasteLinearFn = Coercion +fromLinearFn :: Coercion (LinearFunction s (LinearFunction s u v) w) + (Tensor s (LinearFunction s v u) w) +fromLinearFn = Coercion +asLinearFn :: Coercion (Tensor s (LinearFunction s u v) w) + (LinearFunction s (LinearFunction s v u) w) +asLinearFn = Coercion -instance (LSpace u, LSpace v, Scalar u ~ s, Scalar v ~ s) - => TensorSpace (LinearFunction s u v) where - type TensorProduct (LinearFunction s u v) w = LinearFunction s u (Tensor s v w) - zeroTensor = deferLinearFn $ const0 - toFlatTensor = arr deferLinearFn . fmap toFlatTensor - fromFlatTensor = fmap fromFlatTensor . arr hasteLinearFn - addTensors t s = deferLinearFn $ (hasteLinearFn$t)^+^(hasteLinearFn$s) - subtractTensors t s = deferLinearFn $ (hasteLinearFn$t)^-^(hasteLinearFn$s) - scaleTensor = LinearFunction $ \μ -> arr deferLinearFn . scaleWith μ . arr hasteLinearFn - negateTensor = arr deferLinearFn . lNegateV . arr hasteLinearFn - tensorProduct = flipBilin $ LinearFunction $ - \w -> arr deferLinearFn . fmap (flipBilin tensorProduct $ w) - transposeTensor = arr hasteLinearFn >>> LinearFunction tp - where tp f = fmap (LinearFunction $ \dw -> (flipBilin contractTensorWith$dw) . f) - $ idTensor - fmapTensor = bilinearFunction $ \f g - -> deferLinearFn $ fmap f . (hasteLinearFn$g) - fzipTensorWith = bilinearFunction $ \f (g,h) - -> deferLinearFn $ fzipWith f - <<< (hasteLinearFn$g)&&&(hasteLinearFn$h) - coerceFmapTensorProduct = cftpLf - where cftpLf :: ∀ s u v a b p . TensorSpace v - => p (LinearFunction s u v) -> Coercion a b - -> Coercion (LinearFunction s u (Tensor s v a)) - (LinearFunction s u (Tensor s v b)) - cftpLf p c = case coerceFmapTensorProduct ([]::[v]) c of - Coercion -> Coercion -coCurryLinearFn :: Coercion (LinearMap s (LinearFunction s u v) w) - (LinearFunction s v (Tensor s u w)) -coCurryLinearFn = Coercion +instance ∀ s u v . (LinearSpace u, LinearSpace v, Scalar u ~ s, Scalar v ~ s) + => TensorSpace (LinearFunction s u v) where + type TensorProduct (LinearFunction s u v) w = LinearFunction s (LinearFunction s v u) w + scalarSpaceWitness = case ( scalarSpaceWitness :: ScalarSpaceWitness u + , scalarSpaceWitness :: ScalarSpaceWitness v ) of + (ScalarSpaceWitness, ScalarSpaceWitness) -> ScalarSpaceWitness + zeroTensor = fromLinearFn $ const0 + toFlatTensor = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> fmap fromLinearFn $ applyDualVector + fromFlatTensor = case ( scalarSpaceWitness :: ScalarSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness u ) of + (ScalarSpaceWitness, DualSpaceWitness) + -> arr asLinearFn >>> LinearFunction`id` + \f -> let t = transposeTensor . (fmapTensor-+$>fromLinearForm) + -+$> coCurryLinearMap + $ sampleLinearFunction-+$> f . applyLinear + in applyLinear $ fromTensor $ t + addTensors t s = fromLinearFn $ (asLinearFn$t)^+^(asLinearFn$s) + subtractTensors t s = fromLinearFn $ (asLinearFn$t)^-^(asLinearFn$s) + scaleTensor = bilinearFunction $ \μ (Tensor f) -> Tensor $ μ *^ f + negateTensor = LinearFunction $ \(Tensor f) -> Tensor $ negateV f + tensorProduct = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> bilinearFunction $ \uv w -> Tensor $ + (applyDualVector-+$>uv) >>> scaleV w + transposeTensor = tt scalarSpaceWitness dualSpaceWitness + where tt :: ∀ w . (TensorSpace w, Scalar w ~ s) + => ScalarSpaceWitness u -> DualSpaceWitness u + -> Tensor s (LinearFunction s u v) w + -+> Tensor s w (LinearFunction s u v) + tt ScalarSpaceWitness DualSpaceWitness + = LinearFunction $ arr asLinearFn >>> \f + -> (fmapTensor-+$>applyLinear) + -+$> fmap fromTensor . rassocTensor + $ transposeTensor . fmap transposeTensor + -+$> fmap asTensor . coCurryLinearMap + $ sampleLinearFunctionFn -+$> f + fmapTensor = bilinearFunction $ \f -> arr asLinearFn + >>> \g -> fromLinearFn $ f . g + fzipTensorWith = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> bilinearFunction $ \f (g,h) + -> fromLinearFn $ f . ((asLinearFn$g)&&&(asLinearFn$h)) + coerceFmapTensorProduct _ Coercion = Coercion -coUncurryLinearFn :: Coercion (LinearFunction s u (Tensor s v w)) - (LinearMap s (LinearFunction s v u) w) -coUncurryLinearFn = Coercion +exposeLinearFn :: Coercion (LinearMap s (LinearFunction s u v) w) + (LinearFunction s (LinearFunction s u v) w) +exposeLinearFn = Coercion -instance (LSpace u, LSpace v, Scalar u ~ s, Scalar v ~ s) +instance (LinearSpace u, LinearSpace v, Scalar u ~ s, Scalar v ~ s) => LinearSpace (LinearFunction s u v) where type DualVector (LinearFunction s u v) = LinearFunction s v u - linearId = coUncurryLinearFn $ LinearFunction $ - \v -> fmap (fmap (scaleV v) . applyDualVector) $ idTensor + dualSpaceWitness = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> DualSpaceWitness + linearId = sym exposeLinearFn $ id + tensorId = uncurryLinearMap . sym exposeLinearFn + $ LinearFunction $ \f -> sampleLinearFunction-+$>tensorProduct-+$>f coerceDoubleDual = Coercion - blockVectSpan = arr deferLinearFn . bilinearFunction (\w u - -> fmap ( arr coUncurryLinearFn - . fmap (flipBilin tensorProduct$w) . applyLinear ) - $ (blockVectSpan$u) ) - contractTensorMap = arr coCurryLinearFn - >>> arr (fmap (fmap hasteLinearFn)) - >>> sampleLinearFunction - >>> fmap contractFnTensor - >>> contractTensorMap - contractMapTensor = arr hasteLinearFn - >>> arr (fmap (fmap coCurryLinearFn)) - >>> sampleLinearFunction - >>> fmap contractFnTensor - >>> contractTensorMap - contractLinearMapAgainst = arr coCurryLinearFn - >>> bilinearFunction (\v2uw w2uv - -> trace . fmap (contractTensorFn . fmap v2uw) - . sampleLinearFunction $ w2uv ) - applyDualVector = sampleLinearFunction >>> contractLinearMapAgainst - applyLinear = arr coCurryLinearFn >>> LinearFunction (\f - -> contractTensorFn . fmap f) - composeLinear = LinearFunction $ \f - -> arr coCurryLinearFn >>> fmap (fmap $ applyLinear $ f) - >>> arr coUncurryLinearFn + sampleLinearFunction = LinearFunction . arr $ sym exposeLinearFn +--contractLinearMapAgainst = arr coCurryLinearFn +-- >>> bilinearFunction (\v2uw w2uv +-- -> trace . fmap (contractTensorFn . fmap v2uw) +-- . sampleLinearFunction $ w2uv ) + applyDualVector = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> bilinearFunction $ + \f g -> trace . sampleLinearFunction -+$> f . g + applyLinear = bilinearFunction $ \f g -> (exposeLinearFn $ f) -+$> g + applyTensorFunctional = atf scalarSpaceWitness dualSpaceWitness + where atf :: ∀ w . (LinearSpace w, Scalar w ~ s) + => ScalarSpaceWitness u -> DualSpaceWitness w + -> LinearFunction s + (LinearMap s (LinearFunction s u v) (DualVector w)) + (LinearFunction s (Tensor s (LinearFunction s u v) w) s) + atf ScalarSpaceWitness DualSpaceWitness = bilinearFunction $ \f g + -> trace -+$> fromTensor $ transposeTensor + -+$> fmap ((exposeLinearFn $ f) . applyLinear) + -+$> ( transposeTensor + -+$> deferLinearMap + $ fmap transposeTensor + -+$> hasteLinearMap + $ transposeTensor + -+$> coCurryLinearMap + $ sampleLinearFunctionFn + -+$> asLinearFn $ g ) + applyTensorLinMap = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> bilinearFunction $ \f g + -> contractMapTensor . transposeTensor + -+$> fmap ((asLinearFn $ g) . applyLinear) + -+$> ( transposeTensor + -+$> deferLinearMap + $ fmap transposeTensor + -+$> hasteLinearMap + $ transposeTensor + -+$> coCurryLinearMap + $ sampleLinearFunctionFn + -+$> exposeLinearFn . curryLinearMap $ f ) 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 @@ -48,12 +48,20 @@ import Math.VectorSpace.ZeroDimensional +(<.>^) :: LinearSpace v => DualVector v -> v -> Scalar v +f<.>^v = (applyDualVector-+$>f)-+$>v + + type ℝ = Double +instance Num' ℝ where + closedScalarWitness = ClosedScalarWitness + instance TensorSpace ℝ where type TensorProduct ℝ w = w + scalarSpaceWitness = ScalarSpaceWitness zeroTensor = Tensor zeroV - scaleTensor = LinearFunction (pretendLike Tensor) . scale + scaleTensor = bilinearFunction $ \μ (Tensor t) -> Tensor $ μ*^t addTensors (Tensor v) (Tensor w) = Tensor $ v ^+^ w subtractTensors (Tensor v) (Tensor w) = Tensor $ v ^-^ w negateTensor = pretendLike Tensor lNegateV @@ -67,30 +75,38 @@ coerceFmapTensorProduct _ Coercion = Coercion instance LinearSpace ℝ where type DualVector ℝ = ℝ + dualSpaceWitness = DualSpaceWitness linearId = LinearMap 1 + tensorId = uncurryLinearMap $ LinearMap $ fmap (follow Tensor) -+$> id idTensor = Tensor 1 fromLinearForm = flout LinearMap coerceDoubleDual = Coercion contractTensorMap = flout Tensor . flout LinearMap contractMapTensor = flout LinearMap . flout Tensor - contractTensorWith = flout Tensor >>> applyDualVector - contractLinearMapAgainst = flout LinearMap >>> flipBilin lApply - blockVectSpan = follow Tensor . follow LinearMap applyDualVector = scale - applyLinear = elacs . flout LinearMap - composeLinear = LinearFunction $ \f -> follow LinearMap . arr f . flout LinearMap + applyLinear = LinearFunction $ \(LinearMap w) -> scaleV w + applyTensorFunctional = bilinearFunction $ \(LinearMap du) (Tensor u) -> du<.>^u + applyTensorLinMap = bilinearFunction $ \fℝuw (Tensor u) + -> let LinearMap fuw = curryLinearMap $ fℝuw + in (applyLinear-+$>fuw) -+$> u + composeLinear = bilinearFunction $ \f (LinearMap g) + -> LinearMap $ (applyLinear-+$>f)-+$>g -#define FreeLinearSpace(V, LV, tp, bspan, tenspl, dspan, contraction, contraaction) \ -instance Num''' s => TensorSpace (V s) where { \ +#define FreeLinearSpace(V, LV, tp, tenspl, tenid, dspan, contraction, contraaction) \ +instance ∀ s . Num' s => TensorSpace (V s) where { \ type TensorProduct (V s) w = V w; \ + scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of{ \ + ClosedScalarWitness -> ScalarSpaceWitness}; \ zeroTensor = Tensor $ pure zeroV; \ addTensors (Tensor m) (Tensor n) = Tensor $ liftA2 (^+^) m n; \ subtractTensors (Tensor m) (Tensor n) = Tensor $ liftA2 (^-^) m n; \ negateTensor = LinearFunction $ Tensor . fmap negateV . getTensorProduct; \ scaleTensor = bilinearFunction \ $ \μ -> Tensor . fmap (μ*^) . getTensorProduct; \ - toFlatTensor = follow Tensor; \ - fromFlatTensor = flout Tensor; \ + toFlatTensor = case closedScalarWitness :: ClosedScalarWitness s of{ \ + ClosedScalarWitness -> follow Tensor}; \ + fromFlatTensor = case closedScalarWitness :: ClosedScalarWitness s of{ \ + ClosedScalarWitness -> flout Tensor}; \ tensorProduct = bilinearFunction $ \w v -> Tensor $ fmap (*^v) w; \ transposeTensor = LinearFunction (tp); \ fmapTensor = bilinearFunction $ \ @@ -99,36 +115,48 @@ \(LinearFunction f) (Tensor vw, Tensor vx) \ -> Tensor $ liftA2 (curry f) vw vx; \ coerceFmapTensorProduct _ Coercion = Coercion }; \ -instance Num''' s => LinearSpace (V s) where { \ +instance ∀ s . Num' s => LinearSpace (V s) where { \ type DualVector (V s) = V s; \ + dualSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of \ + {ClosedScalarWitness -> DualSpaceWitness}; \ linearId = LV Mat.identity; \ idTensor = Tensor Mat.identity; \ + tensorId = ti dualSpaceWitness where \ + { ti :: ∀ w . (LinearSpace w, Scalar w ~ s) => DualSpaceWitness w -> (V s⊗w)+>(V s⊗w) \ + ; ti DualSpaceWitness = LinearMap $ \ + fmap (\f -> fmap (LinearFunction $ Tensor . f)-+$>asTensor $ id) \ + (tenid :: V (w -> V w)) }; \ coerceDoubleDual = Coercion; \ - fromLinearForm = flout LinearMap; \ - blockVectSpan = LinearFunction $ Tensor . (bspan); \ + fromLinearForm = case closedScalarWitness :: ClosedScalarWitness s of{ \ + ClosedScalarWitness -> flout LinearMap}; \ contractTensorMap = LinearFunction $ (contraction) . coerce . getLinearMap; \ contractMapTensor = LinearFunction $ (contraction) . coerce . getTensorProduct; \ - contractTensorWith = bilinearFunction $ \ - \(Tensor wv) dw -> fmap (arr $ applyDualVector $ dw) wv; \ +{-contractTensorWith = bilinearFunction $ \ + \(Tensor wv) dw -> fmap (arr $ applyDualVector $ dw) wv; -} \ contractLinearMapAgainst = bilinearFunction $ getLinearMap >>> (contraaction); \ applyDualVector = bilinearFunction Mat.dot; \ applyLinear = bilinearFunction $ \(LV m) \ -> foldl' (^+^) zeroV . liftA2 (^*) m; \ + applyTensorFunctional = bilinearFunction $ \(LinearMap f) (Tensor t) \ + -> sum $ liftA2 (<.>^) f t; \ + applyTensorLinMap = bilinearFunction $ \(LinearMap f) (Tensor t) \ + -> foldl' (^+^) zeroV $ liftA2 (arr fromTensor >>> \ + getLinearFunction . getLinearFunction applyLinear) f t; \ composeLinear = bilinearFunction $ \ - \f (LinearMap g) -> LinearMap $ fmap (f$) g } + \f (LinearMap g) -> LinearMap $ fmap ((applyLinear-+$>f)-+$>) g } FreeLinearSpace( V0 , LinearMap , \(Tensor V0) -> zeroV - , \_ -> V0 , \_ -> LinearMap V0 + , V0 , LinearMap V0 , \V0 -> zeroV , \V0 _ -> 0 ) FreeLinearSpace( V1 , LinearMap , \(Tensor (V1 w₀)) -> w₀⊗V1 1 - , \w -> V1 (LinearMap $ V1 w) , \w -> LinearMap $ V1 (Tensor $ V1 w) + , V1 V1 , LinearMap . V1 . blockVectSpan $ V1 1 , \(V1 (V1 w)) -> w , \(V1 x) f -> (f$x)^._x ) @@ -136,10 +164,9 @@ , LinearMap , \(Tensor (V2 w₀ w₁)) -> w₀⊗V2 1 0 ^+^ w₁⊗V2 0 1 - , \w -> V2 (LinearMap $ V2 w zeroV) - (LinearMap $ V2 zeroV w) , \w -> LinearMap $ V2 (Tensor $ V2 w zeroV) (Tensor $ V2 zeroV w) + , V2 (`V2`zeroV) (V2 zeroV) , LinearMap $ V2 (blockVectSpan $ V2 1 0) (blockVectSpan $ V2 0 1) , \(V2 (V2 w₀ _) @@ -150,12 +177,12 @@ , \(Tensor (V3 w₀ w₁ w₂)) -> w₀⊗V3 1 0 0 ^+^ w₁⊗V3 0 1 0 ^+^ w₂⊗V3 0 0 1 - , \w -> V3 (LinearMap $ V3 w zeroV zeroV) - (LinearMap $ V3 zeroV w zeroV) - (LinearMap $ V3 zeroV zeroV w) , \w -> LinearMap $ V3 (Tensor $ V3 w zeroV zeroV) (Tensor $ V3 zeroV w zeroV) (Tensor $ V3 zeroV zeroV w) + , V3 (\w -> V3 w zeroV zeroV) + (\w -> V3 zeroV w zeroV) + (\w -> V3 zeroV zeroV w) , LinearMap $ V3 (blockVectSpan $ V3 1 0 0) (blockVectSpan $ V3 0 1 0) (blockVectSpan $ V3 0 0 1) @@ -173,10 +200,10 @@ (LinearMap $ V4 zeroV w zeroV zeroV) (LinearMap $ V4 zeroV zeroV w zeroV) (LinearMap $ V4 zeroV zeroV zeroV w) - , \w -> LinearMap $ V4 (Tensor $ V4 w zeroV zeroV zeroV) - (Tensor $ V4 zeroV w zeroV zeroV) - (Tensor $ V4 zeroV zeroV w zeroV) - (Tensor $ V4 zeroV zeroV zeroV w) + , V4 (\w -> V4 w zeroV zeroV zeroV) + (\w -> V4 zeroV w zeroV zeroV) + (\w -> V4 zeroV zeroV w zeroV) + (\w -> V4 zeroV zeroV zeroV w) , LinearMap $ V4 (blockVectSpan $ V4 1 0 0 0) (blockVectSpan $ V4 0 1 0 0) (blockVectSpan $ V4 0 0 1 0) @@ -189,7 +216,7 @@ -instance (Num''' n, TensorProduct (DualVector n) n ~ n) => Num (LinearMap n n n) where +instance (Num' n, TensorProduct (DualVector n) n ~ n) => Num (LinearMap n n n) where LinearMap n + LinearMap m = LinearMap $ n + m LinearMap n - LinearMap m = LinearMap $ n - m LinearMap n * LinearMap m = LinearMap $ n * m @@ -197,7 +224,7 @@ signum (LinearMap n) = LinearMap $ signum n fromInteger = LinearMap . fromInteger -instance (Fractional'' n, TensorProduct (DualVector n) n ~ n) +instance (Fractional' n, TensorProduct (DualVector n) n ~ n) => Fractional (LinearMap n n n) where LinearMap n / LinearMap m = LinearMap $ n / m recip (LinearMap n) = LinearMap $ recip n diff --git a/Math/VectorSpace/Docile.hs b/Math/VectorSpace/Docile.hs --- a/Math/VectorSpace/Docile.hs +++ b/Math/VectorSpace/Docile.hs @@ -48,7 +48,7 @@ import Control.Arrow.Constrained import Linear ( V0(V0), V1(V1), V2(V2), V3(V3), V4(V4) - , _x, _y, _z, _w ) + , _x, _y, _z, _w, ex, ey, ez, ew ) import qualified Data.Vector.Unboxed as UArr import Data.VectorSpace.Free import Math.VectorSpace.ZeroDimensional @@ -70,7 +70,7 @@ -- infinite-dimensional space. -- -- Of course, this also works for spaces which are already finite-dimensional themselves. -class LSpace v => SemiInner v where +class LinearSpace v => SemiInner v where -- | Lazily enumerate choices of a basis of functionals that can be made dual -- to the given vectors, in order of preference (which roughly means, large in -- the normal direction.) I.e., if the vector @𝑣@ is assigned early to the @@ -84,6 +84,9 @@ -- For simple finite-dimensional array-vectors, you can easily define this -- method using 'cartesianDualBasisCandidates'. dualBasisCandidates :: [(Int,v)] -> Forest (Int, DualVector v) + + tensorDualBasisCandidates :: (SemiInner w, Scalar w ~ Scalar v) + => [(Int, v⊗w)] -> Forest (Int, DualVector (v⊗w)) cartesianDualBasisCandidates :: [DualVector v] -- ^ Set of canonical basis functionals. @@ -115,41 +118,90 @@ (_,[]) -> [] (f,s:l') -> s : f++l' -instance (Fractional'' s, SemiInner s) => SemiInner (ZeroDim s) where +instance (Fractional' s, SemiInner s) => SemiInner (ZeroDim s) where dualBasisCandidates _ = [] -instance (Fractional'' s, SemiInner s) => SemiInner (V0 s) where + tensorDualBasisCandidates _ = [] +instance (Fractional' s, SemiInner s) => SemiInner (V0 s) where dualBasisCandidates _ = [] - -(<.>^) :: LSpace v => DualVector v -> v -> Scalar v -f<.>^v = (applyDualVector$f)$v + tensorDualBasisCandidates _ = [] -orthonormaliseDuals :: (SemiInner v, LSpace v, RealFrac' (Scalar v)) - => Scalar v -> [(v, DualVector v)] -> [(v,DualVector v)] -orthonormaliseDuals _ [] = [] -orthonormaliseDuals ε ((v,v'₀):ws) - | abs ovl > ε = (v,v') : [(w, w' ^-^ (w'<.>^v)*^v') | (w,w')<-wssys] - | otherwise = (v,zeroV) : wssys - where wssys = orthonormaliseDuals ε ws - v'₁ = foldl' (\v'i (w,w') -> v'i ^-^ (v'i<.>^w)*^w') (v'₀ ^/ (v'₀<.>^v)) wssys - v' = v'₁ ^/ ovl - ovl = v'₁<.>^v +orthonormaliseDuals :: ∀ v . (SemiInner v, LSpace v, RealFrac' (Scalar v)) + => Scalar v -> [(v, DualVector v)] + -> [(v,Maybe (DualVector v))] +orthonormaliseDuals = od dualSpaceWitness + where od _ _ [] = [] + od (DualSpaceWitness :: DualSpaceWitness v) ε ((v,v'₀):ws) + | abs ovl₀ > 0, abs ovl₁ > ε + = (v,Just v') + : [ (w, fmap (\w' -> w' ^-^ (w'<.>^v)*^v') w's) + | (w,w's)<-wssys ] + | otherwise = (v,Nothing) : wssys + where wssys = orthonormaliseDuals ε ws + v'₁ = foldl' (\v'i₀ (w,w's) + -> foldl' (\v'i w' -> v'i ^-^ (v'i<.>^w)*^w') v'i₀ w's) + (v'₀ ^/ ovl₀) wssys + v' = v'₁ ^/ ovl₁ + ovl₀ = v'₀<.>^v + ovl₁ = v'₁<.>^v -dualBasis :: (SemiInner v, LSpace v, RealFrac' (Scalar v)) => [v] -> [DualVector v] -dualBasis vs = snd <$> orthonormaliseDuals epsilon (zip' vsIxed candidates) +dualBasis :: ∀ v . (SemiInner v, LSpace v, RealFrac' (Scalar v)) + => [v] -> [Maybe (DualVector v)] +dualBasis vs = snd <$> result where zip' ((i,v):vs) ((j,v'):ds) | i^(lookupArr Arr.! i) /= 0 - , Just best' <- findBest (n-1) bv' - = Just $ (i,v') : best' - | otherwise = findBest n alts + result :: [(v, Maybe (DualVector v))] + result = case findBest n n $ dualBasisCandidates vsIxed of + Right bestCandidates + -> orthonormaliseDuals epsilon + (zip' vsIxed $ sortBy (comparing fst) bestCandidates) + Left (_, bestCompromise) + -> let survivors :: [(Int, DualVector v)] + casualties :: [Int] + (casualties, survivors) + = second (sortBy $ comparing fst) + $ mapEither (\case + (i,Nothing) -> Left i + (i,Just v') -> Right (i,v') + ) bestCompromise + bestEffort = orthonormaliseDuals epsilon + [ (lookupArr Arr.! i, v') + | (i,v') <- survivors ] + in map snd . sortBy (comparing fst) + $ zipWith ((,) . fst) survivors bestEffort + ++ [ (i,(lookupArr Arr.! i, Nothing)) + | i <- casualties ] + where findBest :: Int -- ^ Dual vectors needed for complete dual basis + -> Int -- ^ Maximum numbers of alternatives to consider + -- (to prevent exponential blowup of possibilities) + -> Forest (Int, DualVector v) + -> Either (Int, [(Int, Maybe (DualVector v))]) + [(Int, DualVector v)] + findBest 0 _ _ = Right [] + findBest nMissing _ [] = Left (nMissing, []) + findBest n maxCompromises (Node (i,v') bv' : alts) + | Just _ <- guardedv' + , Right best' <- straightContinue = Right $ (i,v') : best' + | maxCompromises > 0 + , Right goodAlt <- alternative = Right goodAlt + | otherwise = case straightContinue of + Right goodOtherwise -> Left (1, second Just <$> goodOtherwise) + Left (nBad, badAnyway) + | maxCompromises > 0 + , Left (nBadAlt, badAlt) <- alternative + , nBadAlt < nBad + myBadness + -> Left (nBadAlt, badAlt) + | otherwise -> Left ( nBad + myBadness + , (i, guardedv') : badAnyway ) + where guardedv' = case v'<.>^(lookupArr Arr.! i) of + 0 -> Nothing + _ -> Just v' + myBadness = case guardedv' of + Nothing -> 1 + Just _ -> 0 + straightContinue = findBest (n-1) (maxCompromises-1) bv' + alternative = findBest n (maxCompromises-1) alts vsIxed = zip [0..] vs lookupArr = Arr.fromList vs n = Arr.length lookupArr @@ -158,57 +210,84 @@ dualBasisCandidates = fmap ((`Node`[]) . second recip) . sortBy (comparing $ negate . abs . snd) . filter ((/=0) . snd) + tensorDualBasisCandidates = map (second getTensorProduct) + >>> dualBasisCandidates + >>> fmap (fmap $ second LinearMap) -instance (Fractional'' s, Ord s, SemiInner s) => SemiInner (V1 s) where +instance (Fractional' s, Ord s, SemiInner s) => SemiInner (V1 s) where dualBasisCandidates = fmap ((`Node`[]) . second recip) . sortBy (comparing $ negate . abs . snd) . filter ((/=0) . snd) + tensorDualBasisCandidates = map (second $ \(Tensor (V1 w)) -> w) + >>> dualBasisCandidates + >>> fmap (fmap . second $ LinearMap . V1) -#define FreeSemiInner(V, sabs) \ -instance SemiInner (V) where { \ - dualBasisCandidates \ - = cartesianDualBasisCandidates Mat.basis (fmap sabs . toList) } -FreeSemiInner(V2 ℝ, abs) -FreeSemiInner(V3 ℝ, abs) -FreeSemiInner(V4 ℝ, abs) +instance SemiInner (V2 ℝ) where + dualBasisCandidates = cartesianDualBasisCandidates Mat.basis (toList . fmap abs) + tensorDualBasisCandidates = map (second $ \(Tensor (V2 x y)) -> (x,y)) + >>> dualBasisCandidates + >>> map (fmap . second $ LinearMap . \(dx,dy) -> V2 dx dy) +instance SemiInner (V3 ℝ) where + dualBasisCandidates = cartesianDualBasisCandidates Mat.basis (toList . fmap abs) + tensorDualBasisCandidates = map (second $ \(Tensor (V3 x y z)) -> (x,(y,z))) + >>> dualBasisCandidates + >>> map (fmap . second $ LinearMap . \(dx,(dy,dz)) -> V3 dx dy dz) +instance SemiInner (V4 ℝ) where + dualBasisCandidates = cartesianDualBasisCandidates Mat.basis (toList . fmap abs) + tensorDualBasisCandidates = map (second $ \(Tensor (V4 x y z w)) -> ((x,y),(z,w))) + >>> dualBasisCandidates + >>> map (fmap . second $ LinearMap . \((dx,dy),(dz,dw)) -> V4 dx dy dz dw) instance ∀ u v . ( SemiInner u, SemiInner v, Scalar u ~ Scalar v ) => SemiInner (u,v) where dualBasisCandidates = fmap (\(i,(u,v))->((i,u),(i,v))) >>> unzip >>> dualBasisCandidates *** dualBasisCandidates - >>> combineBaseis False mempty - where combineBaseis :: Bool -> Set Int + >>> combineBaseis (dualSpaceWitness,dualSpaceWitness) False mempty + where combineBaseis :: (DualSpaceWitness u, DualSpaceWitness v) -> Bool -> Set Int -> ( Forest (Int, DualVector u) , Forest (Int, DualVector v) ) -> Forest (Int, (DualVector u, DualVector v)) - combineBaseis _ _ ([], []) = [] - combineBaseis False forbidden (Node (i,du) bu' : abu, bv) - | i`Set.member`forbidden = combineBaseis False forbidden (abu, bv) + combineBaseis _ _ _ ([], []) = [] + combineBaseis wit@(DualSpaceWitness,DualSpaceWitness) + False forbidden (Node (i,du) bu' : abu, bv) + | i`Set.member`forbidden = combineBaseis wit False forbidden (abu, bv) | otherwise = Node (i, (du, zeroV)) - (combineBaseis True (Set.insert i forbidden) (bu', bv)) - : combineBaseis False forbidden (abu, bv) - combineBaseis True forbidden (bu, Node (i,dv) bv' : abv) - | i`Set.member`forbidden = combineBaseis True forbidden (bu, abv) + (combineBaseis wit True (Set.insert i forbidden) (bu', bv)) + : combineBaseis wit False forbidden (abu, bv) + combineBaseis wit@(DualSpaceWitness,DualSpaceWitness) + True forbidden (bu, Node (i,dv) bv' : abv) + | i`Set.member`forbidden = combineBaseis wit True forbidden (bu, abv) | otherwise = Node (i, (zeroV, dv)) - (combineBaseis False (Set.insert i forbidden) (bu, bv')) - : combineBaseis True forbidden (bu, abv) - combineBaseis _ forbidden (bu, []) = combineBaseis False forbidden (bu,[]) - combineBaseis _ forbidden ([], bv) = combineBaseis True forbidden ([],bv) + (combineBaseis wit False (Set.insert i forbidden) (bu, bv')) + : combineBaseis wit True forbidden (bu, abv) + combineBaseis wit _ forbidden (bu, []) = combineBaseis wit False forbidden (bu,[]) + combineBaseis wit _ forbidden ([], bv) = combineBaseis wit True forbidden ([],bv) + tensorDualBasisCandidates = case scalarSpaceWitness :: ScalarSpaceWitness u of + ScalarSpaceWitness -> map (second $ \(Tensor (tu, tv)) -> (tu, tv)) + >>> dualBasisCandidates + >>> map (fmap . second $ \(LinearMap lu, LinearMap lv) + -> LinearMap $ (Tensor lu, Tensor lv) ) -instance ∀ s u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ s, Scalar v ~ s ) +instance ∀ s u v . ( SemiInner u, SemiInner v, Scalar u ~ s, Scalar v ~ s ) => SemiInner (Tensor s u v) where - dualBasisCandidates = map (fmap (second $ arr transposeTensor . arr asTensor)) - . dualBasisCandidates - . map (second $ arr asLinearMap) + dualBasisCandidates = tensorDualBasisCandidates + tensorDualBasisCandidates = map (second $ arr rassocTensor) + >>> tensorDualBasisCandidates + >>> map (fmap . second $ arr uncurryLinearMap) -instance ∀ s u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ s, Scalar v ~ s ) +instance ∀ s u v . ( LinearSpace u, SemiInner (DualVector u), SemiInner v + , Scalar u ~ s, Scalar v ~ s ) => SemiInner (LinearMap s u v) where - dualBasisCandidates = sequenceForest - . map (second pseudoInverse) -- this is not efficient - where sequenceForest [] = [] - sequenceForest (x:xs) = [Node x $ sequenceForest xs] + dualBasisCandidates = case dualSpaceWitness :: DualSpaceWitness u of + DualSpaceWitness -> (coerce :: [(Int, LinearMap s u v)] + -> [(Int, Tensor s (DualVector u) v)]) + >>> tensorDualBasisCandidates + >>> coerce + tensorDualBasisCandidates = map (second $ arr hasteLinearMap) + >>> dualBasisCandidates + >>> map (fmap . second $ arr coUncurryLinearMap) (^/^) :: (InnerSpace v, Eq (Scalar v), Fractional (Scalar v)) => v -> v -> Scalar v v^/^w = case (v<.>w) of @@ -217,7 +296,7 @@ type DList x = [x]->[x] -class (LSpace v, LSpace (Scalar v)) => FiniteDimensional v where +class (LSpace 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. @@ -259,7 +338,7 @@ recomposeContraLinMapTensor :: ( FiniteDimensional u, LinearSpace w , Scalar u ~ Scalar v, Scalar w ~ Scalar v, Hask.Functor f ) - => (f (Scalar w) -> w) -> f (DualVector v⊗DualVector u) -> (v⊗u)+>w + => (f (Scalar w) -> w) -> f (v+>DualVector u) -> (v⊗u)+>w -- | The existance of a finite basis gives us an isomorphism between a space -- and its dual space. Note that this isomorphism is not natural (i.e. it @@ -269,7 +348,7 @@ uncanonicallyToDual :: v -+> DualVector v -instance (Num''' s) => FiniteDimensional (ZeroDim s) where +instance (Num' s) => FiniteDimensional (ZeroDim s) where data SubBasis (ZeroDim s) = ZeroBasis entireBasis = ZeroBasis enumerateSubBasis ZeroBasis = [] @@ -284,7 +363,7 @@ uncanonicallyFromDual = id uncanonicallyToDual = id -instance (Num''' s, LinearSpace s) => FiniteDimensional (V0 s) where +instance (Num' s, LinearSpace s) => FiniteDimensional (V0 s) where data SubBasis (V0 s) = V0Basis entireBasis = V0Basis enumerateSubBasis V0Basis = [] @@ -312,12 +391,12 @@ decomposeLinMapWithin RealsBasis (LinearMap v) = pure (v:) recomposeContraLinMap fw = LinearMap . fw recomposeContraLinMapTensor fw = arr uncurryLinearMap . LinearMap - . recomposeContraLinMap fw . fmap getTensorProduct + . recomposeContraLinMap fw . fmap getLinearMap uncanonicallyFromDual = id uncanonicallyToDual = id #define FreeFiniteDimensional(V, VB, dimens, take, give) \ -instance (Num''' s, LSpace s) \ +instance (Num' s, LSpace s) \ => FiniteDimensional (V s) where { \ data SubBasis (V s) = VB deriving (Show); \ entireBasis = VB; \ @@ -334,10 +413,14 @@ decomposeLinMapWithin VB (LinearMap m) = pure (toList m ++); \ recomposeContraLinMap fw mv \ = LinearMap $ (\v -> fw $ fmap (<.>^v) mv) <$> Mat.identity; \ - recomposeContraLinMapTensor fw mv = LinearMap $ \ + recomposeContraLinMapTensor = rclmt dualSpaceWitness \ + where {rclmt :: ∀ u w f . ( FiniteDimensional u, LinearSpace w \ + , Scalar u ~ s, Scalar w ~ s, Hask.Functor f ) => DualSpaceWitness u \ + -> (f (Scalar w) -> w) -> f (V s+>DualVector u) -> (V s⊗u)+>w \ + ; rclmt DualSpaceWitness fw mv = LinearMap $ \ (\v -> fromLinearMap $ recomposeContraLinMap fw \ - $ fmap (\(Tensor q) -> foldl' (^+^) zeroV $ liftA2 (*^) v q) mv) \ - <$> Mat.identity } + $ fmap (\(LinearMap q) -> foldl' (^+^) zeroV $ liftA2 (*^) v q) mv) \ + <$> Mat.identity } } 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₂ ) @@ -352,24 +435,34 @@ deriving instance Show (SubBasis ℝ) -instance ( FiniteDimensional u, FiniteDimensional v - , Scalar u ~ Scalar v ) +instance ∀ u v . ( FiniteDimensional u, FiniteDimensional v + , Scalar u ~ Scalar v, Scalar (DualVector u) ~ Scalar (DualVector v) ) => FiniteDimensional (u,v) where data SubBasis (u,v) = TupleBasis !(SubBasis u) !(SubBasis v) entireBasis = TupleBasis entireBasis entireBasis enumerateSubBasis (TupleBasis bu bv) = ((,zeroV)<$>enumerateSubBasis bu) ++ ((zeroV,)<$>enumerateSubBasis bv) subbasisDimension (TupleBasis bu bv) = subbasisDimension bu + subbasisDimension bv - decomposeLinMap (LinearMap (fu, fv)) - = case (decomposeLinMap (asLinearMap$fu), decomposeLinMap (asLinearMap$fv)) of - ((bu, du), (bv, dv)) -> (TupleBasis bu bv, du . dv) - decomposeLinMapWithin (TupleBasis bu bv) (LinearMap (fu, fv)) - = case ( decomposeLinMapWithin bu (asLinearMap$fu) - , decomposeLinMapWithin bv (asLinearMap$fv) ) of - (Left (bu', du), Left (bv', dv)) -> Left (TupleBasis bu' bv', du . dv) - (Left (bu', du), Right dv) -> Left (TupleBasis bu' bv, du . dv) - (Right du, Left (bv', dv)) -> Left (TupleBasis bu bv', du . dv) - (Right du, Right dv) -> Right $ du . dv + decomposeLinMap = dclm dualSpaceWitness dualSpaceWitness dualSpaceWitness + where dclm :: ∀ w . (LinearSpace w, Scalar w ~ Scalar u) + => DualSpaceWitness u -> DualSpaceWitness v -> DualSpaceWitness w + -> ((u,v)+>w) -> (SubBasis (u,v), DList w) + dclm DualSpaceWitness DualSpaceWitness DualSpaceWitness (LinearMap (fu, fv)) + = case (decomposeLinMap (asLinearMap$fu), decomposeLinMap (asLinearMap$fv)) of + ((bu, du), (bv, dv)) -> (TupleBasis bu bv, du . dv) + decomposeLinMapWithin = dclm dualSpaceWitness dualSpaceWitness dualSpaceWitness + where dclm :: ∀ w . (LinearSpace w, Scalar w ~ Scalar u) + => DualSpaceWitness u -> DualSpaceWitness v -> DualSpaceWitness w + -> SubBasis (u,v) -> ((u,v)+>w) + -> Either (SubBasis (u,v), DList w) (DList w) + dclm DualSpaceWitness DualSpaceWitness DualSpaceWitness + (TupleBasis bu bv) (LinearMap (fu, fv)) + = case ( decomposeLinMapWithin bu (asLinearMap$fu) + , decomposeLinMapWithin bv (asLinearMap$fv) ) of + (Left (bu', du), Left (bv', dv)) -> Left (TupleBasis bu' bv', du . dv) + (Left (bu', du), Right dv) -> Left (TupleBasis bu' bv, du . dv) + (Right du, Left (bv', dv)) -> Left (TupleBasis bu bv', du . dv) + (Right du, Right dv) -> Right $ du . dv recomposeSB (TupleBasis bu bv) coefs = case recomposeSB bu coefs of (u, coefs') -> case recomposeSB bv coefs' of (v, coefs'') -> ((u,v), coefs'') @@ -381,14 +474,24 @@ recomposeContraLinMap fw dds = recomposeContraLinMap fw (fst<$>dds) ⊕ recomposeContraLinMap fw (snd<$>dds) - recomposeContraLinMapTensor fw dds - = uncurryLinearMap + recomposeContraLinMapTensor fw dds = case ( scalarSpaceWitness :: ScalarSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness,DualSpaceWitness,DualSpaceWitness) -> uncurryLinearMap $ LinearMap ( fromLinearMap . curryLinearMap - $ recomposeContraLinMapTensor fw (fmap (\(Tensor(tu,_))->tu) dds) + $ recomposeContraLinMapTensor fw + (fmap (\(LinearMap(Tensor tu,_))->LinearMap tu) dds) , fromLinearMap . curryLinearMap - $ recomposeContraLinMapTensor fw (fmap (\(Tensor(_,tv))->tv) dds) ) - uncanonicallyFromDual = uncanonicallyFromDual *** uncanonicallyFromDual - uncanonicallyToDual = uncanonicallyToDual *** uncanonicallyToDual + $ recomposeContraLinMapTensor fw + (fmap (\(LinearMap(_,Tensor tv))->LinearMap tv) dds) ) + uncanonicallyFromDual = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness,DualSpaceWitness) + -> uncanonicallyFromDual *** uncanonicallyFromDual + uncanonicallyToDual = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness,DualSpaceWitness) + -> uncanonicallyToDual *** uncanonicallyToDual deriving instance (Show (SubBasis u), Show (SubBasis v)) => Show (SubBasis (u,v)) @@ -396,35 +499,67 @@ instance ∀ s u v . ( FiniteDimensional u, FiniteDimensional v - , Scalar u~s, Scalar v~s, Fractional' (Scalar v) ) + , Scalar u~s, Scalar v~s, Scalar (DualVector u)~s, Scalar (DualVector v)~s + , Fractional' (Scalar v) ) => FiniteDimensional (Tensor s u v) where data SubBasis (Tensor s u v) = TensorBasis !(SubBasis u) !(SubBasis v) entireBasis = TensorBasis entireBasis entireBasis enumerateSubBasis (TensorBasis bu bv) = [ u⊗v | u <- enumerateSubBasis bu, v <- enumerateSubBasis bv ] subbasisDimension (TensorBasis bu bv) = subbasisDimension bu * subbasisDimension bv - decomposeLinMap muvw = case decomposeLinMap $ curryLinearMap $ muvw of - (bu, mvwsg) -> first (TensorBasis bu) . go $ mvwsg [] - where (go, _) = tensorLinmapDecompositionhelpers - decomposeLinMapWithin (TensorBasis bu bv) muvw + decomposeLinMap = dlm dualSpaceWitness + where dlm :: ∀ w . (LSpace w, Scalar w ~ Scalar v) + => DualSpaceWitness w -> ((u⊗v)+>w) -> (SubBasis (u⊗v), DList w) + dlm DualSpaceWitness muvw = case decomposeLinMap $ curryLinearMap $ muvw of + (bu, mvwsg) -> first (TensorBasis bu) . go $ mvwsg [] + where (go, _) = tensorLinmapDecompositionhelpers + decomposeLinMapWithin = dlm dualSpaceWitness + where dlm :: ∀ w . (LSpace w, Scalar w ~ Scalar v) + => DualSpaceWitness w -> SubBasis (u⊗v) + -> ((u⊗v)+>w) -> Either (SubBasis (u⊗v), DList w) (DList w) + dlm DualSpaceWitness (TensorBasis bu bv) muvw = case decomposeLinMapWithin bu $ curryLinearMap $ muvw of - Left (bu', mvwsg) -> let (_, (bv', ws)) = goWith bv id (mvwsg []) id - in Left (TensorBasis bu' bv', ws) - where (_, goWith) = tensorLinmapDecompositionhelpers + Left (bu', mvwsg) -> let (_, (bv', ws)) = goWith bv id (mvwsg []) id + in Left (TensorBasis bu' bv', ws) + where (_, goWith) = tensorLinmapDecompositionhelpers recomposeSB (TensorBasis bu bv) = recomposeSBTensor bu bv - recomposeSBTensor (TensorBasis bu bv) bw + recomposeSBTensor = rst dualSpaceWitness + where rst :: ∀ w . (FiniteDimensional w, Scalar w ~ s) + => DualSpaceWitness w -> SubBasis (u⊗v) + -> SubBasis w -> [s] -> ((u⊗v)⊗w, [s]) + rst DualSpaceWitness (TensorBasis bu bv) bw = first (arr lassocTensor) . recomposeSBTensor bu (TensorBasis bv bw) - recomposeLinMap (TensorBasis bu bv) ws = - ( uncurryLinearMap $ fst . recomposeLinMap bu $ unfoldr (pure . recomposeLinMap bv) ws - , drop (subbasisDimension bu * subbasisDimension bv) ws ) - recomposeContraLinMap = recomposeContraLinMapTensor - recomposeContraLinMapTensor fw dds - = uncurryLinearMap . uncurryLinearMap . fmap (curryLinearMap) . curryLinearMap - $ recomposeContraLinMapTensor fw $ fmap (arr rassocTensor) dds - uncanonicallyToDual = fmap uncanonicallyToDual + recomposeLinMap = rlm dualSpaceWitness + where rlm :: ∀ w . (LSpace w, Scalar w ~ Scalar v) + => DualSpaceWitness w -> SubBasis (u⊗v) -> [w] + -> ((u⊗v)+>w, [w]) + rlm DualSpaceWitness (TensorBasis bu bv) ws + = ( uncurryLinearMap $ fst . recomposeLinMap bu + $ unfoldr (pure . recomposeLinMap bv) ws + , drop (subbasisDimension bu * subbasisDimension bv) ws ) + recomposeContraLinMap = case dualSpaceWitness :: DualSpaceWitness u of + DualSpaceWitness -> recomposeContraLinMapTensor + recomposeContraLinMapTensor = rclt dualSpaceWitness dualSpaceWitness + where rclt :: ∀ u' w f . ( FiniteDimensional u', Scalar u' ~ s + , LinearSpace w, Scalar w ~ s + , Hask.Functor f ) + => DualSpaceWitness u -> DualSpaceWitness u' + -> (f (Scalar w) -> w) + -> f (Tensor s u v +> DualVector u') + -> (Tensor s u v ⊗ u') +> w + rclt DualSpaceWitness DualSpaceWitness fw dds + = uncurryLinearMap . uncurryLinearMap + . fmap (curryLinearMap) . curryLinearMap + $ recomposeContraLinMapTensor fw $ fmap (arr curryLinearMap) dds + uncanonicallyToDual = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> fmap uncanonicallyToDual >>> transposeTensor >>> fmap uncanonicallyToDual - >>> transposeTensor - uncanonicallyFromDual = fmap uncanonicallyFromDual + >>> transposeTensor >>> arr fromTensor + uncanonicallyFromDual = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> arr asTensor + >>> fmap uncanonicallyFromDual >>> transposeTensor >>> fmap uncanonicallyFromDual >>> transposeTensor @@ -458,37 +593,72 @@ instance ∀ s u v . ( LSpace u, FiniteDimensional (DualVector u), FiniteDimensional v - , Scalar u~s, Scalar v~s, Fractional' (Scalar 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 entireBasis of TensorBasis bu bv -> LinMapBasis bu bv - enumerateSubBasis (LinMapBasis bu bv) - = arr (fmap asLinearMap) . enumerateSubBasis $ TensorBasis bu bv + entireBasis = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) + -> case entireBasis of TensorBasis bu bv -> LinMapBasis bu bv + enumerateSubBasis + = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> \(LinMapBasis bu bv) + -> arr (fmap asLinearMap) . enumerateSubBasis $ TensorBasis bu bv subbasisDimension (LinMapBasis bu bv) = subbasisDimension bu * subbasisDimension bv - decomposeLinMap = first (\(TensorBasis bv bu)->LinMapBasis bu bv) + decomposeLinMap = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) + -> first (\(TensorBasis bu bv)->LinMapBasis bu bv) . decomposeLinMap . coerce - decomposeLinMapWithin (LinMapBasis bu bv) m - = case decomposeLinMapWithin (TensorBasis bv bu) (coerce m) of + decomposeLinMapWithin = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) + -> \(LinMapBasis bu bv) m + -> case decomposeLinMapWithin (TensorBasis bu bv) (coerce m) of Right ws -> Right ws - Left (TensorBasis bv' bu', ws) -> Left (LinMapBasis bu' bv', ws) - recomposeSB (LinMapBasis bu bv) - = recomposeSB (TensorBasis bu bv) >>> first (arr fromTensor) - recomposeSBTensor (LinMapBasis bu bv) bw - = recomposeSBTensor (TensorBasis bu bv) bw >>> first coerce - recomposeLinMap (LinMapBasis bu bv) ws = - ( coUncurryLinearMap . fmap asTensor $ fst . recomposeLinMap bv - $ unfoldr (pure . recomposeLinMap bu) ws - , drop (subbasisDimension bu * subbasisDimension bv) ws ) - recomposeContraLinMap fw dds = coUncurryLinearMap . fmap fromLinearMap . curryLinearMap - $ recomposeContraLinMapTensor fw $ fmap (arr asTensor) dds - recomposeContraLinMapTensor fw dds - = uncurryLinearMap . coUncurryLinearMap - . fmap (fromLinearMap . curryLinearMap) . curryLinearMap - $ recomposeContraLinMapTensor fw $ fmap (arr $ asTensor . hasteLinearMap) dds - uncanonicallyToDual = fmap uncanonicallyToDual >>> arr asTensor - >>> transposeTensor >>> arr fromTensor >>> fmap uncanonicallyToDual - uncanonicallyFromDual = fmap uncanonicallyFromDual >>> arr asTensor - >>> transposeTensor >>> arr fromTensor >>> fmap uncanonicallyFromDual + Left (TensorBasis bu' bv', ws) -> Left (LinMapBasis bu' bv', ws) + recomposeSB = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> \(LinMapBasis bu bv) + -> recomposeSB (TensorBasis bu bv) >>> first (arr fromTensor) + recomposeSBTensor = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> \(LinMapBasis bu bv) bw + -> recomposeSBTensor (TensorBasis bu bv) bw >>> first coerce + recomposeLinMap = rlm dualSpaceWitness dualSpaceWitness + where rlm :: ∀ w . (LSpace w, Scalar w ~ Scalar v) + => DualSpaceWitness u -> DualSpaceWitness w -> SubBasis (u+>v) -> [w] + -> ((u+>v)+>w, [w]) + rlm 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 + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) -> \fw dds + -> argFromTensor $ recomposeContraLinMapTensor fw $ fmap (arr asLinearMap) dds + recomposeContraLinMapTensor = rclmt dualSpaceWitness 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' + -> (f (Scalar w) -> w) -> f ((u+>v)+>DualVector u') -> ((u+>v)⊗u')+>w + rclmt DualSpaceWitness DualSpaceWitness DualSpaceWitness fw dds + = uncurryLinearMap . coUncurryLinearMap + . fmap curryLinearMap . coCurryLinearMap . argFromTensor + $ recomposeContraLinMapTensor fw + $ fmap (arr $ asLinearMap . coCurryLinearMap) dds + uncanonicallyToDual = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) + -> arr asTensor >>> fmap uncanonicallyToDual >>> transposeTensor + >>> fmap uncanonicallyToDual >>> transposeTensor + uncanonicallyFromDual = case ( dualSpaceWitness :: DualSpaceWitness u + , dualSpaceWitness :: DualSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness) + -> arr fromTensor <<< fmap uncanonicallyFromDual <<< transposeTensor + <<< fmap uncanonicallyFromDual <<< transposeTensor deriving instance (Show (SubBasis (DualVector u)), Show (SubBasis v)) => Show (SubBasis (LinearMap s u v)) @@ -525,11 +695,11 @@ (\$) :: ∀ u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v ) => (u+>v) -> v -> u (\$) m - | du > dv = (unsafeRightInverse m $) - | du < dv = (unsafeLeftInverse m $) + | du > dv = ((applyLinear-+$>unsafeRightInverse m)-+$>) + | du < dv = ((applyLinear-+$>unsafeLeftInverse m)-+$>) | otherwise = let v's = dualBasis $ mdecomp [] (mbas, mdecomp) = decomposeLinMap m - in fst . \v -> recomposeSB mbas [v'<.>^v | v' <- v's] + in fst . \v -> recomposeSB mbas [ maybe 0 (<.>^v) v' | v' <- v's ] where du = subbasisDimension (entireBasis :: SubBasis u) dv = subbasisDimension (entireBasis :: SubBasis v) @@ -550,9 +720,12 @@ -- @ unsafeLeftInverse :: ∀ u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v ) => (u+>v) -> v+>u -unsafeLeftInverse m = unsafeInverse (m' . (fmap uncanonicallyToDual $ m)) +unsafeLeftInverse = uli dualSpaceWitness dualSpaceWitness + where uli :: DualSpaceWitness u -> DualSpaceWitness v -> (u+>v) -> v+>u + uli DualSpaceWitness DualSpaceWitness m + = unsafeInverse (m' . (fmap uncanonicallyToDual $ m)) . m' . arr uncanonicallyToDual - where m' = adjoint $ m :: DualVector v +> DualVector u + where m' = adjoint $ m :: DualVector v +> DualVector u -- | If @f@ is surjective, then -- @@ -561,41 +734,54 @@ -- @ unsafeRightInverse :: ∀ u v . ( SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v ) => (u+>v) -> v+>u -unsafeRightInverse m = (fmap uncanonicallyToDual $ m') +unsafeRightInverse = uri dualSpaceWitness dualSpaceWitness + where uri :: DualSpaceWitness u -> DualSpaceWitness v -> (u+>v) -> v+>u + uri DualSpaceWitness DualSpaceWitness m + = (fmap uncanonicallyToDual $ m') . unsafeInverse (m . (fmap uncanonicallyToDual $ m')) - where m' = adjoint $ m :: DualVector v +> DualVector u + 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 ) => (u+>v) -> v+>u -unsafeInverse m = recomposeContraLinMap (fst . recomposeSB mbas) v's +unsafeInverse m = recomposeContraLinMap (fst . recomposeSB mbas) + $ [maybe zeroV id v' | v'<-v's] where v's = dualBasis $ mdecomp [] (mbas, mdecomp) = decomposeLinMap m -- | The -- provides an isomorphism between a Hilbert space and its (continuous) dual space. -riesz :: (FiniteDimensional v, InnerSpace v) => DualVector v -+> v -riesz = LinearFunction $ \dv -> +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 [] -sRiesz :: FiniteDimensional v => DualSpace v -+> v -sRiesz = LinearFunction $ \dv -> +sRiesz :: ∀ v . FiniteDimensional v => DualSpace v -+> v +sRiesz = case ( scalarSpaceWitness :: ScalarSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness,DualSpaceWitness) -> LinearFunction $ \dv -> let (bas, compos) = decomposeLinMap $ dv in fst . recomposeSB bas $ compos [] -coRiesz :: (LSpace v, Num''' (Scalar v), InnerSpace v) => v -+> DualVector v -coRiesz = fromFlatTensor . arr asTensor . sampleLinearFunction . inner +coRiesz :: ∀ v . (LSpace v, InnerSpace v) => v -+> DualVector v +coRiesz = case ( scalarSpaceWitness :: ScalarSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness,DualSpaceWitness) + -> fromFlatTensor . arr asTensor . sampleLinearFunction . inner -- | Functions are generally a pain to display, but since linear functionals -- in a Hilbert space can be represented by /vectors/ in that space, -- this can be used for implementing a 'Show' instance. -showsPrecAsRiesz :: ( FiniteDimensional v, InnerSpace v, Show v - , HasBasis (Scalar v), Basis (Scalar v) ~ () ) +showsPrecAsRiesz :: ∀ v . ( FiniteDimensional v, InnerSpace v, Show v + , HasBasis (Scalar v), Basis (Scalar v) ~ () ) => Int -> DualSpace v -> ShowS -showsPrecAsRiesz p dv = showParen (p>0) $ ("().<"++) - . showsPrec 7 (sRiesz$dv) +showsPrecAsRiesz = case ( scalarSpaceWitness :: ScalarSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness v ) of + (ScalarSpaceWitness,DualSpaceWitness) + -> \p dv -> showParen (p>0) $ ("().<"++) . showsPrec 7 (sRiesz$dv) instance Show (LinearMap ℝ (V0 ℝ) ℝ) where showsPrec = showsPrecAsRiesz instance Show (LinearMap ℝ ℝ ℝ) where showsPrec = showsPrecAsRiesz @@ -604,53 +790,168 @@ instance Show (LinearMap ℝ (V3 ℝ) ℝ) where showsPrec = showsPrecAsRiesz instance Show (LinearMap ℝ (V4 ℝ) ℝ) where showsPrec = showsPrecAsRiesz +class TensorDecomposable u => RieszDecomposable u where + rieszDecomposition :: (FiniteDimensional v, v ~ DualVector v, Scalar v ~ Scalar u) + => (v +> u) -> [(Basis u, v)] +instance RieszDecomposable ℝ where + rieszDecomposition (LinearMap r) = [((), fromFlatTensor $ Tensor r)] +instance ( RieszDecomposable x, RieszDecomposable y + , Scalar x ~ Scalar y, Scalar (DualVector x) ~ Scalar (DualVector y) ) + => RieszDecomposable (x,y) where + rieszDecomposition m = map (first Left) (rieszDecomposition $ fst . m) + ++ map (first Right) (rieszDecomposition $ snd . m) + +instance RieszDecomposable (V0 ℝ) where + rieszDecomposition _ = [] +instance RieszDecomposable (V1 ℝ) where + rieszDecomposition m = [(ex, sRiesz $ fmap (LinearFunction (^._x)) $ m)] +instance RieszDecomposable (V2 ℝ) where + rieszDecomposition m = [ (ex, sRiesz $ fmap (LinearFunction (^._x)) $ m) + , (ey, sRiesz $ fmap (LinearFunction (^._y)) $ m) ] +instance RieszDecomposable (V3 ℝ) where + rieszDecomposition m = [ (ex, sRiesz $ fmap (LinearFunction (^._x)) $ m) + , (ey, sRiesz $ fmap (LinearFunction (^._y)) $ m) + , (ez, sRiesz $ fmap (LinearFunction (^._z)) $ m) ] +instance RieszDecomposable (V4 ℝ) where + rieszDecomposition m = [ (ex, sRiesz $ fmap (LinearFunction (^._x)) $ m) + , (ey, sRiesz $ fmap (LinearFunction (^._y)) $ m) + , (ez, sRiesz $ fmap (LinearFunction (^._z)) $ m) + , (ew, sRiesz $ fmap (LinearFunction (^._w)) $ m) ] + infixl 7 .< -- | Outer product of a general @v@-vector and a basis element from @w@. -- Note that this operation is in general pretty inefficient; it is -- provided mostly to lay out matrix definitions neatly. -(.<) :: ( FiniteDimensional v, Num''' (Scalar v) +(.<) :: ( FiniteDimensional v, Num' (Scalar v) , InnerSpace v, LSpace w, HasBasis w, Scalar v ~ Scalar w ) => Basis w -> v -> v+>w bw .< v = sampleLinearFunction $ LinearFunction $ \v' -> recompose [(bw, v<.>v')] + +rieszDecomposeShowsPrec :: ∀ u v s . ( RieszDecomposable u + , FiniteDimensional v, v ~ DualVector v, Show v + , Scalar u ~ s, Scalar v ~ s ) + => Int -> LinearMap s v u -> ShowS +rieszDecomposeShowsPrec p m = case rieszDecomposition m of + [] -> ("zeroV"++) + ((b₀,dv₀):dvs) -> showParen (p>6) + $ \s -> showsPrecBasis ([]::[u]) 7 b₀ + . (".<"++) . showsPrec 7 dv₀ + $ foldr (\(b,dv) + -> (" ^+^ "++) . showsPrecBasis ([]::[u]) 7 b + . (".<"++) . showsPrec 7 dv) s dvs + instance Show (LinearMap s v (V0 s)) where show _ = "zeroV" instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) => Show (LinearMap ℝ v (V1 ℝ)) where - showsPrec p m = showParen (p>6) $ ("ex .< "++) - . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m) + showsPrec = rieszDecomposeShowsPrec instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) => Show (LinearMap ℝ v (V2 ℝ)) where - showsPrec p m = showParen (p>6) - $ ("ex.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m) - . (" ^+^ ey.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._y)) $ m) + showsPrec = rieszDecomposeShowsPrec instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) => Show (LinearMap ℝ v (V3 ℝ)) where - showsPrec p m = showParen (p>6) - $ ("ex.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m) - . (" ^+^ ey.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._y)) $ m) - . (" ^+^ ez.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._z)) $ m) + showsPrec = rieszDecomposeShowsPrec instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) => Show (LinearMap ℝ v (V4 ℝ)) where - showsPrec p m = showParen (p>6) - $ ("ex.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._x)) $ m) - . (" ^+^ ey.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._y)) $ m) - . (" ^+^ ez.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._z)) $ m) - . (" ^+^ ew.<"++) . showsPrec 7 (sRiesz $ fmap (LinearFunction (^._w)) $ m) + showsPrec = rieszDecomposeShowsPrec +instance ( FiniteDimensional v, v ~ DualVector v, Show v + , RieszDecomposable x, RieszDecomposable y + , Scalar x ~ s, Scalar y ~ s, Scalar v ~ s + , Scalar (DualVector x) ~ s, Scalar (DualVector y) ~ s ) + => Show (LinearMap s v (x,y)) where + showsPrec = case + (dualSpaceWitness::DualSpaceWitness x, dualSpaceWitness::DualSpaceWitness y) of + (DualSpaceWitness, DualSpaceWitness) -> rieszDecomposeShowsPrec +infixr 7 .⊗ +(.⊗) :: ( TensorSpace v, HasBasis v, TensorSpace w + , Num' (Scalar v), Scalar v ~ Scalar w ) + => Basis v -> w -> v⊗w +b .⊗ w = basisValue b ⊗ w +class (FiniteDimensional v, HasBasis v) => TensorDecomposable v where + tensorDecomposition :: v⊗w -> [(Basis v, w)] + showsPrecBasis :: Hask.Functor p => p v -> 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) ) + => 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 + +instance TensorDecomposable (V0 ℝ) where + tensorDecomposition _ = [] + showsPrecBasis _ _ (Mat.E q) = (V0^.q ++) +instance TensorDecomposable (V1 ℝ) where + tensorDecomposition (Tensor (V1 w)) = [(ex, w)] + 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 ++) +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 ++) +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 ++) + +tensorDecomposeShowsPrec :: ∀ u v s + . ( TensorDecomposable u, FiniteDimensional v, Show v, Scalar u ~ s, Scalar v ~ s ) + => Int -> Tensor s u v -> ShowS +tensorDecomposeShowsPrec p t = case tensorDecomposition t of + [] -> ("zeroV"++) + ((b₀,dv₀):dvs) -> showParen (p>6) + $ \s -> showsPrecBasis ([]::[u]) 7 b₀ + . (".⊗"++) . showsPrec 7 dv₀ + $ foldr (\(b,dv) + -> (" ^+^ "++) . showsPrecBasis ([]::[u]) 7 b + . (".⊗"++) . showsPrec 7 dv) s dvs + +instance Show (Tensor s (V0 s) v) where + show _ = "zeroV" +instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) + => Show (Tensor ℝ (V1 ℝ) v) where + showsPrec = tensorDecomposeShowsPrec +instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) + => Show (Tensor ℝ (V2 ℝ) v) where + showsPrec = tensorDecomposeShowsPrec +instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) + => Show (Tensor ℝ (V3 ℝ) v) where + showsPrec = tensorDecomposeShowsPrec +instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) + => Show (Tensor ℝ (V4 ℝ) v) where + showsPrec = tensorDecomposeShowsPrec + +instance ( FiniteDimensional v, v ~ DualVector v, Show v + , TensorDecomposable x, TensorDecomposable y + , Scalar x ~ s, Scalar y ~ s, Scalar v ~ s ) + => Show (Tensor s (x,y) v) where + showsPrec = case + (dualSpaceWitness::DualSpaceWitness x, dualSpaceWitness::DualSpaceWitness y) of + (DualSpaceWitness, DualSpaceWitness) -> tensorDecomposeShowsPrec + + (^) :: Num a => a -> Int -> a (^) = (Hask.^) type HilbertSpace v = (LSpace v, InnerSpace v, DualVector v ~ v) -type RealFrac' s = (IEEE s, HilbertSpace s, Scalar s ~ s) +type RealFrac' s = (Fractional' s, IEEE s, InnerSpace s) type RealFloat' s = (RealFrac' s, Floating s) type SimpleSpace v = ( FiniteDimensional v, FiniteDimensional (DualVector v) @@ -672,7 +973,8 @@ instance ∀ s u v . ( LSpace u, FiniteDimensional (DualVector u), LSpace v, FiniteFreeSpace v - , Scalar u~s, Scalar v~s ) => FiniteFreeSpace (Tensor s u v) where + , Scalar u~s, Scalar v~s, Scalar (DualVector u)~s, Scalar (DualVector v)~s ) + => FiniteFreeSpace (Tensor s u v) where freeDimension _ = subbasisDimension (entireBasis :: SubBasis (DualVector u)) * freeDimension ([]::[v]) toFullUnboxVect = arr asLinearMap >>> decomposeLinMapWithin entireBasis >>> \case @@ -709,6 +1011,9 @@ -- -- But /not/ @(v+>w) -> (w+>v)@, in general (though in a Hilbert space, this too is -- equivalent, via 'riesz' isomorphism). -adjoint :: (LSpace v, LSpace w, Scalar v ~ Scalar w) +adjoint :: ∀ v w . (LSpace v, LSpace w, Scalar v ~ Scalar w) => (v +> DualVector w) -+> (w +> DualVector v) -adjoint = arr fromTensor . transposeTensor . arr asTensor +adjoint = case ( dualSpaceWitness :: DualSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness w ) of + (DualSpaceWitness, DualSpaceWitness) + -> arr fromTensor . transposeTensor . arr asTensor 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.1.0.1 +version: 0.2.0.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 @@ -45,7 +45,7 @@ Math.LinearMap.Category.Instances Math.VectorSpace.Docile other-extensions: FlexibleInstances, UndecidableInstances, FunctionalDependencies, TypeOperators, TypeFamilies - build-depends: base >=4.8 && <4.9, + build-depends: base >=4.8 && <5, vector-space >=0.10 && <0.11, constrained-categories >=0.3 && <0.4, containers, vector,