diff --git a/Math/LinearMap/Category.hs b/Math/LinearMap/Category.hs --- a/Math/LinearMap/Category.hs +++ b/Math/LinearMap/Category.hs @@ -28,6 +28,7 @@ -- ** Function implementation LinearFunction (..), type (-+>)(), Bilinear + , lfun -- ** Tensor implementation , LinearMap (..), type (+>)() , (⊕), (>+<) @@ -37,6 +38,9 @@ , (<.>^), (-+|>) -- * Tensor spaces , Tensor (..), type (⊗)(), (⊗) + -- ** Symmetric + , SymmetricTensor(..), squareV, squareVs + , type (⊗〃+>)(), currySymBilin -- * Norms -- $metricIntro , Norm(..), Seminorm @@ -49,12 +53,13 @@ , normSpanningSystem , normSpanningSystem' -- ** Variances - , Variance, spanVariance, varianceSpanningSystem + , Variance, spanVariance, (|&>), varianceSpanningSystem , dualNorm, dualNorm', dependence -- ** Utility - , densifyNorm + , densifyNorm, wellDefinedNorm -- * Solving linear equations , (\$), pseudoInverse, roughDet + , linearRegressionW, linearRegressionWVar -- * Eigenvalue problems , eigen , constructEigenSystem @@ -66,7 +71,7 @@ , TensorSpace (..) , LinearSpace (..) -- ** Orthonormal systems - , SemiInner (..), cartesianDualBasisCandidates + , SemiInner (..), cartesianDualBasisCandidates, embedFreeSubspace -- ** Finite baseis , FiniteDimensional (..) -- * Utility @@ -80,7 +85,7 @@ , HilbertSpace, SimpleSpace , Num'(..) , Fractional' - , RealFrac', RealFloat' + , RealFrac', RealFloat', LinearShowable -- ** Double-dual, scalar-scalar etc. identity , ClosedScalarWitness(..), ScalarSpaceWitness(..), DualSpaceWitness(..) , LinearManifoldWitness(..) @@ -90,6 +95,8 @@ , summandSpaceNorms, sumSubspaceNorms , sharedNormSpanningSystem, sharedSeminormSpanningSystem , sharedSeminormSpanningSystem' + , convexPolytopeHull + , convexPolytopeRepresentatives ) where import Math.LinearMap.Category.Class @@ -124,8 +131,13 @@ import qualified Linear.Vector as Mat import Control.Lens ((^.)) +import qualified Data.Vector.Unboxed as UArr + import Numeric.IEEE +import qualified GHC.Exts as GHC +import qualified Data.Type.Coercion as GHC + -- $linmapIntro -- This library deals with linear functions, i.e. functions @f :: v -> w@ -- that fulfill @@ -218,7 +230,6 @@ - -- $metricIntro -- A norm is a way to quantify the magnitude/length of different vectors, -- even if they point in different directions. @@ -366,6 +377,8 @@ -- @ -- ('euclideanNorm' '<$|' v) '<.>^' w ≡ v '<.>' w -- @ +-- +-- See also '|&>'. (<$|) :: LSpace v => Norm v -> v -> DualVector v Norm m <$| v = m-+$>v @@ -382,6 +395,16 @@ (|$|) :: (LSpace v, Floating (Scalar v)) => Seminorm v -> v -> Scalar v (|$|) m = sqrt . normSq m +infixl 1 |&> +-- | Flipped, “ket” version of '<$|'. +-- +-- @ +-- v '<.>^' (w |&> 'euclideanNorm') ≡ v '<.>' w +-- @ +(|&>) :: LSpace v => DualVector v -> Variance v -> v +dv |&> Norm m = GHC.sym coerceDoubleDual $ m-+$>dv + + -- | 'spanNorm' / 'spanVariance' are inefficient if the number of vectors -- is similar to the dimension of the space, or even larger than it. -- Use this function to optimise the underlying operator to a dense @@ -391,6 +414,12 @@ DualSpaceWitness -> \(Norm m) -> Norm . arr $ sampleLinearFunction $ m +-- | Like 'densifyNorm', but also perform a “sanity check” to eliminate NaN etc. problems. +wellDefinedNorm :: ∀ v . LinearSpace v => Norm v -> Maybe (Norm v) +wellDefinedNorm = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness + -> \(Norm m) -> Norm <$> wellDefinedVector m + data OrthonormalSystem v = OrthonormalSystem { orthonormalityNorm :: Norm v , orthonormalVectors :: [v] @@ -631,7 +660,7 @@ 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) + | μn^2 > epsilon = (v'^*μn, Just $ sqrt (max 0 $ 1 - μn^2)/μn) | otherwise = (v', Nothing) where v' = combined<$|v @@ -679,3 +708,64 @@ instance (SimpleSpace v, Show (DualVector v)) => Show (Norm v) where showsPrec p n = showParen (p>9) $ ("spanNorm "++) . shows (normSpanningSystem n) + +type LinearShowable v = (Show v, RieszDecomposable v) + + + +convexPolytopeHull :: ∀ v . SimpleSpace v => [v] -> [DualVector v] +convexPolytopeHull vs = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness + -> [dv^/η | (dv,η) <- candidates, all ((<=η) . (dv<.>^)) vs] + where vrv = spanVariance vs + nmv = dualNorm' vrv + candidates :: [(DualVector v, Scalar v)] + candidates = [ (dv, dv<.>^v) | v <- vs + , let dv = nmv<$|v ] + +convexPolytopeRepresentatives :: ∀ v . SimpleSpace v => [DualVector v] -> [v] +convexPolytopeRepresentatives dvs + = [v^/η | ((v,η),dv) <- zip candidates dvs + , all (\(w,ψ) -> dv<.>^w <= ψ) candidates] + where nmv :: Norm v + nmv = spanNorm dvs + vrv = dualNorm nmv + candidates :: [(v, Scalar v)] + candidates = [ (v, dv<.>^v) | dv <- dvs + , let v = dv|&>vrv ] + +linearRegressionW :: ∀ s x m y + . ( LinearSpace x, FiniteDimensional y, SimpleSpace m + , Scalar x ~ s, Scalar y ~ s, Scalar m ~ s, RealFrac' s ) + => Norm y -> (x -> (m +> y)) -> [(x,y)] -> m +linearRegressionW σy modelMap = fst . linearRegressionWVar modelMap . map (second (,σy)) + +linearRegressionWVar :: ∀ s x m y + . ( LinearSpace x, FiniteDimensional y, SimpleSpace m + , Scalar x ~ s, Scalar y ~ s, Scalar m ~ s, RealFrac' s ) + => (x -> (m +> y)) -> [(x, (y, Norm y))] -> (m, [DualVector m]) +linearRegressionWVar = lrw (dualSpaceWitness, dualSpaceWitness) + where lrw :: (DualSpaceWitness y, DualSpaceWitness m) + -> (x -> (m +> y)) -> [(x, (y, Norm y))] -> (m, [DualVector m]) + lrw (DualSpaceWitness, DualSpaceWitness) modelMap dataxy + = ( leastSquareSol, deviations ) + where leastSquareSol = (lfun $ forward' . zipWith ((<$|) . snd . snd) dataxy + . forward) + \$ forward' [σy<$|y | (_,(y,σy)) <- dataxy] + forward :: m -> [y] + forward m = [modelMap x $ m | (x,_)<-dataxy] + forward' :: [DualVector y] -> DualVector m + forward' = sumV . zipWith ($) modelGens + modelGens :: [DualVector y +> DualVector m] + modelGens = ((adjoint$) . modelMap . fst)<$>dataxy + deviations = [ m $ dy ^/ ψ + | (m,(dy,ψ)) <- zip modelGens ddys + , ψ > 0 + ] + ddys = [ (dy, ψ) | (x,(yd,σy)) <- dataxy + , let ym = modelMap x $ leastSquareSol + δy = yd ^-^ ym + dy = σy<$|δy + ψ = dy<.>^δy + ] + 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 @@ -22,6 +22,7 @@ {-# LANGUAGE TupleSections #-} {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE GADTs #-} +{-# LANGUAGE DefaultSignatures #-} module Math.LinearMap.Category.Class where @@ -79,6 +80,11 @@ => (v ⊗ w) -+> (v ⊗ w) tensorProduct :: (TensorSpace w, Scalar w ~ Scalar v) => Bilinear v w (v ⊗ w) + tensorProducts :: (TensorSpace w, Scalar w ~ Scalar v) + => [(v,w)] -> (v ⊗ w) + tensorProducts vws = sumV [ getLinearFunction ( + getLinearFunction tensorProduct v) w + | (v,w) <- vws ] transposeTensor :: (TensorSpace w, Scalar w ~ Scalar v) => (v ⊗ w) -+> (w ⊗ v) fmapTensor :: (TensorSpace w, TensorSpace x, Scalar w ~ Scalar v, Scalar x ~ Scalar v) @@ -88,6 +94,14 @@ => Bilinear ((w,x) -+> u) (v⊗w, v⊗x) (v⊗u) coerceFmapTensorProduct :: Hask.Functor p => p v -> Coercion a b -> Coercion (TensorProduct v a) (TensorProduct v b) + -- | “Sanity-check” a vector. This typically amounts to detecting any NaN components, + -- which should trigger a @Nothing@ result. Otherwise, the result should be @Just@ + -- the input, but may also be optimised / memoised if applicable (i.e. for + -- function spaces). + wellDefinedVector :: v -> Maybe v + default wellDefinedVector :: Eq v => v -> Maybe v + wellDefinedVector v = if v==v then Just v else Nothing + wellDefinedTensor :: (TensorSpace w, Scalar w ~ Scalar v) => v⊗w -> Maybe (v⊗w) infixl 7 ⊗ @@ -215,6 +229,8 @@ fmapTensor = biConst0 fzipTensorWith = biConst0 coerceFmapTensorProduct _ Coercion = Coercion + wellDefinedVector Origin = Just Origin + wellDefinedTensor (Tensor Origin) = Just (Tensor Origin) instance Num' s => LinearSpace (ZeroDim s) where type DualVector (ZeroDim s) = ZeroDim s dualSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of @@ -445,6 +461,9 @@ ( coerceFmapTensorProduct (fst<$>p) cab , coerceFmapTensorProduct (snd<$>p) cab ) of (Coercion, Coercion) -> Coercion + wellDefinedVector (u,v) = liftA2 (,) (wellDefinedVector u) (wellDefinedVector v) + wellDefinedTensor (Tensor (u,v)) + = liftA2 ((Tensor.) . (,)) (wellDefinedTensor u) (wellDefinedTensor v) instance ∀ u v . ( LinearSpace u, LinearSpace v, Scalar u ~ Scalar v ) => LinearSpace (u,v) where type DualVector (u,v) = (DualVector u, DualVector v) @@ -612,6 +631,10 @@ cftlp DualSpaceWitness _ c = coerceFmapTensorProduct ([]::[DualVector u]) (fmap c :: Coercion (v⊗a) (v⊗b)) + wellDefinedVector = case dualSpaceWitness :: DualSpaceWitness u of + DualSpaceWitness -> arr asTensor >>> wellDefinedTensor >>> arr (fmap fromTensor) + wellDefinedTensor + = arr hasteLinearMap >>> wellDefinedVector >>> arr (fmap deferLinearMap) -- | @((u+>v)+>w) -> u⊗(v+>w)@ coCurryLinearMap :: ∀ s u v w . ( LinearSpace u, Scalar u ~ s @@ -739,6 +762,8 @@ (TensorProduct u (Tensor s v b)) cftlp _ c = coerceFmapTensorProduct ([]::[u]) (fmap c :: Coercion (v⊗a) (v⊗b)) + wellDefinedVector = wellDefinedTensor + wellDefinedTensor = arr rassocTensor >>> wellDefinedTensor >>> arr (fmap lassocTensor) 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) = LinearMap s u (DualVector v) @@ -919,6 +944,14 @@ ScalarSpaceWitness -> bilinearFunction $ \f (g,h) -> fromLinearFn $ f . ((asLinearFn$g)&&&(asLinearFn$h)) coerceFmapTensorProduct _ Coercion = Coercion + wellDefinedVector = arr sampleLinearFunction >>> wellDefinedVector + >>> fmap (arr applyLinear) + wellDefinedTensor = arr asLinearFn >>> (. applyLinear) + >>> getLinearFunction sampleLinearFunction + >>> wellDefinedVector + >>> fmap (arr fromLinearFn <<< \m + -> sampleLinearFunction + >>> getLinearFunction applyLinear m) exposeLinearFn :: Coercion (LinearMap s (LinearFunction s u v) w) (LinearFunction s (LinearFunction s u v) w) @@ -991,3 +1024,11 @@ (.+^) = (^+^) +-- | Use a function as a linear map. This is only well-defined if the function /is/ +-- linear (this condition is not checked). +lfun :: ( EnhancedCat f (LinearFunction s) + , LinearSpace u, TensorSpace v, Scalar u ~ s, Scalar v ~ s + , Object f u, Object f v ) => (u->v) -> f u v +lfun = arr . LinearFunction + + diff --git a/Math/LinearMap/Category/Derivatives.hs b/Math/LinearMap/Category/Derivatives.hs new file mode 100644 --- /dev/null +++ b/Math/LinearMap/Category/Derivatives.hs @@ -0,0 +1,74 @@ +-- | +-- Module : Math.LinearMap.Category.Derivatives +-- Copyright : (c) Justus Sagemüller 2016 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE PatternSynonyms #-} +{-# LANGUAGE ViewPatterns #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE DefaultSignatures #-} + +module Math.LinearMap.Category.Derivatives + {-# WARNING "These lenses will probably change their domain in the future." #-} where + +import Data.VectorSpace +import Data.VectorSpace.Free + +import Prelude () +import qualified Prelude as Hask + +import Control.Category.Constrained.Prelude +import Control.Arrow.Constrained + +import Data.Type.Coercion +import Data.Tagged + +import Math.Manifold.Core.PseudoAffine +import Math.LinearMap.Asserted +import Math.LinearMap.Category.Instances +import Math.LinearMap.Category.Class + +import Control.Lens + +infixr 7 *∂, /∂, .∂ +(/∂) :: ∀ s x y v q + . ( Num' s, LinearSpace x, LinearSpace y, LinearSpace v, LinearSpace q + , s ~ Scalar x, s ~ Scalar y, s ~ Scalar v, s ~ Scalar q ) + => Lens' y v -> Lens' x q -> Lens' (LinearMap s x y) (LinearMap s q v) +𝑣/∂𝑞 = lens (\m -> fmap (LinearFunction (^.𝑣)) + $ m . arr (LinearFunction $ \q -> zeroV & 𝑞.~q)) + (\m u -> arr.LinearFunction + $ \x -> (m $ x & 𝑞.~zeroV) + ^+^ (𝑣.~(u $ x^.𝑞) $ m $ zeroV & 𝑞.~(x^.𝑞)) ) + +(*∂) :: ∀ s a q v . ( Num' s, OneDimensional q, LinearSpace q, LinearSpace v + , s ~ Scalar a, s ~ Scalar q, s ~ Scalar v ) + => q -> Lens' a (LinearMap s q v) -> Lens' a v +q*∂𝑚 = lens (\a -> a^.𝑚 $ q) + (\a v -> (a & 𝑚 .~ arr (LinearFunction $ \q' -> v ^* (q'^/!q))) ) + +(.∂) :: ∀ s x z . ( Fractional' s, LinearSpace x, s ~ Scalar x, LinearSpace z, s ~ Scalar z ) + => (∀ w . (LinearSpace w, Scalar w ~ s) => Lens' (TensorProduct x w) w) + -> Lens' x z -> Lens' (SymmetricTensor s x) z +𝑤.∂𝑦 = case closedScalarWitness :: ClosedScalarWitness s of + ClosedScalarWitness -> lens + (\(SymTensor t) + -> (getTensorProduct $ fmap (LinearFunction (^.𝑦)) $ t)^.𝑤 ^* 0.5) + (\(SymTensor (Tensor t)) z -> SymTensor . Tensor $ (𝑤.𝑦.~z^*2) t) + 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 @@ -14,6 +14,7 @@ {-# LANGUAGE TypeOperators #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE CPP #-} {-# LANGUAGE TupleSections #-} @@ -40,6 +41,8 @@ import Data.Foldable (foldl') import Data.VectorSpace.Free +import Data.VectorSpace.Free.FiniteSupportedSequence +import Data.VectorSpace.Free.Sequence as Seq import qualified Linear.Matrix as Mat import qualified Linear.Vector as Mat import qualified Linear.Metric as Mat @@ -47,9 +50,13 @@ , _x, _y, _z, _w ) import Control.Lens ((^.)) +import qualified Data.Vector as Arr +import qualified Data.Vector.Unboxed as UArr + import Math.LinearMap.Asserted import Math.VectorSpace.ZeroDimensional +import qualified GHC.Exts as GHC infixr 7 <.>^ (<.>^) :: LinearSpace v => DualVector v -> v -> Scalar v @@ -78,6 +85,7 @@ fzipTensorWith = LinearFunction $ \f -> follow Tensor <<< f <<< flout Tensor *** flout Tensor coerceFmapTensorProduct _ Coercion = Coercion + wellDefinedTensor (Tensor w) = Tensor <$> wellDefinedVector w instance LinearSpace ℝ where type DualVector ℝ = ℝ dualSpaceWitness = DualSpaceWitness @@ -105,7 +113,7 @@ translateP = Tagged (^+^) }; \ instance Num s => PseudoAffine (V s) where { \ v.-~.w = pure (v^-^w); (.-~!) = (^-^) }; \ -instance ∀ s . Num' s => TensorSpace (V s) where { \ +instance ∀ s . (Num' s, Eq s) => TensorSpace (V s) where { \ type TensorProduct (V s) w = V w; \ scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of{ \ ClosedScalarWitness -> ScalarSpaceWitness}; \ @@ -127,8 +135,10 @@ fzipTensorWith = bilinearFunction $ \ \(LinearFunction f) (Tensor vw, Tensor vx) \ -> Tensor $ liftA2 (curry f) vw vx; \ - coerceFmapTensorProduct _ Coercion = Coercion }; \ -instance ∀ s . Num' s => LinearSpace (V s) where { \ + coerceFmapTensorProduct _ Coercion = Coercion; \ + wellDefinedTensor = getTensorProduct >>> Hask.traverse wellDefinedVector \ + >>> fmap Tensor }; \ +instance ∀ s . (Num' s, Eq s) => LinearSpace (V s) where { \ type DualVector (V s) = V s; \ dualSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of \ {ClosedScalarWitness -> DualSpaceWitness}; \ @@ -246,5 +256,238 @@ +instance (Num' n, UArr.Unbox n) => Semimanifold (FinSuppSeq n) where + type Needle (FinSuppSeq n) = FinSuppSeq n + (.+~^) = (.+^); translateP = Tagged (.+^) + toInterior = pure; fromInterior = id +instance (Num' n, UArr.Unbox n) => PseudoAffine (FinSuppSeq n) where + v.-~.w = Just $ v.-.w; (.-~!) = (.-.) +instance (Num' n, UArr.Unbox n) => TensorSpace (FinSuppSeq n) where + type TensorProduct (FinSuppSeq n) v = [v] + wellDefinedVector (FinSuppSeq v) = FinSuppSeq <$> UArr.mapM wellDefinedVector v + scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness n of + ClosedScalarWitness -> ScalarSpaceWitness + linearManifoldWitness = LinearManifoldWitness BoundarylessWitness + zeroTensor = Tensor [] + toFlatTensor = LinearFunction $ Tensor . UArr.toList . getFiniteSeq + fromFlatTensor = LinearFunction $ FinSuppSeq . UArr.fromList . getTensorProduct + addTensors (Tensor s) (Tensor t) = Tensor $ Mat.liftU2 (^+^) s t + scaleTensor = bilinearFunction $ \μ (Tensor t) -> Tensor $ (μ*^)<$>t + negateTensor = LinearFunction $ \(Tensor t) -> Tensor $ negateV<$>t + tensorProduct = bilinearFunction + $ \(FinSuppSeq v) w -> Tensor $ (*^w)<$>UArr.toList v + transposeTensor = LinearFunction $ \(Tensor a) + -> let n = length a + in foldl' (^+^) zeroV + $ zipWith ( \i w -> getLinearFunction tensorProduct w $ basisValue i ) + [0..] a + fmapTensor = bilinearFunction $ \f (Tensor a) -> Tensor $ map (f$) a + fzipTensorWith = bilinearFunction $ \f (Tensor a, Tensor b) + -> Tensor $ zipWith (curry $ arr f) a b + coerceFmapTensorProduct _ Coercion = Coercion + wellDefinedTensor (Tensor a) = Tensor <$> Hask.traverse wellDefinedVector a + + +instance (Num' n, UArr.Unbox n) => Semimanifold (Sequence n) where + type Needle (Sequence n) = Sequence n + (.+~^) = (.+^); translateP = Tagged (.+^) + toInterior = pure; fromInterior = id + +instance (Num' n, UArr.Unbox n) => PseudoAffine (Sequence n) where + v.-~.w = Just $ v.-.w; (.-~!) = (.-.) + +instance (Num' n, UArr.Unbox n) => TensorSpace (Sequence n) where + type TensorProduct (Sequence n) v = [v] + wellDefinedVector (SoloChunk n c) = SoloChunk n <$> UArr.mapM wellDefinedVector c + wellDefinedVector (Sequence h r) = Sequence <$> UArr.mapM wellDefinedVector h + <*> wellDefinedVector r + wellDefinedTensor (Tensor a) = Tensor <$> Hask.traverse wellDefinedVector a + scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness n of + ClosedScalarWitness -> ScalarSpaceWitness + linearManifoldWitness = LinearManifoldWitness BoundarylessWitness + zeroTensor = Tensor [] + toFlatTensor = LinearFunction $ Tensor . GHC.toList + fromFlatTensor = LinearFunction $ GHC.fromList . getTensorProduct + addTensors (Tensor s) (Tensor t) = Tensor $ Mat.liftU2 (^+^) s t + scaleTensor = bilinearFunction $ \μ (Tensor t) -> Tensor $ (μ*^)<$>t + negateTensor = LinearFunction $ \(Tensor t) -> Tensor $ negateV<$>t + tensorProduct = bilinearFunction + $ \v w -> Tensor $ (*^w)<$>GHC.toList v + transposeTensor = LinearFunction $ \(Tensor a) + -> let n = length a + in foldl' (^+^) zeroV + $ zipWith (\i w -> (getLinearFunction tensorProduct w) $ basisValue i) + [0..] a + fmapTensor = bilinearFunction $ \f (Tensor a) -> Tensor $ map (f$) a + fzipTensorWith = bilinearFunction $ \f (Tensor a, Tensor b) + -> Tensor $ zipWith (curry $ arr f) a b + coerceFmapTensorProduct _ Coercion = Coercion + +instance (Num' n, UArr.Unbox n) => LinearSpace (Sequence n) where + type DualVector (Sequence n) = FinSuppSeq n + dualSpaceWitness = case closedScalarWitness :: ClosedScalarWitness n of + ClosedScalarWitness -> DualSpaceWitness + linearId = LinearMap [basisValue i | i<-[0..]] + tensorId = LinearMap [asTensor $ fmap (LinearFunction $ + \w -> Tensor $ replicate (i-1) zeroV ++ [w]) $ id | i<-[0..]] + applyDualVector = bilinearFunction $ adv Seq.minimumChunkSize + where adv _ (FinSuppSeq v) (Seq.SoloChunk o q) + = UArr.sum $ UArr.zipWith (*) (UArr.drop o v) q + adv chunkSize (FinSuppSeq v) (Sequence c r) + | UArr.length v > chunkSize + = UArr.sum (UArr.zipWith (*) v c) + + adv (chunkSize*2) (FinSuppSeq $ UArr.drop chunkSize v) r + | otherwise = UArr.sum $ UArr.zipWith (*) v c + applyLinear = bilinearFunction $ apl Seq.minimumChunkSize + where apl _ (LinearMap m) (Seq.SoloChunk o q) + = sumV $ zipWith (*^) (UArr.toList q) (drop o m) + apl chunkSize (LinearMap m) (Sequence c r) + | null mr = sumV $ zipWith (*^) (UArr.toList c) mc + | otherwise = foldl' (^+^) (apl (chunkSize*2) (LinearMap mr) r) + (zipWith (*^) (UArr.toList c) mc) + where (mc, mr) = splitAt chunkSize m + applyTensorFunctional = bilinearFunction + $ \(LinearMap m) (Tensor t) -> sum $ zipWith (<.>^) m t + applyTensorLinMap = bilinearFunction $ arr curryLinearMap >>> + \(LinearMap m) (Tensor t) + -> sumV $ zipWith (getLinearFunction . getLinearFunction applyLinear) m t +instance (Num' n, UArr.Unbox n) => LinearSpace (FinSuppSeq n) where + type DualVector (FinSuppSeq n) = Sequence n + dualSpaceWitness = case closedScalarWitness :: ClosedScalarWitness n of + ClosedScalarWitness -> DualSpaceWitness + linearId = LinearMap [basisValue i | i<-[0..]] + tensorId = LinearMap [asTensor $ fmap (LinearFunction $ + \w -> Tensor $ replicate (i-1) zeroV ++ [w]) $ id | i<-[0..]] + applyDualVector = bilinearFunction $ adv Seq.minimumChunkSize + where adv _ (Seq.SoloChunk o q) (FinSuppSeq v) + = UArr.sum $ UArr.zipWith (*) q (UArr.drop o v) + adv chunkSize (Sequence c r) (FinSuppSeq v) + | UArr.length v > chunkSize + = UArr.sum (UArr.zipWith (*) c v) + + adv (chunkSize*2) r (FinSuppSeq $ UArr.drop chunkSize v) + | otherwise = UArr.sum $ UArr.zipWith (*) c v + applyLinear = bilinearFunction $ \(LinearMap m) (FinSuppSeq v) + -> foldl' (^+^) zeroV $ zipWith (*^) (UArr.toList v) m + applyTensorFunctional = bilinearFunction + $ \(LinearMap m) (Tensor t) -> sum $ zipWith (<.>^) m t + applyTensorLinMap = bilinearFunction $ arr curryLinearMap >>> + \(LinearMap m) (Tensor t) + -> sumV $ zipWith (getLinearFunction . getLinearFunction applyLinear) m t + + + +instance GHC.IsList (Tensor s (Sequence s) v) where + type Item (Tensor s (Sequence s) v) = v + fromList = Tensor + toList = getTensorProduct + +instance GHC.IsList (Tensor s (FinSuppSeq s) v) where + type Item (Tensor s (FinSuppSeq s) v) = v + fromList = Tensor + toList = getTensorProduct + + + +newtype SymmetricTensor s v + = SymTensor { getSymmetricTensor :: Tensor s v v } +deriving instance (Show (Tensor s v v)) => Show (SymmetricTensor s v) + +instance (TensorSpace v, Scalar v ~ s) => AffineSpace (SymmetricTensor s v) where + type Diff (SymmetricTensor s v) = SymmetricTensor s v + (.+^) = (^+^) + (.-.) = (^-^) +instance (TensorSpace v, Scalar v ~ s) => AdditiveGroup (SymmetricTensor s v) where + SymTensor s ^+^ SymTensor t = SymTensor $ s ^+^ t + zeroV = SymTensor zeroV + negateV (SymTensor t) = SymTensor $ negateV t + +instance (TensorSpace v, Scalar v ~ s) + => VectorSpace (SymmetricTensor s v) where + type Scalar (SymmetricTensor s v) = s + μ *^ SymTensor f = SymTensor $ μ*^f + +instance (TensorSpace v, Scalar v ~ s) => Semimanifold (SymmetricTensor s v) where + type Needle (SymmetricTensor s v) = SymmetricTensor s v + (.+~^) = (^+^) + fromInterior = id + toInterior = pure + translateP = Tagged (^+^) +instance (TensorSpace v, Scalar v ~ s) => PseudoAffine (SymmetricTensor s v) where + (.-~!) = (^-^) +instance (Num' s, TensorSpace v, Scalar v ~ s) => TensorSpace (SymmetricTensor s v) where + type TensorProduct (SymmetricTensor s v) x = Tensor s v (Tensor s v x) + wellDefinedVector (SymTensor t) = SymTensor <$> wellDefinedVector t + scalarSpaceWitness = case closedScalarWitness :: ClosedScalarWitness s of + ClosedScalarWitness -> ScalarSpaceWitness + linearManifoldWitness = LinearManifoldWitness BoundarylessWitness + zeroTensor = Tensor zeroV + toFlatTensor = case closedScalarWitness :: ClosedScalarWitness s of + ClosedScalarWitness -> LinearFunction $ \(SymTensor t) + -> Tensor $ fmap toFlatTensor $ t + fromFlatTensor = case closedScalarWitness :: ClosedScalarWitness s of + ClosedScalarWitness -> LinearFunction $ \(Tensor t) + -> SymTensor $ fmap fromFlatTensor $ t + addTensors (Tensor f) (Tensor g) = Tensor $ f^+^g + subtractTensors (Tensor f) (Tensor g) = Tensor $ f^-^g + negateTensor = LinearFunction $ \(Tensor f) -> Tensor $ negateV f + scaleTensor = bilinearFunction $ \μ (Tensor f) -> Tensor $ μ *^ f + tensorProduct = bilinearFunction $ \(SymTensor t) g + -> Tensor $ fmap (LinearFunction (⊗g)) $ t + transposeTensor = LinearFunction $ \(Tensor f) -> getLinearFunction ( + arr (fmap Coercion) . transposeTensor . arr lassocTensor) f + fmapTensor = bilinearFunction $ \f (Tensor t) -> Tensor $ fmap (fmap f) $ t + fzipTensorWith = bilinearFunction $ \f (Tensor s, Tensor t) + -> Tensor $ fzipWith (fzipWith f) $ (s,t) + coerceFmapTensorProduct _ crc = fmap (fmap crc) + wellDefinedTensor (Tensor t) = Tensor <$> wellDefinedVector t + +instance (Num' s, LinearSpace v, Scalar v ~ s) => LinearSpace (SymmetricTensor s v) where + type DualVector (SymmetricTensor s v) = SymmetricTensor s (DualVector v) + dualSpaceWitness = case ( closedScalarWitness :: ClosedScalarWitness s + , dualSpaceWitness :: DualSpaceWitness v ) of + (ClosedScalarWitness, DualSpaceWitness) -> DualSpaceWitness + linearId = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> LinearMap $ rassocTensor . asTensor + . fmap (follow SymTensor . asTensor) $ id + tensorId = LinearMap $ asTensor . fmap asTensor . curryLinearMap + . fmap asTensor + . curryLinearMap + . fmap (follow $ \t -> Tensor $ rassocTensor $ t) + $ id + applyLinear = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> bilinearFunction $ \(LinearMap f) (SymTensor t) + -> (getLinearFunction applyLinear + $ fromTensor . deferLinearMap . asLinearMap $ f) $ t + applyDualVector = bilinearFunction $ \(SymTensor f) (SymTensor v) + -> getLinearFunction + (getLinearFunction applyDualVector $ fromTensor $ f) v + applyTensorFunctional = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> bilinearFunction $ \(LinearMap f) (Tensor t) + -> getLinearFunction + (getLinearFunction applyTensorFunctional + $ fromTensor . fmap fromTensor $ f) t + applyTensorLinMap = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> bilinearFunction $ \(LinearMap (Tensor f)) (Tensor t) + -> getLinearFunction (getLinearFunction applyTensorLinMap + $ uncurryLinearMap + . fmap (uncurryLinearMap . fromTensor . fmap fromTensor) + $ LinearMap f) t + + + + +squareV :: (Num' s, s ~ Scalar v) + => TensorSpace v => v -> SymmetricTensor s v +squareV v = SymTensor $ v⊗v + +squareVs :: (Num' s, s ~ Scalar v) + => TensorSpace v => [v] -> SymmetricTensor s v +squareVs vs = SymTensor $ tensorProducts [(v,v) | v<-vs] + + +type v⊗〃+>w = LinearMap (Scalar v) (SymmetricTensor (Scalar v) v) w + +currySymBilin :: LinearSpace v => (v⊗〃+>w) -+> (v+>(v+>w)) +currySymBilin = LinearFunction . arr $ fmap fromTensor . fromTensor . flout LinearMap diff --git a/Math/VectorSpace/Docile.hs b/Math/VectorSpace/Docile.hs --- a/Math/VectorSpace/Docile.hs +++ b/Math/VectorSpace/Docile.hs @@ -21,6 +21,8 @@ {-# LANGUAGE TupleSections #-} {-# LANGUAGE LambdaCase #-} {-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE EmptyCase #-} module Math.VectorSpace.Docile where @@ -29,23 +31,27 @@ import Math.LinearMap.Asserted import Data.Tree (Tree(..), Forest) -import Data.List (sortBy, foldl') +import Data.List (sortBy, foldl', tails) import qualified Data.Set as Set import Data.Set (Set) import Data.Ord (comparing) import Data.List (maximumBy, unfoldr) import qualified Data.Vector as Arr import Data.Foldable (toList) +import Data.List (transpose) import Data.Semigroup import Data.VectorSpace import Data.Basis +import Data.Void + import Prelude () import qualified Prelude as Hask import Control.Category.Constrained.Prelude hiding ((^)) import Control.Arrow.Constrained +import Control.Monad.Trans.State import Linear ( V0(V0), V1(V1), V2(V2), V3(V3), V4(V4) , _x, _y, _z, _w, ex, ey, ez, ew ) @@ -54,7 +60,7 @@ import Math.VectorSpace.ZeroDimensional import qualified Linear.Matrix as Mat import qualified Linear.Vector as Mat -import Control.Lens ((^.)) +import Control.Lens ((^.), Lens', lens, ReifiedLens', ReifiedLens(..)) import Data.Coerce import Numeric.IEEE @@ -87,6 +93,25 @@ tensorDualBasisCandidates :: (SemiInner w, Scalar w ~ Scalar v) => [(Int, v⊗w)] -> Forest (Int, DualVector (v⊗w)) + + symTensorDualBasisCandidates + :: [(Int, SymmetricTensor (Scalar v) v)] + -> Forest (Int, SymmetricTensor (Scalar v) (DualVector v)) + + symTensorTensorDualBasisCandidates :: ∀ w . (SemiInner w, Scalar w ~ Scalar v) + => [(Int, SymmetricTensor (Scalar v) v ⊗ w)] + -> Forest (Int, SymmetricTensor (Scalar v) v +> DualVector w) + -- Delegate to the transposed tensor. This is a hack that will sooner or + -- later catch up with us. TODO: make a proper implementation. + symTensorTensorDualBasisCandidates + = case ( dualSpaceWitness :: DualSpaceWitness v + , dualSpaceWitness :: DualSpaceWitness w + , scalarSpaceWitness :: ScalarSpaceWitness v ) of + (DualSpaceWitness, DualSpaceWitness, ScalarSpaceWitness) + -> map (second $ getLinearFunction transposeTensor) + >>> dualBasisCandidates + >>> fmap (fmap . second $ + arr asTensor >>> arr transposeTensor >>> arr fromTensor) cartesianDualBasisCandidates :: [DualVector v] -- ^ Set of canonical basis functionals. @@ -121,11 +146,13 @@ instance (Fractional' s, SemiInner s) => SemiInner (ZeroDim s) where dualBasisCandidates _ = [] tensorDualBasisCandidates _ = [] + symTensorDualBasisCandidates _ = [] instance (Fractional' s, SemiInner s) => SemiInner (V0 s) where dualBasisCandidates _ = [] tensorDualBasisCandidates _ = [] + symTensorDualBasisCandidates _ = [] -orthonormaliseDuals :: ∀ v . (SemiInner v, LSpace v, RealFrac' (Scalar v)) +orthonormaliseDuals :: ∀ v . (SemiInner v, RealFrac' (Scalar v)) => Scalar v -> [(v, DualVector v)] -> [(v,Maybe (DualVector v))] orthonormaliseDuals = od dualSpaceWitness @@ -144,7 +171,7 @@ ovl₀ = v'₀<.>^v ovl₁ = v'₁<.>^v -dualBasis :: ∀ v . (SemiInner v, LSpace v, RealFrac' (Scalar v)) +dualBasis :: ∀ v . (SemiInner v, RealFrac' (Scalar v)) => [v] -> [Maybe (DualVector v)] dualBasis vs = snd <$> result where zip' ((i,v):vs) ((j,v'):ds) @@ -206,6 +233,27 @@ lookupArr = Arr.fromList vs n = Arr.length lookupArr + +zipTravWith :: Hask.Traversable t => (a->b->c) -> t a -> [b] -> Maybe (t c) +zipTravWith f = evalStateT . Hask.traverse zp + where zp a = do + bs <- get + case bs of + [] -> StateT $ const Nothing + (b:bs') -> put bs' >> return (f a b) + +embedFreeSubspace :: ∀ v t r . (SemiInner v, RealFrac' (Scalar v), Hask.Traversable t) + => t v -> Maybe (ReifiedLens' v (t (Scalar v))) +embedFreeSubspace vs = fmap (\(g,s) -> Lens (lens g s)) result + where vsList = toList vs + result = fmap (genGet&&&genSet) . sequenceA $ dualBasis vsList + genGet vsDuals u = case zipTravWith (\_v dv -> dv<.>^u) vs vsDuals of + Just cs -> cs + genSet vsDuals u coefs = case zipTravWith (,) coefs $ zip vsList vsDuals of + Just updators -> foldl' (\ur (c,(v,v')) -> ur ^+^ v^*(c - v'<.>^ur)) + u updators + + instance SemiInner ℝ where dualBasisCandidates = fmap ((`Node`[]) . second recip) . sortBy (comparing $ negate . abs . snd) @@ -213,6 +261,9 @@ tensorDualBasisCandidates = map (second getTensorProduct) >>> dualBasisCandidates >>> fmap (fmap $ second LinearMap) + symTensorDualBasisCandidates = map (second getSymmetricTensor) + >>> dualBasisCandidates + >>> fmap (fmap $ second (arr asTensor >>> SymTensor)) instance (Fractional' s, Ord s, SemiInner s) => SemiInner (V1 s) where dualBasisCandidates = fmap ((`Node`[]) . second recip) @@ -221,28 +272,84 @@ tensorDualBasisCandidates = map (second $ \(Tensor (V1 w)) -> w) >>> dualBasisCandidates >>> fmap (fmap . second $ LinearMap . V1) + symTensorDualBasisCandidates = map (second getSymmetricTensor) + >>> dualBasisCandidates + >>> fmap (fmap $ second (arr asTensor >>> SymTensor)) 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) + symTensorDualBasisCandidates = cartesianDualBasisCandidates + (SymTensor . Tensor<$>[ V2 (V2 1 0) zeroV + , V2 (V2 0 sqrt¹₂) (V2 sqrt¹₂ 0) + , V2 zeroV (V2 0 1)]) + (\(SymTensor (Tensor (V2 (V2 xx xy) + (V2 yx yy)))) + -> abs <$> [xx, (xy+yx)*sqrt¹₂, yy]) + where sqrt¹₂ = sqrt 0.5 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) + symTensorDualBasisCandidates = cartesianDualBasisCandidates + (SymTensor . Tensor<$>[ V3 (V3 1 0 0) zeroV zeroV + , V3 (V3 0 sqrt¹₂ 0) (V3 sqrt¹₂ 0 0) zeroV + , V3 (V3 0 0 sqrt¹₂) zeroV (V3 sqrt¹₂ 0 0) + , V3 zeroV (V3 0 1 0) zeroV + , V3 zeroV (V3 0 0 sqrt¹₂) (V3 0 sqrt¹₂ 0) + , V3 zeroV zeroV (V3 0 0 1)]) + (\(SymTensor (Tensor (V3 (V3 xx xy xz) + (V3 yx yy yz) + (V3 zx zy zz)))) + -> abs <$> [ xx, (xy+yx)*sqrt¹₂, (xz+zx)*sqrt¹₂ + , yy , (yz+zy)*sqrt¹₂ + , zz ]) + where sqrt¹₂ = sqrt 0.5 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) + symTensorDualBasisCandidates = cartesianDualBasisCandidates + (SymTensor . Tensor<$>[ V4 (V4 1 0 0 0) zeroV zeroV zeroV + , V4 (V4 0 sqrt¹₂ 0 0) (V4 sqrt¹₂ 0 0 0) zeroV zeroV + , V4 (V4 0 0 sqrt¹₂ 0) zeroV (V4 sqrt¹₂ 0 0 0) zeroV + , V4 (V4 0 0 0 sqrt¹₂) zeroV zeroV (V4 sqrt¹₂ 0 0 0) + , V4 zeroV (V4 0 1 0 0) zeroV zeroV + , V4 zeroV (V4 0 0 sqrt¹₂ 0) (V4 0 sqrt¹₂ 0 0) zeroV + , V4 zeroV (V4 0 0 0 sqrt¹₂) zeroV (V4 0 sqrt¹₂ 0 0) + , V4 zeroV zeroV (V4 0 0 1 0) zeroV + , V4 zeroV zeroV (V4 0 0 0 sqrt¹₂) (V4 0 0 sqrt¹₂ 0) + , V4 zeroV zeroV zeroV (V4 0 0 0 1)]) + (\(SymTensor (Tensor (V4 (V4 xx xy xz xw) + (V4 yx yy yz yw) + (V4 zx zy zz zw) + (V4 wx wy wz ww)))) + -> abs <$> [ xx, (xy+yx)*sqrt¹₂, (xz+zx)*sqrt¹₂, (xw+wx)*sqrt¹₂ + , yy , (yz+zy)*sqrt¹₂, (yw+wy)*sqrt¹₂ + , zz , (zw+wz)*sqrt¹₂ + , ww ]) + where sqrt¹₂ = sqrt 0.5 -instance ∀ u v . ( SemiInner u, SemiInner v, Scalar u ~ Scalar v ) => SemiInner (u,v) where +infixl 4 ⊗<$> +(⊗<$>) :: ( Num' s + , Object (LinearFunction s) u + , Object (LinearFunction s) v + , Object (LinearFunction s) w ) + => LinearFunction s v w -> Tensor s u v -> Tensor s u w +f⊗<$>t = fmap f $ t + +instance ∀ u v . ( SemiInner u, SemiInner v, Scalar u ~ Scalar v, Num' (Scalar u) ) + => SemiInner (u,v) where dualBasisCandidates = fmap (\(i,(u,v))->((i,u),(i,v))) >>> unzip >>> dualBasisCandidates *** dualBasisCandidates >>> combineBaseis (dualSpaceWitness,dualSpaceWitness) False mempty - where combineBaseis :: (DualSpaceWitness u, DualSpaceWitness v) -> Bool -> Set Int + where combineBaseis :: (DualSpaceWitness u, DualSpaceWitness v) + -> Bool -- ^ “Bias flag”: iff True, v will be preferred. + -> Set Int -- ^ Set of already-assigned basis indices. -> ( Forest (Int, DualVector u) , Forest (Int, DualVector v) ) -> Forest (Int, (DualVector u, DualVector v)) @@ -263,6 +370,60 @@ : combineBaseis wit True forbidden (bu, abv) combineBaseis wit _ forbidden (bu, []) = combineBaseis wit False forbidden (bu,[]) combineBaseis wit _ forbidden ([], bv) = combineBaseis wit True forbidden ([],bv) + symTensorDualBasisCandidates = fmap (\(i,SymTensor (Tensor (u_uv, v_uv))) + -> ( (i, snd ⊗<$> u_uv) + ,((i, SymTensor $ fst ⊗<$> u_uv) + , (i, SymTensor $ snd ⊗<$> v_uv))) ) + >>> unzip >>> second unzip + >>> dualBasisCandidates *** dualBasisCandidates *** dualBasisCandidates + >>> combineBaseis (dualSpaceWitness,dualSpaceWitness) (Just False) mempty + where combineBaseis :: (DualSpaceWitness u, DualSpaceWitness v) + -> Maybe Bool -- ^ @Just True@: prefer v⊗v, @Nothing@: prefer u⊗v + -> Set Int + -> ( Forest (Int, LinearMap (Scalar u) u (DualVector v)) + ,(Forest (Int, SymmetricTensor (Scalar u) (DualVector u)) + , Forest (Int, SymmetricTensor (Scalar v) (DualVector v))) ) + -> Forest (Int, SymmetricTensor (Scalar u) (DualVector u, DualVector v)) + combineBaseis _ _ _ ([], ([],[])) = [] + combineBaseis wit@(DualSpaceWitness,DualSpaceWitness) + Nothing forbidden + (Node (i, duv) buv' : abuv, (bu, bv)) + | i`Set.member`forbidden + = combineBaseis wit Nothing forbidden (abuv, (bu, bv)) + | otherwise + = Node (i, SymTensor $ Tensor + ( (zeroV&&&id)⊗<$>(asTensor$duv) + , (id&&&zeroV)⊗<$>(transposeTensor$asTensor$duv) ) ) + (combineBaseis wit (Just False) + (Set.insert i forbidden) (buv', (bu, bv))) + : combineBaseis wit Nothing forbidden (abuv, (bu, bv)) + combineBaseis wit Nothing forbidden ([], (bu, bv)) + = combineBaseis wit (Just False) forbidden ([], (bu, bv)) + combineBaseis wit@(DualSpaceWitness,DualSpaceWitness) + (Just False) forbidden + (buv, (Node (i,SymTensor du) bu' : abu, bv)) + | i`Set.member`forbidden + = combineBaseis wit (Just False) forbidden (buv, (abu, bv)) + | otherwise + = Node (i, SymTensor $ Tensor ((id&&&zeroV)⊗<$> du, zeroV)) + (combineBaseis wit (Just True) + (Set.insert i forbidden) (buv, (bu', bv))) + : combineBaseis wit (Just False) forbidden (buv, (abu, bv)) + combineBaseis wit (Just False) forbidden (buv, ([], bv)) + = combineBaseis wit (Just True) forbidden (buv, ([], bv)) + combineBaseis wit@(DualSpaceWitness,DualSpaceWitness) + (Just True) forbidden + (buv, (bu, Node (i,SymTensor dv) bv' : abv)) + | i`Set.member`forbidden + = combineBaseis wit (Just True) forbidden (buv, (bu, abv)) + | otherwise + = Node (i, SymTensor $ Tensor (zeroV, (zeroV&&&id)⊗<$> dv)) + (combineBaseis wit Nothing + (Set.insert i forbidden) (buv, (bu, bv'))) + : combineBaseis wit (Just True) forbidden (buv, (bu, abv)) + combineBaseis wit (Just True) forbidden (buv, (bu, [])) + = combineBaseis wit Nothing forbidden (buv, (bu, [])) + tensorDualBasisCandidates = case scalarSpaceWitness :: ScalarSpaceWitness u of ScalarSpaceWitness -> map (second $ \(Tensor (tu, tv)) -> (tu, tv)) >>> dualBasisCandidates @@ -277,6 +438,12 @@ >>> tensorDualBasisCandidates >>> map (fmap . second $ arr uncurryLinearMap) +instance ∀ s v . ( Num' s, SemiInner v, Scalar v ~ s ) + => SemiInner (SymmetricTensor s v) where + dualBasisCandidates = symTensorDualBasisCandidates + tensorDualBasisCandidates = symTensorTensorDualBasisCandidates + symTensorTensorDualBasisCandidates = case () of {} + instance ∀ s u v . ( LinearSpace u, SemiInner (DualVector u), SemiInner v , Scalar u ~ s, Scalar v ~ s ) => SemiInner (LinearMap s u v) where @@ -363,7 +530,7 @@ uncanonicallyFromDual = id uncanonicallyToDual = id -instance (Num' s, LinearSpace s) => FiniteDimensional (V0 s) where +instance (Num' s, Eq s, LinearSpace s) => FiniteDimensional (V0 s) where data SubBasis (V0 s) = V0Basis entireBasis = V0Basis enumerateSubBasis V0Basis = [] @@ -396,7 +563,7 @@ uncanonicallyToDual = id #define FreeFiniteDimensional(V, VB, dimens, take, give) \ -instance (Num' s, LSpace s) \ +instance (Num' s, Eq s, LSpace s) \ => FiniteDimensional (V s) where { \ data SubBasis (V s) = VB deriving (Show); \ entireBasis = VB; \ @@ -521,6 +688,10 @@ = case decomposeLinMapWithin bu $ curryLinearMap $ muvw of Left (bu', mvwsg) -> let (_, (bv', ws)) = goWith bv id (mvwsg []) id in Left (TensorBasis bu' bv', ws) + Right mvwsg -> let (changed, (bv', ws)) = goWith bv id (mvwsg []) id + in if changed + then Left (TensorBasis bu bv', ws) + else Right ws where (_, goWith) = tensorLinmapDecompositionhelpers recomposeSB (TensorBasis bu bv) = recomposeSBTensor bu bv recomposeSBTensor = rst dualSpaceWitness @@ -590,7 +761,115 @@ deriving instance (Show (SubBasis u), Show (SubBasis v)) => Show (SubBasis (Tensor s u v)) +instance ∀ s v . + ( FiniteDimensional v, Scalar v~s, Scalar (DualVector v)~s + , RealFloat' s ) + => FiniteDimensional (SymmetricTensor s v) where + newtype SubBasis (SymmetricTensor s v) = SymTensBasis (SubBasis v) + entireBasis = SymTensBasis entireBasis + enumerateSubBasis (SymTensBasis b) = do + v:vs <- tails $ enumerateSubBasis b + squareV v + : [ (squareV (v^+^w) ^-^ squareV v ^-^ squareV w) ^* sqrt¹₂ | w <- vs ] + where sqrt¹₂ = sqrt 0.5 + subbasisDimension (SymTensBasis b) = ((n-1)*n)`quot`2 + where n = subbasisDimension b + decomposeLinMap = dclm dualSpaceWitness + where dclm (DualSpaceWitness :: DualSpaceWitness v) (LinearMap f) + = (SymTensBasis bf, rmRedundant 0 . symmetrise $ dlw []) + where rmRedundant _ [] = id + rmRedundant k (row:rest) + = (sclOffdiag (drop k row)++) . rmRedundant (k+1) rest + symmetrise l = zipWith (zipWith (^+^)) lm $ transpose lm + where lm = matr l + matr [] = [] + matr l = case splitAt n l of + (row,rest) -> row : matr rest + n = case subbasisDimension bf of + nbf | nbf == subbasisDimension bf' -> nbf + (LinMapBasis bf bf', dlw) + = decomposeLinMap $ asLinearMap . lassocTensor $ f + sclOffdiag (d:o) = 0.5*^d : ((^*sqrt¹₂)<$>o) + sqrt¹₂ = sqrt 0.5 :: s + recomposeSB = rclm dualSpaceWitness + where rclm (DualSpaceWitness :: DualSpaceWitness v) (SymTensBasis b) ws + = case recomposeSB (TensorBasis b b) + $ mkSym (subbasisDimension b) (repeat id) ws of + (t, remws) -> (SymTensor t, remws) + mkSym _ _ [] = [] + mkSym 0 _ ws = ws + mkSym n (sd₀:sds) ws = let (d:o,rest) = splitAt n ws + oscld = (sqrt 0.5*)<$>o + in sd₀ [] ++ [d] ++ oscld + ++ mkSym (n-1) (zipWith (.) sds $ (:)<$>oscld) rest + recomposeLinMap = rclm dualSpaceWitness + where rclm (DualSpaceWitness :: DualSpaceWitness v) (SymTensBasis b) ws + = case recomposeLinMap (LinMapBasis b b) + $ mkSym (subbasisDimension b) (repeat id) ws of + (f, remws) -> (LinearMap $ rassocTensor . asTensor $ f, remws) + mkSym _ _ [] = [] + mkSym 0 _ ws = ws + mkSym n (sd₀:sds) ws = let (d:o,rest) = splitAt n ws + oscld = (sqrt 0.5*^)<$>o + in sd₀ [] ++ [d] ++ oscld + ++ mkSym (n-1) (zipWith (.) sds $ (:)<$>oscld) rest + recomposeSBTensor = rcst + where rcst :: ∀ w . (FiniteDimensional w, Scalar w ~ s) + => SubBasis (SymmetricTensor s v) -> SubBasis w + -> [s] -> (Tensor s (SymmetricTensor s v) w, [s]) + rcst (SymTensBasis b) bw μs + = case recomposeSBTensor (TensorBasis b b) bw + $ mkSym (subbasisDimension bw) (subbasisDimension b) (repeat id) μs of + (Tensor t, remws) -> ( Tensor $ Tensor t + :: Tensor s (SymmetricTensor s v) w + , remws ) + mkSym _ _ _ [] = [] + mkSym _ 0 _ ws = ws + mkSym nw n (sd₀:sds) ws = let (d:o,rest) = multiSplit nw n ws + oscld = map (sqrt 0.5*)<$>o + in concat (sd₀ []) ++ d ++ concat oscld + ++ mkSym nw (n-1) (zipWith (.) sds $ (:)<$>oscld) rest + recomposeContraLinMap f tenss + = LinearMap . arr (rassocTensor . asTensor) . rcCLM dualSpaceWitness f + $ fmap getSymmetricTensor tenss + where rcCLM :: (Hask.Functor f, LinearSpace w, s~Scalar w) + => DualSpaceWitness v + -> (f s->w) -> f (Tensor s (DualVector v) (DualVector v)) + -> LinearMap s (LinearMap s (DualVector v) v) w + rcCLM DualSpaceWitness f = recomposeContraLinMap f + recomposeContraLinMapTensor = rcCLMT' + where rcCLMT' :: ∀ f u w . (Hask.Functor f, LinearSpace w, s~Scalar w + , FiniteDimensional u, s~Scalar u) + => (f s->w) -> f (SymmetricTensor s v +> DualVector u) + -> (SymmetricTensor s v ⊗ u) +> w + rcCLMT' f tenss + = LinearMap . arr (fmap rassocTensor . rassocTensor . asTensor) + . rcCLMT (dualSpaceWitness, dualSpaceWitness) f + $ fmap getLinearMap tenss + where rcCLMT :: (DualSpaceWitness v, DualSpaceWitness u) + -> (f s->w) -> f (Tensor s (DualVector v) + (Tensor s (DualVector v) (DualVector u))) + -- -> LinearMap s (Tensor s (SymmetricTensor s v) u) w + -- ∼ TensorProduct (LinearMap s (SymmetricTensor s v) (DualVector u)) w + -- ⩵ TensorProduct (SymmetricTensor s (DualVector v)) (DualVector u ⊗ w) + -- ⩵ Tensor s (DualVector v) (DualVector v ⊗ (DualVector u ⊗ w)) + -> LinearMap s (LinearMap s (DualVector v) + (LinearMap s (DualVector v) u)) w + -- ∼ Tensor s (Tensor s (DualVector v) + -- (DualVector v ⊗ DualVector u)) w + -- ∼ Tensor s (DualVector v) + -- (Tensor s (DualVector v ⊗ DualVector u) w) + rcCLMT (DualSpaceWitness, DualSpaceWitness) f = recomposeContraLinMap f + uncanonicallyFromDual = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> LinearFunction + $ \(SymTensor t) -> SymTensor $ arr fromLinearMap . uncanonicallyFromDual $ t + uncanonicallyToDual = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> LinearFunction + $ \(SymTensor t) -> SymTensor $ uncanonicallyToDual . arr asLinearMap $ t + +deriving instance (Show (SubBasis v)) => Show (SubBasis (SymmetricTensor s v)) + instance ∀ s u v . ( LSpace u, FiniteDimensional (DualVector u), FiniteDimensional v , Scalar u~s, Scalar v~s, Scalar (DualVector v)~s, Fractional' (Scalar v) ) @@ -783,6 +1062,7 @@ (ScalarSpaceWitness,DualSpaceWitness) -> \p dv -> showParen (p>0) $ ("().<"++) . showsPrec 7 (sRiesz$dv) +instance Show (LinearMap ℝ (ZeroDim ℝ) ℝ) where showsPrec = showsPrecAsRiesz instance Show (LinearMap ℝ (V0 ℝ) ℝ) where showsPrec = showsPrecAsRiesz instance Show (LinearMap ℝ ℝ ℝ) where showsPrec = showsPrecAsRiesz instance Show (LinearMap ℝ (V1 ℝ) ℝ) where showsPrec = showsPrecAsRiesz @@ -802,6 +1082,8 @@ rieszDecomposition m = map (first Left) (rieszDecomposition $ fst . m) ++ map (first Right) (rieszDecomposition $ snd . m) +instance RieszDecomposable (ZeroDim ℝ) where + rieszDecomposition _ = [] instance RieszDecomposable (V0 ℝ) where rieszDecomposition _ = [] instance RieszDecomposable (V1 ℝ) where @@ -842,7 +1124,9 @@ $ foldr (\(b,dv) -> (" ^+^ "++) . showsPrecBasis ([]::[u]) 7 b . (".<"++) . showsPrec 7 dv) s dvs - + +instance Show (LinearMap s v (ZeroDim s)) where + show _ = "zeroV" instance Show (LinearMap s v (V0 s)) where show _ = "zeroV" instance (FiniteDimensional v, v ~ DualVector v, Scalar v ~ ℝ, Show v) @@ -893,6 +1177,9 @@ showsPrecBasis proxy p (Right by) = showParen (p>9) $ ("Right "++) . showsPrecBasis (snd<$>proxy) 10 by +instance TensorDecomposable (ZeroDim ℝ) where + tensorDecomposition _ = [] + showsPrecBasis _ _ = absurd instance TensorDecomposable (V0 ℝ) where tensorDecomposition _ = [] showsPrecBasis _ _ (Mat.E q) = (V0^.q ++) @@ -1011,9 +1298,17 @@ -- -- But /not/ @(v+>w) -> (w+>v)@, in general (though in a Hilbert space, this too is -- equivalent, via 'riesz' isomorphism). -adjoint :: ∀ v w . (LSpace v, LSpace w, Scalar v ~ Scalar w) +adjoint :: ∀ v w . (LinearSpace v, LinearSpace w, Scalar v ~ Scalar w) => (v +> DualVector w) -+> (w +> DualVector v) adjoint = case ( dualSpaceWitness :: DualSpaceWitness v , dualSpaceWitness :: DualSpaceWitness w ) of (DualSpaceWitness, DualSpaceWitness) -> arr fromTensor . transposeTensor . arr asTensor + + + + +multiSplit :: Int -> Int -> [a] -> ([[a]], [a]) +multiSplit chunkSize 0 l = ([],l) +multiSplit chunkSize nChunks l = case splitAt chunkSize l of + (chunk, rest) -> first (chunk:) $ multiSplit chunkSize (nChunks-1) rest 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.3.0.1 +version: 0.3.2.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 @@ -40,6 +40,7 @@ library exposed-modules: Math.LinearMap.Category Math.VectorSpace.ZeroDimensional + Math.LinearMap.Category.Derivatives other-modules: Math.LinearMap.Category.Class Math.LinearMap.Asserted Math.LinearMap.Category.Instances @@ -50,8 +51,8 @@ constrained-categories >=0.3 && <0.4, containers, vector, tagged, - free-vector-spaces >= 0.1.1 && < 0.2, - linear, lens, + free-vector-spaces >= 0.1.2 && < 0.2, + linear, lens, transformers, manifolds-core >= 0.4 && < 0.5, semigroups, ieee754 >= 0.7 && < 0.9