diff --git a/Data/CoNat.hs b/Data/CoNat.hs deleted file mode 100644 --- a/Data/CoNat.hs +++ /dev/null @@ -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 diff --git a/Data/Manifold/Cone.hs b/Data/Manifold/Cone.hs --- a/Data/Manifold/Cone.hs +++ b/Data/Manifold/Cone.hs @@ -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 diff --git a/Data/Manifold/FibreBundle.hs b/Data/Manifold/FibreBundle.hs --- a/Data/Manifold/FibreBundle.hs +++ b/Data/Manifold/FibreBundle.hs @@ -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 + diff --git a/Data/Manifold/Function/Interpolation.hs b/Data/Manifold/Function/Interpolation.hs --- a/Data/Manifold/Function/Interpolation.hs +++ b/Data/Manifold/Function/Interpolation.hs @@ -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 diff --git a/Data/Manifold/Mesh.hs b/Data/Manifold/Mesh.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/Mesh.hs @@ -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 + diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs --- a/Data/Manifold/PseudoAffine.hs +++ b/Data/Manifold/PseudoAffine.hs @@ -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) + diff --git a/Data/Manifold/Riemannian.hs b/Data/Manifold/Riemannian.hs --- a/Data/Manifold/Riemannian.hs +++ b/Data/Manifold/Riemannian.hs @@ -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) diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs --- a/Data/Manifold/TreeCover.hs +++ b/Data/Manifold/TreeCover.hs @@ -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' diff --git a/Data/Manifold/Types/Primitive.hs b/Data/Manifold/Types/Primitive.hs --- a/Data/Manifold/Types/Primitive.hs +++ b/Data/Manifold/Types/Primitive.hs @@ -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 diff --git a/Data/Manifold/Web.hs b/Data/Manifold/Web.hs --- a/Data/Manifold/Web.hs +++ b/Data/Manifold/Web.hs @@ -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)) diff --git a/Data/Simplex/Abstract.hs b/Data/Simplex/Abstract.hs new file mode 100644 --- /dev/null +++ b/Data/Simplex/Abstract.hs @@ -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 -> [] diff --git a/Data/SimplicialComplex.hs b/Data/SimplicialComplex.hs deleted file mode 100644 --- a/Data/SimplicialComplex.hs +++ /dev/null @@ -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 - - diff --git a/Math/Manifold/Real/Coordinates.hs b/Math/Manifold/Real/Coordinates.hs new file mode 100644 --- /dev/null +++ b/Math/Manifold/Real/Coordinates.hs @@ -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 && θ 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 diff --git a/manifolds.cabal b/manifolds.cabal --- a/manifolds.cabal +++ b/manifolds.cabal @@ -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 diff --git a/test/tasty/test.hs b/test/tasty/test.hs --- a/test/tasty/test.hs +++ b/test/tasty/test.hs @@ -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 +-- +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