manifolds 0.4.4.0 → 0.4.5.0

raw patch · 15 files changed

+1194/−142 lines, 15 filesdep +QuickCheckdep ~constrained-categoriesdep ~free-vector-spacesdep ~manifolds-corePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: QuickCheck

Dependency ranges changed: constrained-categories, free-vector-spaces, manifolds-core

API changes (from Hackage documentation)

Files

• Data/Manifold/Atlas.hs +11/−11
• Data/Manifold/DifferentialEquation.hs +2/−2
• Data/Manifold/FibreBundle.hs +360/−0
• Data/Manifold/Function/Interpolation.hs +106/−0
• Data/Manifold/Function/LocalModel.hs +17/−0
• Data/Manifold/PseudoAffine.hs +55/−31
• Data/Manifold/Riemannian.hs +8/−8
• Data/Manifold/TreeCover.hs +26/−13
• Data/Manifold/Types.hs +6/−4
• Data/Manifold/Types/Primitive.hs +74/−56
• Data/Manifold/Web.hs +61/−2
• Data/Manifold/Web/Internal.hs +11/−0
• Math/Manifold/Embedding/Simple/Class.hs +26/−0
• manifolds.cabal +9/−4
• test/tasty/test.hs +422/−11

@@ -73,20 +73,20 @@   lookupAtlas = id instance Atlas S¹ where   type ChartIndex S¹ = S⁰-  chartReferencePoint NegativeHalfSphere = S¹ $ -pi/2-  chartReferencePoint PositiveHalfSphere = S¹ $ pi/2-  interiorChartReferencePoint _ NegativeHalfSphere = S¹ $ -pi/2-  interiorChartReferencePoint _ PositiveHalfSphere = S¹ $ pi/2-  lookupAtlas (S¹ φ) | φ<0        = NegativeHalfSphere+  chartReferencePoint NegativeHalfSphere = S¹Polar $ -pi/2+  chartReferencePoint PositiveHalfSphere = S¹Polar $ pi/2+  interiorChartReferencePoint _ NegativeHalfSphere = S¹Polar $ -pi/2+  interiorChartReferencePoint _ PositiveHalfSphere = S¹Polar $ pi/2+  lookupAtlas (S¹Polar φ) | φ<0        = NegativeHalfSphere                      | otherwise  = PositiveHalfSphere instance Atlas S² where   type ChartIndex S² = S⁰-  chartReferencePoint PositiveHalfSphere = S² 0 0-  chartReferencePoint NegativeHalfSphere = S² pi 0-  interiorChartReferencePoint _ PositiveHalfSphere = S² 0 0-  interiorChartReferencePoint _ NegativeHalfSphere = S² pi 0-  lookupAtlas (S² ϑ _) | ϑ<pi/2     = PositiveHalfSphere-                       | otherwise  = NegativeHalfSphere+  chartReferencePoint PositiveHalfSphere = S²Polar 0 0+  chartReferencePoint NegativeHalfSphere = S²Polar pi 0+  interiorChartReferencePoint _ PositiveHalfSphere = S²Polar 0 0+  interiorChartReferencePoint _ NegativeHalfSphere = S²Polar pi 0+  lookupAtlas (S²Polar ϑ _) | ϑ<pi/2     = PositiveHalfSphere+                            | otherwise  = NegativeHalfSphere  instance (LinearSpace (a n), Needle (a n) ~ a n, Interior (a n) ~ a n)               => Atlas (LinAff.Point a n) where

