manifolds 0.1.0.2 → 0.1.3.0
raw patch · 14 files changed
+3265/−206 lines, 14 filesdep +deepseqdep +hmatrixdep −MonadRandomdep −randombinary-added
Dependencies added: deepseq, hmatrix
Dependencies removed: MonadRandom, random
Files
- Data/CoNat.hs +314/−0
- Data/Embedding.hs +167/−0
- Data/LinearMap/Category.hs +190/−0
- Data/LinearMap/HerMetric.hs +263/−53
- Data/List/FastNub.hs +29/−7
- Data/Manifold.hs +3/−2
- Data/Manifold/PseudoAffine.hs +241/−70
- Data/Manifold/TreeCover.hs +910/−0
- Data/Manifold/Types.hs +214/−69
- Data/Manifold/Types/Primitive.hs +253/−0
- Data/SimplicialComplex.hs +500/−0
- Data/VectorSpace/FiniteDimensional.hs +163/−0
- images/examples/simple-2d-ShadeTree.png binary
- manifolds.cabal +18/−5
+ Data/CoNat.hs view
@@ -0,0 +1,314 @@+-- |+-- Module : Data.CoNat+-- 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 MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE UnicodeSyntax #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE PolyKinds #-}++module Data.CoNat where++import Data.Tagged+import Data.Semigroup++import Data.MemoTrie+import Data.VectorSpace+import Data.AffineSpace+import Data.Basis+import Data.AdditiveGroup+import qualified Data.List as List+ +import qualified Prelude as Hask hiding(foldl)+import qualified Control.Applicative as Hask+import qualified Control.Monad as Hask+import qualified Data.Foldable as Hask+import qualified Data.Traversable as Hask+++import Control.Category.Constrained.Prelude hiding ((^))+import Data.Traversable.Constrained+++import qualified Data.Vector as Arr+import qualified Numeric.LinearAlgebra.HMatrix as HMat++import Unsafe.Coerce++ +-- | Mainly intended to be used as a data kind.+-- Of course, we'd rather use "GHC.TypeLits" naturals, but they aren't mature enough yet.+data Nat = Z | S Nat deriving (Eq)++natToInt :: Nat -> Int+natToInt Z = 0; natToInt (S n) = 1 + natToInt n++fromNat :: Num a => Nat -> a+fromNat = fromIntegral . natToInt++natTagLast :: forall n f n' . (KnownNat n, Num n') => Tagged (f n) n'+natTagLast = retag (theNatN :: Tagged n n')+natTagPænultimate :: forall n f n' x . (KnownNat n, Num n') => Tagged (f n x) n'+natTagPænultimate = retag (theNatN :: Tagged n n')+natTagAntepænultimate :: forall n f n' x y . (KnownNat n, Num n') => Tagged (f n x y) n'+natTagAntepænultimate = retag (theNatN :: Tagged n n')++natSelfSucc :: forall n . KnownNat n => Tagged (S n) Nat+natSelfSucc = Tagged $ S n+ where (Tagged n) = theNat :: Tagged n Nat+natSelfSuccN :: forall n a . (KnownNat n, Num a) => Tagged (S n) a+natSelfSuccN = Tagged $ n + 1+ where (Tagged n) = theNatN :: Tagged n a++class KnownNat (n :: Nat) where+ theNat :: Tagged n Nat+ theNatN :: Num n' => Tagged n n'+ + cozero :: s Z -> Option (s n)+ cozeroT :: c Z x -> Option (c n x)+ + cosucc :: (forall k . KnownNat k => s (S k)) -> Option (s n)+ fCosucc :: Hask.Alternative f => (forall k . KnownNat k => f (s (S k))) -> f (s n)+ cosuccT :: (forall k . KnownNat k => s (S k) x) -> Option (s n x)+ fCosuccT :: Hask.Alternative f => (forall k . KnownNat k => f (s (S k) x)) -> f (s n x)+ + coNat :: (s Z->r) -> ( forall k . KnownNat k => s (S k) -> r ) -> s n -> r+ coNatT :: (c Z x->r) -> ( forall k . KnownNat k => c (S k) x -> r ) -> c n x -> r+ + coInduce :: s Z -> (forall k . KnownNat k => s k -> s (S k)) -> s n+ coInduceT :: c Z x -> (forall k . KnownNat k => c k x -> c (S k) x) -> c n x+ + ftorCoInduce :: f (s Z) -> (forall k . KnownNat k => f (s k) -> f (s (S k))) -> f (s n)+ ftorCoInduceT :: f (c Z x) -> (forall k . KnownNat k => f (c k x) -> f (c (S k) x))+ -> f (c n x)++ tryToMatch :: KnownNat k => (∀ j . KnownNat j => b j -> b (S j)) -> b k -> Option (b n)+++instance KnownNat Z where+ theNat = Tagged Z+ theNatN = Tagged 0+ cozero = pure; cosucc _ = Hask.empty; fCosucc _ = Hask.empty+ cozeroT = pure; cosuccT _ = Hask.empty; fCosuccT _ = Hask.empty+ coNat f _ = f; coNatT f _ = f+ coInduce s _ = s+ coInduceT s _ = s+ ftorCoInduce s _ = s+ ftorCoInduceT s _ = s+ tryToMatch = ttmZ+ where ttmZ :: ∀ b k . KnownNat k+ => (∀ j . KnownNat j => b j -> b (S j)) -> b k -> Option (b Z)+ ttmZ sc nf = case k of+ Z -> return $ unsafeCoerce nf+ S _ -> Hask.empty+ where (Tagged k) = theNat :: Tagged k Nat+instance (KnownNat n) => KnownNat (S n) where+ theNat = natSelfSucc+ theNatN = natSelfSuccN+ cozero _ = Hask.empty; cosucc v = pure v; fCosucc v = v+ cozeroT _ = Hask.empty; cosuccT v = pure v; fCosuccT v = v+ coNat _ f = f; coNatT _ f = f+ coInduce s f = f $ coInduce s f+ coInduceT s f = f $ coInduceT s f+ ftorCoInduce s f = f $ ftorCoInduce s f+ ftorCoInduceT s f = f $ ftorCoInduceT s f+ tryToMatch = ttmS+ where ttmS :: ∀ b k n . (KnownNat k, KnownNat n)+ => (∀ j . KnownNat j => b j -> b (S j)) -> b k -> Option (b (S n))+ ttmS sc nf | k == sn = return $ unsafeCoerce nf+ | k <= sn = tryToMatch sc $ sc nf+ | otherwise = Hask.empty+ where (Tagged k) = theNatN :: Tagged k Int+ (Tagged sn) = theNatN :: Tagged (S n) Int+ +++newtype NatTagAtPænultimate t x n+ = NatTagAtPænultimate { getNatTagAtPænultimate :: t n x }+mapNatTagAtPænultimate :: (s n x -> t m y)+ -> NatTagAtPænultimate s x n -> NatTagAtPænultimate t y m+mapNatTagAtPænultimate f (NatTagAtPænultimate x) = NatTagAtPænultimate $ f x++newtype NatTagAtAntepænultimate t x y n+ = NatTagAtAntepænultimate { getNatTagAtAntepænultimate :: t n x y }+mapNatTagAtAntepænultimate :: (s n w x -> t m y z)+ -> NatTagAtAntepænultimate s w x n -> NatTagAtAntepænultimate t y z m+mapNatTagAtAntepænultimate f (NatTagAtAntepænultimate x) = NatTagAtAntepænultimate $ f x++newtype NatTagAtPreantepænultimate t x y z n+ = NatTagAtPreantepænultimate { getNatTagAtPreantepænultimate :: t n x y z }+mapNatTagAtPreantepænultimate :: (s n u v w -> t m x y z)+ -> NatTagAtPreantepænultimate s u v w n -> NatTagAtPreantepænultimate t x y z m+mapNatTagAtPreantepænultimate f (NatTagAtPreantepænultimate x) = NatTagAtPreantepænultimate $ f x++newtype NatTagAtFtorUltimate f t n+ = NatTagAtFtorUltimate { getNatTagAtFtorUltimate :: f (t n) }+mapNatTagAtFtorUltimate :: (f (s n) -> f (t m))+ -> NatTagAtFtorUltimate f s n -> NatTagAtFtorUltimate f t m+mapNatTagAtFtorUltimate f (NatTagAtFtorUltimate x) = NatTagAtFtorUltimate $ f x++newtype NatTagAtFtorPænultimate f t x n+ = NatTagAtFtorPænultimate { getNatTagAtFtorPænultimate :: f (t n x) }+mapNatTagAtFtorPænultimate :: (f (s n x) -> f (t m y))+ -> NatTagAtFtorPænultimate f s x n -> NatTagAtFtorPænultimate f t y m+mapNatTagAtFtorPænultimate f (NatTagAtFtorPænultimate x) = NatTagAtFtorPænultimate $ f x++newtype NatTagAtFtorAntepænultimate f t x y n+ = NatTagAtFtorAntepænultimate { getNatTagAtFtorAntepænultimate :: f (t n x y) }+mapNatTagAtFtorAntepænultimate :: (f (s n w x) -> f (t m y z))+ -> NatTagAtFtorAntepænultimate f s w x n -> NatTagAtFtorAntepænultimate f t y z m+mapNatTagAtFtorAntepænultimate f (NatTagAtFtorAntepænultimate x) = NatTagAtFtorAntepænultimate $ f x+++tryToMatchT :: (KnownNat k, KnownNat j)+ => (∀ n . KnownNat n => c n x -> c (S n) x) -> c k x -> Option (c j x)+tryToMatchT f = fmap getNatTagAtPænultimate+ . tryToMatch (mapNatTagAtPænultimate f) . NatTagAtPænultimate+tryToMatchTT ::(KnownNat k, KnownNat j) => (∀ n . KnownNat n => d n x y -> d (S n) x y) -> d k x y -> Option (d j x y)+tryToMatchTT f = fmap getNatTagAtAntepænultimate+ . tryToMatch (mapNatTagAtAntepænultimate f) . NatTagAtAntepænultimate+tryToMatchTTT :: (KnownNat k, KnownNat j) => (∀ n . KnownNat n => e n x y z -> e (S n) x y z)+ -> e k x y z -> Option (e j x y z)+tryToMatchTTT f = fmap getNatTagAtPreantepænultimate+ . tryToMatch (mapNatTagAtPreantepænultimate f) . NatTagAtPreantepænultimate++ftorTryToMatch :: (KnownNat k, KnownNat j) =>+ (∀ n . KnownNat n => f (b n) -> f (b (S n))) -> f (b k) -> Option (f (b j))+ftorTryToMatch f = fmap getNatTagAtFtorUltimate+ . tryToMatch (mapNatTagAtFtorUltimate f) . NatTagAtFtorUltimate+ftorTryToMatchT :: (KnownNat k, KnownNat j) => (∀ n . KnownNat n => f (c n x) -> f (c (S n) x)) -> f (c k x) -> Option (f (c j x))+ftorTryToMatchT f = fmap getNatTagAtFtorPænultimate+ . tryToMatch (mapNatTagAtFtorPænultimate f) . NatTagAtFtorPænultimate+ftorTryToMatchTT :: (KnownNat k, KnownNat j) => (∀ n . KnownNat n => f (d n x y) -> f (d (S n) x y)) -> f (d k x y) -> Option (f (d j x y))+ftorTryToMatchTT f = fmap getNatTagAtFtorAntepænultimate+ . tryToMatch (mapNatTagAtFtorAntepænultimate f) . NatTagAtFtorAntepænultimate+++++++newtype Range (n::Nat) = InRange { getInRange :: Int -- ^ MUST be between 0 and @n-1@.+ }++clipToRange :: forall n . KnownNat n => Int -> Option (Range n)+clipToRange = \k -> if k < n then Hask.pure $ InRange n+ else Hask.empty+ where (Tagged n) = theNatN :: Tagged n Int+ +instance KnownNat n => HasTrie (Range n) where+ data Range n :->: x = RangeTrie (FreeVect n x)+ trie = RangeTrie . \f -> fmap f ids+ where ids = fmap InRange indices+ untrie (RangeTrie (FreeVect arr)) = \(InRange i) -> arr Arr.! i+ enumerate (RangeTrie (FreeVect arr)) = Arr.ifoldr (\i x l -> (InRange i, x) : l) [] arr+++newtype FreeVect (n::Nat) x = FreeVect+ { getFreeVect :: Arr.Vector x -- ^ MUST have length @n@.+ } deriving (Hask.Functor, Hask.Foldable)++instance (KnownNat n) => Hask.Applicative (FreeVect n) where+ pure = replicVector+ (<*>) = perfectZipWith ($)+instance (KnownNat n) => Traversable (FreeVect n) (FreeVect n) (->) (->) where+ traverse f (FreeVect v) = fmap FreeVect . runAsHaskFunctor+ $ Hask.traverse (AsHaskFunctor . f) v+instance (KnownNat n, Show x) => Show (FreeVect n x) where+ show (FreeVect v) = "(freeTuple $->$ ("+ ++ List.intercalate "," [show x | x<-Arr.toList v] ++ "))"++type x ^ n = FreeVect n x++instance (Num x, KnownNat n) => AffineSpace (FreeVect n x) where+ type Diff (FreeVect n x) = FreeVect n x+ (.-.) = perfectZipWith (-)+ (.+^) = perfectZipWith (+)+instance (Num x, KnownNat n) => AdditiveGroup (FreeVect n x) where+ zeroV = replicVector 0+ negateV = fmap negate+ (^+^) = perfectZipWith (+)+instance (Num x, KnownNat n) => VectorSpace (FreeVect n x) where+ type Scalar (FreeVect n x) = x+ (*^) = fmap . (*)+instance (Num x, AdditiveGroup x, KnownNat n) => InnerSpace (FreeVect n x) where+ FreeVect v<.>FreeVect w = Arr.sum $ Arr.zipWith (*) v w+instance (Num x, KnownNat n) => HasBasis (FreeVect n x) where+ type Basis (FreeVect n x) = Range n+ basisValue = \(InRange i) -> fmap (\k -> if i==k then 1 else 0) ids+ where ids = indices+ decompose (FreeVect arr) = Arr.ifoldr (\i x l -> (InRange i, x) : l) [] arr+ decompose' (FreeVect arr) (InRange i) = arr Arr.! i+++replicVector :: forall n x . KnownNat n => x -> FreeVect n x+replicVector = FreeVect . Arr.replicate n+ where (Tagged n) = theNatN :: Tagged n Int+++freeVector :: forall l n x . (KnownNat n, Hask.Foldable l) => l x -> Option (FreeVect n x)+freeVector c'+ | length c == n = pure . FreeVect $ Arr.fromList c+ | otherwise = Hask.empty+ where (Tagged n) = theNatN :: Tagged n Int+ c = Hask.toList c'++-- | Free vector containing the (0-based) indices of its fields as entries.+indices :: forall n n' . (KnownNat n, Num n') => FreeVect n n'+indices = FreeVect $ Arr.enumFromN 0 n+ where (Tagged n) = theNatN :: Tagged n Int+++perfectZipWith :: forall n a b c . KnownNat n+ => (a->b->c) -> FreeVect n a -> FreeVect n b -> FreeVect n c+perfectZipWith f (FreeVect va) (FreeVect vb) = FreeVect $ Arr.zipWith f va vb++freeSortBy :: forall n a . KnownNat n+ => (a->a->Ordering) -> a^n -> a^n+freeSortBy cmp (FreeVect xs) = FreeVect $ Arr.fromList (List.sortBy cmp $ Arr.toList xs)++freeRotate :: ∀ n a . KnownNat n => Int -> a^n -> a^n+freeRotate j' = \(FreeVect v) -> FreeVect $ Arr.unsafeBackpermute v rot+ where (Tagged n) = theNatN :: Tagged n Int+ rot = Arr.enumFromN j (n-j) Arr.++ Arr.enumFromN 0 j+ j = j'`mod`n++++freeCons :: a -> FreeVect n a -> FreeVect (S n) a+freeCons x (FreeVect xs) = FreeVect $ Arr.cons x xs++freeSnoc :: FreeVect n a -> a -> FreeVect (S n) a+freeSnoc (FreeVect xs) x = FreeVect $ Arr.snoc xs x+++++newtype AsHaskFunctor f x = AsHaskFunctor { runAsHaskFunctor :: f x }++instance (Functor f (->) (->)) => Hask.Functor (AsHaskFunctor f) where+ fmap f (AsHaskFunctor c) = AsHaskFunctor $ fmap f c+instance (Monoidal f (->) (->)) => Hask.Applicative (AsHaskFunctor f) where+ pure x = fmap (const x) . AsHaskFunctor $ pureUnit ()+ AsHaskFunctor fs <*> AsHaskFunctor xs = AsHaskFunctor . fmap (uncurry ($)) $ fzip (fs, xs)
+ Data/Embedding.hs view
@@ -0,0 +1,167 @@+-- |+-- Module : Data.Embedding+-- 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 MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE UnicodeSyntax #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE DataKinds #-}++module Data.Embedding where++import Data.Tagged+import Data.Semigroup++import qualified Prelude as Hask hiding(foldl)+import qualified Control.Applicative as Hask+import qualified Control.Monad as Hask+import qualified Data.Foldable as Hask+++import Control.Category.Constrained.Prelude hiding ((^))+import Control.Arrow.Constrained+++++data Isomorphism c a b = Isomorphism { forwardIso :: c a b+ , backwardIso :: c b a }++infixr 0 $->$, $<-$+($->$) :: (Function c, Object c a, Object c b) => Isomorphism c a b -> a -> b+Isomorphism f _ $->$ x = f $ x++($<-$) :: (Function c, Object c b, Object c a) => Isomorphism c a b -> b -> a+Isomorphism _ p $<-$ y = p $ y++fromInversePair :: c a b -> c b a -> Isomorphism c a b+fromInversePair = Isomorphism++perfectInvert :: Isomorphism c a b -> Isomorphism c b a+perfectInvert (Isomorphism f b) = Isomorphism b f++instance (Category c) => Category (Isomorphism c) where+ type Object (Isomorphism c) a = Object c a+ id = Isomorphism id id+ Isomorphism e p . Isomorphism f q = Isomorphism (e.f) (q.p)++instance (Cartesian c) => Cartesian (Isomorphism c) where+ type UnitObject (Isomorphism c) = UnitObject c+ type PairObjects (Isomorphism c) a b = PairObjects c a b+ swap = Isomorphism swap swap+ attachUnit = Isomorphism attachUnit detachUnit+ detachUnit = Isomorphism detachUnit attachUnit+ regroup = Isomorphism regroup regroup'+ regroup' = Isomorphism regroup' regroup++instance (CoCartesian c) => CoCartesian (Isomorphism c) where+ type ZeroObject (Isomorphism c) = ZeroObject c+ type SumObjects (Isomorphism c) a b = SumObjects c a b+ coSwap = Isomorphism coSwap coSwap+ attachZero = Isomorphism attachZero detachZero+ detachZero = Isomorphism detachZero attachZero+ coRegroup = Isomorphism coRegroup coRegroup'+ coRegroup' = Isomorphism coRegroup' coRegroup+ maybeAsSum = Isomorphism maybeAsSum maybeFromSum+ maybeFromSum = Isomorphism maybeFromSum maybeAsSum+ boolAsSum = Isomorphism boolAsSum boolFromSum+ boolFromSum = Isomorphism boolFromSum boolAsSum++instance (Morphism c) => Morphism (Isomorphism c) where+ Isomorphism e p *** Isomorphism f q = Isomorphism (e***f) (p***q)+ +instance (MorphChoice c) => MorphChoice (Isomorphism c) where+ Isomorphism e p +++ Isomorphism f q = Isomorphism (e+++f) (p+++q)++instance (Category c) => EnhancedCat c (Isomorphism c) where + arr = forwardIso++instance (Category c) => EnhancedCat (Embedding c) (Isomorphism c) where + arr (Isomorphism f b) = Embedding f b+++ +-- | A pair of matching functions. The projection must be a left (but not necessarily right)+-- inverse of the embedding,+-- i.e. the cardinality of @a@ will have to be less or equal than the cardinality+-- of @b@.+data Embedding c a b = Embedding { embedding :: c a b+ , projection :: c b a+ }++fromEmbedProject :: c a b -- ^ Injective embedding+ -> c b a -- ^ Surjective projection. Must be left inverse of embedding.+ -> Embedding c a b+fromEmbedProject = Embedding+++infixr 0 $->, >-$+($->) :: (Function c, Object c a, Object c b) => Embedding c a b -> a -> b+Embedding f _ $-> x = f $ x++(>-$) :: (Function c, Object c b, Object c a) => Embedding c a b -> b -> a+Embedding _ p >-$ y = p $ y+++instance (Category c) => Category (Embedding c) where+ type Object (Embedding c) a = Object c a+ id = Embedding id id+ Embedding e p . Embedding f q = Embedding (e.f) (q.p)++instance (Cartesian c) => Cartesian (Embedding c) where+ type UnitObject (Embedding c) = UnitObject c+ type PairObjects (Embedding c) a b = PairObjects c a b+ swap = Embedding swap swap+ attachUnit = Embedding attachUnit detachUnit+ detachUnit = Embedding detachUnit attachUnit+ regroup = Embedding regroup regroup'+ regroup' = Embedding regroup' regroup++instance (CoCartesian c) => CoCartesian (Embedding c) where+ type ZeroObject (Embedding c) = ZeroObject c+ type SumObjects (Embedding c) a b = SumObjects c a b+ coSwap = Embedding coSwap coSwap+ attachZero = Embedding attachZero detachZero+ detachZero = Embedding detachZero attachZero+ coRegroup = Embedding coRegroup coRegroup'+ coRegroup' = Embedding coRegroup' coRegroup+ maybeAsSum = Embedding maybeAsSum maybeFromSum+ maybeFromSum = Embedding maybeFromSum maybeAsSum+ boolAsSum = Embedding boolAsSum boolFromSum+ boolFromSum = Embedding boolFromSum boolAsSum++instance (Morphism c) => Morphism (Embedding c) where+ Embedding e p *** Embedding f q = Embedding (e***f) (p***q)+ +instance (MorphChoice c) => MorphChoice (Embedding c) where+ Embedding e p +++ Embedding f q = Embedding (e+++f) (p+++q)++instance (Category c) => EnhancedCat c (Embedding c) where + arr = embedding+++ ++++
+ Data/LinearMap/Category.hs view
@@ -0,0 +1,190 @@+-- |+-- Module : Data.LinearMap.Category+-- 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 MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE UnicodeSyntax #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE DataKinds #-}++module Data.LinearMap.Category where++import Data.Tagged+import Data.Semigroup++import Data.MemoTrie+import Data.VectorSpace+import Data.VectorSpace.FiniteDimensional+import Data.AffineSpace+import Data.Basis+import Data.AdditiveGroup+ +import qualified Prelude as Hask hiding(foldl)+import qualified Control.Applicative as Hask+import qualified Control.Monad as Hask+import qualified Data.Foldable as Hask+++import Control.Category.Constrained.Prelude hiding ((^))+import Control.Arrow.Constrained++import Data.Manifold.Types.Primitive+import Data.CoNat+import Data.Embedding++import qualified Data.Vector as Arr+import qualified Numeric.LinearAlgebra.HMatrix as HMat+++ +-- | A linear mapping between finite-dimensional spaces, implemeted as a dense matrix.+data Linear s a b = DenseLinear { getDenseMatrix :: HMat.Matrix s }++identMat :: forall v w . FiniteDimensional v => Linear (Scalar v) w v+identMat = DenseLinear $ HMat.ident n+ where (Tagged n) = dimension :: Tagged v Int++instance (SmoothScalar s) => Category (Linear s) where+ type Object (Linear s) v = (FiniteDimensional v, Scalar v~s)+ id = identMat+ DenseLinear f . DenseLinear g = DenseLinear $ HMat.mul f g++instance (SmoothScalar s) => Cartesian (Linear s) where+ type UnitObject (Linear s) = ZeroDim s+ swap = lSwap+ where lSwap :: forall v w s+ . (FiniteDimensional v, FiniteDimensional w, Scalar v~s, Scalar w~s)+ => Linear s (v,w) (w,v)+ lSwap = DenseLinear $ HMat.assoc (n,n) 0 l+ where l = [ ((i,i+nv), 1) | i<-[0.. nw-1] ] ++ [ ((i+nw,i), 1) | i<-[0.. nv-1] ] + (Tagged nv) = dimension :: Tagged v Int+ (Tagged nw) = dimension :: Tagged w Int+ n = nv + nw+ attachUnit = identMat+ detachUnit = identMat+ regroup = identMat+ regroup' = identMat++instance (SmoothScalar s) => Morphism (Linear s) where+ DenseLinear f *** DenseLinear g = DenseLinear $ HMat.diagBlock [f,g]++instance (SmoothScalar s) => PreArrow (Linear s) where+ DenseLinear f &&& DenseLinear g = DenseLinear $ HMat.fromBlocks [[f], [g]]+ fst = lFst+ where lFst :: forall v w s+ . (FiniteDimensional v, FiniteDimensional w, Scalar v~s, Scalar w~s)+ => Linear s (v,w) v+ lFst = DenseLinear $ HMat.assoc (nv,n) 0 l+ where l = [ ((i,i), 1) | i<-[0.. nv-1] ]+ (Tagged nv) = dimension :: Tagged v Int+ (Tagged nw) = dimension :: Tagged w Int+ n = nv + nw+ snd = lSnd+ where lSnd :: forall v w s+ . (FiniteDimensional v, FiniteDimensional w, Scalar v~s, Scalar w~s)+ => Linear s (v,w) w+ lSnd = DenseLinear $ HMat.assoc (nw,n) 0 l+ where l = [ ((i,i+nv), 1) | i<-[0.. nw-1] ]+ (Tagged nv) = dimension :: Tagged v Int+ (Tagged nw) = dimension :: Tagged w Int+ n = nv + nw+ terminal = lTerminal+ where lTerminal :: forall v s . (FiniteDimensional v, Scalar v~s)+ => Linear s v (ZeroDim s)+ lTerminal = DenseLinear $ (0 HMat.>< n) []+ where (Tagged n) = dimension :: Tagged v Int++instance (SmoothScalar s) => EnhancedCat (->) (Linear s) where+ arr (DenseLinear mat) = fromPackedVector . HMat.app mat . asPackedVector+++++canonicalIdentityMatrix :: forall n v s+ . (KnownNat n, IsFreeSpace v, FreeDimension v ~ n, Scalar v ~ s)+ => Linear s v (FreeVect n s)+canonicalIdentityMatrix = DenseLinear $ HMat.ident n+ where (Tagged n) = theNatN :: Tagged n Int++-- | Class of spaces that directly represent a free vector space, i.e. that are simply+-- @n@-fold products of the base field.+-- This class basically contains 'ℝ', 'ℝ²', 'ℝ³' etc., in future also the complex and+-- probably integral versions.+class (FiniteDimensional v, KnownNat (FreeDimension v)) => IsFreeSpace v where+ type FreeDimension v :: Nat+ identityMatrix :: Isomorphism (Linear (Scalar v))+ v+ (FreeVect (FreeDimension v) (Scalar v))+ identityMatrix = fromInversePair emb proj+ where emb@(DenseLinear i) = canonicalIdentityMatrix+ proj = DenseLinear i++instance (KnownNat n, Num s, SmoothScalar s) => IsFreeSpace (FreeVect n s) where + type FreeDimension (FreeVect n s) = n+ identityMatrix = fromInversePair id id++instance IsFreeSpace ℝ where+ type FreeDimension ℝ = S Z+ +instance ( SmoothScalar s, IsFreeSpace v, Scalar v ~ s, FiniteDimensional s, s ~ Scalar s )+ => IsFreeSpace (v,s) where+ type FreeDimension (v,s) = S (FreeDimension v)++++class VectorSpace v => FreeTuple v where+ type Tuplity v :: Nat+ freeTuple :: Isomorphism (->) v (FreeVect (Tuplity v) (Scalar v))++#define FreeScalar(s) \+instance FreeTuple (s) where { \+ type Tuplity (s) = S Z; \+ freeTuple = fromInversePair (FreeVect . pure) (\(FreeVect v) -> v Arr.! 0); }++#define FreePair(s) \+FreeScalar(s); \+instance FreeTuple (s,s) where { \+ type Tuplity (s,s) = S(S Z); \+ freeTuple = fromInversePair (\(a,b) -> FreeVect $ Arr.fromList[a,b]) \+ (\(FreeVect v) -> (v Arr.! 0, v Arr.! 1)); }++#define FreeTriple(s) \+FreePair(s); \+instance FreeTuple (s,s,s) where { \+ type Tuplity (s,s,s) = S(S(S Z)); \+ freeTuple = fromInversePair (\(a,b,c) -> FreeVect $ Arr.fromList[a,b,c]) \+ (\(FreeVect v) -> (v Arr.! 0, v Arr.! 1, v Arr.! 2)); };\+instance FreeTuple (s,(s,s)) where { \+ type Tuplity (s,(s,s)) = S(S(S Z)); \+ freeTuple = fromInversePair (\(a,(b,c)) -> FreeVect $ Arr.fromList[a,b,c]) \+ (\(FreeVect v) -> (v Arr.! 0, (v Arr.! 1, v Arr.! 2))); };\+instance FreeTuple ((s,s),s) where { \+ type Tuplity ((s,s),s) = S(S(S Z)); \+ freeTuple = fromInversePair (\((a,b),c) -> FreeVect $ Arr.fromList[a,b,c]) \+ (\(FreeVect v) -> ((v Arr.! 0, v Arr.! 1), v Arr.! 2)); }++FreeTriple(ℝ)+FreeTriple(Int)++
Data/LinearMap/HerMetric.hs view
@@ -1,5 +1,6 @@ {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE TupleSections #-} {-# LANGUAGE TypeFamilies #-}@@ -16,32 +17,62 @@ HerMetric, HerMetric' -- * Evaluating metrics , metricSq, metricSq', metric, metric', metrics, metrics'- -- * Defining metrics by projectors+ -- * Defining metrics , projector, projector'- -- * Utility- , adjoint+ , euclideanMetric'+ -- * Metrics induce inner products+ , spanHilbertSubspace+ , spanSubHilbertSpace+ , IsFreeSpace+ -- * Utility for metrics , transformMetric, transformMetric' , dualiseMetric, dualiseMetric'- , HasMetric(..)+ , recipMetric, recipMetric'+ , eigenSpan, eigenSpan'+ , eigenCoSpan, eigenCoSpan'+ , metriScale', metriScale+ , adjoint+ -- * The dual-space class+ , HasMetric+ , HasMetric'(..) , (^<.>)- , metriScale, metriScale'+-- , riesz, riesz'+ -- * Fundamental requirements+ , MetricScalar+ , FiniteDimensional(..) ) where -import Prelude hiding ((^))- 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 Control.Applicative+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+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 <.>^, ^<.> @@ -56,30 +87,45 @@ -- -- Yet other possible interpretations of this type include /density matrix/ (as in -- quantum mechanics), /standard range of statistical fluctuations/, and /volume element/.-newtype HerMetric v = HerMetric { getHerMetric :: v :-* DualSpace v }+newtype HerMetric v = HerMetric {+ -- morally: @getHerMetric :: v :-* DualSpace v@.+ metricMatrix :: Maybe (HMat.Matrix (Scalar v)) -- @Nothing@ for zero metric.+ } +matrixMetric :: HasMetric v => HMat.Matrix (Scalar v) -> HerMetric v+matrixMetric = HerMetric . Just -instance HasMetric v => AdditiveGroup (HerMetric v) where- zeroV = HerMetric zeroV- negateV (HerMetric m) = HerMetric $ negateV m- HerMetric m ^+^ HerMetric n = HerMetric $ m ^+^ n+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 $ s *^ m + s *^ (HerMetric m) = HerMetric $ HMat.scale s <$> m -- | A metric on the dual space; equivalent to a linear mapping from the dual space -- to the original vector space. -- -- Prime-versions of the functions in this module target those dual-space metrics, so -- we can avoid some explicit handling of double-dual spaces.-newtype HerMetric' v = HerMetric' { dualMetric :: DualSpace v :-* v }+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' zeroV- negateV (HerMetric' m) = HerMetric' $ negateV m- HerMetric' m ^+^ HerMetric' n = HerMetric' $ m ^+^ n-instance (HasMetric v) => VectorSpace (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' $ s *^ m + s *^ (HerMetric' m) = HerMetric' $ HMat.scale s <$> m -- | A metric on @v@ that simply yields the squared overlap of a vector with the@@ -92,31 +138,52 @@ -- Metrics generated this way are positive definite if no negative coefficients have -- been introduced with the '*^' scaling operator or with '^-^'. projector :: HasMetric v => DualSpace v -> HerMetric v-projector u = HerMetric (linear $ \v -> u ^* (u<.>^v))+projector u = matrixMetric $ HMat.outer uDecomp uDecomp+ where uDecomp = asPackedVector u projector' :: HasMetric v => v -> HerMetric' v-projector' v = HerMetric' . linear $ \u -> v ^* (v^<.>u)+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 ++ -- | Evaluate a vector through a metric. For the canonical metric on a Hilbert space, -- this will be simply 'magnitudeSq'. metricSq :: HasMetric v => HerMetric v -> v -> Scalar v-metricSq (HerMetric m) v = lapply m v <.>^ 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' m) u = lapply m u ^<.> u+metricSq' (HerMetric' Nothing) _ = 0+metricSq' (HerMetric' (Just m)) u = uDecomp `HMat.dot` HMat.app m uDecomp+ where uDecomp = asPackedVector u -- | Evaluate a vector's “magnitude” through a metric. This assumes an actual -- mathematical metric, i.e. positive definite – otherwise the internally used -- square root may get negative arguments (though it can still produce results if the -- scalars are complex; however, complex spaces aren't supported yet). metric :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> v -> Scalar v-metric (HerMetric m) v = sqrt $ lapply m v <.>^ v+metric m = sqrt . metricSq m metric' :: (HasMetric v, Floating (Scalar v)) => HerMetric' v -> DualSpace v -> Scalar v-metric' (HerMetric' m) u = sqrt $ lapply m u ^<.> u+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++-- | “Anti-normalise” a vector: /multiply/ with its own norm, according to metric. metriScale :: (HasMetric v, Floating (Scalar v)) => HerMetric v -> v -> v metriScale m v = metric m v *^ v @@ -139,28 +206,101 @@ transformMetric :: (HasMetric v, HasMetric w, Scalar v ~ Scalar w) => (w :-* v) -> HerMetric v -> HerMetric w-transformMetric t (HerMetric m) = HerMetric $ adjoint t *.* m *.* t+transformMetric _ (HerMetric Nothing) = HerMetric Nothing+transformMetric t (HerMetric (Just m)) = matrixMetric $ HMat.tr tmat HMat.<> m HMat.<> tmat+ where tmat = asPackedMatrix t transformMetric' :: ( HasMetric v, HasMetric w, Scalar v ~ Scalar w ) => (v :-* w) -> HerMetric' v -> HerMetric' w-transformMetric' t (HerMetric' m)- = HerMetric' $ t *.* m *.* adjoint t+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, HasMetric (DualSpace v))- => HerMetric (DualSpace v) -> HerMetric' v-dualiseMetric (HerMetric m) = HerMetric' $ linear doubleDual' *.* m+-- | This doesn't really do anything at all, since @'HerMetric' v@ is essentially a+-- synonym for @'HerMetric' ('DualSpace' v)@.+dualiseMetric :: HasMetric v => HerMetric (DualSpace v) -> HerMetric' v+dualiseMetric (HerMetric m) = HerMetric' m -dualiseMetric' :: (HasMetric v, HasMetric (DualSpace v))- => HerMetric' v -> HerMetric (DualSpace v)-dualiseMetric' (HerMetric' m) = HerMetric $ linear doubleDual *.* m+dualiseMetric' :: HasMetric v => HerMetric' v -> HerMetric (DualSpace v)+dualiseMetric' (HerMetric' m) = HerMetric m +-- | The inverse mapping of a metric tensor. Since a metric maps from+-- a space to its dual, the inverse maps from the dual into the+-- (double-dual) space – i.e., it is a metric on the dual space.+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++++-- | The eigenbasis of a /positive definite/ metric, with each eigenvector scaled+-- to the square root of the eigenvalue.+-- +-- This constitutes, in a sense,+-- a decomposition of a metric into a set of 'projector'' vectors. If those+-- are 'sumV'ed again, the original metric is obtained. (This holds even for+-- non-Hilbert/Banach spaces, even though the concept of eigenbasis and+-- “scaled length” doesn't really makes sense then in the usual way!)+eigenSpan :: (HasMetric v, Scalar v ~ ℝ) => HerMetric' v -> [v]+eigenSpan (HerMetric' Nothing) = []+eigenSpan (HerMetric' (Just m)) = map fromPackedVector eigSpan+ where (μs,vsm) = HMat.eigSH m -- TODO: replace with `eigSH'`, which is unchecked+ -- (`HerMetric` is always Hermitian!)+ 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 -- TODO: replace with `eigSH'`, which is unchecked+ -- (`HerMetric` is always Hermitian!)+ 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 -- TODO: replace with `eigSH'`, which is unchecked+ -- (`HerMetric` is always Hermitian!)+ 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 -- TODO: replace with `eigSH'`, which is unchecked+ -- (`HerMetric` is always Hermitian!)+ eigSpan = zipWith (HMat.scale . recip . sqrt) (HMat.toList μs) (HMat.toColumns vsm)+++-- | Constraint that a space's scalars need to fulfill so it can be used for 'HerMetric'.+type MetricScalar s = ( SmoothScalar s+ , Ord s -- We really rather wouldn't require this...+ )+++type HasMetric v = (HasMetric' v, HasMetric' (DualSpace v), DualSpace (DualSpace v) ~ v)++ -- | While the main purpose of this class is to express 'HerMetric', it's actually -- all about dual spaces.-class ( HasBasis v, VectorSpace (Scalar v), HasTrie (Basis v)+class ( FiniteDimensional v, FiniteDimensional (DualSpace v) , VectorSpace (DualSpace v), HasBasis (DualSpace v)- , Scalar v ~ Scalar (DualSpace v), Basis v ~ Basis (DualSpace v) )- => HasMetric v where+ , MetricScalar (Scalar v), Scalar v ~ Scalar (DualSpace v)+ , Basis v ~ Basis (DualSpace v) )+ => HasMetric' v where -- | @'DualSpace' v@ is isomorphic to the space of linear functionals on @v@, i.e. -- @v ':-*' 'Scalar' v@.@@ -183,10 +323,10 @@ -- | While isomorphism between a space and its dual isn't generally canonical, -- the /double-dual/ space should be canonically isomorphic in pretty much -- all relevant cases. Indeed, it is recommended that they are the very same type;- -- the tuple instance actually assumes this to be able to offer an efficient- -- implementation (namely, 'id') of the isomorphisms.- doubleDual :: HasMetric (DualSpace v) => v -> DualSpace (DualSpace v)- doubleDual' :: HasMetric (DualSpace v) => DualSpace (DualSpace v) -> v+ -- this condition is enforced by the 'HasMetric' constraint (which is recommended+ -- over using 'HasMetric'' itself in signatures).+ doubleDual :: HasMetric' (DualSpace v) => v -> DualSpace (DualSpace v)+ doubleDual' :: HasMetric' (DualSpace v) => DualSpace (DualSpace v) -> v @@ -194,20 +334,36 @@ (^<.>) :: HasMetric v => v -> DualSpace v -> Scalar v ket ^<.> bra = bra <.>^ ket -instance (VectorSpace k) => HasMetric (ZeroDim k) where++euclideanMetric' :: forall v . (HasMetric v, InnerSpace v) => HerMetric v+euclideanMetric' = HerMetric . pure $ HMat.ident n+ where (Tagged n) = dimension :: Tagged v Int++-- -- | Associate a Hilbert space vector canonically with its dual-space counterpart,+-- -- as by the Riesz representation theorem.+-- -- +-- -- Note that usually, Hilbert spaces should just implement @DualSpace v ~ v@,+-- -- according to that same correspondence, so 'riesz' is essentially just a more explicit+-- -- (and less efficient) way of writing @'id' :: v -> DualSpace v'.+-- riesz :: (HasMetric v, InnerSpace v) => v -> DualSpace v+-- riesz v = functional (v<.>)+-- +-- riesz' :: (HasMetric v, InnerSpace v) => DualSpace v -> v+-- riesz' f = doubleDual' . functional (f<.>^)+++instance (MetricScalar k) => HasMetric' (ZeroDim k) where Origin<.>^Origin = zeroV functional _ = Origin doubleDual = id; doubleDual'= id-instance HasMetric Double where+instance HasMetric' Double where (<.>^) = (<.>) functional f = f 1 doubleDual = id; doubleDual'= id instance ( HasMetric v, HasMetric w, Scalar v ~ Scalar w- , HasMetric (DualSpace v), DualSpace (DualSpace v) ~ v- , HasMetric (DualSpace w), DualSpace (DualSpace w) ~ w- ) => HasMetric (v,w) where+ ) => HasMetric' (v,w) where type DualSpace (v,w) = (DualSpace v, DualSpace w)- (v,w)<.>^(v',w') = v<.>^v' ^+^ w<.>^w'+ (v,w)<.>^(v',w') = v<.>^v' + w<.>^w' functional f = (functional $ f . (,zeroV), functional $ f . (zeroV,)) doubleDual = id; doubleDual'= id @@ -225,8 +381,10 @@ -metrConst :: (HasMetric v, v ~ DualSpace v, Num (Scalar v)) => Scalar v -> HerMetric v-metrConst = HerMetric . linear . (*^)+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@@ -234,7 +392,7 @@ negate = negateV -- | This does /not/ work correctly if the metrics don't share an eigenbasis!- HerMetric m * HerMetric n = HerMetric $ m *.* n+ HerMetric m * HerMetric n = HerMetric $ liftA2 (HMat.<>) m n -- | Undefined, though it could actually be done. abs = error "abs undefined for HerMetric"@@ -243,7 +401,8 @@ metrNumFun :: (HasMetric v, v ~ Scalar v, v ~ DualSpace v, Num v) => (v -> v) -> HerMetric v -> HerMetric v-metrNumFun f (HerMetric m) = HerMetric . linear . (*^) . f $ lapply m 1+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@@ -267,3 +426,54 @@ 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 -> Hask.empty+ μ -> pure (v ^/ μ)++orthonormalPairsWith :: forall v . HasMetric v => HerMetric v -> [v] -> [(v, DualSpace v)]+orthonormalPairsWith met = mkON+ where mkON :: [v] -> [(v, DualSpace v)] -- Generalised Gram-Schmidt process+ 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+ +++spanHilbertSubspace :: forall 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@.+ -> Option (Embedding (Linear s) w v)+ -- ^ An embedding of the @n@-dimensional free subspace @w@ (if the given+ -- vectors actually span such a space) into the main space @v@.+ -- Regardless of the structure of @v@ (which doesn't need to have an+ -- inner product at all!), @w@ will be an 'InnerSpace' with the scalar+ -- product defined by the given metric.+spanHilbertSubspace met = emb . orthonormalPairsWith met+ where emb onb'+ | n'==n = return $ Embedding emb prj . arr identityMatrix+ | otherwise = Hask.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+++-- | Same as 'spanHilbertSubspace', but with the standard 'euclideanMetric' (i.e., the+-- basis vectors will be orthonormal in the usual sense, in both @w@ and @v@).+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'+
Data/List/FastNub.hs view
@@ -16,15 +16,17 @@ fastNubBy :: (a->a->Ordering) -> [a] -> [a] fastNubBy _ [] = [] fastNubBy _ [e] = [e]-fastNubBy cmp es = merge(fastNubBy cmp lhs)(fastNubBy cmp rhs)+fastNubBy cmp es = fnubMergeBy cmp (fastNubBy cmp lhs) (fastNubBy cmp rhs) where (lhs,rhs) = splitAt (length es `quot` 2) es- merge [] rs = rs- merge ls [] = ls- merge (l:ls) (r:rs) = case cmp l r of- LT -> l : merge ls (r:rs)- GT -> r : merge (l:ls) rs- EQ -> merge (l:ls) rs +fnubMergeBy :: (a->a->Ordering) -> [a] -> [a] -> [a]+fnubMergeBy _ [] rs = rs+fnubMergeBy _ ls [] = ls+fnubMergeBy cmp (l:ls) (r:rs) = case cmp l r of+ LT -> l : fnubMergeBy cmp ls (r:rs)+ GT -> r : fnubMergeBy cmp (l:ls) rs+ EQ -> fnubMergeBy cmp (l:ls) rs+ -- | Like 'fastNubBy', but doesn't just discard duplicates but \"merges\" them. -- @'fastNubBy' cmp = cmp `'fastNubByWith'` 'const'@. fastNubByWith :: (a->a->Ordering) -> (a->a->a) -> [a] -> [a]@@ -41,3 +43,23 @@ sfGroupBy :: (a->a->Ordering) -> [a] -> [[a]] sfGroupBy cmp = fastNubByWith (cmp`on`head) (++) . map(:[])+++++fnubConcatBy :: (a->a->Ordering) -> [[a]] -> [a]+fnubConcatBy cmp = foldr (fnubMergeBy cmp) [] . map (fastNubBy cmp)++fnubConcat :: FastNub a => [[a]] -> [a]+fnubConcat = foldr (fnubMergeBy compare) [] . map fastNub++fnubConcatMap :: FastNub b => (a -> [b]) -> [a] -> [b]+fnubConcatMap f = fnubConcat . map f++fnubIntersect :: FastNub a => [a] -> [a] -> [a]+fnubIntersect xs ys = fis (fastNub xs) (fastNub ys)+ where fis [] _ = []+ fis _ [] = []+ fis (x:xs) (y:ys) | x<y = fis xs (y:ys)+ | x>y = fis (x:xs) ys+ | otherwise = x : fis xs ys
Data/Manifold.hs view
@@ -29,7 +29,7 @@ {-# LANGUAGE RecordWildCards #-} -module Data.Manifold (module Data.Manifold, module Data.Manifold.Types) where+module Data.Manifold (module Data.Manifold, module Data.Manifold.Types.Primitive) where import Data.List import Data.Maybe@@ -41,7 +41,7 @@ import Data.Basis import Data.Complex hiding (magnitude) import Data.Void-import Data.Manifold.Types+import Data.Manifold.Types.Primitive import qualified Prelude @@ -204,6 +204,7 @@ +type EuclidSpace v = (HasBasis v, EqFloating(Scalar v), Eq v) isInUpperHemi :: EuclidSpace v => v -> Bool isInUpperHemi v = (snd . head) (decompose v) >= 0
Data/Manifold/PseudoAffine.hs view
@@ -23,6 +23,7 @@ {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE LiberalTypeSynonyms #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE TupleSections #-}@@ -36,17 +37,34 @@ module Data.Manifold.PseudoAffine ( -- * Manifold class- PseudoAffine(..)+ Manifold+ , Semimanifold(..)+ , PseudoAffine(..)+ , Metric, Metric', euclideanMetric -- * Regions within a manifold , Region -- * Hierarchy of manifold-categories+ -- ** Everywhere differentiable functions , Differentiable- , PWDiffable, RWDiffable+ -- ** Almost everywhere diff'able funcs+ , PWDiffable+ -- ** Region-wise defined diff'able funcs+ , RWDiffable+ -- * Helper constraints+ , RealDimension, AffineManifold+ , LinearManifold+ , WithField+ , HilbertSpace+ , EuclidSpace+ -- * Misc+ , palerp ) 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)@@ -55,14 +73,18 @@ import Data.VectorSpace import Data.LinearMap import Data.LinearMap.HerMetric-import Data.MemoTrie (HasTrie)+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+import Data.Manifold.Types.Primitive +import Data.CoNat++import qualified Numeric.LinearAlgebra.HMatrix as HMat+ import qualified Prelude import Control.Category.Constrained.Prelude hiding ((^))@@ -74,9 +96,42 @@ infix 6 .-~.-infixl 6 .+~^+infixl 6 .+~^, .-~^ --- | 'PseudoAffine' is intended as an alternative class for 'Data.Manifold.Manifold's.+class (AdditiveGroup (Needle x)) => Semimanifold x where+ -- | 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+ -- line segments (aka vectors) between two points, i.e. the same as 'Diff'.+ -- The 'AffineManifold' constraint makes that requirement explicit.+ -- + -- This space should be isomorphic to the tangent space (and is in fact+ -- used somewhat synonymously).+ type Needle x :: *+ + -- | Generalised translation operation.+ (.+~^) :: x -> Needle x -> x+ + -- | Shorthand for @\\p v -> p .+~^ 'negateV' v@, which should obey the /asymptotic/ law+ -- + -- @+ -- p .-~^ v .+~^ v ≅ p+ -- @+ -- + -- Meaning: if @v@ is scaled down with sufficiently small factors /η/, then+ -- the difference @(p.-~^v.+~^v) .-~. p@ should scale down even faster:+ -- 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+ 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.)+-- -- 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.@@ -86,83 +141,158 @@ -- designated origin, a pseudo-affine space can have nontrivial topology on the global -- scale, and yet be used in practically the same way as an affine space. At least the -- usual spheres and tori make good instances, perhaps the class is in fact equivalent to--- /parallelisable manifolds/.-class PseudoAffine x where- type PseudoDiff x :: *- (.-~.) :: x -> x -> Option (PseudoDiff x)- (.+~^) :: x -> PseudoDiff x -> x+-- 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+ -- | 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+ -- + -- @+ -- p .+~^ (q.-~.p) ≡ q+ -- @+ -- + -- should hold, at least save for floating-point precision limits etc..+ (.-~.) :: x -> 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 -type LocallyScalable s x = ( PseudoAffine x, (PseudoDiff x) ~ PseudoDiff x- , HasMetric (PseudoDiff x)- , DualSpace (PseudoDiff x) ~ DualSpace (PseudoDiff x)- , HasMetric (DualSpace (PseudoDiff x))- , PseudoDiff x ~ DualSpace (DualSpace (PseudoDiff x))- , s ~ Scalar (PseudoDiff x)- , s ~ Scalar (DualSpace (PseudoDiff x)) )-type LinearManifold s x = ( PseudoAffine x, PseudoDiff x ~ x- , HasMetric x, HasMetric (DualSpace x)- , DualSpace (DualSpace x) ~ x- , s ~ Scalar x, s ~ Scalar (DualSpace x) )-type RealDimension r = ( PseudoAffine r, PseudoDiff r ~ r+type LocallyScalable s x = ( PseudoAffine x, (Needle x) ~ Needle x+ , HasMetric (Needle x)+ , DualSpace (Needle x) ~ DualSpace (Needle x)+ , s ~ Scalar (Needle x) )++-- | Basically just an “updated” version of the 'VectorSpace' class.+-- Every vector space is a manifold, this constraint makes it explicit.+-- +-- (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 )++-- | Require some constraint on a manifold, and also fix the type of the manifold's+-- underlying field. For example, @WithField ℝ 'HilbertSpace' v@ constrains+-- @v@ to be a real (i.e., 'Double'-) Hilbert space.+-- Note that for this to compile, you will in+-- general need the @-XLiberalTypeSynonyms@ extension (except if the constraint+-- is an actual type class (like 'Manifold'): only those can always be partially+-- applied, for @type@ constraints this is by default not allowed).+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 , HasMetric r, DualSpace r ~ r, Scalar r ~ r , RealFloat r ) +-- | The 'AffineSpace' class plus manifold constraints.+type AffineManifold m = ( PseudoAffine 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+-- also a manifold.+-- +-- (Stricly speaking, that doesn't have much to do with the completeness criterion;+-- 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) ) -palerp :: (PseudoAffine x, VectorSpace (PseudoDiff x))- => x -> x -> Option (Scalar (PseudoDiff x) -> x)+-- | An euclidean space is a real affine space whose tangent space is a Hilbert space.+type EuclidSpace x = ( AffineManifold x, InnerSpace (Diff x)+ , DualSpace (Diff x) ~ Diff x, Floating (Scalar (Diff x)) )++euclideanMetric :: EuclidSpace x => Tagged x (Metric x)+euclideanMetric = Tagged euclideanMetric'+++-- | The word “metric” is used in the sense as in general relativity. Cf. 'HerMetric'.+type Metric x = HerMetric (Needle x)+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 #define deriveAffine(t) \-instance PseudoAffine t where { \- type PseudoDiff t = Diff t; \- a.-~.b = pure (a.-.b); \- (.+~^) = (.+^) }+instance Semimanifold (t) where { \+ type Needle (t) = Diff (t); \+ (.+~^) = (.+^) }; \+instance PseudoAffine (t) where { \+ a.-~.b = pure (a.-.b); } deriveAffine(Double) deriveAffine(Rational) +instance Semimanifold (ZeroDim k) where+ type Needle (ZeroDim k) = ZeroDim k+ Origin .+~^ Origin = Origin+ Origin .-~^ Origin = Origin instance PseudoAffine (ZeroDim k) where- type PseudoDiff (ZeroDim k) = ZeroDim k Origin .-~. Origin = pure Origin- Origin .+~^ Origin = Origin++instance (Semimanifold a, Semimanifold b) => Semimanifold (a,b) where+ type Needle (a,b) = (Needle a, Needle b)+ (a,b).+~^(v,w) = (a.+~^v, b.+~^w)+ (a,b).-~^(v,w) = (a.-~^v, b.-~^w) instance (PseudoAffine a, PseudoAffine b) => PseudoAffine (a,b) where- type PseudoDiff (a,b) = (PseudoDiff a, PseudoDiff b) (a,b).-~.(c,d) = liftA2 (,) (a.-~.c) (b.-~.d)- (a,b).+~^(v,w) = (a.+~^v, b.+~^w)++instance (Semimanifold a, Semimanifold b, Semimanifold c) => Semimanifold (a,b,c) where+ type Needle (a,b,c) = (Needle a, Needle b, Needle 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) instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) => PseudoAffine (a,b,c) where- type PseudoDiff (a,b,c) = (PseudoDiff a, PseudoDiff b, PseudoDiff c) (a,b,c).-~.(d,e,f) = liftA3 (,,) (a.-~.d) (b.-~.e) (c.-~.f)- (a,b,c).+~^(v,w,x) = (a.+~^v, b.+~^w, c.+~^x) +instance (MetricScalar a, KnownNat n) => Semimanifold (FreeVect n a) where+ type Needle (FreeVect n a) = FreeVect n a+ (.+~^) = (.+^)+instance (MetricScalar a, KnownNat n) => PseudoAffine (FreeVect n a) where+ a.-~.b = pure (a.-.b) ++instance Semimanifold S⁰ where+ type Needle S⁰ = ℝ⁰+ p .+~^ Origin = p+ p .-~^ Origin = p+instance PseudoAffine S⁰ where+ PositiveHalfSphere .-~. PositiveHalfSphere = pure Origin+ NegativeHalfSphere .-~. NegativeHalfSphere = pure Origin+ _ .-~. _ = Option Nothing++instance Semimanifold S¹ where+ type Needle S¹ = ℝ+ S¹ φ₀ .+~^ δφ+ | φ' < 0 = S¹ $ φ' + tau+ | otherwise = S¹ $ φ'+ where φ' = toS¹range $ φ₀ + δφ instance PseudoAffine S¹ where- type PseudoDiff S¹ = ℝ S¹ φ₁ .-~. S¹ φ₀ | δφ > pi = pure (δφ - 2*pi) | δφ < (-pi) = pure (δφ + 2*pi) | otherwise = pure δφ where δφ = φ₁ - φ₀- S¹ φ₀ .+~^ δφ- | φ' < 0 = S¹ $ φ' + tau- | otherwise = S¹ $ φ'- where φ' = (φ₀ + δφ)`mod'`tau +instance Semimanifold S² where+ type Needle S² = ℝ²+ S² ϑ₀ φ₀ .+~^ δv+ | ϑ₀ < pi/2 = sphereFold PositiveHalfSphere $ ϑ₀*^embed(S¹ φ₀) ^+^ δv+ | otherwise = sphereFold NegativeHalfSphere $ (pi-ϑ₀)*^embed(S¹ φ₀) ^+^ δv instance PseudoAffine S² where- type PseudoDiff S² = ℝ² S² ϑ₁ φ₁ .-~. S² ϑ₀ φ₀ | ϑ₀ < pi/2 = pure ( ϑ₁*^embed(S¹ φ₁) ^-^ ϑ₀*^embed(S¹ φ₀) ) | otherwise = pure ( (pi-ϑ₁)*^embed(S¹ φ₁) ^-^ (pi-ϑ₀)*^embed(S¹ φ₀) )- S² ϑ₀ φ₀ .+~^ δv- | ϑ₀ < pi/2 = sphereFold PositiveHalfSphere $ ϑ₀*^embed(S¹ φ₀) ^+^ δv- | otherwise = sphereFold NegativeHalfSphere $ (pi-ϑ₀)*^embed(S¹ φ₀) ^+^ δv sphereFold :: S⁰ -> ℝ² -> S² sphereFold hfSphere v- | ϑ₀ > pi = S² (inv $ tau - ϑ₀) ((φ₀+pi)`mod'`tau)+ | ϑ₀ > pi = S² (inv $ tau - ϑ₀) (toS¹range $ φ₀+pi) | otherwise = S² (inv ϑ₀) φ₀ where S¹ φ₀ = coEmbed v ϑ₀ = magnitude v `mod'` tau@@ -170,15 +300,51 @@ NegativeHalfSphere -> pi - ϑ +instance Semimanifold ℝP² where+ type Needle ℝP² = ℝ²+ ℝP² r₀ φ₀ .+~^ (δr, δφ)+ | r₀ > 1/2 = case r₀ + δr of+ r₁ | r₁ > 1 -> ℝP² (2-r₁) (toS¹range $ φ₀+δφ+pi)+ | otherwise -> ℝP² r₁ (toS¹range $ φ₀+δφ)+ ℝP² r₀ φ₀ .+~^ δxy = let v = r₀*^embed(S¹ φ₀) ^+^ δxy+ S¹ φ₁ = coEmbed v+ r₁ = magnitude v `mod'` 1+ in ℝP² r₁ φ₁ +instance PseudoAffine ℝP² where+ ℝP² r₁ φ₁ .-~. ℝP² r₀ φ₀+ | r₀ > 1/2 = pure `id` case φ₁-φ₀ of+ δφ | δφ > 3*pi/2 -> ( r₁ - r₀, δφ - 2*pi)+ | δφ < -3*pi/2 -> ( r₁ - r₀, δφ + 2*pi)+ | δφ > pi/2 -> (2-r₁ - r₀, δφ - pi )+ | δφ < -pi/2 -> (2-r₁ - r₀, δφ + pi )+ | otherwise -> ( r₁ - r₀, δφ )+ | otherwise = pure ( r₁*^embed(S¹ φ₁) ^-^ r₀*^embed(S¹ φ₀) ) -tau :: Double++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₀+ ++++tau :: ℝ tau = 2 * pi +toS¹range :: ℝ -> ℝ+toS¹range φ = (φ+pi)`mod'`tau - pi -type LinDevPropag d c = HerMetric (PseudoDiff c) -> HerMetric (PseudoDiff d)+type LinDevPropag d c = Metric c -> Metric d dev_ε_δ :: RealDimension a => (a -> a) -> LinDevPropag a a@@ -219,11 +385,11 @@ -- overlap from exceeding one; this makes the concept actually work on general manifolds.) newtype Differentiable s d c = Differentiable { runDifferentiable ::- d -> ( c, PseudoDiff d :-* PseudoDiff c, LinDevPropag d c ) }+ d -> ( c, Needle d :-* Needle c, LinDevPropag d c ) } type (-->) = Differentiable ℝ -instance (VectorSpace s) => Category (Differentiable s) where+instance (MetricScalar s) => Category (Differentiable s) where type Object (Differentiable s) o = LocallyScalable s o id = Differentiable $ \x -> (x, idL, const zeroV) Differentiable f . Differentiable g = Differentiable $@@ -235,7 +401,7 @@ in (z, f'*.*g', devfg) -instance (VectorSpace s) => Cartesian (Differentiable s) where+instance (MetricScalar s) => Cartesian (Differentiable s) where type UnitObject (Differentiable s) = ZeroDim s swap = Differentiable $ \(x,y) -> ((y,x), lSwap, const zeroV) where lSwap = linear swap@@ -249,7 +415,7 @@ where lRegroup = linear regroup' -instance (VectorSpace s) => Morphism (Differentiable s) where+instance (MetricScalar s) => Morphism (Differentiable s) where Differentiable f *** Differentiable g = Differentiable h where h (x,y) = ((fx, gy), lPar, devfg) where (fx, f', devf) = f x@@ -263,7 +429,7 @@ lcofst = linear (,zeroV); lcosnd = linear (zeroV,) -instance (VectorSpace s) => PreArrow (Differentiable s) where+instance (MetricScalar s) => PreArrow (Differentiable s) where terminal = Differentiable $ \_ -> (Origin, zeroV, const zeroV) fst = Differentiable $ \(x,_) -> (x, lfst, const zeroV) where lfst = linear fst@@ -279,7 +445,7 @@ lcofst = linear (,zeroV); lcosnd = linear (zeroV,) -instance (VectorSpace s) => WellPointed (Differentiable s) where+instance (MetricScalar s) => WellPointed (Differentiable s) where unit = Tagged Origin globalElement x = Differentiable $ \Origin -> (x, zeroV, const zeroV) const x = Differentiable $ \_ -> (x, zeroV, const zeroV)@@ -288,30 +454,30 @@ type DfblFuncValue s = GenericAgent (Differentiable s) -instance (VectorSpace s) => HasAgent (Differentiable s) where+instance (MetricScalar s) => HasAgent (Differentiable s) where alg = genericAlg ($~) = genericAgentMap-instance (VectorSpace s) => CartesianAgent (Differentiable s) where+instance (MetricScalar s) => CartesianAgent (Differentiable s) where alg1to2 = genericAlg1to2 alg2to1 = genericAlg2to1 alg2to2 = genericAlg2to2-instance (VectorSpace s)+instance (MetricScalar s) => PointAgent (DfblFuncValue s) (Differentiable s) a x where point = genericPoint -actuallyLinear :: ( LinearManifold s x, LinearManifold s y )+actuallyLinear :: ( WithField s LinearManifold x, WithField s LinearManifold y ) => (x:-*y) -> Differentiable s x y actuallyLinear f = Differentiable $ \x -> (lapply f x, f, const zeroV) -actuallyAffine :: ( LinearManifold s x, LinearManifold s y )+actuallyAffine :: ( WithField s LinearManifold x, WithField s LinearManifold y ) => y -> (x:-*y) -> Differentiable s x y actuallyAffine y₀ f = Differentiable $ \x -> (y₀ ^+^ lapply f x, f, const zeroV) dfblFnValsFunc :: ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s d- , v ~ PseudoDiff c, v' ~ PseudoDiff c'+ , v ~ Needle c, v' ~ Needle c' , ε ~ HerMetric v, ε ~ HerMetric v' ) => (c' -> (c, v':-*v, ε->ε)) -> DfblFuncValue s d c' -> DfblFuncValue s d c dfblFnValsFunc f = (Differentiable f $~)@@ -319,7 +485,7 @@ dfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' , LocallyScalable s d- , v ~ PseudoDiff c, v' ~ PseudoDiff c', v'' ~ PseudoDiff c''+ , v ~ Needle c, v' ~ Needle c', v'' ~ Needle c'' , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) -> DfblFuncValue s d c' -> DfblFuncValue s d c'' -> DfblFuncValue s d c@@ -346,7 +512,7 @@ -instance (LinearManifold s v, LocallyScalable s a, Floating s)+instance (WithField s LinearManifold v, LocallyScalable s a, Floating s) => AdditiveGroup (DfblFuncValue s a v) where zeroV = point zeroV (^+^) = dfblFnValsCombine $ \a b -> (a^+^b, lPlus, const zeroV)@@ -387,7 +553,7 @@ -- roots, but the square root of a nontrivial-vector-space metric requires -- an eigenbasis transform, which we have not implemented yet. -- --- instance (LinearManifold s v, LocallyScalable s a, Floating s)+-- instance (WithField s LinearManifold v, LocallyScalable s a, Floating s) -- => VectorSpace (DfblFuncValue s a v) where -- type Scalar (DfblFuncValue s a v) = DfblFuncValue s a (Scalar v) -- (*^) = dfblFnValsCombine $ \μ v -> (μ*^v, lScl, \ε -> (ε ^* sqrt 2, ε ^* sqrt 2))@@ -569,7 +735,7 @@ gpwDfblFnValsFunc :: ( RealDimension s , LocallyScalable s c, LocallyScalable s c', LocallyScalable s d- , v ~ PseudoDiff c, v' ~ PseudoDiff c'+ , v ~ Needle c, v' ~ Needle c' , ε ~ HerMetric v, ε ~ HerMetric v' ) => (c' -> (c, v':-*v, ε->ε)) -> PWDfblFuncValue s d c' -> PWDfblFuncValue s d c gpwDfblFnValsFunc f = (PWDiffable (\_ -> (GlobalRegion, Differentiable f)) $~)@@ -577,7 +743,7 @@ gpwDfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' , LocallyScalable s d, RealDimension s- , v ~ PseudoDiff c, v' ~ PseudoDiff c', v'' ~ PseudoDiff c''+ , v ~ Needle c, v' ~ Needle c', v'' ~ Needle c'' , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) -> PWDfblFuncValue s d c' -> PWDfblFuncValue s d c'' -> PWDfblFuncValue s d c@@ -604,7 +770,7 @@ lcosnd = linear(zeroV,) -instance (LinearManifold s v, LocallyScalable s a, RealDimension s)+instance (WithField s LinearManifold v, LocallyScalable s a, RealDimension s) => AdditiveGroup (PWDfblFuncValue s a v) where zeroV = point zeroV (^+^) = gpwDfblFnValsCombine $ \a b -> (a^+^b, lPlus, const zeroV)@@ -677,7 +843,7 @@ -- need to exhaustively 'isNaN'-check all results...) -- -- @--- hb :: RWDiffable R R R+-- hb :: RWDiffable ℝ ℝ ℝ -- hb = alg (\\p -> - p * logBase 2 p - (1-p) * logBase 2 (1-p) ) -- @ newtype RWDiffable s d c@@ -784,7 +950,7 @@ grwDfblFnValsFunc :: ( RealDimension s , LocallyScalable s c, LocallyScalable s c', LocallyScalable s d- , v ~ PseudoDiff c, v' ~ PseudoDiff c'+ , v ~ Needle c, v' ~ Needle c' , ε ~ HerMetric v, ε ~ HerMetric v' ) => (c' -> (c, v':-*v, ε->ε)) -> RWDfblFuncValue s d c' -> RWDfblFuncValue s d c grwDfblFnValsFunc f = (RWDiffable (\_ -> (GlobalRegion, pure (Differentiable f))) $~)@@ -792,7 +958,7 @@ grwDfblFnValsCombine :: forall d c c' c'' v v' v'' ε ε' ε'' s. ( LocallyScalable s c, LocallyScalable s c', LocallyScalable s c'' , LocallyScalable s d, RealDimension s- , v ~ PseudoDiff c, v' ~ PseudoDiff c', v'' ~ PseudoDiff c''+ , v ~ Needle c, v' ~ Needle c', v'' ~ Needle c'' , ε ~ HerMetric v , ε' ~ HerMetric v' , ε'' ~ HerMetric v'', ε~ε', ε~ε'' ) => ( c' -> c'' -> (c, (v',v''):-*v, ε -> (ε',ε'')) ) -> RWDfblFuncValue s d c' -> RWDfblFuncValue s d c'' -> RWDfblFuncValue s d c@@ -824,7 +990,7 @@ -instance (LinearManifold s v, LocallyScalable s a, RealDimension s)+instance (WithField s LinearManifold v, LocallyScalable s a, RealDimension s) => AdditiveGroup (RWDfblFuncValue s a v) where zeroV = point zeroV (^+^) = grwDfblFnValsCombine $ \a b -> (a^+^b, lPlus, const zeroV)@@ -1011,5 +1177,10 @@ -- Empirical, with epsEst upper bound. - ++++++
+ Data/Manifold/TreeCover.hs view
@@ -0,0 +1,910 @@+-- |+-- Module : Data.Manifold.TreeCover+-- 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 RecordWildCards #-}+{-# LANGUAGE DataKinds #-}+++module Data.Manifold.TreeCover (+ -- * Shades + Shade, shadeCtr, shadeExpanse, fullShade, pointsShades+ -- * Shade trees+ , ShadeTree(..), fromLeafPoints+ -- * Simple view helpers+ , onlyNodes, onlyLeaves+ -- ** Auxiliary types+ , SimpleTree, Trees, NonEmptyTree, GenericTree(..)+ -- * Misc+ , sShSaw, chainsaw, HasFlatView(..)+ -- ** Triangulation-builders+ , TriangBuild, doTriangBuild, singleFullSimplex, autoglueTriangulation+ , AutoTriang, elementaryTriang, breakdownAutoTriang+ ) 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.SimplicialComplex+import Data.Manifold.Types+import Data.Manifold.Types.Primitive ((^))+import Data.Manifold.PseudoAffine+ +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)+import Control.Arrow.Constrained+import Control.Monad.Constrained hiding (forM)+import Data.Foldable.Constrained++import GHC.Generics (Generic)+++-- | Possibly / Partially / asymPtotically singular metric.+data PSM x = PSM {+ psmExpanse :: !(Metric' x)+ , relevantEigenspan :: ![DualSpace (Needle x)]+ }+ ++-- | A 'Shade' is a very crude description of a region within a manifold. It+-- can be interpreted as either an ellipsoid shape, or as the Gaussian peak+-- of a normal distribution (use <http://hackage.haskell.org/package/manifold-random>+-- for actually sampling from that distribution).+-- +-- 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) }++instance (AffineManifold x) => Semimanifold (Shade x) where+ type Needle (Shade x) = Diff x+ 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++subshadeId' :: WithField ℝ Manifold x+ => x -> NonEmpty (DualSpace (Needle x)) -> x -> (Int, HourglassBulb)+subshadeId' c expvs x = case x .-~. c of+ Option (Just v) -> let (iu,vl) = maximumBy (comparing $ abs . snd)+ $ zip [0..] (map (v <.>^) $ NE.toList expvs)+ in (iu, if vl>0 then UpperBulb else LowerBulb)+ _ -> (-1, error "Trying to obtain the subshadeId of a point not actually included in the shade.")++subshadeId :: WithField ℝ Manifold x => Shade x -> x -> (Int, HourglassBulb)+subshadeId (Shade c expa) = subshadeId' c . NE.fromList $ eigenCoSpan expa+ +++-- | Attempt to find a 'Shade' that “covers” the given points.+-- At least in an affine space (and thus locally in any manifold), this can be used to+-- estimate the parameters of a normal distribution from which some points were+-- sampled.+-- +-- For /nonconnected/ manifolds it will be necessary to yield separate shades+-- for each connected component. And for an empty input list, there is no shade!+-- Hence the list result.+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 + Option (Just δ) -> (acc .+~^ δ^/i, (p:rb, nr))+ _ -> (acc, (rb, p:nr)) )+ (p₀, mempty)+ ( zip [1..] $ p₀:psr )++pointsShades' :: WithField ℝ Manifold x => Metric' x -> [x] -> [([x], Shade x)]+pointsShades' _ [] = []+pointsShades' minExt ps = case expa of + Option (Just e) -> (ps, fullShade ctr e)+ : pointsShades' minExt unreachable+ _ -> pointsShades' minExt inc'd+ ++ pointsShades' minExt unreachable+ where (ctr,(inc'd,unreachable)) = pseudoECM $ NE.fromList ps+ expa = ( (^+^minExt) . (^/ fromIntegral(length ps)) . sumV . map projector' )+ <$> mapM (.-~.ctr) ps+ ++minusLogOcclusion :: (PseudoAffine x, HasMetric (Needle 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+ _ -> 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 )+ => Shade x -> x -> s+occlusion (Shade p₀ δ) = occ+ where occ p = case p .-~. p₀ of+ Option(Just vd) -> exp . negate $ metricSq δinv vd+ _ -> zeroV+ δinv = recipMetric δ++++-- | Hourglass as the geometric shape (two opposing ~conical volumes, sharing+-- only a single point in the middle); has nothing to do with time.+data Hourglass s = Hourglass { upperBulb, lowerBulb :: !s }+ deriving (Generic, Hask.Functor, Hask.Foldable)+instance (NFData s) => NFData (Hourglass s)+instance (Semigroup s) => Semigroup (Hourglass s) where+ Hourglass u l <> Hourglass u' l' = Hourglass (u<>u') (l<>l')+ sconcat hgs = let (us,ls) = NE.unzip $ (upperBulb&&&lowerBulb) <$> hgs+ in Hourglass (sconcat us) (sconcat ls)+instance (Monoid s, Semigroup s) => Monoid (Hourglass s) where+ mempty = Hourglass mempty mempty; mappend = (<>)+ mconcat hgs = let (us,ls) = unzip $ (upperBulb&&&lowerBulb) <$> hgs+ in Hourglass (mconcat us) (mconcat ls)+instance Hask.Applicative Hourglass where+ pure x = Hourglass x x+ Hourglass f g <*> Hourglass x y = Hourglass (f x) (g y)+instance Foldable Hourglass (->) (->) where+ ffoldl f (x, Hourglass a b) = f (f(x,a), b)+ foldMap f (Hourglass a b) = f a `mappend` f b++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+oneBulb LowerBulb f (Hourglass u l) = Hourglass u (f l)++++data ShadeTree x = PlainLeaves [x]+ | DisjointBranches !Int (NonEmpty (ShadeTree x))+ | OverlappingBranches !Int !(Shade x) (NonEmpty (DBranch x))+ deriving (Generic)+ +data DBranch' x c = DBranch { boughDirection :: !(DualSpace (Needle x))+ , boughContents :: !(Hourglass c) }+ deriving (Generic, Hask.Functor, Hask.Foldable)+type DBranch x = DBranch' x (ShadeTree x)++newtype DBranches' x c = DBranches (NonEmpty (DBranch' x c))+ deriving (Generic, Hask.Functor, Hask.Foldable)++-- ^ /Unsafe/: this assumes the direction information of both containers to be equivalent.+instance (Semigroup c) => Semigroup (DBranches' x c) where+ DBranches b1 <> DBranches b2 = DBranches $ NE.zipWith (\(DBranch d1 c1) (DBranch _ c2)+ -> DBranch d1 $ c1<>c2 ) b1 b2+ +++instance (NFData x) => NFData (ShadeTree x) where+ rnf (PlainLeaves xs) = rnf xs+ rnf (DisjointBranches n bs) = n `seq` rnf (NE.toList bs)+ rnf (OverlappingBranches n sh bs) = n `seq` sh `seq` rnf (NE.toList bs)+instance (NFData x) => NFData (DBranch x)+ +-- | Experimental. There might be a more powerful instance possible.+instance (AffineManifold x) => Semimanifold (ShadeTree x) where+ type Needle (ShadeTree x) = Diff x+ PlainLeaves xs .+~^ v = PlainLeaves $ (.+^v)<$>xs + OverlappingBranches n sh br .+~^ v+ = OverlappingBranches n (sh.+~^v)+ $ fmap (\(DBranch d c) -> DBranch d $ (.+~^v)<$>c) br+ DisjointBranches n br .+~^ v = DisjointBranches n $ (.+~^v)<$>br++-- | WRT union.+instance WithField ℝ Manifold x => Semigroup (ShadeTree x) where+ PlainLeaves [] <> t = t+ t <> PlainLeaves [] = t+ t <> s = fromLeafPoints $ onlyLeaves t ++ onlyLeaves s+ -- Could probably be done more efficiently+ sconcat = mconcat . NE.toList+instance WithField ℝ Manifold x => Monoid (ShadeTree x) where+ mempty = PlainLeaves []+ mappend = (<>)+ mconcat l = case filter ne l of+ [] -> mempty+ [t] -> t+ l' -> fromLeafPoints $ onlyLeaves =<< l'+ where ne (PlainLeaves []) = False; ne _ = True+++-- | Build a really 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+-- +-- > plotWindow [ plot . onlyNodes $ fromLeafPoints testPts4+-- > , plot [circle 0.06 & moveTo p & fc green :: PlainGraphics | p <- testPts4] ]+-- @+-- +-- <<images/examples/simple-2d-ShadeTree.png>>+fromLeafPoints :: forall x. WithField ℝ Manifold x => [x] -> ShadeTree x+fromLeafPoints = go zeroV+ where go :: Metric' x -> [x] -> ShadeTree x+ go preShExpa = \xs -> case pointsShades' (preShExpa^/10) xs of+ [] -> mempty+ [(_,rShade)] -> let trials = sShIdPartition rShade xs+ in case reduce rShade trials of+ Just redBrchs+ -> OverlappingBranches+ (length xs) rShade+ (branchProc (shadeExpanse rShade) redBrchs)+ _ -> PlainLeaves xs+ partitions -> DisjointBranches (length xs)+ . NE.fromList+ $ map (\(xs',pShade) -> go zeroV xs') partitions+ where + branchProc redSh = fmap (fmap $ go redSh)+ + reduce :: Shade x -> NonEmpty (DBranch' x [x])+ -> Maybe (NonEmpty (DBranch' x [x]))+ reduce sh@(Shade c _) brCandidates+ = case findIndex deficient cards of+ Just i | (DBranch _ reBr, o:ok)+ <- amputateId i (NE.toList brCandidates)+ -> reduce sh+ $ sShIdPartition' c (fold reBr) (o:|ok)+ | otherwise -> Nothing+ _ -> Just brCandidates+ where (cards, maxCard) = (NE.toList &&& maximum')+ $ fmap (fmap length . boughContents) brCandidates+ deficient (Hourglass u l) = any (\c -> c^2 <= maxCard + 1) [u,l]+ maximum' = maximum . NE.toList . fmap (\(Hourglass u l) -> max u l)+++sShIdPartition' :: WithField ℝ Manifold x+ => x -> [x] -> NonEmpty (DBranch' x [x])->NonEmpty (DBranch' x [x])+sShIdPartition' c xs st+ = foldr (\p -> let (i,h) = ssi p+ in asList $ update_nth (\(DBranch d c)+ -> DBranch d (oneBulb h (p:) c))+ i )+ st xs+ where ssi = subshadeId' c (boughDirection<$>st)+sShIdPartition :: WithField ℝ Manifold x => Shade x -> [x] -> NonEmpty (DBranch' x [x])+sShIdPartition (Shade c expa) xs+ | b:bs <- [DBranch v mempty | v <- eigenCoSpan expa]+ = sShIdPartition' c xs $ b:|bs+ ++asList :: ([a]->[b]) -> NonEmpty a->NonEmpty b+asList f = NE.fromList . f . NE.toList++update_nth :: (a->a) -> Int -> [a] -> [a]+update_nth _ n l | n<0 = l+update_nth f 0 (c:r) = f c : r+update_nth f n [] = []+update_nth f n (l:r) = l : update_nth f (n-1) r+++amputateId :: Int -> [a] -> (a,[a])+amputateId i l = let ([a],bs) = amputateIds [i] l in (a, bs)++deleteIds :: [Int] -> [a] -> [a]+deleteIds kids = snd . amputateIds kids++amputateIds :: [Int] -- ^ Sorted list of non-negative indices to extract+ -> [a] -- ^ Input list+ -> ([a],[a]) -- ^ (Extracted elements, remaining elements)+amputateIds = go 0+ where go _ _ [] = ([],[])+ go _ [] l = ([],l)+ go i (k:ks) (x:xs)+ | i==k = first (x:) $ go (i+1) ks xs+ | otherwise = second (x:) $ go (i+1) (k:ks) xs+++++sortByKey :: Ord a => [(a,b)] -> [b]+sortByKey = map snd . sortBy (comparing fst)++++++ ++-- simplexFaces :: forall n x . Simplex (S n) x -> Triangulation n x+-- simplexFaces (Simplex p (ZeroSimplex q)) = TriangVertices $ Arr.fromList [p, q]+-- simplexFaces splx = carpent splx $ TriangVertices ps+-- where ps = Arr.fromList $ p : splxVertices qs+-- where carpent (ZeroSimplex (Simplex p qs@(Simplex _ _))+-- | Triangulation es <- simplexFaces qs = TriangSkeleton $ Simplex p <$> es+++++newtype BaryCoords n = BaryCoords { getBaryCoordsTail :: FreeVect n ℝ }++instance (KnownNat n) => AffineSpace (BaryCoords n) where+ type Diff (BaryCoords n) = FreeVect n ℝ+ BaryCoords v .-. BaryCoords w = v ^-^ w+ BaryCoords v .+^ w = BaryCoords $ v ^+^ w+instance (KnownNat n) => Semimanifold (BaryCoords n) where+ type Needle (BaryCoords n) = FreeVect n ℝ+ (.+~^) = (.+^)+instance (KnownNat n) => PseudoAffine (BaryCoords n) where+ (.-~.) = pure .: (.-.)++getBaryCoords :: BaryCoords n -> ℝ ^ S n+getBaryCoords (BaryCoords (FreeVect bcs)) = FreeVect $ (1 - Arr.sum bcs) `Arr.cons` bcs+ +getBaryCoords' :: BaryCoords n -> [ℝ]+getBaryCoords' (BaryCoords (FreeVect bcs)) = 1 - Arr.sum bcs : Arr.toList bcs++getBaryCoord :: BaryCoords n -> Int -> ℝ+getBaryCoord (BaryCoords (FreeVect bcs)) 0 = 1 - Arr.sum bcs+getBaryCoord (BaryCoords (FreeVect bcs)) i = case bcs Arr.!? i of+ Just a -> a+ _ -> 0++mkBaryCoords :: KnownNat n => ℝ ^ S n -> BaryCoords n+mkBaryCoords (FreeVect bcs) = BaryCoords $ FreeVect (Arr.tail bcs) ^/ Arr.sum bcs++mkBaryCoords' :: KnownNat n => [ℝ] -> Option (BaryCoords n)+mkBaryCoords' bcs = fmap (BaryCoords . (^/sum bcs)) . freeVector . Arr.fromList $ tail bcs++newtype ISimplex n x = ISimplex { iSimplexBCCordEmbed :: Embedding (->) (BaryCoords n) x }+++++data TriangBuilder n x where+ TriangVerticesSt :: [x] -> TriangBuilder Z x+ TriangBuilder :: Triangulation (S n) x+ -> [x]+ -> [(Simplex n x, [x] -> Option x)]+ -> TriangBuilder (S n) x++++-- startTriangulation :: forall n x . (KnownNat n, WithField ℝ Manifold x)+-- => ISimplex n x -> TriangBuilder n x+-- startTriangulation ispl@(ISimplex emb) = startWith $ fromISimplex ispl+-- where startWith (ZeroSimplex p) = TriangVerticesSt [p]+-- startWith s@(Simplex _ _)+-- = TriangBuilder (Triangulation [s])+-- (splxVertices s)+-- [ (s', expandInDir j)+-- | j<-[0..n]+-- | s' <- getTriangulation $ simplexFaces s ]+-- where expandInDir j xs = case sortBy (comparing snd) $ filter ((> -1) . snd) xs_bc of+-- ((x, q) : _) | q<0 -> pure x+-- _ -> Hask.empty+-- where xs_bc = map (\x -> (x, getBaryCoord (emb >-$ x) j)) xs+-- (Tagged n) = theNatN :: Tagged n Int++-- extendTriangulation :: forall n x . (KnownNat n, WithField ℝ Manifold x)+-- => [x] -> TriangBuilder n x -> TriangBuilder n x+-- extendTriangulation xs (TriangBuilder tr tb te) = foldr tryex (TriangBuilder tr tb []) te+-- where tryex (bspl, expd) (TriangBuilder (Triangulation tr') tb' te')+-- | Option (Just fav) <- expd xs+-- = let snew = Simplex fav bspl+-- in TriangBuilder (Triangulation $ snew:tr') (fav:tb') undefined++ +bottomExtendSuitability :: (KnownNat n, WithField ℝ Manifold x)+ => ISimplex (S n) x -> x -> ℝ+bottomExtendSuitability (ISimplex emb) x = case getBaryCoord (emb >-$ x) 0 of+ 0 -> 0+ r -> - recip r++optimalBottomExtension :: (KnownNat n, WithField ℝ Manifold x)+ => ISimplex (S n) x -> [x] -> Option Int+optimalBottomExtension s xs+ = case filter ((>0).snd)+ $ zipWith ((. bottomExtendSuitability s) . (,)) [0..] xs of+ [] -> Hask.empty+ qs -> pure . fst . maximumBy (comparing snd) $ qs+++simplexPlane :: forall n x . (KnownNat n, WithField ℝ Manifold x)+ => Metric x -> Simplex n x -> Embedding (Linear ℝ) (FreeVect n ℝ) (Needle x)+simplexPlane m s = embedding+ where bc = barycenter s+ spread = init . map ((.-~.bc) >>> \(Option (Just v)) -> v) $ splxVertices s+ embedding = case spanHilbertSubspace m spread of+ (Option (Just e)) -> e+ _ -> error "Trying to obtain simplexPlane from zero-volume\+ \ simplex (which cannot span sufficient basis vectors)."++++-- simplexShade :: forall x n . (KnownNat n, WithField ℝ Manifold x)+barycenter :: forall x n . (KnownNat n, WithField ℝ Manifold x) => Simplex n x -> x+barycenter = bc + where bc (ZS x) = x+ bc (x :<| xs') = x .+~^ sumV [x'–x | x'<-splxVertices xs'] ^/ (n+1)+ + Tagged n = theNatN :: Tagged n ℝ+ x' – x = case x'.-~.x of {Option(Just v)->v}++toISimplex :: forall x n . (KnownNat n, WithField ℝ Manifold x)+ => Metric x -> Simplex n x -> ISimplex n x+toISimplex m s = ISimplex $ fromEmbedProject fromBrc toBrc+ where bc = barycenter s+ (Embedding emb (DenseLinear prj))+ = simplexPlane m s+ (r₀:rs) = [ prj HMat.#> asPackedVector v+ | x <- splxVertices s, let (Option (Just v)) = x.-~.bc ]+ tmat = HMat.inv $ HMat.fromColumns [ r - r₀ | r<-rs ] + toBrc x = case x.-~.bc of+ Option (Just v) -> let rx = prj HMat.#> asPackedVector v - r₀+ in finalise $ tmat HMat.#> rx+ finalise v = case freeVector $ HMat.toList v of+ Option (Just bv) -> BaryCoords bv+ fromBrc bccs = bc .+~^ (emb $ v)+ where v = linearCombo $ (fromPackedVector r₀, b₀) : zip (fromPackedVector<$>rs) bs+ (b₀:bs) = getBaryCoords' bccs++fromISimplex :: forall x n . (KnownNat n, WithField ℝ Manifold x)+ => ISimplex n x -> Simplex n x+fromISimplex (ISimplex emb) = s+ where (Option (Just s))+ = makeSimplex' [ emb $-> jOnly+ | j <- [0..n]+ , let (Option (Just jOnly)) = mkBaryCoords' [ if k==j then 1 else 0+ | k<-[0..n] ]+ ]+ (Tagged n) = theNatN :: Tagged n Int++iSimplexSideViews :: ∀ n x . KnownNat n => ISimplex n x -> [ISimplex n x]+iSimplexSideViews = \(ISimplex is)+ -> take (n+1) $ [ISimplex $ rot j is | j<-[0..] ]+ where rot j (Embedding emb proj)+ = Embedding ( emb . mkBaryCoords . freeRotate j . getBaryCoords )+ ( mkBaryCoords . freeRotate (n-j) . getBaryCoords . proj )+ (Tagged n) = theNatN :: Tagged n Int+++type FullTriang t n x = TriangT t n x+ (State (Map.Map (SimplexIT t n x) (ISimplex n x)))++type TriangBuild t n x = TriangT t (S n) x+ ( State (Map.Map (SimplexIT t n x) (Metric x, ISimplex (S n) x) ))++doTriangBuild :: KnownNat n => (∀ t . TriangBuild t n x ()) -> [Simplex (S n) x]+doTriangBuild t = runIdentity (fst <$>+ doTriangT (unliftInTriangT (`evalStateT`mempty) t >> simplexITList >>= mapM lookSimplex))++singleFullSimplex :: ∀ t n x . (KnownNat n, WithField ℝ Manifold x)+ => ISimplex n x -> FullTriang t n x (SimplexIT t n x)+singleFullSimplex is = do+ frame <- disjointSimplex (fromISimplex is)+ lift . modify' $ Map.insert frame is+ return frame+ +fullOpenSimplex :: ∀ t n x . (KnownNat n, WithField ℝ Manifold x)+ => Metric x -> Simplex (S n) x -> TriangBuild t n x [SimplexIT t n x]+fullOpenSimplex m s = do+ let is = toISimplex m s+ frame <- disjointSimplex (fromISimplex is)+ fsides <- toList <$> lookSplxFacesIT frame+ lift . forM (zip fsides $ iSimplexSideViews is)+ $ \(fside,is') -> modify' $ Map.insert fside (m,is')+ return fsides+++hypotheticalSimplexScore :: ∀ t n n' x . (KnownNat n', WithField ℝ Manifold x, n~S n')+ => SimplexIT t Z x+ -> SimplexIT t n x+ -> TriangBuild t n x ( Option Double )+hypotheticalSimplexScore p b = do+ altViews :: [(SimplexIT t Z x, SimplexIT t n x)] <- do+ pSups <- lookSupersimplicesIT p+ nOpts <- forM pSups $ \psup -> fmap (fmap $ \((bq,_p), _b') -> (bq,psup))+ $ distinctSimplices b psup+ return $ catOptions nOpts+ scores <- forM ((p,b) :| altViews) $ \(p',b') -> do+ x <- lookVertexIT p'+ q <- lift $ Map.lookup b' <$> get+ return $ case q of+ Just(_,is) | s<-bottomExtendSuitability is x, s>0+ -> pure s+ _ -> Hask.empty+ return . fmap sum $ Hask.sequence scores++spanSemiOpenSimplex :: ∀ t n n' x . (KnownNat n', WithField ℝ Manifold x, n~S n')+ => SimplexIT t Z x -- ^ Tip of the desired simplex.+ -> SimplexIT t n x -- ^ Base of the desired simplex.+ -> TriangBuild t n x [SimplexIT t n x]+ -- ^ Return the exposed faces of the new simplices.+spanSemiOpenSimplex p b = do+ m <- lift $ fst <$> (Map.!b) <$> get+ neighbours <- filterM isAdjacent =<< lookSupersimplicesIT p+ let bs = b:|neighbours+ frame <- webinateTriang p b+ backSplx <- lookSimplex frame+ let iSplx = toISimplex m backSplx+ fsides <- toList <$> lookSplxFacesIT frame+ let sviews = filter (not . (`elem`bs) . fst) $ zip fsides (iSimplexSideViews iSplx)+ lift . forM sviews $ \(fside,is') -> modify' $ Map.insert fside (m,is')+ lift . Hask.forM_ bs $ \fside -> modify' $ Map.delete fside+ return $ fst <$> sviews+ where isAdjacent = fmap (isJust . getOption) . sharedBoundary b++multiextendTriang :: ∀ t n n' x . (KnownNat n', WithField ℝ Manifold x, n~S n')+ => [SimplexIT t Z x] -> TriangBuild t n x ()+multiextendTriang vs = do+ ps <- mapM lookVertexIT vs+ sides <- lift $ Map.toList <$> get+ forM_ sides $ \(f,(m,s)) ->+ case optimalBottomExtension s ps of+ Option (Just c) -> spanSemiOpenSimplex (vs !! c) f+ _ -> return []++-- | BUGGY: this does connect the supplied triangulations, but it doesn't choose+-- the right boundary simplices yet. Probable cause: inconsistent internal+-- numbering of the subsimplices.+autoglueTriangulation :: ∀ t n n' n'' x+ . (KnownNat n'', WithField ℝ Manifold x, n~S n', n'~S n'')+ => (∀ t' . TriangBuild t' n' x ()) -> TriangBuild t n' x ()+autoglueTriangulation tb = do+ mbBounds <- Map.toList <$> lift get+ mps <- pointsOfSurf mbBounds+ + WriterT gbBounds <- liftInTriangT $ mixinTriangulation tb'+ lift . forM_ gbBounds $ \(i,ms) -> do+ modify' $ Map.insert i ms+ gps <- pointsOfSurf gbBounds+ + autoglue mps gbBounds+ autoglue gps mbBounds+ + where tb' :: ∀ s . TriangT s n x Identity+ (WriterT (Metric x, ISimplex n x) [] (SimplexIT s n' x))+ tb' = unliftInTriangT (`evalStateT`mempty) $+ tb >> (WriterT . Map.toList) <$> lift get+ + pointsOfSurf s = fnubConcatMap Hask.toList <$> forM s (lookSplxVerticesIT . fst)+ + autoglue :: [SimplexIT t Z x] -> [(SimplexIT t n' x, (Metric x, ISimplex n x))]+ -> TriangBuild t n' x ()+ autoglue vs sides = do+ forM_ sides $ \(f,_) -> do+ possibs <- forM vs $ \p -> fmap(p,) <$> hypotheticalSimplexScore p f+ case catOptions possibs of+ [] -> return ()+ qs -> do+ spanSemiOpenSimplex (fst `id` maximumBy (comparing $ snd) qs) f+ return ()+++data AutoTriang n x where+ AutoTriang :: { getAutoTriang :: ∀ t . TriangBuild t n x () } -> AutoTriang (S n) x++instance (KnownNat n, WithField ℝ Manifold x) => Semigroup (AutoTriang (S (S n)) x) where+ (<>) = autoTriangMappend++autoTriangMappend :: ∀ n n' n'' x . ( KnownNat n'', n ~ S n', n' ~ S n''+ , WithField ℝ Manifold x )+ => AutoTriang n x -> AutoTriang n x -> AutoTriang n x+AutoTriang a `autoTriangMappend` AutoTriang b = AutoTriang c+ where c :: ∀ t . TriangBuild t n' x ()+ c = a >> autoglueTriangulation b++elementaryTriang :: ∀ n n' x . (KnownNat n', n~S n', WithField ℝ EuclidSpace x)+ => Simplex n x -> AutoTriang n x+elementaryTriang t = AutoTriang (fullOpenSimplex m t >> return ())+ where (Tagged m) = euclideanMetric :: Tagged x (Metric x)++breakdownAutoTriang :: ∀ n n' x . (KnownNat n', n ~ S n') => AutoTriang n x -> [Simplex n x]+breakdownAutoTriang (AutoTriang t) = doTriangBuild t+ + +-- where tr :: Triangulation n x+-- outfc :: Map.Map (SimplexIT t n' x) (Metric x, ISimplex n x)+-- (((), tr), outfc) = runState (doTriangT tb') mempty+-- tb' :: ∀ t' . TriangT t' n x +-- ( State ( Map.Map (SimplexIT t' n' x)+-- (Metric x, ISimplex n x) ) ) ()+-- tb' = tb+ + + + ++-- primitiveTriangulation :: forall x n . (KnownNat n,WithField ℝ Manifold x)+-- => [x] -> Triangulation n x+-- primitiveTriangulation xs = head $ build <$> buildOpts+-- where build :: ([x], [x]) -> Triangulation n x+-- build (mainVerts, sideVerts) = Triangulation [mainSplx]+-- where (Option (Just mainSplx)) = makeSimplex mainVerts+-- -- mainFaces = Map.fromAscList . zip [0..] . getTriangulation+-- -- $ simplexFaces mainSplx+-- buildOpts = partitionsOfFstLength n xs+-- (Tagged n) = theNatN :: Tagged n Int+ +partitionsOfFstLength :: Int -> [a] -> [([a],[a])]+partitionsOfFstLength 0 l = [([],l)]+partitionsOfFstLength n [] = []+partitionsOfFstLength n (x:xs) = first (x:) <$> partitionsOfFstLength (n-1) xs+ ++ second (x:) <$> partitionsOfFstLength n xs++splxVertices :: Simplex n x -> [x]+splxVertices (ZS x) = [x]+splxVertices (x :<| s') = x : splxVertices s'++++-- triangulate :: forall x n . (KnownNat n, WithField ℝ Manifold x)+-- => ShadeTree x -> Triangulation n x+-- triangulate (DisjointBranches _ brs)+-- = Triangulation $ Hask.foldMap (getTriangulation . triangulate) brs+-- triangulate (PlainLeaves xs) = primitiveTriangulation xs++-- triangBranches :: WithField ℝ Manifold x+-- => ShadeTree x -> Branchwise x (Triangulation x) n+-- triangBranches _ = undefined+-- +-- tringComplete :: WithField ℝ Manifold x+-- => Triangulation x (n-1) -> Triangulation x n -> Triangulation x n+-- tringComplete (Triangulation trr) (Triangulation tr) = undefined+-- where +-- bbSimplices = Map.fromList [(i, Left s) | s <- tr | i <- [0::Int ..] ]+-- bbVertices = [(i, splxVertices s) | s <- tr | i <- [0::Int ..] ]+-- + +++++-- |+-- @+-- 'SimpleTree' x ≅ Maybe (x, 'Trees' x)+-- @+type SimpleTree = GenericTree Maybe []+-- |+-- @+-- 'Trees' x ≅ [(x, 'Trees' x)]+-- @+type Trees = GenericTree [] []+-- |+-- @+-- 'NonEmptyTree' x ≅ (x, 'Trees' x)+-- @+type NonEmptyTree = GenericTree NonEmpty []+ +newtype GenericTree c b x = GenericTree { treeBranches :: c (x,GenericTree b b x) }+ deriving (Hask.Functor)+instance (Hask.MonadPlus c) => Semigroup (GenericTree c b x) where+ GenericTree b1 <> GenericTree b2 = GenericTree $ Hask.mplus b1 b2+instance (Hask.MonadPlus c) => Monoid (GenericTree c b x) where+ mempty = GenericTree Hask.mzero+ mappend = (<>)+deriving instance Show (c (x, GenericTree b b x)) => Show (GenericTree c b x)++-- | Imitate the specialised 'ShadeTree' structure with a simpler, generic tree.+onlyNodes :: WithField ℝ Manifold x => ShadeTree x -> Trees x+onlyNodes (PlainLeaves []) = GenericTree []+onlyNodes (PlainLeaves ps) = let (ctr,_) = pseudoECM $ NE.fromList ps+ in GenericTree [ (ctr, GenericTree $ (,mempty) <$> ps) ]+onlyNodes (DisjointBranches _ brs) = Hask.foldMap onlyNodes brs+onlyNodes (OverlappingBranches _ (Shade ctr _) brs)+ = GenericTree [ (ctr, Hask.foldMap (Hask.foldMap onlyNodes) brs) ]+++-- | Left (and, typically, also right) inverse of 'fromLeafNodes'.+onlyLeaves :: WithField ℝ Manifold x => ShadeTree x -> [x]+onlyLeaves tree = dismantle tree []+ where dismantle (PlainLeaves xs) = (xs++)+ dismantle (OverlappingBranches _ _ brs)+ = foldr ((.) . dismantle) id $ Hask.foldMap (Hask.toList) brs+ dismantle (DisjointBranches _ brs) = foldr ((.) . dismantle) id $ NE.toList brs+++++++++data Sawbones x = Sawbones { sawnTrunk1, sawnTrunk2 :: [x]->[x]+ , sawdust1, sawdust2 :: [x] }+instance Semigroup (Sawbones x) where+ Sawbones st11 st12 sd11 sd12 <> Sawbones st21 st22 sd21 sd22+ = Sawbones (st11.st21) (st12.st22) (sd11<>sd21) (sd12<>sd22)+instance Monoid (Sawbones x) where+ mempty = Sawbones id id [] []+ mappend = (<>)+++chainsaw :: WithField ℝ Manifold x => Cutplane x -> ShadeTree x -> Sawbones x+chainsaw cpln (PlainLeaves xs) = Sawbones (sd1++) (sd2++) sd2 sd1+ where (sd1,sd2) = partition (\x -> sideOfCut cpln x == Option(Just PositiveHalfSphere)) xs+chainsaw cpln (DisjointBranches _ brs) = Hask.foldMap (chainsaw cpln) brs+chainsaw cpln (OverlappingBranches _ (Shade _ bexpa) brs) = Sawbones t1 t2 d1 d2+ where (Sawbones t1 t2 subD1 subD2)+ = Hask.foldMap (Hask.foldMap (chainsaw cpln) . boughContents) brs+ [d1,d2] = map (foldl' go [] . foci) [subD1, subD2]+ where go d' (dp,dqs) = case fathomCD dp of+ Option (Just dpCD) | not $ any (shelter dpCD) dqs+ -> dp:d' -- dp is close enough to cut plane to make dust.+ _ -> d' -- some dq is actually closer than the cut plane => discard dp.+ where shelter dpCutDist dq = case ptsDist dp dq of+ Option (Just d) -> d < abs dpCutDist+ _ -> False+ ptsDist = fmap (metric $ recipMetric bexpa) .: (.-~.)+ fathomCD = fathomCutDistance cpln bexpa+ ++type DList x = [x]->[x]+ +data DustyEdges x = DustyEdges { sawChunk :: DList x, chunkDust :: DBranches' x [x] }+instance Semigroup (DustyEdges x) where+ DustyEdges c1 d1 <> DustyEdges c2 d2 = DustyEdges (c1.c2) (d1<>d2)++data Sawboneses x = SingleCut (Sawbones x)+ | Sawboneses (DBranches' x (DustyEdges x))+ deriving (Generic)+instance Semigroup (Sawboneses x) where+ SingleCut c <> SingleCut d = SingleCut $ c<>d+ Sawboneses c <> Sawboneses d = Sawboneses $ c<>d++++-- | Saw a tree into the domains covered by the respective branches of another tree.+sShSaw :: WithField ℝ Manifold x+ => ShadeTree x -- ^ “Reference tree”, defines the cut regions.+ -- Must be at least one level of 'OverlappingBranches' deep.+ -> ShadeTree x -- ^ Tree to take the actual contents from.+ -> Sawboneses x -- ^ All points within each region, plus those from the+ -- boundaries of each neighbouring region.+sShSaw (OverlappingBranches _ (Shade sh _) (DBranch dir _ :| [])) src+ = SingleCut $ chainsaw (Cutplane sh $ stiefel1Project dir) src+sShSaw (OverlappingBranches _ (Shade cctr _) cbrs) (PlainLeaves xs)+ = Sawboneses . DBranches $ NE.fromList ngbsAdded+ where brsEmpty = fmap (\(DBranch dir _)-> DBranch dir mempty) cbrs+ srcDistrib = sShIdPartition' cctr xs brsEmpty+ ngbsAdded = fmap (\(DBranch dir (Hourglass u l), othrs)+ -> let [allOthr,allOthr']+ = map (DBranches . NE.fromList)+ [othrs, fmap (\(DBranch d' o)+ ->DBranch(negateV d') o) othrs]+ in DBranch dir $ Hourglass (DustyEdges (u++) allOthr)+ (DustyEdges (l++) allOthr')+ ) $ foci (NE.toList srcDistrib)+sShSaw cuts@(OverlappingBranches _ (Shade sh _) cbrs)+ (OverlappingBranches _ (Shade _ bexpa) brs)+ = Sawboneses . DBranches $ ftr'd+ where Option (Just (Sawboneses (DBranches recursed)))+ = Hask.foldMap (Hask.foldMap (pure . sShSaw cuts) . boughContents) brs+ ftr'd = fmap (\(DBranch dir1 ds) -> DBranch dir1 $ fmap (+ \(DustyEdges bk (DBranches dds))+ -> DustyEdges bk . DBranches $ fmap (obsFilter dir1) dds+ ) ds ) recursed+ obsFilter dir1 (DBranch dir2 (Hourglass pd2 md2))+ = DBranch dir2 $ Hourglass pd2' md2'+ where cpln cpSgn = Cutplane sh . stiefel1Project $ dir1 ^+^ cpSgn*^dir2+ [pd2', md2'] = zipWith (occl . cpln) [-1, 1] [pd2, md2] + occl cpl = foldl' go [] . foci+ where go d' (dp,dqs) = case fathomCD dp of+ Option (Just dpCD) | not $ any (shelter dpCD) dqs+ -> dp:d'+ _ -> d'+ where shelter dpCutDist dq = case ptsDist dp dq of+ Option (Just d) -> d < abs dpCutDist+ _ -> False+ ptsDist = fmap (metric $ recipMetric bexpa) .: (.-~.)+ fathomCD = fathomCutDistance cpl bexpa+sShSaw _ _ = error "`sShSaw` is not supposed to cut anything else but `OverlappingBranches`"+++foci :: [a] -> [(a,[a])]+foci [] = []+foci (x:xs) = (x,xs) : fmap (second (x:)) (foci xs)+ ++(.:) :: (c->d) -> (a->b->c) -> a->b->d +(.:) = (.) . (.)+++catOptions :: [Option a] -> [a]+catOptions = catMaybes . map getOption++++class HasFlatView f where+ type FlatView f x+ flatView :: f x -> FlatView f x+ superFlatView :: f x -> [[x]]+ +instance HasFlatView Sawbones where+ type FlatView Sawbones x = [([x],[[x]])]+ flatView (Sawbones t1 t2 d1 d2) = [(t1[],[d1]), (t2[],[d2])]+ superFlatView = foldMap go . flatView+ where go (t,ds) = t : ds++instance HasFlatView Sawboneses where+ type FlatView Sawboneses x = [([x],[[x]])]+ flatView (SingleCut (Sawbones t1 t2 d1 d2)) = [(t1[],[d1]), (t2[],[d2])]+ flatView (Sawboneses (DBranches bs)) = + [ (m[], NE.toList ds >>= \(DBranch _ (Hourglass u' l')) -> [u',l'])+ | (DBranch _ (Hourglass u l)) <- NE.toList bs+ , (DustyEdges m (DBranches ds)) <- [u,l]+ ]+ superFlatView = foldMap go . flatView+ where go (t,ds) = t : ds+
Data/Manifold/Types.hs view
@@ -7,17 +7,21 @@ -- Stability : experimental -- Portability : portable -- +-- Several commonly-used manifolds, represented in some simple way as Haskell+-- data types. All these are in the 'PseudoAffine' class. {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-}--- {-# LANGUAGE OverlappingInstances #-}+{-# LANGUAGE CPP #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE LiberalTypeSynonyms #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE TupleSections #-}+{-# LANGUAGE MultiWayIf #-} {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE PatternGuards #-} {-# LANGUAGE TypeOperators #-}@@ -25,16 +29,52 @@ {-# LANGUAGE RecordWildCards #-} -module Data.Manifold.Types where+module Data.Manifold.Types (+ -- * Index / ASCII names+ Real0, Real1, RealPlus, Real2, Real3+ , Sphere0, Sphere1, Sphere2+ , Projective1, Projective2+ , Disk1, Disk2, Cone, OpenCone+ -- * Linear manifolds+ , ZeroDim(..)+ , ℝ⁰, ℝ, ℝ², ℝ³+ -- * Hyperspheres+ -- ** General form: Stiefel manifolds+ , Stiefel1, stiefel1Project, stiefel1Embed+ -- ** Specific examples+ , HasUnitSphere(..)+ , S⁰(..), S¹(..), S²(..)+ -- * Projective spaces+ , ℝP¹, ℝP²(..)+ -- * Intervals\/disks\/cones+ , D¹(..), D²(..)+ , ℝay+ , CD¹(..), Cℝay(..)+ -- * Misc+ -- * Cut-planes+ , Cutplane(..)+ , fathomCutDistance, sideOfCut+ ) where import Data.VectorSpace import Data.AffineSpace+import Data.MemoTrie (HasTrie(..)) import Data.Basis-import Data.Complex hiding (magnitude)+import Data.Fixed import Data.Void+import Data.Tagged import Data.Monoid+import Data.Semigroup+import qualified Numeric.LinearAlgebra.HMatrix as HMat+import qualified Data.Vector.Generic as Arr+import qualified Data.Vector +import Data.Manifold.Types.Primitive+import Data.Manifold.PseudoAffine+import Data.LinearMap.HerMetric+import Data.VectorSpace.FiniteDimensional+ import qualified Prelude import Control.Category.Constrained.Prelude hiding ((^))@@ -42,99 +82,204 @@ import Control.Monad.Constrained import Data.Foldable.Constrained +#define deriveAffine(c,t) \+instance (c) => Semimanifold (t) where { \+ type Needle (t) = Diff (t); \+ (.+~^) = (.+^) }; \+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+s1bTrie = \f -> St1BTrie $ fmap (f . Stiefel1Basis) allIs+ where (Tagged d) = dimension :: Tagged v Int+ allIs = Arr.fromList [0 .. d-2] +instance FiniteDimensional v => HasTrie (Stiefel1Basis v) where+ data (Stiefel1Basis v :->: a) = St1BTrie ( Array a )+ trie = s1bTrie; untrie (St1BTrie a) (Stiefel1Basis i) = a Arr.! i+ enumerate (St1BTrie a) = Arr.ifoldr (\i x l -> (Stiefel1Basis i,x):l) [] a -type EuclidSpace v = (HasBasis v, EqFloating(Scalar v), Eq v)-type EqFloating f = (Eq f, Ord f, Floating f)+type Array = Data.Vector.Vector +instance(SmoothScalar(Scalar v),FiniteDimensional v)=>AdditiveGroup(Stiefel1Needle v) where+ Stiefel1Needle v ^+^ Stiefel1Needle w = Stiefel1Needle $ v + w+ zeroV = s1nZ; negateV (Stiefel1Needle v) = Stiefel1Needle $ negate v+s1nZ :: forall v. FiniteDimensional v => Stiefel1Needle v+s1nZ=Stiefel1Needle .HMat.fromList$replicate(d-1)0 where(Tagged d)=dimension::Tagged v Int +instance (SmoothScalar(Scalar v),FiniteDimensional v)=>VectorSpace(Stiefel1Needle v) where+ type Scalar (Stiefel1Needle v) = Scalar v+ μ *^ Stiefel1Needle v = Stiefel1Needle $ HMat.scale μ v -data GraphWindowSpec = GraphWindowSpec {- lBound, rBound, bBound, tBound :: Double- , xResolution, yResolution :: Int- }+instance (SmoothScalar (Scalar v), FiniteDimensional v)=>HasBasis (Stiefel1Needle v) where+ type Basis (Stiefel1Needle v) = Stiefel1Basis v+ basisValue = s1bV+ decompose (Stiefel1Needle v) = zipWith ((,).Stiefel1Basis) [0..] $ HMat.toList v+ decompose' (Stiefel1Needle v) (Stiefel1Basis i) = v HMat.! i+s1bV :: forall v b. FiniteDimensional v => Stiefel1Basis v -> Stiefel1Needle v+s1bV = \(Stiefel1Basis i) -> Stiefel1Needle+ $ HMat.fromList [ if k==i then 1 else 0 | k<-[0..d-2] ]+ where (Tagged d) = dimension :: Tagged v Int +instance (SmoothScalar (Scalar v), FiniteDimensional v)+ => FiniteDimensional (Stiefel1Needle v) where+ dimension = s1nD+ basisIndex = Tagged $ \(Stiefel1Basis i) -> i+ indexBasis = Tagged Stiefel1Basis+ fromPackedVector = Stiefel1Needle+ asPackedVector = getStiefel1Tangent+s1nD :: forall v. FiniteDimensional v => Tagged (Stiefel1Needle v) Int+s1nD = Tagged (d - 1) where (Tagged d) = dimension :: Tagged v Int +instance (SmoothScalar (Scalar v), FiniteDimensional v)+ => AffineSpace (Stiefel1Needle v) where+ type Diff (Stiefel1Needle v) = Stiefel1Needle v+ (.+^) = (^+^)+ (.-.) = (^-^) +deriveAffine((SmoothScalar (Scalar v), FiniteDimensional v), Stiefel1Needle v) -data ZeroDim k = Origin deriving(Eq, Show)-instance Monoid (ZeroDim k) where- mempty = Origin- mappend Origin Origin = Origin-instance AdditiveGroup (ZeroDim k) where- zeroV = Origin- Origin ^+^ Origin = Origin- negateV Origin = Origin-instance VectorSpace (ZeroDim k) where- type Scalar (ZeroDim k) = k- _ *^ Origin = Origin-instance HasBasis (ZeroDim k) where- type Basis (ZeroDim k) = Void- basisValue = absurd- decompose Origin = []- decompose' Origin = absurd+instance (MetricScalar (Scalar v), FiniteDimensional v)+ => HasMetric' (Stiefel1Needle v) where+ type DualSpace (Stiefel1Needle v) = Stiefel1Needle v+ Stiefel1Needle v <.>^ Stiefel1Needle w = HMat.dot v w + functional = s1nF+ doubleDual = id; doubleDual' = id+s1nF :: forall v. FiniteDimensional v => (Stiefel1Needle v->Scalar v)->Stiefel1Needle v+s1nF = \f -> Stiefel1Needle $ HMat.fromList [f $ basisValue b | b <- cb]+ where (Tagged cb) = completeBasis :: Tagged (Stiefel1Needle v) [Stiefel1Basis v] -data S⁰ = PositiveHalfSphere | NegativeHalfSphere deriving(Eq, Show)-newtype S¹ = S¹ { φParamS¹ :: Double -- [-π, π[- } deriving (Show)-data S² = S² { ϑParamS² :: !Double -- [0, π[- , φParamS² :: !Double -- [-π, π[- } deriving (Show)+instance (WithField k LinearManifold v, Real k) => Semimanifold (Stiefel1 v) where + type Needle (Stiefel1 v) = Stiefel1Needle v+ 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+ ιmν = 1-ν + in insi ιmν m+ | otherwise -> let m = HMat.scale ιmν spro + HMat.scale ((abs ιmν-1)/ν') n+ ιmν = ν-3+ in insi ιmν m+ where d = HMat.size s'+ s'= asPackedVector s+ ν' = l2norm n+ quop = signum s'i / ν'+ ν = ν' `mod'` 4+ im = HMat.maxIndex $ HMat.cmap abs s'+ s'i = s' HMat.! im+ spro = let v = deli s' in HMat.scale (recip s'i) v+ deli v = Arr.take im v Arr.++ Arr.drop (im+1) v+ insi ti v = Arr.generate d $ \i -> if | i<im -> v Arr.! i+ | i>im -> v Arr.! (i-1) + | otherwise -> ti+instance (WithField k LinearManifold v, Real k) => PseudoAffine (Stiefel1 v) where + Stiefel1 s .-~. Stiefel1 t = pure . Stiefel1Needle $ case s' HMat.! im of+ 0 -> HMat.scale (recip $ l2norm delis) delis+ 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+ in HMat.scale μ v+ | t'i/s'i > 0 -> samePoint+ | otherwise -> antipode+ where d = HMat.size t'+ s'= asPackedVector s; t' = asPackedVector t+ im = HMat.maxIndex $ HMat.cmap abs t'+ t'i = t' HMat.! im+ tpro = let v = deli t' in HMat.scale (recip t'i) v+ delis = deli s'+ deli v = Arr.take im v Arr.++ Arr.drop (im+1) v+ 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 -class NaturallyEmbedded m v where- embed :: m -> v- coEmbed :: v -> m++stiefel1Project :: LinearManifold v =>+ DualSpace v -- ^ Must be nonzero.+ -> Stiefel1 v+stiefel1Project = Stiefel1++stiefel1Embed :: HilbertSpace v => Stiefel1 v -> v+stiefel1Embed (Stiefel1 n) = normalized n -instance (VectorSpace y) => NaturallyEmbedded x (x,y) where- embed x = (x, zeroV)- coEmbed (x,_) = x-instance (VectorSpace y, VectorSpace z) => NaturallyEmbedded x ((x,y),z) where- embed x = (embed x, zeroV)- coEmbed (x,_) = coEmbed x+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 NaturallyEmbedded S⁰ ℝ where- embed PositiveHalfSphere = 1- embed NegativeHalfSphere = -1- coEmbed x | x>=0 = PositiveHalfSphere- | otherwise = NegativeHalfSphere-instance NaturallyEmbedded S¹ ℝ² where- embed (S¹ φ) = (cos φ, sin φ)- coEmbed (x,y) = S¹ $ atan2 y x-instance NaturallyEmbedded S² ℝ³ where- embed (S² ϑ φ) = ((cos φ * sin ϑ, sin φ * sin ϑ), cos ϑ)- coEmbed ((x,y),z) = S² (acos $ z/r) (atan2 y x)- where r = sqrt $ x^2 + y^2 + z^2- +instance HasUnitSphere ℝ where type UnitSphere ℝ = S⁰+instance HasUnitSphere ℝ² where type UnitSphere ℝ² = S¹+instance HasUnitSphere ℝ³ where type UnitSphere ℝ³ = S² +instance (HasUnitSphere v, v ~ DualSpace v) => NaturallyEmbedded (Stiefel1 v) v where+ embed = embed . unstiefel+ coEmbed = stiefel . coEmbed -type Endomorphism a = a->a -type ℝ = Double-type ℝ² = (ℝ,ℝ)-type ℝ³ = (ℝ²,ℝ) -instance VectorSpace () where- type Scalar () = ℝ- _ *^ () = ()+-- | Oriented hyperplanes, naïvely generalised to 'PseudoAffine' manifolds:+-- @'Cutplane' p w@ represents the set of all points 'q' such that+-- @(q.-~.p) ^\<.\> w ≡ 0@.+-- +-- In vector spaces this is indeed a hyperplane; for general manifolds it should+-- behave locally as a plane, globally as an (/n/−1)-dimensional submanifold.+data Cutplane x = Cutplane { sawHandle :: x+ , cutNormal :: Stiefel1 (Needle x) } -instance HasBasis () where- type Basis () = Void- basisValue = absurd- decompose () = []- decompose' () = absurd-instance InnerSpace () where- () <.> () = 0 +sideOfCut :: WithField ℝ Manifold x => Cutplane x -> x -> Option S⁰+sideOfCut (Cutplane sh (Stiefel1 cn)) p = decideSide . (cn<.>^) =<< p .-~. sh+ where decideSide 0 = mzero+ decideSide μ | μ > 0 = pure PositiveHalfSphere+ | otherwise = pure NegativeHalfSphere -(^) :: Num a => a -> Int -> a-(^) = (Prelude.^)++fathomCutDistance :: WithField ℝ Manifold x+ => Cutplane x -- ^ Hyperplane to measure the distance from.+ -> HerMetric'(Needle x) -- ^ Metric to use for measuring that distance.+ -- This can only be accurate if the metric+ -- is valid both around the cut-plane's 'sawHandle', and+ -- around the points you measure.+ -- (Strictly speaking, we would need /parallel transport/+ -- to ensure this).+ -> x -- ^ Point to measure the distance to.+ -> Option ℝ -- ^ A signed number, giving the distance from plane+ -- to point with indication on which side the point lies.+ -- 'Nothing' if the point isn't reachable from the plane.+fathomCutDistance (Cutplane sh (Stiefel1 cn)) met = \x -> fmap fathom $ x .-~. sh+ where fathom v = (cn <.>^ v) / scaleDist+ scaleDist = metric' met cn+
+ Data/Manifold/Types/Primitive.hs view
@@ -0,0 +1,253 @@+-- |+-- Module : Data.Manifold.Types.Primitive+-- Copyright : (c) Justus Sagemüller 2015+-- License : GPL v3+-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de+-- Stability : experimental+-- Portability : portable+-- +-- Several low-dimensional manifolds, represented in some simple way as Haskell+-- data types. All these are in the 'PseudoAffine' class.+-- +-- Also included in this module are some misc helper constraints etc., which don't really+-- belong here.+++{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE UndecidableInstances #-}+-- {-# LANGUAGE OverlappingInstances #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE RecordWildCards #-}+++module Data.Manifold.Types.Primitive (+ -- * Index / ASCII names+ Real0, Real1, RealPlus, Real2, Real3+ , Sphere0, Sphere1, Sphere2+ , Projective1, Projective2+ , Disk1, Disk2, Cone, OpenCone+ -- * Linear manifolds+ , ZeroDim(..)+ , ℝ⁰, ℝ, ℝ², ℝ³+ -- * Hyperspheres+ , S⁰(..), S¹(..), S²(..)+ -- * Projective spaces+ , ℝP¹, ℝP²(..)+ -- * Intervals\/disks\/cones+ , D¹(..), D²(..)+ , ℝay+ , CD¹(..), Cℝay(..)+ -- * Utility (deprecated)+ , NaturallyEmbedded(..)+ , GraphWindowSpec(..), Endomorphism, (^), EqFloating+ ) where+++import Data.VectorSpace+import Data.AffineSpace+import Data.Basis+import Data.Complex hiding (magnitude)+import Data.Void+import Data.Monoid++import qualified Prelude++import Control.Category.Constrained.Prelude hiding ((^))+import Control.Arrow.Constrained+import Control.Monad.Constrained+import Data.Foldable.Constrained+++++++type EqFloating f = (Eq f, Ord f, Floating f)++++data GraphWindowSpec = GraphWindowSpec {+ lBound, rBound, bBound, tBound :: Double+ , xResolution, yResolution :: Int+ }++++-- | A single point. Can be considered a zero-dimensional vector space, WRT any scalar.+data ZeroDim k = Origin deriving(Eq, Show)+instance Monoid (ZeroDim k) where+ mempty = Origin+ mappend Origin Origin = Origin+instance AdditiveGroup (ZeroDim k) where+ zeroV = Origin+ Origin ^+^ Origin = Origin+ negateV Origin = Origin+instance VectorSpace (ZeroDim k) where+ type Scalar (ZeroDim k) = k+ _ *^ Origin = Origin+instance HasBasis (ZeroDim k) where+ type Basis (ZeroDim k) = Void+ basisValue = absurd+ decompose Origin = []+ decompose' Origin = absurd++-- | 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.+data S⁰ = PositiveHalfSphere | NegativeHalfSphere deriving(Eq, Show)+-- | The unit circle.+newtype S¹ = S¹ { φParamS¹ :: Double -- ^ Must be in range @[-π, π[@.+ } deriving (Show)+-- | The ordinary unit sphere.+data S² = S² { ϑParamS² :: !Double -- ^ Range @[0, π[@.+ , φParamS² :: !Double -- ^ Range @[-π, π[@.+ } deriving (Show)+++++type ℝP¹ = S¹++-- | The two-dimensional real projective space, implemented as a unit disk with+-- opposing points on the rim glued together.+data ℝP² = ℝP² { rParamℝP² :: !Double -- ^ Range @[0, 1]@.+ , φParamℝP² :: !Double -- ^ Range @[-π, π[@.+ } deriving (Show)++++-- | The “one-dimensional disk” – really just the line segment between+-- the two points -1 and 1 of 'S⁰', i.e. this is simply a closed interval.+newtype D¹ = D¹ { xParamD¹ :: Double -- ^ Range @[-1, 1]@.+ }++-- | The standard, closed unit disk. Homeomorphic to the cone over 'S¹', but not in the+-- the obvious, “flat” way. (And not at all, despite+-- the identical ADT definition, to the projective space 'ℝP²'!)+data D² = D² { rParamD² :: !Double -- ^ Range @[0, 1]@.+ , φParamD² :: !Double -- ^ Range @[-π, π[@.+ } deriving (Show)++-- | A (closed) cone over a space @x@ is the product of @x@ with the closed interval 'D¹'+-- of “heights”,+-- except on its “tip”: here, @x@ is smashed to a single point.+-- +-- This construct becomes (homeomorphic-to-) an actual geometric cone (and to 'D²') in the+-- special case @x = 'S¹'@.+data CD¹ x = CD¹ { hParamCD¹ :: !Double -- ^ Range @[0, 1]@+ , pParamCD¹ :: !x -- ^ Irrelevant at @h = 0@.+ }+++-- | An open cone is homeomorphic to a closed cone without the “lid”,+-- i.e. without the “last copy” of @x@, at the far end of the height+-- interval. Since that means the height does not include its supremum, it is actually+-- more natural to express it as the entire real ray, hence the name.+data Cℝay x = Cℝay { hParamCℝay :: !Double -- ^ Range @[0, ∞[@+ , pParamCℝay :: !x -- ^ Irrelevant at @h = 0@.+ }++class NaturallyEmbedded m v where+ embed :: m -> v+ coEmbed :: v -> m+ ++instance (VectorSpace y) => NaturallyEmbedded x (x,y) where+ embed x = (x, zeroV)+ coEmbed (x,_) = x+instance (VectorSpace y, VectorSpace z) => NaturallyEmbedded x ((x,y),z) where+ embed x = (embed x, zeroV)+ coEmbed (x,_) = coEmbed x++instance NaturallyEmbedded S⁰ ℝ where+ embed PositiveHalfSphere = 1+ embed NegativeHalfSphere = -1+ coEmbed x | x>=0 = PositiveHalfSphere+ | otherwise = NegativeHalfSphere+instance NaturallyEmbedded S¹ ℝ² where+ embed (S¹ φ) = (cos φ, sin φ)+ coEmbed (x,y) = S¹ $ atan2 y x+instance NaturallyEmbedded S² ℝ³ where+ embed (S² ϑ φ) = ((cos φ * sin ϑ, sin φ * sin ϑ), cos ϑ)+ coEmbed ((x,y),z) = S² (acos $ z/r) (atan2 y x)+ where r = sqrt $ x^2 + y^2 + z^2+ +instance NaturallyEmbedded ℝP² ℝ³ where+ embed (ℝP² r φ) = ((r * cos φ, r * sin φ), sqrt $ 1-r^2)+ coEmbed ((x,y),z) = ℝP² (sqrt $ 1-(z/r)^2) (atan2 (y/r) (x/r))+ where r = sqrt $ x^2 + y^2 + z^2++instance NaturallyEmbedded D¹ ℝ where+ embed = xParamD¹+ coEmbed = D¹ . max (-1) . min 1++instance (NaturallyEmbedded x p) => NaturallyEmbedded (Cℝay x) (p,ℝ) where+ embed (Cℝay h p) = (embed p, h)+ coEmbed (v,z) = Cℝay (max 0 z) (coEmbed v)++++type Endomorphism a = a->a+++type ℝ⁰ = ZeroDim ℝ+type ℝ = Double+type ℝ² = (ℝ,ℝ)+type ℝ³ = (ℝ²,ℝ)+++-- | Better known as ℝ⁺ (which is not a legal Haskell name), the ray+-- of positive numbers (including zero, i.e. closed on one end).+type ℝay = Cℝay ℝ⁰+++++type Real0 = ℝ⁰+type Real1 = ℝ+type RealPlus = ℝay+type Real2 = ℝ²+type Real3 = ℝ³++type Sphere0 = S⁰+type Sphere1 = S¹+type Sphere2 = S²++type Projective1 = ℝP¹+type Projective2 = ℝP²++type Disk1 = D¹+type Disk2 = D²++type Cone = CD¹ +type OpenCone = Cℝay++++instance VectorSpace () where+ type Scalar () = ℝ+ _ *^ () = ()++instance HasBasis () where+ type Basis () = Void+ basisValue = absurd+ decompose () = []+ decompose' () = absurd+instance InnerSpace () where+ () <.> () = 0++++(^) :: Num a => a -> Int -> a+(^) = (Prelude.^)+
+ Data/SimplicialComplex.hs view
@@ -0,0 +1,500 @@+-- |+-- Module : Data.SimplicialComplex+-- 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 TypeOperators #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DataKinds #-}+++module Data.SimplicialComplex (+ -- * Simplices+ Simplex(..)+ -- ** Construction+ , (.<.), makeSimplex, makeSimplex'+ -- ** Deconstruction+ , simplexVertices, simplexVertices'+ -- * Simplicial complexes+ , Triangulation+ , singleSimplex+ -- * Triangulation-builder monad+ , TriangT+ , evalTriangT, runTriangT, doTriangT, getTriang+ -- ** Subsimplex-references+ , SimplexIT, simplexITList, lookSimplex+ , lookSplxFacesIT, lookSupersimplicesIT, tgetSimplexIT+ , lookVertexIT, lookSplxVerticesIT+ , sharedBoundary+ , distinctSimplices, NeighbouringSimplices+ -- ** Building triangulations+ , disjointTriangulation+ , disjointSimplex+ , mixinTriangulation+ , introVertToTriang+ , webinateTriang+ -- * Misc util+ , HaskMonad, liftInTriangT, unliftInTriangT+ , Nat, Zero, One, Two, Three, Succ+ ) where++++import Data.List hiding (filter, all, elem)+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 Data.VectorSpace+import Data.LinearMap+import Data.LinearMap.Category+import Data.Void+import Data.Tagged+import Data.Proxy++import Data.Manifold.Types+import Data.Manifold.Types.Primitive ((^))+import Data.Manifold.PseudoAffine+ +import Data.Embedding+import Data.CoNat++import qualified Prelude as Hask hiding(foldl)+import qualified Control.Applicative as Hask+import qualified Control.Monad as Hask+import Control.Monad.Trans.List+import Control.Monad.Trans.Class+import qualified Data.Foldable as Hask+import Data.Foldable (all, elem)++import Data.Functor.Identity (Identity, runIdentity)++import Control.Category.Constrained.Prelude hiding ((^), all, elem)+import Control.Arrow.Constrained+import Control.Monad.Constrained+import Data.Foldable.Constrained++import GHC.Generics (Generic)++infixr 5 :<|, .<.++-- | An /n/-simplex is a connection of /n/+1 points in a simply connected region of a manifold.+data Simplex :: Nat -> * -> * where+ ZS :: !x -> Simplex Z x+ (:<|) :: KnownNat n => !x -> !(Simplex n x) -> Simplex (S n) x++deriving instance (Show x) => Show (Simplex n x)+instance Hask.Functor (Simplex n) where+ fmap f (ZS x) = ZS (f x)+ fmap f (x:<|xs) = f x :<| fmap f xs++-- | Use this together with ':<|' to easily build simplices, like you might construct lists.+-- E.g. @(0,0) ':<|' (1,0) '.<.' (0,1) :: 'Simplex' 'Two' ℝ²@.+(.<.) :: x -> x -> Simplex One x+x .<. y = x :<| ZS y+++makeSimplex :: ∀ x n . KnownNat n => x ^ S n -> Simplex n x+makeSimplex xs = case makeSimplex' $ Hask.toList xs of+ Option (Just s) -> s++makeSimplex' :: ∀ x n . KnownNat n => [x] -> Option (Simplex n x)+makeSimplex' [] = Option Nothing+makeSimplex' [x] = cozeroT $ ZS x+makeSimplex' (x:xs) = fCosuccT ((x:<|) <$> makeSimplex' xs)++simplexVertices :: ∀ x n . Simplex n x -> x ^ S n+simplexVertices (ZS x) = pure x+simplexVertices (x :<| s) = freeCons x (simplexVertices s)++simplexVertices' :: ∀ x n . Simplex n x -> [x]+simplexVertices' (ZS x) = [x]+simplexVertices' (x :<| s) = x : simplexVertices' s+++type Array = Arr.Vector++-- | An /n/-dimensional /abstract simplicial complex/ is a collection of /n/-simplices+-- which are “glued together” in some way. The preferred way to construct+-- such complexes is to run a 'TriangT' builder.+data Triangulation (n :: Nat) (x :: *) where+ TriangSkeleton :: KnownNat n+ => Triangulation n x -- The lower-dimensional skeleton.+ -> Array -- Array of `S n`-simplices in this triangulation.+ ( Int ^ S (S n) -- “down link” – the subsimplices+ , [Int] -- “up link” – what higher simplices have+ ) -- this one as a subsimplex?+ -> Triangulation (S n) x+ TriangVertices :: Array (x, [Int]) -> Triangulation Z x+instance Hask.Functor (Triangulation n) where+ fmap f (TriangVertices vs) = TriangVertices $ first f <$> vs+ fmap f (TriangSkeleton sk vs) = TriangSkeleton (f<$>sk) vs+deriving instance (Show x) => Show (Triangulation n x)++-- | Consider a single simplex as a simplicial complex, consisting only of+-- this simplex and its faces.+singleSimplex :: ∀ n x . KnownNat n => Simplex n x -> Triangulation n x+singleSimplex (ZS x) = TriangVertices $ pure (x, [])+singleSimplex (x :<| s)+ = runIdentity . execTriangT insX $ TriangSkeleton (singleSimplex s) mempty+ where insX :: ∀ t . TriangT t n x Identity ()+ insX = introVertToTriang x [SimplexIT 0] >> return()++nTopSplxs :: Triangulation n' x -> Int+nTopSplxs (TriangVertices vs) = Arr.length vs+nTopSplxs (TriangSkeleton _ vs) = Arr.length vs++nSplxs :: ∀ k n x . (KnownNat k, KnownNat n) => Triangulation n x -> Tagged k Int+nSplxs t = case t of+ TriangVertices vs | n == k -> Tagged $ Arr.length vs+ TriangSkeleton _ vs | n == k -> Tagged $ Arr.length vs+ TriangSkeleton sk _ | n > k -> nSplxs sk+ _ -> Tagged 0+ where (Tagged k) = theNatN :: Tagged k Int+ (Tagged n) = theNatN :: Tagged n Int++-- | Combine two triangulations (assumed as disjoint) to a single, non-connected complex.+instance (KnownNat n) => Semigroup (Triangulation n x) where+ TriangVertices vs₁ <> TriangVertices vs₂ = TriangVertices $ vs₁ Arr.++ vs₂+ TriangSkeleton sk₁ sp₁ <> TriangSkeleton sk₂ sp₂+ = TriangSkeleton (sk₁ <> shiftUprefs (Arr.length sp₁) sk₂)+ (sp₁ Arr.++ fmap (first $ fmap (+ nTopSplxs sk₁)) sp₂)+ where shiftUprefs :: Int -> Triangulation n' x -> Triangulation n' x+ shiftUprefs δn (TriangVertices vs)+ = TriangVertices $ fmap (second $ fmap (+δn)) vs+ shiftUprefs δn (TriangSkeleton sk' vs)+ = TriangSkeleton sk' $ fmap (second $ fmap (+δn)) vs+instance (KnownNat n) => Monoid (Triangulation n x) where+ mappend = (<>)+ mempty = coInduceT (TriangVertices mempty) (`TriangSkeleton`mempty)+++++ +-- | A “conservative” state monad containing a 'Triangulation'. It+-- can be extended by new simplices, which can then be indexed using 'SimplexIT'.+-- The universally-quantified @t@ argument ensures you can't index simplices that+-- don't actually exist in this triangulation.+newtype TriangT t n x m y = TriangT {+ unsafeRunTriangT :: Triangulation n x -> m (y, Triangulation n x) }+ deriving (Hask.Functor)+instance (Hask.Functor m, Monad m (->))+ => Hask.Applicative (TriangT t n x m) where+ pure x = TriangT $ pure . (x,)+ TriangT fs <*> TriangT xs = TriangT $+ fs >=> \(f, t') -> fmap (first f) $ xs t'+instance (Hask.Functor m, Monad m (->)) => Hask.Monad (TriangT t n x m) where+ return x = TriangT $ pure . (x,)+ TriangT xs >>= f = TriangT $+ \t -> xs t >>= \(y,t') -> let (TriangT zs) = f y in zs t'++instance MonadTrans (TriangT t n x) where+ lift m = TriangT $ \tr -> Hask.liftM (,tr) m++type HaskMonad m = (Hask.Applicative m, Hask.Monad m)++triangReadT :: ∀ t n x m y . HaskMonad m => (Triangulation n x -> m y) -> TriangT t n x m y+triangReadT f = TriangT $ \t -> fmap (,t) $ f t++unsafeEvalTriangT :: ∀ n t x m y . HaskMonad m+ => TriangT t n x m y -> Triangulation n x -> m y+unsafeEvalTriangT t = fmap fst . unsafeRunTriangT t++execTriangT :: ∀ n x m y . HaskMonad m => (∀ t . TriangT t n x m y)+ -> Triangulation n x -> m (Triangulation n x)+execTriangT t = fmap snd . unsafeRunTriangT (t :: TriangT () n x m y)++evalTriangT :: ∀ n x m y . (KnownNat n, HaskMonad m) => (∀ t . TriangT t n x m y) -> m y+evalTriangT t = fmap fst (unsafeRunTriangT (t :: TriangT () n x m y) mempty)++runTriangT :: ∀ n x m y . (∀ t . TriangT t n x m y)+ -> Triangulation n x -> m (y, Triangulation n x)+runTriangT t = unsafeRunTriangT (t :: TriangT () n x m y)++doTriangT :: ∀ n x m y . KnownNat n => (∀ t . TriangT t n x m y) -> m (y, Triangulation n x)+doTriangT t = runTriangT t mempty++getEntireTriang :: ∀ t n x m . HaskMonad m => TriangT t n x m (Triangulation n x)+getEntireTriang = TriangT $ \t -> pure (t, t)++getTriang :: ∀ t n k x m . (HaskMonad m, KnownNat k, KnownNat n)+ => TriangT t n x m (Option (Triangulation k x))+getTriang = onSkeleton getEntireTriang++liftInTriangT :: ∀ t n x m μ y . (HaskMonad m, MonadTrans μ)+ => TriangT t n x m y -> TriangT t n x (μ m) y+liftInTriangT (TriangT b) = TriangT $ lift . b++unliftInTriangT :: ∀ t n x m μ y . (HaskMonad m, MonadTrans μ)+ => (∀ m' a . μ m a -> m a) -> TriangT t n x (μ m) y -> TriangT t n x m y+unliftInTriangT unlift (TriangT b) = TriangT $ \t -> unlift (b t)++++forgetVolumes :: ∀ n x t m y . (KnownNat n, HaskMonad m)+ => TriangT t n x m y -> TriangT t (S n) x m y+forgetVolumes (TriangT f) = TriangT $ \(TriangSkeleton l bk)+ -> fmap (\(y, l') -> (y, TriangSkeleton l' bk)) $ f l++onSkeleton :: ∀ n k x t m y . (KnownNat k, KnownNat n, HaskMonad m)+ => TriangT t k x m y -> TriangT t n x m (Option y)+onSkeleton q@(TriangT qf) = case tryToMatchTTT forgetVolumes q of+ Option (Just q') -> pure <$> q'+ _ -> return Hask.empty+++newtype SimplexIT (t :: *) (n :: Nat) (x :: *) = SimplexIT { tgetSimplexIT' :: Int }+ deriving (Eq, Ord, Show)++-- | A unique (for the given dimension) ID of a triagulation's simplex. It is the index+-- where that simplex can be found in the 'simplexITList'.+tgetSimplexIT :: SimplexIT t n x -> Int+tgetSimplexIT = tgetSimplexIT'++-- | Reference the /k/-faces of a given simplex in a triangulation.+lookSplxFacesIT :: ∀ t m n k x . (HaskMonad m, KnownNat k, KnownNat n)+ => SimplexIT t (S k) x -> TriangT t n x m (SimplexIT t k x ^ S(S k))+lookSplxFacesIT = fmap (\(Option(Just r))->r) . onSkeleton . lookSplxFacesIT'++lookSplxFacesIT' :: ∀ t m n x . (HaskMonad m, KnownNat n)+ => SimplexIT t (S n) x -> TriangT t (S n) x m (SimplexIT t n x ^ S(S n))+lookSplxFacesIT' (SimplexIT i) = triangReadT rc+ where rc (TriangSkeleton _ ssb) = return . fmap SimplexIT . fst $ ssb Arr.! i++lookSplxVerticesIT :: ∀ t m n k x . (HaskMonad m, KnownNat k, KnownNat n)+ => SimplexIT t k x -> TriangT t n x m (SimplexIT t Z x ^ S k)+lookSplxVerticesIT = fmap (\(Option(Just r))->r) . onSkeleton . lookSplxVerticesIT'++lookSplxVerticesIT' :: ∀ t m n x . (HaskMonad m, KnownNat n)+ => SimplexIT t n x -> TriangT t n x m (SimplexIT t Z x ^ S n)+lookSplxVerticesIT' i = fmap + (\vs -> case freeVector vs of+ Option (Just vs') -> vs'+ _ -> error $ "Impossible number " ++ show (length vs) ++ " of vertices for "+ ++ show n ++ "-simplex in `lookSplxVerticesIT'`."+ ) $ lookSplxsVerticesIT [i]+ where (Tagged n) = theNatN :: Tagged n Int+ ++lookSplxsVerticesIT :: ∀ t m n x . HaskMonad m+ => [SimplexIT t n x] -> TriangT t n x m [SimplexIT t Z x]+lookSplxsVerticesIT is = triangReadT rc+ where rc (TriangVertices _) = return is+ rc (TriangSkeleton sk up) = unsafeEvalTriangT+ ( lookSplxsVerticesIT+ $ SimplexIT <$> fastNub [ j | SimplexIT i <- is+ , j <- Hask.toList . fst $ up Arr.! i ]+ ) sk++lookVertexIT :: ∀ t m n x . (HaskMonad m, KnownNat n)+ => SimplexIT t Z x -> TriangT t n x m x+lookVertexIT = fmap (\(Option(Just r))->r) . onSkeleton . lookVertexIT'++lookVertexIT' :: ∀ t m x . HaskMonad m => SimplexIT t Z x -> TriangT t Z x m x+lookVertexIT' (SimplexIT i) = triangReadT $ \(TriangVertices vs) -> return.fst $ vs Arr.! i++lookSimplex :: ∀ t m n k x . (HaskMonad m, KnownNat k, KnownNat n)+ => SimplexIT t k x -> TriangT t n x m (Simplex k x)+lookSimplex s = do + vis <- lookSplxVerticesIT s+ fmap makeSimplex $ mapM lookVertexIT vis++simplexITList :: ∀ t m n k x . (HaskMonad m, KnownNat k, KnownNat n)+ => TriangT t n x m [SimplexIT t k x]+simplexITList = fmap (\(Option(Just r))->r) $ onSkeleton simplexITList'++simplexITList' :: ∀ t m n x . (HaskMonad m, KnownNat n)+ => TriangT t n x m [SimplexIT t n x]+simplexITList' = triangReadT $ return . sil+ where sil :: Triangulation n x -> [SimplexIT t n x]+ sil (TriangVertices vs) = [ SimplexIT i | i <- [0 .. Arr.length vs - 1] ]+ sil (TriangSkeleton _ bk) = [ SimplexIT i | i <- [0 .. Arr.length bk - 1] ]+++lookSupersimplicesIT :: ∀ t m n k j x . (HaskMonad m, KnownNat k, KnownNat j, KnownNat n)+ => SimplexIT t k x -> TriangT t n x m [SimplexIT t j x]+lookSupersimplicesIT = runListT . defLstt . matchLevel . pure+ where lvlIt :: ∀ i . (KnownNat i, KnownNat n) => ListT (TriangT t n x m) (SimplexIT t i x)+ -> ListT (TriangT t n x m) (SimplexIT t (S i) x)+ lvlIt (ListT m) = ListT . fmap (fnubConcatBy $ comparing tgetSimplexIT)+ $ mapM lookSupersimplicesIT' =<< m+ matchLevel = ftorTryToMatchT lvlIt+ defLstt (Option (Just lt)) = lt+ defLstt _ = ListT $ return []++lookSupersimplicesIT' :: ∀ t m n k x . (HaskMonad m, KnownNat k, KnownNat n)+ => SimplexIT t k x -> TriangT t n x m [SimplexIT t (S k) x]+lookSupersimplicesIT' = fmap (\(Option(Just r))->r) . onSkeleton . lookSupersimplicesIT''++lookSupersimplicesIT'' :: ∀ t m n x . (HaskMonad m, KnownNat n)+ => SimplexIT t n x -> TriangT t (S n) x m [SimplexIT t (S n) x]+lookSupersimplicesIT'' (SimplexIT i) =+ fmap ( \tr -> SimplexIT <$> case tr of+ TriangSkeleton (TriangSkeleton _ tsps) _ -> snd (tsps Arr.! i)+ TriangSkeleton (TriangVertices tsps) _ -> snd (tsps Arr.! i)+ ) getEntireTriang++sharedBoundary :: ∀ t m n k x . (HaskMonad m, KnownNat k, KnownNat n)+ => SimplexIT t (S k) x -> SimplexIT t (S k) x+ -> TriangT t n x m (Option (SimplexIT t k x))+sharedBoundary i j = fmap snd <$> distinctSimplices i j++type NeighbouringSimplices t n x = ((SimplexIT t Z x, SimplexIT t Z x), SimplexIT t n x)++distinctSimplices :: ∀ t m n k x . (HaskMonad m, KnownNat k, KnownNat n)+ => SimplexIT t (S k) x -> SimplexIT t (S k) x+ -> TriangT t n x m (Option (NeighbouringSimplices t k x))+distinctSimplices i j = do+ [iSubs,jSubs] <- mapM lookSplxFacesIT [i,j]+ case fnubIntersect (Hask.toList iSubs) (Hask.toList jSubs) of+ [shBound] -> do+ shVerts <- lookSplxVerticesIT shBound+ [[iIVert], [jIVert]] <- forM [i,j]+ $ fmap (filter (not . (`elem` shVerts)) . Hask.toList) . lookSplxVerticesIT+ return $ pure ((iIVert, jIVert), shBound)+ _ -> return Hask.empty+++triangulationBulk :: ∀ t m n k x . (HaskMonad m, KnownNat k, KnownNat n) => TriangT t n x m [Simplex k x]+triangulationBulk = simplexITList >>= mapM lookSimplex++withThisSubsimplex :: ∀ t m n k j x . (HaskMonad m, KnownNat j, KnownNat k, KnownNat n)+ => SimplexIT t j x -> TriangT t n x m [SimplexIT t k x]+withThisSubsimplex s = do+ svs <- lookSplxVerticesIT s+ simplexITList >>= filterM (lookSplxVerticesIT >>> fmap`id`+ \s'vs -> all (`elem`s'vs) svs )++lookupSimplexCone :: ∀ t m n k x . ( HaskMonad m, KnownNat k, KnownNat n )+ => SimplexIT t Z x -> SimplexIT t k x -> TriangT t n x m (Option (SimplexIT t (S k) x))+lookupSimplexCone tip base = do+ tipSups :: [SimplexIT t (S k) x] <- lookSupersimplicesIT tip+ baseSups :: [SimplexIT t (S k) x] <- lookSupersimplicesIT base+ return $ case intersect tipSups baseSups of+ (res:_) -> pure res+ _ -> Hask.empty+ +++-- | Import an entire triangulation, as disjoint from everything already in the monad.+disjointTriangulation :: ∀ t m n x . (KnownNat n, HaskMonad m)+ => Triangulation n x -> TriangT t n x m [SimplexIT t n x]+disjointTriangulation t = TriangT $+ \tr -> return ( [ SimplexIT k+ | k <- take (nTopSplxs t) [nTopSplxs tr ..] ]+ , tr <> t )++disjointSimplex :: ∀ t m n x . (KnownNat n, HaskMonad m)+ => Simplex n x -> TriangT t n x m (SimplexIT t n x)+disjointSimplex s = TriangT $ \tr -> return ( SimplexIT $ nTopSplxs tr+ , tr <> singleSimplex s )+++-- | Import a triangulation like with 'disjointTriangulation',+-- together with references to some of its subsimplices.+mixinTriangulation :: ∀ t m f k n x . ( KnownNat n, KnownNat k+ , HaskMonad m, Functor f (->) (->) )+ => (∀ s . TriangT s n x m (f (SimplexIT s k x)))+ -> TriangT t n x m (f (SimplexIT t k x))+mixinTriangulation t+ = TriangT $ \tr -> do+ (sqs, tr') <- doTriangT t'+ let (Tagged n) = nSplxs tr :: Tagged k Int+ return ( fmap (\k -> SimplexIT $ n + k) sqs, tr <> tr' )+ where t' :: ∀ s . TriangT s n x m (f Int)+ t' = fmap (fmap tgetSimplexIT) t+++webinateTriang :: ∀ t m n x . (HaskMonad m, KnownNat n)+ => SimplexIT t Z x -> SimplexIT t n x -> TriangT t (S n) x m (SimplexIT t (S n) x)+webinateTriang ptt@(SimplexIT pt) bst@(SimplexIT bs) = do+ existsReady <- lookupSimplexCone ptt bst+ case existsReady of+ Option (Just ext) -> return ext+ _ -> TriangT $ \(TriangSkeleton sk cnn)+ -> let resi = Arr.length cnn+ res = SimplexIT $ Arr.length cnn :: SimplexIT t (S n) x+ in case sk of+ TriangVertices vs -> return+ $ ( res+ , TriangSkeleton (TriangVertices+ $ vs Arr.// [ (pt, second (resi:) $ vs Arr.! pt)+ , (bs, second (resi:) $ vs Arr.! bs) ]+ ) $ Arr.snoc cnn (freeTuple$->$(pt, bs), []) )+ TriangSkeleton _ cnn'+ -> let (cnbs,_) = cnn' Arr.! bs+ in do (cnws,sk') <- unsafeRunTriangT ( do+ cnws <- forM cnbs $ \j -> do+ kt@(SimplexIT k) <- webinateTriang ptt (SimplexIT j)+ addUplink' res kt+ return k+ addUplink' res bst+ return cnws+ ) sk+ let snocer = (freeSnoc cnws bs, [])+ return $ (res, TriangSkeleton sk' $ Arr.snoc cnn snocer)+ where addUplink' :: SimplexIT t (S n) x -> SimplexIT t n x -> TriangT t n x m ()+ addUplink' (SimplexIT i) (SimplexIT j) = TriangT+ $ \sk -> pure ((), case sk of+ TriangVertices vs+ -> let (v,ul) = vs Arr.! j+ in TriangVertices $ vs Arr.// [(j, (v, i:ul))]+ TriangSkeleton skd us+ -> let (b,tl) = us Arr.! j+ in TriangSkeleton skd $ us Arr.// [(j, (b, i:tl))]+ )+ ++++introVertToTriang :: ∀ t m n x . (HaskMonad m, KnownNat n)+ => x -> [SimplexIT t n x] -> TriangT t (S n) x m (SimplexIT t Z x)+introVertToTriang v glues = do+ j <- fmap (\(Option(Just k)) -> SimplexIT k) . onSkeleton . TriangT+ $ return . tVertSnoc+ mapM_ (webinateTriang j) glues+ return j+ where tVertSnoc :: Triangulation Z x -> (Int, Triangulation Z x)+ tVertSnoc (TriangVertices vs)+ = (Arr.length vs, TriangVertices $ vs `Arr.snoc` (v,[]))+ +++++-- | Type-level zero of kind 'Nat'.+type Zero = Z+type One = S Zero+type Two = S One+type Three = S Two+type Succ = S++
+ Data/VectorSpace/FiniteDimensional.hs view
@@ -0,0 +1,163 @@+-- |+-- Module : Data.VectorSpace.FiniteDimensional+-- Copyright : (c) Justus Sagemüller 2015+-- License : GPL v3+-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de+-- Stability : experimental+-- Portability : portable+-- +{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE ScopedTypeVariables #-}+++++module Data.VectorSpace.FiniteDimensional (+ FiniteDimensional(..)+ , SmoothScalar + ) where+ ++ ++import Prelude hiding ((^))++import Data.VectorSpace+import Data.LinearMap+import Data.Basis+import Data.MemoTrie+import Data.Tagged+import Data.Void++import Control.Applicative+ +import Data.Manifold.Types.Primitive+import Data.CoNat++import qualified Data.Vector as Arr+import qualified Numeric.LinearAlgebra.HMatrix as HMat+++++-- | Constraint that a space's scalars need to fulfill so it can be used for efficient linear algebra.+-- Fulfilled pretty much only by the basic real and complex floating-point types.+type SmoothScalar s = ( VectorSpace s, HMat.Numeric s, HMat.Field s+ , 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.+-- More importantly, @hmatrix@ offers linear facilities+-- such as inverse and eigenbasis transformations, which aren't available in the+-- @vector-space@ library yet. But the classes from that library are strongly preferrable+-- to plain matrices and arrays, conceptually.+-- +-- The 'FiniteDimensional' class is used to convert between both representations.+-- It would be nice not to have the requirement of finite dimension on 'HerMetric',+-- but it's probably not feasible to get rid of it in forseeable time.+-- +-- Instead of the run-time 'dimension' information, we would rather have a compile-time+-- @type Dimension v :: Nat@, but type-level naturals are not mature enough yet. This+-- will almost certainly change in the future.+class (HasBasis v, HasTrie (Basis v), SmoothScalar (Scalar v)) => FiniteDimensional v where+ dimension :: Tagged v Int+ basisIndex :: Tagged v (Basis v -> Int)+ -- | Index must be in @[0 .. dimension-1]@, otherwise this is undefined.+ indexBasis :: Tagged v (Int -> Basis v)+ completeBasis :: Tagged v [Basis v]+ completeBasis = liftA2 (\dim f -> f <$> [0 .. dim - 1]) dimension indexBasis+ + asPackedVector :: v -> HMat.Vector (Scalar v)+ asPackedVector v = HMat.fromList $ snd <$> decompose v+ + asPackedMatrix :: (FiniteDimensional w, Scalar w ~ Scalar v)+ => (v :-* w) -> HMat.Matrix (Scalar v)+ asPackedMatrix = defaultAsPackedMatrix+ where defaultAsPackedMatrix :: forall v w s .+ (FiniteDimensional v, FiniteDimensional w, s~Scalar v, s~Scalar w)+ => (v :-* w) -> HMat.Matrix s+ defaultAsPackedMatrix m = HMat.fromRows $ asPackedVector . atBasis m <$> cb+ where (Tagged cb) = completeBasis :: Tagged v [Basis v]+ + fromPackedVector :: HMat.Vector (Scalar v) -> v+ fromPackedVector v = result+ where result = recompose $ zip cb (HMat.toList v)+ cb = witness completeBasis result++instance (SmoothScalar k) => FiniteDimensional (ZeroDim k) where+ dimension = Tagged 0+ basisIndex = Tagged absurd+ indexBasis = Tagged $ const undefined+ completeBasis = Tagged []+ asPackedVector Origin = HMat.fromList []+ fromPackedVector _ = Origin+instance FiniteDimensional ℝ where+ dimension = Tagged 1+ basisIndex = Tagged $ \() -> 0+ indexBasis = Tagged $ \0 -> ()+ completeBasis = Tagged [()]+ asPackedVector x = HMat.fromList [x]+ asPackedMatrix f = HMat.asRow . asPackedVector $ atBasis f ()+ fromPackedVector v = v HMat.! 0+instance (FiniteDimensional a, FiniteDimensional b, Scalar a~Scalar b)+ => FiniteDimensional (a,b) where+ dimension = tupDim+ where tupDim :: forall a b.(FiniteDimensional a,FiniteDimensional b)=>Tagged(a,b)Int+ tupDim = Tagged $ da+db+ where (Tagged da)=dimension::Tagged a Int; (Tagged db)=dimension::Tagged b Int+ basisIndex = basId+ where basId :: forall a b . (FiniteDimensional a, FiniteDimensional b)+ => Tagged (a,b) (Either (Basis a) (Basis b) -> Int)+ basId = Tagged basId'+ where basId' (Left ba) = basIda ba+ basId' (Right bb) = da + basIdb bb+ (Tagged da) = dimension :: Tagged a Int+ (Tagged basIda) = basisIndex :: Tagged a (Basis a->Int)+ (Tagged basIdb) = basisIndex :: Tagged b (Basis b->Int)+ indexBasis = basId+ where basId :: forall a b . (FiniteDimensional a, FiniteDimensional b)+ => Tagged (a,b) (Int -> Either (Basis a) (Basis b))+ basId = Tagged basId'+ where basId' i | i < da = Left $ basIda i+ | otherwise = Right . basIdb $ i - da+ (Tagged da) = dimension :: Tagged a Int+ (Tagged basIda) = indexBasis :: Tagged a (Int->Basis a)+ (Tagged basIdb) = indexBasis :: Tagged b (Int->Basis b)+ completeBasis = cb+ where cb :: forall a b . (FiniteDimensional a, FiniteDimensional b)+ => Tagged (a,b) [Either (Basis a) (Basis b)]+ cb = Tagged $ map Left cba ++ map Right cbb+ where (Tagged cba) = completeBasis :: Tagged a [Basis a]+ (Tagged cbb) = completeBasis :: Tagged b [Basis b]+ asPackedVector (a,b) = HMat.vjoin [asPackedVector a, asPackedVector b]+ fromPackedVector = fPV+ where fPV :: forall a b . (FiniteDimensional a, FiniteDimensional b, Scalar a~Scalar b)+ => HMat.Vector (Scalar a) -> (a,b)+ fPV v = (fromPackedVector l, fromPackedVector r)+ where (Tagged da) = dimension :: Tagged a Int+ (Tagged db) = dimension :: Tagged b Int+ l = HMat.subVector 0 da v+ r = HMat.subVector da db v+ + +instance (SmoothScalar x, KnownNat n) => FiniteDimensional (FreeVect n x) where+ dimension = natTagPænultimate+ basisIndex = Tagged getInRange+ indexBasis = Tagged InRange+ asPackedVector (FreeVect arr) = Arr.convert arr+ fromPackedVector arr = FreeVect (Arr.convert arr)+ -- asPackedMatrix = _ -- could be done quite efficiently here!+ +
+ images/examples/simple-2d-ShadeTree.png view
binary file changed (absent → 107507 bytes)
manifolds.cabal view
@@ -1,5 +1,5 @@ Name: manifolds-Version: 0.1.0.2+Version: 0.1.3.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,@@ -25,10 +25,16 @@ License: GPL-3 License-file: COPYING Author: Justus Sagemüller+Homepage: https://github.com/leftaroundabout/manifolds Maintainer: (@) sagemueller $ geo.uni-koeln.de Build-Type: Simple Cabal-Version: >=1.10+Extra-Doc-Files: images/examples/*.png +Source-Repository head+ type: git+ location: git://github.com/leftaroundabout/manifolds.git+ Library Build-Depends: base>=4.5 && < 6 , transformers@@ -36,13 +42,13 @@ , MemoTrie , vector , vector-algorithms+ , hmatrix >= 0.16 && < 0.18 , containers- , random- , MonadRandom , comonad , semigroups , void , tagged+ , deepseq , constrained-categories >= 0.2 && < 0.3 other-extensions: FlexibleInstances , TypeFamilies@@ -57,10 +63,17 @@ ghc-options: -O2 Exposed-modules: Data.Manifold Data.Manifold.PseudoAffine+ Data.Manifold.TreeCover+ Data.SimplicialComplex Data.LinearMap.HerMetric -- Data.Manifold.Visualisation.R3.GLUT- Other-modules: Data.Manifold.Types- Data.List.FastNub+ Data.Manifold.Types+ Other-modules: Data.List.FastNub+ Data.Manifold.Types.Primitive+ Data.CoNat+ Data.Embedding+ Data.LinearMap.Category+ Data.VectorSpace.FiniteDimensional Util.Associate Util.LtdShow default-language: Haskell2010