module Data.LinearMap.HerMetric (
HerMetric, HerMetric'
, metricSq, metricSq', metric, metric', metrics, metrics'
, projector, projector'
, euclideanMetric'
, spanHilbertSubspace
, spanSubHilbertSpace
, IsFreeSpace
, factoriseMetric, factoriseMetric'
, productMetric, productMetric'
, metricAsLength, metricFromLength, metric'AsLength
, transformMetric, transformMetric'
, dualiseMetric, dualiseMetric'
, recipMetric, recipMetric'
, eigenSpan, eigenSpan'
, eigenCoSpan, eigenCoSpan'
, metriScale', metriScale
, adjoint
, extendMetric
, HasMetric
, HasMetric'(..)
, (^<.>)
, MetricScalar
, FiniteDimensional(..)
, Stiefel1(..)
) where
import Data.VectorSpace
import Data.LinearMap
import Data.Basis
import Data.MemoTrie
import Data.Semigroup
import Data.Tagged
import Data.Void
import qualified Data.List as List
import qualified Prelude as Hask
import qualified Control.Applicative as Hask
import qualified Control.Monad as Hask
import Control.Category.Constrained.Prelude hiding ((^))
import Control.Arrow.Constrained
import Data.Manifold.Types.Primitive
import Data.CoNat
import qualified Data.Vector as Arr
import qualified Numeric.LinearAlgebra.HMatrix as HMat
import Data.VectorSpace.FiniteDimensional
import Data.LinearMap.Category
import Data.Embedding
infixr 7 <.>^, ^<.>
newtype HerMetric v = HerMetric {
metricMatrix :: Maybe (HMat.Matrix (Scalar v))
}
matrixMetric :: HasMetric v => HMat.Matrix (Scalar v) -> HerMetric v
matrixMetric = HerMetric . Just
instance (HasMetric v) => AdditiveGroup (HerMetric v) where
zeroV = HerMetric Nothing
negateV (HerMetric m) = HerMetric $ negate <$> m
HerMetric Nothing ^+^ HerMetric n = HerMetric n
HerMetric m ^+^ HerMetric Nothing = HerMetric m
HerMetric (Just m) ^+^ HerMetric (Just n) = HerMetric . Just $ m + n
instance HasMetric v => VectorSpace (HerMetric v) where
type Scalar (HerMetric v) = Scalar v
s *^ (HerMetric m) = HerMetric $ HMat.scale s <$> m
newtype HerMetric' v = HerMetric' {
metricMatrix' :: Maybe (HMat.Matrix (Scalar v))
}
extendMetric :: (HasMetric v, Scalar v~ℝ) => HerMetric v -> v -> HerMetric v
extendMetric (HerMetric Nothing) _ = HerMetric Nothing
extendMetric (HerMetric (Just m)) v
| isInfinite' detm = HerMetric $ Just m
| isInfinite' detmninv = singularMetric
| otherwise = HerMetric $ Just mn
where
(minv, (detm, _)) = HMat.invlndet m
(mn, (detmninv, _)) = HMat.invlndet (minv + HMat.outer vv vv)
vv = asPackedVector v
matrixMetric' :: HasMetric v => HMat.Matrix (Scalar v) -> HerMetric' v
matrixMetric' = HerMetric' . Just
instance (HasMetric v) => AdditiveGroup (HerMetric' v) where
zeroV = HerMetric' Nothing
negateV (HerMetric' m) = HerMetric' $ negate <$> m
HerMetric' Nothing ^+^ HerMetric' n = HerMetric' n
HerMetric' m ^+^ HerMetric' Nothing = HerMetric' m
HerMetric' (Just m) ^+^ HerMetric' (Just n) = matrixMetric' $ m + n
instance HasMetric v => VectorSpace (HerMetric' v) where
type Scalar (HerMetric' v) = Scalar v
s *^ (HerMetric' m) = HerMetric' $ HMat.scale s <$> m
projector :: HasMetric v => DualSpace v -> HerMetric v
projector u = matrixMetric $ HMat.outer uDecomp uDecomp
where uDecomp = asPackedVector u
projector' :: HasMetric v => v -> HerMetric' v
projector' v = matrixMetric' $ HMat.outer vDecomp vDecomp
where vDecomp = asPackedVector v
singularMetric :: forall v . HasMetric v => HerMetric v
singularMetric = matrixMetric $ HMat.scale (1/0) (HMat.ident dim)
where (Tagged dim) = dimension :: Tagged v Int
singularMetric' :: forall v . HasMetric v => HerMetric' v
singularMetric' = matrixMetric' $ HMat.scale (1/0) (HMat.ident dim)
where (Tagged dim) = dimension :: Tagged v Int
metricSq :: HasMetric v => HerMetric v -> v -> Scalar v
metricSq (HerMetric Nothing) _ = 0
metricSq (HerMetric (Just m)) v = vDecomp `HMat.dot` HMat.app m vDecomp
where vDecomp = asPackedVector v
metricSq' :: HasMetric v => HerMetric' v -> DualSpace v -> Scalar v
metricSq' (HerMetric' Nothing) _ = 0
metricSq' (HerMetric' (Just m)) u = uDecomp `HMat.dot` HMat.app m uDecomp
where uDecomp = asPackedVector u
metric :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> v -> Scalar v
metric m = sqrt . metricSq m
metric' :: (HasMetric v, Floating (Scalar v)) => HerMetric' v -> DualSpace v -> Scalar v
metric' m = sqrt . metricSq' m
toDualWith :: HasMetric v => HerMetric v -> v -> DualSpace v
toDualWith (HerMetric Nothing) = const zeroV
toDualWith (HerMetric (Just m)) = fromPackedVector . HMat.app m . asPackedVector
metriScale :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> v -> v
metriScale m v = metric m v *^ v
metriScale' :: (HasMetric v, Floating (Scalar v))
=> HerMetric' v -> DualSpace v -> DualSpace v
metriScale' m v = metric' m v *^ v
metrics :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> [v] -> Scalar v
metrics m vs = sqrt . sum $ metricSq m <$> vs
metrics' :: (HasMetric v, Floating (Scalar v)) => HerMetric' v -> [DualSpace v] -> Scalar v
metrics' m vs = sqrt . sum $ metricSq' m <$> vs
transformMetric :: (HasMetric v, HasMetric w, Scalar v ~ Scalar w)
=> (w :-* v) -> HerMetric v -> HerMetric w
transformMetric _ (HerMetric Nothing) = HerMetric Nothing
transformMetric t (HerMetric (Just m)) = matrixMetric $ tmat HMat.<> m HMat.<> HMat.tr tmat
where tmat = asPackedMatrix t
transformMetric' :: ( HasMetric v, HasMetric w, Scalar v ~ Scalar w )
=> (v :-* w) -> HerMetric' v -> HerMetric' w
transformMetric' _ (HerMetric' Nothing) = HerMetric' Nothing
transformMetric' t (HerMetric' (Just m))
= matrixMetric' $ HMat.tr tmat HMat.<> m HMat.<> tmat
where tmat = asPackedMatrix t
dualiseMetric :: HasMetric v => HerMetric (DualSpace v) -> HerMetric' v
dualiseMetric (HerMetric m) = HerMetric' m
dualiseMetric' :: HasMetric v => HerMetric' v -> HerMetric (DualSpace v)
dualiseMetric' (HerMetric' m) = HerMetric m
recipMetric' :: HasMetric v => HerMetric v -> HerMetric' v
recipMetric' (HerMetric Nothing) = singularMetric'
recipMetric' (HerMetric (Just m))
| isInfinite' detm = singularMetric'
| otherwise = matrixMetric' minv
where (minv, (detm, _)) = HMat.invlndet m
recipMetric :: HasMetric v => HerMetric' v -> HerMetric v
recipMetric (HerMetric' Nothing) = singularMetric
recipMetric (HerMetric' (Just m))
| isInfinite' detm = singularMetric
| otherwise = matrixMetric minv
where (minv, (detm, _)) = HMat.invlndet m
isInfinite' :: (Eq a, Num a) => a -> Bool
isInfinite' x = x==x*2
eigenSpan :: (HasMetric v, Scalar v ~ ℝ) => HerMetric' v -> [v]
eigenSpan (HerMetric' Nothing) = []
eigenSpan (HerMetric' (Just m)) = map fromPackedVector eigSpan
where (μs,vsm) = HMat.eigSH' m
eigSpan = zipWith (HMat.scale . sqrt) (HMat.toList μs) (HMat.toColumns vsm)
eigenSpan' :: (HasMetric v, Scalar v ~ ℝ) => HerMetric v -> [DualSpace v]
eigenSpan' (HerMetric Nothing) = []
eigenSpan' (HerMetric (Just m)) = map fromPackedVector eigSpan
where (μs,vsm) = HMat.eigSH' m
eigSpan = zipWith (HMat.scale . sqrt) (HMat.toList μs) (HMat.toColumns vsm)
eigenCoSpan :: (HasMetric v, Scalar v ~ ℝ) => HerMetric' v -> [DualSpace v]
eigenCoSpan (HerMetric' Nothing) = []
eigenCoSpan (HerMetric' (Just m)) = map fromPackedVector eigSpan
where (μs,vsm) = HMat.eigSH' m
eigSpan = zipWith (HMat.scale . recip . sqrt) (HMat.toList μs) (HMat.toColumns vsm)
eigenCoSpan' :: (HasMetric v, Scalar v ~ ℝ) => HerMetric v -> [v]
eigenCoSpan' (HerMetric Nothing) = []
eigenCoSpan' (HerMetric (Just m)) = map fromPackedVector eigSpan
where (μs,vsm) = HMat.eigSH' m
eigSpan = zipWith (HMat.scale . recip . sqrt) (HMat.toList μs) (HMat.toColumns vsm)
type MetricScalar s = ( SmoothScalar s
, Ord s
)
type HasMetric v = (HasMetric' v, HasMetric' (DualSpace v), DualSpace (DualSpace v) ~ v)
class ( FiniteDimensional v, FiniteDimensional (DualSpace v)
, VectorSpace (DualSpace v), HasBasis (DualSpace v)
, MetricScalar (Scalar v), Scalar v ~ Scalar (DualSpace v)
, Basis v ~ Basis (DualSpace v) )
=> HasMetric' v where
type DualSpace v :: *
type DualSpace v = v
(<.>^) :: DualSpace v -> v -> Scalar v
functional :: (v -> Scalar v) -> DualSpace v
doubleDual :: HasMetric' (DualSpace v) => v -> DualSpace (DualSpace v)
doubleDual' :: HasMetric' (DualSpace v) => DualSpace (DualSpace v) -> v
(^<.>) :: HasMetric v => v -> DualSpace v -> Scalar v
ket ^<.> bra = bra <.>^ ket
euclideanMetric' :: forall v . (HasMetric v, InnerSpace v) => HerMetric v
euclideanMetric' = HerMetric . pure $ HMat.ident n
where (Tagged n) = dimension :: Tagged v Int
instance (MetricScalar k) => HasMetric' (ZeroDim k) where
Origin<.>^Origin = zeroV
functional _ = Origin
doubleDual = id; doubleDual'= id
instance HasMetric' Double where
(<.>^) = (<.>)
functional f = f 1
doubleDual = id; doubleDual'= id
instance ( HasMetric v, HasMetric w, Scalar v ~ Scalar w
) => HasMetric' (v,w) where
type DualSpace (v,w) = (DualSpace v, DualSpace w)
(v,w)<.>^(v',w') = v<.>^v' + w<.>^w'
functional f = (functional $ f . (,zeroV), functional $ f . (zeroV,))
doubleDual = id; doubleDual'= id
instance (SmoothScalar s, Ord s, KnownNat n) => HasMetric' (s^n) where
type DualSpace (s^n) = s^n
(<.>^) = (<.>)
functional = fnal
where fnal :: ∀ s n . (SmoothScalar s, KnownNat n) => (s^n -> s) -> s^n
fnal f = FreeVect . Arr.generate n $
\i -> f . FreeVect . Arr.generate n $ \j -> if i==j then 1 else 0
where Tagged n = theNatN :: Tagged n Int
doubleDual = id; doubleDual'= id
instance (HasMetric v, s~Scalar v) => HasMetric' (FinVecArrRep t v s) where
type DualSpace (FinVecArrRep t v s) = FinVecArrRep t (DualSpace v) s
FinVecArrRep v <.>^ FinVecArrRep w = HMat.dot v w
functional = fnal
where fnal :: ∀ v . HasMetric v =>
(FinVecArrRep t v (Scalar v) -> Scalar v)
-> FinVecArrRep t (DualSpace v) (Scalar v)
fnal f = FinVecArrRep . (n HMat.|>)
$ (f . FinVecArrRep) <$> HMat.toRows (HMat.ident n)
Tagged n = dimension :: Tagged v Int
doubleDual = id; doubleDual'= id
adjoint :: (HasMetric v, HasMetric w, Scalar w ~ Scalar v)
=> (v :-* w) -> DualSpace w :-* DualSpace v
adjoint m = linear $ \w -> functional $ \v
-> w <.>^lapply m v
metrConst :: forall v. (HasMetric v, v ~ DualSpace v, Num (Scalar v))
=> Scalar v -> HerMetric v
metrConst μ = matrixMetric $ HMat.scale μ (HMat.ident dim)
where (Tagged dim) = dimension :: Tagged v Int
instance (HasMetric v, v ~ DualSpace v, Num (Scalar v)) => Num (HerMetric v) where
fromInteger = metrConst . fromInteger
(+) = (^+^)
negate = negateV
HerMetric m * HerMetric n = HerMetric $ liftA2 (HMat.<>) m n
abs = error "abs undefined for HerMetric"
signum = error "signum undefined for HerMetric"
metrNumFun :: (HasMetric v, v ~ Scalar v, v ~ DualSpace v, Num v)
=> (v -> v) -> HerMetric v -> HerMetric v
metrNumFun f (HerMetric Nothing) = matrixMetric . HMat.scalar $ f 0
metrNumFun f (HerMetric (Just m)) = matrixMetric . HMat.scalar . f $ m HMat.! 0 HMat.! 0
instance (HasMetric v, v ~ Scalar v, v ~ DualSpace v, Fractional v)
=> Fractional (HerMetric v) where
fromRational = metrConst . fromRational
recip = metrNumFun recip
instance (HasMetric v, v ~ Scalar v, v ~ DualSpace v, Floating v)
=> Floating (HerMetric v) where
pi = metrConst pi
sqrt = metrNumFun sqrt
exp = metrNumFun exp
log = metrNumFun log
sin = metrNumFun sin
cos = metrNumFun cos
tan = metrNumFun tan
asin = metrNumFun asin
acos = metrNumFun acos
atan = metrNumFun atan
sinh = metrNumFun sinh
cosh = metrNumFun cosh
asinh = metrNumFun asinh
atanh = metrNumFun atanh
acosh = metrNumFun acosh
normaliseWith :: HasMetric v => HerMetric v -> v -> Option v
normaliseWith m v = case metric m v of
0 -> empty
μ -> pure (v ^/ μ)
orthonormalPairsWith :: forall v . HasMetric v => HerMetric v -> [v] -> [(v, DualSpace v)]
orthonormalPairsWith met = mkON
where mkON :: [v] -> [(v, DualSpace v)]
mkON [] = []
mkON (v:vs) = let onvs = mkON vs
v' = List.foldl' (\va (vb,pb) -> va ^-^ vb ^* (pb <.>^ va)) v onvs
p' = toDualWith met v'
in case sqrt (p' <.>^ v') of
0 -> onvs
μ -> (v'^/μ, p'^/μ) : onvs
factoriseMetric :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)
=> HerMetric (v,w) -> (HerMetric v, HerMetric w)
factoriseMetric (HerMetric Nothing) = (HerMetric Nothing, HerMetric Nothing)
factoriseMetric met = (sumV *** sumV) . unzip
$ (projector.fst &&& projector.snd) <$> eigenSpan' met
factoriseMetric' :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)
=> HerMetric' (v,w) -> (HerMetric' v, HerMetric' w)
factoriseMetric' met = (sumV *** sumV) . unzip
$ (projector'.fst &&& projector'.snd) <$> eigenSpan met
productMetric :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)
=> HerMetric v -> HerMetric w -> HerMetric (v,w)
productMetric (HerMetric Nothing) (HerMetric Nothing) = HerMetric Nothing
productMetric (HerMetric (Just mv)) (HerMetric (Just mw))
= HerMetric . Just $ HMat.diagBlock [mv, mw]
productMetric (HerMetric Nothing) (HerMetric (Just mw))
= HerMetric . Just $ HMat.diagBlock [HMat.konst 0 (dv,dv), mw]
where (Tagged dv) = dimension :: Tagged v Int
productMetric (HerMetric (Just mv)) (HerMetric Nothing)
= HerMetric . Just $ HMat.diagBlock [mv, HMat.konst 0 (dw,dw)]
where (Tagged dw) = dimension :: Tagged w Int
productMetric' :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)
=> HerMetric' v -> HerMetric' w -> HerMetric' (v,w)
productMetric' (HerMetric' Nothing) (HerMetric' Nothing) = HerMetric' Nothing
productMetric' (HerMetric' (Just mv)) (HerMetric' (Just mw))
= HerMetric' . Just $ HMat.diagBlock [mv, mw]
productMetric' (HerMetric' Nothing) (HerMetric' (Just mw))
= HerMetric' . Just $ HMat.diagBlock [HMat.konst 0 (dv,dv), mw]
where (Tagged dv) = dimension :: Tagged v Int
productMetric' (HerMetric' (Just mv)) (HerMetric' Nothing)
= HerMetric' . Just $ HMat.diagBlock [mv, HMat.konst 0 (dw,dw)]
where (Tagged dw) = dimension :: Tagged w Int
metricAsLength :: HerMetric ℝ -> ℝ
metricAsLength m = case metricSq m 1 of
o | o > 0 -> recip o
| o < 0 -> error "Metric fails to be positive definite!"
| o == 0 -> error "Trying to use zero metric as length."
metricFromLength :: ℝ -> HerMetric ℝ
metricFromLength = projector . recip
metric'AsLength :: HerMetric' ℝ -> ℝ
metric'AsLength = recip . (`metric'`1)
spanHilbertSubspace :: ∀ s v w
. (HasMetric v, Scalar v ~ s, IsFreeSpace w, Scalar w ~ s)
=> HerMetric v
-> [v]
-> Option (Embedding (Linear s) w v)
spanHilbertSubspace met = emb . orthonormalPairsWith met
where emb onb'
| n'==n = return $ Embedding emb prj . arr identityMatrix
| otherwise = empty
where emb = DenseLinear . HMat.fromColumns $ (asPackedVector . fst) <$> onb
prj = DenseLinear . HMat.fromRows $ (asPackedVector . snd) <$> onb
n' = length onb'
onb = take n onb'
(Tagged n) = theNatN :: Tagged (FreeDimension w) Int
spanSubHilbertSpace :: forall s v w
. (HasMetric v, InnerSpace v, Scalar v ~ s, IsFreeSpace w, Scalar w ~ s)
=> [v]
-> Option (Embedding (Linear s) w v)
spanSubHilbertSpace = spanHilbertSubspace euclideanMetric'
newtype Stiefel1 v = Stiefel1 { getStiefel1N :: DualSpace v }
instance (HasMetric v, Scalar v ~ Double, Show (DualSpace v)) => Show (HerMetric v) where
showsPrec p m
| null eigSp = showString "zeroV"
| otherwise = showParen (p>5)
. foldr1 ((.) . (.(" ^+^ "++)))
$ ((("projector "++).).showsPrec 6)<$>eigSp
where eigSp = eigenSpan' m
instance (HasMetric v, Scalar v ~ Double, Show v) => Show (HerMetric' v) where
showsPrec p m
| null eigSp = showString "zeroV"
| otherwise = showParen (p>5)
. foldr1 ((.) . (.(" ^+^ "++)))
$ ((("projector' "++).).showsPrec 6)<$>eigSp
where eigSp = eigenSpan m