manifolds 0.1.3.1 → 0.1.5.0
raw patch · 11 files changed
+1382/−126 lines, 11 filesdep −vector-algorithms
Dependencies removed: vector-algorithms
Files
- Data/LinearMap/HerMetric.hs +81/−1
- Data/Manifold/Cone.hs +332/−0
- Data/Manifold/Griddable.hs +182/−0
- Data/Manifold/PseudoAffine.hs +183/−38
- Data/Manifold/Riemannian.hs +244/−0
- Data/Manifold/TreeCover.hs +209/−33
- Data/Manifold/Types.hs +22/−48
- Data/Manifold/Types/Primitive.hs +20/−2
- Data/VectorSpace/FiniteDimensional.hs +104/−1
- images/examples/simple-2d-ShadeTree.png binary
- manifolds.cabal +5/−3
Data/LinearMap/HerMetric.hs view
@@ -8,6 +8,7 @@ {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UnicodeSyntax #-} @@ -24,6 +25,10 @@ , spanHilbertSubspace , spanSubHilbertSpace , IsFreeSpace+ -- * One-dimensional axes and product spaces+ , factoriseMetric, factoriseMetric'+ , productMetric, productMetric'+ , metricAsLength, metric'AsLength -- * Utility for metrics , transformMetric, transformMetric' , dualiseMetric, dualiseMetric'@@ -40,6 +45,8 @@ -- * Fundamental requirements , MetricScalar , FiniteDimensional(..)+ -- * Misc+ , Stiefel1(..) ) where @@ -366,6 +373,26 @@ (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 @@ -448,7 +475,51 @@ -spanHilbertSubspace :: forall s v w+-- | Project a metric on each of the factors of a product space. This works by+-- projecting the eigenvectors into both subspaces.+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 = recip . (`metric`1)++metric'AsLength :: HerMetric' ℝ -> ℝ+metric'AsLength = recip . (`metric'`1)+++spanHilbertSubspace :: ∀ s v w . (HasMetric v, Scalar v ~ s, IsFreeSpace w, Scalar w ~ s) => HerMetric v -- ^ Metric to induce the inner product on the Hilbert space. -> [v] -- ^ @n@ linearly independent vectors, to span the subspace @w@.@@ -477,3 +548,12 @@ -> Option (Embedding (Linear s) w v) spanSubHilbertSpace = spanHilbertSubspace euclideanMetric' ++-- | The /n/-th Stiefel manifold is the space of all possible configurations of+-- /n/ orthonormal vectors. In the case /n/ = 1, simply the subspace of normalised+-- vectors, i.e. equivalent to the 'UnitSphere'. Even so, it strictly speaking+-- requires the containing space to be at least metric (if not Hilbert); we would+-- however like to be able to use this concept also in spaces with no inner product,+-- therefore we define this space not as normalised vectors, but rather as all+-- vectors modulo scaling by positive factors.+newtype Stiefel1 v = Stiefel1 { getStiefel1N :: DualSpace v }
+ Data/Manifold/Cone.hs view
@@ -0,0 +1,332 @@+-- |+-- Module : Data.Manifold.Cone+-- Copyright : (c) Justus Sagemüller 2015+-- License : GPL v3+-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de+-- Stability : experimental+-- Portability : portable+-- ++{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE LiberalTypeSynonyms #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE UnicodeSyntax #-}+{-# LANGUAGE MultiWayIf #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE CPP #-}+++module Data.Manifold.Cone where+ +++import Data.List+import qualified Data.Vector.Generic as Arr+import qualified Data.Vector+import Data.Maybe+import Data.Semigroup+import Data.Function (on)+import Data.Fixed++import Data.VectorSpace+import Data.LinearMap+import Data.LinearMap.HerMetric+import Data.MemoTrie (HasTrie(..))+import Data.AffineSpace+import Data.Basis+import Data.Complex hiding (magnitude)+import Data.Void+import Data.Tagged+import Data.Manifold.Types.Primitive++import Data.CoNat+import Data.VectorSpace.FiniteDimensional++import qualified Numeric.LinearAlgebra.HMatrix as HMat++import qualified Prelude+import qualified Control.Applicative as Hask++import Control.Category.Constrained.Prelude hiding ((^))+import Control.Arrow.Constrained+import Control.Monad.Constrained+import Data.Foldable.Constrained++import Data.Manifold.PseudoAffine+import Data.Embedding++++type ConeVecArr m = FinVecArrRep Cℝay (CℝayInterior m) (Scalar (Needle m))+type ConeNeedle m = Needle (ConeVecArr m)+type SConn'dConeVecArr m = FinVecArrRep Cℝay (ℝ, Interior m) ℝ+++class ( Semimanifold m, Semimanifold (Interior (Interior m))+ , Semimanifold (ConeVecArr m)+ , Interior (ConeVecArr m) ~ ConeVecArr m )+ => ConeSemimfd m where+ {-# MINIMAL (fromCℝayInterior | fromCD¹Interior)+ , (toCℝayInterior | toCD¹Interior) #-}+ type CℝayInterior m :: *+ + fromCℝayInterior :: ConeVecArr m -> Cℝay m+ fromCℝayInterior = projCD¹ToCℝay . fromCD¹Interior+ fromCD¹Interior :: ConeVecArr m -> CD¹ m+ fromCD¹Interior = embCℝayToCD¹ . fromCℝayInterior+ + toCℝayInterior :: Cℝay m -> Option (ConeVecArr m)+ toCℝayInterior = toCD¹Interior . embCℝayToCD¹+ toCD¹Interior :: CD¹ m -> Option (ConeVecArr m)+ toCD¹Interior = toCℝayInterior . projCD¹ToCℝay++ ++++instance (ConeSemimfd m) => Semimanifold (Cℝay m) where+ type Needle (Cℝay m) = ConeNeedle m+ type Interior (Cℝay m) = ConeVecArr m+ fromInterior = fromCℝayInterior+ toInterior = toCℝayInterior+ translateP = ctp+ where ctp :: Tagged (Cℝay m) (ConeVecArr m -> ConeNeedle m -> ConeVecArr m)+ ctp = Tagged ctp'+ where Tagged ctp' = translateP+ :: Tagged (ConeVecArr m) (ConeVecArr m -> ConeNeedle m -> ConeVecArr m)+ +instance (ConeSemimfd m) => Semimanifold (CD¹ m) where+ type Needle (CD¹ m) = ConeNeedle m+ type Interior (CD¹ m) = ConeVecArr m+ fromInterior = fromCD¹Interior+ toInterior = toCD¹Interior+ translateP = ctp+ where ctp :: Tagged (CD¹ m) (ConeVecArr m -> ConeNeedle m -> ConeVecArr m)+ ctp = Tagged ctp'+ where Tagged ctp' = translateP+ :: Tagged (ConeVecArr m) (ConeVecArr m -> ConeNeedle m -> ConeVecArr m)++instance (ConeSemimfd m, SmoothScalar (Scalar (Needle m))) => PseudoAffine (Cℝay m) where+ p.-~.i = (.-~.i) =<< toInterior p+instance (ConeSemimfd m, SmoothScalar (Scalar (Needle m))) => PseudoAffine (CD¹ m) where+ p.-~.i = (.-~.i) =<< toInterior p+++instance ConeSemimfd (ZeroDim ℝ) where+ type CℝayInterior (ZeroDim ℝ) = ℝ+ fromCℝayInterior (FinVecArrRep qb) | HMat.size qb == 0 = Cℝay 1 Origin+ | x <- qb HMat.! 0 = Cℝay (bijectℝtoℝplus x) Origin + toCℝayInterior (Cℝay 0 Origin) = Hask.empty+ toCℝayInterior (Cℝay y Origin) = pure . FinVecArrRep $ 1 HMat.|>[bijectℝplustoℝ y]+instance ConeSemimfd ℝ where+ type CℝayInterior ℝ = ℝ²+ fromCℝayInterior (FinVecArrRep qb) = Cℝay (q'+b') (q'-b')+ where [q', b'] = HMat.toList $ HMat.cmap ((/2) . bijectℝtoℝplus) qb+ toCℝayInterior (Cℝay 0 _) = Hask.empty+ toCℝayInterior (Cℝay h x) = pure . FinVecArrRep + . HMat.cmap bijectℝplustoℝ $ HMat.fromList [h+x, h-x]+ fromCD¹Interior (FinVecArrRep qb) = CD¹ (bijectℝplustoIntv $ q'+b') (q'-b')+ where [q', b'] = HMat.toList $ HMat.cmap ((/2) . bijectℝtoℝplus) qb+ toCD¹Interior (CD¹ h x) = pure . FinVecArrRep+ . HMat.cmap bijectℝplustoℝ $ HMat.fromList [h'+x, h'-x]+ where h' = bijectIntvtoℝplus h++instance ConeSemimfd S⁰ where+ type CℝayInterior S⁰ = ℝ+ fromCℝayInterior xa | x>0 = Cℝay x PositiveHalfSphere+ | otherwise = Cℝay (-x) NegativeHalfSphere+ where x = getFinVecArrRep xa HMat.! 0+ toCℝayInterior (Cℝay x PositiveHalfSphere) = return . FinVecArrRep $ HMat.scalar x+ toCℝayInterior (Cℝay x NegativeHalfSphere) = return . FinVecArrRep . HMat.scalar $ -x+ fromCD¹Interior xa | x>0 = CD¹ (bijectℝtoIntv x) PositiveHalfSphere+ | otherwise = CD¹ (-bijectℝtoIntv x) NegativeHalfSphere+ where x = getFinVecArrRep xa HMat.! 0+ toCD¹Interior (CD¹ 1 _) = Hask.empty+ toCD¹Interior (CD¹ x PositiveHalfSphere)+ = return . FinVecArrRep . HMat.scalar $ bijectIntvtoℝ x+ toCD¹Interior (CD¹ x NegativeHalfSphere)+ = return . FinVecArrRep . HMat.scalar $ -bijectℝtoIntv x+++instance ConeSemimfd S¹ where+ type CℝayInterior S¹ = ℝ²+ fromCℝayInterior (FinVecArrRep xy) = Cℝay r (S¹ $ atan2 y x)+ where r = HMat.norm_2 xy+ [x,y] = HMat.toList xy+ toCℝayInterior (Cℝay r (S¹ φ)) = return . FinVecArrRep+ . HMat.scale r $ HMat.fromList [cos φ, sin φ]+ fromCD¹Interior (FinVecArrRep xy) = CD¹ (bijectℝtoIntv r) (S¹ $ atan2 y x)+ where r = HMat.norm_2 xy+ [x,y] = HMat.toList xy+ toCD¹Interior (CD¹ 1 _) = Hask.empty+ toCD¹Interior (CD¹ r (S¹ φ)) = return . FinVecArrRep+ . HMat.scale r' $ HMat.fromList [cos φ, sin φ]+ where r' = bijectIntvtoℝ r+++instance ConeSemimfd S² where+ type CℝayInterior S² = ℝ³+ fromCℝayInterior (FinVecArrRep xyz) = Cℝay r (S² (acos $ z/r) (atan2 y x))+ where r = HMat.norm_2 xyz+ [x,y,z] = HMat.toList xyz+ toCℝayInterior (Cℝay r (S² ϑ φ)) = return . FinVecArrRep+ . HMat.scale r $ HMat.fromList [w*x₀, w*y₀, z₀]+ where x₀ = cos φ; y₀ = sin φ; z₀ = cos ϑ; w = sin ϑ++ +++-- | Products of simply connected spaces.+instance ( PseudoAffine x, PseudoAffine y+ , WithField ℝ HilbertSpace (Interior x), WithField ℝ HilbertSpace (Interior y)+ , LinearManifold (FinVecArrRep Cℝay (ℝ, (Interior x, Interior y)) ℝ)+ ) => ConeSemimfd (x,y) where+ type CℝayInterior (x,y) = (ℝ, (Interior x, Interior y))+ fromCℝayInterior = simplyCncted_fromCℝayInterior+ toCℝayInterior = simplyCncted_toCℝayInterior++instance ( KnownNat n ) => ConeSemimfd (ℝ^n) where+ type CℝayInterior (ℝ^n) = (ℝ, ℝ^n)+ fromCℝayInterior = simplyCncted_fromCℝayInterior+ toCℝayInterior = simplyCncted_toCℝayInterior++instance ( HilbertSpace (FinVecArrRep t v ℝ) ) => ConeSemimfd (FinVecArrRep t v ℝ) where+ type CℝayInterior (FinVecArrRep t v ℝ) = (ℝ, FinVecArrRep t v ℝ)+ fromCℝayInterior = simplyCncted_fromCℝayInterior+ toCℝayInterior = simplyCncted_toCℝayInterior+++ +instance ( WithField ℝ ConeSemimfd x, PseudoAffine (Cℝay x)+ , HilbertSpace (CℝayInterior x)+ , HilbertSpace (FinVecArrRep Cℝay (CℝayInterior x) ℝ)+ ) => ConeSemimfd (CD¹ x) where+ type CℝayInterior (CD¹ x) = (ℝ, ConeVecArr x)+ fromCℝayInterior i = Cℝay h (embCℝayToCD¹ o)+ where (Cℝay h o) = simplyCncted_fromCℝayInterior i+ toCℝayInterior (Cℝay _ (CD¹ 1 _)) = Hask.empty+ toCℝayInterior (Cℝay h p) = simplyCncted_toCℝayInterior $ Cℝay h (projCD¹ToCℝay p)+ + +instance ( WithField ℝ ConeSemimfd x, PseudoAffine (Cℝay x)+ , HilbertSpace (CℝayInterior x)+ , HilbertSpace (FinVecArrRep Cℝay (CℝayInterior x) ℝ)+ ) => ConeSemimfd (Cℝay x) where+ type CℝayInterior (Cℝay x) = (ℝ, ConeVecArr x)+ fromCℝayInterior = simplyCncted_fromCℝayInterior+ toCℝayInterior = simplyCncted_toCℝayInterior+ + +simplyCncted_fromCℝayInterior :: (PseudoAffine x, WithField ℝ HilbertSpace (Interior x))+ => SConn'dConeVecArr x -> Cℝay x+simplyCncted_fromCℝayInterior (FinVecArrRep ri) = Cℝay h . fromInterior . fromPackedVector+ $ subtract (h/n) `Arr.map` Arr.tail cmps+ where h = Arr.sum cmps+ cmps = bijectℝtoℝplus `HMat.cmap` ri+ n = fromIntegral $ Arr.length cmps+ +simplyCncted_toCℝayInterior :: (PseudoAffine x, WithField ℝ HilbertSpace (Interior x))+ => Cℝay x -> Option (SConn'dConeVecArr x)+simplyCncted_toCℝayInterior (Cℝay h v) | h/=0, Option (Just vi) <- toInterior v + = let cmps'' = asPackedVector vi+ cmps' = (+ h/n) `HMat.cmap` cmps''+ cmps = (h - Arr.sum cmps') `Arr.cons` cmps+ n = fromIntegral $ Arr.length cmps+ in return $ FinVecArrRep (bijectℝplustoℝ `Arr.map` cmps)+simplyCncted_toCℝayInterior (Cℝay _ _) = Hask.empty+++-- Some essential homeomorphisms+bijectℝtoℝplus , bijectℝplustoℝ+ , bijectIntvtoℝplus, bijectℝplustoIntv+ , bijectIntvtoℝ, bijectℝtoIntv+ :: ℝ -> ℝ++bijectℝplustoℝ x = x - 1/x+bijectℝtoℝplus y = y/2 + sqrt(y^2/4 + 1)++-- [0, 1[ ↔ ℝ⁺+bijectℝplustoIntv y = 1 - recip (y+1)+bijectIntvtoℝplus x = recip(1-x) - 1++-- ]-1, 1[ ↔ ℝ (Similar to 'tanh', but converges less quickly towards ±1.)+bijectℝtoIntv y | y>0 = -1/(2*y) + sqrt(1/(4*y^2) + 1)+ | y<0 = -1/(2*y) - sqrt(1/(4*y^2) + 1)+ | otherwise = 0+ -- 0 = x² + x/y - 1+ -- x = -1/2y ± sqrt(1/4y² + 1)+bijectIntvtoℝ x = x / (1-x^2)++embCℝayToCD¹ :: Cℝay m -> CD¹ m+embCℝayToCD¹ (Cℝay h m) = CD¹ (bijectℝplustoIntv h) m++projCD¹ToCℝay :: CD¹ m -> Cℝay m+projCD¹ToCℝay (CD¹ h m) = Cℝay (bijectIntvtoℝplus h) m++-- instance (WithScalar ℝ PseudoAffine m) => Semimanifold (Cℝay m) where+-- type Needle (Cℝay m) = (Needle m, ℝ)+-- type Interior (Cℝay m) = (Interior m, ℝ)+-- +-- fromInterior (im, d)+-- | d>38 = Cℝay m d -- from 38 on, the +1 is numerically+-- -- insignificant against the exponential.+-- | otherwise = cℝay m (log $ exp d + 1)+-- -- note that (for the same reason we can shortcut above 38)+-- -- such negative arguments will actually yield the value zero.+-- -- This means we're actually reaching the “infinitely far”+-- -- rim rather quickly. This might be a problem, but normally+-- -- shouldn't really matter much.+-- -- It would perhaps be better to have homeomorphism that+-- -- approaches -1/x in the negative limit, but such a+-- -- function doesn't seem as easy to come by.+-- where m = fromInterior im+-- toInterior (Cℝay m q)+-- | q>38 = fmap (,q) im+-- | q>0 = fmap (, log $ exp d - 1) im+-- | otherwise = Hask.empty+-- where im = toInterior m++stiefel1Project :: LinearManifold v =>+ DualSpace v -- ^ Must be nonzero.+ -> Stiefel1 v+stiefel1Project = Stiefel1++stiefel1Embed :: HilbertSpace v => Stiefel1 v -> v+stiefel1Embed (Stiefel1 n) = normalized n+ ++class (PseudoAffine v, InnerSpace v, NaturallyEmbedded (UnitSphere v) (DualSpace v))+ => HasUnitSphere v where+ type UnitSphere v :: *+ stiefel :: UnitSphere v -> Stiefel1 v+ stiefel = Stiefel1 . embed+ unstiefel :: Stiefel1 v -> UnitSphere v+ unstiefel = coEmbed . getStiefel1N++instance HasUnitSphere ℝ where type UnitSphere ℝ = S⁰+instance HasUnitSphere (FinVecArrRep t ℝ ℝ) where type UnitSphere (FinVecArrRep t ℝ ℝ) = S⁰++instance HasUnitSphere ℝ² where type UnitSphere ℝ² = S¹+instance HasUnitSphere (FinVecArrRep t ℝ² ℝ) where type UnitSphere (FinVecArrRep t ℝ² ℝ) = S¹++instance HasUnitSphere ℝ³ where type UnitSphere ℝ³ = S²+instance HasUnitSphere (FinVecArrRep t ℝ³ ℝ) where type UnitSphere (FinVecArrRep t ℝ³ ℝ) = S²++-- instance (HasUnitSphere v, v ~ DualSpace v) => NaturallyEmbedded (Stiefel1 v) v where+-- embed = embed . unstiefel+-- coEmbed = stiefel . coEmbed++++
+ Data/Manifold/Griddable.hs view
@@ -0,0 +1,182 @@+-- |+-- Module : Data.Manifold.Griddable+-- Copyright : (c) Justus Sagemüller 2015+-- License : GPL v3+-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de+-- Stability : experimental+-- Portability : portable+-- +{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE DeriveFoldable #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE ParallelListComp #-}+{-# LANGUAGE UnicodeSyntax #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE LiberalTypeSynonyms #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DataKinds #-}+++module Data.Manifold.Griddable (GridAxis(..), Griddable(..)) where+++import Data.List hiding (filter, all, elem, sum)+import Data.Maybe+import qualified Data.Map as Map+import qualified Data.Vector as Arr+import Data.List.NonEmpty (NonEmpty(..))+import Data.List.FastNub+import qualified Data.List.NonEmpty as NE+import Data.Semigroup++import Data.VectorSpace+import Data.LinearMap+import Data.LinearMap.HerMetric+import Data.LinearMap.Category+import Data.AffineSpace+import Data.Basis+import Data.Complex hiding (magnitude)+import Data.Void+import Data.Tagged+import Data.Proxy++import Data.SimplicialComplex+import Data.Manifold.Types+import Data.Manifold.Types.Primitive ((^), (^.))+import Data.Manifold.PseudoAffine+import Data.Manifold.TreeCover (Shade(..), fullShade, shadeCtr, shadeExpanse)+ +import Data.Embedding+import Data.CoNat++import qualified Prelude as Hask hiding(foldl, sum, sequence)+import qualified Control.Applicative as Hask+import qualified Control.Monad as Hask hiding(forM_, sequence)+import Data.Functor.Identity+import Control.Monad.Trans.State+import Control.Monad.Trans.Writer+import Control.Monad.Trans.Class+import qualified Data.Foldable as Hask+import Data.Foldable (all, elem, toList, sum)+import qualified Data.Traversable as Hask+import Data.Traversable (forM)++import qualified Numeric.LinearAlgebra.HMatrix as HMat++import Control.Category.Constrained.Prelude hiding+ ((^), all, elem, sum, forM, Foldable(..), Traversable)+import Control.Arrow.Constrained+import Control.Monad.Constrained hiding (forM)+import Data.Foldable.Constrained++import Text.Printf+import GHC.Generics (Generic)+++data GridAxis m g = GridAxInterval (Shade m)+ | GridAxCons (Shade m) g (GridAxis m g)+ | GridAxisClosed g (GridAxis m g)+ deriving (Hask.Functor)++gshmap :: (Shade m -> Shade n) -> GridAxis m g -> GridAxis n g+gshmap f (GridAxInterval i) = GridAxInterval $ f i+gshmap f (GridAxCons i g ax) = GridAxCons (f i) g $ gshmap f ax+gshmap f (GridAxisClosed g ax) = GridAxisClosed g $ gshmap f ax++axisEnumFromStepTo :: (ℝ->a) -> ℝ -> ℝ -> ℝ -> GridAxis ℝ a+axisEnumFromStepTo f l st r+ | l' > r = GridAxInterval $ intvl2Shade (Interval l l')+ | otherwise = GridAxCons (intvl2Shade $ Interval l l')+ (f l') $ axisEnumFromStepTo f l' st r+ where l' = l+st++axisGrLength :: GridAxis m a -> Int+axisGrLength (GridAxInterval _) = 0+axisGrLength (GridAxCons _ _ ax) = 1 + axisGrLength ax+axisGrLength (GridAxisClosed _ ax) = axisGrLength ax++class (WithField ℝ Manifold m) => Griddable m g where+ data GriddingParameters m g :: *+ mkGridding :: GriddingParameters m g -> Int -> Shade m -> [GridAxis m g]+++instance Griddable ℝ String where+ data GriddingParameters ℝ String = ℝGridParam+ mkGridding ℝGridParam n (Shade c expa') = [ax]+ where l = c - expa+ r = c + expa+ + expa = metric'AsLength expa'+ + (Just ax) = find ((>=n) . axisGrLength)+ $ [ let qe = 10^^lqe' * nb+ in axisEnumFromStepTo (prettyFloatShow lqe')+ ( qe * fromIntegral (floor $ l / qe) ) qe r+ | lqe' <- [lqe - 1, lqe - 2 ..], nb <- [5, 2, 1] ]+ + lqe = lqef expa :: Int+ lqef n | n > 0 = floor $ lg n+ | n < 0 = floor $ lg (-n)+++instance (Griddable m a, Griddable n a) => Griddable (m,n) a where+ data GriddingParameters (m,n) a = PairGriddingParameters {+ fstGriddingParams :: GriddingParameters m a+ , sndGriddingParams :: GriddingParameters n a }+ mkGridding (PairGriddingParameters p₁ p₂) n (Shade (c₁,c₂) e₁e₂)+ = gshmap ( uncurry fullShade . ( (,c₂).(^.shadeCtr)+ &&& (`productMetric'`e₂).(^.shadeExpanse)) )+ <$> g₁s+ ++ gshmap ( uncurry fullShade . ( (c₁,).(^.shadeCtr)+ &&& ( productMetric' e₁).(^.shadeExpanse)) )+ <$> g₂s+ where g₁s = mkGridding p₁ n $ fullShade c₁ e₁+ g₂s = mkGridding p₂ n $ fullShade c₂ e₂+ (e₁,e₂) = factoriseMetric' e₁e₂ ++prettyFloatShow :: Int -> Double -> String+prettyFloatShow _ 0 = "0"+prettyFloatShow preci x+ | preci >= 0, preci < 4 = show $ round x+ | preci < 0, preci > -2 = printf "%.1f" x+ | otherwise = case ceiling (0.01 + lg (abs x/10^^(preci+1))) + preci of+ 0 | preci < 0 -> printf ("%."++show(-preci)++"f") x+ expn | expn>preci -> printf ("%."++show(expn-preci)++"f*10^%i")+ (x/10^^expn) expn+ | otherwise -> printf ("%i*10^%i")+ (round $ x/10^^expn :: Int) expn++++data Interval = Interval { ivLBound, ivRBound :: ℝ }++shade2Intvl :: Shade ℝ -> Interval+shade2Intvl sh = Interval l r+ where c = sh ^. shadeCtr+ expa = metric'AsLength $ sh ^. shadeExpanse+ l = c - expa; r = c + expa++intvl2Shade :: Interval -> Shade ℝ+intvl2Shade (Interval l r) = fullShade c (projector' expa)+ where c = (l+r) / 2+ expa = (r-l) / 2+ ++lg :: Floating a => a -> a+lg = logBase 10+
Data/Manifold/PseudoAffine.hs view
@@ -17,6 +17,14 @@ -- diffeomorphic. At the moment, we mainly focus on /region-wise differentiable functions/, -- which are a promising compromise between flexibility of definition and provability of -- analytic properties. In particular, they are well-suited for visualisation purposes.+-- +-- The classes in this module are mostly aimed at manifolds /without boundary/.+-- Manifolds with boundary (which we call @MWBound@, never /manifold/!)+-- are more or less treated as a disjoint sum of the interior and the boundary.+-- To understand how this module works, best first forget about boundaries – in this case,+-- @'Interior' x ~ x@, 'fromInterior' and 'toInterior' are trivial, and+-- '.+~|', '|-~.' and 'betweenBounds' are irrelevant.+-- The manifold structure of the boundary itself is not considered at all here. {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-}@@ -84,10 +92,12 @@ import Data.Manifold.Types.Primitive import Data.CoNat+import Data.VectorSpace.FiniteDimensional import qualified Numeric.LinearAlgebra.HMatrix as HMat import qualified Prelude+import qualified Control.Applicative as Hask import Control.Category.Constrained.Prelude hiding ((^)) import Control.Arrow.Constrained@@ -100,11 +110,13 @@ infix 6 .-~. infixl 6 .+~^, .-~^ -class (AdditiveGroup (Needle x)) => Semimanifold x where+class ( AdditiveGroup (Needle x), Interior (Interior x) ~ Interior x )+ => Semimanifold x where+ {-# MINIMAL ((.+~^) | fromInterior), toInterior, translateP #-} -- | The space of “natural” ways starting from some reference point -- and going to some particular target point. Hence, -- the name: like a compass needle, but also with an actual length.- -- For affine space, 'Needle' is simply the space of+ -- For affine spaces, 'Needle' is simply the space of -- line segments (aka vectors) between two points, i.e. the same as 'Diff'. -- The 'AffineManifold' constraint makes that requirement explicit. -- @@ -112,9 +124,37 @@ -- used somewhat synonymously). type Needle x :: * - -- | Generalised translation operation.- (.+~^) :: x -> Needle x -> x+ -- | Manifolds with boundary are a bit tricky. We support such manifolds,+ -- but carry out most calculations only in “the fleshy part” – the+ -- interior, which is an “infinite space”, so you can arbitrarily scale paths.+ -- + -- The default implementation is @'Interior' x = x@, which corresponds+ -- to a manifold that has no boundary to begin with.+ type Interior x :: *+ type Interior x = x + -- | Generalised translation operation. Note that the result will always also+ -- be in the interior; scaling up the needle can only get you ever /closer/+ -- to a boundary.+ (.+~^) :: Interior x -> Needle x -> x+ (.+~^) = addvp+ where addvp :: ∀ x . Semimanifold x => Interior x -> Needle x -> x+ addvp p = fromInterior . tp p+ where (Tagged tp) = translateP :: Tagged x (Interior x -> Needle x -> Interior x)+ + -- | 'id' sans boundary.+ fromInterior :: Interior x -> x+ fromInterior p = p .+~^ zeroV + + toInterior :: x -> Option (Interior x)+ + -- | The signature of '.+~^' should really be @'Interior' x -> 'Needle' x -> 'Interior' x@,+ -- only, this is not possible because it only consists of non-injective type families.+ -- The solution is this tagged signature, which is of course rather unwieldy. That's+ -- why '.+~^' has the stronger, but easier usable signature. Without boundary, these+ -- functions should be equivalent, i.e. @translateP = Tagged (.+~^)@.+ translateP :: Tagged x (Interior x -> Needle x -> Interior x)+ -- | Shorthand for @\\p v -> p .+~^ 'negateV' v@, which should obey the /asymptotic/ law -- -- @@@ -126,17 +166,19 @@ -- as /O/ (/η/²). For large vectors, it will however behave differently, -- except in flat spaces (where all this should be equivalent to the 'AffineSpace' -- instance).- (.-~^) :: x -> Needle x -> x+ (.-~^) :: Interior x -> Needle x -> x p .-~^ v = p .+~^ negateV v --- | This is the class underlying manifolds. ('Manifold' only adds an extra constraint that--- would be circular if it was in a single class. You can always just use 'Manifold'--- as a constraint in your signatures, but you must /define/ only 'PseudoAffine' for--- manifold types – the 'Manifold' instance follows universally from this.)+ +-- | This is the class underlying manifolds. ('Manifold' only precludes boundaries+-- and adds an extra constraint that would be circular if it was in a single+-- class. You can always just use 'Manifold' as a constraint in your signatures,+-- but you must /define/ only 'PseudoAffine' for manifold types –+-- the 'Manifold' instance follows universally from this, if @'Interior x ~ x@.) -- --- The interface is almost identical to the better-known 'AffineSpace' class, but unlike--- in the mathematical definition of affine spaces we don't require associativity --- of '.+~^' with '^+^' – except in an asymptotic sense for small vectors.+-- The interface is (boundaries aside) almost identical to the better-known+-- 'AffineSpace' class, but we don't require associativity of '.+~^' with '^+^'+-- – except in an /asymptotic sense/ for small vectors. -- -- That innocent-looking change makes the class applicable to vastly more general types: -- while an affine space is basically nothing but a vector space without particularly@@ -146,22 +188,38 @@ -- manifolds in their usual maths definition (with an atlas of charts: a family of -- overlapping regions of the topological space, each homeomorphic to the 'Needle' -- vector space or some simply-connected subset thereof).-class Semimanifold x => PseudoAffine x where+class ( Semimanifold x, Semimanifold (Interior x)+ , Needle (Interior x) ~ Needle x, Interior (Interior x) ~ Interior x)+ => PseudoAffine x where -- | The path reaching from one point to another.- -- Should only yield 'Nothing' if the points are on disjoint segments of a- -- non–path-connected manifold. Otherwise, the identity+ -- Should only yield 'Nothing' if+ -- + -- * The points are on disjoint segments of a non–path-connected space.+ -- + -- * Either of the points is on the boundary. Use '|-~.' to deal with this.+ -- + -- On manifolds, the identity -- -- @ -- p .+~^ (q.-~.p) ≡ q -- @ -- -- should hold, at least save for floating-point precision limits etc..- (.-~.) :: x -> x -> Option (Needle x)+ -- + -- '.-~.' and '.+~^' only really work in manifolds without boundary. If you consider+ -- the path between two points, one of which lies on the boundary, it can't really+ -- be possible to scale this path any longer – it would have to reach “out of the+ -- manifold”. To adress this problem, these functions basically consider only the+ -- /interior/ of the space.+ (.-~.) :: x -> Interior x -> Option (Needle x)+ + + -- | See 'Semimanifold' and 'PseudoAffine' for the methods.-class (PseudoAffine m, LinearManifold (Needle m)) => Manifold m-instance (PseudoAffine m, LinearManifold (Needle m)) => Manifold m+class (PseudoAffine m, LinearManifold (Needle m), Interior m ~ m) => Manifold m+instance (PseudoAffine m, LinearManifold (Needle m), Interior m ~ m) => Manifold m type LocallyScalable s x = ( PseudoAffine x, (Needle x) ~ Needle x , HasMetric (Needle x)@@ -173,7 +231,7 @@ -- -- (Actually, 'LinearManifold' is stronger than 'VectorSpace' at the moment, since -- 'HasMetric' requires 'FiniteDimensional'. This might be lifted in the future.)-type LinearManifold x = ( PseudoAffine x, Needle x ~ x, HasMetric x )+type LinearManifold x = ( PseudoAffine x, Interior x ~ x, Needle x ~ x, HasMetric x ) -- | Require some constraint on a manifold, and also fix the type of the manifold's -- underlying field. For example, @WithField ℝ 'HilbertSpace' v@ constrains@@ -185,12 +243,12 @@ type WithField s c x = ( c x, s ~ Scalar (Needle x) ) -- | The 'RealFloat' class plus manifold constraints.-type RealDimension r = ( PseudoAffine r, Needle r ~ r+type RealDimension r = ( PseudoAffine r, Interior r ~ r, Needle r ~ r , HasMetric r, DualSpace r ~ r, Scalar r ~ r , RealFloat r ) -- | The 'AffineSpace' class plus manifold constraints.-type AffineManifold m = ( PseudoAffine m, AffineSpace m+type AffineManifold m = ( PseudoAffine m, Interior m ~ m, AffineSpace m , Needle m ~ Diff m, LinearManifold (Diff m) ) -- | A Hilbert space is a /complete/ inner product space. Being a vector space, it is@@ -200,7 +258,8 @@ -- but since 'Manifold's are at the moment confined to finite dimension, they are in -- fact (trivially) complete.) type HilbertSpace x = ( LinearManifold x, InnerSpace x- , Needle x ~ x, DualSpace x ~ x, Floating (Scalar x) )+ , Interior x ~ x, Needle x ~ x, DualSpace x ~ x+ , Floating (Scalar x) ) -- | An euclidean space is a real affine space whose tangent space is a Hilbert space. type EuclidSpace x = ( AffineManifold x, InnerSpace (Diff x)@@ -215,46 +274,98 @@ type Metric' x = HerMetric' (Needle x) --- | Interpolate between points, approximately linearly.-palerp :: (PseudoAffine x, VectorSpace (Needle x))- => x -> x -> Option (Scalar (Needle x) -> x)-palerp p1 p2 = fmap (\v t -> p1 .+~^ t *^ v) $ p2 .-~. p1+-- | Interpolate between points, approximately linearly. For+-- points that aren't close neighbours (i.e. lie in an almost+-- flat region), the pathway is basically undefined – save for+-- its end points.+-- +-- A proper, really well-defined (on global scales) interpolation+-- only makes sense on a Riemannian manifold, as geodesics.+-- This is a task to be tackled in the future.+palerp :: ∀ x. Manifold x+ => Interior x -> Interior x -> Option (Scalar (Needle x) -> x)+palerp p1 p2 = case (fromInterior p2 :: x) .-~. p1 of+ Option (Just v) -> return $ \t -> p1 .+~^ t *^ v+ _ -> Hask.empty #define deriveAffine(t) \ instance Semimanifold (t) where { \ type Needle (t) = Diff (t); \- (.+~^) = (.+^) }; \-instance PseudoAffine (t) where { \+ fromInterior = id; \+ toInterior = pure; \+ translateP = Tagged (.+^); \+ (.+~^) = (.+^) }; \+instance PseudoAffine (t) where { \ a.-~.b = pure (a.-.b); } deriveAffine(Double) deriveAffine(Rational) +instance SmoothScalar s => Semimanifold (FinVecArrRep t b s) where+ type Needle (FinVecArrRep t b s) = FinVecArrRep t b s+ type Interior (FinVecArrRep t b s) = FinVecArrRep t b s+ fromInterior = id+ toInterior = pure+ translateP = Tagged (.+^)+ (.+~^) = (.+^)+instance SmoothScalar s => PseudoAffine (FinVecArrRep t b s) where+ a.-~.b = pure (a.-.b)+ + instance Semimanifold (ZeroDim k) where type Needle (ZeroDim k) = ZeroDim k+ fromInterior = id+ toInterior = pure Origin .+~^ Origin = Origin Origin .-~^ Origin = Origin+ translateP = Tagged (.+~^) instance PseudoAffine (ZeroDim k) where Origin .-~. Origin = pure Origin instance (Semimanifold a, Semimanifold b) => Semimanifold (a,b) where type Needle (a,b) = (Needle a, Needle b)+ type Interior (a,b) = (Interior a, Interior b) (a,b).+~^(v,w) = (a.+~^v, b.+~^w) (a,b).-~^(v,w) = (a.-~^v, b.-~^w)+ fromInterior (i,j) = (fromInterior i, fromInterior j)+ toInterior (a,b) = fzip (toInterior a, toInterior b)+ translateP = tp+ where tp :: ∀ a b . (Semimanifold a, Semimanifold b)+ => Tagged (a,b) ( (Interior a, Interior b) + -> (Needle a, Needle b)+ -> (Interior a, Interior b) )+ tp = Tagged $ \(a,b) (v,w) -> (ta a v, tb b w)+ where Tagged ta = translateP :: Tagged a (Interior a -> Needle a -> Interior a)+ Tagged tb = translateP :: Tagged b (Interior b -> Needle b -> Interior b) instance (PseudoAffine a, PseudoAffine b) => PseudoAffine (a,b) where (a,b).-~.(c,d) = liftA2 (,) (a.-~.c) (b.-~.d) instance (Semimanifold a, Semimanifold b, Semimanifold c) => Semimanifold (a,b,c) where type Needle (a,b,c) = (Needle a, Needle b, Needle c)+ type Interior (a,b,c) = (Interior a, Interior b, Interior c) (a,b,c).+~^(v,w,x) = (a.+~^v, b.+~^w, c.+~^x) (a,b,c).-~^(v,w,x) = (a.-~^v, b.-~^w, c.-~^x)+ fromInterior (i,j,k) = (fromInterior i, fromInterior j, fromInterior k)+ toInterior (a,b,c) = liftA3 (,,) (toInterior a) (toInterior b) (toInterior c)+ translateP = tp+ where tp :: ∀ a b v . (Semimanifold a, Semimanifold b, Semimanifold c)+ => Tagged (a,b,c) ( (Interior a, Interior b, Interior c) + -> (Needle a, Needle b, Needle c)+ -> (Interior a, Interior b, Interior c) )+ tp = Tagged $ \(a,b,c) (v,w,x) -> (ta a v, tb b w, tc c x)+ where Tagged ta = translateP :: Tagged a (Interior a -> Needle a -> Interior a)+ Tagged tb = translateP :: Tagged b (Interior b -> Needle b -> Interior b)+ Tagged tc = translateP :: Tagged c (Interior c -> Needle c -> Interior c) instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) => PseudoAffine (a,b,c) where (a,b,c).-~.(d,e,f) = liftA3 (,,) (a.-~.d) (b.-~.e) (c.-~.f) instance (MetricScalar a, KnownNat n) => Semimanifold (FreeVect n a) where type Needle (FreeVect n a) = FreeVect n a+ fromInterior = id+ toInterior = pure+ translateP = Tagged (.+~^) (.+~^) = (.+^) instance (MetricScalar a, KnownNat n) => PseudoAffine (FreeVect n a) where a.-~.b = pure (a.-.b)@@ -262,6 +373,9 @@ instance Semimanifold S⁰ where type Needle S⁰ = ℝ⁰+ fromInterior = id+ toInterior = pure+ translateP = Tagged (.+~^) p .+~^ Origin = p p .-~^ Origin = p instance PseudoAffine S⁰ where@@ -271,6 +385,9 @@ instance Semimanifold S¹ where type Needle S¹ = ℝ+ fromInterior = id+ toInterior = pure+ translateP = Tagged (.+~^) S¹ φ₀ .+~^ δφ | φ' < 0 = S¹ $ φ' + tau | otherwise = S¹ $ φ'@@ -282,8 +399,25 @@ | otherwise = pure δφ where δφ = φ₁ - φ₀ +instance Semimanifold D¹ where+ type Needle D¹ = ℝ+ type Interior D¹ = ℝ+ fromInterior = D¹ . tanh+ toInterior (D¹ x) | abs x < 1 = return $ atanh x+ | otherwise = Hask.empty+ translateP = Tagged (+)+instance PseudoAffine D¹ where+ D¹ 1 .-~. _ = Hask.empty+ D¹ (-1) .-~. _ = Hask.empty+ D¹ x .-~. y+ | abs x < 1 = return $ atanh x - y+ | otherwise = Hask.empty+ instance Semimanifold S² where type Needle S² = ℝ²+ fromInterior = id+ toInterior = pure+ translateP = Tagged (.+~^) S² ϑ₀ φ₀ .+~^ δv | ϑ₀ < pi/2 = sphereFold PositiveHalfSphere $ ϑ₀*^embed(S¹ φ₀) ^+^ δv | otherwise = sphereFold NegativeHalfSphere $ (pi-ϑ₀)*^embed(S¹ φ₀) ^+^ δv@@ -304,6 +438,9 @@ instance Semimanifold ℝP² where type Needle ℝP² = ℝ²+ fromInterior = id+ toInterior = pure+ translateP = Tagged (.+~^) ℝP² r₀ φ₀ .+~^ (δr, δφ) | r₀ > 1/2 = case r₀ + δr of r₁ | r₁ > 1 -> ℝP² (2-r₁) (toS¹range $ φ₀+δφ+pi)@@ -323,16 +460,16 @@ | otherwise = pure ( r₁*^embed(S¹ φ₁) ^-^ r₀*^embed(S¹ φ₀) ) -instance (PseudoAffine m, VectorSpace (Needle m), Scalar (Needle m) ~ ℝ)- => Semimanifold (CD¹ m) where- type Needle (CD¹ m) = (Needle m, ℝ)- CD¹ h₀ m₀ .+~^ (h₁δm, δh)- = let h₁ = min 1 . max 1e-300 $ h₀+δh; δm = h₁δm^/h₁- in CD¹ h₁ (m₀.+~^δm)-instance (PseudoAffine m, VectorSpace (Needle m), Scalar (Needle m) ~ ℝ)- => PseudoAffine (CD¹ m) where- CD¹ h₁ m₁ .-~. CD¹ h₀ m₀- = fmap ( \δm -> (h₁*^δm, h₁-h₀) ) $ m₁.-~.m₀+-- instance (PseudoAffine m, VectorSpace (Needle m), Scalar (Needle m) ~ ℝ)+-- => Semimanifold (CD¹ m) where+-- type Needle (CD¹ m) = (Needle m, ℝ)+-- CD¹ h₀ m₀ .+~^ (h₁δm, δh)+-- = let h₁ = min 1 . max 1e-300 $ h₀+δh; δm = h₁δm^/h₁+-- in CD¹ h₁ (m₀.+~^δm)+-- instance (PseudoAffine m, VectorSpace (Needle m), Scalar (Needle m) ~ ℝ)+-- => PseudoAffine (CD¹ m) where+-- CD¹ h₁ m₁ .-~. CD¹ h₀ m₀+-- = fmap ( \δm -> (h₁*^δm, h₁-h₀) ) $ m₁.-~.m₀ @@ -403,6 +540,9 @@ in (z, f'*.*g', devfg) +instance (RealDimension s) => EnhancedCat (->) (Differentiable s) where+ arr (Differentiable f) x = let (y,_,_) = f x in y+ instance (MetricScalar s) => Cartesian (Differentiable s) where type UnitObject (Differentiable s) = ZeroDim s swap = Differentiable $ \(x,y) -> ((y,x), lSwap, const zeroV)@@ -694,6 +834,11 @@ instance (RealDimension s) => EnhancedCat (PWDiffable s) (Differentiable s) where arr = globalDiffable+instance (RealDimension s) => EnhancedCat (->) (PWDiffable s) where+ arr (PWDiffable g) x = let (_,Differentiable f) = g x+ (y,_,_) = f x + in y+ instance (RealDimension s) => Cartesian (PWDiffable s) where type UnitObject (PWDiffable s) = ZeroDim s
+ Data/Manifold/Riemannian.hs view
@@ -0,0 +1,244 @@+-- |+-- Module : Data.Manifold.Riemannian+-- Copyright : (c) Justus Sagemüller 2015+-- License : GPL v3+-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de+-- Stability : experimental+-- Portability : portable+-- +-- Riemannian manifolds are manifolds equipped with a 'Metric' at each point.+-- That means, these manifolds aren't merely topological objects anymore, but+-- have a geometry as well. This gives, in particular, a notion of distance+-- and shortest paths (geodesics) along which you can interpolate.+-- +-- Keep in mind that the types in this library are+-- generally defined in an abstract-mathematical spirit, which may not always+-- match the intuition if you think about manifolds as embedded in ℝ³.+-- (For instance, the torus inherits its geometry from the decomposition as+-- @'S¹' × 'S¹'@, not from the “doughnut” embedding; the cone over @S¹@ is+-- simply treated as the unit disk, etc..)++{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE DeriveFoldable #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE ParallelListComp #-}+{-# LANGUAGE UnicodeSyntax #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE LiberalTypeSynonyms #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds #-}+++module Data.Manifold.Riemannian where+++import Data.List hiding (filter, all, elem, sum)+import Data.Maybe+import qualified Data.Map as Map+import qualified Data.Vector as Arr+import Data.List.NonEmpty (NonEmpty(..))+import Data.List.FastNub+import qualified Data.List.NonEmpty as NE+import Data.Semigroup+import Data.Ord (comparing)+import Control.DeepSeq++import Data.VectorSpace+import Data.LinearMap+import Data.LinearMap.HerMetric+import Data.LinearMap.Category+import Data.AffineSpace+import Data.Basis+import Data.Complex hiding (magnitude)+import Data.Void+import Data.Tagged+import Data.Proxy++import Data.Manifold.Types+import Data.Manifold.Types.Primitive ((^), embed, coEmbed)+import Data.Manifold.PseudoAffine+import Data.VectorSpace.FiniteDimensional+ +import Data.Embedding+import Data.CoNat++import qualified Prelude as Hask hiding(foldl, sum, sequence)+import qualified Control.Applicative as Hask+import qualified Control.Monad as Hask hiding(forM_, sequence)+import Data.Functor.Identity+import Control.Monad.Trans.State+import Control.Monad.Trans.Writer+import Control.Monad.Trans.Class+import qualified Data.Foldable as Hask+import Data.Foldable (all, elem, toList, sum)+import qualified Data.Traversable as Hask+import Data.Traversable (forM)++import qualified Numeric.LinearAlgebra.HMatrix as HMat++import Control.Category.Constrained.Prelude hiding+ ((^), all, elem, sum, forM, Foldable(..), Traversable)+import Control.Arrow.Constrained+import Control.Monad.Constrained hiding (forM)+import Data.Foldable.Constrained++import GHC.Generics (Generic)+++class PseudoAffine x => Geodesic x where+ geodesicBetween ::+ x -- ^ Starting point; the interpolation will yield this at -1.+ -> x -- ^ End point, for +1.+ -- + -- If the two points are actually connected by a path...+ -> Option (D¹ -> x) -- ^ ...then this is the interpolation function. Attention: + -- the type will change to 'Differentiable' in the future.++interpolate :: (Geodesic x, IntervalLike i) => x -> x -> Option (i -> x)+interpolate a b = (. toClosedInterval) <$> geodesicBetween a b+++++#define deriveAffineGD(x) \+instance Geodesic x where { \+ geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹ \+ }++deriveAffineGD (ℝ)++instance Geodesic (ZeroDim ℝ) where+ geodesicBetween Origin Origin = return $ \_ -> Origin++instance (Geodesic a, Geodesic b) => Geodesic (a,b) where+ geodesicBetween (a,b) (α,β) = liftA2 (&&&) (geodesicBetween a α) (geodesicBetween b β)++instance (Geodesic a, Geodesic b, Geodesic c) => Geodesic (a,b,c) where+ geodesicBetween (a,b,c) (α,β,γ)+ = liftA3 (\ia ib ic t -> (ia t, ib t, ic t))+ (geodesicBetween a α) (geodesicBetween b β) (geodesicBetween c γ)++instance (KnownNat n) => Geodesic (FreeVect n ℝ) where+ geodesicBetween (FreeVect v) (FreeVect w)+ = return $ \(D¹ t) -> let μv = (1-t)/2; μw = (t+1)/2+ in FreeVect $ Arr.zipWith (\vi wi -> μv*vi + μw*wi) v w++instance (PseudoAffine v) => Geodesic (FinVecArrRep t v ℝ) where+ geodesicBetween (FinVecArrRep v) (FinVecArrRep w)+ | HMat.size v>0 && HMat.size w>0+ = return $ \(D¹ t) -> let μv = (1-t)/2; μw = (t+1)/2+ in FinVecArrRep $ HMat.scale μv v + HMat.scale μw w++instance (Geodesic v, WithField ℝ HilbertSpace v)+ => Geodesic (Stiefel1 v) where+ geodesicBetween (Stiefel1 p') (Stiefel1 q')+ = (\f -> \(D¹ t) -> Stiefel1 . f . D¹ $ g * tan (ϑ*t))+ <$> geodesicBetween p q+ where p = normalized p'; q = normalized q'+ l = magnitude $ p^-^q+ ϑ = asin $ l/2+ g = sqrt $ 4/l^2 - 1+++instance Geodesic S⁰ where+ geodesicBetween PositiveHalfSphere PositiveHalfSphere = return $ const PositiveHalfSphere+ geodesicBetween NegativeHalfSphere NegativeHalfSphere = return $ const NegativeHalfSphere+ geodesicBetween _ _ = Hask.empty++instance Geodesic S¹ where+ geodesicBetween (S¹ φ) (S¹ ϕ)+ | abs (φ-ϕ) < pi = (>>> S¹) <$> geodesicBetween φ ϕ+ | φ > 0 = (>>> S¹ . \ψ -> signum ψ*pi - ψ)+ <$> geodesicBetween (pi-φ) (-ϕ-pi)+ | otherwise = (>>> S¹ . \ψ -> signum ψ*pi - ψ)+ <$> geodesicBetween (-pi-φ) (pi-ϕ)+++instance Geodesic (Cℝay S⁰) where+ geodesicBetween p q = (>>> fromℝ) <$> geodesicBetween (toℝ p) (toℝ q)+ where toℝ (Cℝay h PositiveHalfSphere) = h+ toℝ (Cℝay h NegativeHalfSphere) = -h+ fromℝ x | x>0 = Cℝay x PositiveHalfSphere+ | otherwise = Cℝay (-x) NegativeHalfSphere++instance Geodesic (CD¹ S⁰) where+ geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q)+ where toI (CD¹ h PositiveHalfSphere) = h+ toI (CD¹ h NegativeHalfSphere) = -h+ fromI x | x>0 = CD¹ x PositiveHalfSphere+ | otherwise = CD¹ (-x) NegativeHalfSphere++instance Geodesic (Cℝay S¹) where+ geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q)+ where fromP = fromInterior+ toP w = case toInterior w of {Option (Just i) -> i}++instance Geodesic (CD¹ S¹) where+ geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q)+ where toI (CD¹ h (S¹ φ)) = (h*cos φ, h*sin φ)+ fromI (x,y) = CD¹ (sqrt $ x^2+y^2) (S¹ $ atan2 y x)++instance Geodesic (Cℝay S²) where+ geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q)+ where fromP = fromInterior+ toP w = case toInterior w of {Option (Just i) -> i}++instance Geodesic (CD¹ S²) where+ geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q :: ℝ³)+ where toI (CD¹ h sph) = h *^ embed sph+ fromI v = CD¹ (magnitude v) (coEmbed v)++#define geoVSpCone(c,t) \+instance (c) => Geodesic (Cℝay (t)) where { \+ geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q) \+ where { fromP (x,0) = Cℝay 0 x \+ ; fromP (x,h) = Cℝay h (x^/h) \+ ; toP (Cℝay h w) = ( h*^w, h ) } } ; \+instance (c) => Geodesic (CD¹ (t)) where { \+ geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q) \+ where { fromP (x,0) = CD¹ 0 x \+ ; fromP (x,h) = CD¹ h (x^/h) \+ ; toP (CD¹ h w) = ( h*^w, h ) } }++geoVSpCone ((), ℝ)+geoVSpCone ((), ℝ⁰)+geoVSpCone ((WithField ℝ HilbertSpace a, WithField ℝ HilbertSpace b, Geodesic (a,b)), (a,b))+geoVSpCone (KnownNat n, FreeVect n ℝ)+geoVSpCone ((Geodesic v, WithField ℝ HilbertSpace v), FinVecArrRep t v ℝ)+++++-- | One-dimensional manifolds, whose closure is homeomorpic to the unit interval.+class WithField ℝ PseudoAffine i => IntervalLike i where+ toClosedInterval :: i -> D¹ -- Differentiable ℝ i D¹++instance IntervalLike D¹ where+ toClosedInterval = id+instance IntervalLike (CD¹ S⁰) where+ toClosedInterval (CD¹ h PositiveHalfSphere) = D¹ h+ toClosedInterval (CD¹ h NegativeHalfSphere) = D¹ (-h)+instance IntervalLike (Cℝay S⁰) where+ toClosedInterval (Cℝay h PositiveHalfSphere) = D¹ $ tanh h+ toClosedInterval (Cℝay h NegativeHalfSphere) = D¹ $ -tanh h+instance IntervalLike (CD¹ ℝ⁰) where+ toClosedInterval (CD¹ h Origin) = D¹ $ h*2 - 1+instance IntervalLike (Cℝay ℝ⁰) where+ toClosedInterval (Cℝay h Origin) = D¹ $ 1 - 2/(h+1)+instance IntervalLike ℝ where+ toClosedInterval x = D¹ $ tanh x+
Data/Manifold/TreeCover.hs view
@@ -34,7 +34,9 @@ module Data.Manifold.TreeCover ( -- * Shades - Shade, shadeCtr, shadeExpanse, fullShade, pointsShades+ Shade(..), Shade'(..)+ -- ** Lenses and constructors+ , shadeCtr, shadeExpanse, shadeNarrowness, fullShade, fullShade', pointsShades -- * Shade trees , ShadeTree(..), fromLeafPoints -- * Simple view helpers@@ -42,14 +44,14 @@ -- ** Auxiliary types , SimpleTree, Trees, NonEmptyTree, GenericTree(..) -- * Misc- , sShSaw, chainsaw, HasFlatView(..)+ , sShSaw, chainsaw, HasFlatView(..), shadesMerge, smoothInterpolate -- ** Triangulation-builders , TriangBuild, doTriangBuild, singleFullSimplex, autoglueTriangulation , AutoTriang, elementaryTriang, breakdownAutoTriang ) where -import Data.List hiding (filter, all, elem, sum)+import Data.List hiding (filter, all, elem, sum, foldr1) import Data.Maybe import qualified Data.Map as Map import qualified Data.Vector as Arr@@ -87,14 +89,14 @@ import Control.Monad.Trans.Writer import Control.Monad.Trans.Class import qualified Data.Foldable as Hask-import Data.Foldable (all, elem, toList, sum)+import Data.Foldable (all, elem, toList, sum, foldr1) import qualified Data.Traversable as Hask import Data.Traversable (forM) import qualified Numeric.LinearAlgebra.HMatrix as HMat import Control.Category.Constrained.Prelude hiding- ((^), all, elem, sum, forM, Foldable(..), Traversable)+ ((^), all, elem, sum, forM, Foldable(..), foldr1, Traversable) import Control.Arrow.Constrained import Control.Monad.Constrained hiding (forM) import Data.Foldable.Constrained@@ -116,17 +118,48 @@ -- -- For a /precise/ description of an arbitrarily-shaped connected subset of a manifold, -- there is 'Region', whose implementation is vastly more complex.-data Shade x = Shade { shadeCtr :: !x- , shadeExpanse :: !(Metric' x) }+data Shade x = Shade { _shadeCtr :: !(Interior x)+ , _shadeExpanse :: !(Metric' x) } +-- | A “co-shade” can describe ellipsoid regions as well, but unlike+-- 'Shade' it can be unlimited / infinitely wide in some directions.+-- It does OTOH need to have nonzero thickness, which 'Shade' needs not.+data Shade' x = Shade' { _shade'Ctr :: !(Interior x)+ , _shade'Narrowness :: !(Metric x) }++class IsShade shade where+-- type (*) shade :: *->*+ -- | Access the center of a 'Shade' or a 'Shade''.+ shadeCtr :: Functor f (->) (->) => (Interior x->f (Interior x)) -> shade x -> f (shade x)+-- -- | Convert between 'Shade' and 'Shade' (which must be neither singular nor infinite).+-- unsafeDualShade :: WithField ℝ Manifold x => shade x -> shade* x++instance IsShade Shade where+ shadeCtr f (Shade c e) = fmap (`Shade`e) $ f c++shadeExpanse :: Functor f (->) (->) => (Metric' x -> f (Metric' x)) -> Shade x -> f (Shade x)+shadeExpanse f (Shade c e) = fmap (Shade c) $ f e++instance IsShade Shade' where+ shadeCtr f (Shade' c e) = fmap (`Shade'`e) $ f c++shadeNarrowness :: Functor f (->) (->) => (Metric x -> f (Metric x)) -> Shade' x -> f (Shade' x)+shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e+ instance (AffineManifold x) => Semimanifold (Shade x) where type Needle (Shade x) = Diff x+ fromInterior = id+ toInterior = pure+ translateP = Tagged (.+~^) Shade c e .+~^ v = Shade (c.+^v) e Shade c e .-~^ v = Shade (c.-^v) e fullShade :: WithField ℝ Manifold x => x -> Metric' x -> Shade x fullShade ctr expa = Shade ctr expa +fullShade' :: WithField ℝ Manifold x => x -> Metric x -> Shade' x+fullShade' ctr expa = Shade' ctr expa+ subshadeId' :: WithField ℝ Manifold x => x -> NonEmpty (DualSpace (Needle x)) -> x -> (Int, HourglassBulb) subshadeId' c expvs x = case x .-~. c of@@ -151,6 +184,7 @@ pointsShades :: WithField ℝ Manifold x => [x] -> [Shade x] pointsShades = map snd . pointsShades' zeroV + pseudoECM :: WithField ℝ Manifold x => NonEmpty x -> (x, ([x],[x])) pseudoECM (p₀ NE.:| psr) = foldl' ( \(acc, (rb,nr)) (i,p) -> case p.-~.acc of @@ -171,22 +205,52 @@ <$> mapM (.-~.ctr) ps -minusLogOcclusion :: (PseudoAffine x, HasMetric (Needle x)- , s ~ (Scalar (Needle x)), RealDimension s )+-- | Attempt to reduce the number of shades to fewer (ideally, a single one).+-- In the simplest cases these should guaranteed cover the same area;+-- for non-flat manifolds it only works in a heuristic sense.+shadesMerge :: WithField ℝ Manifold x+ => ℝ -- ^ How near (inverse normalised distance, relative to shade expanse)+ -- two shades must be to be merged. If this is zero, any shades+ -- in the same connected region of a manifold are merged.+ -> [Shade x] -- ^ A list of /n/ shades.+ -> [Shade x] -- ^ /m/ ≤ /n/ shades which cover at least the same area.+shadesMerge fuzz (sh₁@(Shade c₁ e₁) : shs) = case extractJust tryMerge shs of+ (Just mg₁, shs') -> shadesMerge fuzz+ $ shs'++[mg₁] -- Append to end to prevent undue weighting+ -- of first shade and its mergers.+ (_, shs') -> sh₁ : shadesMerge fuzz shs' + where tryMerge (Shade c₂ e₂)+ | Option (Just v) <- c₁.-~.c₂+ , Option (Just v') <- c₂.-~.c₁+ , [e₁',e₂'] <- recipMetric<$>[e₁, e₂] + , b₁ <- metric e₂' v+ , b₂ <- metric e₁' v+ , fuzz*b₁*b₂ <= b₁ + b₂+ = Just $ let cc = c₂ .+~^ v ^/ 2+ Option (Just cv₁) = c₁.-~.cc+ Option (Just cv₂) = c₂.-~.cc+ in Shade cc . sumV $ [e₁, e₂] ++ projector'<$>[cv₁, cv₂] + | otherwise = Nothing+shadesMerge _ shs = shs++minusLogOcclusion :: ( Manifold x, s ~ (Scalar (Needle x)), RealDimension s ) => Shade x -> x -> s minusLogOcclusion (Shade p₀ δ) = occ where occ p = case p .-~. p₀ of- Option(Just vd) -> metricSq δinv vd+ Option(Just vd) | mSq <- metricSq δinv vd+ , mSq == mSq -- avoid NaN+ -> mSq _ -> 1/0 δinv = recipMetric δ -- | Check the statistical likelyhood of a point being within a shade.-occlusion :: (PseudoAffine x, HasMetric (Needle x)- , s ~ (Scalar (Needle x)), RealDimension s )+occlusion :: ( Manifold x, s ~ (Scalar (Needle x)), RealDimension s ) => Shade x -> x -> s occlusion (Shade p₀ δ) = occ where occ p = case p .-~. p₀ of- Option(Just vd) -> exp . negate $ metricSq δinv vd+ Option(Just vd) | mSq <- metricSq δinv vd+ , mSq == mSq -- avoid NaN+ -> exp (negate mSq) _ -> zeroV δinv = recipMetric δ @@ -215,11 +279,6 @@ flipHour :: Hourglass s -> Hourglass s flipHour (Hourglass u l) = Hourglass l u -newtype Hourglasses s = Hourglasses {- getHourglasses :: NonEmpty (Hourglass s) }- deriving (Generic, Hask.Functor, Hask.Foldable)-instance (NFData s) => NFData (Hourglasses s)- data HourglassBulb = UpperBulb | LowerBulb oneBulb :: HourglassBulb -> (a->a) -> Hourglass a->Hourglass a oneBulb UpperBulb f (Hourglass u l) = Hourglass (f u) l@@ -256,6 +315,9 @@ -- | Experimental. There might be a more powerful instance possible. instance (AffineManifold x) => Semimanifold (ShadeTree x) where type Needle (ShadeTree x) = Diff x+ fromInterior = id+ toInterior = pure+ translateP = Tagged (.+~^) PlainLeaves xs .+~^ v = PlainLeaves $ (.+^v)<$>xs OverlappingBranches n sh br .+~^ v = OverlappingBranches n (sh.+~^v)@@ -279,25 +341,23 @@ where ne (PlainLeaves []) = False; ne _ = True --- | Build a really quite nicely balanced tree from a cloud of points, on--- any real manifold.+-- | Build a quite nicely balanced tree from a cloud of points, on any real manifold. -- -- Example: -- -- @--- > :m +Graphics.Dynamic.Plot.R2 Data.Manifold.TreeCover Data.VectorSpace Data.AffineSpace--- > import Diagrams.Prelude ((^&), P2, R2, circle, fc, (&), moveTo, green)--- --- > let testPts0 = [0^&0, 0^&1, 1^&1, 1^&2, 2^&2] :: [P2] -- Generate sort-of–random point cloud--- > let testPts1 = [p .+^ v^/3 | p<-testPts0, v <- [0^&0, (-1)^&1, 1^&2]]--- > let testPts2 = [p .+^ v^/4 | p<-testPts1, v <- [0^&0, (-1)^&1, 1^&2]]--- > let testPts3 = [p .+^ v^/5 | p<-testPts2, v <- [0^&0, (-2)^&1, 1^&2]]--- > let testPts4 = [p .+^ v^/7 | p<-testPts3, v <- [0^&1, (-2)^&1, 1^&2]]--- > length testPts4--- 405+-- > :m +Graphics.Dynamic.Plot.R2 Data.Manifold.TreeCover Data.VectorSpace Data.AffineSpace +-- > import Diagrams.Prelude ((^&), p2, r2, P2, circle, fc, (&), moveTo, opacity) -- --- > plotWindow [ plot . onlyNodes $ fromLeafPoints testPts4--- > , plot [circle 0.06 & moveTo p & fc green :: PlainGraphics | p <- testPts4] ]+-- > -- Generate sort-of–random cloud of lots of points+-- > let testPts0 = p2 \<$\> [(0,0), (0,1), (1,1), (1,2), (2,2)] :: [P2 Double]+-- > let testPts1 = [p .+^ v^/3 | p\<-testPts0, v \<- r2\<$\>[(0,0), (-1,1), (1,2)]]+-- > let testPts2 = [p .+^ v^/4 | p\<-testPts1, v \<- r2\<$\>[(0,0), (-1,1), (1,2)]]+-- > let testPts3 = [p .+^ v^/5 | p\<-testPts2, v \<- r2\<$\>[(0,0), (-2,1), (1,2)]]+-- > let testPts4 = [p .+^ v^/7 | p\<-testPts3, v \<- r2\<$\>[(0,1), (-1,1), (1,2)]]+-- +-- > plotWindow [ plot [ shapePlot $ circle 0.06 & moveTo p & opacity 0.3 | p <- testPts4 ]+-- > , plot . onlyNodes $ 'fromLeafPoints' testPts4 ] -- @ -- -- <<images/examples/simple-2d-ShadeTree.png>>@@ -311,7 +371,7 @@ Just redBrchs -> OverlappingBranches (length xs) rShade- (branchProc (shadeExpanse rShade) redBrchs)+ (branchProc (_shadeExpanse rShade) redBrchs) _ -> PlainLeaves xs partitions -> DisjointBranches (length xs) . NE.fromList@@ -406,6 +466,9 @@ BaryCoords v .+^ w = BaryCoords $ v ^+^ w instance (KnownNat n) => Semimanifold (BaryCoords n) where type Needle (BaryCoords n) = FreeVect n ℝ+ fromInterior = id+ toInterior = pure+ translateP = Tagged (.+~^) (.+~^) = (.+^) instance (KnownNat n) => PseudoAffine (BaryCoords n) where (.-~.) = pure .: (.-.)@@ -874,6 +937,109 @@ sShSaw _ _ = error "`sShSaw` is not supposed to cut anything else but `OverlappingBranches`" ++-- | Essentially the same as @(x,y)@, but not considered as a product topology.+-- The 'Semimanifold' etc. instances just copy the topology of @x@, ignoring @y@.+data x`WithAny`y+ = WithAny { _untopological :: y+ , _topological :: !x }+ deriving (Hask.Functor)++instance (Semimanifold x) => Semimanifold (x`WithAny`y) where+ type Needle (WithAny x y) = Needle x+ type Interior (WithAny x y) = Interior x `WithAny` y+ WithAny y x .+~^ δx = WithAny y $ x.+~^δx+ fromInterior (WithAny y x) = WithAny y $ fromInterior x+ toInterior (WithAny y x) = fmap (WithAny y) $ toInterior x+ translateP = tpWD+ where tpWD :: ∀ x y . Semimanifold x => Tagged (WithAny x y)+ (Interior x`WithAny`y -> Needle x -> Interior x`WithAny`y)+ tpWD = Tagged `id` \(WithAny y x) δx -> WithAny y $ tpx x δx+ where Tagged tpx = translateP :: Tagged x (Interior x -> Needle x -> Interior x)+ +instance (PseudoAffine x) => PseudoAffine (x`WithAny`y) where+ WithAny _ x .-~. WithAny _ ξ = x.-~.ξ++instance (AffineSpace x) => AffineSpace (x`WithAny`y) where+ type Diff (WithAny x y) = Diff x+ WithAny _ x .-. WithAny _ ξ = x.-.ξ+ WithAny y x .+^ δx = WithAny y $ x.+^δx ++instance (VectorSpace x, Monoid y) => VectorSpace (x`WithAny`y) where+ type Scalar (WithAny x y) = Scalar x+ μ *^ WithAny y x = WithAny y $ μ*^x ++instance (AdditiveGroup x, Monoid y) => AdditiveGroup (x`WithAny`y) where+ zeroV = WithAny mempty zeroV+ negateV (WithAny y x) = WithAny y $ negateV x+ WithAny y x ^+^ WithAny υ ξ = WithAny (mappend y υ) (x^+^ξ)++instance (AdditiveGroup x) => Hask.Applicative (WithAny x) where+ pure x = WithAny x zeroV+ WithAny f x <*> WithAny t ξ = WithAny (f t) (x^+^ξ)+ +instance (AdditiveGroup x) => Hask.Monad (WithAny x) where+ return x = WithAny x zeroV+ WithAny y x >>= f = WithAny r $ x^+^q+ where WithAny r q = f y++shadeWithoutAnything :: Shade (x`WithAny`y) -> Shade x+shadeWithoutAnything (Shade (WithAny _ b) e) = Shade b e++-- | This is to 'ShadeTree' as 'Data.Map.Map' is to 'Data.Set.Set'.+type x`Shaded`y = ShadeTree (x`WithAny`y)++stiWithDensity :: (WithField ℝ Manifold x, WithField ℝ LinearManifold y)+ => x`Shaded`y -> x -> Cℝay y+stiWithDensity (PlainLeaves lvs)+ | [locShape@(Shade baryc expa)] <- pointsShades $ _topological <$> lvs+ = let nlvs = fromIntegral $ length lvs :: ℝ+ indiShapes = [(Shade p expa, y) | WithAny y p <- lvs]+ in \x -> let lcCoeffs = [ occlusion psh x | (psh, _) <- indiShapes ]+ dens = sum lcCoeffs+ in mkCone dens . linearCombo . zip (snd<$>indiShapes)+ $ (/dens)<$>lcCoeffs+stiWithDensity (DisjointBranches _ lvs)+ = \x -> foldr1 qGather $ (`stiWithDensity`x)<$>lvs+ where qGather (Cℝay 0 _) o = o+ qGather o _ = o+stiWithDensity (OverlappingBranches n (Shade (WithAny _ bc) extend) brs) = ovbSWD+ where ovbSWD x = case x .-~. bc of+ Option (Just v)+ | dist² <- metricSq ε v+ , dist² < 9+ , att <- exp(1/(dist²-9)+1/9)+ -> qGather att $ fmap ($x) downPrepared+ _ -> coneTip+ ε = recipMetric extend+ downPrepared = dp =<< brs+ where dp (DBranch _ (Hourglass up dn))+ = fmap stiWithDensity $ up:|[dn]+ qGather att contribs = mkCone (att*dens)+ $ linearCombo [(v, d/dens) | Cℝay d v <- NE.toList contribs]+ where dens = sum (hParamCℝay <$> contribs)+++smoothInterpolate :: (WithField ℝ Manifold x, WithField ℝ LinearManifold y)+ => NonEmpty (x,y) -> x -> y+smoothInterpolate l = \x ->+ case ltr x of+ Cℝay 0 _ -> defy+ Cℝay _ y -> y+ where defy = linearCombo [(y, 1/n) | WithAny y _ <- l']+ n = fromIntegral $ length l'+ l' = (uncurry WithAny . swap) <$> NE.toList l+ ltr = stiWithDensity $ fromLeafPoints l'+++coneTip :: (AdditiveGroup v) => Cℝay v+coneTip = Cℝay 0 zeroV++mkCone :: AdditiveGroup v => ℝ -> v -> Cℝay v+mkCone 0 _ = coneTip+mkCone h v = Cℝay h v++ foci :: [a] -> [(a,[a])] foci [] = [] foci (x:xs) = (x,xs) : fmap (second (x:)) (foci xs)@@ -909,4 +1075,14 @@ ] superFlatView = foldMap go . flatView where go (t,ds) = t : ds+++++++extractJust :: (a->Maybe b) -> [a] -> (Maybe b, [a])+extractJust f [] = (Nothing,[])+extractJust f (x:xs) | Just r <- f x = (Just r, xs)+ | otherwise = second (x:) $ extractJust f xs
Data/Manifold/Types.hs view
@@ -72,6 +72,7 @@ import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine+import Data.Manifold.Cone import Data.LinearMap.HerMetric import Data.VectorSpace.FiniteDimensional @@ -85,20 +86,14 @@ #define deriveAffine(c,t) \ instance (c) => Semimanifold (t) where { \ type Needle (t) = Diff (t); \- (.+~^) = (.+^) }; \-instance (c) => PseudoAffine (t) where { \+ fromInterior = id; \+ toInterior = pure; \+ translateP = Tagged (.+~^); \+ (.+~^) = (.+^) }; \+instance (c) => PseudoAffine (t) where { \ a.-~.b = pure (a.-.b); } --- | The /n/-th Stiefel manifold is the space of all possible configurations of--- /n/ orthonormal vectors. In the case /n/ = 1, simply the subspace of normalised--- vectors, i.e. equivalent to the 'UnitSphere'. Even so, it strictly speaking--- requires the containing space to be at least metric (if not Hilbert); we would--- however like to be able to use this concept also in spaces with no inner product,--- therefore we define this space not as normalised vectors, but rather as all--- vectors modulo scaling by positive factors.-newtype Stiefel1 v = Stiefel1 { getStiefel1N :: DualSpace v }- newtype Stiefel1Needle v = Stiefel1Needle { getStiefel1Tangent :: HMat.Vector (Scalar v) } newtype Stiefel1Basis v = Stiefel1Basis { getStiefel1Basis :: Int } s1bTrie :: forall v b. FiniteDimensional v => (Stiefel1Basis v->b) -> Stiefel1Basis v:->:b@@ -163,19 +158,12 @@ instance (WithField k LinearManifold v, Real k) => Semimanifold (Stiefel1 v) where type Needle (Stiefel1 v) = Stiefel1Needle v+ fromInterior = id+ toInterior = pure+ translateP = Tagged (.+~^) Stiefel1 s .+~^ Stiefel1Needle n = Stiefel1 . fromPackedVector . HMat.scale (signum s'i) $ if| ν==0 -> s' -- ν'≡0 is a special case of this, so we can otherwise assume ν'>0.--- -- | ν<=1 -> let -- κ = (-1 − 1/(ν−1)) / ν'--- -- m ∝ spro + κ · n--- -- ∝ (1−ν) · spro + (1−ν) · κ · n--- -- = (1−ν) · spro + (-(1−ν) − -1)/ν' · n--- m = HMat.scale (1-ν) spro + HMat.scale (ν/ν') n--- in insi (1-ν) m- | ν<=2 -> let -- κ = (1/(ν−1) − 1) / ν'- -- m ∝ - spro + κ · n- -- ∝ (1−ν) · spro + (ν−1) · κ · n- -- = (1−ν) · spro + (1 − (ν−1))/ν' · n- m = HMat.scale ιmν spro + HMat.scale ((1-abs ιmν)/ν') n+ | ν<=2 -> let m = HMat.scale ιmν spro + HMat.scale ((1-abs ιmν)/ν') n ιmν = 1-ν in insi ιmν m | otherwise -> let m = HMat.scale ιmν spro + HMat.scale ((abs ιmν-1)/ν') n@@ -199,8 +187,7 @@ s'i | v <- HMat.scale (recip s'i) delis - tpro , absv <- l2norm v , absv > 0- -> let μ -- = (1 − recip (|v| + 1)) / |v| for sgn sᵢ = sgn tᵢ- = (signum (t'i/s'i) - recip(absv + 1)) / absv+ -> let μ = (signum (t'i/s'i) - recip(absv + 1)) / absv in HMat.scale μ v | t'i/s'i > 0 -> samePoint | otherwise -> antipode@@ -214,34 +201,21 @@ samePoint = (d-1) HMat.|> repeat 0 antipode = (d-1) HMat.|> (2 : repeat 0) -l2norm :: MetricScalar s => HMat.Vector s -> s-l2norm = realToFrac . HMat.norm_2 +instance ( WithField ℝ HilbertSpace x ) => ConeSemimfd (Stiefel1 x) where+ type CℝayInterior (Stiefel1 x) = x+ fromCℝayInterior (FinVecArrRep v) = case HMat.size v of+ 0 -> Cℝay 0 $ Stiefel1 zeroV+ _ -> Cℝay (HMat.norm_2 v) $ Stiefel1 (fromPackedVector v)+ toCℝayInterior (Cℝay 0 _) = pure zeroV+ toCℝayInterior (Cℝay l (Stiefel1 v))+ = pure.FinVecArrRep $ HMat.scale (l/HMat.norm_2 v') v'+ where v' = asPackedVector v -stiefel1Project :: LinearManifold v =>- DualSpace v -- ^ Must be nonzero.- -> Stiefel1 v-stiefel1Project = Stiefel1 -stiefel1Embed :: HilbertSpace v => Stiefel1 v -> v-stiefel1Embed (Stiefel1 n) = normalized n- --class (PseudoAffine v, InnerSpace v, NaturallyEmbedded (UnitSphere v) (DualSpace v))- => HasUnitSphere v where- type UnitSphere v :: *- stiefel :: UnitSphere v -> Stiefel1 v- stiefel = Stiefel1 . embed- unstiefel :: Stiefel1 v -> UnitSphere v- unstiefel = coEmbed . getStiefel1N--instance HasUnitSphere ℝ where type UnitSphere ℝ = S⁰-instance HasUnitSphere ℝ² where type UnitSphere ℝ² = S¹-instance HasUnitSphere ℝ³ where type UnitSphere ℝ³ = S²+l2norm :: MetricScalar s => HMat.Vector s -> s+l2norm = realToFrac . HMat.norm_2 -instance (HasUnitSphere v, v ~ DualSpace v) => NaturallyEmbedded (Stiefel1 v) v where- embed = embed . unstiefel- coEmbed = stiefel . coEmbed
Data/Manifold/Types/Primitive.hs view
@@ -37,7 +37,7 @@ , Projective1, Projective2 , Disk1, Disk2, Cone, OpenCone -- * Linear manifolds- , ZeroDim(..)+ , ZeroDim(..), isoAttachZeroDim , ℝ⁰, ℝ, ℝ², ℝ³ -- * Hyperspheres , S⁰(..), S¹(..), S²(..)@@ -49,7 +49,7 @@ , CD¹(..), Cℝay(..) -- * Utility (deprecated) , NaturallyEmbedded(..)- , GraphWindowSpec(..), Endomorphism, (^), EqFloating+ , GraphWindowSpec(..), Endomorphism, (^), (^.), EqFloating ) where @@ -60,6 +60,8 @@ import Data.Void import Data.Monoid +import Control.Applicative (Const(..))+ import qualified Prelude import Control.Category.Constrained.Prelude hiding ((^))@@ -67,11 +69,13 @@ import Control.Monad.Constrained import Data.Foldable.Constrained +import Data.Embedding + type EqFloating f = (Eq f, Ord f, Floating f) @@ -101,6 +105,13 @@ decompose Origin = [] decompose' Origin = absurd +{-# INLINE isoAttachZeroDim #-}+isoAttachZeroDim :: ( WellPointed c, UnitObject c ~ (), ObjectPair c a ()+ , Object c (ZeroDim k), ObjectPair c a (ZeroDim k)+ , PointObject c (ZeroDim k) )+ => Isomorphism c a (a, ZeroDim k)+isoAttachZeroDim = second (Isomorphism (const Origin) terminal) . attachUnit+ -- | The zero-dimensional sphere is actually just two points. Implementation might -- therefore change to @ℝ⁰ 'Control.Category.Constrained.+' ℝ⁰@: the disjoint sum of two -- single-point spaces.@@ -247,7 +258,14 @@ () <.> () = 0 +infixr 8 ^ (^) :: Num a => a -> Int -> a (^) = (Prelude.^)+++infixl 8 ^.+{-# INLINE (^.) #-}+(^.) :: s -> (forall f . Prelude.Functor f => (a->f a) -> s->f s) -> a+o ^. g = getConst (g Const o)
Data/VectorSpace/FiniteDimensional.hs view
@@ -10,13 +10,16 @@ {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE TupleSections #-} {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE PolyKinds #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UnicodeSyntax #-} @@ -24,6 +27,7 @@ module Data.VectorSpace.FiniteDimensional ( FiniteDimensional(..) , SmoothScalar + , FinVecArrRep(..), concreteArrRep, (⊗), splitArrRep ) where @@ -31,6 +35,7 @@ import Prelude hiding ((^)) +import Data.AffineSpace import Data.VectorSpace import Data.LinearMap import Data.Basis@@ -42,7 +47,10 @@ import Data.Manifold.Types.Primitive import Data.CoNat+import Data.Embedding +import Control.Arrow+ import qualified Data.Vector as Arr import qualified Numeric.LinearAlgebra.HMatrix as HMat @@ -55,7 +63,6 @@ , Num(HMat.Vector s), HMat.Indexable(HMat.Vector s)s , HMat.Normed(HMat.Vector s) ) - -- | Many linear algebra operations are best implemented via packed, dense 'HMat.Matrix'es. -- For one thing, that makes common general vector operations quite efficient, -- in particular on high-dimensional spaces.@@ -160,4 +167,100 @@ fromPackedVector arr = FreeVect (Arr.convert arr) -- asPackedMatrix = _ -- could be done quite efficiently here! +++-- | Semantically the same as @'Tagged' tag refvs@, but directly uses the+-- packed-vector array representation.+-- +-- The tag should really be kind-polymorphic, but at least GHC-7.8 doesn't quite+-- handle the associated types of the manifold classes then.+newtype FinVecArrRep (tag :: * -> *) refvs scalar+ = FinVecArrRep { getFinVecArrRep :: HMat.Vector scalar }++instance (SmoothScalar s) => AffineSpace (FinVecArrRep t b s) where+ type Diff (FinVecArrRep t b s) = FinVecArrRep t b s+ (.-.) = (^-^)+ (.+^) = (^+^)+ +instance (SmoothScalar s) => AdditiveGroup (FinVecArrRep t b s) where+ zeroV = FinVecArrRep $ 0 HMat.|> []+ negateV (FinVecArrRep v) = FinVecArrRep $ negate v+ FinVecArrRep v ^+^ FinVecArrRep w+ | HMat.size v == 0 = FinVecArrRep w+ | HMat.size w == 0 = FinVecArrRep w+ | otherwise = FinVecArrRep $ v + w++instance (SmoothScalar s) => VectorSpace (FinVecArrRep t b s) where+ type Scalar (FinVecArrRep t b s) = s+ μ *^ FinVecArrRep v = FinVecArrRep $ HMat.scale μ v++instance (SmoothScalar s) => InnerSpace (FinVecArrRep t b s) where+ FinVecArrRep v <.> FinVecArrRep w+ | HMat.size v == 0 = 0+ | HMat.size w == 0 = 0+ | otherwise = v`HMat.dot`w++concreteArrRep :: (SmoothScalar s, FiniteDimensional r, Scalar r ~ s)+ => Isomorphism (->) r (FinVecArrRep t r s)+concreteArrRep = Isomorphism (FinVecArrRep . asPackedVector)+ (fromPackedVector . getFinVecArrRep)++(⊗) :: ∀ t s v w . ( SmoothScalar s, FiniteDimensional v, FiniteDimensional w+ , Scalar v ~ s, Scalar w ~ s )+ => FinVecArrRep t v s -> FinVecArrRep t w s -> FinVecArrRep t (v,w) s+FinVecArrRep v ⊗ FinVecArrRep w+ | HMat.size v + HMat.size w == 0 = FinVecArrRep v+ | HMat.size v == 0 = FinVecArrRep $ HMat.vjoin [HMat.konst 0 nv, w]+ | HMat.size w == 0 = FinVecArrRep $ HMat.vjoin [v, HMat.konst 0 nw]+ | otherwise = FinVecArrRep $ HMat.vjoin [v,w]+ where Tagged nv = dimension :: Tagged v Int+ Tagged nw = dimension :: Tagged w Int++splitArrRep :: ∀ t s v w . ( SmoothScalar s, FiniteDimensional v, FiniteDimensional w+ , Scalar v ~ s, Scalar w ~ s )+ => FinVecArrRep t (v,w) s -> (FinVecArrRep t v s, FinVecArrRep t w s)+splitArrRep (FinVecArrRep vw)+ | HMat.size vw == 0 = (FinVecArrRep vw, FinVecArrRep vw)+ | otherwise = ( FinVecArrRep $ HMat.subVector 0 nv vw+ , FinVecArrRep $ HMat.subVector nv nw vw )+ where Tagged nv = dimension :: Tagged v Int+ Tagged nw = dimension :: Tagged w Int+ ++instance (SmoothScalar s, FiniteDimensional r, Scalar r ~ s)+ => HasBasis (FinVecArrRep t r s) where+ type Basis (FinVecArrRep t r s) = Basis r+ basisValue = (concreteArrRep$->$) . basisValue+ decompose = decompose . (concreteArrRep$<-$)+ decompose' = decompose' . (concreteArrRep$<-$)++instance (SmoothScalar s, FiniteDimensional r, Scalar r ~ s)+ => FiniteDimensional (FinVecArrRep t r s) where+ dimension = d+ where d :: ∀ t r s . FiniteDimensional r => Tagged (FinVecArrRep t r s) Int+ d = Tagged n+ where Tagged n = dimension :: Tagged r Int+ indexBasis = d+ where d :: ∀ t r s . FiniteDimensional r => Tagged (FinVecArrRep t r s) (Int -> Basis r)+ d = Tagged n+ where Tagged n = indexBasis :: Tagged r (Int -> Basis r)+ basisIndex = d+ where d :: ∀ t r s . FiniteDimensional r => Tagged (FinVecArrRep t r s) (Basis r -> Int)+ d = Tagged n+ where Tagged n = basisIndex :: Tagged r (Basis r -> Int)+ asPackedVector = apv+ where apv :: ∀ t r s . (FiniteDimensional r, SmoothScalar s)+ => FinVecArrRep t r s -> HMat.Vector s+ apv (FinVecArrRep v)+ | HMat.size v == 0 = HMat.konst 0 n+ | otherwise = v+ where Tagged n = dimension :: Tagged r Int+ fromPackedVector = FinVecArrRep+++instance (NaturallyEmbedded m r, FiniteDimensional r, s ~ Scalar r)+ => NaturallyEmbedded m (FinVecArrRep t r s) where+ embed = (concreteArrRep$<-$) . embed+ coEmbed = coEmbed . (concreteArrRep$->$)+
images/examples/simple-2d-ShadeTree.png view
binary file changed (107507 → 103235 bytes)
manifolds.cabal view
@@ -1,8 +1,8 @@ Name: manifolds-Version: 0.1.3.1+Version: 0.1.5.0 Category: Math Synopsis: Working with manifolds in a direct, embedding-free way.-Description: Manifolds, a generalisation of the notion of \"smooth curves\" or sufaces,+Description: Manifolds, a generalisation of the notion of “smooth curves” or surfaces, are topological spaces /locally homeomorphic to a vector space/. This gives rise to what is actually the most natural / mathematically elegant way of dealing with them: calculations can be carried out locally, in connection with Riemannian@@ -41,7 +41,6 @@ , vector-space>=0.8 , MemoTrie , vector- , vector-algorithms , hmatrix >= 0.16 && < 0.18 , containers , comonad@@ -68,8 +67,11 @@ Data.LinearMap.HerMetric -- Data.Manifold.Visualisation.R3.GLUT Data.Manifold.Types+ Data.Manifold.Griddable+ Data.Manifold.Riemannian Other-modules: Data.List.FastNub Data.Manifold.Types.Primitive+ Data.Manifold.Cone Data.CoNat Data.Embedding Data.LinearMap.Category