@@ -92,7 +92,7 @@ --   be an arbitrary one-dimensional space (i.e. basically real intervals or 'S¹'). --   In these cases, there is always only one partial derivative: that which we --   integrate over, in the only possible direction for propagation.-type ODE x y = DifferentialEqn AffineModel x y+type ODE x y = DifferentialEqn QuadraticModel x y  constLinearDEqn :: ∀ x y . ( SimpleSpace x                            , SimpleSpace y, AffineManifold y@@ -132,7 +132,7 @@     ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ) -> \bwt' ->     let bwt'inv = pseudoInverse bwt'     in \(Shade (_x,y) δxy) -> LocalDifferentialEqn-            (\(AffineModel shy' _) ->+            (\(QuadraticModel shy' _ _) ->                     let shy = dualShade shy'                     in ( return $ shy & shadeNarrowness %~ scaleNorm 0.01                        , return $ projectShade (Embedding (arr bwt')

@@ -0,0 +1,360 @@+-- |+-- Module      : Data.Manifold.FibreBundle+-- Copyright   : (c) Justus Sagemüller 2018+-- License     : GPL v3+-- +-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de+-- Stability   : experimental+-- Portability : portable+-- ++{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE MultiParamTypeClasses      #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE UndecidableInstances       #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE UnicodeSyntax              #-}+{-# LANGUAGE GADTs                      #-}+{-# LANGUAGE DefaultSignatures          #-}+{-# LANGUAGE CPP                        #-}+#if __GLASGOW_HASKELL__ >= 800+{-# LANGUAGE UndecidableSuperClasses    #-}+#endif+++module Data.Manifold.FibreBundle where+++import Data.AdditiveGroup+import Data.VectorSpace+import Math.LinearMap.Category++import Data.Manifold.Types.Primitive+import Data.Manifold.PseudoAffine+    +import qualified Prelude as Hask++import Control.Category.Constrained.Prelude hiding ((^))+import Control.Category.Discrete+import Control.Arrow.Constrained++import Linear.V2 (V2(V2))+import Linear.V3 (V3(V3))++import Data.Tagged+++data TransportOnNeedleWitness k m f where+  TransportOnNeedle :: (ParallelTransporting (LinearFunction (Scalar (Needle m)))+                                             (Needle m) (Needle f))+                     => TransportOnNeedleWitness k m f++data ForgetTransportProperties k m f where+  ForgetTransportProperties :: ParallelTransporting (->) m f+                     => ForgetTransportProperties k m f++class (PseudoAffine m, m ~ Interior m, Category k, Object k f)+           => ParallelTransporting k m f where+  transportOnNeedleWitness :: TransportOnNeedleWitness k m f+  default transportOnNeedleWitness+      :: ParallelTransporting (LinearFunction (Scalar (Needle m))) (Needle m) (Needle f)+           => TransportOnNeedleWitness k m f+  transportOnNeedleWitness = TransportOnNeedle+  forgetTransportProperties :: ForgetTransportProperties k m f+  default forgetTransportProperties :: ParallelTransporting (->) m f+           => ForgetTransportProperties k m f+  forgetTransportProperties = ForgetTransportProperties+  +  parallelTransport :: m -> Needle m -> k f f+  translateAndInvblyParTransport+        :: m -> Needle m -> (m, (k f f, k f f))+  translateAndInvblyParTransport p v+              = (q, ( parallelTransport p v+                    , parallelTransport q $ p.-~!q ))+   where q = p.+~^v++instance ∀ m s . (PseudoAffine m, m ~ Interior m, s ~ (Scalar (Needle m)), Num' s)+      => ParallelTransporting Discrete m (ZeroDim s) where+  transportOnNeedleWitness = case (pseudoAffineWitness :: PseudoAffineWitness m) of+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) -> TransportOnNeedle+  forgetTransportProperties = case (pseudoAffineWitness :: PseudoAffineWitness m) of+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))+        -> ForgetTransportProperties+  parallelTransport _ _ = id+instance ∀ m s . (PseudoAffine m, m ~ Interior m, s ~ (Scalar (Needle m)), Num' s)+      => ParallelTransporting (LinearFunction s) m (ZeroDim s) where+  transportOnNeedleWitness = case (pseudoAffineWitness :: PseudoAffineWitness m) of+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) -> TransportOnNeedle+  forgetTransportProperties = case (pseudoAffineWitness :: PseudoAffineWitness m) of+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))+        -> ForgetTransportProperties+  parallelTransport _ _ = id+instance ∀ m s . (PseudoAffine m, m ~ Interior m, s ~ (Scalar (Needle m)), Num' s)+      => ParallelTransporting (->) m (ZeroDim s) where+  transportOnNeedleWitness = case (pseudoAffineWitness :: PseudoAffineWitness m) of+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) -> TransportOnNeedle+  forgetTransportProperties = case (pseudoAffineWitness :: PseudoAffineWitness m) of+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))+        -> ForgetTransportProperties+  parallelTransport _ _ = id++instance (Category k, Object k ℝ) => ParallelTransporting k ℝ ℝ where+  parallelTransport _ _ = id+instance (Category k, Object k ℝ²) => ParallelTransporting k ℝ² ℝ² where+  parallelTransport _ _ = id+instance (Category k, Object k ℝ³) => ParallelTransporting k ℝ³ ℝ³ where+  parallelTransport _ _ = id+instance (Category k, Object k ℝ⁴) => ParallelTransporting k ℝ⁴ ℝ⁴ where+  parallelTransport _ _ = id++instance (Category k, Object k ℝ) => ParallelTransporting k S¹ ℝ where+  parallelTransport _ _ = id++instance (EnhancedCat k (LinearMap ℝ), Object k ℝ²)+             => ParallelTransporting k S² ℝ² where+  parallelTransport p v = (fst . snd) (translateAndInvblyParTransport p v)+  translateAndInvblyParTransport (S²Polar θ₀ φ₀) 𝐯+     | d < pi     = (S²Polar θ₁ φ₁, (arr fwd, arr bwd))+     | d < 2*pi   = translateAndInvblyParTransport (S²Polar θ₀ φ₀)+                      $ 𝐯^*(-(2*pi-d)/d)+     | otherwise  = translateAndInvblyParTransport (S²Polar θ₀ φ₀)+                      $ let revolutions = floor $ d/(2*pi)+                        in 𝐯^*((d - 2*pi*fromIntegral revolutions)/d)+   where -- See images/constructions/sphericoords-needles.svg. Translation as in+         -- "Data.Manifold.PseudoAffine" instance.+         S¹Polar γc₀ = coEmbed 𝐯+         γ₀ | θ₀ < pi/2   = γc₀ - φ₀+            | otherwise   = γc₀ + φ₀+         d = magnitude 𝐯+         S¹Polar φ₁ = S¹Polar φ₀ .+~^ δφ+         +         -- Cartesian coordinates of p₁ in the system whose north pole is p₀+         -- with φ₀ as the zero meridian+         V3 bx by bz = embed $ S²Polar d γ₀+         +         sθ₀ = sin θ₀; cθ₀ = cos θ₀+         -- Cartesian coordinates of p₁ in the system with the standard north pole,+         -- but still φ₀ as the zero meridian+         (qx,qz) = ( cθ₀ * bx + sθ₀ * bz+                   ,-sθ₀ * bx + cθ₀ * bz )+         qy      = by+         +         S²Polar θ₁ δφ = coEmbed $ V3 qx qy qz+         +         sθ₁ = sin θ₁; cθ₁ = cos θ₁+         +         γ₁+          | sθ₀<=sθ₁  = let+              -- Cartesian coordinates of the standard north pole in the system whose north+              -- pole is p₀ with 𝐯 along the zero meridian+              V3 nbx nby nbz = embed $ S²Polar θ₀ (pi-γ₀)+              +              sd = sin d; cd = cos d+              -- Cartesian coordinates of the standard north pole in the system whose north+              -- pole is p₁ with 𝐯 along the zero meridian+              (ox,oz) = ( cd * nbx - sd * nbz+                        , sd * nbx + cd * nbz )+              oy      = nby++           in atan2 oy (-ox)++          | otherwise = let+              -- Cartesian coordinates of p₀ in the system with the standard north pole,+              -- with p₁ on the zero meridian+              V3 gx gy gz = embed $ S²Polar θ₀ (-δφ)+              +              -- Cartesian coordinates of p₀ in the system whose north+              -- pole is p₁ and the standard north pole on the zero meridian+              (ux,uz) = ( cθ₁ * gx - sθ₁ * gz+                        , sθ₁ * gx + cθ₁ * gz )+              uy      = gy++           in atan2 (-uy) (-ux)++         γc₁ | θ₁ < pi/2  = γ₁ + φ₁+             | otherwise  = γ₁ - φ₁++         (sδγc, cδγc) = sin &&& cos $ γc₁ - γc₀++         fwd = LinearMap (V2 (V2   cδγc  sδγc)+                             (V2 (-sδγc) cδγc)) :: LinearMap ℝ ℝ² ℝ²+         bwd = LinearMap (V2 (V2 cδγc (-sδγc))+                             (V2 sδγc   cδγc )) :: LinearMap ℝ ℝ² ℝ²+++instance {-# OVERLAPS #-} ∀ k a b fa fb s .+         ( ParallelTransporting k a fa, ParallelTransporting k b fb+         , PseudoAffine fa, PseudoAffine fb+         , Scalar (Needle a) ~ s, Scalar (Needle b) ~ s+         , Scalar (Needle fa) ~ s, Scalar (Needle fb) ~ s+         , Num' s+         , Morphism k, ObjectPair k fa fb )+              => ParallelTransporting k (a,b) (fa,fb) where+  transportOnNeedleWitness = case+         ( pseudoAffineWitness :: PseudoAffineWitness a+         , pseudoAffineWitness :: PseudoAffineWitness b+         , pseudoAffineWitness :: PseudoAffineWitness fa+         , pseudoAffineWitness :: PseudoAffineWitness fb+         , transportOnNeedleWitness :: TransportOnNeedleWitness k a fa+         , transportOnNeedleWitness :: TransportOnNeedleWitness k b fb ) of+     ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+      ,TransportOnNeedle, TransportOnNeedle)+         -> TransportOnNeedle+  forgetTransportProperties = case+    ( forgetTransportProperties :: ForgetTransportProperties k a fa+    , forgetTransportProperties :: ForgetTransportProperties k b fb ) of+     (ForgetTransportProperties, ForgetTransportProperties) -> ForgetTransportProperties+  parallelTransport (pa,pb) (va,vb)+       = parallelTransport pa va  *** parallelTransport pb vb++instance ∀ k a f g s .+         ( ParallelTransporting k a f, ParallelTransporting k a g+         , ParallelTransporting (LinearFunction s) (Needle a) (Needle f, Needle g)+         , PseudoAffine f, PseudoAffine g+         , Morphism k, ObjectPair k f g )+              => ParallelTransporting k a (f,g) where+  transportOnNeedleWitness = case+         ( pseudoAffineWitness :: PseudoAffineWitness a+         , pseudoAffineWitness :: PseudoAffineWitness f+         , pseudoAffineWitness :: PseudoAffineWitness g+         , transportOnNeedleWitness :: TransportOnNeedleWitness k a f+         , transportOnNeedleWitness :: TransportOnNeedleWitness k a g ) of+     ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+      ,TransportOnNeedle, TransportOnNeedle)+         -> TransportOnNeedle+  forgetTransportProperties = case+    ( forgetTransportProperties :: ForgetTransportProperties k a f+    , forgetTransportProperties :: ForgetTransportProperties k a g ) of+     (ForgetTransportProperties, ForgetTransportProperties) -> ForgetTransportProperties+  parallelTransport p v+       = parallelTransport p v *** parallelTransport p v+++instance ( ParallelTransporting (LinearFunction (Scalar f)) m f, AdditiveGroup m+         , VectorSpace f )+                => AdditiveGroup (FibreBundle m f) where+  zeroV = FibreBundle zeroV zeroV+  FibreBundle p v ^+^ FibreBundle q w = FibreBundle (p^+^q) (v^+^w)+  negateV (FibreBundle p v) = FibreBundle (negateV p) (negateV v)++instance ∀ m f s .+         ( ParallelTransporting (->) m (Interior f), Semimanifold f+         , ParallelTransporting (LinearFunction s) (Needle m) (Needle f)+         , s ~ Scalar (Needle m) )+                => Semimanifold (FibreBundle m f) where+  type Interior (FibreBundle m f) = FibreBundle m (Interior f)+  type Needle (FibreBundle m f) = FibreBundle (Needle m) (Needle f)+  toInterior (FibreBundle p f) = FibreBundle p <$> toInterior f+  translateP = Tagged $ case ( translateP :: Tagged m (Interior m -> Needle m -> Interior m)+                             , semimanifoldWitness :: SemimanifoldWitness f) of+      (Tagged tpm, SemimanifoldWitness BoundarylessWitness)+           -> \(FibreBundle p f) (FibreBundle v δf)+                   -> FibreBundle (tpm p v) (parallelTransport p v f.+~^δf)+  semimanifoldWitness = case ( semimanifoldWitness :: SemimanifoldWitness m+                             , semimanifoldWitness :: SemimanifoldWitness f+                             , forgetTransportProperties+                               :: ForgetTransportProperties (LinearFunction s) (Needle m) (Needle f)+                             ) of+         (SemimanifoldWitness BoundarylessWitness, SemimanifoldWitness BoundarylessWitness+          ,ForgetTransportProperties)+           -> SemimanifoldWitness BoundarylessWitness+  FibreBundle p f .+~^ FibreBundle v δf+      = FibreBundle (p.+~^v) (parallelTransport p v f.+~^δf)++instance ∀ m f s .+         ( ParallelTransporting (->) m f, ParallelTransporting (->) m (Interior f)+         , PseudoAffine f+         , ParallelTransporting (LinearFunction s) (Needle m) (Needle f)+         , s ~ Scalar (Needle m) )+                => PseudoAffine (FibreBundle m f) where+  pseudoAffineWitness = case ( pseudoAffineWitness :: PseudoAffineWitness m+                             , pseudoAffineWitness :: PseudoAffineWitness f+                             , forgetTransportProperties+                               :: ForgetTransportProperties (LinearFunction s) (Needle m) (Needle f)+                             ) of+     ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+      ,ForgetTransportProperties)+         -> PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+  FibreBundle p f .-~. FibreBundle q g = case p.-~.q of+      Nothing -> Nothing+      Just v  -> FibreBundle v <$> f .-~. parallelTransport p v g+++instance (AdditiveGroup f, x ~ Interior x) => NaturallyEmbedded x (FibreBundle x f) where+  embed x = FibreBundle x zeroV+  coEmbed (FibreBundle x _) = x++instance (NaturallyEmbedded (Interior m) (Interior 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)+           => 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++instance NaturallyEmbedded v w+      => NaturallyEmbedded (FibreBundle ℝ v) (FibreBundle ℝ w) where+  embed (FibreBundle p v) = FibreBundle p $ embed v+  coEmbed (FibreBundle p w) = FibreBundle p $ coEmbed w+instance NaturallyEmbedded v w+      => NaturallyEmbedded (FibreBundle ℝ² v) (FibreBundle ℝ² w) where+  embed (FibreBundle p v) = FibreBundle p $ embed v+  coEmbed (FibreBundle p w) = FibreBundle p $ coEmbed w+instance NaturallyEmbedded v w+      => NaturallyEmbedded (FibreBundle ℝ³ v) (FibreBundle ℝ³ w) where+  embed (FibreBundle p v) = FibreBundle p $ embed v+  coEmbed (FibreBundle p w) = FibreBundle p $ coEmbed w+instance NaturallyEmbedded v w+      => NaturallyEmbedded (FibreBundle ℝ⁴ v) (FibreBundle ℝ⁴ w) where+  embed (FibreBundle p v) = FibreBundle p $ embed v+  coEmbed (FibreBundle p w) = FibreBundle p $ coEmbed w++instance NaturallyEmbedded (FibreBundle S¹ ℝ) (FibreBundle ℝ² ℝ²) where+  embed (FibreBundle (S¹Polar φ) l) = FibreBundle (V2 cφ sφ) $ l*^(V2 (-sφ) cφ)+   where (cφ, sφ) = (cos &&& sin) φ+  coEmbed (FibreBundle (V2 0 0) (V2 _ δy)) = FibreBundle (S¹Polar 0) δy+  coEmbed (FibreBundle p (V2 δx δy)) = FibreBundle (S¹Polar $ atan2 sφ cφ) $ cφ*δy - sφ*δx+   where V2 cφ sφ = p^/r+         r = magnitude p++instance NaturallyEmbedded (FibreBundle S² ℝ²) (FibreBundle ℝ³ ℝ³) where+  embed (FibreBundle (S²Polar θ φ) 𝐯@(V2 δξ δυ))+       = FibreBundle (V3 (sθ*cφ) (sθ*sφ) cθ) 𝐯r+   where [V2 cθ sθ, V2 cφ sφ] = embed . S¹Polar <$> [θ,φ]+         S¹Polar γc = coEmbed 𝐯+         γ | θ < pi/2   = γc - φ+           | otherwise  = γc + φ+         d = magnitude 𝐯++         V2 δθ δφ = d *^ embed (S¹Polar γ)+         +         𝐞φ = V3 (-sφ) cφ 0+         𝐞θ = V3 (cθ*cφ) (cθ*sφ) (-sθ)+         𝐯r = δθ*^𝐞θ ^+^ δφ*^𝐞φ+  +  coEmbed (FibreBundle (V3 x y z) 𝐯r)+           = FibreBundle (S²Polar θ φ) (magnitude (δθ,δφ) *^ embed (S¹Polar γc))+   where r = sqrt $ x^2 + y^2 + z^2+         rxy = sqrt $ x^2 + y^2+         θ = atan2 rxy z+         φ = atan2 y x+         cθ = z / r+         sθ = rxy / r+         (cφ,sφ) | rxy>0      = (x,y)^/rxy+                 | otherwise  = (1,0)+         𝐞φ = V3 (-sφ) cφ 0+         𝐞θ = V3 (cθ*cφ) (cθ*sφ) (-sθ)+         δθ = 𝐞θ <.> 𝐯r+         δφ = 𝐞φ <.> 𝐯r+         γ = atan2 δφ δθ+         γc | θ < pi/2   = γ + φ+            | otherwise  = γ - φ

@@ -0,0 +1,106 @@+-- |+-- Module      : Data.Manifold.Function.Interpolation+-- Copyright   : (c) Justus Sagemüller 2017+-- License     : GPL v3+-- +-- Maintainer  : (@) jsagemue $ uni-koeln.de+-- Stability   : experimental+-- Portability : portable+-- ++{-# LANGUAGE ScopedTypeVariables      #-}+{-# LANGUAGE UnicodeSyntax            #-}+{-# LANGUAGE TypeOperators            #-}+{-# LANGUAGE TupleSections            #-}+{-# LANGUAGE TypeFamilies             #-}+{-# LANGUAGE UndecidableInstances     #-}+{-# LANGUAGE FlexibleContexts         #-}+{-# LANGUAGE StandaloneDeriving       #-}+{-# LANGUAGE TemplateHaskell          #-}+{-# LANGUAGE ConstraintKinds          #-}++module Data.Manifold.Function.Interpolation (+      InterpolationFunction+    ) where+++import Data.Manifold.Types+import Data.Manifold.Types.Primitive ((^))+import Data.Manifold.PseudoAffine+import Data.Manifold.Shade+import Data.Manifold.Web+import Data.Manifold.Web.Internal+import Data.Manifold.Function.LocalModel++import Data.VectorSpace+import Math.LinearMap.Category++import Data.Foldable (toList)+import Data.List.NonEmpty (NonEmpty (..))+import qualified Data.List.NonEmpty as NE++import qualified Prelude as Hask++import Control.Category.Constrained.Prelude+import Control.Arrow.Constrained+import Control.Monad.Constrained++import Control.Lens+import Control.Lens.TH+++newtype InterpolationFunction ㄇ x y = InterpolationFunction {+      _interpWeb :: PointsWeb x (ㄇ x y)+    }+makeLenses ''InterpolationFunction+++fromPointsWeb :: (ModellableRelation x y, LocalModel ㄇ)+                 => PointsWeb x (Shade' y) -> InterpolationFunction ㄇ x y+fromPointsWeb = InterpolationFunction . localFmapWeb (+                 \locInfo -> case fitLocally $+                                    (zeroV, locInfo^.thisNodeData)+                                  : [ (ngbx, ngb^.thisNodeData)+                                    | (ngbx,ngb) <- concat $ localOnion locInfo []] of+                                 Just locModl -> locModl )+++adjustMetricToModel :: ∀ x y ㄇ . (ModellableRelation x y, LocalModel ㄇ)+                 => InterpolationFunction ㄇ x y -> InterpolationFunction ㄇ x y+adjustMetricToModel = _interpWeb >>> webLocalInfo+    >>> \(PointsWeb w) -> InterpolationFunction . PointsWeb $ fmap remetricise w+ where remetricise :: Neighbourhood x (WebLocally x (ㄇ x y))+             -> Neighbourhood x (ㄇ x y)+       remetricise nd = nd & dataAtNode .~ localModel+                           & localScalarProduct .~ newNorm+        where localModel = nd^.dataAtNode.thisNodeData+              newNorm = spanNorm+                  [ dx ^/ ((0.1 + occlusion (ngb^.thisNodeData.tweakLocalOffset)+                                            (fromInterior ySynth))+                           * (dx<.>^δx))+                  | (δx,ngb) <- concat . take 2 $ localOnion (nd^.dataAtNode) []+                  , let dx = nd^.localScalarProduct<$|δx+                        Shade' ySynth _ = evalLocalModel localModel δx ]+                      :: Metric x+++upsampleAtLargeDist :: (ModellableRelation x y, LocalModel ㄇ)+        => ℝ -> (x -> ㄇ x y -> Needle x -> Shade' y)+            -> InterpolationFunction ㄇ x y -> PointsWeb x (Shade' y)+upsampleAtLargeDist dmax gapFillFunc (InterpolationFunction web)+     = fromWebNodes (\(Shade x _) -> case nearestNeighbour webI (fromInterior x) of+                         Just (_,nearest) -> nearest ^. nodeLocalScalarProduct) $ do+          local <- toList webI+          (local^.thisNodeCoord, evalLocalModel (local^.thisNodeData) zeroV) : do +             (ngId, (δx, ngb)) <- local^.nodeNeighbours+             guard (ngId > local^.thisNodeId+                   && (local^.nodeLocalScalarProduct|$|δx) > dmax)+             return ( local^.thisNodeCoord !+~^ δx^/2+                    , gapFillFunc (local^.thisNodeCoord)+                                  (local^.thisNodeData)+                                  (δx^/2) )+ where webI = webLocalInfo web++autoUpsampleAtLargeDist :: (ModellableRelation x y, LocalModel ㄇ)+        => ℝ -> InterpolationFunction ㄇ x y -> PointsWeb x (Shade' y)+autoUpsampleAtLargeDist dmax = upsampleAtLargeDist dmax $ const evalLocalModel

@@ -260,6 +260,7 @@                   => [(Needle x, Shade' y)] -> Maybe (ㄇ x y)   tweakLocalOffset :: ModellableRelation x y                   => Lens' (ㄇ x y) (Shade y)+  evalLocalModel :: ModellableRelation x y => ㄇ x y -> Needle x -> Shade' y  modelParametersOverdetMargin :: Int -> Int modelParametersOverdetMargin n = n + round (sqrt $ fromIntegral n) - 1@@ -296,6 +297,12 @@                      $ (p₀:|ps++[pω])           | otherwise  = Nothing   tweakLocalOffset = affineModelOffset+  evalLocalModel = aEvL pseudoAffineWitness+   where aEvL :: ∀ x y . ModellableRelation x y+                => PseudoAffineWitness y -> AffineModel x y -> Needle x -> Shade' y+         aEvL (PseudoAffineWitness (SemimanifoldWitness _)) (AffineModel shy₀ shj) δx+          = convolveShade' (dualShade shy₀)+                           (dualShade . linearProjectShade (lfun ($ δx)) $ shj)  instance LocalModel QuadraticModel where   fitLocally = qFitL@@ -308,3 +315,13 @@                      $ (p₀:|ps++[pω])           | otherwise  = Nothing   tweakLocalOffset = quadraticModelOffset+  evalLocalModel = aEvL pseudoAffineWitness+   where aEvL :: ∀ x y . ModellableRelation x y+                => PseudoAffineWitness y -> QuadraticModel x y -> Needle x -> Shade' y+         aEvL (PseudoAffineWitness (SemimanifoldWitness _))+              (QuadraticModel shy₀ shj shjj) δx+          = (dualShade shy₀)+           `convolveShade'`+            (dualShade . linearProjectShade (lfun ($ δx)) $ shj)+           `convolveShade'`+            (dualShade . linearProjectShade (lfun ($ squareV δx)) $ shjj)

@@ -55,7 +55,7 @@             , PseudoAffine(..)             -- * Type definitions             -- ** Needles-            , Local(..), (⊙+^)+            , Local(..), (⊙+^), (!+~^)             -- ** Metrics             , Metric, Metric'             , RieMetric, RieMetric'@@ -107,6 +107,7 @@  import Control.Lens (Lens', lens, (^.), (&), (%~), (.~)) +import Data.CallStack (HasCallStack) import GHC.Exts (Constraint)  @@ -385,46 +386,62 @@   fromInterior = id   toInterior = pure   translateP = Tagged (.+~^)-  S² ϑ₀ φ₀ .+~^ δv-     | ϑ₀ < pi/2  = sphereFold PositiveHalfSphere $ ϑ₀*^embed(S¹ φ₀) ^+^ δv-     | otherwise  = sphereFold NegativeHalfSphere $ (pi-ϑ₀)*^embed(S¹ φ₀) ^+^ δv+  S²Polar θ₀ φ₀ .+~^ 𝐯 = S²Polar θ₁ φ₁+   where -- See images/constructions/sphericoords-needles.svg.+         S¹Polar γc = coEmbed 𝐯+         γ | θ₀ < pi/2   = γc - φ₀+           | otherwise   = γc + φ₀+         d = magnitude 𝐯+         S¹Polar φ₁ = S¹Polar φ₀ .+~^ δφ+         +         -- Cartesian coordinates of p₁ in the system whose north pole is p₀+         -- with φ₀ as the zero meridian+         V3 bx by bz = embed $ S²Polar d γ+         +         sθ₀ = sin θ₀; cθ₀ = cos θ₀+         -- Cartesian coordinates of p₁ in the system with the standard north pole,+         -- but still φ₀ as the zero meridian+         (qx,qz) = ( cθ₀ * bx + sθ₀ * bz+                   ,-sθ₀ * bx + cθ₀ * bz )+         qy      = by+         +         S²Polar θ₁ δφ = coEmbed $ V3 qx qy qz+ instance PseudoAffine S² where-  S² ϑ₁ φ₁ .-~. S² ϑ₀ φ₀-     | ϑ₀ < pi/2  = pure ( ϑ₁*^embed(S¹ φ₁) ^-^ ϑ₀*^embed(S¹ φ₀) )-     | otherwise  = pure ( (pi-ϑ₁)*^embed(S¹ φ₁) ^-^ (pi-ϑ₀)*^embed(S¹ φ₀) )+  S²Polar θ₁ φ₁ .-~! S²Polar θ₀ φ₀ = d *^ embed(S¹Polar γc)+   where -- See images/constructions/sphericoords-needles.svg.+         V3 qx qy qz = embed $ S²Polar θ₁ (φ₁-φ₀) -sphereFold :: S⁰ -> ℝ² -> S²-sphereFold hfSphere v-   | ϑ₀ > pi     = S² (inv $ tau - ϑ₀) (toS¹range $ φ₀+pi)-   | otherwise  = S² (inv ϑ₀) φ₀- where S¹ φ₀ = coEmbed v-       ϑ₀ = magnitude v `mod'` tau-       inv ϑ = case hfSphere of PositiveHalfSphere -> ϑ-                                NegativeHalfSphere -> pi - ϑ+         sθ₀ = sin θ₀; cθ₀ = cos θ₀+         (bx,bz) = ( cθ₀ * qx - sθ₀ * qz+                   , sθ₀ * qx + cθ₀ * qz )+         by      = qy +         S²Polar d γ = coEmbed $ V3 bx by bz+         +         γc | θ₀ < pi/2   = γ + φ₀+            | otherwise   = γ - φ₀ +++ instance Semimanifold ℝP² where   type Needle ℝP² = ℝ²   fromInterior = id   toInterior = pure   translateP = Tagged (.+~^)-  ℝP² r₀ φ₀ .+~^ V2 δr δφ-   | r₀ > 1/2   = case r₀ + δr of-                   r₁ | r₁ > 1     -> ℝP² (2-r₁) (toS¹range $ φ₀+δφ+pi)-                      | otherwise  -> ℝP²    r₁  (toS¹range $ φ₀+δφ)-  ℝP² r₀ φ₀ .+~^ δxy = let v = r₀*^embed(S¹ φ₀) ^+^ δxy-                           S¹ φ₁ = coEmbed v-                           r₁ = magnitude v `mod'` 1-                       in ℝP² r₁ φ₁  +  HemisphereℝP²Polar θ₀ φ₀ .+~^ v+      = case S²Polar θ₀ φ₀ .+~^ v of+          S²Polar θ₁ φ₁+           | θ₁ > pi/2   -> HemisphereℝP²Polar (pi-θ₁) (-φ₁)+           | otherwise   -> HemisphereℝP²Polar θ₁        φ₁ instance PseudoAffine ℝP² where-  ℝP² r₁ φ₁ .-~. ℝP² r₀ φ₀-   | r₀ > 1/2   = pure `id` case φ₁-φ₀ of-                          δφ | δφ > 3*pi/2  -> V2 (  r₁ - r₀) (δφ - 2*pi)-                             | δφ < -3*pi/2 -> V2 (  r₁ - r₀) (δφ + 2*pi)-                             | δφ > pi/2    -> V2 (2-r₁ - r₀) (δφ - pi  )-                             | δφ < -pi/2   -> V2 (2-r₁ - r₀) (δφ + pi  )-                             | otherwise    -> V2 (  r₁ - r₀) (δφ       )-   | otherwise  = pure ( r₁*^embed(S¹ φ₁) ^-^ r₀*^embed(S¹ φ₀) )+  HemisphereℝP²Polar θ₁ φ₁ .-~! HemisphereℝP²Polar θ₀ φ₀+      = case S²Polar θ₁ φ₁ .-~! S²Polar θ₀ φ₀ of+          v -> let r² = magnitudeSq v+               in if r²>pi^2/4+                   then S²Polar (pi-θ₁) (-φ₁) .-~! S²Polar θ₀ φ₀+                   else v   -- instance (PseudoAffine m, VectorSpace (Needle m), Scalar (Needle m) ~ ℝ)@@ -460,6 +477,13 @@  type DualNeedleWitness x = DualSpaceWitness (Needle x) +++infixl 6 !+~^+-- | Boundary-unsafe version of `.+~^`.+(!+~^) :: ∀ x . (Semimanifold x, HasCallStack) => x -> Needle x -> x+p!+~^v = case toInterior p of+           Just p' -> p'.+~^v   infix 6 ⊙+^

@@ -164,16 +164,16 @@   geodesicBetween _ _ = empty  instance Geodesic S¹ where-  geodesicBetween (S¹ φ) (S¹ ϕ)-    | abs (φ-ϕ) < pi  = (>>> S¹) <$> geodesicBetween φ ϕ-    | φ > 0           = (>>> S¹ . \ψ -> signum ψ*pi - ψ)+  geodesicBetween (S¹Polar φ) (S¹Polar ϕ)+    | abs (φ-ϕ) < pi  = (>>> S¹Polar) <$> geodesicBetween φ ϕ+    | φ > 0           = (>>> S¹Polar . \ψ -> signum ψ*pi - ψ)                         <$> geodesicBetween (pi-φ) (-ϕ-pi)-    | otherwise       = (>>> S¹ . \ψ -> signum ψ*pi - ψ)+    | otherwise       = (>>> S¹Polar . \ψ -> signum ψ*pi - ψ)                         <$> geodesicBetween (-pi-φ) (pi-ϕ)-  middleBetween (S¹ φ) (S¹ ϕ)-    | abs (φ-ϕ) < pi  = S¹ <$> middleBetween φ ϕ-    | φ > 0           = S¹ <$> middleBetween (pi-φ) (-ϕ-pi)-    | otherwise       = S¹ <$> middleBetween (-pi-φ) (pi-ϕ)+  middleBetween (S¹Polar φ) (S¹Polar ϕ)+    | abs (φ-ϕ) < pi  = S¹Polar <$> middleBetween φ ϕ+    | φ > 0           = S¹Polar <$> middleBetween (pi-φ) (-ϕ-pi)+    | otherwise       = S¹Polar <$> middleBetween (-pi-φ) (pi-ϕ)   -- instance Geodesic (Cℝay S⁰) where

@@ -365,26 +365,39 @@     , (pre, (x,node):post) <- splitAt i lvs               -> Right . fmap (PlainLeaves . (pre++) . (:post) . (x,)) $ f node     | otherwise -> Left $ i-n-treeLeaf i f (DisjointBranches n brs)-    | i<n        = foldl (\case -                             Left i' -> (treeLeaf i' f)-                             result  -> return result-                         ) (Left i) brs-    | otherwise  = Left $ i-n-treeLeaf i f sh@(OverlappingBranches n _ brs)-    | i<n        = foldl (\case -                             Left i' -> (treeLeaf i' f)-                             result  -> return result-                         ) (Left i) (toList brs>>=toList)-    | otherwise  = Left $ i-n+treeLeaf i f (DisjointBranches n _)+    | i>=n   = Left $ i-n+treeLeaf i f (DisjointBranches n (br:|[]))+        = fmap (DisjointBranches n . pure) <$> treeLeaf i f br+treeLeaf i f (DisjointBranches n (br:|br':brs))+        = case treeLeaf i f br of+            Left overshoot -> fmap (\(DisjointBranches _ (br'':|brs'))+                                   -> DisjointBranches n (br:|br'':brs'))+                  <$> treeLeaf overshoot f+                     (DisjointBranches (n-nLeaves br) $ br':|brs)+            Right done -> Right $ DisjointBranches n . (:|br':brs) <$> done+treeLeaf i f (OverlappingBranches n extend (br@(DBranch dir (Hourglass t b)):|brs))+    | i<nt       = fmap (OverlappingBranches n extend+                         . (:|brs) . DBranch dir . (`Hourglass`b))+                    <$> treeLeaf i f t+    | i<nt+nb    = fmap (OverlappingBranches n extend+                         . (:|brs) . DBranch dir . ( Hourglass t))+                    <$> treeLeaf (i-nt) f b+    | br':brs' <- brs+                 = fmap (\(OverlappingBranches _ _ (br'':|brs''))+                         -> OverlappingBranches n extend $ br:|br'':brs'')+                    <$> treeLeaf (i-nt-nb) f (OverlappingBranches n extend $ br':|brs')+    | otherwise  = Left $ i - nt - nb+ where [nt,nb] = nLeaves<$>[t,b] + -- | “Inverse indexing” of a tree. This is roughly a nearest-neighbour search, --   but not guaranteed to give the correct result unless evaluated at the --   precise position of a tree leaf. positionIndex :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x))        => Maybe (Metric x)   -- ^ For deciding (at the lowest level) what “close” means;                              --   this is optional for any tree of depth >1.-        -> x`Shaded`y        -- ^ The tree to index into+        -> (x`Shaded`y)      -- ^ The tree to index into         -> x                 -- ^ Position to look up         -> Maybe (Int, ([x`Shaded`y], (x,y)))                    -- ^ Index of the leaf near to the query point, the “path” of

@@ -28,14 +28,16 @@ {-# LANGUAGE TypeOperators            #-} {-# LANGUAGE ScopedTypeVariables      #-} {-# LANGUAGE UnicodeSyntax            #-}+{-# LANGUAGE PatternSynonyms          #-}   module Data.Manifold.Types (         -- * Index / ASCII names           Real0, Real1, RealPlus, Real2, Real3         , Sphere0, Sphere1, Sphere2-        , Projective1, Projective2+        , Projective0, Projective1, Projective2         , Disk1, Disk2, Cone, OpenCone+        , FibreBundle(..), TangentBundle         -- * Linear manifolds         , ZeroDim(..)         , ℝ, ℝ⁰, ℝ¹, ℝ², ℝ³, ℝ⁴@@ -44,11 +46,11 @@         , Stiefel1(..), stiefel1Project, stiefel1Embed         -- ** Specific examples         , HasUnitSphere(..)-        , S⁰(..), S¹(..), S²(..)+        , S⁰(..), S¹(..), pattern S¹, S²(..), pattern S²         -- * Projective spaces-        , ℝP¹,  ℝP²(..)+        , ℝP⁰(..), ℝP¹(..), pattern ℝP¹,  ℝP²(..), pattern ℝP²         -- * Intervals\/disks\/cones-        , D¹(..), D²(..)+        , D¹(..), D²(..), pattern D²         , ℝay         , CD¹(..), Cℝay(..)         -- * Affine subspaces

@@ -28,23 +28,25 @@ {-# LANGUAGE TypeOperators            #-} {-# LANGUAGE ScopedTypeVariables      #-} {-# LANGUAGE RecordWildCards          #-}+{-# LANGUAGE PatternSynonyms          #-}   module Data.Manifold.Types.Primitive (         -- * Index / ASCII names           Real0, Real1, RealPlus, Real2, Real3         , Sphere0, Sphere1, Sphere2-        , Projective1, Projective2+        , Projective0, Projective1, Projective2         , Disk1, Disk2, Cone, OpenCone+        , FibreBundle(..), TangentBundle         -- * Linear manifolds         , ZeroDim(..)         , ℝ, ℝ⁰, ℝ¹, ℝ², ℝ³, ℝ⁴         -- * Hyperspheres-        , S⁰(..), otherHalfSphere, S¹(..), S²(..)+        , S⁰(..), otherHalfSphere, S¹(..), pattern S¹, S²(..), pattern S²         -- * Projective spaces-        , ℝP¹,  ℝP²(..)+        , ℝP⁰(..), ℝP¹(..), pattern ℝP¹,  ℝP²(..), pattern ℝP²         -- * Intervals\/disks\/cones-        , D¹(..), fromIntv0to1, D²(..)+        , D¹(..), fromIntv0to1, D²(..), pattern D²         , ℝay         , CD¹(..), Cℝay(..)         -- * Tensor products@@ -57,6 +59,7 @@   import Math.Manifold.Core.Types+import Math.Manifold.Core.PseudoAffine (FibreBundle(..), TangentBundle, Interior)  import Data.VectorSpace import Data.VectorSpace.Free@@ -67,6 +70,7 @@ import Data.Basis import Data.Void import Data.Monoid+import Data.Fixed (mod') import Math.LinearMap.Category (type (⊗)())  import Control.Applicative (Const(..), Alternative(..))@@ -80,7 +84,8 @@  import Data.Embedding -+import qualified Test.QuickCheck as QC+import qualified Text.Show.Pragmatic as SP   @@ -98,52 +103,9 @@   --- | The ordinary unit sphere.-data S² = S² { ϑParamS² :: !Double -- ^ Range @[0, π[@.-             , φParamS² :: !Double -- ^ Range @[-π, π[@.-             } deriving (Show)   --- | The two-dimensional real projective space, implemented as a unit disk with---   opposing points on the rim glued together.-data ℝP² = ℝP² { rParamℝP² :: !Double -- ^ Range @[0, 1]@.-               , φParamℝP² :: !Double -- ^ Range @[-π, π[@.-               } deriving (Show)------ | The standard, closed unit disk. Homeomorphic to the cone over 'S¹', but not in the---   the obvious, &#x201c;flat&#x201d; way. (And not at all, despite---   the identical ADT definition, to the projective space 'ℝP²'!)-data D² = D² { rParamD² :: !Double -- ^ Range @[0, 1]@.-             , φParamD² :: !Double -- ^ Range @[-π, π[@.-             } deriving (Show)---- | A (closed) cone over a space @x@ is the product of @x@ with the closed interval 'D¹'---   of &#x201c;heights&#x201d;,---   except on its &#x201c;tip&#x201d;: here, @x@ is smashed to a single point.---   ---   This construct becomes (homeomorphic-to-) an actual geometric cone (and to 'D²') in the---   special case @x = 'S¹'@.-data CD¹ x = CD¹ { hParamCD¹ :: !Double -- ^ Range @[0, 1]@-                 , pParamCD¹ :: !x      -- ^ Irrelevant at @h = 0@.-                 } deriving (Show)----- | An open cone is homeomorphic to a closed cone without the &#x201c;lid&#x201d;,---   i.e. without the &#x201c;last copy&#x201d; of @x@, at the far end of the height---   interval. Since that means the height does not include its supremum, it is actually---   more natural to express it as the entire real ray, hence the name.-data Cℝay x = Cℝay { hParamCℝay :: !Double -- ^ Range @[0, &#x221e;[@-                   , pParamCℝay :: !x      -- ^ Irrelevant at @h = 0@.-                   } deriving (Show)------ class NaturallyEmbedded m v where   embed :: m -> v   coEmbed :: v -> m@@ -156,23 +118,33 @@   embed x = (embed x, zeroV)   coEmbed (x,_) = coEmbed x +instance NaturallyEmbedded ℝ⁰ ℝ⁰ where embed = id; coEmbed = id+instance NaturallyEmbedded ℝ  ℝ  where embed = id; coEmbed = id+instance NaturallyEmbedded ℝ² ℝ² where embed = id; coEmbed = id+instance NaturallyEmbedded ℝ³ ℝ³ where embed = id; coEmbed = id+instance NaturallyEmbedded ℝ⁴ ℝ⁴ where embed = id; coEmbed = id+ instance NaturallyEmbedded S⁰ ℝ where   embed PositiveHalfSphere = 1   embed NegativeHalfSphere = -1   coEmbed x | x>=0       = PositiveHalfSphere             | otherwise  = NegativeHalfSphere instance NaturallyEmbedded S¹ ℝ² where-  embed (S¹ φ) = V2 (cos φ) (sin φ)-  coEmbed (V2 x y) = S¹ $ atan2 y x+  embed (S¹Polar φ) = V2 (cos φ) (sin φ)+  coEmbed (V2 x y) = S¹Polar $ atan2 y x instance NaturallyEmbedded S² ℝ³ where-  embed (S² ϑ φ) = V3 (cos φ * sin ϑ) (sin φ * sin ϑ) (cos ϑ)-  coEmbed (V3 x y z) = S² (acos $ z/r) (atan2 y x)-   where r = sqrt $ x^2 + y^2 + z^2+  embed (S²Polar ϑ φ) = V3 (cos φ * sϑ) (sin φ * sϑ) (cos ϑ)+   where sϑ = sin ϑ+  {-# INLINE embed #-}+  coEmbed (V3 x y z) = S²Polar (atan2 rxy z) (atan2 y x)+   where rxy = sqrt $ x^2 + y^2+  {-# INLINE coEmbed #-}   instance NaturallyEmbedded ℝP² ℝ³ where-  embed (ℝP² r φ) = V3 (r * cos φ) (r * sin φ) (sqrt $ 1-r^2)-  coEmbed (V3 x y z) = ℝP² (sqrt $ 1-(z/r)^2) (atan2 (y/r) (x/r))-   where r = sqrt $ x^2 + y^2 + z^2+  embed (HemisphereℝP²Polar θ φ) = V3 (cθ * cos φ) (cθ * sin φ) (sin θ)+   where cθ = cos θ+  coEmbed (V3 x y z) = HemisphereℝP²Polar (atan2 rxy z) (atan2 y x)+   where rxy = sqrt $ x^2 + y^2  instance NaturallyEmbedded D¹ ℝ where   embed = xParamD¹@@ -210,6 +182,7 @@ type Sphere1 = S¹ type Sphere2 = S² +type Projective0 = ℝP⁰ type Projective1 = ℝP¹ type Projective2 = ℝP² @@ -229,3 +202,48 @@   +instance QC.Arbitrary S⁰ where+  arbitrary = (\hsph -> if hsph then PositiveHalfSphere else NegativeHalfSphere)+               <$> QC.arbitrary+instance SP.Show S⁰ where+  showsPrec = showsPrec++instance QC.Arbitrary S¹ where+  arbitrary = S¹Polar . (pi-) . (`mod'`(2*pi))+               <$> QC.arbitrary+  shrink (S¹Polar φ) = S¹Polar . (pi/12*) <$> QC.shrink (φ*12/pi)+instance SP.Show S¹ where+  showsPrec p (S¹Polar φ) = showParen (p>9) $ ("S¹Polar "++) . SP.showsPrec 10 φ++instance QC.Arbitrary S² where+  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 SP.Show S² where+  showsPrec p (S²Polar θ φ) = showParen (p>9) $ ("S²Polar "++)+                           . SP.showsPrec 10 θ . (' ':) . SP.showsPrec 10 φ++instance QC.Arbitrary ℝP⁰ where+  arbitrary = pure ℝPZero++instance QC.Arbitrary ℝP¹ where+  arbitrary = ( \θ -> HemisphereℝP¹Polar (pi/2 - (θ`mod'`pi)) ) <$> QC.arbitrary+  shrink (HemisphereℝP¹Polar θ) = HemisphereℝP¹Polar . (pi/6*) <$> QC.shrink (θ*6/pi)++instance QC.Arbitrary ℝP² where+  arbitrary = ( \θ φ -> HemisphereℝP²Polar (θ`mod'`pi/2) (pi - (φ`mod'`(2*pi))) )+               <$> QC.arbitrary<*>QC.arbitrary+  shrink (HemisphereℝP²Polar θ φ) = [ HemisphereℝP²Polar (θ'*pi/6) (φ'*pi/12)+                                    | θ' <- QC.shrink (θ*6/pi)+                                    , φ' <- QC.shrink (φ*12/pi) ]+++instance (SP.Show (Interior 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+  arbitrary = FibreBundle <$> QC.arbitrary <*> QC.arbitrary+  shrink (FibreBundle m v) = [ FibreBundle m' v'+                             | m' <- QC.shrink m+                             , v' <- QC.shrink v ]

@@ -49,7 +49,8 @@             , localModels_CGrid               -- * Differential equations               -- ** Fixed resolution-            , iterateFilterDEqn_static, iterateFilterDEqn_static_selective+            , iterateFilterDEqn_static, iterateFilterDEqn_pathwise+            , iterateFilterDEqn_static_selective               -- ** Automatic resolution             , filterDEqnSolutions_adaptive, iterateFilterDEqn_adaptive               -- ** Configuration@@ -59,7 +60,7 @@             , PropagationInconsistency(..)               -- * Misc             , ConvexSet(..), ellipsoid, ellipsoidSet, coerceWebDomain-            , rescanPDELocally, webOnions, knitShortcuts+            , rescanPDELocally, localOnion, webOnions, knitShortcuts             ) where  @@ -818,6 +819,28 @@                            . fmap (shading $->)  +iterateFilterDEqn_pathwise+     :: ( ModellableRelation x y, Hask.MonadPlus m, Hask.Traversable m, LocalModel ㄇ )+       => InformationMergeStrategy [] m (x,Shade' y) iy+           -> Embedding (->) (Shade' y) iy+           -> DifferentialEqn ㄇ x y+                 -> PointsWeb x (Shade' y) -> Cofree m (PointsWeb x (Shade' y))+iterateFilterDEqn_pathwise strategy shading f+            = fmap (fmap (shading >-$))+            . (`evalState`7438)+            . unfoldM (\oldWeb -> do+               r <- get+               let i = r `mod` nLeaves (webNodeRsc oldWeb)+                   m = 2^31 - 1+                   a = 963345    :: Int  -- revised Park-Miller+               put $ (a*r)`mod`m+               return ( oldWeb+                      , filterDEqnSolutions_static strategy shading f+                       =<<filterDEqnSolutions_pathsTowards i strategy shading f oldWeb+                      ))+            . fmap (shading $->)++ iterateFilterDEqn_static_selective :: ( ModellableRelation x y                                       , Hask.MonadPlus m, badness ~ ℝ                                       , LocalModel ㄇ )@@ -869,6 +892,42 @@               _ -> mergeInformation strategy oldValue empty         ) ++filterDEqnSolutions_pathsTowards :: ∀ x y ㄇ iy m .+                     ( ModellableRelation x y, Hask.MonadPlus m, LocalModel ㄇ )+       => WebNodeId+          -> InformationMergeStrategy [] m  (x,Shade' y) iy+          -> Embedding (->) (Shade' y) iy+          -> DifferentialEqn ㄇ x y -> PointsWeb x iy -> m (PointsWeb x iy)+filterDEqnSolutions_pathsTowards = case ( geodesicWitness :: GeodesicWitness y+                                        , boundarylessWitness :: BoundarylessWitness x ) of+   (GeodesicWitness _, BoundarylessWitness) -> \targetNode strategy shading f+       -> traversePathsTowards targetNode+            (\(PathStep stepStart stepEnd) -> StateT $+              \odeState ->+                let apriori = shading >-$ stepEnd^.thisNodeData+                in case propagateDEqnSolution_loc+                                f+                                (LocalDataPropPlan{+                                   _sourcePosition = stepStart^.thisNodeCoord+                                 , _targetPosOffset = (stepEnd^.thisNodeCoord)+                                                        .-~! (stepStart^.thisNodeCoord)+                                 , _sourceData = odeState+                                 , _targetAPrioriData = apriori+                                 , _relatedData+                                     = (fmap (second ((shading>-$) . _thisNodeData))+                                               . concat . tail $ localOnion stepStart+                                                                  [stepEnd^.thisNodeId])+                                 }) of+                          Nothing -> undefined+                              <$> mergeInformation strategy (stepEnd^.thisNodeData) []+                          Just propd -> (, propd)+                                  <$> mergeInformation strategy+                                        (stepEnd^.thisNodeData)+                                        [ ( stepEnd^.thisNodeCoord, apriori )+                                        , ( stepStart^.thisNodeCoord, propd ) ] )+            (\startPoint pathTrav+               -> evalStateT pathTrav $ shading >-$ startPoint^.thisNodeData)   data Average a = Average { weight :: Int

@@ -63,6 +63,7 @@ import Control.Monad (guard, forM_) import Control.Comonad import Control.Monad.Trans.State+import Control.Monad.Trans.Writer  import Control.DeepSeq @@ -511,6 +512,16 @@  type WNIPath = [WebNodeId] type NodeSet = ℤSet.IntSet+++pathsTowards :: ∀ x y . (WithField ℝ Manifold x, HasCallStack)+     => WebNodeId -> PointsWeb x y -> [[y]]+pathsTowards target web = execWriter $ traversePathsTowards+       target+       (\(PathStep _ y) -> tell [y^.thisNodeData] >> return (y^.thisNodeData))+       (\startNode (WriterT (Identity (ν, pathTrav)))+            -> tell [startNode^.thisNodeData : pathTrav] >> return ν)+       web  traversePathInIWeb :: ∀ φ x y . (WithField ℝ Manifold x, Monad φ, HasCallStack)      => [WebNodeId] -> (PathStep x y -> φ y)

@@ -0,0 +1,26 @@+-- |+-- Module      : Math.Manifold.Embedding.Simple.Class+-- Copyright   : (c) Justus Sagemüller 2018+-- License     : GPL v3+-- +-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de+-- Stability   : experimental+-- Portability : portable+-- +-- Some manifolds are “naturally” embedded within some bigger space. For instance,+-- the topological spheres are readily identified with the geometric unit spheres in+-- real vector spaces.+--+-- An embedding is a pretty strong relationship, but often all that's needed is being+-- able to map single points from the manifold to the enclosing space. This module offers+-- a class which does just that.+++++module Math.Manifold.Embedding.Simple.Class (+          NaturallyEmbedded(..)+   ) where+++import Data.Manifold.Types.Primitive

@@ -1,5 +1,5 @@ Name:                manifolds-Version:             0.4.4.0+Version:             0.4.5.0 Category:            Math Synopsis:            Coordinate-free hypersurfaces Description:         Manifolds, a generalisation of the notion of &#x201c;smooth curves&#x201d; or surfaces,@@ -40,10 +40,10 @@  Library   Build-Depends:     base>=4.5 && < 6-                     , manifolds-core == 0.4.4.0+                     , manifolds-core == 0.4.5.0                      , transformers                      , vector-space>=0.8-                     , free-vector-spaces>=0.1.1+                     , free-vector-spaces>=0.1.5                      , linear                      , MemoTrie                      , vector@@ -59,8 +59,9 @@                      , placeholders                      , lens                      , call-stack-                     , constrained-categories >= 0.2.3 && < 0.3.1+                     , constrained-categories >= 0.3.1 && < 0.4                      , pragmatic-show+                     , QuickCheck >=2.8 && <3   other-extensions:  FlexibleInstances                      , TypeFamilies                      , FlexibleContexts@@ -80,6 +81,7 @@                      Data.Manifold.Web.Internal                      Data.Manifold.DifferentialEquation                      Data.Manifold.Function.LocalModel+                     Data.Manifold.Function.Interpolation                      Data.SimplicialComplex                      Data.Function.Differentiable                      Data.Function.Affine@@ -87,7 +89,9 @@                      Data.Manifold.Types.Stiefel                      Data.Manifold.Griddable                      Data.Manifold.Atlas+                     Data.Manifold.FibreBundle                      Data.Manifold.Riemannian+                     Math.Manifold.Embedding.Simple.Class   Other-modules:   Data.List.FastNub                    Data.Manifold.Types.Primitive                    Data.SetLike.Intersection@@ -119,6 +123,7 @@     , pragmatic-show     , containers     , vector-space+    , linear     , constrained-categories     , linearmap-category     , lens

@@ -9,17 +9,23 @@ --   {-# LANGUAGE OverloadedLists, TypeFamilies, FlexibleContexts, UndecidableInstances #-}-{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE FlexibleInstances, AllowAmbiguousTypes  #-}+{-# LANGUAGE TypeOperators, TypeApplications, ScopedTypeVariables, UnicodeSyntax #-}  module Main where  import Data.Manifold.Types import Data.Manifold.PseudoAffine+import Data.Manifold.FibreBundle import Data.Manifold.TreeCover import Data.Manifold.Web import Data.Manifold.Web.Internal import Data.Manifold.Function.LocalModel+import Math.Manifold.Embedding.Simple.Class import Data.VectorSpace+import Data.Cross (cross3)+import Linear.V2 (V2(V2))+import Linear.V3 (V3(V3)) import Math.LinearMap.Category import Prelude hiding (id, fst, snd) import Control.Category.Constrained (id)@@ -35,7 +41,7 @@ import qualified Data.Graph as Graph import qualified Data.Set as Set import Control.Arrow-import Control.Lens+import Control.Lens hiding ((<.>))  import qualified Text.Show.Pragmatic as SP @@ -44,7 +50,262 @@  tests :: TestTree tests = testGroup "Tests"- [ testGroup "Graph structure of webs"+ [ testGroup "Semimanifold laws"+  [ testGroup "Asymptotic associativity"+   [ QC.testProperty "Real vector space" (nearlyAssociative @(ℝ,ℝ))+   , QC.testProperty "1-sphere" (nearlyAssociative @S¹)+   , QC.testProperty "Projective line" (nearlyAssociative @ℝP¹)+   , QC.testProperty "2-sphere" (QC.expectFailure $ nearlyAssociative @S²)+   , QC.testProperty "Projective plane" (QC.expectFailure $ nearlyAssociative @ℝP²)+   ]+  ]+ , testGroup "Pseudoaffine laws"+  [ testGroup "Displacement cancellation"+   [ QC.testProperty "Real vector space" (originCancellation @(ℝ,ℝ))+   , QC.testProperty "1-sphere" (originCancellation @S¹)+   , QC.testProperty "Projective line" (originCancellation @ℝP¹)+   , QC.testProperty "2-sphere" (originCancellation @S²)+   , testGroup "2-sphere corner cases"+    [ QC.testProperty "To north pole"+        $ \(S¹Polar φ) p -> originCancellation (S²Polar 0 φ) p+    , QC.testProperty "From north pole"+        $ \(S¹Polar φ) p -> originCancellation p (S²Polar 0 φ)+    , QC.testProperty "To south pole"+        $ \(S¹Polar φ) p -> originCancellation (S²Polar pi φ) p+    , QC.testProperty "From south pole"+        $ \(S¹Polar φ) p -> originCancellation p (S²Polar pi φ)+    , QC.testProperty "South- to north pole"+        $ \(S¹Polar φ) (S¹Polar ψ) -> originCancellation (S²Polar 0 φ) (S²Polar pi ψ)+    , QC.testProperty "North- to south pole"+        $ \(S¹Polar φ) (S¹Polar ψ) -> originCancellation (S²Polar pi ψ) (S²Polar 0 φ)+    , QC.testProperty "Along equator"+        $ \(S¹Polar φ) (S¹Polar ψ) -> originCancellation (S²Polar (pi/2) ψ) (S²Polar (pi/2) φ)+    , QC.testProperty "Just south of equator"+        $ \(S¹Polar φ) (S¹Polar ψ) -> originCancellation (S²Polar (pi/2 + 1e-10) ψ) (S²Polar (pi/2 + 1e-10) φ)+    , QC.testProperty "Just across the equator"+        $ \(S¹Polar φ) (S¹Polar ψ) -> originCancellation (S²Polar (pi/2) ψ) (S²Polar (pi/2 + 1e-10) φ)+    , QC.testProperty "To equator"+        $ \(S¹Polar φ) p -> originCancellation (S²Polar (pi/2) φ) p+    , QC.testProperty "From equator"+        $ \(S¹Polar φ) p -> originCancellation p (S²Polar (pi/2) φ)+    ]+   , QC.testProperty "Projective plane" (originCancellation @ℝP²)+   ]+  ]+ , testGroup "Natural embeddings"+  [ testGroup "1-sphere"+     [ testCase "North pole" $ embed (S¹Polar $ pi/2) @?≈ (V2 0 1 :: ℝ²)+     , testCase "South pole" $ embed (S¹Polar $ -pi/2) @?≈ (V2 0 (-1) :: ℝ²)+     ]+  , testGroup "2-sphere"+     [ testCase "North pole" $ embed (S²Polar 0 0) @?≈ (V3 0 0 1 :: ℝ³)+     , testCase "South pole" $ embed (S²Polar pi 0) @?≈ (V3 0 0 (-1) :: ℝ³)+     ]+  , testGroup "1-sphere tangent bundle"+     [ testCase "North pole"+           $ embed (FibreBundle (S¹Polar $  pi/2) 1 :: TangentBundle S¹)+               @?≈ (FibreBundle (V2 0 1) (V2 (-1) 0) :: TangentBundle ℝ²)+     , testCase "South pole"+           $ embed (FibreBundle (S¹Polar $ -pi/2) 1 :: TangentBundle S¹)+               @?≈ (FibreBundle (V2 0 (-1)) (V2 1 0) :: TangentBundle ℝ²)+     , testCase "45°"+           $ embed (FibreBundle (S¹Polar $ pi/4) 1 :: TangentBundle S¹)+               @?≈ (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²)+               @?≈ (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²)+               @?≈ (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²)+               @?≈ (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²)+               @?≈ (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²)+               @?≈ (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²)+               @?≈ (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²)+               @?≈ (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²)+               @?≈ (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²)+               @?≈ (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²)+               @?≈ (FibreBundle (V3 1 0 0) (V3 0 0 (-1)) :: TangentBundle ℝ³)+     ]+  ]+ , testGroup "Embedding tangent bundles"+  [ QC.testProperty "Real vector space" (embeddingTangentiality @ℝ² @ℝ² 1)+  , QC.testProperty "1-sphere (unlimited)" (QC.expectFailure+                                       $ embeddingTangentiality @ℝ² @S¹ 1)+  , QC.testProperty "1-sphere" (embeddingTangentiality @ℝ² @S¹ 1e-5)+  , QC.testProperty "2-sphere" (embeddingTangentiality @ℝ³ @S² 1e-5)+  ]+ , testGroup "Embedding back-projection"+  [ QC.testProperty "Real vector space" (embeddingBackProject @(ℝ,ℝ) @ℝ)+  , QC.testProperty "1-sphere" (embeddingBackProject @ℝ² @S¹)+  , QC.testProperty "2-sphere" (embeddingBackProject @ℝ³ @S²)+  , QC.testProperty "Vector space tangent bundle"+       (embeddingBackProject @(TangentBundle (ℝ,ℝ)) @(TangentBundle ℝ) )+  , QC.testProperty "S¹ tangent bundle"+       (embeddingBackProject @(TangentBundle ℝ²) @(TangentBundle S¹) )+  , QC.testProperty "S² tangent bundle"+       (embeddingBackProject @(TangentBundle ℝ³) @(TangentBundle S²) )+  ]+ , 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²)+                    q = p .+~^ v :: S²+                    qCart = embed q :: ℝ³+                    axis = pCart `cross3` qCart+                    FibreBundle _ axisProj :: TangentBundle S²+                        = coEmbed (FibreBundle pCart axis :: TangentBundle ℝ³)+                in vCart <.> axis + 1 ≈ 1    -- i.e. the movement vector is always+                  && v <.> axisProj + 1 ≈ 1  -- orthogonal to the rotation axis.+   ]+  ]+ , testGroup "Parallel transport"+  [ testGroup "Displacement cancellation"+   [ QC.testProperty "Real vector space" (parTransportAssociativity @(ℝ,ℝ))+   , QC.testProperty "1-sphere" (parTransportAssociativity @S¹)+   ]+  , testGroup "Nearby tangent spaces of embedding"+   [ QC.testProperty "Real vector space" (nearbyTangentSpaceEmbedding @(ℝ,ℝ) @ℝ 1)+   , QC.testProperty "1-sphere (unlimited)"+         $ QC.expectFailure (nearbyTangentSpaceEmbedding @ℝ² @S¹ 1)+   , QC.testProperty "1-sphere" (nearbyTangentSpaceEmbedding @ℝ² @S¹ 1e-5)+   , QC.testProperty "2-sphere" (nearbyTangentSpaceEmbedding @ℝ³ @S² 1e-5)+   ]+  , testGroup "2-sphere"+   [ testCase "Non-movement on the equator"+        $ sphereParallelTransportTest+            (S²Polar (pi/2) 0) (S²Polar (pi/2) 0) [V3 0 0 1] [V3 0 0 1]+   , testCase "Micro-movement on the equator"+        $ sphereParallelTransportTest+            (S²Polar (pi/2) 0) (S²Polar (pi/2) 1e-3) [V3 0 0 1] [V3 0 0 1]+   , testCase "Small movement on the equator (ez)"+        $ sphereParallelTransportTest+            (S²Polar (pi/2) 0) (S²Polar (pi/2) (pi/2)) [V3 0 0 1, V3   0  1 0]+                                                       [V3 0 0 1, V3 (-1) 0 0]+   , testCase "Big movement on the equator"+        $ sphereParallelTransportTest+            (S²Polar (pi/2) 0) (S²Polar (pi/2) 3) [V3 0 0 1] [V3 0 0 1]+   , testCase "Big negative movement on the equator"+        $ sphereParallelTransportTest+            (S²Polar (pi/2) 0) (S²Polar (pi/2) (-3)) [V3 0 0 1] [V3 0 0 1]+   , testCase "Movement on the zero meridian from north pole"+        $ sphereParallelTransportTest+            (S²Polar 0 0) (S²Polar (pi/2) 0) [V3 0 1 0] [V3 0 1 0]+   , testCase "Movement on the zero meridian to north pole"+        $ sphereParallelTransportTest+            (S²Polar (pi/2) 0) (S²Polar 0 0) [V3 0 1 0, V3   0  0 1]+                                             [V3 0 1 0, V3 (-1) 0 0]+   , testCase "Crossing the equator on the zero meridian"+        $ sphereParallelTransportTest+            (S²Polar (pi/4) 0) (S²Polar (3*pi/4) 0) [V3 0 1 0, V3 (-1) 0 1] +                                                    [V3 0 1 0, V3   1  0 1]+   , testCase "Crossing the equator on the 90° meridian"+        $ sphereParallelTransportTest+            (S²Polar (pi/4) (pi/2)) (S²Polar (3*pi/4) (pi/2)) [V3 1 0 0, V3 0 (-1) 1]+                                                              [V3 1 0 0, V3 0   1  1]+   , testCase "Crossing the equator on the 180° meridian"+        $ sphereParallelTransportTest+            (S²Polar (pi/4) pi) (S²Polar (3*pi/4) pi) [V3 0 1 0, V3   1  0 1]+                                                      [V3 0 1 0, V3 (-1) 0 1]+   , testCase "Crossing the equator on the -90° meridian"+        $ sphereParallelTransportTest+            (S²Polar (pi/4) (-pi/2)) (S²Polar (3*pi/4) (-pi/2)) [V3 1 0 0, V3 0   1  1]+                                                                [V3 1 0 0, V3 0 (-1) 1]+   , QC.testProperty "Movement on the equator" . QC.expectFailure+        $ \(S¹Polar φ₀) (S¹Polar φ₁) -> assertParTransportNeedleTargetFixpoint+                 (S²Polar 0 0, Just "north pole")+                 (S²Polar (pi/2) φ₀)+                 (S²Polar (pi/2) φ₁)+   , QC.testProperty "Just north of the equator"+        $ \p@(S¹Polar φ₀) q@(S¹Polar φ₁) -> abs (p.-~!q) < 2+            ==> assertParTransportNeedleTargetFixpoint+                 (S²Polar 0 0, Just "north pole")+                 (S²Polar (pi/2-1e-13) φ₀)+                 (S²Polar (pi/2-1e-13) φ₁)+   , QC.testProperty "Just slightly crossing the equator"+        $ \(S¹Polar φ₀) (S¹Polar φ₁) -> assertParTransportNeedleTargetFixpoint+                 (S²Polar 0 0, Just "north pole")+                 (S²Polar (pi/2-1e-13) φ₀)+                 (S²Polar (pi/2+1e-13) φ₁)+   , QC.testProperty "Just south of the equator"+        $ \p@(S¹Polar φ₀) q@(S¹Polar φ₁) -> abs (p.-~!q) < 2+            ==> assertParTransportNeedleTargetFixpoint+                 (S²Polar pi 0, Just "south pole")+                 (S²Polar (pi/2+1e-13) φ₀)+                 (S²Polar (pi/2+1e-13) φ₁)+   , QC.testProperty "Movement on the zero meridian"+        $ \(S¹Polar θ₀) (S¹Polar θ₁) -> assertParTransportNeedleTargetFixpoint+                 (S²Polar (pi/2) (pi/2), Nothing)+                 (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²+                      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 ℝ³+                      pxv = pCart`cross3`vCart+                      qxw = qCart`cross3`wCart+                    in QC.counterexample+                           ("  𝑝 = "++SP.show p++"\t ≃ "++SP.show pCart+                        ++"\n  𝑞 = "++SP.show q++"\t ≃ "++SP.show qCart+                        ++"\n  𝑣 = "++SP.show v++"\t = "++SP.show vCart++" @ 𝑝"+                        ++"\n  𝑤 = "++SP.show w++"\t = "++SP.show wCart++" @ 𝑞"+                        ++"\n𝑝×𝑣 = "++SP.show pxv    -- rotation axis+                        ++"\n𝑞×𝑤 = "++SP.show qxw    -- rotation axis+                             )+                       $ 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 ℝ³+                        infix 7 ×+                        (×) = cross3+                        pxq = pCart×qCart+                        fㄧg = fCart ^-^ gCart+                        ㄍ = magnitudeSq+                    in QC.counterexample+                           ("              𝑝 = "++SP.show p+                        ++"\n              𝑞 = "++SP.show q+                        ++"\n              𝑓 = "++SP.show f+                        ++"\n              𝑔 = "++SP.show g+                        ++"\n            𝑝×𝑞 = "++SP.show pxq  -- rotation axis+                        ++"\n          𝑓 − 𝑔 = "++SP.show fㄧg -- movement in the rot.-plane+                        ++"\n    (𝑝×𝑞)×(𝑓−𝑔) = "++SP.show (pxq × fㄧg)+                        ++"\n    (𝑝×𝑞)·(𝑓−𝑔) = "++SP.show (pxq <.> fㄧg)+                        ++"\n ‖(𝑝×𝑞)×(𝑓−𝑔)‖² = "++SP.show (ㄍ $ pxq × fㄧg)+                        ++"\n         ‖𝑝×𝑞‖² = "++SP.show (ㄍ pxq)+                        ++"\n         ‖𝑓−𝑔‖² = "++SP.show (ㄍ fㄧg)+                        ++"\n  ‖𝑝×𝑞‖²·‖𝑓−𝑔‖² = "++SP.show (ㄍ pxq*ㄍ fㄧg)+                             )+                       $ ㄍ (pxq × fㄧg)      -- Check that 𝑝×𝑞 and 𝑓−𝑔 are orthogonal.+                          ≈ ㄍ pxq * ㄍ fㄧg  -- (For orthogonal 𝐚 and 𝐛, we have+                                              -- ‖𝐚×𝐛‖ = ‖𝐚‖·‖𝐛‖.)+   ]+  ]+ , testGroup "Graph structure of webs"   [ testCase "Manually-defined empty web."     $ toList (fst $ toGraph emptyWeb) @?= []   , testCase "Manually-defined single-point web."@@ -439,33 +700,80 @@  infix 4 ≈ class AEq e where+  fuzzyEq :: ℝ -> e -> e -> Bool+  unitEpsilon :: ℝ+  unitEpsilon = 1e-9   (≈) :: e -> e -> Bool+  (≈) = fuzzyEq (unitEpsilon @e)++instance AEq Double where+  fuzzyEq η x y  = x + abs x*η >= y+          && x - abs x*η <= y instance (SimpleSpace v, Needle v~v, Interior v~v, Floating (Scalar v))              => AEq (Shade' v) where-  Shade' c₀ σ₀ ≈ Shade' c₁ σ₁+  fuzzyEq η (Shade' c₀ σ₀) (Shade' c₁ σ₁)     = (σ₀|$|δ) < ε && (σ₀|$|δ) < ε      && all (is1 . (σ₀|$|)) (normSpanningSystem' σ₁)      && all (is1 . (σ₁|$|)) (normSpanningSystem' σ₀)    where δ = c₁ ^-^ c₀-         ε = 1e-2+         ε = 1e-2 + realToFrac η          is1 x = abs (x-1) < ε instance ( SimpleSpace v, DualVector (Needle' v) ~ v, Interior v ~ v          , InnerSpace (Scalar v), Scalar (Needle' v) ~ Scalar v )               => AEq (Shade v) where-  Shade c₀ σ₀ ≈ Shade c₁ σ₁+  fuzzyEq η (Shade c₀ σ₀) (Shade c₁ σ₁)     = (dualNorm σ₀|$|δ) < ε && (dualNorm σ₀|$|δ) < ε      && all (is1 . (dualNorm σ₀|$|)) (normSpanningSystem σ₁)      && all (is1 . (dualNorm σ₁|$|)) (normSpanningSystem σ₀)    where δ = c₁ ^-^ c₀-         ε = 1e-2+         ε = 1e-2 + realToFrac η          is1 x = abs (x-1) < ε instance AEq a => AEq (Maybe a) where-  Just x ≈ Just y = x ≈ y-  Nothing ≈ Nothing = True-  _ ≈ _ = False+  fuzzyEq η (Just x) (Just y) = fuzzyEq η x y+  fuzzyEq _ Nothing Nothing = True+  fuzzyEq _ _ _ = False instance (AEq (Shade y), AEq (Shade (Needle x +> Needle y)))               => AEq (AffineModel x y) where-  AffineModel b₀ a₀ ≈ AffineModel b₁ a₁ = b₀ ≈ b₁ && a₀ ≈ a₁+  fuzzyEq η (AffineModel b₀ a₀) (AffineModel b₁ a₁) = fuzzyEq η b₀ b₁ && fuzzyEq η a₀ a₁++instance (AEq a, AEq b) => (AEq (a,b)) where+  fuzzyEq η (x,y) (ξ,υ) = fuzzyEq η x ξ && fuzzyEq η y υ+instance AEq S¹ where+  fuzzyEq η (S¹Polar φ) (S¹Polar ϕ)+   | φ > pi/2, ϕ < -pi/2  = fuzzyEq η (S¹Polar $ φ - 2*pi) (S¹Polar ϕ)+   | ϕ > pi/2, φ < -pi/2  = fuzzyEq η (S¹Polar φ) (S¹Polar $ ϕ - 2*pi)+   | otherwise            = abs (φ - ϕ) < η+instance AEq S² where+  fuzzyEq η (S²Polar θ φ) (S²Polar ϑ ϕ)+   | φ > pi/2, ϕ < -pi/2  = fuzzyEq η (S²Polar θ $ φ - 2*pi) (S²Polar ϑ ϕ)+   | ϕ > pi/2, φ < -pi/2  = fuzzyEq η (S²Polar θ φ) (S²Polar ϑ $ ϕ - 2*pi)+   | otherwise            = abs (θ - ϑ) < η && abs (φ - ϕ) * sin θ < η++instance AEq ℝ² where+  fuzzyEq η (V2 x y) (V2 ξ υ) = abs (x - ξ) <= ε && abs (y - υ) <= ε+   where ε = (maximum @[]) (abs<$>[x,y,ξ,υ]) * η+instance AEq ℝ³ where+  fuzzyEq η (V3 x y z) (V3 ξ υ ζ) = (all @[]) ((ε>=) . abs) $ [x-ξ, y-υ, z-ζ]+   where ε = (maximum @[]) (abs<$>[x,y,z,ξ,υ,ζ]) * η++instance AEq ℝP⁰ where+  fuzzyEq _ ℝPZero ℝPZero  = True+instance AEq ℝP¹ where+  fuzzyEq η (HemisphereℝP¹Polar θ) (HemisphereℝP¹Polar ϑ)+   = fuzzyEq η (S¹Polar $ θ*2) (S¹Polar $ ϑ*2)+instance AEq ℝP² where+  fuzzyEq η (HemisphereℝP²Polar θ φ) (HemisphereℝP²Polar ϑ ϕ)+   | φ > pi/2, ϕ < -pi/2  = fuzzyEq η (HemisphereℝP²Polar θ $ φ - 2*pi) (HemisphereℝP²Polar ϑ ϕ)+   | ϕ > pi/2, φ < -pi/2  = fuzzyEq η (HemisphereℝP²Polar θ φ) (HemisphereℝP²Polar ϑ $ ϕ - 2*pi)+   | θ < pi/2             = abs (θ - ϑ) < η && abs (φ - ϕ) * θ < η+   | φ > pi/4, ϕ < -pi/4  = fuzzyEq η (HemisphereℝP²Polar (pi/2) $ φ - pi)+                                      (HemisphereℝP²Polar (pi/2) ϕ)+   | ϕ > pi/4, φ < -pi/4  = fuzzyEq η (HemisphereℝP²Polar (pi/2) φ)+                                      (HemisphereℝP²Polar (pi/2) $ ϕ - pi)+   | otherwise            = abs (φ - ϕ) < η++instance (AEq (Interior m), AEq f) => AEq (FibreBundle m f) where+  fuzzyEq η (FibreBundle p v) (FibreBundle q w) = fuzzyEq η p q && fuzzyEq η v w                                          infix 1 @?≈        (@?≈) :: (AEq e, Show e) => e -> e -> Assertion@@ -473,4 +781,107 @@  | a≈b        = return ()  | otherwise  = assertFailure $ "Expected "++show b++", but got "++show a +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)) +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)++originCancellation :: ∀ m . (AEq m, Manifold m, Show m, Show (Needle m))+                         => m -> m -> QC.Property+originCancellation p q = case ( boundarylessWitness :: BoundarylessWitness m+                              , p.-~.q ) of+      (BoundarylessWitness, Just v)+          -> let p' = q.+~^v+             in QC.counterexample ("v = "++show v++", q+v = "++show p') $ p' ≈ p++embeddingBackProject :: ∀ m n . ( NaturallyEmbedded n m, AEq n, SP.Show m, SP.Show n )+       => n -> QC.Property+embeddingBackProject p = QC.counterexample ("Embedded: "++SP.show ep+                                          ++", back-projected: "++SP.show p')+                           $ p' ≈ p+ where ep = embed p :: m+       p' = coEmbed ep++embeddingTangentiality :: ∀ m n . ( Semimanifold m, Semimanifold 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+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'+ where rvub = realToFrac $ magnitude vub+       v = vub ^* consistRadius+       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++nearbyTangentSpaceEmbedding :: ∀ m n+                     . ( Semimanifold m, Semimanifold n+                       , m ~ Interior m, n ~ Interior n+                       , NaturallyEmbedded n m+                       , NaturallyEmbedded (TangentBundle n) (TangentBundle m)+                       , ParallelTransporting (->) n (Needle n)+                       , SP.Show n, SP.Show (Needle n), AEq (Needle n)+                       , InnerSpace (Needle n), RealFloat (Scalar (Needle n)) )+       => Scalar (Needle n) -> Interior n -> Needle n -> Needle n -> QC.Property+nearbyTangentSpaceEmbedding consistRadius p vub f+         = QC.counterexample ("𝑓 embd. at 𝑝, then proj. at 𝑝+𝑣 = "++SP.show fReProj+                              ++", 𝑓 moved by 𝑣 = "++SP.show g)+            $ fuzzyEq (unitEpsilon @(Needle n) * (1+rvub^2)) g fReProj+ where rvub = realToFrac $ magnitude vub+       v = vub ^* consistRadius+       q :: n+       q = p .+~^ v :: n+       qEmbd = embed q :: m+       FibreBundle _ fReProj :: TangentBundle n+               = coEmbed (FibreBundle qEmbd fEmbd :: TangentBundle m)+       g = parallelTransport p v f+       o :: TangentBundle n+       o = FibreBundle p f+       FibreBundle pEmbd fEmbd = embed o :: TangentBundle m++parTransportAssociativity :: ∀ m+           . ( AEq m, Manifold m, SP.Show m+             , ParallelTransporting (->) m (Needle m)+             , InnerSpace (Needle m), RealFloat (Scalar (Needle m)) )+                         => m -> Needle m -> Needle m -> QC.Property+parTransportAssociativity p v w+ = maximum (map magnitude [v,w]) < 1000+       -- Very vast vectors incur inevitable floating-point uncertainty+  ==> let q, q' :: m+          q = (p .+~^ v) .+~^ parallelTransport p v w+          q' = p .+~^ (v^+^w)+      in QC.counterexample ("(p+v) + 〔pTp. v〕w = "++SP.show q++", p+(v+w) = "++SP.show q')+          $ q ≈ q'++assertParTransportNeedleTargetFixpoint :: ∀ m+     . ( AEq m, Manifold m, SP.Show m, Show (Needle m)+       , ParallelTransporting (->) m (Needle m) )+    => (m, Maybe String) -> m -> m -> QC.Property+assertParTransportNeedleTargetFixpoint (q, qName) p₀ p₁+         = let q'= p₁ .+~^ parallelTransport p₀ (p₁ .-~! p₀) (q .-~! p₀)+           in QC.counterexample+                 ("Should keep pointing on "++qShw++", but got "++ SP.show q')+               $ q' ≈ q+ where qShw = case qName of+        Just s  -> s+        Nothing -> SP.show q+++sphereParallelTransportTest :: S² -> S² -> [ℝ³] -> [ℝ³] -> Assertion+sphereParallelTransportTest p q [] [] = assert True+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²+            | (o,u) <- [(p,v), (q,w)] ]