diff --git a/Data/CoNat.hs b/Data/CoNat.hs new file mode 100644 --- /dev/null +++ b/Data/CoNat.hs @@ -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) diff --git a/Data/Embedding.hs b/Data/Embedding.hs new file mode 100644 --- /dev/null +++ b/Data/Embedding.hs @@ -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 + + + + + + + diff --git a/Data/LinearMap/Category.hs b/Data/LinearMap/Category.hs new file mode 100644 --- /dev/null +++ b/Data/LinearMap/Category.hs @@ -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) + + diff --git a/Data/LinearMap/HerMetric.hs b/Data/LinearMap/HerMetric.hs --- a/Data/LinearMap/HerMetric.hs +++ b/Data/LinearMap/HerMetric.hs @@ -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' + diff --git a/Data/List/FastNub.hs b/Data/List/FastNub.hs --- a/Data/List/FastNub.hs +++ b/Data/List/FastNub.hs @@ -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) | xy = fis (x:xs) ys + | otherwise = x : fis xs ys diff --git a/Data/Manifold.hs b/Data/Manifold.hs --- a/Data/Manifold.hs +++ b/Data/Manifold.hs @@ -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 diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs --- a/Data/Manifold/PseudoAffine.hs +++ b/Data/Manifold/PseudoAffine.hs @@ -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. - + + + + + + diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/TreeCover.hs @@ -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 +-- 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] ] +-- @ +-- +-- <> +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 + diff --git a/Data/Manifold/Types.hs b/Data/Manifold/Types.hs --- a/Data/Manifold/Types.hs +++ b/Data/Manifold/Types.hs @@ -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 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 + diff --git a/Data/Manifold/Types/Primitive.hs b/Data/Manifold/Types/Primitive.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/Types/Primitive.hs @@ -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.^) + diff --git a/Data/SimplicialComplex.hs b/Data/SimplicialComplex.hs new file mode 100644 --- /dev/null +++ b/Data/SimplicialComplex.hs @@ -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 + + diff --git a/Data/VectorSpace/FiniteDimensional.hs b/Data/VectorSpace/FiniteDimensional.hs new file mode 100644 --- /dev/null +++ b/Data/VectorSpace/FiniteDimensional.hs @@ -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! + + diff --git a/images/examples/simple-2d-ShadeTree.png b/images/examples/simple-2d-ShadeTree.png new file mode 100644 Binary files /dev/null and b/images/examples/simple-2d-ShadeTree.png differ diff --git a/manifolds.cabal b/manifolds.cabal --- a/manifolds.cabal +++ b/manifolds.cabal @@ -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