linearmap-category 0.3.0.1 → 0.3.2.0
raw patch · 6 files changed
+766/−22 lines, 6 filesdep +transformersdep ~free-vector-spaces
Dependencies added: transformers
Dependency ranges changed: free-vector-spaces
Files
- Math/LinearMap/Category.hs +96/−6
- Math/LinearMap/Category/Class.hs +41/−0
- Math/LinearMap/Category/Derivatives.hs +74/−0
- Math/LinearMap/Category/Instances.hs +246/−3
- Math/VectorSpace/Docile.hs +305/−10
- linearmap-category.cabal +4/−3
Math/LinearMap/Category.hs view
@@ -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+ ]+
Math/LinearMap/Category/Class.hs view
@@ -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++
+ Math/LinearMap/Category/Derivatives.hs view
@@ -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)+
Math/LinearMap/Category/Instances.hs view
@@ -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
Math/VectorSpace/Docile.hs view
@@ -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
linearmap-category.cabal view
@@ -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