-- | -- Module : Main -- Copyright : (c) Justus Sagemüller 2017 -- License : GPL v3 -- -- Maintainer : (@) jsag $ hvl.no -- Stability : experimental -- Portability : portable -- {-# LANGUAGE OverloadedLists, TypeFamilies, FlexibleContexts, UndecidableInstances #-} {-# 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 Math.Manifold.Real.Coordinates 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, 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, sort) import qualified Data.Graph as Graph import qualified Data.Set as Set import Control.Arrow import Control.Lens hiding ((<.>)) import Data.Fixed (mod') import qualified Text.Show.Pragmatic as SP main = defaultMain tests tests :: TestTree tests = testGroup "Tests" [ 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 (TangentBundle (S¹Polar $ pi/2) 1) @?≈ (FibreBundle (V2 0 1) (V2 (-1) 0) :: TangentBundle ℝ²) , testCase "South pole" $ embed (TangentBundle (S¹Polar $ -pi/2) 1) @?≈ (FibreBundle (V2 0 (-1)) (V2 1 0) :: TangentBundle ℝ²) , testCase "45°" $ 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 (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 (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 (TangentBundle (S²Polar 0 0) (V2 0 1)) @?≈ (FibreBundle (V3 0 0 1) (V3 0 1 0) :: TangentBundle ℝ³) , testCase "Close to north pole" $ 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 (TangentBundle (S²Polar pi 0) (V2 1 0)) @?≈ (FibreBundle (V3 0 0 (-1)) (V3 (-1) 0 0) :: TangentBundle ℝ³) , testCase "South pole, y-dir" $ 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 (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 (TangentBundle (S²Polar (pi/2) 0) (V2 0 1)) @?≈ (FibreBundle (V3 1 0 0) (V3 0 1 0) :: TangentBundle ℝ³) , testCase "Equator, x-dir" $ 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 (TangentBundle (S²Polar (pi/2) 0) (V2 1 0)) @?≈ (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 -> magnitude v < 1e6 ==> let TangentBundle pCart vCart :: TangentBundle ℝ³ = embed $ TangentBundle p v q = p .+~^ v :: S² qCart = embed q :: ℝ³ axis = pCart `cross3` qCart 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 @(ℝ,ℝ)) , 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 -> magnitude v < 1e6 ==> let q = p .+~^ v :: S² w = parallelTransport p v v vCart :@. pCart = embed (v:@.p) :: TangentBundle ℝ³ wCart :@. qCart = embed (w:@.q) :: 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 :: Needle S² fCart :@. pCart = embed (f :@. p) :: TangentBundle ℝ³ gCart :@. qCart = embed (g :@. q) :: 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 "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) @?= [] , testCase "Manually-defined single-point web." $ toList (fst $ toGraph singletonWeb) @?= [[]] , testCase "Manually-defined simple triangular web." $ toList (fst $ toGraph triangularWeb) @?= [[1,2],[0,2],[0,1]] , testCase "Manually-defined simple quadratic web." $ toList (fst $ toGraph quadraticWeb) @?= [[1,2],[0,3],[0,3],[1,2]] , testCase "Envi-aware traversal over simple quadratic web." $ toList (fst . toGraph $ dummyWebFmap quadraticWeb) @?= [[1,2],[0,3],[0,3],[1,2]] , testCase "Direct neighbours in empty web." $ toList (directNeighbours emptyWeb) @?= [] , testCase "Direct neighbours in single-point web." $ toList (directNeighbours singletonWeb) @?= [[]] , testCase "Direct neighbours in simple triangular web." $ toList (directNeighbours triangularWeb) @?= [[1,2],[0,2],[0,1]] , testCase "Direct neighbours in simple quadratic web." $ toList (directNeighbours quadraticWeb) @?= [[1,2],[0,3],[0,3],[1,2]] , testCase "Direct neighbours in quadratic web with one-direction diagonals." $ toList (directNeighbours unidirDiagonalLinkedWeb) @?= [[1,2,3],[0,3],[0,1,3],[1,2]] , testCase "Direct neighbours in 1-dir diag quadratic web after bidirectionalisation." $ toList (directNeighbours $ bidirectionaliseWebLinks unidirDiagonalLinkedWeb) @?= [[1,2,3],[0,2,3],[0,1,3],[0,1,2]] , testCase "Direct neighbours in unsymmetric web." $ toList (directNeighbours unsymmetricWeb) @?= [[5],[2,3,0],[4,3],[4,2,5,1],[5],[0,1,6],[5],[4,6]] , testCase "Next-neighbours in simple quadratic web." $ toList (nextNeighbours quadraticWeb) @?= [ [(1,[0,3]),(2,[0,3])] , [(0,[1,2]),(3,[1,2])] , [(0,[1,2]),(3,[1,2])] , [(1,[0,3]),(2,[0,3])] ] , testCase "Next-neighbours in triangular web (after scrambling)" $ toList (nextNeighbours $ scrambleKnitting triangularWeb) @?= [ [(2,[1,0]),(1,[2,0])] , [(2,[1,0]),(0,[2,1])] , [(1,[2,0]),(0,[2,1])] ] , testCase "Layers in a nested web" $ toList (pointsLocInEnvi nestedWeb) @?= [ [((1, 朳[(o,朳[ {-LEAF-} (o,朳[]) ])]), 0) ,((2, 朳[(o,朳[ {- {- -} {- -}-} (o,朳[(o,朳[]),(o,朳[])]) ])]), 0) ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 0) ] , [((1, 朳[(o,朳[ (o,朳[]) {-LEAF-} ])]), 1) ,((2, 朳[(o,朳[ {- {- -} {- -}-} (o,朳[(o,朳[]),(o,朳[])]) ])]), 0) ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 0) ] , [((1, 朳[(o,朳[ {-LEAF-} (o,朳[]) ])]), 0) ,((2, 朳[(o,朳[ (o,朳[(o,朳[]),(o,朳[])]) {- {- -} {- -}-} ])]), 2) ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 0) ] , [((1, 朳[(o,朳[ (o,朳[]) {-LEAF-} ])]), 1) ,((2, 朳[(o,朳[ (o,朳[(o,朳[]),(o,朳[])]) {- {- -} {- -}-} ])]), 2) ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 0) ] , [((1, 朳[(o,朳[ {-LEAF-} (o,朳[]) ])]), 0) ,((2, 朳[(o,朳[ {- {- -} {- -}-} (o,朳[(o,朳[]),(o,朳[])]) ])]), 0) ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 4) ] , [((1, 朳[(o,朳[ (o,朳[]) {-LEAF-} ])]), 1) ,((2, 朳[(o,朳[ {- {- -} {- -}-} (o,朳[(o,朳[]),(o,朳[])]) ])]), 0) ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 4) ] , [((1, 朳[(o,朳[ {-LEAF-} (o,朳[]) ])]), 0) ,((2, 朳[(o,朳[ (o,朳[(o,朳[]),(o,朳[])]) {- {- -} {- -}-} ])]), 2) ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 4) ] , [((1, 朳[(o,朳[ (o,朳[]) {-LEAF-} ])]), 1) ,((2, 朳[(o,朳[ (o,朳[(o,朳[]),(o,朳[])]) {- {- -} {- -}-} ])]), 2) ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 4) ] ] , testCase "Next-neighbours in nested web." $ toList (nextNeighbours nestedWeb) @?= [ [ (1,[0,3,4]), (2,[0,3]) ] , [ (0,[1,2]) , (3,[1,6]) , (4,[1,5,6]) ] , [ (0,[1,2]) , (3,[1,6]) ] , [ (1,[0,3,4]), (6,[3,4,7]) ] , [ (1,[0,3,4]), (5,[4,7]) , (6,[3,4,7]) ] , [ (4,[1,5,6]), (7,[5,6]) ] , [ (3,[1,6]) , (4,[1,5,6]), (7,[5,6]) ] , [ (5,[4,7]) , (6,[3,4,7]) ] ] , testCase "Next-neighbours in unsymmetric web." $ toList (nextNeighbours unsymmetricWeb) @?= [ [ (5,[0,1,6]) ] , [ (2,[4,3]) , (3,[4,2,5,1]), (0,[5]) ] , [ (4,[5]) , (3,[4,2,5,1]) ] , [ (4,[5]) , (2,[4,3]) , (5,[0,1,6]), (1,[2,3,0]) ] , [ (5,[0,1,6]) ] , [ (0,[5]) , (1,[2,3,0]) , (6,[5]) ] , [ (5,[0,1,6]) ] , [ (4,[5]) , (6,[5]) ] ] , testCase "Neighbours in unsymmetric web after scrambling." $ toList (directNeighbours $ scrambleKnitting unsymmetricWeb) @?= [ [1,6], [4,3,2,5], [5,4,1], [5,4,0,1,6,2], [0,1,6], [2,3,0], [0,1], [5] ] ] , testGroup "Adjacency layers around points in a web" [ testCase "Onions in nested web" $ toList (webOnions $ localFmapWeb _thisNodeId nestedWeb) @?= [ [[(o,0)],[(o,1),(o,2)],[(o,3),(o,4)],[(o,6),(o,5)],[(o,7)]] , [[(o,1)],[(o,0),(o,3),(o,4)],[(o,6),(o,2),(o,5)],[(o,7)]] , [[(o,2)],[(o,0),(o,3)],[(o,1),(o,6)],[(o,4),(o,7)],[(o,5)]] , [[(o,3)],[(o,1),(o,6)],[(o,4),(o,0),(o,7)],[(o,5),(o,2)]] , [[(o,4)],[(o,1),(o,5),(o,6)],[(o,3),(o,7),(o,0)],[(o,2)]] , [[(o,5)],[(o,4),(o,7)],[(o,6),(o,1)],[(o,3),(o,0)],[(o,2)]] , [[(o,6)],[(o,3),(o,4),(o,7)],[(o,1),(o,5)],[(o,0)],[(o,2)]] , [[(o,7)],[(o,5),(o,6)],[(o,4),(o,3)],[(o,1)],[(o,0)],[(o,2)]] ] ] , testGroup "Neighbour-search for web knitting." [ testCase "1D line of points" $ bestNeighbours (euclideanNorm :: Norm ℝ) (zip [0..] [-1, -0.7 .. 1]) @?= ([3,4], Nothing) , testCase "Origin-boundary excluding two points on the x- and y-axes" $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) [(0, (1,0)), (1, (0,1))] @?= ([0,1], Just (sqrt 2/2, sqrt 2/2)) , testCase "Origin-boundary excluding points in the x≥0 half plane" $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) [(0, (1,0)), (1, (0,1)), (2, (0,-1))] @?= ([0,1,2], Just (1, -1.922877998462862e-16)) , testCase "Best neighbours in a quadratic surrounding" $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) [ (1, (0,-1)), (2, (1,-1)) , (3, (-1,0)), (4, (1,0)) , (5, (-1,1)), (6, (0,1)), (7, (1,1)) ] @?= ([1,3,4,6], Nothing) , testCase "Best neighbours to the corner of a rectangular grid" $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) [ ( 1,(1,0)), ( 2,(2,0)), ( 3,(3,0)) , (10,(0,1)), (11,(1,1)), (12,(2,1)), (13,(3,1)) , (20,(0,2)), (21,(1,2)), (22,(2,2)), (23,(3,2)) ] @?= ([1,10], Just (sqrt 2/2, sqrt 2/2)) , testCase "Best neighbours in a rectangular grid" $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) ((id&&&id) <$> [ (-2,-1), (-1,-1), ( 0,-1), ( 1,-1), ( 2,-1) , (-2, 0), (-1, 0),{-ORIGIN-}( 1, 0), ( 2, 0) , (-2, 1), (-1, 1), ( 0, 1), ( 1, 1), ( 2, 1) ]) @?= ([(0,-1), (-1,0), (1,0), (0,1)], Nothing) , testCase "Best neighbours in a big rectangular grid" $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) ((id&&&id) <$> [ (-3,-3), (-2,-3), (-1,-3), ( 0,-3), ( 1,-3), ( 2,-3), ( 3,-3) , (-3,-2), (-2,-2), (-1,-2), ( 0,-2), ( 1,-2), ( 2,-2), ( 3,-2) , (-3,-1), (-2,-1), (-1,-1), ( 0,-1), ( 1,-1), ( 2,-1), ( 3,-1) , (-3, 0), (-2, 0), (-1, 0),{-ORIGIN-}( 1, 0), ( 2, 0), ( 3, 0) , (-3, 1), (-2, 1), (-1, 1), ( 0, 1), ( 1, 1), ( 2, 1), ( 3, 1) , (-3, 2), (-2, 2), (-1, 2), ( 0, 2), ( 1, 2), ( 2, 2), ( 3, 2) ]) @?= ([(0,-1), (-1,0), (1,0), (0,1)], Nothing) , testCase "Best neighbours in an irregular point-cloud" $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) ((id&&&id) <$> [ (-1,-6) , (0,-5), (4,-5),(5,-5) , (-1,-4) ,(-6,-3), (-4,-3),(2,-3) , (-2,-2), (0,-2) , (1,-1), (4,-1),(5,-1) {-ORIGIN-} ,(-5,1), (-3,1),(-2,1), (2,1), (4,1), (5,1) , (-3,2),(-2,3), (1,3),(2,3) ,(-5,4), (-1,4),(3,4) , (-3,5), (3,5) , (2,6), (5,6),(6,6) ]) @?= ([(1,-1), (-2,-2), (2,1), (-6,-3), (-2,1)], Nothing) , testCase "Best neighbours in degenerate near-boundary constellation" $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) ((id &&& (^-^(3.6, 3.0))) <$> reverse [ (3.15,3.6) , (3.29,3.4) , (3.3,3.2), (3.45,3.2), (3.6,3.2) , (3.15,3.0), (3.3,3.0) {-ORIGIN-} , (3.15,2.8), (3.3,2.8), (3.45,2.8), (3.6,2.8), (3.75,2.8) , (3.45,2.6), (3.6,2.6), (3.75,2.6), (3.9,2.6) , (3.15,2.2) ]) @?= ([(3.6,2.8), (3.3,3.0), (3.6,3.2), (3.75,2.8)], Nothing) , testCase "Best neighbours in point selection from almost-rectangular grid" $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) ([ (235,(0.0,-0.2)) , (248,(-0.7499999999999996,0.0)) , (267,(0.0,0.2)) , (268,(0.15,0.0)) , (271,(-0.14999,0.0)) ]) @?= ([271,267,268,235], Nothing) , testCase "Best neighbours in point selection of 1D web test" $ bestNeighbours (euclideanNorm :: Norm ℝ) ((id &&& (^-^467)) <$> [ 565.5193483520385, 254.62827644949562 , 203.3896874080876, 214.87356399193985 ]) @?= ([565.5193483520385, 254.62827644949562], Nothing) ] , testGroup "Automatically building webs" [ testCase "Minimal, 3-point 1D “web”" $ let web = fromWebNodes euclideanMetric [(x, ()) | x<-[0,1,2]] :: PointsWeb ℝ () in toList (localFmapWeb (\info -> ( fst <$> info^.nodeNeighbours , info^.webBoundingPlane ) ) web) @?= [([1], Just 1), ([0,2], Nothing), ([1], Just $ -1)] , testCase "Linear 1D “web”" $ map sort (toList $ directNeighbours (fromWebNodes euclideanMetric [(x, ()) | x<-[0, 0.1 .. 2]] :: PointsWeb ℝ () )) @?= [ [1,9], [0,2], [1,3], [2,4], [3], [6,12], [5,7], [6,8], [7,9], [0,8], [11,15] , [10,12],[5,11],[14,20],[13,15],[10,14],[17],[16,18],[17,19],[18,20],[13,19] ] , testCase "Small linear 1D web with nonuniform spacing" $ toList (directNeighbours (fromWebNodes euclideanMetric [ (x, ()) | x<-[ 203.3896874080876 , 214.87356399193985 , 254.62827644949562 , 467.0 , 565.5193483520385 ] ] :: PointsWeb ℝ () )) @?= [ [1], [0,2], [1,3], [4,2], [3] ] , 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" $ \ps -> length ps >= 2 ==> length (webBoundary (fromWebNodes euclideanMetric [(x, ()) | x<-Set.toList ps] :: PointsWeb ℝ () ) ) == 2 ] ] , testGroup "Shades" [ testCase "Equality of `Shade`s" $ (1 :± [1]) @?≈ (1 :± [1] :: Shade ℝ) , testCase "Equality of `Shade'`s" $ ((1,0)|±|[(1,-2),(3,4)]) @?≈ ((1,0)|±|[(1,-2),(3,4)] :: Shade' (ℝ,ℝ)) , testCase "Pragmatically showing" $ SP.show ((1,0)|±|[(1,-2),(3,4)] :: Shade' (ℝ,ℝ)) @?= "(1,0)|±|[(5,2),(0,2)]" , testCase "Pragmatically showing (with orthogonal span)" $ SP.show ((1,0)|±|[(6,0),(0,2)] :: Shade' (ℝ,ℝ)) @?= "(1,0)|±|[(6,0),(0,2)]" ] , testGroup "Function models for uncertain data" [ testCase "Fitting a 1D affine model to constant data" $ fitLocally [ (-1, 5|±|[1]), (0, 5|±|[1]), (1, 5|±|[1]) ] @?≈ Just ( AffineModel (5:±[1.15]) (zeroV:±[id^/sqrt 2]) :: AffineModel ℝ ℝ ) , testCase "Fitting a 2D affine model to constant data" $ fitLocally [ ((0,1), 5|±|[1]) , ((-1,0), 5|±|[1]), ((0,0), 5|±|[1]), ((1,0), 5|±|[1]) , ((0,-1), 5|±|[1]) ] @?≈ Just ( AffineModel (5:±[0.9]) (zeroV:±((^/sqrt 2)<$>[fst, snd])) :: AffineModel (ℝ,ℝ) ℝ ) , testCase "Fitting a 1D affine model to rising-uncertainty data" $ fitLocally [ (-1, 3|±|[0.1]), (0, 4|±|[0.5]), (1, 5|±|[1]) ] @?≈ Just ( AffineModel (4:±[1/sqrt 2]) (id:±[id^*0.36]) :: AffineModel ℝ ℝ ) , testCase "Fitting a 1D affine model to quadratic data" $ fitLocally [ (-1, 3|±|[0.1]), (0, 0|±|[0.1]), (1, 3|±|[0.1]) ] @?≈ Just ( AffineModel (2:±[2.94]) (zeroV:±[id^*1.8]) :: AffineModel ℝ ℝ ) ] ] emptyWeb, singletonWeb, triangularWeb, quadraticWeb, nestedWeb, unsymmetricWeb , unidirDiagonalLinkedWeb :: PointsWeb ℝ⁰ () emptyWeb = PointsWeb $ PlainLeaves [] singletonWeb = PointsWeb $ PlainLeaves [ (o, Neighbourhood () mempty euclideanNorm Nothing) ] triangularWeb = PointsWeb $ PlainLeaves [ (o, Neighbourhood () [1,2] euclideanNorm Nothing) , (o, Neighbourhood () [-1,1] euclideanNorm Nothing) , (o, Neighbourhood () [-2,-1] euclideanNorm Nothing) ] quadraticWeb = PointsWeb $ OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass (PlainLeaves [ (o, Neighbourhood () [1,2] euclideanNorm Nothing) , (o, Neighbourhood () [-1,2] euclideanNorm Nothing) ]) (PlainLeaves [ (o, Neighbourhood () [-2,1] euclideanNorm Nothing) , (o, Neighbourhood () [-2,-1] euclideanNorm Nothing) ]) ) nestedWeb = PointsWeb $ OverlappingBranches 8 (Shade o mempty) (pure . DBranch o $ Hourglass (OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass (PlainLeaves [ (o, Neighbourhood () [1,2] euclideanNorm Nothing) , (o, Neighbourhood () [-1,2,3] euclideanNorm Nothing) ]) (PlainLeaves [ (o, Neighbourhood () [-2,1] euclideanNorm Nothing) , (o, Neighbourhood () [-2,3] euclideanNorm Nothing) ]) )) (OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass (PlainLeaves [ (o, Neighbourhood () [-3,1,2] euclideanNorm Nothing) , (o, Neighbourhood () [-1,2] euclideanNorm Nothing) ]) (PlainLeaves [ (o, Neighbourhood () [-3,-2,1] euclideanNorm Nothing) , (o, Neighbourhood () [-2,-1] euclideanNorm Nothing) ]) )) ) unsymmetricWeb = PointsWeb $ OverlappingBranches 8 (Shade o mempty) (pure . DBranch o $ Hourglass (OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass (PlainLeaves [ (o, Neighbourhood () [5] euclideanNorm Nothing) , (o, Neighbourhood () [1,2,-1] euclideanNorm Nothing) ]) (PlainLeaves [ (o, Neighbourhood () [2,1] euclideanNorm Nothing) , (o, Neighbourhood () [1,-1,2,-2] euclideanNorm Nothing) ]) )) (OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass (PlainLeaves [ (o, Neighbourhood () [1] euclideanNorm Nothing) , (o, Neighbourhood () [-5,-4,1] euclideanNorm Nothing) ]) (PlainLeaves [ (o, Neighbourhood () [-1] euclideanNorm Nothing) , (o, Neighbourhood () [-3,-1] euclideanNorm Nothing) ]) )) ) unidirDiagonalLinkedWeb = PointsWeb $ OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass (PlainLeaves [ (o, Neighbourhood () [1,2,3] euclideanNorm Nothing) , (o, Neighbourhood () [-1,2] euclideanNorm Nothing) ]) (PlainLeaves [ (o, Neighbourhood () [-2,-1,1] euclideanNorm Nothing) , (o, Neighbourhood () [-2,-1] euclideanNorm Nothing) ]) ) o = zeroV :: ℝ⁰ dummyWebFmap :: PointsWeb ℝ⁰ a -> PointsWeb ℝ⁰ a dummyWebFmap = localFmapWeb $ \info -> info^.thisNodeData directNeighbours :: WithField ℝ Manifold v => PointsWeb v () -> PointsWeb v [WebNodeId] directNeighbours = localFmapWeb $ \info -> fst <$> info^.nodeNeighbours nextNeighbours :: PointsWeb ℝ⁰ a -> PointsWeb ℝ⁰ [(WebNodeId, [WebNodeId])] nextNeighbours = webLocalInfo >>> localFmapWeb `id` \info -> [ ( nId ≡! nId' ≡! (nInfo^.thisNodeId) ≡! (nInfo'^.thisNodeId) , (fst<$>nInfo^.nodeNeighbours) ≡! (fst<$>nInfo'^.nodeNeighbours) ) | ((nId,(_,nInfo)),(nId',(_,nInfo'))) <- zip (info^.nodeNeighbours) (info^.thisNodeData.nodeNeighbours) , all (==Origin) [ nInfo''^.thisNodeCoord | (_,(_,nInfo''))<-nInfo'^.nodeNeighbours ] ] pointsLocInEnvi :: PointsWeb ℝ⁰ a -> PointsWeb ℝ⁰ [((Int, Trees ℝ⁰), WebNodeId)] pointsLocInEnvi = fmapNodesInEnvi $ \(NodeInWeb (_, orig) env) -> fmap (const $ first ((nLeaves&&&onlyNodes) . fmap (const ())) <$> env) orig scrambleKnitting :: PointsWeb ℝ⁰ a -> PointsWeb ℝ⁰ a scrambleKnitting = tweakWebGeometry euclideanMetric $ \info -> nub [ i' | (_, (_, nInfo)) <- info^.nodeNeighbours , (i',_) <- nInfo^.nodeNeighbours , i' /= info^.thisNodeId ] infixl 4 ≡! (≡!) :: (Eq a, Show a) => a -> a -> a x ≡! y | x==y = x | otherwise = error $ show x++" ≠ "++show y 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 fuzzyEq η (Shade' c₀ σ₀) (Shade' c₁ σ₁) = (σ₀|$|δ) < ε && (σ₀|$|δ) < ε && all (is1 . (σ₀|$|)) (normSpanningSystem' σ₁) && all (is1 . (σ₁|$|)) (normSpanningSystem' σ₀) where δ = c₁ ^-^ c₀ ε = 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 fuzzyEq η (Shade c₀ σ₀) (Shade c₁ σ₁) = (dualNorm σ₀|$|δ) < ε && (dualNorm σ₀|$|δ) < ε && all (is1 . (dualNorm σ₀|$|)) (normSpanningSystem σ₁) && all (is1 . (dualNorm σ₁|$|)) (normSpanningSystem σ₀) where δ = c₁ ^-^ c₀ ε = 1e-2 + realToFrac η is1 x = abs (x-1) < ε instance AEq a => AEq (Maybe a) where 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 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 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 a@?≈b | 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 , 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 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 , 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) -> 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) TangentBundle pEmbd vEmbd = embed (TangentBundle p v) 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 fReProj :@. _= coEmbed (fEmbd :@. qEmbd) :: TangentBundle n g = parallelTransport p v f fEmbd :@. pEmbd = embed (f:@.p) :: 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 [] [] = mempty sphereParallelTransportTest p q (v:vs) (w:ws) = (parallelTransport p (q.-~!p) vSph @?≈ wSph) >> sphereParallelTransportTest p q vs ws 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