manifolds 0.4.5.0 → 0.5.0.0
raw patch · 15 files changed
+929/−933 lines, 15 filesdep +arraydep +ieee754dep +spatial-rotationsdep ~manifolds-core
Dependencies added: array, ieee754, spatial-rotations
Dependency ranges changed: manifolds-core
Files
- Data/CoNat.hs +0/−325
- Data/Manifold/Cone.hs +0/−2
- Data/Manifold/FibreBundle.hs +34/−2
- Data/Manifold/Function/Interpolation.hs +14/−0
- Data/Manifold/Mesh.hs +62/−0
- Data/Manifold/PseudoAffine.hs +28/−4
- Data/Manifold/Riemannian.hs +0/−2
- Data/Manifold/TreeCover.hs +0/−122
- Data/Manifold/Types/Primitive.hs +22/−2
- Data/Manifold/Web.hs +32/−2
- Data/Simplex/Abstract.hs +106/−0
- Data/SimplicialComplex.hs +0/−423
- Math/Manifold/Real/Coordinates.hs +379/−0
- manifolds.cabal +9/−4
- test/tasty/test.hs +243/−45
− Data/CoNat.hs
@@ -1,325 +0,0 @@--- |--- 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 ExplicitNamespaces #-}-{-# LANGUAGE DataKinds #-}-{-# LANGUAGE PolyKinds #-}--module Data.CoNat ( Nat(..), natToInt, fromNat- , natTagLast, natTagPænultimate, natTagAntepænultimate- , tryToMatchT, tryToMatchTT, tryToMatchTTT- , ftorTryToMatch, ftorTryToMatchT, ftorTryToMatchTT- , KnownNat(..)- , Range(..)- , FreeVect(..), type (^)(), freeVector, freeCons, freeSnoc- , replicVector, indices, perfectZipWith, freeRotate- , ) where--import Data.Tagged-import Data.Semigroup--import Data.MemoTrie-import Data.VectorSpace-import Data.AffineSpace-import Data.Basis-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 ((^), Foldable(..), Traversable(..))-import Data.Traversable.Constrained---import qualified Data.Vector as Arr--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 _ = empty; fCosucc _ = empty- cozeroT = pure; cosuccT _ = empty; fCosuccT _ = 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 _ -> empty- where (Tagged k) = theNat :: Tagged k Nat-instance (KnownNat n) => KnownNat (S n) where- theNat = natSelfSucc- theNatN = natSelfSuccN- cozero _ = empty; cosucc v = pure v; fCosucc v = v- cozeroT _ = 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 = 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 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'- | List.length c == n = pure . FreeVect $ Arr.fromList c- | otherwise = 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)---empty :: Hask.Alternative m => m a-empty = Hask.empty
Data/Manifold/Cone.hs view
@@ -40,8 +40,6 @@ import Data.Manifold.Types.Stiefel import Math.LinearMap.Category -import Data.CoNat- import qualified Prelude import qualified Control.Applicative as Hask
Data/Manifold/FibreBundle.hs view
@@ -18,6 +18,7 @@ {-# LANGUAGE GADTs #-} {-# LANGUAGE DefaultSignatures #-} {-# LANGUAGE CPP #-}+{-# LANGUAGE PatternSynonyms #-} #if __GLASGOW_HASKELL__ >= 800 {-# LANGUAGE UndecidableSuperClasses #-} #endif@@ -32,6 +33,8 @@ import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine++import Math.Rotations.Class import qualified Prelude as Hask @@ -45,6 +48,20 @@ import Data.Tagged +pattern TangentBundle :: m -> Needle m -> FibreBundle m (Needle m)+pattern TangentBundle p v = FibreBundle p v++infixr 5 :@.+-- | Provided for convenience. Flipped synonym of 'FibreBundle', restricted to manifolds+-- without boundary (so the type of the whole can be inferred from its interior).+pattern (:@.) :: f -> m -> FibreBundle m f+pattern f :@. p = FibreBundle p f++-- | A zero vector in the fibre bundle at the given position. Intended to be used+-- with tangent-modifying lenses such as 'Math.Manifold.Real.Coordinates.delta'.+tangentAt :: (AdditiveGroup (Needle m), m ~ Interior m) => m -> TangentBundle m+tangentAt p = zeroV :@. p+ data TransportOnNeedleWitness k m f where TransportOnNeedle :: (ParallelTransporting (LinearFunction (Scalar (Needle m))) (Needle m) (Needle f))@@ -291,12 +308,12 @@ embed x = FibreBundle x zeroV coEmbed (FibreBundle x _) = x -instance (NaturallyEmbedded (Interior m) (Interior v), VectorSpace f)+instance (NaturallyEmbedded m v, VectorSpace f) => NaturallyEmbedded (FibreBundle m ℝ⁰) (FibreBundle v f) where embed (FibreBundle x Origin) = FibreBundle (embed x) zeroV coEmbed (FibreBundle u _) = FibreBundle (coEmbed u) Origin -instance (AdditiveGroup (Interior y), AdditiveGroup g)+instance (AdditiveGroup y, AdditiveGroup g) => NaturallyEmbedded (FibreBundle x f) (FibreBundle (x,y) (f,g)) where embed (FibreBundle x δx) = FibreBundle (x,zeroV) (δx,zeroV) coEmbed (FibreBundle (x,_) (δx,_)) = FibreBundle x δx@@ -358,3 +375,18 @@ γ = atan2 δφ δθ γc | θ < pi/2 = γ + φ | otherwise = γ - φ+++-- | @ex -> ey@, @ey -> ez@, @ez -> ex@+transformEmbeddedTangents+ :: ∀ x f v . ( NaturallyEmbedded (FibreBundle x f) (FibreBundle v v)+ , v ~ Interior v )+ => (v -> v) -> FibreBundle x f -> FibreBundle x f+transformEmbeddedTangents f p = case embed p :: FibreBundle v v of+ FibreBundle v δv -> coEmbed (FibreBundle (f v) (f δv) :: FibreBundle v v)+++instance Rotatable (FibreBundle S² ℝ²) where+ type AxisSpace (FibreBundle S² ℝ²) = ℝP²+ rotateAbout axis angle = transformEmbeddedTangents $ rotateℝ³AboutCenteredAxis axis angle+
Data/Manifold/Function/Interpolation.hs view
@@ -20,7 +20,13 @@ {-# LANGUAGE ConstraintKinds #-} module Data.Manifold.Function.Interpolation (+ -- * Interpolation functions InterpolationFunction+ , interpWeb+ , fromUncertainFunction+ -- * Specialised implementations+ -- * Local models+ , AffineModel, QuadraticModel ) where @@ -104,3 +110,11 @@ autoUpsampleAtLargeDist :: (ModellableRelation x y, LocalModel ㄇ) => ℝ -> InterpolationFunction ㄇ x y -> PointsWeb x (Shade' y) autoUpsampleAtLargeDist dmax = upsampleAtLargeDist dmax $ const evalLocalModel+++fromUncertainFunction :: (ModellableRelation x y, LocalModel ㄇ)+ => (x -> Shade' y) -- ^ Function to sample.+ -> PointsWeb x () -- ^ Minimum-resolution domain coverage.+ -> InterpolationFunction ㄇ x y+fromUncertainFunction f domain = fromPointsWeb+ $ localFmapWeb (f . (^.thisNodeCoord)) domain
+ Data/Manifold/Mesh.hs view
@@ -0,0 +1,62 @@+-- |+-- Module : Data.Manifold.Mesh+-- Copyright : (c) Justus Sagemüller 2018+-- License : GPL v3+-- +-- Maintainer : (@) jsagemue $ uni-koeln.de+-- Stability : experimental+-- Portability : portable+-- ++{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE ConstraintKinds #-}++module Data.Manifold.Mesh where++import Data.Manifold.Types.Primitive+import Math.Manifold.Core.PseudoAffine+import Data.Manifold.PseudoAffine+import Data.Simplex.Abstract++import Data.Manifold.Web+import Data.Manifold.Web.Internal+import Data.Manifold.FibreBundle++import GHC.Exts (Constraint)++-- | A mesh is a container data structure whose nodes are in some way located+-- distributed over a manifold, and are aware of the topology by way of having+-- access to their neighbours. Any such grid can be seen as a 'PointsWeb', but it+-- may have extra structure (e.g. rectangular) in addition to that.+class SimplexSpanning (MeshDomainSpace メ) => Mesh メ where+ type MeshDomainSpace メ :: *+ type MeshGridDataConstraint メ y :: Constraint+ type MeshGridDataConstraint メ y = ()+ + asWeb :: MeshGridDataConstraint メ y+ => メ y -> PointsWeb (MeshDomainSpace メ) y+ + meshSimplicesInWeb :: メ y -> [AbstractSimplex (Needle (MeshDomainSpace メ)) WebNodeId]+ + meshSimplices :: MeshGridDataConstraint メ y+ => メ y -> [SimplexF (MeshDomainSpace メ) y]+ meshSimplices mesh+ = map (fmap $ \i -> case indexWeb web i of+ Just (x,info) -> (info^.thisNodeData):@.x+ Nothing -> error $ "Faulty `Mesh` instance: node #"++show i+ ++" not in web." )+ nodeRefs+ where web = webLocalInfo $ asWeb mesh+ nodeRefs = meshSimplicesInWeb mesh+ + extrapolateGrid :: (WithField ℝ Manifold y, Connected y, MeshGridDataConstraint メ y)+ => メ y -> MeshDomainSpace メ -> y++-- | A mesh that “covers” the entire manifold, i.e. any point lies between some nodes+-- of the mesh.+class Mesh メ => CoveringMesh メ where+ interpolateGrid :: (WithField ℝ Manifold y, Connected y, MeshGridDataConstraint メ y)+ => メ y -> MeshDomainSpace メ -> y+ interpolateGrid = extrapolateGrid+
Data/Manifold/PseudoAffine.hs view
@@ -72,6 +72,7 @@ -- * Misc , alerpB, palerp, palerpB, LocallyCoercible(..), CanonicalDiffeomorphism(..) , ImpliesMetric(..), coerceMetric, coerceMetric'+ , Connected (..) ) where @@ -95,8 +96,6 @@ import Data.Tagged import Data.Manifold.Types.Primitive -import Data.CoNat- import qualified Prelude as Hask import qualified Control.Applicative as Hask @@ -278,8 +277,6 @@ instance (c) => PseudoAffine (t) where { \ a.-~.b = pure (a.-.b); } -deriveAffine(KnownNat n, FreeVect n ℝ)- instance (NumPrime s) => LocallyCoercible (ZeroDim s) (V0 s) where locallyTrivialDiffeomorphism Origin = V0 coerceNeedle _ = LinearFunction $ \Origin -> V0@@ -491,3 +488,30 @@ (⊙+^) :: ∀ x proxy . Semimanifold x => Interior x -> Needle x -> proxy x -> Interior x (⊙+^) x v _ = tp x v where Tagged tp = translateP :: Tagged x (Interior x -> Needle x -> Interior x)++++infix 6 .−.+-- | A connected manifold is one where any point can be reached by translation from+-- any other point.+class (PseudoAffine m) => Connected m where+ {-# MINIMAL #-}+ -- | Safe version of '(.-~.)'.+ (.−.) :: m -> m -> Needle m+ (.−.) = (.-~!)++instance Connected ℝ⁰+instance Connected ℝ+instance Connected ℝ¹+instance Connected ℝ²+instance Connected ℝ³+instance Connected ℝ⁴+instance Connected S¹+instance Connected S²+instance Connected ℝP⁰+instance Connected ℝP¹+instance Connected ℝP²+instance (Connected x, Connected y) => Connected (x,y)+instance (Connected x, Connected y, PseudoAffine (FibreBundle x y))+ => Connected (FibreBundle x y)+
Data/Manifold/Riemannian.hs view
@@ -64,8 +64,6 @@ import Data.Manifold.PseudoAffine import Data.Manifold.Atlas (AffineManifold) -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)
Data/Manifold/TreeCover.hs view
@@ -60,9 +60,6 @@ , stiAsIntervalMapping, spanShading , DBranch, DBranch'(..), Hourglass(..) , unsafeFmapTree- -- ** Triangulation-builders- , TriangBuild, doTriangBuild- , AutoTriang, breakdownAutoTriang -- ** External , AffineManifold, euclideanMetric ) where@@ -84,7 +81,6 @@ import Math.LinearMap.Category import Data.Tagged -import Data.SimplicialComplex import Data.Manifold.Shade import Data.Manifold.Types import Data.Manifold.Types.Primitive ((^), empty)@@ -94,7 +90,6 @@ import Data.Function.Affine import Data.Embedding-import Data.CoNat import Control.Lens (Lens', (^.), (.~), (%~), (&), _2, swapped) import Control.Lens.TH@@ -793,123 +788,6 @@ , wall <- takeWhile ((==depth) . fst . _wallID) gsc ] -----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 ℝ- fromInterior = id- toInterior = pure- translateP = Tagged (.+~^)- (.+~^) = (.+^)- semimanifoldWitness = undefined-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--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] -> Maybe x)]- -> TriangBuilder (S n) x---- -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] -> Maybe Int-optimalBottomExtension s xs- = case filter ((>0).snd)- $ zipWith ((. bottomExtendSuitability s) . (,)) [0..] xs of- [] -> empty- qs -> pure . fst . maximumBy (comparing snd) $ qs-----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))---------data AutoTriang n x where- AutoTriang :: { getAutoTriang :: ∀ t . TriangBuild t n x () } -> AutoTriang (S n) x----breakdownAutoTriang :: ∀ n n' x . (KnownNat n', n ~ S n') => AutoTriang n x -> [Simplex n x]-breakdownAutoTriang (AutoTriang t) = doTriangBuild t- - - - - - -- -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'
Data/Manifold/Types/Primitive.hs view
@@ -29,6 +29,7 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE RecordWildCards #-} {-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE LambdaCase #-} module Data.Manifold.Types.Primitive (@@ -85,6 +86,7 @@ import Data.Embedding import qualified Test.QuickCheck as QC+import qualified Test.QuickCheck.Function as QC (Function (..), functionMap) import qualified Text.Show.Pragmatic as SP @@ -205,6 +207,12 @@ instance QC.Arbitrary S⁰ where arbitrary = (\hsph -> if hsph then PositiveHalfSphere else NegativeHalfSphere) <$> QC.arbitrary+instance QC.CoArbitrary S⁰ where+ coarbitrary PositiveHalfSphere = QC.coarbitrary (2255841931547 :: Int)+ coarbitrary NegativeHalfSphere = QC.coarbitrary (1710032008738 :: Int)+instance QC.Function S⁰ where+ function = QC.functionMap (\case {PositiveHalfSphere->True; NegativeHalfSphere->False})+ (\case {True->PositiveHalfSphere; False->NegativeHalfSphere}) instance SP.Show S⁰ where showsPrec = showsPrec @@ -212,6 +220,10 @@ arbitrary = S¹Polar . (pi-) . (`mod'`(2*pi)) <$> QC.arbitrary shrink (S¹Polar φ) = S¹Polar . (pi/12*) <$> QC.shrink (φ*12/pi)+instance QC.CoArbitrary S¹ where+ coarbitrary (S¹Polar φ) = QC.coarbitrary φ+instance QC.Function S¹ where+ function = QC.functionMap (\(S¹Polar φ) -> tan $ φ/2) (S¹Polar . (*2) . atan) instance SP.Show S¹ where showsPrec p (S¹Polar φ) = showParen (p>9) $ ("S¹Polar "++) . SP.showsPrec 10 φ @@ -219,6 +231,14 @@ arbitrary = ( \θ φ -> S²Polar (θ`mod'`pi) (pi - (φ`mod'`(2*pi))) ) <$> QC.arbitrary<*>QC.arbitrary shrink (S²Polar θ φ) = uncurry S²Polar . (pi/12*^) <$> QC.shrink (θ*12/pi, φ*12/pi)+instance QC.CoArbitrary S² where+ coarbitrary (S²Polar 0 φ) = QC.coarbitrary (544317577041 :: Int)+ coarbitrary (S²Polar θ φ)+ | θ < pi = QC.coarbitrary (θ,φ)+ | otherwise = QC.coarbitrary (1771964485166 :: Int)+instance QC.Function S² where+ function = QC.functionMap (\(S²Polar θ φ) -> (cos φ, sin φ)^*tan (θ/2))+ (\(x,y) -> S²Polar (2 * (atan . sqrt $ x^2 + y^2)) (atan2 y x)) instance SP.Show S² where showsPrec p (S²Polar θ φ) = showParen (p>9) $ ("S²Polar "++) . SP.showsPrec 10 θ . (' ':) . SP.showsPrec 10 φ@@ -238,11 +258,11 @@ , φ' <- QC.shrink (φ*12/pi) ] -instance (SP.Show (Interior m), SP.Show f) => SP.Show (FibreBundle m f) where+instance (SP.Show m, SP.Show f) => SP.Show (FibreBundle m f) where showsPrec p (FibreBundle m v) = showParen (p>9) $ ("FibreBundle "++) . SP.showsPrec 10 m . (' ':) . SP.showsPrec 10 v-instance (QC.Arbitrary (Interior m), QC.Arbitrary f) => QC.Arbitrary (FibreBundle m f) where+instance (QC.Arbitrary m, QC.Arbitrary f) => QC.Arbitrary (FibreBundle m f) where arbitrary = FibreBundle <$> QC.arbitrary <*> QC.arbitrary shrink (FibreBundle m v) = [ FibreBundle m' v' | m' <- QC.shrink m
Data/Manifold/Web.hs view
@@ -36,7 +36,7 @@ -- * The web data type PointsWeb -- ** Construction- , fromWebNodes, fromShadeTree_auto, fromShadeTree, fromShaded+ , fromWebNodes, fromShadeTree_auto, fromShadeTree, fromShaded, fromGraph -- ** Lookup , nearestNeighbour, indexWeb, toGraph, webBoundary -- ** Decomposition@@ -70,6 +70,7 @@ import qualified Data.Vector as Arr import qualified Data.Vector.Mutable as MArr import qualified Data.Vector.Unboxed as UArr+import qualified Data.Array as PArr import Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NE import Data.List.FastNub (fastNub,fastNubBy)@@ -123,7 +124,7 @@ import Control.Comonad (Comonad(..)) import Control.Comonad.Cofree-import Control.Lens ((&), (%~), (^.), (.~), (+~), ix)+import Control.Lens ((&), (%~), (^.), (.~), (+~), ix, iover, indexing) import Control.Lens.TH import GHC.Generics (Generic)@@ -131,6 +132,12 @@ import Development.Placeholders +unlinkedFromWebNodes :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x))+ => (MetricChoice x) -> [(x,y)] -> PointsWeb x y+unlinkedFromWebNodes = case boundarylessWitness :: BoundarylessWitness x of+ BoundarylessWitness ->+ \mf -> unlinkedFromShaded mf . fromLeafPoints_+ fromWebNodes :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (MetricChoice x) -> [(x,y)] -> PointsWeb x y fromWebNodes = case boundarylessWitness :: BoundarylessWitness x of@@ -710,6 +717,29 @@ =<< (localOnion info []) of Just ㄇ -> ㄇ) >>= intersectShade's . (:|[info^.thisNodeData])++fromGraph :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x))+ => MetricChoice x -> Graph -> (Vertex -> (x, y)) -> PointsWeb x y+fromGraph metricf gr dataLookup+ = introduceLinks $ unlinkedFromWebNodes metricf+ [(fst (dataLookup v), v) | v <- vertices gr]+ where introduceLinks :: PointsWeb x Vertex -> PointsWeb x y+ introduceLinks (PointsWeb w) = PointsWeb $+ iover (indexing Hask.traverse)+ (\wi (Neighbourhood vert _ sclPr bound)+ -> let neighbours = gr PArr.! wi+ neighbourwis = (vertToWebNode Map.!) <$> neighbours+ (x, y) = dataLookup vert+ in Neighbourhood y+ (UArr.fromList $ subtract wi<$>neighbourwis)+ sclPr+ (snd (bestNeighbours sclPr+ [ ((), fst (dataLookup ni).-~!x)+ | ni<-neighbours ])) )+ w+ where webNodeToVert = Map.fromList assocs+ vertToWebNode = Map.fromList $ swap<$>assocs+ assocs = zip [0..] [vert | Neighbourhood vert _ _ _ <- toList w] toGraph :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => PointsWeb x y -> (Graph, Vertex -> (x, y))
+ Data/Simplex/Abstract.hs view
@@ -0,0 +1,106 @@+-- |+-- Module : Data.Simplex.Abstract+-- Copyright : (c) Justus Sagemüller 2018+-- License : GPL v3+-- +-- Maintainer : (@) jsagemue $ uni-koeln.de+-- Stability : experimental+-- Portability : portable+-- ++{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE DeriveFoldable #-}+{-# LANGUAGE DeriveTraversable #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ConstraintKinds #-}++module Data.Simplex.Abstract where++import Data.Manifold.Types.Primitive+import Math.Manifold.Core.PseudoAffine+import Data.Manifold.PseudoAffine++import Math.LinearMap.Category (spanVariance, dualNorm', (<$|), (<.>^), SimpleSpace)+import Data.VectorSpace (VectorSpace, Scalar)++import Data.Foldable (toList)+import Data.Traversable (Traversable)++import GHC.Generics++data family AbstractSimplex v x++data instance AbstractSimplex ℝ⁰ x = ℝ⁰Simplex !x+ deriving (Functor, Foldable, Traversable)+instance Applicative (AbstractSimplex ℝ⁰) where+ pure = ℝ⁰Simplex+ ℝ⁰Simplex p <*> ℝ⁰Simplex q = ℝ⁰Simplex $ p q++data instance AbstractSimplex ℝ x = ℝSimplex !x !x+ deriving (Functor, Foldable, Traversable)+data instance AbstractSimplex ℝ¹ x = ℝ¹Simplex !x !x+ deriving (Functor, Foldable, Traversable)++data instance AbstractSimplex ℝ² x = ℝ²Simplex !x !x !x+ deriving (Functor, Foldable, Traversable)++data instance AbstractSimplex ℝ³ x = ℝ³Simplex !x !x !x !x+ deriving (Functor, Foldable, Traversable)++data instance AbstractSimplex ℝ⁴ x = ℝ⁴Simplex !x !x !x !x !x+ deriving (Functor, Foldable, Traversable)++data instance AbstractSimplex (ℝ, v) x = ConeSimplex !x !(AbstractSimplex v x)+deriving instance (Functor (AbstractSimplex v)) => Functor (AbstractSimplex (ℝ,v))+deriving instance (Foldable (AbstractSimplex v)) => Foldable (AbstractSimplex (ℝ,v))+deriving instance (Traversable (AbstractSimplex v)) => Traversable (AbstractSimplex (ℝ,v))++newtype instance AbstractSimplex (GenericNeedle m) x+ = GenericSimplex (AbstractSimplex (Rep m ()) x)+deriving instance (Functor (AbstractSimplex (Rep m ())))+ => Functor (AbstractSimplex (GenericNeedle m))+deriving instance (Foldable (AbstractSimplex (Rep m ())))+ => Foldable (AbstractSimplex (GenericNeedle m))+deriving instance (Traversable (AbstractSimplex (Rep m ())))+ => Traversable (AbstractSimplex (GenericNeedle m))++newtype instance AbstractSimplex (NeedleProductSpace f g p) x+ = GenProdSimplex (AbstractSimplex (Needle (f p), Needle (g p)) x)+deriving instance (Functor (AbstractSimplex (Needle (f p), Needle (g p))))+ => Functor (AbstractSimplex (NeedleProductSpace f g p))+deriving instance (Foldable (AbstractSimplex (Needle (f p), Needle (g p))))+ => Foldable (AbstractSimplex (NeedleProductSpace f g p))+deriving instance (Traversable (AbstractSimplex (Needle (f p), Needle (g p))))+ => Traversable (AbstractSimplex (NeedleProductSpace f g p))+++type Simplex m = AbstractSimplex (Needle m) m+type SimplexF m y = AbstractSimplex (Needle m) (FibreBundle m y)++type SimplexSpanning m+ = ( WithField ℝ Manifold m, VectorSpace (Needle m)+ , Traversable (AbstractSimplex (Needle m)) )++seenFromOneVertex :: (WithField ℝ Manifold m, Foldable (AbstractSimplex (Needle m)))+ => Simplex m -> (m, [Needle m])+seenFromOneVertex s = case toList s of+ (p₀:ps) -> (p₀, [ case p.-~.p₀ of+ Just v -> v+ Nothing -> error "A simplex must always be path-connected."+ | p <- ps ])+ [] -> error "A simplex type must contain at least one value!" ++toBarycentric :: ( WithField ℝ Manifold m+ , Foldable (AbstractSimplex (Needle m))+ , SimpleSpace (Needle m) )+ => Simplex m -> m -> [ℝ]+toBarycentric s = case seenFromOneVertex s of+ (p₀, vs) -> let v's = (dualNorm' (spanVariance vs)<$|) <$> vs+ in \q -> case q.-~.p₀ of+ Just w -> let vws = (<.>^w) <$> v's+ in (1 - sum vws) : vws+ Nothing -> []
− Data/SimplicialComplex.hs
@@ -1,423 +0,0 @@--- |--- 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- -- * Triangulation-builder monad- , TriangT- , evalTriangT, runTriangT, doTriangT, getTriang- -- ** Subsimplex-references- , SimplexIT, simplexITList, lookSimplex- , lookSplxFacesIT, lookSupersimplicesIT, tgetSimplexIT- , lookVertexIT, lookSplxVerticesIT- , sharedBoundary- , distinctSimplices, NeighbouringSimplices- -- ** Building triangulations- , disjointTriangulation- , mixinTriangulation- -- * 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.Vector as Arr-import Data.List.FastNub-import qualified Data.List.NonEmpty as NE-import Data.Semigroup-import Data.Ord (comparing)--import Math.LinearMap.Category-import Data.Tagged--import Data.Manifold.Types.Primitive ((^), empty)-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---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)--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 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 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- _ -> 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 )----- | 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--- -------- | Type-level zero of kind 'Nat'.-type Zero = Z-type One = S Zero-type Two = S One-type Three = S Two-type Succ = S--
+ Math/Manifold/Real/Coordinates.hs view
@@ -0,0 +1,379 @@+-- |+-- Module : Math.Manifold.Real.Coordinates+-- Copyright : (c) Justus Sagemüller 2018+-- License : GPL v3+-- +-- Maintainer : (@) jsagemue $ uni-koeln.de+-- Stability : experimental+-- Portability : portable+-- ++{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UnicodeSyntax #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE EmptyCase #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE ScopedTypeVariables #-}++++module Math.Manifold.Real.Coordinates+ ( Coordinate, coordinate+ , HasCoordinates(..)+ -- * Vector space axes+ , HasXCoord(..), HasYCoord(..), HasZCoord(..)+ -- * Fibre bundle / tangent space diffs+ , location's+ , CoordDifferential(..)+ -- * Spherical coordinates+ , HasAzimuth(..)+ , HasZenithDistance(..)+ ) where+++import Data.Manifold.Types.Primitive+import Data.Manifold.Types.Stiefel+import Data.Manifold.PseudoAffine+import Math.LinearMap.Category+import Data.VectorSpace++import Control.Lens hiding ((<.>))+import Data.List (intercalate, transpose)++import qualified Linear as Lin++import qualified Test.QuickCheck as QC+import qualified Test.QuickCheck.Gen as QC (unGen)+import qualified Test.QuickCheck.Random as QC (mkQCGen)+import Data.Maybe (fromJust, isJust)++import qualified Numeric.IEEE as IEEE++-- | To give a custom type coordinate axes, first define an instance of this class.+class HasCoordinates m where+ -- | A unique description of a coordinate axis.+ data CoordinateIdentifier m :: *+ -- | How to use a coordinate axis for points in the containing space.+ -- This is what 'coordinate' calls under the hood.+ coordinateAsLens :: CoordinateIdentifier m -> Lens' m ℝ+ -- | Delimiters for the possible values one may choose for a given coordinate,+ -- around a point on the manifold.+ -- For example, in spherical coordinates, the 'azimuth' generally has a range+ -- of @(-'pi', 'pi')@, except at the poles where it's @(0,0)@.+ validCoordinateRange :: CoordinateIdentifier m -> m -> (ℝ,ℝ)+ validCoordinateRange _ _ = (-1/0, 1/0)++class CoordinateIsh q m | q -> m where+ useCoordinate :: CoordinateIdentifier m -> q++instance CoordinateIsh (CoordinateIdentifier m) m where+ useCoordinate = id+instance (Functor f, HasCoordinates m, a ~ (ℝ -> f ℝ), b ~ (m -> f m))+ => CoordinateIsh (a -> b) m where+ useCoordinate = coordinateAsLens++coordinate :: CoordinateIdentifier m -> Coordinate m+coordinate = useCoordinate++-- | A coordinate is a function that can be used both to determine the position+-- of a point on a manifold along the one of some family of (possibly curved) axes on+-- which it lies, and for moving the point along that axis.+-- Basically, this is a 'Lens' and can indeed be used with the '^.', '.~' and '%~'+-- operators.+-- +-- @+-- 'Coordinate' m ~ 'Lens'' m 'ℝ'+-- @+-- +-- In addition, each type may also have a way of identifying particular coordinate+-- axes. This is done with 'CoordinateIdentifier', which is what should be used+-- for /defining/ given coordinate axes.+type Coordinate m = ∀ q . CoordinateIsh q m => q++instance HasCoordinates ℝ⁰ where+ data CoordinateIdentifier ℝ⁰+ coordinateAsLens b = case b of {}++instance HasCoordinates ℝ where+ data CoordinateIdentifier ℝ = RealCoord { realAxisTfmStretch :: !ℝ }+ deriving (Eq,Show)+ coordinateAsLens (RealCoord μ) = iso (/μ) (*μ)+ {-# INLINE coordinateAsLens #-}++instance QC.Arbitrary (CoordinateIdentifier ℝ) where+ arbitrary = RealCoord . QC.getNonZero <$> QC.arbitrary+ shrink (RealCoord μ) = [ RealCoord ν | ν <- QC.shrink μ, ν/=0 ]++data OriginAxisCoord v = OriginAxisCoord+ { coordHeading :: !v -- ^ Must be conjugate to heading, i.e.+ , coordSensor :: !(DualVector v) -- ^ @'coordSensor' <.>^ 'coordHeading' = 1@.+ }+deriving instance (Show v, Show (DualVector v)) => Show (OriginAxisCoord v)+deriving instance (Eq v, Eq (DualVector v)) => Eq (OriginAxisCoord v)++originAxisCoordAsLens :: LinearSpace v => OriginAxisCoord v -> Lens' v (Scalar v)+originAxisCoordAsLens (OriginAxisCoord v dv)+ = lens (dv<.>^)+ (\w c' -> w ^+^ (c' - dv<.>^w)*^v)+{-# INLINE originAxisCoordAsLens #-}++instance (QC.Arbitrary v, InnerSpace v, v ~ DualVector v, Scalar v ~ ℝ)+ => QC.Arbitrary (OriginAxisCoord v) where+ arbitrary = QC.arbitrary `suchThatMap` \v+ -> case magnitudeSq v of+ 0 -> Nothing+ v² -> Just $ OriginAxisCoord v (v^/v²)+ shrink (OriginAxisCoord v _) = [ OriginAxisCoord w (w^/w²)+ | w <- QC.shrink v+ , let w² = magnitudeSq w+ , w² > 0 ]++instance HasCoordinates ℝ² where+ data CoordinateIdentifier ℝ² = ℝ²Coord !(OriginAxisCoord ℝ²) deriving (Eq,Show)+ coordinateAsLens (ℝ²Coord b) = originAxisCoordAsLens b+ {-# INLINE coordinateAsLens #-}++instance QC.Arbitrary ℝ² => QC.Arbitrary (CoordinateIdentifier ℝ²) where+ arbitrary = ℝ²Coord <$> QC.arbitrary+ shrink (ℝ²Coord q) = ℝ²Coord <$> QC.shrink q++instance HasCoordinates ℝ³ where+ data CoordinateIdentifier ℝ³ = ℝ³Coord !(OriginAxisCoord ℝ³) deriving (Eq,Show)+ coordinateAsLens (ℝ³Coord b) = originAxisCoordAsLens b+ {-# INLINE coordinateAsLens #-}++instance QC.Arbitrary ℝ³ => QC.Arbitrary (CoordinateIdentifier ℝ³) where+ arbitrary = ℝ³Coord <$> QC.arbitrary+ shrink (ℝ³Coord q) = ℝ³Coord <$> QC.shrink q++instance (HasCoordinates a, HasCoordinates b) => HasCoordinates (a,b) where+ data CoordinateIdentifier (a,b) = LSubspaceCoord (CoordinateIdentifier a)+ | RSubspaceCoord (CoordinateIdentifier b)+ coordinateAsLens (LSubspaceCoord ca) = _1 . coordinateAsLens ca+ coordinateAsLens (RSubspaceCoord cb) = _2 . coordinateAsLens cb+ {-# INLINE coordinateAsLens #-}++deriving instance (Eq (CoordinateIdentifier a), Eq (CoordinateIdentifier b))+ => Eq (CoordinateIdentifier (a,b))+deriving instance (Show (CoordinateIdentifier a), Show (CoordinateIdentifier b))+ => Show (CoordinateIdentifier (a,b))++instance (QC.Arbitrary (CoordinateIdentifier a), QC.Arbitrary (CoordinateIdentifier b))+ => QC.Arbitrary (CoordinateIdentifier (a,b)) where+ arbitrary = QC.oneof [LSubspaceCoord<$>QC.arbitrary, RSubspaceCoord<$>QC.arbitrary]+ shrink (LSubspaceCoord ba) = LSubspaceCoord <$> QC.shrink ba+ shrink (RSubspaceCoord bb) = RSubspaceCoord <$> QC.shrink bb++class HasCoordinates m => HasXCoord m where+ xCoord :: Coordinate m++instance HasXCoord ℝ where+ xCoord = coordinate (RealCoord 1)+ {-# INLINE xCoord #-}+instance HasXCoord ℝ² where+ xCoord = coordinate (ℝ²Coord $ OriginAxisCoord (Lin.V2 1 0) (Lin.V2 1 0))+ {-# INLINE xCoord #-}+instance HasXCoord ℝ³ where+ xCoord = coordinate (ℝ³Coord $ OriginAxisCoord (Lin.V3 1 0 0) (Lin.V3 1 0 0))+ {-# INLINE xCoord #-}+instance (HasXCoord v, HasCoordinates w) => HasXCoord (v,w) where+ xCoord = coordinate $ LSubspaceCoord xCoord++class HasYCoord m where+ yCoord :: Coordinate m++instance HasYCoord ℝ² where+ yCoord = coordinate (ℝ²Coord $ OriginAxisCoord (Lin.V2 0 1) (Lin.V2 0 1))+ {-# INLINE yCoord #-}+instance HasYCoord ℝ³ where+ yCoord = coordinate (ℝ³Coord $ OriginAxisCoord (Lin.V3 0 1 0) (Lin.V3 0 1 0))+ {-# INLINE yCoord #-}+instance HasCoordinates w => HasYCoord ((ℝ,ℝ),w) where+ yCoord = coordinate $ LSubspaceCoord yCoord+instance (HasXCoord w) => HasYCoord (ℝ,w) where+ yCoord = coordinate $ RSubspaceCoord xCoord++class HasZCoord m where+ zCoord :: Coordinate m++instance HasZCoord ℝ³ where+ zCoord = coordinate (ℝ³Coord $ OriginAxisCoord (Lin.V3 0 0 1) (Lin.V3 0 0 1))+ {-# INLINE zCoord #-}+instance HasXCoord w => HasZCoord ((ℝ,ℝ),w) where+ zCoord = coordinate $ RSubspaceCoord xCoord+instance (HasYCoord w) => HasZCoord (ℝ,w) where+ zCoord = coordinate $ RSubspaceCoord yCoord++instance (HasCoordinates b, HasCoordinates f)+ => HasCoordinates (FibreBundle b f) where+ data CoordinateIdentifier (FibreBundle b f)+ = BaseSpaceCoordinate (CoordinateIdentifier b)+ | FibreSpaceCoordinate (b -> CoordinateIdentifier f)+ coordinateAsLens (BaseSpaceCoordinate b)+ = lens (\(FibreBundle p _) -> p)+ (\(FibreBundle _ f) p -> FibreBundle p f)+ . coordinateAsLens b+ coordinateAsLens (FibreSpaceCoordinate b)+ = \φ pf@(FibreBundle p f) -> case coordinateAsLens $ b p of+ fLens -> FibreBundle p <$> fLens φ f+ validCoordinateRange (BaseSpaceCoordinate b) (FibreBundle p _) = validCoordinateRange b p+ validCoordinateRange (FibreSpaceCoordinate bf) (FibreBundle p f)+ = validCoordinateRange (bf p) f+ +instance ∀ b f . ( Show (CoordinateIdentifier b)+ , Show (CoordinateIdentifier f)+ , Eq b, Eq (CoordinateIdentifier f)+ , QC.Arbitrary b, Show b )+ => Show (CoordinateIdentifier (FibreBundle b f)) where+ showsPrec p (BaseSpaceCoordinate b)+ = showParen (p>9) $ ("BaseSpaceCoordinate "++) . showsPrec 10 b+ showsPrec p (FibreSpaceCoordinate bf)+ = showParen (p>0) $ \cont ->+ "BaseSpaceCoordinate $ \\case {"+ ++ intercalate "; " [ showsPrec 5 p . (" -> "++) . shows (bf p) $ ""+ | p <- exampleArgs ]+ ++ "... }" ++ cont+ where exampleArgs :: [b]+ exampleArgs = head $ go 1 0 2384148716156+ where go :: Int -> Int -> Int -> [[b]]+ go n tries seed+ | length candidate == n, allDifferent candidate+ , (shrunk:_) <- filter (allDifferent . map bf)+ $ shrinkElems candidate ++ [candidate]+ , [] <- take (5-n) $ go (n+1) 0 seed'+ = [shrunk]+ | tries*(n-1) > 15 = []+ | otherwise = go n (tries+1) seed'+ where candidate = take n $ generateFrom seed QC.arbitrary+ seed' = generateFrom seed QC.arbitrary+ allDifferent (x:ys) = all (x/=) ys && allDifferent ys+ allDifferent [] = True++generateFrom :: QC.CoArbitrary s => s -> QC.Gen a -> a+generateFrom seed val = QC.unGen (QC.coarbitrary seed val) (QC.mkQCGen 256592) 110818++-- | Keep length of the list, but shrink the individual elements.+shrinkElems :: QC.Arbitrary a => [a] -> [[a]]+shrinkElems l = filter ((==length l) . length) . transpose $ map QC.shrink l+++location's :: (HasCoordinates b, Interior b ~ b, HasCoordinates f)+ => CoordinateIdentifier b -> Coordinate (FibreBundle b f)+location's = coordinate . BaseSpaceCoordinate++class HasCoordinates m => CoordDifferential m where+ -- | Observe local, small variations (in the tangent space) of a coordinate.+ -- The idea is that @((p & coord+~δc) − p) ^. delta coord ≈ δc@, thus the name+ -- “'delta'”. Note however that this only holds exactly for flat spaces;+ -- in most manifolds it can (by design) only be understood in an asymptotic+ -- sense, i.e. used for evaluating directional derivatives of some function.+ -- In particular, @delta 'azimuth'@ is unstable near the poles of a sphere,+ -- because it has to compensate for the sensitive rotation of the @eφ@ unit vector.+ delta :: CoordinateIdentifier m -> Coordinate (TangentBundle m)++instance ( CoordDifferential m, f ~ Needle m, m ~ Interior m+ , QC.Arbitrary m+ , QC.Arbitrary (CoordinateIdentifier m)+ , QC.Arbitrary (CoordinateIdentifier f) )+ => QC.Arbitrary (CoordinateIdentifier (FibreBundle m f)) where+ arbitrary = QC.oneof [ BaseSpaceCoordinate <$> QC.arbitrary+ , delta <$> QC.arbitrary ]+ shrink (BaseSpaceCoordinate b) = BaseSpaceCoordinate <$> QC.shrink b+ shrink (FibreSpaceCoordinate bf) = FibreSpaceCoordinate . const+ <$> QC.shrink (bf cRef)+ where cRef₀ = QC.unGen QC.arbitrary (QC.mkQCGen 534373) 57314+ cRef = head $ QC.shrink cRef₀ ++ [cRef₀]++instance CoordDifferential ℝ where+ delta ζ = coordinate . FibreSpaceCoordinate $ const ζ+instance CoordDifferential ℝ² where+ delta ζ = coordinate . FibreSpaceCoordinate $ const ζ+instance CoordDifferential ℝ³ where+ delta ζ = coordinate . FibreSpaceCoordinate $ const ζ++instance (CoordDifferential a, CoordDifferential b) => CoordDifferential (a,b) where+ delta (LSubspaceCoord ba) = coordinate $ case delta ba of+ FibreSpaceCoordinate bf -> FibreSpaceCoordinate $ \(δa,_) -> LSubspaceCoord $ bf δa+ delta (RSubspaceCoord bb) = coordinate $ case delta bb of+ FibreSpaceCoordinate bf -> FibreSpaceCoordinate $ \(_,δb) -> RSubspaceCoord $ bf δb++instance HasCoordinates S¹ where+ data CoordinateIdentifier S¹ = S¹Azimuth deriving (Eq,Show)+ coordinateAsLens S¹Azimuth = lens φParamS¹ (const S¹Polar)+ validCoordinateRange S¹Azimuth _ = (-pi, pi)++instance QC.Arbitrary (CoordinateIdentifier S¹) where+ arbitrary = return S¹Azimuth++class HasAzimuth m where+ azimuth :: Coordinate m++instance HasAzimuth S¹ where+ azimuth = coordinate S¹Azimuth++instance CoordDifferential S¹ where+ delta S¹Azimuth = coordinate . FibreSpaceCoordinate $ const xCoord+ +instance HasCoordinates S² where+ data CoordinateIdentifier S² = S²ZenithAngle | S²Azimuth deriving (Eq,Show)+ coordinateAsLens S²ZenithAngle = lens ϑParamS² (\(S²Polar _ φ) θ -> S²Polar θ φ)+ coordinateAsLens S²Azimuth = lens φParamS² (\(S²Polar θ _) φ -> S²Polar θ φ)+ validCoordinateRange S²ZenithAngle _ = (0, pi)+ validCoordinateRange S²Azimuth (S²Polar θ _)+ | θ>0 && θ<pi = (-pi, pi)+ | otherwise = (0, 0)++instance QC.Arbitrary (CoordinateIdentifier S²) where+ arbitrary = QC.elements [S²Azimuth, S²ZenithAngle]++instance HasAzimuth S² where+ azimuth = coordinate S²Azimuth+ +class HasZenithDistance m where+ zenithAngle :: Coordinate m++instance HasZenithDistance S² where+ zenithAngle = coordinate S²ZenithAngle++instance CoordDifferential S² where+ delta S²ZenithAngle = coordinate . FibreSpaceCoordinate+ $ \(S²Polar θ φ) -> let eθ+ | θ < pi/2 = embed . S¹Polar $ φ+ | otherwise = embed . S¹Polar $ -φ+ in ℝ²Coord $ OriginAxisCoord eθ eθ+ delta S²Azimuth = coordinate . FibreSpaceCoordinate+ $ \(S²Polar θ φ) -> let eφ+ | θ < pi/2 = embed . S¹Polar $ φ + pi/2+ | otherwise = embed . S¹Polar $ pi/2 - φ+ sθ = sin θ + tiny+ -- ^ Right at the poles, azimuthal movements+ -- become inexpressible, which manifests itself+ -- in giving infinite diffs. Moreover,+ -- we also can't retrieve tangent diffs we put+ -- in anymore. Arguably, this just expresses+ -- the fact that azimuthal changes are meaningless+ -- at the poles, however it violates the lens+ -- laws, so prevent the infinity by keeping+ -- sin θ very slightly above 0.+ in ℝ²Coord $ OriginAxisCoord (eφ^*sθ) (eφ^/sθ)++-- | @2e-162@. A value that's so small that it can't notably disturb any nonzero value+-- you might realistically encounter (i.e. @x + tiny == x@), but still large enough+-- that ratios can reliably be represented (i.e. @x * tiny / tiny == x@).+tiny :: ℝ+tiny = IEEE.bisectIEEE IEEE.minNormal IEEE.epsilon+ ++suchThatMap :: QC.Gen a -> (a -> Maybe b) -> QC.Gen b+#if !MIN_VERSION_QuickCheck(2,11,0)+gen `suchThatMap` f =+ fmap fromJust $ fmap f gen `QC.suchThat` isJust+#else+suchThatMap = QC.suchThatMap+#endif
manifolds.cabal view
@@ -1,5 +1,5 @@ Name: manifolds-Version: 0.4.5.0+Version: 0.5.0.0 Category: Math Synopsis: Coordinate-free hypersurfaces Description: Manifolds, a generalisation of the notion of “smooth curves” or surfaces,@@ -40,7 +40,7 @@ Library Build-Depends: base>=4.5 && < 6- , manifolds-core == 0.4.5.0+ , manifolds-core == 0.5.0.0 , transformers , vector-space>=0.8 , free-vector-spaces>=0.1.5@@ -48,12 +48,15 @@ , MemoTrie , vector , linearmap-category >= 0.3.4 && < 0.4+ , spatial-rotations >= 0.1 && < 0.2 , containers+ , array , comonad , free , semigroups , void , number-show >= 0.1 && < 0.2+ , ieee754 >= 0.8 && < 1 , tagged , deepseq , placeholders@@ -79,10 +82,11 @@ Data.Manifold.Shade Data.Manifold.Web Data.Manifold.Web.Internal+ Data.Manifold.Mesh Data.Manifold.DifferentialEquation Data.Manifold.Function.LocalModel Data.Manifold.Function.Interpolation- Data.SimplicialComplex+ Data.Simplex.Abstract Data.Function.Differentiable Data.Function.Affine Data.Manifold.Types@@ -91,12 +95,12 @@ Data.Manifold.Atlas Data.Manifold.FibreBundle Data.Manifold.Riemannian+ Math.Manifold.Real.Coordinates Math.Manifold.Embedding.Simple.Class Other-modules: Data.List.FastNub Data.Manifold.Types.Primitive Data.SetLike.Intersection Data.Manifold.Cone- Data.CoNat Data.Embedding Data.Manifold.Function.Quadratic Data.Function.Differentiable.Data@@ -124,6 +128,7 @@ , containers , vector-space , linear+ , spatial-rotations , constrained-categories , linearmap-category , lens
test/tasty/test.hs view
@@ -18,6 +18,7 @@ import Data.Manifold.PseudoAffine import Data.Manifold.FibreBundle import Data.Manifold.TreeCover+import Math.Manifold.Real.Coordinates import Data.Manifold.Web import Data.Manifold.Web.Internal import Data.Manifold.Function.LocalModel@@ -27,14 +28,18 @@ import Linear.V2 (V2(V2)) import Linear.V3 (V3(V3)) import Math.LinearMap.Category-import Prelude hiding (id, fst, snd)+import Prelude hiding (id, fst, snd, asinh) import Control.Category.Constrained (id) import Control.Arrow.Constrained (fst,snd) +import Math.Rotations.Class+import Data.Simplex.Abstract+ import Test.Tasty import Test.Tasty.HUnit import qualified Test.Tasty.QuickCheck as QC import Test.Tasty.QuickCheck ((==>))+import Data.Typeable import Data.Foldable (toList) import Data.List (nub)@@ -43,6 +48,8 @@ import Control.Arrow import Control.Lens hiding ((<.>)) +import Data.Fixed (mod')+ import qualified Text.Show.Pragmatic as SP @@ -103,45 +110,45 @@ ] , testGroup "1-sphere tangent bundle" [ testCase "North pole"- $ embed (FibreBundle (S¹Polar $ pi/2) 1 :: TangentBundle S¹)+ $ embed (TangentBundle (S¹Polar $ pi/2) 1) @?≈ (FibreBundle (V2 0 1) (V2 (-1) 0) :: TangentBundle ℝ²) , testCase "South pole"- $ embed (FibreBundle (S¹Polar $ -pi/2) 1 :: TangentBundle S¹)+ $ embed (TangentBundle (S¹Polar $ -pi/2) 1) @?≈ (FibreBundle (V2 0 (-1)) (V2 1 0) :: TangentBundle ℝ²) , testCase "45°"- $ embed (FibreBundle (S¹Polar $ pi/4) 1 :: TangentBundle S¹)+ $ embed (TangentBundle (S¹Polar $ pi/4) 1) @?≈ (FibreBundle (V2 1 1^/sqrt 2) (V2 (-1) 1^/sqrt 2) :: TangentBundle ℝ²) ] , testGroup "2-sphere tangent bundle" [ testCase "North pole, x-dir"- $ embed (FibreBundle (S²Polar 0 0) (V2 1 0) :: TangentBundle S²)+ $ embed (TangentBundle (S²Polar 0 0) (V2 1 0)) @?≈ (FibreBundle (V3 0 0 1) (V3 1 0 0) :: TangentBundle ℝ³) , testCase "North pole (alternative φ), x-dir"- $ embed (FibreBundle (S²Polar 0 1.524) (V2 1 0) :: TangentBundle S²)+ $ embed (TangentBundle (S²Polar 0 1.524) (V2 1 0)) @?≈ (FibreBundle (V3 0 0 1) (V3 1 0 0) :: TangentBundle ℝ³) , testCase "North pole, y-dir"- $ embed (FibreBundle (S²Polar 0 0) (V2 0 1) :: TangentBundle S²)+ $ embed (TangentBundle (S²Polar 0 0) (V2 0 1)) @?≈ (FibreBundle (V3 0 0 1) (V3 0 1 0) :: TangentBundle ℝ³) , testCase "Close to north pole"- $ embed (FibreBundle (S²Polar 1e-11 0.602) (V2 3.7 1.1) :: TangentBundle S²)+ $ embed (TangentBundle (S²Polar 1e-11 0.602) (V2 3.7 1.1)) @?≈ (FibreBundle (V3 0 0 1) (V3 3.7 1.1 0) :: TangentBundle ℝ³) , testCase "South pole, x-dir"- $ embed (FibreBundle (S²Polar pi 0) (V2 1 0) :: TangentBundle S²)+ $ embed (TangentBundle (S²Polar pi 0) (V2 1 0)) @?≈ (FibreBundle (V3 0 0 (-1)) (V3 (-1) 0 0) :: TangentBundle ℝ³) , testCase "South pole, y-dir"- $ embed (FibreBundle (S²Polar pi 0) (V2 0 1) :: TangentBundle S²)+ $ embed (TangentBundle (S²Polar pi 0) (V2 0 1)) @?≈ (FibreBundle (V3 0 0 (-1)) (V3 0 1 0) :: TangentBundle ℝ³) , testCase "Close to south pole"- $ embed (FibreBundle (S²Polar (pi-1e-11) 0.602) (V2 3.7 1.1) :: TangentBundle S²)+ $ embed (TangentBundle (S²Polar (pi-1e-11) 0.602) (V2 3.7 1.1)) @?≈ (FibreBundle (V3 0 0 (-1)) (V3 (-3.7) 1.1 0) :: TangentBundle ℝ³) , testCase "Equator, y-dir"- $ embed (FibreBundle (S²Polar (pi/2) 0) (V2 0 1) :: TangentBundle S²)+ $ embed (TangentBundle (S²Polar (pi/2) 0) (V2 0 1)) @?≈ (FibreBundle (V3 1 0 0) (V3 0 1 0) :: TangentBundle ℝ³) , testCase "Equator, x-dir"- $ embed (FibreBundle (S²Polar (pi/2) (pi/2)) (V2 1 0) :: TangentBundle S²)+ $ embed (TangentBundle (S²Polar (pi/2) (pi/2)) (V2 1 0)) @?≈ (FibreBundle (V3 0 1 0) (V3 (-1) 0 0) :: TangentBundle ℝ³) , testCase "Equator, z-dir"- $ embed (FibreBundle (S²Polar (pi/2) 0) (V2 1 0) :: TangentBundle S²)+ $ embed (TangentBundle (S²Polar (pi/2) 0) (V2 1 0)) @?≈ (FibreBundle (V3 1 0 0) (V3 0 0 (-1)) :: TangentBundle ℝ³) ] ]@@ -166,17 +173,111 @@ , testGroup "Special properties of translations" [ testGroup "2-sphere" [ QC.testProperty "S²-movement as rotation in ℝ³"- $ \p v -> let FibreBundle pCart vCart :: TangentBundle ℝ³- = embed (FibreBundle p v :: TangentBundle S²)+ $ \p v -> magnitude v < 1e6+ ==> let TangentBundle pCart vCart :: TangentBundle ℝ³+ = embed $ TangentBundle p v q = p .+~^ v :: S² qCart = embed q :: ℝ³ axis = pCart `cross3` qCart- FibreBundle _ axisProj :: TangentBundle S²- = coEmbed (FibreBundle pCart axis :: TangentBundle ℝ³)+ TangentBundle _ axisProj :: TangentBundle S²+ = coEmbed $ TangentBundle pCart axis in vCart <.> axis + 1 ≈ 1 -- i.e. the movement vector is always && v <.> axisProj + 1 ≈ 1 -- orthogonal to the rotation axis. ] ]+ , testGroup "Rotation"+ [ testCase "Pole to eqt / prime meridian"+ $ let rotated = 90° yAxis $ V2 1 0 :@. S²Polar 0 0+ in V2 (rotated ^. delta zenithAngle) (rotated ^. delta azimuth)+ @?≈ V2 1 0+ , testCase "Pole to eqt / 90°E"+ $ let rotated = 90° xAxis $ V2 1 0 :@. S²Polar 0 0+ in V2 (rotated ^. delta zenithAngle) (rotated ^. delta azimuth)+ @?≈ V2 0 1+ , QC.testProperty "Undo – arbitrary axis / angle and points in 𝑇S²."+ $ \ax ψ p -> rotateAboutThenUndo @(TangentBundle S²) ax ψ p ≈ p+ ]+ , testGroup "Coordinates"+ [ testGroup "Single dimension"+ [ QC.testProperty "Access" $ \x -> x^.xCoord ≈ x+ , QC.testProperty "Update" $ \x₀ x₁ -> (xCoord.~x₁) x₀ ≈ (x₁ :: ℝ) ]+ , testGroup "x-coordinate"+ [ QC.testProperty "Access" $ \x y -> V2 x y^.xCoord ≈ x+ , QC.testProperty "Update" $ \x₀ y x₁ -> (xCoord.~x₁) (V2 x₀ y) ≈ V2 x₁ y ]+ , testGroup "y-coordinate"+ [ QC.testProperty "Access" $ \x y -> V2 x y^.yCoord ≈ y+ , QC.testProperty "Update" $ \x y₀ y₁ -> (yCoord.~y₁) (V2 x y₀) ≈ V2 x y₁ ]+ , testGroup "z-coordinate"+ [ QC.testProperty "Access" $ \x y z -> V3 x y z^.zCoord ≈ z+ , QC.testProperty "Update" $ \x y z₀ z₁ -> (zCoord.~z₁) (V3 x y z₀) ≈ V3 x y z₁ ]+ , testGroup "Lens laws"+ [ coordinateLensLaws @ℝ+ , coordinateLensLaws @ℝ²+ , coordinateLensLaws @ℝ³+ , coordinateLensLaws @S¹+ , coordinateLensLaws @S²+ , coordinateLensLaws @(TangentBundle ℝ)+ , coordinateLensLaws @(TangentBundle ℝ²)+ , coordinateLensLaws @(TangentBundle ℝ³)+ , coordinateLensLaws @(TangentBundle S¹)+ , coordinateLensLaws @(TangentBundle S²)+ ]+ , testGroup "Finite differences"+ [ QC.testProperty "ℝ" $ coordinateFiniteDifference @ℝ 1 1e6 1e100+ , QC.testProperty "ℝ²" $ coordinateFiniteDifference @ℝ² 1 1e6 1e100+ , QC.testProperty "ℝ³" $ coordinateFiniteDifference @ℝ³ 1 1e6 1e100+ , QC.testProperty "(ℝ,ℝ)" $ coordinateFiniteDifference @(ℝ,ℝ) 1 1e6 1e100+ , QC.testProperty "S¹" $ coordinateFiniteDifference @S¹ 1 1e6 (2*pi)+ , QC.testProperty "S² (unlimited)"+ . QC.expectFailure $ coordinateFiniteDifference @S² 0.5 pi (2*pi)+ , QC.testProperty "S²" $ \p@(S²Polar θ _)+ -> let poleDist = sin θ+ in poleDist > 0.1+ ==> coordinateFiniteDifference @S² (poleDist^2 * 1e-6)+ (poleDist/2)+ (2*pi) p+ ]+ , testGroup "Location"+ [ QC.testProperty "S²" $ \p v+ -> TangentBundle @S² p v ^. location's azimuth ≈ p^.azimuth+ ]+ , testGroup "x-coordinate diff"+ [ QC.testProperty "Access" $ \x y δx δy+ -> (TangentBundle (V2 x y) (V2 δx δy))+ ^.delta xCoord ≈ δx+ , QC.testProperty "Update" $ \x y δx₀ δx₁ δy+ -> (delta xCoord.~δx₁)+ (TangentBundle (V2 x y) (V2 δx₀ δy))+ ≈ TangentBundle (V2 x y) (V2 δx₁ δy) ]+ , testGroup "Spheres"+ [ testGroup "S¹"+ [ QC.testProperty "Azimuth access" $ \φ -> S¹Polar φ^.azimuth ≈ φ+ , QC.testProperty "Azimuth update" $ \p φ -> (azimuth .~ φ) p ≈ S¹Polar φ+ ]+ , testGroup "S²"+ [ QC.testProperty "Azimuth access" $ \θ φ -> S²Polar θ φ^.azimuth ≈ φ+ , QC.testProperty "Azimuth update" $ \θ φ₀ φ₁+ -> (azimuth .~ φ₁) (S²Polar θ φ₀) ≈ S²Polar θ φ₁+ , QC.testProperty "Zenith-distance access" $ \θ φ -> S²Polar θ φ^.zenithAngle ≈ θ+ , QC.testProperty "Zenith-distance update" $ \θ₀ θ₁ φ+ -> (zenithAngle .~ θ₁) (S²Polar θ₀ φ) ≈ S²Polar θ₁ φ+ , testGroup "Tangent space examples"+ [ testCase "Zenith-angle at equator | prime meridian"+ $ (TangentBundle (S²Polar (pi/2-1e-6) 0) (V2 1 0))+ ^. delta zenithAngle @?≈ 1+ , testCase "Azimuth at just north of equator | prime meridian"+ $ (TangentBundle (S²Polar (pi/2-1e-6) 0) (V2 0 1))+ ^. delta azimuth @?≈ 1+ , testCase "Azimuth at just north of equator | 90°E"+ $ (TangentBundle (S²Polar (pi/2-1e-6) (pi/2)) (V2 1 0))+ ^. delta azimuth @?≈ -1+ , testCase "Azimuth at 45°N | prime meridian"+ $ (TangentBundle (S²Polar (pi/4) 0) (V2 0 1))+ ^. delta azimuth @?≈ sqrt 2+ ]+ ]+ ]+ ] , testGroup "Parallel transport" [ testGroup "Displacement cancellation" [ QC.testProperty "Real vector space" (parTransportAssociativity @(ℝ,ℝ))@@ -257,12 +358,11 @@ (S²Polar (abs θ₀) (if θ₀>0 then 0 else pi)) (S²Polar (abs θ₁) (if θ₁>0 then 0 else pi)) , QC.testProperty "Rotation axis – heading-vector"- $ \p v -> let q = p .+~^ v :: S²+ $ \p v -> magnitude v < 1e6+ ==> let q = p .+~^ v :: S² w = parallelTransport p v v- FibreBundle pCart vCart- = embed (FibreBundle p v :: TangentBundle S²) :: TangentBundle ℝ³- FibreBundle qCart wCart- = embed (FibreBundle q w :: TangentBundle S²) :: TangentBundle ℝ³+ vCart :@. pCart = embed (v:@.p) :: TangentBundle ℝ³+ wCart :@. qCart = embed (w:@.q) :: TangentBundle ℝ³ pxv = pCart`cross3`vCart qxw = qCart`cross3`wCart in QC.counterexample@@ -276,11 +376,9 @@ $ pxv ≈ qxw , QC.testProperty "Rotation axis – arbitrary vectors" $ \p v f -> let q = p .+~^ v :: S²- g = parallelTransport p v f- FibreBundle pCart fCart- = embed (FibreBundle p f :: TangentBundle S²) :: TangentBundle ℝ³- FibreBundle qCart gCart- = embed (FibreBundle q g :: TangentBundle S²) :: TangentBundle ℝ³+ g = parallelTransport p v f :: Needle S²+ fCart :@. pCart = embed (f :@. p) :: TangentBundle ℝ³+ gCart :@. qCart = embed (g :@. q) :: TangentBundle ℝ³ infix 7 × (×) = cross3 pxq = pCart×qCart@@ -305,6 +403,16 @@ -- ‖𝐚×𝐛‖ = ‖𝐚‖·‖𝐛‖.) ] ]+ , testGroup "Simplices"+ [ testGroup "Barycentric coordinates"+ [ QC.testProperty "In ℝ²"+ $ \p q r μ ν -> not (p≈q || q≈r || r≈p)+ ==> let λ = 1-μ-ν+ in toBarycentric (ℝ²Simplex p q r :: Simplex ℝ²)+ (p^*λ ^+^ q^*μ ^+^ r^*ν)+ ?≈! [ λ, μ, ν]+ ]+ ] , testGroup "Graph structure of webs" [ testCase "Manually-defined empty web." $ toList (fst $ toGraph emptyWeb) @?= []@@ -535,17 +643,21 @@ , 565.5193483520385 ] ] :: PointsWeb ℝ () )) @?= [ [1], [0,2], [1,3], [4,2], [3] ]- , QC.testProperty "Random 1D web should be strongly connected"+ , adjustOption (\(QC.QuickCheckTests n)+ -> QC.QuickCheckTests (ceiling . sqrt $ fromIntegral n))+ $ testGroup "QuickCheck"+ [ QC.testProperty "Random 1D web should be strongly connected" $ \ps -> length ps >= 2 ==> length (Graph.scc . fst $ toGraph ( fromWebNodes euclideanMetric [(x, ()) | x<-Set.toList ps] :: PointsWeb ℝ () ) ) == 1- , QC.testProperty "Random 1D web should have only 2 boundary-points"+ , QC.testProperty "Random 1D web should have only 2 boundary-points" $ \ps -> length ps >= 2 ==> length (webBoundary (fromWebNodes euclideanMetric [(x, ()) | x<-Set.toList ps] :: PointsWeb ℝ () ) ) == 2+ ] ] , testGroup "Shades" [ testCase "Equality of `Shade`s"@@ -772,8 +884,13 @@ (HemisphereℝP²Polar (pi/2) $ ϕ - pi) | otherwise = abs (φ - ϕ) < η -instance (AEq (Interior m), AEq f) => AEq (FibreBundle m f) where+instance (AEq m, AEq f) => AEq (FibreBundle m f) where fuzzyEq η (FibreBundle p v) (FibreBundle q w) = fuzzyEq η p q && fuzzyEq η v w++instance (AEq a) => AEq [a] where+ fuzzyEq _ [] [] = True+ fuzzyEq η (x:xs) (y:ys) = fuzzyEq η x y && fuzzyEq η xs ys+ fuzzyEq _ _ _ = False infix 1 @?≈ (@?≈) :: (AEq e, Show e) => e -> e -> Assertion@@ -781,14 +898,25 @@ | a≈b = return () | otherwise = assertFailure $ "Expected "++show b++", but got "++show a +infix 4 ?≈!+(?≈!) :: (AEq e, SP.Show e) => e -> e -> QC.Property+a?≈!b = QC.counterexample ("Expected "++SP.show b++", but got "++SP.show a) $ a≈b+ instance QC.Arbitrary ℝ² where arbitrary = (\(x,y)->V2 x y) <$> QC.arbitrary shrink (V2 x y) = V2 <$> ((/12)<$>QC.shrink (x*12)) <*> ((/12)<$>QC.shrink (y*12))+instance QC.Arbitrary ℝ³ where+ arbitrary = (\(x,y,z)->V3 x y z) <$> QC.arbitrary+ shrink (V3 x y z) = V3 <$> ((/12)<$>QC.shrink (x*12))+ <*> ((/12)<$>QC.shrink (y*12))+ <*> ((/12)<$>QC.shrink (z*12)) -nearlyAssociative :: ∀ m . (AEq m, Semimanifold m, Interior m ~ m)- => m -> Needle m -> Needle m -> Bool-nearlyAssociative p v w = (p .+~^ v) .+~^ w ≈ (p .+~^ (v^+^w) :: m)+nearlyAssociative :: ∀ m . ( AEq m, Semimanifold m, Interior m ~ m+ , InnerSpace (Needle m), RealFloat (Scalar (Needle m)) )+ => m -> Needle m -> Needle m -> QC.Property+nearlyAssociative p v w = maximum (map magnitude [v,w]) < 1e6+ ==> (p .+~^ v) .+~^ w ≈ (p .+~^ (v^+^w) :: m) originCancellation :: ∀ m . (AEq m, Manifold m, Show m, Show (Needle m)) => m -> m -> QC.Property@@ -807,11 +935,12 @@ p' = coEmbed ep embeddingTangentiality :: ∀ m n . ( Semimanifold m, Semimanifold n+ , Interior m ~ m, Interior n ~ n , NaturallyEmbedded n m , NaturallyEmbedded (TangentBundle n) (TangentBundle m) , SP.Show n, AEq n , InnerSpace (Needle n), RealFloat (Scalar (Needle n)) )- => Scalar (Needle n) -> Interior n -> Needle n -> QC.Property+ => Scalar (Needle n) -> n -> Needle n -> QC.Property embeddingTangentiality consistRadius p vub = QC.counterexample ("p+v = "++SP.show q++", coEmbed (embed p+v) = "++SP.show q') $ fuzzyEq (unitEpsilon @n * (1+rvub^2)) q q'@@ -820,9 +949,7 @@ q, q' :: n q = p .+~^ v q' = coEmbed $ (pEmbd .+~^ vEmbd :: m)- o :: TangentBundle n- o = FibreBundle p v- FibreBundle pEmbd vEmbd = embed o :: TangentBundle m+ TangentBundle pEmbd vEmbd = embed (TangentBundle p v) nearbyTangentSpaceEmbedding :: ∀ m n . ( Semimanifold m, Semimanifold n@@ -842,12 +969,9 @@ q :: n q = p .+~^ v :: n qEmbd = embed q :: m- FibreBundle _ fReProj :: TangentBundle n- = coEmbed (FibreBundle qEmbd fEmbd :: TangentBundle m)+ fReProj :@. _= coEmbed (fEmbd :@. qEmbd) :: TangentBundle n g = parallelTransport p v f- o :: TangentBundle n- o = FibreBundle p f- FibreBundle pEmbd fEmbd = embed o :: TangentBundle m+ fEmbd :@. pEmbd = embed (f:@.p) :: TangentBundle m parTransportAssociativity :: ∀ m . ( AEq m, Manifold m, SP.Show m@@ -882,6 +1006,80 @@ sphereParallelTransportTest p q (v:vs) (w:ws) = (parallelTransport p (q.-~!p) vSph @?≈ wSph) >> sphereParallelTransportTest p q vs ws- where [FibreBundle _ vSph, FibreBundle _ wSph]- = [ coEmbed (FibreBundle (embed o) u :: TangentBundle ℝ³) :: TangentBundle S²+ where [vSph:@._, wSph:@._]+ = [ coEmbed (u :@. embed o :: TangentBundle ℝ³) :: TangentBundle S² | (o,u) <- [(p,v), (q,w)] ]+++coordinateLensLaws :: ∀ m . ( Typeable m, HasCoordinates m+ , Show m, Show (CoordinateIdentifier m)+ , SP.Show m, AEq m+ , QC.Arbitrary m, QC.Arbitrary (CoordinateIdentifier m) )+ => TestTree +coordinateLensLaws = testGroup (show $ typeRep ([]::[m]))+ [ QC.testProperty "Retrieval" retrieval+ , QC.testProperty "Identity-pasting" idPasting+ , QC.testProperty "Putting twice" twicePutting+ ]+ where retrieval :: CoordinateIdentifier m -> m -> ℝ -> QC.Property+ retrieval c p a = (QC.counterexample ("Got back "++SP.show retrieved)+ $ retrieved ≈ x)+ where retrieved = (coordinate c.~x) p ^. coordinate c+ x = constrainToRange (validCoordinateRange c p) a+ idPasting :: CoordinateIdentifier m -> m -> QC.Property+ idPasting c p = (QC.counterexample ("Putting the viewed coordinate back in gives "+ ++ SP.show backPasted)+ $ backPasted ≈ p)+ where backPasted = coordinate c .~ (p^.coordinate c) $ p+ twicePutting :: CoordinateIdentifier m -> m -> ℝ -> QC.Property+ twicePutting c p a = (QC.counterexample ("Second putting made it "++SP.show dubPut)+ $ dubPut ≈ singlyPut)+ where singlyPut = p & coordinate c .~ x+ dubPut = singlyPut & coordinate c .~ x+ x = constrainToRange (validCoordinateRange c p) a++constrainToRange :: (ℝ,ℝ) -> ℝ -> ℝ+constrainToRange (lul,uul) = \x -> sinh $ m + rd * tanh (asinh x / (4 + rd))+ where l = asinh $ max (-huge) lul+ u = asinh $ min huge uul+ rd = (u-l)/2+ m = l + rd+ huge = 1e9++-- | 'Prelude.asinh' is (as of GHC-8.2) unstable for negative arguments, see+-- <https://ghc.haskell.org/trac/ghc/ticket/14927>+asinh :: RealFloat a => a -> a+asinh x+ | x > 1e20 = log 2 + log x+ | x < 0 = -asinh (-x)+ | otherwise = log $ x + sqrt (1 + x^2)++++coordinateFiniteDifference :: ∀ m .+ ( Semimanifold m, HasCoordinates m, m ~ Interior m+ , HasCoordinates (Needle m), CoordDifferential m+ , AEq (Needle m), InnerSpace (Needle m), Scalar (Needle m) ~ ℝ+ , SP.Show m )+ => ℝ -- ^ Radius of consistency (within which we expect order-1 accuracy)+ -> ℝ -- ^ Radius of stability (without we don't expect sensible results at all)+ -> ℝ -- ^ Modularity+ -> m -> CoordinateIdentifier m -> Needle m -> QC.Property+coordinateFiniteDifference consistRadius stabilRadius modl p c vub+ = QC.counterexample ("Fin. diff: "++SP.show finitesimal+ ++", tangential component: "++SP.show infinitesimal+ ++"\n(q = "++SP.show q++")")+ $ rvub * consistRadius < stabilRadius+ ==> fuzzyEq (unitEpsilon @(Needle m) * (1+rvub^2))+ (orthoCorrection + finitesimal) (orthoCorrection + infinitesimal)+ where rvub = realToFrac $ magnitude vub+ v = vub ^* consistRadius+ q = p .+~^ v+ infinitesimal = (FibreBundle p v ^. delta c)`mod'`modl+ finitesimal = (q^.coordinate c - p^.coordinate c)`mod'`modl+ orthoCorrection = signum infinitesimal+++rotateAboutThenUndo :: Rotatable m => AxisSpace m -> S¹ -> m -> m+rotateAboutThenUndo ax g@(S¹Polar w) p+ = rotateAbout ax (S¹Polar $ -w) $ rotateAbout ax g p