manifolds 0.4.1.0 → 0.4.4.0

raw patch · 18 files changed

+2326/−893 lines, 18 filesdep +call-stackdep +manifoldsdep +placeholdersdep ~basedep ~constrained-categoriesdep ~linearmap-categorybinary-addedPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: call-stack, manifolds, placeholders, pragmatic-show, tasty, tasty-hunit, tasty-quickcheck

Dependency ranges changed: base, constrained-categories, linearmap-category, manifolds-core, vector-space

API changes (from Hackage documentation)

Files

• Data/Function/Differentiable.hs +3/−3
• Data/Manifold/DifferentialEquation.hs +28/−58
• Data/Manifold/Function/LocalModel.hs +310/−0
• Data/Manifold/PseudoAffine.hs +29/−5
• Data/Manifold/Shade.hs +73/−189
• Data/Manifold/TreeCover.hs +219/−202
• Data/Manifold/Types.hs +6/−1
• Data/Manifold/Types/Stiefel.hs +2/−0
• Data/Manifold/Web.hs +532/−431
• Data/Manifold/Web/Internal.hs +616/−0
• images/examples/TreesAndWebs/2D-cartesian-strangeaspect.png binary
• images/examples/TreesAndWebs/2D-cartesiandisk.png binary
• images/examples/TreesAndWebs/2D-hexa-honeycomb.png binary
• images/examples/TreesAndWebs/2D-hexagonal.png binary
• images/examples/TreesAndWebs/2D-normaldistrib.png binary
• images/examples/TreesAndWebs/2D-scatter.png binary
• manifolds.cabal +32/−4
• test/tasty/test.hs +476/−0

@@ -676,7 +676,7 @@  where npr (LinearManifoldWitness BoundarylessWitness)            (ClosedScalarWitness :: ClosedScalarWitness s)                   = PreRegion $ ppr . ngt-        where PreRegion ppr = positivePreRegion'+        where PreRegion ppr = positivePreRegion' :: PreRegion s s               ngt = actuallyLinearEndo $ negateV id  preRegionToInfFrom, preRegionFromMinInfTo :: RealDimension s => s -> PreRegion s s@@ -689,14 +689,14 @@  where prif (LinearManifoldWitness BoundarylessWitness)             (ClosedScalarWitness :: ClosedScalarWitness s)             xs = PreRegion $ ppr . trl-        where PreRegion ppr = positivePreRegion'+        where PreRegion ppr = positivePreRegion' :: PreRegion s s               trl = actuallyAffineEndo (-xs) id preRegionFromMinInfTo' = prif (linearManifoldWitness :: LinearManifoldWitness s)                            (closedScalarWitness :: ClosedScalarWitness s)  where prif (LinearManifoldWitness BoundarylessWitness)             (ClosedScalarWitness :: ClosedScalarWitness s)             xe = PreRegion $ ppr . flp-        where PreRegion ppr = positivePreRegion'+        where PreRegion ppr = positivePreRegion' :: PreRegion s s               flp = actuallyAffineEndo xe (negateV id)  intervalPreRegion :: ∀ s . RealDimension s => (s,s) -> PreRegion s s

@@ -36,7 +36,6 @@               DifferentialEqn, ODE             , constLinearDEqn             , constLinearODE-            , constLinearPDE             , iterateFilterDEqn_static             -- * Cost functions for error bounds             , maxDeviationsGoal@@ -59,6 +58,8 @@  import Data.Manifold.Types import Data.Manifold.PseudoAffine+import Data.Manifold.Shade+import Data.Manifold.Function.LocalModel import Data.Function.Differentiable import Data.Function.Differentiable.Data import Data.Manifold.TreeCover@@ -82,6 +83,7 @@ import Data.Foldable.Constrained import Data.Traversable.Constrained (Traversable, traverse) +import Control.Lens  -- | An ordinary differential equation is one that does not need any a-priori --   partial derivatives to compute the derivative for integration in some@@ -90,44 +92,37 @@ --   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 x ℝ⁰ y+type ODE x y = DifferentialEqn AffineModel x y -constLinearDEqn :: ∀ x y ð . ( SimpleSpace x-                             , SimpleSpace y, AffineManifold y-                             , SimpleSpace ð, AffineManifold ð-                             , Scalar x ~ ℝ, Scalar y ~ ℝ, Scalar ð ~ ℝ )-              => ((y,ð) +> (x +> y)) -> ((x +> y) +> (y,ð)) -> DifferentialEqn x ð y+constLinearDEqn :: ∀ x y . ( SimpleSpace x+                           , SimpleSpace y, AffineManifold y+                           , Scalar x ~ ℝ, Scalar y ~ ℝ )+              => (y +> (x +> y)) -> ((x +> y) +> y) -> DifferentialEqn QuadraticModel x y constLinearDEqn = case ( linearManifoldWitness :: LinearManifoldWitness x                        , dualSpaceWitness :: DualSpaceWitness x                        , linearManifoldWitness :: LinearManifoldWitness y-                       , dualSpaceWitness :: DualSpaceWitness y-                       , linearManifoldWitness :: LinearManifoldWitness ð-                       , dualSpaceWitness :: DualSpaceWitness ð ) of+                       , dualSpaceWitness :: DualSpaceWitness y ) of    ( LinearManifoldWitness BoundarylessWitness, DualSpaceWitness-    ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness     ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ) -> \bwt'inv bwt' ->         \(Shade (_x,y) δxy) -> LocalDifferentialEqn-         { _predictDerivatives-            = \(Shade' ð δð) ->-                let j = bwt'inv $ (y,ð)-                    δj = bwt' `transformNorm`-                           sumSubspaceNorms (transformNorm (zeroV&&&id) $ dualNorm δxy) δð-                in return $ Shade' j δj-         , _rescanDerivatives-            = \shy shjApriori _-                -> ( mixShade's $ shy+         { _rescanDifferentialEqn+            = \(QuadraticModel shy' shj'Apriori _) ->+               let shy = dualShade shy'+                   shjApriori = dualShade shj'Apriori+                in ( mixShade's $ shy                              :| [ projectShade-                                   (Embedding (arr bwt'inv <<< id&&&zeroV)-                                              (arr bwt'    >>> fst))+                                   (Embedding (arr bwt'inv)+                                              (arr bwt'))                                    shjApriori ]                    , return $ projectShade-                                   (Embedding (arr bwt'inv <<< zeroV&&&id)-                                              (arr bwt'    >>> snd))-                                   shjApriori+                                   (Embedding (arr bwt')+                                              (arr bwt'inv))+                                   shy                    )          } -constLinearODE :: ∀ x y . ( SimpleSpace x, Scalar x ~ ℝ, SimpleSpace y, Scalar y ~ ℝ )+constLinearODE :: ∀ x y . ( SimpleSpace x, Scalar x ~ ℝ+                          , AffineManifold y, SimpleSpace y, Scalar y ~ ℝ )               => ((x +> y) +> y) -> ODE x y constLinearODE = case ( linearManifoldWitness :: LinearManifoldWitness x                       , dualSpaceWitness :: DualSpaceWitness x@@ -135,40 +130,15 @@                       , dualSpaceWitness :: DualSpaceWitness y ) of    ( LinearManifoldWitness BoundarylessWitness, DualSpaceWitness     ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ) -> \bwt' ->-    let bwt'inv = (bwt'\$)+    let bwt'inv = pseudoInverse bwt'     in \(Shade (_x,y) δxy) -> LocalDifferentialEqn-            (let j = bwt'inv y-                 δj = (bwt'>>>zeroV&&&id) `transformNorm` dualNorm δxy-             in \_ -> return $ Shade' j δj )-            (\shy _ _ -> (pure shy, Just $ Shade' Origin mempty) )+            (\(AffineModel shy' _) ->+                    let shy = dualShade shy'+                    in ( return $ shy & shadeNarrowness %~ scaleNorm 0.01+                       , return $ projectShade (Embedding (arr bwt')+                                                          (arr bwt'inv)) shy )+            ) -constLinearPDE :: ∀ x y ð .-                  ( WithField ℝ SimpleSpace x-                  , WithField ℝ SimpleSpace y-                  , WithField ℝ SimpleSpace ð, AffineManifold ð )-              => ((x +> y) +> ð) -> (ð +> (x +> y)) -> DifferentialEqn x ð y-constLinearPDE = case ( linearManifoldWitness :: LinearManifoldWitness x-                      , dualSpaceWitness :: DualSpaceWitness x-                      , linearManifoldWitness :: LinearManifoldWitness y-                      , dualSpaceWitness :: DualSpaceWitness y-                      , linearManifoldWitness :: LinearManifoldWitness ð-                      , dualSpaceWitness :: DualSpaceWitness ð ) of-   ( LinearManifoldWitness BoundarylessWitness, DualSpaceWitness-    ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness-    ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness )-           -> \bwt' bwt'inv (Shade (_x,y) δxy)-       -> LocalDifferentialEqn-           { _predictDerivatives-              = \(Shade' ð δð) ->-                 let j = bwt'inv $ ð-                     δj = bwt' `transformNorm` δð-                 in return $ Shade' j δj-           , _rescanDerivatives-              = \shy shjApriori _-                -> ( return shy-                   , return $ projectShade (Embedding (arr bwt'inv) (arr bwt')) shjApriori-                   )-           }  -- | A function that variates, relatively speaking, most strongly --   for arguments around 1. In the zero-limit it approaches a constant

@@ -0,0 +1,310 @@+-- |+-- Module      : Data.Manifold.Function.LocalModel+-- 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.LocalModel (+    -- * The model class+      LocalModel (..), ModellableRelation+    -- ** Local data fit models+    , AffineModel(..), QuadraticModel(..)+    , estimateLocalJacobian, estimateLocalHessian+    , propagationCenteredModel+    , propagationCenteredQuadraticModel+    , quadraticModel_derivatives+    -- ** Differential equations+    , DifferentialEqn, LocalDifferentialEqn(..)+    , propagateDEqnSolution_loc, LocalDataPropPlan(..)+    -- ** Range interpolation+    , rangeWithinVertices+    ) where+++import Data.Manifold.Types+import Data.Manifold.PseudoAffine+import Data.Manifold.Types.Primitive ((^))+import Data.Manifold.Shade+import Data.Manifold.Riemannian++import Data.VectorSpace+import Math.LinearMap.Category++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.Lens+import Control.Lens.TH+++newtype LocalDifferentialEqn ㄇ x y = LocalDifferentialEqn {+      _rescanDifferentialEqn :: ㄇ x y+                             -> (Maybe (Shade' y), Maybe (Shade' (LocalLinear x y)))+    }+makeLenses ''LocalDifferentialEqn++type DifferentialEqn ㄇ x y = Shade (x,y) -> LocalDifferentialEqn ㄇ x y++data LocalDataPropPlan x y = LocalDataPropPlan+       { _sourcePosition :: !(Interior x)+       , _targetPosOffset :: !(Needle x)+       , _sourceData, _targetAPrioriData :: !y+       , _relatedData :: [(Needle x, y)]+       }+deriving instance (Show (Interior x), Show y, Show (Needle x))+             => Show (LocalDataPropPlan x y)++makeLenses ''LocalDataPropPlan+++{-# DEPRECATED estimateLocalJacobian "Use `fitLocally`" #-}+estimateLocalJacobian :: ∀ x y . ( WithField ℝ Manifold x, Refinable y+                                 , SimpleSpace (Needle x), SimpleSpace (Needle y) )+            => Metric x -> [(Local x, Shade' y)]+                             -> Maybe (Shade' (LocalLinear x y))+estimateLocalJacobian = elj ( pseudoAffineWitness :: PseudoAffineWitness x+                            , pseudoAffineWitness :: PseudoAffineWitness y )+ where elj ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+           , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) )+        mex [(Local x₁, Shade' y₁ ey₁),(Local x₀, Shade' y₀ ey₀)]+         = return $ Shade' (dx-+|>δy)+                          (Norm . LinearFunction $ \δj -> δx ⊗ (σey<$|δj $ δx))+        where Just δx = x₁.-~.x₀+              δx' = (mex<$|δx)+              dx = δx'^/(δx'<.>^δx)+              Just δy = y₁.-~.y₀+              σey = convolveMetric ([]::[y]) ey₀ ey₁+       elj _ mex (po:ps)+           | DualSpaceWitness <- dualSpaceWitness :: DualNeedleWitness y+           , length ps > 1+               = mixShade's =<< (:|) <$> estimateLocalJacobian mex ps +                             <*> sequenceA [estimateLocalJacobian mex [po,pi] | pi<-ps]+       elj _ _ _ = return $ Shade' zeroV mempty+++data AffineModel x y = AffineModel {+         _affineModelOffset :: Shade                      y+       , _affineModelLCoeff :: Shade ( Needle x  +>Needle y)+       }+deriving instance (Show (Shade y), Show (Shade (Needle x+>Needle y)))+              => Show (AffineModel x y)+makeLenses ''AffineModel+++data QuadraticModel x y = QuadraticModel {+         _quadraticModelOffset :: Shade                      y+       , _quadraticModelLCoeff :: Shade ( Needle x  +>Needle y)+       , _quadraticModelQCoeff :: Shade (Needle x⊗〃+>Needle y)+       }+deriving instance ( Show (Shade y)+                  , Show (Shade (Needle x+>Needle y))+                  , Show (Shade (Needle x⊗〃+>Needle y)) )+              => Show (QuadraticModel x y)+makeLenses ''QuadraticModel++type QModelTup s x y = ( Needle y, (Needle x+>Needle y+                                 , SymmetricTensor s (Needle x)+>(Needle y)) )++++quadratic_linearRegression :: ∀ s x y .+                      ( WithField s PseudoAffine x+                      , WithField s PseudoAffine y, Geodesic y+                      , SimpleSpace (Needle x), SimpleSpace (Needle y) )+            => NE.NonEmpty (Needle x, Shade' y) -> QuadraticModel x y+quadratic_linearRegression = case ( dualSpaceWitness :: DualSpaceWitness (Needle x)+                                  , dualSpaceWitness :: DualSpaceWitness (Needle y) ) of+    (DualSpaceWitness, DualSpaceWitness) -> gLinearRegression+         (\δx -> lfun $ \(c,(b,a)) -> (a $ squareV δx) ^+^ (b $ δx) ^+^ c )+         (\cmy (cBest, (bBest, aBest)) σ+            -> let (σc, (σb, σa)) = second summandSpaceNorms $ summandSpaceNorms σ+               in QuadraticModel (Shade (cmy⊙+^cBest $ ([]::[y])) σc)+                              (Shade bBest σb)+                              (Shade aBest σa) )++gLinearRegression :: ∀ s x y ㄇ ψ.+                      ( WithField s PseudoAffine x+                      , WithField s PseudoAffine y, Geodesic y+                      , SimpleSpace (Needle x), SimpleSpace (Needle y)+                      , SimpleSpace ψ, Scalar ψ ~ s )+            => (Needle x -> ψ -+> Needle y)+               -> (Interior y -> ψ -> Variance ψ -> ㄇ x y)+               -> NE.NonEmpty (Needle x, Shade' y) -> ㄇ x y+gLinearRegression fwdCalc analyse = qlr (pseudoAffineWitness, geodesicWitness)+ where qlr :: (PseudoAffineWitness y, GeodesicWitness y)+                   -> NE.NonEmpty (Needle x, Shade' y) -> ㄇ x y+       qlr (PseudoAffineWitness (SemimanifoldWitness _), GeodesicWitness _) ps+                 = analyse cmy ψ σψ+        where Just cmy = pointsBarycenter $ _shade'Ctr.snd<$>ps+              Just vsxy = Hask.mapM (\(x, Shade' y ey) -> (x,).(,ey)<$>y.-~.cmy) ps+              ψ = linearFit_bestModel regResult+              σψ = dualNorm . (case linearFit_χν² regResult of+                                     χν² | χν² > 0, recip χν² > 0+                                            -> scaleNorm (recip $ 1 + sqrt χν²)+                                     _ -> {-Dbg.trace ("Fit for regression model requires"+               ++" well-defined χν² (which needs positive number of degrees of freedom)."+               ++"\n Data: "++show (length ps+                                * subbasisDimension (entireBasis :: SubBasis (Needle y)))+               ++"\n Model parameters: "++show (subbasisDimension+                                        (entireBasis :: SubBasis ψ)) )-}+                                          id)+                                $ linearFit_modelUncertainty regResult+              regResult = linearRegression (arr . fwdCalc) (NE.toList vsxy)++quadraticModel_derivatives :: ∀ x y .+          ( PseudoAffine x, PseudoAffine y+          , SimpleSpace (Needle x), SimpleSpace (Needle y)+          , Scalar (Needle y) ~ Scalar (Needle x) ) =>+     QuadraticModel x y -> (Shade' y, (Shade' (LocalLinear x y), Shade' (LocalBilinear x y))) +quadraticModel_derivatives (QuadraticModel sh shð shð²)+    | (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))+                                     :: PseudoAffineWitness y <- pseudoAffineWitness+    , DualSpaceWitness :: DualSpaceWitness (Needle x) <- dualSpaceWitness+    , DualSpaceWitness :: DualSpaceWitness (Needle y) <- dualSpaceWitness+             = (dualShade sh, ( dualShade shð+                              , linIsoTransformShade (2*^id) $ dualShade shð² ))++{-# DEPRECATED estimateLocalHessian "Use `fitLocally`" #-}+estimateLocalHessian :: ∀ x y . ( WithField ℝ Manifold x, Refinable y, Geodesic y+                                , FlatSpace (Needle x), FlatSpace (Needle y) )+            => NonEmpty (Local x, Shade' y) -> QuadraticModel x y+estimateLocalHessian pts = quadratic_linearRegression $ first getLocalOffset <$> pts+++propagationCenteredModel :: ∀ x y ㄇ .+                         ( ModellableRelation x y, LocalModel ㄇ )+         => LocalDataPropPlan x (Shade' y) -> ㄇ x y+propagationCenteredModel propPlan = case fitLocally (NE.toList ptsFromCenter) of+                                       Just ㄇ->ㄇ+ where ctrOffset = propPlan^.targetPosOffset^/2+       ptsFromCenter = (negateV ctrOffset, propPlan^.sourceData)+                     :| [(δx^-^ctrOffset, shy)+                        | (δx, shy)+                            <- (propPlan^.targetPosOffset, propPlan^.targetAPrioriData)+                               : propPlan^.relatedData+                        ]+++propagationCenteredQuadraticModel :: ∀ x y .+                         ( ModellableRelation x y )+         => LocalDataPropPlan x (Shade' y) -> QuadraticModel x y+propagationCenteredQuadraticModel = propagationCenteredModel+++propagateDEqnSolution_loc :: ∀ x y ㄇ . (ModellableRelation x y, LocalModel ㄇ)+           => DifferentialEqn ㄇ x y+               -> LocalDataPropPlan x (Shade' y)+               -> Maybe (Shade' y)+propagateDEqnSolution_loc f propPlan+                  = pdesl (dualSpaceWitness :: DualNeedleWitness x)+                          (dualSpaceWitness :: DualNeedleWitness y)+                          (boundarylessWitness :: BoundarylessWitness x)+                          (pseudoAffineWitness :: PseudoAffineWitness y)+                          (geodesicWitness :: GeodesicWitness y)+ where pdesl DualSpaceWitness DualSpaceWitness BoundarylessWitness+             (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))+             (GeodesicWitness _)+          | Nothing <- jacobian  = Nothing+          | otherwise            = pure result+         where (_,jacobian) = f shxy ^. rescanDifferentialEqn+                               $ propagationCenteredModel propPlan+               jacobianSh :: Shade (LocalLinear x y)+               Just jacobianSh = dualShade' <$> jacobian+               mx = propPlan^.sourcePosition .+~^ propPlan^.targetPosOffset ^/ 2 :: x+               (Shade _ expax' :: Shade x)+                    = coverAllAround (propPlan^.sourcePosition)+                                     [δx | (δx,_) <- propPlan^.relatedData]+               shxy = coverAllAround (mx, mυ)+                                     [ (δx ^-^ propPlan^.targetPosOffset ^/ 2, pυ ^+^ v)+                                     | (δx,neυ) <- (zeroV, propPlan^.sourceData)+                                                  : (second id+                                                      <$> propPlan^.relatedData)+                                     , let Just pυ = neυ^.shadeCtr .-~. mυ+                                     , v <- normSpanningSystem' (neυ^.shadeNarrowness)+                                     ]+                where Just mυ = middleBetween (propPlan^.sourceData.shadeCtr)+                                              (propPlan^.targetAPrioriData.shadeCtr)+               expax = dualNorm expax'+               result :: Shade' y+               result = convolveShade' (propPlan^.sourceData)+                             (dualShade . linearProjectShade (lfun ($ δx)) $ jacobianSh)+               δx = propPlan^.targetPosOffset+++type ModellableRelation x y = ( WithField ℝ Manifold x+                              , Refinable y, Geodesic y+                              , FlatSpace (Needle x), FlatSpace (Needle y) )++class LocalModel ㄇ where+  fitLocally :: ModellableRelation x y+                  => [(Needle x, Shade' y)] -> Maybe (ㄇ x y)+  tweakLocalOffset :: ModellableRelation x y+                  => Lens' (ㄇ x y) (Shade y)++modelParametersOverdetMargin :: Int -> Int+modelParametersOverdetMargin n = n + round (sqrt $ fromIntegral n) - 1+++-- | Dimension of the space of affine functions on @v@.+p¹Dimension :: ∀ v p . FiniteDimensional v => p v -> Int+p¹Dimension _ = 1 + d+ where d = subbasisDimension (entireBasis :: SubBasis v)++-- | Dimension of the space of quadratic functions on @v@.+p²Dimension :: ∀ v p . FiniteDimensional v => p v -> Int+p²Dimension _ = 1 + d + (d*(d+1))`div`2+ where d = subbasisDimension (entireBasis :: SubBasis v)++instance LocalModel AffineModel where+  fitLocally = aFitL dualSpaceWitness+   where aFitL :: ∀ x y . ModellableRelation x y+                    => DualSpaceWitness (Needle y)+                      -> [(Needle x, Shade' y)] -> Maybe (AffineModel x y)+         aFitL DualSpaceWitness dataPts+          | (p₀:ps, pω:_) <- splitAt (modelParametersOverdetMargin+                                        $ p¹Dimension ([]::[Needle x])) dataPts+                 = Just . gLinearRegression+                            (\δx -> lfun $ \(b,a) -> (a $ δx) ^+^ b )+                            (\cmy (bBest, aBest) σ+                               -> let (σb, σa) = summandSpaceNorms σ+                                  in AffineModel (Shade (cmy⊙+^bBest $ ([]::[y]))+                                                        $ scaleNorm 2 σb)+                               -- The magic factor 2 seems dubious ↗, but testing indicates+                               -- that this is necessary to not overrate the accuracy.+                               --   TODO:  check the algorithms in linearmap-category.+                                                 (Shade aBest σa) )+                     $ (p₀:|ps++[pω])+          | otherwise  = Nothing+  tweakLocalOffset = affineModelOffset++instance LocalModel QuadraticModel where+  fitLocally = qFitL+   where qFitL :: ∀ x y . ModellableRelation x y+                    => [(Needle x, Shade' y)] -> Maybe (QuadraticModel x y)+         qFitL dataPts+          | (p₀:ps, pω:_) <- splitAt (modelParametersOverdetMargin+                                        $ p²Dimension ([]::[Needle x])) dataPts+                 = Just . quadratic_linearRegression+                     $ (p₀:|ps++[pω])+          | otherwise  = Nothing+  tweakLocalOffset = quadraticModelOffset

@@ -50,12 +50,12 @@  module Data.Manifold.PseudoAffine (             -- * Manifold class-              Manifold+              Manifold(inInterior)             , Semimanifold(..), Needle'             , PseudoAffine(..)             -- * Type definitions             -- ** Needles-            , Local(..)+            , Local(..), (⊙+^)             -- ** Metrics             , Metric, Metric'             , RieMetric, RieMetric'@@ -97,7 +97,7 @@  import Data.CoNat -import qualified Prelude+import qualified Prelude as Hask import qualified Control.Applicative as Hask  import Control.Category.Constrained.Prelude hiding ((^))@@ -118,6 +118,9 @@   boundarylessWitness :: BoundarylessWitness m   default boundarylessWitness :: (m ~ Interior m) => BoundarylessWitness m   boundarylessWitness = BoundarylessWitness+  inInterior :: m -> Interior m+  default inInterior :: (m ~ Interior m) => m -> Interior m+  inInterior = id instance (PseudoAffine m, LSpace (Needle m), Interior m ~ m) => Manifold m  @@ -135,8 +138,22 @@   --   ≡ locallyTrivialDiffeomorphism p .+~^ 'coerceNeedle' v   -- @   locallyTrivialDiffeomorphism :: x -> ξ-  coerceNeedle :: Functor p (->) (->) => p (x,ξ) -> (Needle x -+> Needle ξ)-  coerceNeedle' :: Functor p (->) (->) => p (x,ξ) -> (Needle' x -+> Needle' ξ)+  coerceNeedle :: Hask.Functor p => p (x,ξ) -> (Needle x -+> Needle ξ)+  coerceNeedle' :: Hask.Functor p => p (x,ξ) -> (Needle' x -+> Needle' ξ)+  coerceNorm :: Hask.Functor p => p (x,ξ) -> Metric x -> Metric ξ+  coerceNorm p = case ( oppositeLocalCoercion :: CanonicalDiffeomorphism ξ x+                      , dualSpaceWitness :: DualSpaceWitness (Needle x)+                      , dualSpaceWitness :: DualSpaceWitness (Needle ξ) ) of+    (CanonicalDiffeomorphism, DualSpaceWitness, DualSpaceWitness)+          -> case ( coerceNeedle (swap<$>p), coerceNeedle' p ) of+              (f, f') -> \(Norm n) -> Norm $ f' . n . f+  coerceVariance :: Hask.Functor p => p (x,ξ) -> Metric' x -> Metric' ξ+  coerceVariance p = case ( oppositeLocalCoercion :: CanonicalDiffeomorphism ξ x+                          , dualSpaceWitness :: DualSpaceWitness (Needle x)+                          , dualSpaceWitness :: DualSpaceWitness (Needle ξ) ) of+    (CanonicalDiffeomorphism, DualSpaceWitness, DualSpaceWitness)+          -> case ( coerceNeedle p, coerceNeedle' (swap<$>p) ) of+              (f, f') -> \(Norm n) -> Norm $ f . n . f'   oppositeLocalCoercion :: CanonicalDiffeomorphism ξ x   default oppositeLocalCoercion :: LocallyCoercible ξ x => CanonicalDiffeomorphism ξ x   oppositeLocalCoercion = CanonicalDiffeomorphism@@ -443,3 +460,10 @@  type DualNeedleWitness x = DualSpaceWitness (Needle x) +++infix 6 ⊙+^+-- | Proxy-version of `translateP`.+(⊙+^) :: ∀ x proxy . Semimanifold x => Interior x -> Needle x -> proxy x -> Interior x+(⊙+^) x v _ = tp x v+ where Tagged tp = translateP :: Tagged x (Interior x -> Needle x -> Interior x)

@@ -44,17 +44,10 @@        , factoriseShade, orthoShades, (✠), intersectShade's, linIsoTransformShade        , embedShade, projectShade        , Refinable, subShade', refineShade', convolveShade', coerceShade-       , mixShade's, dualShade-       -- * Misc-       -- ** Shades+       , mixShade's, dualShade, dualShade', wellDefinedShade', linearProjectShade        , shadesMerge, pointsShades', pseudoECM, convolveMetric        , WithAny(..), shadeWithAny, shadeWithoutAnything-       -- ** Local data fit models-       , estimateLocalJacobian, estimateLocalHessian, QuadraticModel(..)-       -- ** Differential equations-       , DifferentialEqn, LocalDifferentialEqn(..)-       , propagateDEqnSolution_loc, LocalDataPropPlan(..)-       -- ** Range interpolation+       -- * Misc        , rangeOnGeodesic, rangeWithinVertices     ) where @@ -99,6 +92,7 @@ import GHC.Generics (Generic)  import Text.Show.Number+import qualified Text.Show.Pragmatic as SP   -- | A 'Shade' is a very crude description of a region within a manifold. It@@ -108,8 +102,10 @@ --  --   For a /precise/ description of an arbitrarily-shaped connected subset of a manifold, --   there is 'Region', whose implementation is vastly more complex.-data Shade x = Shade { _shadeCtr :: !(Interior x)-                     , _shadeExpanse :: !(Metric' x) }+data Shade x where+   Shade :: (Semimanifold x, SimpleSpace (Needle x))+           =>  { _shadeCtr :: !(Interior x)+               , _shadeExpanse :: !(Metric' x) } -> Shade x deriving instance (Show (Interior x), Show (Metric' x), WithField ℝ PseudoAffine x)                 => Show (Shade x) @@ -119,28 +115,7 @@ data Shade' x = Shade' { _shade'Ctr :: !(Interior x)                        , _shade'Narrowness :: !(Metric x) } -data LocalDifferentialEqn x ð y = LocalDifferentialEqn {-      _predictDerivatives :: Shade' ð -> Maybe (Shade' (LocalLinear x y))-    , _rescanDerivatives :: Shade' y -> Shade' (LocalLinear x y)-                             -> Shade' (LocalBilinear x y)-                             -> (Maybe (Shade' y), Maybe (Shade' ð))-    }-makeLenses ''LocalDifferentialEqn -type DifferentialEqn x ð y = Shade (x,y) -> LocalDifferentialEqn x ð y--data LocalDataPropPlan x ym yr = LocalDataPropPlan-       { _sourcePosition :: !(Interior x)-       , _targetPosOffset :: !(Needle x)-       , _sourceData, _targetAPrioriData :: !ym-       , _relatedData :: [(Needle x, yr)]-       }-deriving instance (Show (Interior x), Show ym, Show yr, Show (Needle x))-             => Show (LocalDataPropPlan x ym yr)--makeLenses ''LocalDataPropPlan-- class IsShade shade where --  type (*) shade :: *->*   -- | Access the center of a 'Shade' or a 'Shade''.@@ -156,7 +131,8 @@                     , PseudoAffine y, SimpleSpace (Needle y)                     , Scalar (Needle x) ~ Scalar (Needle y) )                 => shade (x,y) -> (shade x, shade y)-  coerceShade :: (Manifold x, Manifold y, LocallyCoercible x y) => shade x -> shade y+  coerceShade :: ( Manifold x, Manifold y, LocallyCoercible x y+                 , SimpleSpace (Needle y) ) => shade x -> shade y   -- | ASCII version of '✠'.   orthoShades :: ( PseudoAffine x, SimpleSpace (Needle x)            , PseudoAffine y, SimpleSpace (Needle y)@@ -168,7 +144,7 @@   -- | Squash a shade down into a lower dimensional space.   projectShade :: ( Semimanifold x, Semimanifold y                   , Object (Affine s) (Interior x), Object (Affine s) (Interior y)-                  , SemiInner (Needle x), SemiInner (Needle y) )+                  , SimpleSpace (Needle x), SemiInner (Needle y) )                         => Embedding (Affine s) (Interior x) (Interior y)                               -> shade y -> shade x   -- | Include a shade in a higher-dimensional space. Notice that this behaves@@ -177,13 +153,13 @@   --   pillar” pointing in the projection's orthogonal complement.   embedShade :: ( Semimanifold x, Semimanifold y                 , Object (Affine s) (Interior x), Object (Affine s) (Interior y)-                , SemiInner (Needle x), SemiInner (Needle y) )+                , SemiInner (Needle x), SimpleSpace (Needle y) )                         => Embedding (Affine s) (Interior x) (Interior y)                               -> shade x -> shade y     linearProjectShade :: ∀ s x y-          . (Num' s, LinearSpace x, LinearSpace y, Scalar x ~ s, Scalar y ~ s)+          . (Num' s, LinearSpace x, SimpleSpace y, Scalar x ~ s, Scalar y ~ s)                   => (x+>y) -> Shade x -> Shade y linearProjectShade = case ( linearManifoldWitness :: LinearManifoldWitness x                           , linearManifoldWitness :: LinearManifoldWitness y@@ -235,7 +211,7 @@          fs DualSpaceWitness DualSpaceWitness (Shade x δx) (Shade y δy)              = Shade (x,y) $ sumSubspaceNorms δx δy   coerceShade = cS dualSpaceWitness dualSpaceWitness-   where cS :: ∀ x y . (LocallyCoercible x y)+   where cS :: ∀ x y . (LocallyCoercible x y, SimpleSpace (Needle y))                 => DualNeedleWitness x -> DualNeedleWitness y -> Shade x -> Shade y          cS DualSpaceWitness DualSpaceWitness                     = \(Shade x δxym) -> Shade (internCoerce x) (tN δxym)@@ -246,7 +222,7 @@                       CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism   linIsoTransformShade = lits linearManifoldWitness linearManifoldWitness                               dualSpaceWitness dualSpaceWitness-   where lits :: ∀ x y . ( LinearSpace x, LinearSpace y+   where lits :: ∀ x y . ( LinearSpace x, SimpleSpace y                          , Scalar x ~ Scalar y, Num' (Scalar x) )                => LinearManifoldWitness x -> LinearManifoldWitness y                    -> DualSpaceWitness x -> DualSpaceWitness y@@ -257,8 +233,9 @@               f (Shade x δx)                   = Shade (f $ x) (transformNorm (adjoint $ f) δx)   embedShade = ps' (semimanifoldWitness, semimanifoldWitness)-   where ps' :: ∀ s x y . ( Object (Affine s) (Interior x), Object (Affine s) (Interior y)-                          , SemiInner (Needle x), SemiInner (Needle y) )+   where ps' :: ∀ s x y . ( Semimanifold y+                          , Object (Affine s) (Interior x), Object (Affine s) (Interior y)+                          , SemiInner (Needle x), SimpleSpace (Needle y) )                         => (SemimanifoldWitness x, SemimanifoldWitness y)                -> Embedding (Affine s) (Interior x) (Interior y)                               -> Shade x -> Shade y@@ -267,8 +244,9 @@           where y = q $ x                 (_,j) = evalAffine q x   projectShade = ps' (semimanifoldWitness, semimanifoldWitness)-   where ps' :: ∀ s x y . ( Object (Affine s) (Interior x), Object (Affine s) (Interior y)-                          , SemiInner (Needle x), SemiInner (Needle y) )+   where ps' :: ∀ s x y . ( Semimanifold x+                          , Object (Affine s) (Interior x), Object (Affine s) (Interior y)+                          , SimpleSpace (Needle x), SemiInner (Needle y) )                         => (SemimanifoldWitness x, SemimanifoldWitness y)                -> Embedding (Affine s) (Interior x) (Interior y)                               -> Shade y -> Shade x@@ -412,10 +390,11 @@          Just pinterp = case geodesicWitness :: GeodesicWitness x of             GeodesicWitness _ -> geodesicBetween c ζ -fullShade :: WithField ℝ PseudoAffine x => Interior x -> Metric' x -> Shade x+fullShade :: (Semimanifold x, SimpleSpace (Needle x))+                      => Interior x -> Metric' x -> Shade x fullShade ctr expa = Shade ctr expa -fullShade' :: WithField ℝ PseudoAffine x => Interior x -> Metric x -> Shade' x+fullShade' :: WithField ℝ SimpleSpace x => Interior x -> Metric x -> Shade' x fullShade' ctr expa = Shade' ctr expa  @@ -425,11 +404,11 @@ #if GLASGOW_HASKELL < 800 pattern (:±) :: () #else-pattern (:±) :: (WithField ℝ Manifold x, SimpleSpace (Needle x))+pattern (:±) :: (Semimanifold x, SimpleSpace (Needle x)) #endif-             => (WithField ℝ Manifold x, SimpleSpace (Needle x))+             => (Semimanifold x, SimpleSpace (Needle x))                          => Interior x -> [Needle x] -> Shade x-pattern x :± shs <- Shade x (varianceSpanningSystem -> shs)+pattern x :± shs <- (Shade x (varianceSpanningSystem -> shs))  where x :± shs = fullShade x $ spanVariance shs  -- | Similar to ':±', but instead of expanding the shade, each vector /restricts/ it.@@ -459,7 +438,7 @@ --   Hence the result type is a list. pointsShades :: (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))                                  => [Interior x] -> [Shade x]-pointsShades = map snd . pointsShades' mempty . map fromInterior+pointsShades = map snd . pointsShades' mempty . map ((,()) . fromInterior)  coverAllAround :: ∀ x s . ( Fractional' s, WithField s PseudoAffine x                           , SimpleSpace (Needle x) )@@ -488,10 +467,11 @@ pointsCovers = case pseudoAffineWitness :: PseudoAffineWitness x of                  (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) ->                   \ps -> map (\(ps', Shade x₀ _)-                                -> coverAllAround x₀ [v | p<-ps'+                                -> coverAllAround x₀ [v | (p,())<-ps'                                                         , let Just v                                                                  = p.-~.fromInterior x₀])-                             (pointsShades' mempty (fromInterior<$>ps) :: [([x], Shade x)])+                             (pointsShades' mempty ((,()).fromInterior<$>ps)+                                  :: [([(x,())], Shade x)])  pointsShade's :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))                      => [Interior x] -> [Shade' x]@@ -503,20 +483,20 @@ pointsCover's = case dualSpaceWitness :: DualNeedleWitness x of  DualSpaceWitness -> map (\(Shade c e :: Shade x) -> Shade' c $ dualNorm e) . pointsCovers -pseudoECM :: ∀ x p . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x), Hask.Functor p)-                => p x -> NonEmpty x -> (x, ([x],[x]))+pseudoECM :: ∀ x y p . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x), Hask.Functor p)+                => p x -> NonEmpty (x,y) -> (x, ([(x,y)],[(x,y)])) pseudoECM = case semimanifoldWitness :: SemimanifoldWitness x of  SemimanifoldWitness _ ->-   \_ (p₀ NE.:| psr) -> foldl' ( \(acc, (rb,nr)) (i,p)+   \_ ((p₀,y₀) NE.:| psr) -> foldl' ( \(acc, (rb,nr)) (i,(p,y))                                 -> case (p.-~.acc, toInterior acc) of                                        (Just δ, Just acci)-                                        -> (acci .+~^ δ^/i, (p:rb, nr))-                                      _ -> (acc, (rb, p:nr)) )+                                        -> (acci .+~^ δ^/i, ((p,y):rb, nr))+                                      _ -> (acc, (rb, (p,y):nr)) )                              (p₀, mempty)-                             ( zip [1..] $ p₀:psr )+                             ( zip [1..] $ (p₀,y₀):psr ) -pointsShades' :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                                => Metric' x -> [x] -> [([x], Shade x)]+pointsShades' :: ∀ x y . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))+                                => Metric' x -> [(x,y)] -> [([(x,y)], Shade x)] pointsShades' _ [] = [] pointsShades' minExt ps = case (expa, toInterior ctr) of                             (Just e, Just c)@@ -525,7 +505,7 @@                                   ++ pointsShades' minExt unreachable  where (ctr,(inc'd,unreachable)) = pseudoECM ([]::[x]) $ NE.fromList ps        expa = ( (<>minExt) . spanVariance . map (^/ fromIntegral (length ps)) )-              <$> mapM (.-~.ctr) ps+              <$> mapM ((.-~.ctr) . fst) ps          -- | Attempt to reduce the number of shades to fewer (ideally, a single one).@@ -931,143 +911,14 @@ intersectShade's (sh:|shs) = Hask.foldrM refineShade' sh shs  -estimateLocalJacobian :: ∀ x y . ( WithField ℝ Manifold x, Refinable y-                                 , SimpleSpace (Needle x), SimpleSpace (Needle y) )-            => Metric x -> [(Local x, Shade' y)]-                             -> Maybe (Shade' (LocalLinear x y))-estimateLocalJacobian = elj ( pseudoAffineWitness :: PseudoAffineWitness x-                            , pseudoAffineWitness :: PseudoAffineWitness y )- where elj ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)-           , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) )-        mex [(Local x₁, Shade' y₁ ey₁),(Local x₀, Shade' y₀ ey₀)]-         = return $ Shade' (dx-+|>δy)-                          (Norm . LinearFunction $ \δj -> δx ⊗ (σey<$|δj $ δx))-        where Just δx = x₁.-~.x₀-              δx' = (mex<$|δx)-              dx = δx'^/(δx'<.>^δx)-              Just δy = y₁.-~.y₀-              σey = convolveMetric ([]::[y]) ey₀ ey₁-       elj _ mex (po:ps)-           | DualSpaceWitness <- dualSpaceWitness :: DualNeedleWitness y-           , length ps > 1-               = mixShade's =<< (:|) <$> estimateLocalJacobian mex ps -                             <*> sequenceA [estimateLocalJacobian mex [po,pi] | pi<-ps]-       elj _ _ _ = return $ Shade' zeroV mempty   -data QuadraticModel x y = QuadraticModel {-         _quadraticModelOffset :: Interior y-       , _quadraticModel :: Shade (Needle y, (Needle x+>Needle y, Needle x⊗〃+>Needle y))-       } -quadratic_linearRegression :: ∀ s x y .-                      ( WithField s PseudoAffine x-                      , WithField s PseudoAffine y, Geodesic y-                      , SimpleSpace (Needle x), SimpleSpace (Needle y) )-            => NE.NonEmpty (Needle x, Shade' y) -> QuadraticModel x y-quadratic_linearRegression = qlr-                  ( dualSpaceWitness, pseudoAffineWitness-                  , linearManifoldWitness, dualSpaceWitness-                  , geodesicWitness )- where qlr :: ( DualSpaceWitness (Needle x)-              , PseudoAffineWitness y, LinearManifoldWitness (Needle y)-              , DualSpaceWitness (Needle y)-              , GeodesicWitness y )-                   -> NE.NonEmpty (Needle x, Shade' y) -> QuadraticModel x y-       qlr ( DualSpaceWitness-           , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)-           , LinearManifoldWitness BoundarylessWitness, DualSpaceWitness-           , GeodesicWitness _ ) ps-                 = QuadraticModel cmy-                     $ coverAllAround mBest (convexPolytopeRepresentatives dm)-        where Just cmy = pointsBarycenter $ _shade'Ctr.snd<$>ps-              Just vsxy = Hask.mapM (\(x, Shade' y ey) -> (x,).(,ey)<$>y.-~.cmy) ps-              (mBest :: ( Needle y, (Needle x+>Needle y-                              , SymmetricTensor s (Needle x)+>(Needle y))-                            )-               , dm)-                        = linearRegressionWVar-                           (\δx -> lfun $ \(c,(b,a)) -> (a $ squareV δx)-                                                      ^+^ (b $ δx) ^+^ c )-                           (NE.toList vsxy) -estimateLocalHessian :: ∀ x y . ( WithField ℝ Manifold x, Refinable y, Geodesic y-                                , FlatSpace (Needle x), FlatSpace (Needle y) )-            => NonEmpty (Local x, Shade' y) -> QuadraticModel x y-estimateLocalHessian pts = elj ( pseudoAffineWitness :: PseudoAffineWitness x-                               , pseudoAffineWitness :: PseudoAffineWitness y )- where elj ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)-           , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) )-         = theModel-        where localPts :: NonEmpty (Needle x, Shade' y)-              localPts = pts >>= \(Local x, Shade' y ey)-                             -> NE.fromList [ (x, Shade' (y.+~^σ*^δy) ey)-                                            | δy <- normSpanningSystem' ey-                                            , σ <- [-1,1] ]-              theModel = quadratic_linearRegression localPts   -propagateDEqnSolution_loc :: ∀ x y ð . ( WithField ℝ Manifold x-                                       , Refinable y, Geodesic (Interior y)-                                       , WithField ℝ AffineManifold ð, Geodesic ð-                                       , SimpleSpace (Needle x), SimpleSpace (Needle ð) )-           => DifferentialEqn x ð y-               -> LocalDataPropPlan x (Shade' y, Shade' ð) (Shade' y)-               -> Maybe (Shade' y)-propagateDEqnSolution_loc f propPlan-                  = pdesl (dualSpaceWitness :: DualNeedleWitness x)-                          (dualSpaceWitness :: DualNeedleWitness y)-                          (boundarylessWitness :: BoundarylessWitness x)-                          (pseudoAffineWitness :: PseudoAffineWitness y)- where pdesl DualSpaceWitness DualSpaceWitness BoundarylessWitness-             (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))-          | Nothing <- jacobian  = Nothing-          | otherwise            = pure result-         where jacobian = (f shxy ^. predictDerivatives $ shð)-                           >>= \j -> mixShade's $ j:|[aprioriDirDrv]-               Just (Shade' j₀ jExpa) = jacobian-               jacobianSh :: Shade (LocalLinear x y)-               Just jacobianSh = dualShade' <$> jacobian-               aprioriDirDrv :: Shade' (LocalLinear x y)-               Just aprioriDirDrv = estimateLocalJacobian expax-                                 [ (Local zeroV :: Local x, propPlan^.sourceData._1)-                                 , (Local δx,        propPlan^.targetAPrioriData._1) ]-               mx = propPlan^.sourcePosition .+~^ propPlan^.targetPosOffset ^/ 2 :: x-               Just shð = middleBetween (propPlan^.sourceData._2)-                                        (propPlan^.targetAPrioriData._2)-               shxy = coverAllAround (mx, mυ)-                                     [ (δx ^-^ propPlan^.targetPosOffset ^/ 2, pυ ^+^ v)-                                     | (δx,neυ) <- (zeroV, propPlan^.sourceData._1)-                                                  : (second id-                                                      <$> propPlan^.relatedData)-                                     , let Just pυ = neυ^.shadeCtr .-~. mυ-                                     , v <- normSpanningSystem' (neυ^.shadeNarrowness)-                                     ]-                where Just mυ = middleBetween (propPlan^.sourceData._1.shadeCtr)-                                              (propPlan^.targetAPrioriData._1.shadeCtr)-               (Shade _ expax' :: Shade x)-                    = coverAllAround (propPlan^.sourcePosition)-                                     [δx | (δx,_) <- propPlan^.relatedData]-               expax = dualNorm expax'-               result :: Shade' y-               Just result = wellDefinedShade' $ convolveShade'-                        (case wellDefinedShade' $ propPlan^.sourceData._1 of {Just s->s})-                        (case wellDefinedShade' . dualShade-                               . linearProjectShade (lfun ($ δx))-                                $ jacobianSh-                           of {Just s->s})-                where δyb = j₀ $ δx-               δx = propPlan^.targetPosOffset-------- -- | Essentially the same as @(x,y)@, but not considered as a product topology. --   The 'Semimanifold' etc. instances just copy the topology of @x@, ignoring @y@. data x`WithAny`y@@ -1123,7 +974,7 @@ shadeWithAny :: y -> Shade x -> Shade (x`WithAny`y) shadeWithAny y (Shade x xe) = Shade (WithAny y x) xe -shadeWithoutAnything :: Shade (x`WithAny`y) -> Shade x+shadeWithoutAnything :: Semimanifold x => Shade (x`WithAny`y) -> Shade x shadeWithoutAnything (Shade (WithAny _ b) e) = Shade b e                        @@ -1140,6 +991,10 @@ prettyShowShade' :: LtdErrorShow x => Shade' x -> String prettyShowShade' sh = prettyShowsPrecShade' 0 sh [] +instance LtdErrorShow x => SP.Show (Shade' x) where+  showsPrec = prettyShowsPrecShade'+instance LtdErrorShow x => SP.Show (Shade x) where+  showsPrec = prettyShowsPrecShade   wellDefinedShade' :: LinearSpace (Needle x) => Shade' x -> Maybe (Shade' x)@@ -1157,6 +1012,16 @@                          => LtdErrorShowWitness m   ltdErrorShowWitness = LtdErrorShowWitness pseudoAffineWitness   showsPrecShade'_errorLtdC :: Int -> Shade' m -> ShowS+  prettyShowsPrecShade :: Int -> Shade m -> ShowS+  prettyShowsPrecShade p sh@(Shade c e')+              = showParen (p>6) $ v+                   . (":±["++) . flip (foldr id) (intersperse (',':) u) . (']':)+   where v = showsPrecShade'_errorLtdC 6 (Shade' c e :: Shade' m)+         u :: [ShowS] = case ltdErrorShowWitness :: LtdErrorShowWitness m of+           LtdErrorShowWitness (PseudoAffineWitness (SemimanifoldWitness _)) ->+             [ showsPrecShade'_errorLtdC 6 (Shade' δ e :: Shade' (Needle m))+             | δ <- varianceSpanningSystem e']+         e = dualNorm' e'   prettyShowsPrecShade' :: Int -> Shade' m -> ShowS   prettyShowsPrecShade' p sh@(Shade' c e)               = showParen (p>6) $ v@@ -1213,6 +1078,25 @@    where (shx,shy) = factoriseShade sh          shshx = showsPrecShade'_errorLtdC 0 shx           shshy = showsPrecShade'_errorLtdC 0 shy ++instance ∀ v .+    (HilbertSpace v, SemiInner v, FiniteDimensional v, LtdErrorShow v, Scalar v ~ ℝ)+              => LtdErrorShow (LinearMap ℝ v ℝ) where+  showsPrecShade'_errorLtdC p sh = showParen (p>7) $+         ("().<"++) . showsPrecShade'_errorLtdC 7+                        (linIsoTransformShade (arr fromLinearForm) sh :: Shade' v)+instance ∀ v .+    (HilbertSpace v, SemiInner v, FiniteDimensional v, LtdErrorShow v, Scalar v ~ ℝ)+              => LtdErrorShow (LinearMap ℝ v (ℝ,ℝ)) where+  showsPrecShade'_errorLtdC p sh = showParen (p>7) $+         (   "Left ().<"++) . showsPrecShade'_errorLtdC 7 shx+       . ("^+^Right().<"++) . showsPrecShade'_errorLtdC 7 shy+   where (shx,shy) = factoriseShade+                        (linIsoTransformShade (lfun $ \f+                                                -> ( fromLinearForm $ fst . f+                                                   , fromLinearForm $ snd . f ) ) sh+                             :: Shade' (v,v))+                                 instance LtdErrorShow x => Show (Shade' x) where   showsPrec = prettyShowsPrecShade'

@@ -45,22 +45,21 @@        , Refinable, subShade', refineShade', convolveShade', coerceShade        , mixShade's        -- * Shade trees-       , ShadeTree(..), fromLeafPoints, onlyLeaves, indexShadeTree, positionIndex-       -- * View helpers-       , onlyNodes+       , ShadeTree, fromLeafPoints, fromLeafPoints_, onlyLeaves, onlyLeaves_+       , indexShadeTree, treeLeaf, positionIndex+       -- ** View helpers+       , entireTree, onlyNodes, trunkBranches, nLeaves, treeDepth        -- ** Auxiliary types-       , SimpleTree, Trees, NonEmptyTree, GenericTree(..)+       , SimpleTree, Trees, NonEmptyTree, GenericTree(..), 朳        -- * Misc-       , HasFlatView(..), shadesMerge, smoothInterpolate+       , HasFlatView(..), shadesMerge        , allTwigs, twigsWithEnvirons, Twig, TwigEnviron, seekPotentialNeighbours-       , completeTopShading, flexTwigsShading, coerceShadeTree-       , WithAny(..), Shaded, fmapShaded, joinShaded-       , constShaded, zipTreeWithList, stripShadedUntopological+       , completeTopShading, flexTwigsShading, traverseTrunkBranchChoices+       , Shaded(..), fmapShaded+       , constShaded, zipTreeWithList        , stiAsIntervalMapping, spanShading-       , estimateLocalJacobian-       , DifferentialEqn, LocalDifferentialEqn(..)-       , propagateDEqnSolution_loc, LocalDataPropPlan(..)-       , rangeOnGeodesic+       , DBranch, DBranch'(..), Hourglass(..)+       , unsafeFmapTree        -- ** Triangulation-builders        , TriangBuild, doTriangBuild        , AutoTriang, breakdownAutoTriang@@ -106,6 +105,7 @@ import Data.Functor.Identity import Control.Monad.Trans.State import Control.Monad.Trans.Writer+import Control.Monad.Trans.List import Control.Monad.Trans.OuterMaybe import Control.Monad.Trans.Class import qualified Data.Foldable       as Hask@@ -123,6 +123,7 @@ import GHC.Generics (Generic) import Data.Type.Coercion +import Development.Placeholders   type Depth = Int@@ -158,7 +159,7 @@ -- | Hourglass as the geometric shape (two opposing ~conical volumes, sharing --   only a single point in the middle); has nothing to do with time. data Hourglass s = Hourglass { upperBulb, lowerBulb :: !s }-            deriving (Generic, Hask.Functor, Hask.Foldable, Show)+            deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable, Show) instance (NFData s) => NFData (Hourglass s) instance (Semigroup s) => Semigroup (Hourglass s) where   Hourglass u l <> Hourglass u' l' = Hourglass (u<>u') (l<>l')@@ -184,24 +185,28 @@ oneBulb LowerBulb f (Hourglass u l) = Hourglass u (f l)  +type LeafCount = Int+type LeafIndex = Int -data ShadeTree x = PlainLeaves [x]-                 | DisjointBranches !Int (NonEmpty (ShadeTree x))-                 | OverlappingBranches !Int !(Shade x) (NonEmpty (DBranch x))-  deriving (Generic)+type ShadeTree x = x`Shaded`()++data Shaded x y = PlainLeaves [(x,y)]+                | DisjointBranches !LeafCount (NonEmpty (x`Shaded`y))+                | OverlappingBranches !LeafCount !(Shade x) (NonEmpty (DBranch x y))+  deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable) deriving instance ( WithField ℝ PseudoAffine x, Show x                   , Show (Interior x), Show (Needle' x), Show (Metric' x) )              => Show (ShadeTree x)             data DBranch' x c = DBranch { boughDirection :: !(Needle' x)                             , boughContents :: !(Hourglass c) }-  deriving (Generic, Hask.Functor, Hask.Foldable)-type DBranch x = DBranch' x (ShadeTree x)+  deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable)+type DBranch x y = DBranch' x (x`Shaded`y) deriving instance ( WithField ℝ PseudoAffine x, Show (Needle' x), Show c )              => Show (DBranch' x c)  newtype DBranches' x c = DBranches (NonEmpty (DBranch' x c))-  deriving (Generic, Hask.Functor, Hask.Foldable)+  deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable) deriving instance ( WithField ℝ PseudoAffine x, Show (Needle' x), Show c )              => Show (DBranches' x c) @@ -210,17 +215,31 @@   DBranches b1 <> DBranches b2 = DBranches $ NE.zipWith (\(DBranch d1 c1) (DBranch _ c2)                                                               -> DBranch d1 $ c1<>c2 ) b1 b2 +++trunkBranches :: x`Shaded`y -> NonEmpty (LeafIndex, x`Shaded`y)+trunkBranches (OverlappingBranches _ _ brs)+        = (`evalState`0)+            . forM (brs >>= \(DBranch _ (Hourglass t b)) -> t:|[b]) $ \st -> do+               i₀ <- get+               put $ i₀ + nLeaves st+               return (i₀, st)+trunkBranches (DisjointBranches _ brs) = (`evalState`0) . forM brs $ \st -> do+               i₀ <- get+               put $ i₀ + nLeaves st+               return (i₀, st)+trunkBranches t = pure (0,t)    directionChoices :: WithField ℝ Manifold x-               => [DBranch x]-                 -> [ ( (Needle' x, ShadeTree x)-                      ,[(Needle' x, ShadeTree x)] ) ]+               => [DBranch x y]+                 -> [ ( (Needle' x, x`Shaded`y)+                      ,[(Needle' x, x`Shaded`y)] ) ] directionChoices = map (snd *** map snd) . directionIChoices 0  directionIChoices :: (WithField ℝ PseudoAffine x, AdditiveGroup (Needle' x))-               => Int -> [DBranch x]-                 -> [ ( (Int, (Needle' x, ShadeTree x))-                      ,[(Int, (Needle' x, ShadeTree x))] ) ]+               => Int -> [DBranch x y]+                 -> [ ( (Int, (Needle' x, x`Shaded`y))+                      ,[(Int, (Needle' x, x`Shaded`y))] ) ] directionIChoices _ [] = [] directionIChoices i₀ (DBranch ѧ (Hourglass t b) : hs)          =  ( top, bot : map fst uds )@@ -230,13 +249,12 @@        bot = (i₀+1,(negateV ѧ,b))        uds = directionIChoices (i₀+2) hs -traverseDirectionChoices :: ( WithField ℝ PseudoAffine x, LSpace (Needle x)-                            , Hask.Applicative f )-               => (    (Int, (Needle' x, ShadeTree x))-                    -> [(Int, (Needle' x, ShadeTree x))]-                    -> f (ShadeTree x) )-                 -> [DBranch x]-                 -> f [DBranch x]+traverseDirectionChoices :: ( AdditiveGroup (Needle' x), Hask.Applicative f )+               => (    (Int, (Needle' x, x`Shaded`y))+                    -> [(Int, (Needle' x, x`Shaded`y))]+                    -> f (x`Shaded`z) )+                 -> [DBranch x y]+                 -> f [DBranch x z] traverseDirectionChoices f dbs            = td [] . scanLeafNums 0                $ dbs >>= \(DBranch ѧ (Hourglass τ β))@@ -252,7 +270,25 @@        scanLeafNums i₀ ((v,t):vts) = (i₀, (v,t)) : scanLeafNums (i₀ + nLeaves t) vts  -indexDBranches :: NonEmpty (DBranch x) -> NonEmpty (DBranch' x (Int, ShadeTree x))++traverseTrunkBranchChoices :: Hask.Applicative f+               => ( (Int, x`Shaded`y) -> x`Shaded`y -> f (x`Shaded`z) )+                 -> x`Shaded`y -> f (x`Shaded`z)+traverseTrunkBranchChoices f (OverlappingBranches n sh bs)+        = OverlappingBranches n sh . NE.fromList <$> go 0 id (NE.toList bs)+ where go _ _ [] = pure []+       go i₀ prbs (tbs@(DBranch v (Hourglass τ β)) : dbs)+        = (:) . DBranch v <$>+            (Hourglass <$> (f (i₀, τ) . OverlappingBranches (n-nτ) sh+                            . NE.fromList . prbs $ DBranch v (Hourglass hole β) : dbs)+                       <*> (f (i₀+nτ, β) . OverlappingBranches (n-nβ) sh+                            . NE.fromList . prbs $ DBranch v (Hourglass τ hole) : dbs))+            <*> go (i₀+nτ+nβ) (prbs . (tbs:)) dbs+        where [nτ, nβ] = nLeaves<$>[τ,β]+              hole = PlainLeaves []+++indexDBranches :: NonEmpty (DBranch x y) -> NonEmpty (DBranch' x (Int, x`Shaded`y)) indexDBranches (DBranch d (Hourglass t b) :| l) -- this could more concisely be written as a traversal               = DBranch d (Hourglass (0,t) (nt,b)) :| ixDBs (nt + nb) l  where nt = nLeaves t; nb = nLeaves b@@ -261,31 +297,18 @@                = DBranch δ (Hourglass (i₀,τ) (i₀+nτ,β)) : ixDBs (i₀ + nτ + nβ) l         where nτ = nLeaves τ; nβ = nLeaves β -instance (NFData x, NFData (Needle' x)) => NFData (ShadeTree x) where+instance (NFData x, NFData (Needle' x), NFData y) => NFData (x`Shaded`y) where   rnf (PlainLeaves xs) = rnf xs   rnf (DisjointBranches n bs) = n `seq` rnf (NE.toList bs)   rnf (OverlappingBranches n sh bs) = n `seq` sh `seq` rnf (NE.toList bs)-instance (NFData x, NFData (Needle' x)) => NFData (DBranch x)+instance (NFData x, NFData (Needle' x), NFData y) => NFData (DBranch x y)   --- | Experimental. There might be a more powerful instance possible.-instance (AffineManifold x) => Semimanifold (ShadeTree x) where-  type Needle (ShadeTree x) = Diff x-  fromInterior = id-  toInterior = pure-  translateP = Tagged (.+~^)-  PlainLeaves xs .+~^ v = PlainLeaves $ (.+^v)<$>xs -  OverlappingBranches n sh br .+~^ v-        = OverlappingBranches n (sh.+~^v)-                $ fmap (\(DBranch d c) -> DBranch d $ (.+~^v)<$>c) br-  DisjointBranches n br .+~^ v = DisjointBranches n $ (.+~^v)<$>br-  semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of-     SemimanifoldWitness BoundarylessWitness -> SemimanifoldWitness BoundarylessWitness  -- | WRT union. instance (WithField ℝ Manifold x, SimpleSpace (Needle x)) => Semigroup (ShadeTree x) where   PlainLeaves [] <> t = t   t <> PlainLeaves [] = t-  t <> s = fromLeafPoints $ onlyLeaves t ++ onlyLeaves s+  t <> s = fromLeafPoints $ onlyLeaves_ t ++ onlyLeaves_ s            -- Could probably be done more efficiently   sconcat = mconcat . NE.toList instance (WithField ℝ Manifold x, SimpleSpace (Needle x)) => Monoid (ShadeTree x) where@@ -294,7 +317,7 @@   mconcat l = case filter ne l of                [] -> mempty                [t] -> t-               l' -> fromLeafPoints $ onlyLeaves =<< l'+               l' -> fromLeafPoints $ onlyLeaves_ =<< l'    where ne (PlainLeaves []) = False; ne _ = True  @@ -304,14 +327,17 @@ --  -- <<images/examples/simple-2d-ShadeTree.png>> fromLeafPoints :: ∀ x. (WithField ℝ Manifold x, SimpleSpace (Needle x))-                         => [x] -> ShadeTree x-fromLeafPoints = fromLeafPoints' sShIdPartition+                        => [x] -> ShadeTree x+fromLeafPoints = fromLeafPoints_ . map (,()) +fromLeafPoints_ :: ∀ x y. (WithField ℝ Manifold x, SimpleSpace (Needle x))+                        => [(x,y)] -> x`Shaded`y+fromLeafPoints_ = fromLeafPoints' sShIdPartition + -- | The leaves of a shade tree are numbered. For a given index, this function --   attempts to find the leaf with that ID, within its immediate environment.-indexShadeTree :: ∀ x . WithField ℝ Manifold x-       => ShadeTree x -> Int -> Either Int ([ShadeTree x], x)+indexShadeTree :: ∀ x y . x`Shaded`y -> Int -> Either Int ([x`Shaded`y], (x,y)) indexShadeTree _ i     | i<0        = Left i indexShadeTree sh@(PlainLeaves lvs) i = case length lvs of@@ -330,21 +356,42 @@                          ) (Left i) (toList brs>>=toList)     | otherwise  = Left $ i-n +treeLeaf :: ∀ x y f . Hask.Functor f+        => Int -> (y -> f y) -> x`Shaded`y -> Either Int (f (x`Shaded`y))+treeLeaf i _ _+    | i<0        = Left i+treeLeaf i f sh@(PlainLeaves lvs) = case length lvs of+  n | i<n+    , (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  -- | “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 . (WithField ℝ Manifold x, SimpleSpace (Needle x))+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.-        -> ShadeTree x       -- ^ The tree to index into+        -> x`Shaded`y        -- ^ The tree to index into         -> x                 -- ^ Position to look up-        -> Maybe (Int, ([ShadeTree x], x))+        -> Maybe (Int, ([x`Shaded`y], (x,y)))                    -- ^ Index of the leaf near to the query point, the “path” of                    --   environment trees leading down to its position (in decreasing-                   --   order of size), and actual position of the found node.+                   --   order of size), and actual position+info of the found node. positionIndex (Just m) sh@(PlainLeaves lvs) x-        = case catMaybes [ ((i,p),) . normSq m <$> p.-~.x+        = case catMaybes [ ((i,p),) . normSq m <$> fst p.-~.x                             | (i,p) <- zip [0..] lvs] of            [] -> empty            l | ((i,p),_) <- minimumBy (comparing snd) l@@ -371,27 +418,14 @@   -fromFnGraphPoints :: ∀ x y . ( WithField ℝ Manifold x, WithField ℝ Manifold y-                             , SimpleSpace (Needle x), SimpleSpace (Needle y) )-                     => [(x,y)] -> ShadeTree (x,y)-fromFnGraphPoints = case ( dualSpaceWitness :: DualNeedleWitness x-                         , boundarylessWitness :: BoundarylessWitness x-                         , dualSpaceWitness :: DualNeedleWitness y-                         , boundarylessWitness :: BoundarylessWitness y ) of-    (DualSpaceWitness,BoundarylessWitness,DualSpaceWitness,BoundarylessWitness)-        -> fromLeafPoints' $-     \(Shade c expa) xs -> case-            [ DBranch (v, zeroV) mempty-            | v <- normSpanningSystem' (transformNorm (id&&&zeroV) expa :: Metric' x) ] of-         (b:bs) -> sShIdPartition' c xs $ b:|bs -fromLeafPoints' :: ∀ x. (WithField ℝ Manifold x, SimpleSpace (Needle x)) =>-    (Shade x -> [x] -> NonEmpty (DBranch' x [x])) -> [x] -> ShadeTree x+fromLeafPoints' :: ∀ x y. (WithField ℝ Manifold x, SimpleSpace (Needle x)) =>+    (Shade x -> [(x,y)] -> NonEmpty (DBranch' x [(x,y)])) -> [(x,y)] -> x`Shaded`y fromLeafPoints' sShIdPart = go boundarylessWitness mempty- where go :: BoundarylessWitness x -> Metric' x -> [x] -> ShadeTree x+ where go :: BoundarylessWitness x -> Metric' x -> [(x,y)] -> x`Shaded`y        go bw@BoundarylessWitness preShExpa             = \xs -> case pointsShades' (scaleNorm (1/3) preShExpa) xs of-                     [] -> mempty+                     [] -> PlainLeaves []                      [(_,rShade)] -> let trials = sShIdPart rShade xs                                      in case reduce rShade trials of                                          Just redBrchs@@ -405,8 +439,8 @@         where                branchProc redSh = fmap (fmap $ go bw redSh)                                  -              reduce :: Shade x -> NonEmpty (DBranch' x [x])-                                      -> Maybe (NonEmpty (DBranch' x [x]))+              reduce :: Shade x -> NonEmpty (DBranch' x [(x,y)])+                                      -> Maybe (NonEmpty (DBranch' x [(x,y)]))               reduce sh@(Shade c _) brCandidates                         = case findIndex deficient cards of                             Just i | (DBranch _ reBr, o:ok)@@ -422,16 +456,17 @@   sShIdPartition' :: (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-        => Interior x -> [x] -> NonEmpty (DBranch' x [x])->NonEmpty (DBranch' x [x])+        => Interior x -> [(x,y)] -> NonEmpty (DBranch' x [(x,y)])+                                 -> NonEmpty (DBranch' x [(x,y)]) sShIdPartition' c xs st-           = foldr (\p -> let (i,h) = ssi p+           = foldr (\(p,y) -> let (i,h) = ssi p                           in asList $ update_nth (\(DBranch d c)-                                                    -> DBranch d (oneBulb h (p:) c))+                                                    -> DBranch d (oneBulb h ((p,y):) c))                                       i )                    st xs  where ssi = subshadeId' (fromInterior c) (boughDirection<$>st) sShIdPartition :: (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                    => Shade x -> [x] -> NonEmpty (DBranch' x [x])+                    => Shade x -> [(x,y)] -> NonEmpty (DBranch' x [(x,y)]) sShIdPartition (Shade c expa) xs  | b:bs <- [DBranch v mempty | v <- normSpanningSystem' expa]     = sShIdPartition' c xs $ b:|bs@@ -470,59 +505,48 @@ sortByKey = map snd . sortBy (comparing fst)  -trunks :: ∀ x. (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                  => ShadeTree x -> [Shade x]+trunks :: ∀ x y . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))+                  => x`Shaded`y -> [Shade x] trunks t = case (pseudoAffineWitness :: PseudoAffineWitness x, t) of   (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness), PlainLeaves lvs)-                                         -> pointsCovers . catMaybes $ toInterior<$>lvs-  (_, DisjointBranches _ brs)            -> Hask.foldMap trunks brs-  (_, OverlappingBranches _ sh _)        -> [sh]+                                    -> pointsCovers . catMaybes $ toInterior.fst<$>lvs+  (_, DisjointBranches _ brs)       -> Hask.foldMap trunks brs+  (_, OverlappingBranches _ sh _)   -> [sh]  -nLeaves :: ShadeTree x -> Int+nLeaves :: x`Shaded`y -> Int nLeaves (PlainLeaves lvs) = length lvs nLeaves (DisjointBranches n _) = n nLeaves (OverlappingBranches n _ _) = n +treeDepth :: x`Shaded`y -> Int+treeDepth (PlainLeaves lvs) = 0+treeDepth (DisjointBranches _ brs) = 1 + maximum (treeDepth<$>brs)+treeDepth (OverlappingBranches _ _ brs)+     = 1 + maximum (maximum . fmap treeDepth<$>brs) -instance ImpliesMetric ShadeTree where-  type MetricRequirement ShadeTree x = (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-  inferMetric = stInfMet-   where stInfMet :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                                => ShadeTree x -> Metric x-         stInfMet (OverlappingBranches _ (Shade _ e) _) = dualNorm' e-         stInfMet (PlainLeaves lvs)-               = case pointsShades $ Hask.toList . toInterior =<< lvs :: [Shade x] of-             (Shade _ sh:_) -> dualNorm' sh-             _ -> mempty-         stInfMet (DisjointBranches _ (br:|_)) = inferMetric br-  inferMetric' = stInfMet-   where stInfMet :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                                => ShadeTree x -> Metric' x-         stInfMet (OverlappingBranches _ (Shade _ e) _) = e-         stInfMet (PlainLeaves lvs)-               = case pointsShades $ Hask.toList . toInterior =<< lvs :: [Shade x] of-             (Shade _ sh:_) -> sh-             _ -> mempty-         stInfMet (DisjointBranches _ (br:|_)) = inferMetric' br   -overlappingBranches :: Shade x -> NonEmpty (DBranch x) -> ShadeTree x++overlappingBranches :: Shade x -> NonEmpty (DBranch x y) -> x`Shaded`y overlappingBranches shx brs = OverlappingBranches n shx brs  where n = sum $ fmap (sum . fmap nLeaves) brs -unsafeFmapLeaves :: (x -> x) -> ShadeTree x -> ShadeTree x+unsafeFmapLeaves_ :: (x -> x) -> x`Shaded`y -> x`Shaded`y+unsafeFmapLeaves_ = unsafeFmapLeaves . first++unsafeFmapLeaves :: ((x,y) -> (x,y')) -> x`Shaded`y -> x`Shaded`y' unsafeFmapLeaves f (PlainLeaves lvs) = PlainLeaves $ fmap f lvs unsafeFmapLeaves f (DisjointBranches n brs)-                  = DisjointBranches n $ unsafeFmapLeaves f <$> brs+                 = DisjointBranches n $ unsafeFmapLeaves f <$> brs unsafeFmapLeaves f (OverlappingBranches n sh brs)                   = OverlappingBranches n sh $ fmap (unsafeFmapLeaves f) <$> brs -unsafeFmapTree :: (NonEmpty x -> NonEmpty y)-               -> (Needle' x -> Needle' y)-               -> (Shade x -> Shade y)-               -> ShadeTree x -> ShadeTree y+unsafeFmapTree :: (NonEmpty (x,y) -> NonEmpty (ξ,υ))+               -> (Needle' x -> Needle' ξ)+               -> (Shade x -> Shade ξ)+               -> x`Shaded`y -> ξ`Shaded`υ unsafeFmapTree _ _ _ (PlainLeaves []) = PlainLeaves [] unsafeFmapTree f _ _ (PlainLeaves lvs) = PlainLeaves . toList . f $ NE.fromList lvs unsafeFmapTree f fn fs (DisjointBranches n brs)@@ -534,22 +558,13 @@                       ) brs       in overlappingBranches (fs sh) brs' -coerceShadeTree :: ∀ x y . (LocallyCoercible x y, Manifold x, Manifold y)-                       => ShadeTree x -> ShadeTree y-coerceShadeTree = case ( dualSpaceWitness :: DualNeedleWitness x-                       , dualSpaceWitness :: DualNeedleWitness y ) of-   (DualSpaceWitness,DualSpaceWitness)-      -> unsafeFmapTree (fmap locallyTrivialDiffeomorphism)-                                 (coerceNeedle' ([]::[(x,y)]) $)-                                 coerceShade   --type Twig x = (Int, ShadeTree x)-type TwigEnviron x = [Twig x]+type Twig x y = (Int, x`Shaded`y)+type TwigEnviron x y = [Twig x y] -allTwigs :: ∀ x . WithField ℝ PseudoAffine x => ShadeTree x -> [Twig x]+allTwigs :: ∀ x y . WithField ℝ PseudoAffine x => x`Shaded`y -> [Twig x y] allTwigs tree = go 0 tree []  where go n₀ (DisjointBranches _ dp)          = snd (foldl' (\(n₀',prev) br -> (n₀'+nLeaves br, prev . go n₀' br)) (n₀,id) dp)@@ -567,15 +582,15 @@ -- | Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/Trees-and-Webs.ipynb#pseudorandomCloudTree --  --   <<images/examples/TreesAndWebs/2D-scatter_twig-environs.png>>-twigsWithEnvirons :: ∀ x. (WithField ℝ Manifold x, SimpleSpace (Needle x))-    => ShadeTree x -> [(Twig x, TwigEnviron x)]+twigsWithEnvirons :: ∀ x y. (WithField ℝ Manifold x, SimpleSpace (Needle x))+    => x`Shaded`y -> [(Twig x y, TwigEnviron x y)] twigsWithEnvirons = execWriter . traverseTwigsWithEnvirons (writer . (snd.fst&&&pure)) -traverseTwigsWithEnvirons :: ∀ x f .+traverseTwigsWithEnvirons :: ∀ x y f .             (WithField ℝ PseudoAffine x, SimpleSpace (Needle x), Hask.Applicative f)-    => ( (Twig x, TwigEnviron x) -> f (ShadeTree x) ) -> ShadeTree x -> f (ShadeTree x)+    => ( (Twig x y, TwigEnviron x y) -> f (x`Shaded`y) ) -> x`Shaded`y -> f (x`Shaded`y) traverseTwigsWithEnvirons f = fst . go pseudoAffineWitness [] . (0,)- where go :: PseudoAffineWitness x -> TwigEnviron x -> Twig x -> (f (ShadeTree x), Bool)+ where go :: PseudoAffineWitness x -> TwigEnviron x y -> Twig x y -> (f (x`Shaded`y), Bool)        go sw _ (i₀, DisjointBranches nlvs djbs) = ( fmap (DisjointBranches nlvs)                                                    . Hask.traverse (fst . go sw [])                                                    $ NE.zip ioffs djbs@@ -612,7 +627,7 @@        go (PseudoAffineWitness (SemimanifoldWitness _)) envi plvs@(i₀, (PlainLeaves _))                          = (f $ purgeRemotes (plvs, envi), True)        -       twigProximæ :: PseudoAffineWitness x -> Interior x -> ShadeTree x -> TwigEnviron x+       twigProximæ :: PseudoAffineWitness x -> Interior x -> x`Shaded`y -> TwigEnviron x y        twigProximæ sw x₀ (DisjointBranches _ djbs)                = Hask.foldMap (\(i₀,st) -> first (+i₀) <$> twigProximæ sw x₀ st)                     $ NE.zip ioffs djbs@@ -627,7 +642,7 @@                where overlap = bdir<.>^δxb        twigProximæ _ _ plainLeaves = [(0, plainLeaves)]        -       twigsaveTrim :: (DBranch x -> TwigEnviron x) -> ShadeTree x -> TwigEnviron x+       twigsaveTrim :: (DBranch x y -> TwigEnviron x y) -> x`Shaded`y -> TwigEnviron x y        twigsaveTrim f ct@(OverlappingBranches _ _ dbs)                  = case Hask.mapM (\(i₀,dbr) -> noLeaf $ first(+i₀)<$>f dbr)                                  $ NE.zip ioffs dbs of@@ -637,7 +652,7 @@               noLeaf bqs = pure bqs               ioffs = NE.scanl (\i -> (+i) . sum . fmap nLeaves . toList) 0 dbs        -       purgeRemotes :: (Twig x, TwigEnviron x) -> (Twig x, TwigEnviron x)+       purgeRemotes :: (Twig x y, TwigEnviron x y) -> (Twig x y, TwigEnviron x y)        purgeRemotes = id -- See 7d1f3a4 for the implementation; this didn't work reliable.       completeTopShading :: ∀ x y . ( WithField ℝ PseudoAffine x, WithField ℝ PseudoAffine y@@ -647,14 +662,14 @@                                              , dualSpaceWitness :: DualNeedleWitness y ) of        (DualSpaceWitness, DualSpaceWitness)           -> pointsShade's . catMaybes-               $ toInterior . (_topological &&& _untopological) <$> plvs+               $ toInterior <$> plvs completeTopShading (DisjointBranches _ bqs)                      = take 1 . completeTopShading =<< NE.toList bqs completeTopShading t = case ( dualSpaceWitness :: DualNeedleWitness x                             , dualSpaceWitness :: DualNeedleWitness y ) of        (DualSpaceWitness, DualSpaceWitness)           -> pointsCover's . catMaybes-                . map (toInterior <<< _topological &&& _untopological) $ onlyLeaves t+                . map toInterior $ onlyLeaves t   transferAsNormsDo :: ∀ v . LSpace v => Norm v -> Variance v -> v-+>v@@ -682,8 +697,8 @@               fts (xc, (Shade' yc expay, jtg)) = unsafeFmapLeaves applδj t                where Just δyc = yc.-~.yc₀                      tfm = transferAsNormsDo expay₀ (dualNorm expay)-                     applδj (WithAny y x)-                           = WithAny (yc₀ .+~^ ((tfm $ δy) ^+^ (jtg $ δx) ^+^ δyc)) x+                     applδj (x,y)+                           = (x, yc₀ .+~^ ((tfm $ δy) ^+^ (jtg $ δx) ^+^ δyc))                       where Just δx = x.-~.xc                             Just δy = y.-~.(yc₀.+~^(j₀ $ δx))        @@ -705,24 +720,26 @@                   -seekPotentialNeighbours :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                => ShadeTree x -> x`Shaded`[Int]+seekPotentialNeighbours :: ∀ x y . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))+                => x`Shaded`y -> x`Shaded`(y,[Int]) seekPotentialNeighbours tree = zipTreeWithList tree-                     $ snd<$>leavesWithPotentialNeighbours tree+                     $ case snd<$>leavesWithPotentialNeighbours tree of+                         (n:ns) -> n:|ns -leavesWithPotentialNeighbours :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                => ShadeTree x -> [(x, [Int])]+leavesWithPotentialNeighbours :: ∀ x y+            . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))+                => x`Shaded`y -> [((x,y), [Int])] leavesWithPotentialNeighbours = map (second snd) . go pseudoAffineWitness 0 0 []- where go :: PseudoAffineWitness x -> Depth -> Int -> [Wall x] -> ShadeTree x-                -> [(x, ([Wall x], [Int]))]+ where go :: PseudoAffineWitness x -> Depth -> Int -> [Wall x] -> x`Shaded`y+                -> [((x,y), ([Wall x], [Int]))]        go (PseudoAffineWitness (SemimanifoldWitness _)) depth n₀ walls (PlainLeaves lvs)-               = [ (x, ( [ wall & wallDistance .~ d+               = [ ((x,y), ( [ wall & wallDistance .~ d                          | wall <- walls                          , Just vw <- [toInterior x>>=(.-~.wall^.wallAnchor)]                          , let d = (wall^.wallNormal)<.>^vw                          , d < wall^.wallDistance ]                        , [] ))-                 | x <- lvs ]+                 | (x,y) <- lvs ]        go pw depth n₀ walls (DisjointBranches _ dp)          = snd (foldl' (\(n₀',prev) br -> ( n₀'+nLeaves br                                           , prev . (go pw depth n₀' walls br++)))@@ -731,9 +748,9 @@                depth n₀ walls (OverlappingBranches _ (Shade brCtr _) dp)          = reassemble $ snd              (foldl' assignWalls (n₀,id) . directionIChoices 0 $ NE.toList dp) []-        where assignWalls :: (Int, DList (x, ([Wall x],[Int])))-                     -> ((Int,(Needle' x, ShadeTree x)), [(Int,(Needle' x, ShadeTree x))])-                     -> (Int, DList (x, ([Wall x], [Int])))+        where assignWalls :: (Int, DList ((x,y), ([Wall x],[Int])))+                     -> ((Int,(Needle' x, x`Shaded`y)), [(Int,(Needle' x, x`Shaded`y))])+                     -> (Int, DList ((x,y), ([Wall x], [Int])))               assignWalls (n₀',prev) ((iDir,(thisDir,br)),otherDirs)                     = ( n₀'+nLeaves br                       , prev . (go pw (depth+1) n₀'@@ -747,7 +764,7 @@                      updWall wall = wall & wallDistance %~ min bcDist                       where Just vbw = brCtr.-~.wall^.wallAnchor                             bcDist = (wall^.wallNormal)<.>^vbw-              reassemble :: [(x, ([Wall x],[Int]))] -> [(x, ([Wall x],[Int]))]+              reassemble :: [((x,y), ([Wall x],[Int]))] -> [((x,y), ([Wall x],[Int]))]               reassemble pts = [ (x, (higherWalls, newGroups++deeperGroups))                                | (x, (allWalls, deeperGroups)) <- pts                                , let (levelWalls,higherWalls)@@ -902,6 +919,8 @@ -- 'NonEmptyTree' x &#x2245; (x, 'Trees' x) -- @ type NonEmptyTree = GenericTree NonEmpty []++type LeafyTree x y = GenericTree [] (ListT (Either y)) x      newtype GenericTree c b x = GenericTree { treeBranches :: c (x,GenericTree b b x) }  deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable)@@ -912,22 +931,49 @@ instance (Hask.MonadPlus c) => Monoid (GenericTree c b x) where   mempty = GenericTree Hask.mzero   mappend = (<>)-deriving instance Show (c (x, GenericTree b b x)) => Show (GenericTree c b x)+instance Show (c (x, GenericTree b b x)) => Show (GenericTree c b x) where+  showsPrec p (GenericTree t) = showParen (p>9) $ ('朳':) . showsPrec 10 t+deriving instance Eq (c (x, GenericTree b b x)) => Eq (GenericTree c b x) +-- | @U+6733 CJK UNIFIED IDEOGRAPH tree@.+--  The main purpose of this is to give 'GenericTree' a more concise 'Show' instance.+朳 :: c (x, GenericTree b b x) -> GenericTree c b x+朳 = GenericTree+ -- | Imitate the specialised 'ShadeTree' structure with a simpler, generic tree. onlyNodes :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))                 => ShadeTree x -> Trees x onlyNodes (PlainLeaves []) = GenericTree [] onlyNodes (PlainLeaves ps) = let (ctr,_) = pseudoECM ([]::[x]) $ NE.fromList ps-                             in GenericTree [ (ctr, GenericTree $ (,mempty) <$> ps) ]+                             in GenericTree [ (ctr, GenericTree $ (,mempty).fst <$> ps) ] onlyNodes (DisjointBranches _ brs) = Hask.foldMap onlyNodes brs onlyNodes (OverlappingBranches _ (Shade ctr _) brs)               = GenericTree [ ( fromInterior ctr                               , Hask.foldMap (Hask.foldMap onlyNodes) brs ) ] +entireTree :: ∀ x y . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))+              => x`Shaded`y -> LeafyTree x y+entireTree (PlainLeaves lvs)+    = let (ctr,_) = pseudoECM ([]::[x]) $ NE.fromList lvs+      in  GenericTree [ (ctr, GenericTree . ListT $ Right+                                [ (x, GenericTree . lift $ Left y)+                                | (x,y)<-lvs ] )+                      ]+entireTree (DisjointBranches _ brs)+    = GenericTree [ (x, GenericTree subt)+                  | GenericTree sub <- NE.toList $ fmap entireTree brs+                  , (x, GenericTree subt) <- sub ]+entireTree (OverlappingBranches _ (Shade ctr _) brs)+    = GenericTree [ ( fromInterior ctr+                    , GenericTree . ListT . Right+                       $ Hask.foldMap (Hask.foldMap $ treeBranches . entireTree) brs ) ] + -- | Left (and, typically, also right) inverse of 'fromLeafNodes'.-onlyLeaves :: WithField ℝ PseudoAffine x => ShadeTree x -> [x]+onlyLeaves_ :: WithField ℝ PseudoAffine x => ShadeTree x -> [x]+onlyLeaves_ = map fst . onlyLeaves++onlyLeaves :: WithField ℝ PseudoAffine x => x`Shaded`y -> [(x,y)] onlyLeaves tree = dismantle tree []  where dismantle (PlainLeaves xs) = (xs++)        dismantle (OverlappingBranches _ _ brs)@@ -970,33 +1016,23 @@   -constShaded :: y -> ShadeTree x -> x`Shaded`y-constShaded y = unsafeFmapTree (WithAny y<$>) id (shadeWithAny y)--stripShadedUntopological :: x`Shaded`y -> ShadeTree x-stripShadedUntopological = unsafeFmapTree (fmap _topological) id shadeWithoutAnything--fmapShaded :: (y -> υ) -> (x`Shaded`y) -> (x`Shaded`υ)-fmapShaded f = unsafeFmapTree (fmap $ \(WithAny y x) -> WithAny (f y) x)-                              id-                              (\(Shade yx shx) -> Shade (fmap f yx) shx)+constShaded :: y -> x`Shaded`y₀ -> x`Shaded`y+constShaded y = unsafeFmapTree (fmap . second $ const y) id id -joinShaded :: (x`WithAny`y)`Shaded`z -> x`Shaded`(y,z)-joinShaded = unsafeFmapTree (fmap $ \(WithAny z (WithAny y x)) -> WithAny (y,z) x)-                            id-                            (\(Shade (WithAny z (WithAny y x)) shx)-                                  -> Shade (WithAny (y,z) x) shx )+fmapShaded :: (Semimanifold x, SimpleSpace (Needle x))+                   => (y -> υ) -> (x`Shaded`y) -> (x`Shaded`υ)+fmapShaded f = unsafeFmapTree (fmap $ second f) id id -zipTreeWithList :: ShadeTree x -> [y] -> (x`Shaded`y)-zipTreeWithList tree = go tree . cycle- where go (PlainLeaves lvs) ys = PlainLeaves $ zipWith WithAny ys lvs+zipTreeWithList :: x`Shaded`w -> NonEmpty y -> (x`Shaded`(w,y))+zipTreeWithList tree = go tree . NE.toList . NE.cycle+ where go (PlainLeaves lvs) ys = PlainLeaves $ zipWith (\(x,w) y -> (x,(w,y))) lvs ys        go (DisjointBranches n brs) ys              = DisjointBranches n . NE.fromList                   $ snd (foldl (\(ys',prev) br ->                                      (drop (nLeaves br) ys', prev . (go br ys':)) )                            (ys,id) $ NE.toList brs) []-       go (OverlappingBranches n (Shade xoc shx) brs) ys-             = OverlappingBranches n (Shade (WithAny (head ys) xoc) shx) . NE.fromList+       go (OverlappingBranches n shx brs) ys+             = OverlappingBranches n shx . NE.fromList                   $ snd (foldl (\(ys',prev) (DBranch dir (Hourglass top bot))                         -> case drop (nLeaves top) ys' of                               ys'' -> ( drop (nLeaves bot) ys''@@ -1005,17 +1041,14 @@                                       ) )                            (ys,id) $ NE.toList brs) [] --- | This is to 'ShadeTree' as 'Data.Map.Map' is to 'Data.Set.Set'.-type x`Shaded`y = ShadeTree (x`WithAny`y)- stiWithDensity :: ∀ x y . ( WithField ℝ PseudoAffine x, LinearSpace y, Scalar y ~ ℝ                           , SimpleSpace (Needle x) )          => x`Shaded`y -> x -> Cℝay y stiWithDensity (PlainLeaves lvs)   | [Shade baryc expa :: Shade x] <- pointsShades . catMaybes -                                       $ toInterior . _topological <$> lvs+                                       $ toInterior . fst <$> lvs        = let nlvs = fromIntegral $ length lvs :: ℝ-             indiShapes = [(Shade pi expa, y) | WithAny y p <- lvs+             indiShapes = [(Shade pi expa, y) | (p,y) <- lvs                                               , Just pi <- [toInterior p]]          in \x -> let lcCoeffs = [ occlusion psh x | (psh, _) <- indiShapes ]                       dens = sum lcCoeffs@@ -1025,7 +1058,7 @@            = \x -> foldr1 qGather $ (`stiWithDensity`x)<$>lvs  where qGather (Cℝay 0 _) o = o        qGather o _ = o-stiWithDensity (OverlappingBranches n (Shade (WithAny _ bc) extend) brs)+stiWithDensity (OverlappingBranches n (Shade bc extend) brs)            = ovbSWD (dualSpaceWitness, pseudoAffineWitness)  where ovbSWD :: (DualNeedleWitness x, PseudoAffineWitness x) -> x -> Cℝay y        ovbSWD (DualSpaceWitness, PseudoAffineWitness (SemimanifoldWitness _)) x@@ -1051,35 +1084,19 @@                  -> ( xloc, ( (yloc, recip $ shd|$|(0,1))                             , dependence (dualNorm shd) ) ) -smoothInterpolate :: ∀ x y . ( WithField ℝ Manifold x, LinearSpace y, Scalar y ~ ℝ-                             , SimpleSpace (Needle x) )-             => NonEmpty (x,y) -> x -> y-smoothInterpolate = si boundarylessWitness- where si :: BoundarylessWitness x -> NonEmpty (x,y) -> x -> y-       si BoundarylessWitness l = \x ->-             case ltr x of-               Cℝay 0 _ -> defy-               Cℝay _ y -> y-        where defy = linearCombo [(y, 1/n) | WithAny y _ <- l']-              n = fromIntegral $ length l'-              l' = (uncurry WithAny . swap) <$> NE.toList l-              ltr = stiWithDensity $ fromLeafPoints l' - spanShading :: ∀ x y . ( WithField ℝ Manifold x, WithField ℝ Manifold y                        , SimpleSpace (Needle x), SimpleSpace (Needle y) )           => (Shade x -> Shade y) -> ShadeTree x -> x`Shaded`y-spanShading f = unsafeFmapTree addYs id addYSh- where addYs :: NonEmpty x -> NonEmpty (x`WithAny`y)-       addYs l = foldr (NE.<|) (fmap (WithAny $ fromInterior ymid) l     )-                               (fmap (`WithAny` fromInterior xmid) yexamp)+spanShading f = unsafeFmapTree (addYs . fmap fst) id id+ where addYs :: NonEmpty x -> NonEmpty (x,y)+       addYs l = foldr (NE.<|) (fmap (,fromInterior ymid) l     )+                               (fmap (fromInterior xmid,) yexamp)           where [xsh@(Shade xmid _)] = pointsCovers . catMaybes . toList                                            $ toInterior<$>l                 Shade ymid yexpa = f xsh                 yexamp = [ ymid .+~^ σ*^δy                          | δy <- varianceSpanningSystem yexpa, σ <- [-1,1] ]-       addYSh :: Shade x -> Shade (x`WithAny`y)-       addYSh xsh = shadeWithAny (fromInterior . _shadeCtr $ f xsh) xsh                         

@@ -55,7 +55,7 @@         -- ** Lines         , Line(..), lineAsPlaneIntersection         -- ** Hyperplanes-        , Cutplane(..)+        , Cutplane(..), normalPlane         , fathomCutDistance, sideOfCut, cutPosBetween         -- * Linear mappings         , LinearMap, LocalLinear@@ -319,7 +319,12 @@ --   behave locally as a plane, globally as an (/n/−1)-dimensional submanifold. data Cutplane x = Cutplane { sawHandle :: x                            , cutNormal :: Stiefel1 (Needle x) }+deriving instance (Show x, Show (Needle' x)) => Show (Cutplane x) +normalPlane :: x         -- ^ Some point lying in the desired plane.+            -> Needle' x -- ^ Co-vector perpendicular to the plane. Must be nonzero.+            -> Cutplane x+normalPlane x n = Cutplane x $ Stiefel1 n   sideOfCut :: (WithField ℝ PseudoAffine x, LinearSpace (Needle x))

@@ -17,6 +17,7 @@ -- scalings, and we prefer that definition since it doesn't require a notion of -- unit length (which is only defined in inner-product spaces). +{-# LANGUAGE StandaloneDeriving, FlexibleContexts, UndecidableInstances #-}   @@ -45,3 +46,4 @@   newtype Stiefel1 v = Stiefel1 { getStiefel1N :: DualVector v }+deriving instance (Show (DualVector v)) => Show (Stiefel1 v)

@@ -38,7 +38,7 @@               -- ** Construction             , fromWebNodes, fromShadeTree_auto, fromShadeTree, fromShaded               -- ** Lookup-            , nearestNeighbour, indexWeb, webEdges, toGraph+            , nearestNeighbour, indexWeb, toGraph, webBoundary               -- ** Decomposition             , sliceWeb_lin -- , sampleWebAlongLine_lin             , sampleWeb_2Dcartesian_lin, sampleEntireWeb_2Dcartesian_lin@@ -46,25 +46,25 @@             , localFocusWeb               -- * Uncertain functions             , differentiateUncertainWebFunction, differentiate²UncertainWebFunction+            , localModels_CGrid               -- * Differential equations               -- ** Fixed resolution-            , iterateFilterDEqn_static+            , iterateFilterDEqn_static, iterateFilterDEqn_static_selective               -- ** Automatic resolution             , filterDEqnSolutions_adaptive, iterateFilterDEqn_adaptive               -- ** Configuration             , InconsistencyStrategy(..)             , InformationMergeStrategy(..)-            , naïve, inconsistencyAware, indicateInconsistencies+            , naïve, inconsistencyAware, indicateInconsistencies, postponeInconsistencies             , PropagationInconsistency(..)               -- * Misc             , ConvexSet(..), ellipsoid, ellipsoidSet, coerceWebDomain-            , rescanPDEOnWeb, rescanPDELocally, webOnions+            , rescanPDELocally, webOnions, knitShortcuts             ) where   import Data.List hiding (filter, all, foldr1) import Data.Maybe-import qualified Data.Set as Set import qualified Data.Map as Map import qualified Data.Vector as Arr import qualified Data.Vector.Mutable as MArr@@ -77,7 +77,7 @@ import Control.DeepSeq  import Data.VectorSpace-import Math.LinearMap.Category+import Math.LinearMap.Category hiding (trace)  import Data.Tagged import Data.Function (on)@@ -91,8 +91,10 @@ import Data.SetLike.Intersection import Data.Manifold.Riemannian import Data.Manifold.Atlas+import Data.Manifold.Function.LocalModel import Data.Manifold.Function.Quadratic import Data.Function.Affine+import Data.Manifold.Web.Internal import Data.Embedding      import qualified Prelude as Hask hiding(foldl, sum, sequence)@@ -103,6 +105,7 @@ import Control.Monad.Trans.State import Control.Monad.Trans.List import Control.Monad.Trans.Except+import Control.Monad.Trans.Writer hiding (censor) import Data.Functor.Identity (Identity(..)) import qualified Data.Foldable       as Hask import Data.Foldable (all, toList)@@ -119,100 +122,28 @@  import Control.Comonad (Comonad(..)) import Control.Comonad.Cofree-import Control.Lens ((&), (%~), (^.), (.~), (+~))+import Control.Lens ((&), (%~), (^.), (.~), (+~), ix) import Control.Lens.TH  import GHC.Generics (Generic) --type WebNodeId = Int--data Neighbourhood x = Neighbourhood {-     _neighbours :: UArr.Vector WebNodeId-   , _localScalarProduct :: Metric x-   }-  deriving (Generic)-makeLenses ''Neighbourhood--deriving instance ( WithField ℝ PseudoAffine x-                  , SimpleSpace (Needle x), Show (Needle' x) )-             => Show (Neighbourhood x)--data WebLocally x y = LocalWebInfo {-      _thisNodeCoord :: x-    , _thisNodeData :: y-    , _thisNodeId :: WebNodeId-    , _nodeNeighbours :: [(WebNodeId, (Needle x, WebLocally x y))]-    , _nodeLocalScalarProduct :: Metric x-    , _nodeIsOnBoundary :: Bool-    } deriving (Generic)-makeLenses ''WebLocally--data NeighbourhoodVector x = NeighbourhoodVector-          { _nvectId :: Int-          , _theNVect :: Needle x-          , _nvectNormal :: Needle' x-          , _nvectLength :: Scalar (Needle x)-          , _otherNeighboursOverlap :: Scalar (Needle x)-          }-makeLenses ''NeighbourhoodVector--data PropagationInconsistency x υ = PropagationInconsistency {-      _inconsistentPropagatedData :: [(x,υ)]-    , _inconsistentAPrioriData :: υ }-  | PropagationInconsistencies [PropagationInconsistency x υ]- deriving (Show)-makeLenses ''PropagationInconsistency--instance Monoid (PropagationInconsistency x υ) where-  mempty = PropagationInconsistencies []-  mappend p q = mconcat [p,q]-  mconcat = PropagationInconsistencies--instance (NFData x, NFData (Metric x)) => NFData (Neighbourhood x)---- | A 'PointsWeb' is almost, but not quite a mesh. It is a stongly connected†---   directed graph, backed by a tree for fast nearest-neighbour lookup of points.--- ---   †In general, there can be disconnected components, but every connected---   component is strongly connected.-data PointsWeb :: * -> * -> * where-   PointsWeb :: {-       webNodeRsc :: ShadeTree x-     , webNodeAssocData :: Arr.Vector (y, Neighbourhood x)-     } -> PointsWeb x y-  deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable)--instance (NFData x, NFData (Metric x), NFData (Needle' x), NFData y) => NFData (PointsWeb x y)--instance Foldable (PointsWeb x) (->) (->) where-  ffoldl = uncurry . Hask.foldl' . curry-  foldMap = Hask.foldMap-instance Traversable (PointsWeb x) (PointsWeb x) (->) (->) where-  traverse f (PointsWeb rsc asd)-           = fmap (PointsWeb rsc . (`Arr.zip`ngss) . Arr.fromList)-              . traverse f $ Arr.toList ys-   where (ys,ngss) = Arr.unzip asd----type MetricChoice x = Shade x -> Metric x+import Development.Placeholders   fromWebNodes :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x))                     => (MetricChoice x) -> [(x,y)] -> PointsWeb x y fromWebNodes = case boundarylessWitness :: BoundarylessWitness x of    BoundarylessWitness ->-       \mf -> fromShaded mf . fromLeafPoints . map (uncurry WithAny . swap)+       \mf -> fromShaded mf . fromLeafPoints_  fromTopWebNodes :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x))                     => (MetricChoice x) -> [((x,[Int+Needle x]),y)] -> PointsWeb x y fromTopWebNodes = case boundarylessWitness :: BoundarylessWitness x of    BoundarylessWitness ->-       \mf -> fromTopShaded mf . fromLeafPoints-                   . map (uncurry WithAny . swap . regroup')+       \mf -> fromTopShaded mf . fromLeafPoints_ . map regroup' -fromShadeTree_auto :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => ShadeTree x -> PointsWeb x ()+fromShadeTree_auto :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x))+                              => ShadeTree x -> PointsWeb x () fromShadeTree_auto = fromShaded (dualNorm' . _shadeExpanse) . constShaded ()  fromShadeTree :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x))@@ -226,192 +157,224 @@                               --   Riemannian metric).      -> (x`Shaded`y)          -- ^ Source tree.      -> PointsWeb x y-fromShaded metricf = smoothenWebTopology metricf-                   . fromTopShaded metricf . fmapShaded (first (map Left) . swap)-                       . joinShaded . seekPotentialNeighbours+fromShaded metricf = knitShortcuts metricf . autoLinkWeb . unlinkedFromShaded metricf  toShaded :: WithField ℝ PseudoAffine x => PointsWeb x y -> (x`Shaded`y)-toShaded (PointsWeb shd asd) = zipTreeWithList shd $ Arr.toList (fst<$>asd)+toShaded (PointsWeb shd) = fmap _dataAtNode shd +unlinkedFromShaded :: ∀ x y . SimpleSpace (Needle x)+                 => MetricChoice x -> (x`Shaded`y) -> PointsWeb x y+unlinkedFromShaded metricf = PointsWeb<<<fmap `id` \y+                -> Neighbourhood y mempty nm (Just dv)+ where nm = metricf $notImplemented+       dv = head $ normSpanningSystem nm++++autoLinkWeb :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x))+                => PointsWeb x y -> PointsWeb x y+autoLinkWeb = runIdentity . traverseNodesInEnvi ( pure . fetchNgbs []+                                                  . (id &&& findEnviPts (0,1)) )+ where fetchNgbs :: [(WebNodeIdOffset, Needle x)]+                 -> (NodeInWeb x y, [[(WebNodeIdOffset, (x, Neighbourhood x y))]])+                 -> Neighbourhood x y+       fetchNgbs alreadyFound+                 ( NodeInWeb (x, Neighbourhood y aprNgbs locMetr (Just wall))+                             layersAroundThis+                 , enviLayers )+         | (δi, (v, nh)) : _ <- newNgbCandidates+             = fetchNgbs+                ((δi, v) : alreadyFound)+                ( NodeInWeb (x, Neighbourhood y (UArr.cons δi aprNgbs) locMetr+                                  $ if dimension > 1+                                     then pumpHalfspace locMetr v+                                                 (wall, snd<$>alreadyFound)+                                     else case alreadyFound of+                                            [] -> Just $ locMetr<$|v+                                            [_] -> Nothing+                                                 )+                            layersAroundThis+                , enviLayers )+        where newNgbCandidates+                  = [ (δi, (v, nh))+                    | envi <- enviLayers+                    , (δi, ((v,_), nh)) <- sortBy (comparing $ snd . fst . snd)+                                  [ (δi, ((v, if dimension > 1+                                               then gatherDirectionsBadness+                                                 $ linkingUndesirability distSq wallDist+                                               else distSq+                                                 ), nh))+                                  | (δi,(xp,nh)) <- envi+                                  , let Just v = xp.-~.x+                                        distSq = normSq locMetr v+                                        wallDist = walln<.>^v+                                  , wallDist >= 0+                                  , distSq > wallDist^2+                                     || dimension==1 -- in 1D, we must allow linking+                                                     -- to the direct opposite+                                                     -- (there IS no other direction)+                                  , not . any (==δi) $ UArr.toList aprNgbs+                                                        ++ map fst alreadyFound+                                  ] ]+              locMetr' = dualNorm locMetr+              walln = wall ^/ (- (locMetr'|$|wall))+       fetchNgbs _ (NodeInWeb (_, d) _, _) = d+       findEnviPts (iw,wedgeSize) (NodeInWeb tr ((envi,iSpl):envis))+                  = (zip [-iw-iSpl ..] preds ++ zip [wedgeSize-iw ..] succs)+                     : findEnviPts (iw+iSpl, wedgeSize + iSpl + length succs)+                                   (NodeInWeb tr envis)+               where (preds, succs) = splitAt iSpl $ onlyLeaves envi+       findEnviPts _ _ = []+       dimension = subbasisDimension (entireBasis :: SubBasis (Needle x))+ fromTopShaded :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x))      => (MetricChoice x)      -> (x`Shaded`([Int+Needle x], y))                       -- ^ Source tree, with topology information                       --   (IDs of neighbour-candidates, or needles pointing to them)      -> PointsWeb x y-fromTopShaded metricf shd = PointsWeb shd' assocData - where shd' = stripShadedUntopological shd-       assocData = Hask.foldMap locMesh $ allTwigs shd-       -       locMesh :: (Int, ShadeTree (x`WithAny`([Int+Needle x], y)))-                   -> Arr.Vector (y, Neighbourhood x)-       locMesh (i₀, locT) = Arr.map findNeighbours $ Arr.fromList locLeaves-        where locLeaves :: [ (Int, x`WithAny`([Int+Needle x], y)) ]-              locLeaves = map (first (+i₀)) . zip [0..] $ onlyLeaves locT-              findNeighbours :: (Int, x`WithAny`([Int+Needle x], y)) -> (y, Neighbourhood x)-              findNeighbours (i, WithAny (vns,y) x)-                         = (y, cullNeighbours locRieM-                                 (i, WithAny([ (i,v)-                                             | (i,WithAny _ xN) <- locLeaves-                                             , Just v <- [xN.-~.x] ]-                                                ++ aprioriNgbs)-                                             x))-               where aprioriNgbs :: [(Int, Needle x)]-                     aprioriNgbs = catMaybes-                                    [ (second $ const v) <$>-                                          positionIndex (pure locRieM) shd' xN-                                    | Right v <- vns-                                    , let xN = xi.+~^v :: x ]-                                 ++ [ (i,v) | Left i <- vns-                                            , Right (_,xN) <- [indexShadeTree shd' i]-                                            , Just v <- [xN.-~.x] ]-                     Just xi = toInterior x-              -              locRieM :: Metric x-              locRieM = case pointsCovers . catMaybes . map (toInterior . _topological)-                                  $ onlyLeaves locT of-                          [sh₀] -> metricf sh₀+fromTopShaded metricf shd = $notImplemented -cullNeighbours :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-      => Metric x -> (Int, x`WithAny`[(Int,Needle x)]) -> Neighbourhood x-cullNeighbours locRieM (i, WithAny vns x)-           = Neighbourhood (UArr.fromList . sort $ _nvectId<$>execState seek mempty)-                           locRieM- where seek :: State [NeighbourhoodVector x] ()-       seek = do-          Hask.forM_ ( fastNubBy (comparing fst) $ vns )-                    $ \(iNgb, v) ->-             when (iNgb/=i) `id`do-                oldNgbs <- get-                let w₀ = locRieM<$|v-                    l = sqrt $ w₀<.>^v-                    onOverlap = sum [ o^2 | nw<-oldNgbs-                                          , let o = (nw^.nvectNormal)<.>^v-                                          , o > 0 ]-                when (l > onOverlap) `id`do-                   let w = w₀^/sqrt l^3-                       newCandidates-                          = NeighbourhoodVector iNgb v w l 0-                          : [ ongb & otherNeighboursOverlap .~ 0-                            | ongb <- oldNgbs-                            , let o = w<.>^(ongb^.theNVect)-                                  newOverlap = (if o > 0 then (o^2+) else id)-                                                $ ongb^.otherNeighboursOverlap-                            , newOverlap < ongb^.nvectLength ]-                   put $ recalcOverlaps newCandidates-       recalcOverlaps [] = []-       recalcOverlaps (ngb:ngbs)-             = (ngb & otherNeighboursOverlap +~ furtherOvl)-             : recalcOverlaps [ ngb' & otherNeighboursOverlap +~ max 0 o ^ 2-                              | ngb' <- ngbs-                              , let o = (ngb^.nvectNormal)<.>^(ngb'^.theNVect) ]-        where furtherOvl = sum [ o^2 | nw<-ngbs-                                     , let o = (nw^.nvectNormal)<.>^(ngb^.theNVect)-                                     , o > 0 ]-                -- | Re-calculate the links in a web, so as to give each point a satisfyingly --   “complete-spanning” environment. smoothenWebTopology :: (WithField ℝ Manifold x, SimpleSpace (Needle x))              => MetricChoice x -> PointsWeb x y -> PointsWeb x y-smoothenWebTopology mc = swt- where swt (PointsWeb shd net) = PointsWeb shd . go allNodes Set.empty-                                                   . fst $ makeIndexLinksSymmetric net-        where allNodes = Set.fromList . Arr.toList $ fst <$> Arr.indexed net-              go activeSet pastLinks asd-                 | all (isNothing.fst) refined-                 , Set.null (Set.difference symmetryTouched pastLinks)-                               = Arr.imap finalise asd'-                 | otherwise   = go (Set.fromList-                                         [ j | (Just i, (_,Neighbourhood ngbs' _))-                                               <-refined-                                         , j <- i : UArr.toList ngbs' ]-                                      `Set.union` (Set.map fst symmetryTouched))-                                    updtLinks-                                    asd'-               where refined = reseek<$>Set.toList activeSet-                      where reseek i = ( guard isNews >> pure i-                                       , (y, Neighbourhood newNgbs locRieM) )-                             where isNews = newNgbs /= oldNgbs-                                             && or [ not $ Set.member (i,j) pastLinks-                                                   | j <- UArr.toList newNgbs ]-                                   (y,Neighbourhood oldNgbs locRieM) = asd Arr.! i-                                   nextNeighbours = fastNub-                                     $ UArr.toList oldNgbs-                                     ++ (UArr.toList._neighbours.snd.(asd Arr.!)-                                             =<< UArr.toList oldNgbs)-                                   x = xLookup Arr.! i-                                   Neighbourhood newNgbs _-                                     = cullNeighbours locRieM-                                        ( i, WithAny [ (j,v)-                                                     | j <- nextNeighbours-                                                     , Just v-                                                         <- [x .-~. xLookup Arr.! j] ]-                                                     x )-                     (asd', symmetryTouched) = makeIndexLinksSymmetric-                              $ asd Arr.// [(i,n) | (Just i,n) <- refined]-                     updtLinks = Set.unions-                                   [ pastLinks-                                   , Set.fromList-                                      [ (i,j) | (Just i,(_,Neighbourhood n _)) <- refined-                                              , j<-UArr.toList n ]-                                   , symmetryTouched ]-              finalise i (y, Neighbourhood n em)-                  = (y, cullNeighbours em (i, WithAny [ (j,v)-                                                      | j<-UArr.toList n-                                                      , let xN = xLookup Arr.! j-                                                      , Just v <- [xN.-~.x] ]-                                                      x ))-               where x = xLookup Arr.! i-              xLookup = Arr.fromList $ onlyLeaves shd+smoothenWebTopology = knitShortcuts -makeIndexLinksSymmetric-       :: Arr.Vector (y, Neighbourhood x)-       -> (Arr.Vector (y, Neighbourhood x), Set.Set (WebNodeId,WebNodeId))-makeIndexLinksSymmetric orig = runST (do-    result <- Arr.thaw orig-    touched <- newSTRef $ Set.empty-    (`Arr.imapM_`orig) $ \i (_,Neighbourhood ngbs _) -> do-       UArr.forM_ ngbs $ \j -> do-          (yn, Neighbourhood nngbs lsc) <- MArr.read result j-          when (not $ i`UArr.elem`nngbs) `id`do-             MArr.write result j (yn, Neighbourhood (UArr.snoc nngbs i) lsc)-             modifySTRef touched $ Set.insert (j,i)-    final <- Arr.freeze result-    allTouched <- readSTRef touched-    return (final, allTouched)-  ) -indexWeb :: (WithField ℝ Manifold x, SimpleSpace (Needle x))-                => PointsWeb x y -> WebNodeId -> Maybe (x,y)-indexWeb (PointsWeb rsc assocD) i-  | i>=0, i<Arr.length assocD-  , Right (_,x) <- indexShadeTree rsc i  = pure (x, fst (assocD Arr.! i))-  | otherwise                            = empty -unsafeIndexWebData :: PointsWeb x y -> WebNodeId -> y-unsafeIndexWebData (PointsWeb _ asd) i = fst (asd Arr.! i)+type OSNeedle x = (Needle' x, Needle x)+type OSNode x y = (OSNeedle x, WebLocally x y)+type CPCone x = (Needle' x, ℝ) -webEdges :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x))-            => PointsWeb x y -> [((x,y), (x,y))]-webEdges web@(PointsWeb rsc assoc) = (lookId***lookId) <$> toList allEdges- where allEdges :: Set.Set (WebNodeId,WebNodeId)-       allEdges = Hask.foldMap (\(i,(_, Neighbourhood ngbs _))-                    -> Set.fromList [(min i i', max i i')-                                    | i'<-UArr.toList ngbs ]-                               ) $ Arr.indexed assoc-       lookId i | Just xy <- indexWeb web i  = xy +-- | Consider at each node not just the connections to already known neighbours, but+--   also the connections to /their/ neighbours. If these next-neighbours turn out+--   to be actually situated closer, link to them directly.+knitShortcuts :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x))+             => MetricChoice x -> PointsWeb x y -> PointsWeb x y+knitShortcuts metricf w₀ = tweakWebGeometry metricf closeObtuseAngles+                             $ pseudoFixMaximise (rateLinkings w₀) w₀+ where pseudoFixMaximise oldBadness oldSt+         | newBadness < oldBadness  = pseudoFixMaximise newBadness newSt+         | otherwise                = newSt+        where newSt = tweakWebGeometry metricf pickNewNeighbours+                          $ bidirectionaliseWebLinks oldSt+              newBadness = rateLinkings newSt+       rateLinkings :: PointsWeb x y -> Double+       rateLinkings = geometricMeanOf rateNode . webLocalInfo+       rateNode :: WebLocally x y -> Double+       rateNode info = geometricMeanOf+             (\(_, (δx,_)) -> info^.nodeLocalScalarProduct|$|δx)+             $ info^.nodeNeighbours+       +       pickNewNeighbours :: WebLocally x y -> [WebNodeId]+       pickNewNeighbours me = fst <$> go Nothing [] candidates+        where go Nothing prevs (cs:ccs) = case bestNeighbours' lm' cs of+                        (links, Nothing) -> links+                        (links, Just newWall)+                         | Just _ <- me^.webBoundingPlane -> links+                         | otherwise  ->+                             links ++ go (Just newWall) ((snd<$>links) ++ prevs) ccs+              go (Just wall) prevs (cs:ccs) = case gatherGoodNeighbours+                               lm' lm wall prevs [] cs of+                        (links, Nothing) -> links+                        (links, Just newWall)+                         | Nothing <- me^.webBoundingPlane+                         , (_:_) <-ccs ->+                             links ++ go (Just newWall) ((snd<$>links) ++ prevs) ccs+                         | otherwise   -> links+              go _ _ [] = []+              lm' = me^.nodeLocalScalarProduct :: Metric x+              lm = dualNorm lm'+              candidates :: [[(WebNodeId, Needle x)]]+              candidates = preferred : other+               where (preferred, other) = case localOnion me [] of+                       _l₀:l₁:l₂:ls -> ( first _thisNodeId . swap <$> (l₁++l₂)+                                       , map (first _thisNodeId . swap) <$> ls )+                       [_l₀,l₁] -> (first _thisNodeId . swap <$> l₁, [])+       +       closeObtuseAngles :: WebLocally x y -> [WebNodeId]+       closeObtuseAngles me = go [ (dv,v) ^/ sqrt (dv<.>^v)+                                 | (i,(v,_)) <- me^.nodeNeighbours+                                 , let dv = metric<$|v ]+                                 candidates+        where go :: [OSNeedle x] -> [OSNode x y] -> [WebNodeId]+              go existing fillSrc = case constructUninhabitedCone existing of+                    Nothing -> fst <$> me^.nodeNeighbours+                    Just cone -> case findInCone cone fillSrc of+                      Just ((fv,filler),fillSrc')+                              -> (filler^.thisNodeId) : go (fv:existing) fillSrc'+                      Nothing -> fst <$> me^.nodeNeighbours+              constructUninhabitedCone :: [OSNeedle x] -> Maybe (CPCone x)+              constructUninhabitedCone vs = find (not.(`any`vs).includes)+                                              $ coneBetween <$> choices dimension vs+               where coneBetween :: [(Needle' x, a)] -> (Needle' x, ℝ)+                     coneBetween dvs = (coneDir, (coMetric|$|coneDir)/sqrt 2)+                      where coneDir = sumV $ fst <$> dvs+              findInCone :: CPCone x -> [OSNode x y]+                             -> Maybe (OSNode x y, [OSNode x y])+              findInCone cone ((po,pn):ps) | cone`includes`po  = Just ((po,pn), ps)+              findInCone (coneDir, _) ((po,pn):_)+                | Just wall <- pn^.webBoundingPlane+                , BoundarylessWitness <- boundarylessWitness :: BoundarylessWitness x+                , DualSpaceWitness <- dualSpaceWitness :: DualSpaceWitness (Needle x)+                , testp <- pn^.thisNodeCoord .+~^ (coMetric<$|wall)+                , (metric |$| testp.-~!me^.thisNodeCoord) > (metric|$|snd po)+                    = Nothing+              findInCone cone (p:ps) = second (p:) <$> findInCone cone ps+              findInCone _ [] = Nothing+              includes :: CPCone x -> OSNeedle x -> Bool+              (coneDir, narrowing)`includes`(_, v) = coneDir<.>^v >= narrowing+              candidates :: [OSNode x y]+              candidates = case localOnion me [] of+                       _l₀:_l₁:ls -> concat [ snd <$> sortBy (comparing fst)+                                                [ (distSq, ((dv,v) ^/ sqrt distSq, node))+                                                | (v, node) <- layer+                                                , let dv = metric<$|v+                                                      distSq = dv<.>^v ]+                                            | layer <- ls ]+                       _ -> []+              metric = me^.nodeLocalScalarProduct+              coMetric = dualNorm metric+       dimension = subbasisDimension (entireBasis :: SubBasis (Needle x)) -coerceWebDomain :: ∀ a b y . (Manifold a, Manifold b, LocallyCoercible a b)+choices :: Int -> [a] -> [[a]]+choices n l = go n l id []+ where go 0 _ f = (f[]:)+       go _ [] _ = id+       go n (x:xs) f = go n xs f . go (n-1) xs ((x:).f)++meanOf :: (Hask.Foldable f, Fractional n) => (a -> n) -> f a -> n+meanOf f = renormalise . Hask.foldl' accs (0, 0::Int)+ where renormalise (acc,n) = acc/fromIntegral n+       accs (acc,n) x = (acc+f x, succ n)++geometricMeanOf :: (Hask.Foldable f, Floating n) => (a -> n) -> f a -> n+geometricMeanOf f = exp . meanOf (log . f)++++webBoundary :: WithField ℝ Manifold x => PointsWeb x y -> [(Cutplane x, y)]+webBoundary = webLocalInfo >>> Hask.toList >>> Hask.concatMap`id`+        \info -> [ (Cutplane (info^.thisNodeCoord) (Stiefel1 wall), info^.thisNodeData)+                 | Just wall <- [info^.webBoundingPlane] ]+++coerceWebDomain :: ∀ a b y+     . (Manifold a, Manifold b, LocallyCoercible a b, SimpleSpace (Needle b))                                  => PointsWeb a y -> PointsWeb b y-coerceWebDomain (PointsWeb rsc assoc)-         = case oppositeLocalCoercion :: CanonicalDiffeomorphism b a of-   CanonicalDiffeomorphism-       -> PointsWeb ( coerceShadeTree rsc )-                    ( fmap (second $ localScalarProduct-                              %~transformNorm (arr $ coerceNeedle ([]::[(b,a)])))-                         assoc )+coerceWebDomain (PointsWeb web) = PointsWeb+     $ unsafeFmapTree ( fmap $ \(x, Neighbourhood y ngbs lscl bndry)+                            -> ( locallyTrivialDiffeomorphism x+                               , Neighbourhood y ngbs+                                       (coerceNorm ([]::[(a,b)]) lscl)+                                       (fmap crcNeedle' bndry) ) )+                      crcNeedle' coerceShade web+ where crcNeedle' = case ( dualSpaceWitness :: DualSpaceWitness (Needle a)+                         , dualSpaceWitness :: DualSpaceWitness (Needle b) ) of+           (DualSpaceWitness, DualSpaceWitness) -> arr $ coerceNeedle' ([]::[(a,b)])   data InterpolationIv y = InterpolationIv {@@ -441,7 +404,13 @@                         , Geodesic x, Geodesic y )                => PointsWeb x y -> Cutplane x -> [(x,y)] sliceWeb_lin web = sliceEdgs- where edgs = webEdges web+ where edgs :: [((x,y),(x,y))]+       edgs = [ (gnodes i₀, gnodes i₁)+              | (i₀,i₁) <- fastNub [ (i₀,i₁)+                                   | (il,ir) <- edges graph+                                   , let [i₀,i₁] = sort [il,ir] ]+              ]+       (graph, gnodes) = toGraph web        sliceEdgs cp = [ (xi d, yi d)  -- Brute-force search through all edges                       | ((x₀,y₀), (x₁,y₁)) <- edgs                       , Just d <- [cutPosBetween cp (x₀,x₁)]@@ -518,151 +487,189 @@        y₁ = maximum (snd<$>pts)        pts = fst . fst <$> toList (localFocusWeb web) -webLocalInfo :: ∀ x y . WithField ℝ Manifold x-            => PointsWeb x y -> PointsWeb x (WebLocally x y)-webLocalInfo origWeb = result- where result = wli $ localFocusWeb origWeb-       wli (PointsWeb rsc asd) = PointsWeb rsc asd'-        where asd' = Arr.imap localInfo asd-       localInfo i (((x,y), ngbCo), ngbH)-            = ( LocalWebInfo {-                  _thisNodeCoord = x-                , _thisNodeData = y-                , _thisNodeId = i-                , _nodeNeighbours = [ (iNgb, (δx, neighbour))-                                    | iNgb <- UArr.toList $ ngbH^.neighbours-                                    , let neighbour = unsafeIndexWebData result iNgb-                                          Just δx = _thisNodeCoord neighbour.-~.x-                                    ]-                , _nodeLocalScalarProduct = ngbH^.localScalarProduct-                , _nodeIsOnBoundary = anyUnopposed (ngbH^.localScalarProduct) ngbCo-                }, ngbH )-       anyUnopposed rieM ngbCo = (`any`ngbCo) $ \(v,_)-                         -> not $ (`any`ngbCo) $ \(v',_)-                              -> (rieM<$|v) <.>^ v' < 0 ++hardbakeChunk :: WebChunk x y -> PointsWeb x y+hardbakeChunk = _thisChunk++entireWeb :: PointsWeb x y -> WebChunk x y+entireWeb web = WebChunk web []+ localFocusWeb :: WithField ℝ Manifold x                    => PointsWeb x y -> PointsWeb x ((x,y), [(Needle x, y)])-localFocusWeb (PointsWeb rsc asd) = PointsWeb rsc asd''- where asd' = Arr.imap (\i (y,n) -> case indexShadeTree rsc i of-                                         Right (_,x) -> ((x,y),n) ) asd-       asd''= Arr.map (\((x,y),n) ->-                       (((x,y), [ ( case x'.-~.x of-                                     Just v -> v-                                  , y')-                                | j<-UArr.toList (n^.neighbours)-                                , let ((x',y'),_) = asd' Arr.! j-                                ]), n)-                 ) asd'+localFocusWeb = webLocalInfo >>> fmap `id`\n+           -> ( (n^.thisNodeCoord, n^.thisNodeData)+              , [ (δx, ngb^.thisNodeData)+                | (_, (δx, ngb)) <- n^.nodeNeighbours ] ) +++treewiseTraverseLocalWeb :: ∀ f x y . (WithField ℝ Manifold x, Hask.Applicative f)+     => (WebLocally x y -> f y)+       -> (∀ t i w . (Hask.Traversable t, Ord i) => (w -> f w) -> t (i, w) -> f (t w) )+       -> PointsWeb x y -> f (PointsWeb x y)+treewiseTraverseLocalWeb fl ct = fmap hardbakeChunk . twt . entireWeb+ where twt = treewiseTraverseLocalWeb' fl $ ct twt++treewiseTraverseLocalWeb' :: ∀ f x y . (WithField ℝ Manifold x, Hask.Applicative f)+     => (WebLocally x y -> f y)+       -> (NonEmpty (Int, WebChunk x y) -> f (NonEmpty (WebChunk x y)))+       -> WebChunk x y -> f (WebChunk x y)+treewiseTraverseLocalWeb' fl ct domain+                  = $notImplemented{-+ where rezoomed (PlainLeaves _) _ = localTraverseWebChunk fl domain+       rezoomed tree pos+         | pos == i₀, nLeaves tree == lDomain+             = fmap reassemble $ ct (NE.zipWith+                       (\jb (i₀b,t')+                         -> (jb, domain & overrideStart .~ i₀+i₀b+                                        & overriddenData+                                            .~ Arr.slice i₀b (nLeaves t') domainData ))+                       (0:|[1..]) branches)+         | otherwise                     = go branches+        where branches = trunkBranches tree+              go (_:|((i₀nb,nb):brs))+                | pos+i₀nb <= i₀  = go $ (i₀nb,nb):|brs+              go ((i₀b,t):|_) = rezoomed t (pos+i₀b)+              reassemble :: NonEmpty (WebChunk x y) -> WebChunk x y+              reassemble brs = domain & overriddenData+                                  .~ Hask.foldMap _overriddenData brs+       lDomain = Arr.length domainData+   -}+++ localOnion :: ∀ x y . WithField ℝ Manifold x-            => WebLocally x y -> [[WebLocally x y]]-localOnion origin = go Map.empty $ Map.singleton (origin^.thisNodeId) (1, origin)- where go previous next+            => WebLocally x y -> [WebNodeId] -> [[(Needle x, WebLocally x y)]]+localOnion origin directCandidates = map sortBCDistance . go Map.empty . Map.fromList+                      $ (origin^.thisNodeId, (1, origin)) : seeds+ where seeds :: [(WebNodeId, (Int, WebLocally x y))]+       seeds = [ (nid, (1, ninfo))+               | nid <- directCandidates+               , (_,(_,ninfo)) <- origin^.nodeNeighbours+               , ninfo^.thisNodeId == nid ]+       go previous next         | Map.null next = []-        | otherwise  = ( snd <$> sortBy (comparing $ negate . fst)+        | otherwise  = ( computeOffset . snd+                                    <$> sortBy (comparing $ negate . fst)                                                  (Hask.toList next) )                      : go (Map.union previous next)                           (Map.fromListWith (\(n,ninfo) (n',_) -> (n+n'::Int, ninfo))-                                [ (nnid,(1,nneigh))-                                | (nid,(_,ninfo))<-Map.toList next-                                , (nnid,(_,nneigh))<-ninfo^.nodeNeighbours-                                , Map.notMember nnid previous ])+                             [ (nnid,(1,nneigh))+                             | (nid,(_,ninfo))<-Map.toList next+                             , (nnid,(_,nneigh))<-ninfo^.nodeNeighbours+                             , Map.notMember nnid previous && Map.notMember nnid next ])+       computeOffset p = case p^.thisNodeCoord .-~. origin^.thisNodeCoord of+                Just v -> (v,p)+       sortBCDistance = map snd . sortBy (comparing fst) . map (bcDist&&&id)+        where bcDist (v,_)+                = normSq (origin^.nodeLocalScalarProduct) $ v^-^seedBarycenterOffs+       seedBarycenterOffs = sumV ngbOffs ^/ fromIntegral (length directCandidates + 1)+        where ngbOffs = [ v | (_, (_, n)) <- seeds+                            , let Just v = n^.thisNodeCoord .-~. origin^.thisNodeCoord ]  webOnions :: ∀ x y . WithField ℝ Manifold x             => PointsWeb x y -> PointsWeb x [[(x,y)]]-webOnions = localFmapWeb (map (map $ _thisNodeCoord&&&_thisNodeData) . localOnion)+webOnions = localFmapWeb (map (map $ _thisNodeCoord&&&_thisNodeData <<< snd)+                                . (`localOnion`[])) -nearestNeighbour :: (WithField ℝ Manifold x, SimpleSpace (Needle x))+nearestNeighbour :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x))                       => PointsWeb x y -> x -> Maybe (x,y)-nearestNeighbour (PointsWeb rsc asd) x = fmap lkBest $ positionIndex empty rsc x- where lkBest (iEst, (_, xEst)) = (xProx, yProx)-        where (iProx, (xProx, _)) = minimumBy (comparing $ snd . snd)-                                     $ (iEst, (xEst, normSq locMetr vEst))-                                         : neighbours-              (yProx, _) = asd Arr.! iProx-              (_, Neighbourhood neighbourIds locMetr) = asd Arr.! iEst-              neighbours = [ (i, (xNgb, normSq locMetr v))-                           | i <- UArr.toList neighbourIds-                           , let Right (_, xNgb) = indexShadeTree rsc i-                                 Just v = xNgb.-~.x-                           ]-              Just vEst = xEst.-~.x+nearestNeighbour = webLocalInfo >>> \(PointsWeb rsc) x+                 -> fmap (fine x) (positionIndex empty rsc x)+ where fine :: x -> (Int, ( [Shaded x (Neighbourhood x (WebLocally x y))]+                          , (x, Neighbourhood x (WebLocally x y)) ))+                 -> (x,y)+       fine x (_, (_, (xc, (Neighbourhood c _ locMetr _))))+           = snd . minimumBy (comparing fst)+              . map (first $ (c^.nodeLocalScalarProduct|$|)+                           . (^-^vc))+              $ (zeroV, (xc, c^.thisNodeData))+                : [ (δx, (ngb^.thisNodeCoord, ngb^.thisNodeData))+                  | (_, (δx, ngb)) <- c^.nodeNeighbours ]+        where Just vc = x.-~.xc   -instance Hask.Functor (WebLocally x) where-  fmap f (LocalWebInfo co dt id ng sp bn)-       = LocalWebInfo co (f dt) id (map (second . second $ fmap f) ng) sp bn-instance WithField ℝ Manifold x => Comonad (WebLocally x) where-  extract = _thisNodeData-  extend f this@(LocalWebInfo co _ id ng sp bn)-      = LocalWebInfo co (f this) id (map (second . second $ extend f) ng) sp bn-  duplicate this@(LocalWebInfo co _ id ng sp bn)-      = LocalWebInfo co this id (map (second $ second duplicate) ng) sp bn- -- ^ 'fmap' from the co-Kleisli category of 'WebLocally'.-localFmapWeb :: WithField ℝ Manifold x-                => (WebLocally x y -> z) -> PointsWeb x y -> PointsWeb x z-localFmapWeb f = webLocalInfo >>> fmap f+localTraverseWeb :: (WithField ℝ Manifold x, Hask.Applicative m)+                => (WebLocally x y -> m z) -> PointsWeb x y -> m (PointsWeb x z)+localTraverseWeb f = webLocalInfo >>> Hask.traverse f -traverseWebWithStrategy :: ( WithField ℝ Manifold x, Hask.Applicative m )-               => InconsistencyStrategy m x y -> (WebLocally x y -> Maybe y)-                     -> PointsWeb x y -> m (PointsWeb x y)-traverseWebWithStrategy strat f = webLocalInfo-               >>> traverse (\info -> handleInconsistency strat-                                       (info^.thisNodeData) (f info))+-- ^ 'fmap' from the co-Kleisli category of 'WebLocally', restricted to some+--   contiguous part of a web.+localTraverseWebChunk :: (WithField ℝ Manifold x, Hask.Applicative m)+                => (WebLocally x y -> m y) -> WebChunk x y -> m (WebChunk x y)+localTraverseWebChunk f (WebChunk this outlayers)+      = fmap (\c -> WebChunk c outlayers) $ localTraverseWeb f this  differentiateUncertainWebLocally :: ∀ x y-   . ( WithField ℝ Manifold x, SimpleSpace (Needle x)-     , WithField ℝ Refinable y, SimpleSpace (Needle y) )+   . ( ModellableRelation x y )             => WebLocally x (Shade' y)              -> Shade' (LocalLinear x y)-differentiateUncertainWebLocally info-          = case estimateLocalJacobian-                          (info^.nodeLocalScalarProduct)-                          [ ( Local δx :: Local x, ngb^.thisNodeData )-                          | (δx,ngb) <- (zeroV, info)-                                      : (snd<$>info^.nodeNeighbours)-                          ] of-               Just j -> j-               _      -> Shade' zeroV mempty+differentiateUncertainWebLocally = duwl+                ( dualSpaceWitness :: DualSpaceWitness (Needle x)+                , dualSpaceWitness :: DualSpaceWitness (Needle y) )+ where duwl (DualSpaceWitness, DualSpaceWitness) info+          = case fitLocally $+                          (\(δx,ngb) -> (δx, ngb^.thisNodeData) )+                          <$> (zeroV,info) : envi+                          of+               Just (AffineModel _ j :: AffineModel x y) -> dualShade j+        where _:directEnvi:remoteEnvi = localOnion info []+              envi = directEnvi ++ concat remoteEnvi + differentiateUncertainWebFunction :: ∀ x y-   . ( WithField ℝ Manifold x, SimpleSpace (Needle x)-     , WithField ℝ Manifold y, SimpleSpace (Needle y), Refinable y )+   . ( ModellableRelation x y )             => PointsWeb x (Shade' y)              -> PointsWeb x (Shade' (LocalLinear x y)) differentiateUncertainWebFunction = localFmapWeb differentiateUncertainWebLocally  differentiate²UncertainWebLocally :: ∀ x y-   . ( WithField ℝ Manifold x, FlatSpace (Needle x)-     , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y) )+   . ( ModellableRelation x y )             => WebLocally x (Shade' y)              -> Shade' (Needle x ⊗〃+> Needle y) differentiate²UncertainWebLocally = d²uwl-                ( pseudoAffineWitness :: PseudoAffineWitness x-                , pseudoAffineWitness :: PseudoAffineWitness y-                , dualSpaceWitness :: DualSpaceWitness (Needle x)+                ( dualSpaceWitness :: DualSpaceWitness (Needle x)                 , dualSpaceWitness :: DualSpaceWitness (Needle y) )- where d²uwl ( PseudoAffineWitness (SemimanifoldWitness _)-             , PseudoAffineWitness (SemimanifoldWitness _)-             , DualSpaceWitness, DualSpaceWitness ) info-          = case estimateLocalHessian $-                          (\ngb -> case (ngb^.thisNodeCoord .-~. info^.thisNodeCoord) of-                             Just δx -> (Local δx :: Local x, ngb^.thisNodeData) )-                          <$> info :| envi+ where d²uwl (DualSpaceWitness, DualSpaceWitness) info+          = case fitLocally $+                          (\(δx,ngb) -> (δx, ngb^.thisNodeData) )+                          <$> (zeroV,info) : envi                           of-               QuadraticModel _ h -> dualShade $ projectShade-                          (fromEmbedProject (acoSnd.acoSnd ^/ 2)-                                            (snd.snd ^* 2) ) h-        where xVol :: SymmetricTensor ℝ (Needle x)-              xVol = squareVs $ fst.snd<$>info^.nodeNeighbours-              _:directEnvi:remoteEnvi = localOnion info-              envi = directEnvi ++ take (nMinData - length directEnvi) (concat remoteEnvi)-       nMinData = 1 + regular_neighboursCount-                         (subbasisDimension (entireBasis :: SubBasis (Needle x)))+               Just (QuadraticModel _ _ h :: QuadraticModel x y)+                        -> linIsoTransformShade (2*^id) $ dualShade h+        where _:directEnvi:remoteEnvi = localOnion info []+              envi = directEnvi ++ concat remoteEnvi ++-- | Calculate a quadratic fit with uncertainty margin centered around the connection+--   between any two adjacent nodes. In case of a regular grid (which we by no means+--   require here!) this corresponds to the vector quantities of an Arakawa type C/D+--   grid (cf. A. Arakawa, V.R. Lamb (1977):+--   Computational design of the basic dynamical processes of the UCLA general circulation model)+localModels_CGrid :: ∀ x y ㄇ . ( ModellableRelation x y, LocalModel ㄇ )+          => PointsWeb x (Shade' y) -> [(x, ㄇ x y)]+localModels_CGrid = Hask.concatMap theCGrid . Hask.toList . webLocalInfo+ where theCGrid :: WebLocally x (Shade' y) -> [(x, ㄇ x y)]+       theCGrid node = [ ( pn .-~^ δx^/2+                         , propagationCenteredModel+                             ( LocalDataPropPlan+                                    pn+                                    (negateV δx)+                                    (ngbNode^.thisNodeData)+                                    (node^.thisNodeData)+                                    (fmap (second _thisNodeData)+                                      . concat . tail+                                           $ localOnion ngbNode [node^.thisNodeId] )+                                          ) )+                       | (nid, (δx, ngbNode)) <- node^.nodeNeighbours+                       , nid > node^.thisNodeId+                       , Just pn <- [toInterior $ ngbNode^.thisNodeCoord]+                       ]++ acoSnd :: ∀ s v y . ( Object (Affine s) y, Object (Affine s) v                     , LinearSpace v, Scalar v ~ s ) => Affine s y (v,y) acoSnd = case ( linearManifoldWitness :: LinearManifoldWitness v@@ -671,54 +678,37 @@    (LinearManifoldWitness BoundarylessWitness, DualSpaceWitness, DualSpaceWitness)        -> const zeroV &&& id --- | Heuristic formula, matches the number of neighbours each vertex has in a one----   and two-dimensional count-regular_neighboursCount :: Int -> Int-regular_neighboursCount d- | d>0        = (regular_neighboursCount (d-1) + 1)*2- | otherwise  = 0 - differentiate²UncertainWebFunction :: ∀ x y-   . ( WithField ℝ Manifold x, FlatSpace (Needle x)-     , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y) )+   . ( ModellableRelation x y )          => PointsWeb x (Shade' y)           -> PointsWeb x (Shade' (Needle x ⊗〃+> Needle y))  differentiate²UncertainWebFunction = localFmapWeb differentiate²UncertainWebLocally -rescanPDELocally :: ∀ x y ð .-     ( WithField ℝ Manifold x, FlatSpace (Needle x)-     , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y) )-         => DifferentialEqn x ð y -> WebLocally x (Shade' y)-                                -> (Maybe (Shade' y), Maybe (Shade' ð))+rescanPDELocally :: ∀ x y ㄇ .+     ( ModellableRelation x y, LocalModel ㄇ )+         => DifferentialEqn ㄇ x y -> WebLocally x (Shade' y) -> Maybe (Shade' y) rescanPDELocally = case ( dualSpaceWitness :: DualNeedleWitness x                         , dualSpaceWitness :: DualNeedleWitness y                         , boundarylessWitness :: BoundarylessWitness x                         , pseudoAffineWitness :: PseudoAffineWitness y ) of    ( DualSpaceWitness,DualSpaceWitness,BoundarylessWitness     , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) )-     -> \f info -> let xc = info^.thisNodeCoord+     -> \f info+          -> if isJust $ info^.webBoundingPlane+              then return $ info^.thisNodeData+              else let xc = info^.thisNodeCoord                        yc = info^.thisNodeData.shadeCtr                    in case f $ coverAllAround (xc, yc)                                      [ (δx, (ngb^.thisNodeData.shadeCtr.-~!yc) ^+^ v)                                      | (_,(δx,ngb))<-info^.nodeNeighbours                                      , v <- normSpanningSystem'                                               (ngb^.thisNodeData.shadeNarrowness)] of-                        LocalDifferentialEqn _ rescan-                            -> rescan (info^.thisNodeData)-                                      (differentiateUncertainWebLocally info)-                                      (differentiate²UncertainWebLocally info)--rescanPDEOnWeb :: ( WithField ℝ Manifold x, FlatSpace (Needle x)-                  , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y)-                  , Hask.Applicative m )-                => InconsistencyStrategy m x (Shade' y, Shade' ð)-                  -> DifferentialEqn x ð y -> PointsWeb x (Shade' y)-                                   -> m (PointsWeb x (Shade' y, Shade' ð))-rescanPDEOnWeb strat deq = traverseWebWithStrategy strat-                 (fzip . rescanPDELocally deq . fmap fst)-         . fmap (\shy -> (shy, error-                   "No default value for inconsistent PDE-rescanning on web"))+                        LocalDifferentialEqn rescan -> fst+                             ( rescan $ case fitLocally $ map (id *** _thisNodeData)+                                               =<< (localOnion info []) of+                                 Just ㄇ -> ㄇ)+                                 >>= intersectShade's . (:|[info^.thisNodeData])  toGraph :: (WithField ℝ Manifold x, SimpleSpace (Needle x))               => PointsWeb x y -> (Graph, Vertex -> (x, y))@@ -726,8 +716,8 @@                 (graphFromEdges' edgs)  where edgs :: [(Int, Int, [Int])]        edgs = Arr.toList-            . Arr.imap (\i (_, Neighbourhood ngbs _) -> (i, i, UArr.toList ngbs))-                    $ webNodeAssocData wb+            . Arr.imap (\i (Neighbourhood _ ngbs _ _) -> (i, i, (i+) <$> UArr.toList ngbs))+            . Arr.fromList . Hask.toList $ webNodeRsc wb   @@ -798,6 +788,14 @@                Just r  -> pure r                Nothing -> throwE $ PropagationInconsistency n o ) +postponeInconsistencies :: Hask.Monad m => (NonEmpty υ -> Maybe υ)+   -> InformationMergeStrategy [] (WriterT [PropagationInconsistency x υ] m)+                                  (x,υ) υ+postponeInconsistencies merge = InformationMergeStrategy+           (\o n -> case merge $ o :| fmap snd n of+               Just r  -> pure r+               Nothing -> writer (o,[PropagationInconsistency n o]) )+ maybeAlt :: Hask.Alternative f => Maybe a -> f a maybeAlt (Just x) = pure x maybeAlt Nothing = Hask.empty@@ -809,14 +807,10 @@ deriving instance Hask.Functor (InconsistencyStrategy m x)  -iterateFilterDEqn_static :: ( WithField ℝ Manifold x, FlatSpace (Needle x)-                            , Refinable y, Geodesic y, FlatSpace (Needle y)-                            , WithField ℝ AffineManifold ð, Geodesic ð-                            , SimpleSpace (Needle ð)-                            , Hask.MonadPlus m )+iterateFilterDEqn_static :: ( ModellableRelation x y, Hask.MonadPlus m, LocalModel ㄇ )        => InformationMergeStrategy [] m (x,Shade' y) iy            -> Embedding (->) (Shade' y) iy-           -> DifferentialEqn x ð y+           -> DifferentialEqn ㄇ x y                  -> PointsWeb x (Shade' y) -> Cofree m (PointsWeb x (Shade' y)) iterateFilterDEqn_static strategy shading f                            = fmap (fmap (shading >-$))@@ -824,42 +818,154 @@                            . fmap (shading $->)  -filterDEqnSolutions_static :: ∀ x y iy ð m .-                              ( WithField ℝ Manifold x, FlatSpace (Needle x)-                              , Refinable y, Geodesic y, FlatSpace (Needle y)-                              , WithField ℝ AffineManifold ð, Geodesic ð-                              , SimpleSpace (Needle ð)-                              , Hask.MonadPlus m )+iterateFilterDEqn_static_selective :: ( ModellableRelation x y+                                      , Hask.MonadPlus m, badness ~ ℝ+                                      , LocalModel ㄇ )+       => InformationMergeStrategy [] m (x,Shade' y) iy+           -> Embedding (->) (Shade' y) iy+           -> (x -> iy -> badness)+           -> DifferentialEqn ㄇ x y+                 -> PointsWeb x (Shade' y) -> Cofree m (PointsWeb x (Shade' y))+iterateFilterDEqn_static_selective strategy shading badness f+      = fmap (fmap (shading >-$))+      . coiter (filterDEqnSolutions_static_selective strategy shading badness f)+      . fmap (shading $->)+++filterDEqnSolutions_static :: ∀ x y ㄇ iy m .+                     ( ModellableRelation x y, Hask.MonadPlus m, LocalModel ㄇ )        => InformationMergeStrategy [] m  (x,Shade' y) iy -> Embedding (->) (Shade' y) iy-          -> DifferentialEqn x ð y -> PointsWeb x iy -> m (PointsWeb x iy)+          -> DifferentialEqn ㄇ x y -> PointsWeb x iy -> m (PointsWeb x iy) filterDEqnSolutions_static = case geodesicWitness :: GeodesicWitness y of    GeodesicWitness _ -> \strategy shading f        -> webLocalInfo            >>> fmap (id &&& rescanPDELocally f . fmap (shading>-$))            >>> localFocusWeb >>> Hask.traverse ( \((_,(me,updShy)), ngbs)           -> let oldValue = me^.thisNodeData :: iy-             in  case updShy of-              (Just shy, Just shð) -> case ngbs of+             in if isJust $ me ^. webBoundingPlane+                 then return oldValue+                 else case updShy of+              Just shy -> case ngbs of                   []  -> pure oldValue                   _:_ | BoundarylessWitness <- (boundarylessWitness::BoundarylessWitness x)-                    -> maybeAlt-                          ( sequenceA [ fzip sj-                                >>= \ngbShyð -> (ngbInfo^.thisNodeCoord,)<$>+                    -> sequenceA [ maybeAlt sj+                                >>= \ngbShyð -> fmap ((me^.thisNodeCoord .+~^ δx,)+                                                   . (shading>-$))+                                  . mergeInformation strategy oldValue . Hask.toList+                                  $ (ngbInfo^.thisNodeCoord,)<$>                                      propagateDEqnSolution_loc                                        f (LocalDataPropPlan                                              (ngbInfo^.thisNodeCoord)                                              (negateV δx)                                              ngbShyð-                                             (shy, shð)-                                             (fmap (second ((shading>-$) . _thisNodeData)-                                                    . snd) $ ngbInfo^.nodeNeighbours)+                                             shy+                                             (fmap (second ((shading>-$) . _thisNodeData))+                                               . concat . tail $ localOnion ngbInfo+                                                                     [me^.thisNodeId])                                           )                                   | (δx, (ngbInfo,sj)) <- ngbs-                                  ] )+                                  ]                             >>= mergeInformation strategy (shading$->shy)               _ -> mergeInformation strategy oldValue empty         ) +++data Average a = Average { weight :: Int+                         , averageAcc :: a+                         } deriving (Hask.Functor)+instance Num a => Monoid (Average a) where+  mempty = Average 0 0+  mappend (Average w₀ a₀) (Average w₁ a₁) = Average (w₀+w₁) (a₀+a₁)+instance Hask.Applicative Average where+  pure = Average 1+  Average w₀ a₀ <*> Average w₁ a₁ = Average (w₀*w₁) (a₀ a₁)++average :: Fractional a => Average a -> a+average (Average w a) = a / fromIntegral w++averaging :: VectorSpace a => [a] -> Average a+averaging l = Average (length l) (sumV l)++filterDEqnSolutions_static_selective :: ∀ x y ㄇ iy m badness .+                              ( ModellableRelation x y+                              , Hask.MonadPlus m, badness ~ ℝ+                              , LocalModel ㄇ )+       => InformationMergeStrategy [] m  (x,Shade' y) iy -> Embedding (->) (Shade' y) iy+          -> (x -> iy -> badness)+          -> DifferentialEqn ㄇ x y+          -> PointsWeb x iy -> m (PointsWeb x iy)+filterDEqnSolutions_static_selective = case geodesicWitness :: GeodesicWitness y of+   GeodesicWitness _ -> \strategy shading badness f+       ->  -- Integration step: determine at each point from the function values+           -- what the derivatives should be, and use them to propagate the solution+           -- in all directions. We only spend a single computation step on regions+           -- where nothing much changes (indicating the a-priori information is+           -- too weak yet to make any predictions anyway), but multiple steps in+           -- regions where good progress is noticed.+         fmap fst . (runWriterT :: WriterT (Average badness) m (PointsWeb x iy)+                                        -> m (PointsWeb x iy, Average badness))+         . treewiseTraverseLocalWeb ( \me+          -> let oldValue = me^.thisNodeData :: iy+                 badHere = badness $ me^.thisNodeCoord+                 oldBadness = badHere oldValue+             in if isJust $ me ^. webBoundingPlane+                 then return oldValue+                 else case me^.nodeNeighbours of+                  [] -> pure oldValue+                  _:_ | BoundarylessWitness <- (boundarylessWitness::BoundarylessWitness x)+                    -> WriterT . fmap (\updated+                                    -> (updated, pure (oldBadness / badHere updated)))+                       $ sequenceA [ fmap ((me^.thisNodeCoord .+~^ δx,)+                                                   . (shading>-$))+                                  . mergeInformation strategy oldValue . Hask.toList+                                  $ (ngbInfo^.thisNodeCoord,)<$>+                                     propagateDEqnSolution_loc+                                       f (LocalDataPropPlan+                                             (ngbInfo^.thisNodeCoord)+                                             (negateV δx)+                                             (shading >-$ ngbInfo^.thisNodeData)+                                             (shading >-$ oldValue)+                                             (fmap (second ((shading>-$) . _thisNodeData))+                                               . concat . tail $ localOnion+                                                        ngbInfo [me^.thisNodeId] )+                                          )+                                  | (_, (δx, ngbInfo)) <- me^.nodeNeighbours+                                  ]+                            >>= mergeInformation strategy oldValue )+                 (\combiner branchData -> WriterT $ do+                       (branchResults,improvements)+                         <- runWriterT $ Hask.traverse+                                          (\(i,branch) -> fmap (i,)+                                                          . censor (pure . (i,) . average)+                                                          $ combiner branch)+                                          branchData+                       let (best, _) = maximumBy (comparing snd) improvements+                       (branchResults',improvements')+                         <- runWriterT $ Hask.traverse+                                          (\(i,branch) -> if i==best+                                             then censor (pure . (i,) . average)+                                                              $ combiner branch+                                             else WriterT $ return (branch, pure (i,1)) )+                                          branchResults+                       return ( branchResults'+                              , liftA2 (*) (averaging $ snd<$>improvements)+                                           (averaging $ snd<$>improvements') )+                 )+          >=> -- Boundary-condition / differentiation step: update the local values+              -- based on a-priori boundary conditions, possibly dependent on+              -- numerical derivatives of the current solution estimate.+              localTraverseWeb (\me -> fmap (shading$->)+                                         . maybeAlt $ rescanPDELocally f me)+            . fmap (shading>-$)++-- | The <http://hackage.haskell.org/package/transformers-0.5.4.0/docs/Control-Monad-Trans-Writer-Lazy.html#v:censor transformers version of this>+--   is insufficiently polymorphic, requiring @w ~ w'@.+censor :: Functor m (->) (->) => (w -> w') -> WriterT w m a -> WriterT w' m a+censor = mapWriterT . fmap . second+++ handleInconsistency :: InconsistencyStrategy m x a -> a -> Maybe a -> m a handleInconsistency AbortOnInconsistency _ i = i handleInconsistency IgnoreInconsistencies _ (Just v) = Identity v@@ -887,14 +993,13 @@ oldAndNew' (_, l) = (False,) <$> l  -filterDEqnSolutions_adaptive :: ∀ x y ð badness m-        . ( WithField ℝ Manifold x, SimpleSpace (Needle x)-          , WithField ℝ AffineManifold y, Refinable y, Geodesic y-          , WithField ℝ AffineManifold ð, Geodesic ð, SimpleSpace (Needle ð)-          , badness ~ ℝ, Hask.Monad m )+filterDEqnSolutions_adaptive :: ∀ x y ㄇ ð badness m+        . ( ModellableRelation x y, AffineManifold y+          , badness ~ ℝ, Hask.Monad m+          , LocalModel ㄇ )        => MetricChoice x      -- ^ Scalar product on the domain, for regularising the web.        -> InconsistencyStrategy m x (Shade' y)-       -> DifferentialEqn x ð y+       -> DifferentialEqn ㄇ x y        -> (x -> Shade' y -> badness)              -> PointsWeb x (SolverNodeState x y)                         -> m (PointsWeb x (SolverNodeState x y))@@ -905,7 +1010,7 @@                       -> m (PointsWeb x ( (WebLocally x (SolverNodeState x y)                                         , [(Shade' y, badness)]) ))        tryPreproc BoundarylessWitness (GeodesicWitness _)-               = traverse addPropagation $ webLocalInfo oldState+               = Hask.traverse addPropagation $ webLocalInfo oldState         where addPropagation wl                  | null neighbourInfo = pure (wl, [])                  | otherwise           = (wl,) . map (id&&&badness undefined)@@ -917,9 +1022,8 @@                                            (neigh^.thisNodeCoord)                                            (negateV δx)                                            (convexSetHull $ neigh^.thisNodeData-                                                                  .solverNodeStatus-                                           , undefined)-                                           (thisShy, undefined)+                                                                  .solverNodeStatus)+                                           thisShy                                            [ second (convexSetHull                                                      . _solverNodeStatus . _thisNodeData) nn                                            | (_,nn)<-neigh^.nodeNeighbours ] )@@ -1013,9 +1117,8 @@                                                       (n^.thisNodeCoord)                                                       (stepV ^-^ δx)                                                       (convexSetHull $-                                                        n^.thisNodeData.solverNodeStatus-                                                      , undefined)-                                                      (aprioriInterpolate, undefined)+                                                        n^.thisNodeData.solverNodeStatus)+                                                      aprioriInterpolate                                                       (second (convexSetHull                                                                ._solverNodeStatus                                                                ._thisNodeData)@@ -1067,8 +1170,7 @@                                             , let Just vOld = ngb^.thisNodeCoord .-~. xOld                                             ]                               -recomputeJacobian :: ( WithField ℝ Manifold x, SimpleSpace (Needle x)-                     , WithField ℝ Manifold y, SimpleSpace (Needle y), Refinable y )+recomputeJacobian :: ( ModellableRelation x y )              => PointsWeb x (SolverNodeState x y)              -> PointsWeb x (SolverNodeState x y) recomputeJacobian = webLocalInfo@@ -1079,12 +1181,11 @@   iterateFilterDEqn_adaptive-     :: ( WithField ℝ Manifold x, SimpleSpace (Needle x)-        , WithField ℝ AffineManifold y, Refinable y, Geodesic y, Hask.Monad m-        , WithField ℝ AffineManifold ð, Geodesic ð, SimpleSpace (Needle ð) )+     :: ( ModellableRelation x y, AffineManifold y+        , LocalModel ㄇ, Hask.Monad m )        => MetricChoice x      -- ^ Scalar product on the domain, for regularising the web.        -> InconsistencyStrategy m x (Shade' y)-       -> DifferentialEqn x ð y+       -> DifferentialEqn ㄇ x y        -> (x -> Shade' y -> ℝ) -- ^ Badness function for local results.              -> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)] iterateFilterDEqn_adaptive mf strategy f badness

@@ -0,0 +1,616 @@+-- |+-- Module      : Data.Manifold.Web.Internal+-- Copyright   : (c) Justus Sagemüller 2017+-- License     : GPL v3+-- +-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de+-- Stability   : experimental+-- Portability : portable+-- +{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE UndecidableInstances       #-}+{-# LANGUAGE StandaloneDeriving         #-}+{-# LANGUAGE DeriveGeneric              #-}+{-# LANGUAGE DeriveFunctor              #-}+{-# LANGUAGE DeriveFoldable             #-}+{-# LANGUAGE DeriveTraversable          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE GADTs                      #-}+{-# LANGUAGE Rank2Types                 #-}+{-# LANGUAGE MultiParamTypeClasses      #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE TypeOperators              #-}+{-# LANGUAGE TupleSections              #-}+{-# LANGUAGE KindSignatures             #-}+{-# LANGUAGE TemplateHaskell            #-}+{-# LANGUAGE UnicodeSyntax              #-}+++module Data.Manifold.Web.Internal where+++import Prelude hiding ((^))++import qualified Data.Vector.Unboxed as UArr++import Data.Manifold.Types+import Data.Manifold.Types.Primitive+import Data.Manifold.PseudoAffine+import Data.Manifold.Shade+import Data.Manifold.TreeCover+import Data.Function.Affine+import Data.VectorSpace (Scalar, (^+^), (^/), (^*), sumV)+import Math.LinearMap.Category ( SimpleSpace, LSpace, DualVector, Norm, Variance+                               , (<.>^), dualNorm, (<$|), (|$|), normSq+                               , dualSpaceWitness, DualSpaceWitness(..)+                               , FiniteDimensional (..) )+    +import qualified Data.Foldable       as Hask+import qualified Data.Traversable as Hask+import Data.List (sortBy)+import qualified Data.List.NonEmpty as NE+import Data.List.NonEmpty (NonEmpty((:|)))+import qualified Data.Map as Map+import qualified Data.Foldable.Constrained as CCt+import Data.Functor.Identity+import Data.Function ((&))+import Data.Ord (comparing)+import Data.List.FastNub (fastNub)+import qualified Data.IntSet as ℤSet+import Data.IntSet (IntSet)+import Data.Maybe (isNothing)+import Control.Arrow+import Control.Monad (guard, forM_)+import Control.Comonad+import Control.Monad.Trans.State++import Control.DeepSeq++import GHC.Generics (Generic)++import Control.Lens+import Control.Lens.TH++import Data.CallStack (HasCallStack)+++type WebNodeId = Int+type WebNodeIdOffset = Int++data Neighbourhood x y = Neighbourhood {+     _dataAtNode :: y+   , _neighbours :: UArr.Vector WebNodeIdOffset+   , _localScalarProduct :: Metric x+   , _webBoundaryAtNode :: Maybe (Needle' x)+   }+  deriving (Generic, Functor, Foldable, Traversable)+makeLenses ''Neighbourhood++deriving instance ( WithField ℝ PseudoAffine x+                  , SimpleSpace (Needle x), Show (Needle' x), Show y )+             => Show (Neighbourhood x y)++data WebLocally x y = LocalWebInfo {+      _thisNodeCoord :: x+    , _thisNodeData :: y+    , _thisNodeId :: WebNodeId+    , _nodeNeighbours :: [(WebNodeId, (Needle x, WebLocally x y))]+    , _nodeLocalScalarProduct :: Metric x+    , _webBoundingPlane :: Maybe (Needle' x)+    } deriving (Generic)+makeLenses ''WebLocally++data NeighbourhoodVector x = NeighbourhoodVector+          { _nvectId :: Int+          , _theNVect :: Needle x+          , _nvectNormal :: Needle' x+          , _nvectLength :: Scalar (Needle x)+          , _otherNeighboursOverlap :: Scalar (Needle x)+          }+makeLenses ''NeighbourhoodVector++data PropagationInconsistency x υ = PropagationInconsistency {+      _inconsistentPropagatedData :: [(x,υ)]+    , _inconsistentAPrioriData :: υ }+  | PropagationInconsistencies [PropagationInconsistency x υ]+ deriving (Show)+makeLenses ''PropagationInconsistency++instance Monoid (PropagationInconsistency x υ) where+  mempty = PropagationInconsistencies []+  mappend p q = mconcat [p,q]+  mconcat = PropagationInconsistencies++instance (NFData x, NFData (Metric x), NFData (Needle' x), NFData y)+           => NFData (Neighbourhood x y)++-- | A 'PointsWeb' is almost, but not quite a mesh. It is a stongly connected†+--   directed graph, backed by a tree for fast nearest-neighbour lookup of points.+-- +--   †In general, there can be disconnected components, but every connected+--   component is strongly connected.+newtype PointsWeb :: * -> * -> * where+   PointsWeb :: {+       webNodeRsc :: x`Shaded`Neighbourhood x y+     } -> PointsWeb x y+  deriving (Generic, Functor, Foldable, Traversable)++instance (NFData x, NFData (Metric x), NFData (Needle' x), NFData y) => NFData (PointsWeb x y)++instance CCt.Foldable (PointsWeb x) (->) (->) where+  ffoldl = uncurry . Hask.foldl' . curry+  foldMap = Hask.foldMap+++data WebChunk x y = WebChunk {+     _thisChunk :: PointsWeb x y+   , _layersAroundChunk :: [(x`Shaded`Neighbourhood x y, WebNodeId)]+   }++makeLenses ''WebChunk++data NodeInWeb x y = NodeInWeb {+     _thisNodeOnly :: (x, Neighbourhood x y)+   , _layersAroundNode :: [(x`Shaded`Neighbourhood x y, WebNodeId)]+   }+makeLenses ''NodeInWeb++data PathStep x y = PathStep {+     _pathStepStart :: WebLocally x y+   , _pathStepEnd :: WebLocally x y+   }+makeLenses ''PathStep+++type MetricChoice x = Shade x -> Metric x+++traverseInnermostChunks :: ∀ f x y z . Applicative f+          => (WebChunk x y -> f (PointsWeb x z)) -> PointsWeb x y -> f (PointsWeb x z)+traverseInnermostChunks f = go []+ where go :: [(x`Shaded`Neighbourhood x y, WebNodeId)] -> PointsWeb x y -> f (PointsWeb x z)+       go outlayers (w@(PointsWeb (PlainLeaves _)))+         = f (WebChunk w outlayers) +       go outlayers (PointsWeb w) = PointsWeb <$> traverseTrunkBranchChoices travel w+        where travel :: (Int, (Shaded x (Neighbourhood x y)))+                 -> Shaded x (Neighbourhood x y)+                 -> f (Shaded x (Neighbourhood x z))+              travel (i₀, br) obrs+                  = webNodeRsc <$> go ((obrs,i₀) : outlayers) (PointsWeb br)++traverseNodesInEnvi :: ∀ f x y z . Applicative f+           => (NodeInWeb x y -> f (Neighbourhood x z))+             -> PointsWeb x y -> f (PointsWeb x z)+traverseNodesInEnvi f = traverseInnermostChunks fc+ where fc :: WebChunk x y -> f (PointsWeb x z)+       fc (WebChunk (PointsWeb (PlainLeaves lvs)) outlayers)+            = PointsWeb . PlainLeaves <$> Hask.traverse fn (ixedFoci lvs)+        where fn ((i, (x, ngbh)), nearbyLeaves)+               = (x,) <$> f (NodeInWeb (x,ngbh)+                                     $ (PlainLeaves nearbyLeaves, i) : outlayers)++fmapNodesInEnvi :: (NodeInWeb x y -> Neighbourhood x z) -> PointsWeb x y -> PointsWeb x z+fmapNodesInEnvi f = runIdentity . traverseNodesInEnvi (Identity . f)+++ixedFoci :: [a] -> [((Int, a), [a])]+ixedFoci = go 0+ where go _ [] = []+       go i (x:xs) = ((i,x), xs) : map (second (x:)) (go (i+1) xs)+ ++indexWeb :: PointsWeb x y -> WebNodeId -> Maybe (x,y)+indexWeb (PointsWeb rsc) i = case indexShadeTree rsc i of+       Right (_, (x, Neighbourhood y _ _ _)) -> Just (x, y)+       _ -> Nothing++unsafeIndexWebData :: PointsWeb x y -> WebNodeId -> y+unsafeIndexWebData web i = case indexWeb web i of+              Just (x,y) -> y+++jumpNodeOffset :: WebNodeIdOffset -> NodeInWeb x y -> NodeInWeb x y+jumpNodeOffset 0 node = node+jumpNodeOffset δi (NodeInWeb x environment)+   = case zoomoutWebChunk δie $ WebChunk (PointsWeb $ PlainLeaves [x]) environment of+       (WebChunk bigChunk envi', δi')+           -> case pickNodeInWeb bigChunk δi' of+              NodeInWeb x' envi'' -> NodeInWeb x' $ envi'' ++ envi'+ where δie | δi < 0     = δi+           | otherwise  = δi - 1++webAroundChunk :: WebChunk x y -> PointsWeb x y+webAroundChunk (WebChunk chunk []) = chunk+webAroundChunk (WebChunk (PointsWeb (PlainLeaves lvs))+                         ((PlainLeaves lvsAround, i) : envi))+   = webAroundChunk $ WebChunk (PointsWeb . PlainLeaves $ lvsBefore++lvs++lvsAfter) envi+ where (lvsBefore, lvsAfter) = splitAt i lvsAround+webAroundChunk (WebChunk (PointsWeb chunk)+                         ((OverlappingBranches nw ew (DBranch dir+                            (Hourglass (PlainLeaves[]) d) :| brs), 0) : envi))+   = webAroundChunk $ WebChunk (PointsWeb $ OverlappingBranches (nw+nLeaves chunk) ew+                                          (DBranch dir (Hourglass chunk d) :| brs))+                               envi+webAroundChunk (WebChunk (PointsWeb chunk)+                         ((OverlappingBranches nw ew (DBranch dir+                            (Hourglass u (PlainLeaves[])) :| brs), i) : envi))+ | i==nLeaves u+   = webAroundChunk $ WebChunk (PointsWeb $ OverlappingBranches (nw+nLeaves chunk) ew+                                          (DBranch dir (Hourglass u chunk) :| brs))+                               envi+webAroundChunk (WebChunk chunk+                         (( OverlappingBranches nw ew (br₀@(DBranch _ (Hourglass u d))+                                                          :|br₁:brs)+                          , i) : envi))+  = case webAroundChunk (WebChunk chunk [(OverlappingBranches nw ew (br₁:|brs), i')])+      of PointsWeb (OverlappingBranches nw' ew' (br₁':|brs'))+           -> webAroundChunk $ WebChunk+                    (PointsWeb $ OverlappingBranches nw' ew' (br₀:|br₁':brs'))+                    envi+ where i' = i - nLeaves u - nLeaves d+webAroundChunk (WebChunk _ ((OverlappingBranches nw ew branches, i):_))+    = error $ "Environment with branch sizes "++show (fmap nLeaves . Hask.toList<$>(Hask.toList branches))+                ++" does not have a gap at #"++show i+webAroundChunk (WebChunk _ ((PlainLeaves _, _):_))+    = error "Encountered non-PlainLeaves chunk in a PlainLeaves environment."+++zoomoutWebChunk :: WebNodeIdOffset -> WebChunk x y -> (WebChunk x y, WebNodeId)+zoomoutWebChunk δi (WebChunk chunk ((outlayer, olp) : outlayers))+  | δi < -olp || δi >= nLeaves outlayer - olp+      = zoomoutWebChunk δiOut $ WebChunk widerChunk outlayers+  | otherwise  = (WebChunk widerChunk outlayers, δiIn)+ where δiOut | δi < 0     = δi + olp+             | otherwise  = δi + olp - nLeaves outlayer+       δiIn | δi < 0     = δi + olp+            | otherwise  = δi + olp + nLeaves (webNodeRsc chunk)+       widerChunk = webAroundChunk $ WebChunk chunk [(outlayer,olp)]+zoomoutWebChunk δi (WebChunk _ e)+    = error $ "Can't zoom out δ"++show δi+       ++" from a chunk with "++show (length e)++" environment layers."++pickNodeInWeb :: PointsWeb x y -> WebNodeId -> NodeInWeb x y+pickNodeInWeb = go [] id+ where go _ _ (PointsWeb w) i+        | i<0 || i>=n  = error+           $ "Trying to pick node #"++show i++" in web with "++show n++" nodes."+        where n = nLeaves w+       go lyrsAcc lMod (PointsWeb (PlainLeaves lvs)) i+        | (preds, node:succs)<-splitAt i lvs+                   = NodeInWeb node $ lMod (PlainLeaves $ preds++succs, i) : lyrsAcc+       go lyrsAcc lMod+            (PointsWeb (OverlappingBranches nw ew (DBranch dir (Hourglass u d):|brs))) i+        | i < nu     = go (lMod (OverlappingBranches (nw-nu) ew+                                      (DBranch dir (Hourglass gap d):|brs), 0) : lyrsAcc)+                          id (PointsWeb u) i+        | i < nu+nd  = go (lMod (OverlappingBranches (nw-nd) ew+                                      (DBranch dir (Hourglass u gap):|brs), nu) : lyrsAcc)+                          id (PointsWeb d) (i-nu)+        | (b:rs)<-brs  = go+                          lyrsAcc+                          (lMod . \(OverlappingBranches nwe ewe brse, ne)+                                   -> ( OverlappingBranches (nwe+nu+nd) ewe+                                         $ NE.cons (DBranch dir (Hourglass u d)) brse+                                      , ne+nu+nd ) )+                          (PointsWeb $ OverlappingBranches (nw-nu-nd) ew (b:|rs))+                          (i-nu-nd)+        where gap = PlainLeaves []+              [nu,nd] = nLeaves<$>[u,d]+++webLocalInfo :: ∀ x y . WithField ℝ Manifold x+            => PointsWeb x y -> PointsWeb x (WebLocally x y)+webLocalInfo = fmapNodesInEnvi linkln+ where linkln :: NodeInWeb x y -> Neighbourhood x (WebLocally x y)+       linkln node@(NodeInWeb (x, locloc@(Neighbourhood y ngbs metric nBoundary)) envis)+           = locloc & dataAtNode .~ LocalWebInfo {+                  _thisNodeCoord = x+                , _thisNodeData = y+                , _thisNodeId = i+                , _nodeNeighbours = [ (i + δi, (δx, ngb))+                                    | δi <- UArr.toList ngbs+                                    , let ngbNode@(NodeInWeb (xn, _) _)+                                              = jumpNodeOffset δi node+                                          Just δx = xn .-~. x+                                          Neighbourhood ngb _ _ _ = linkln ngbNode ]+                , _nodeLocalScalarProduct = metric+                , _webBoundingPlane = nBoundary+                }+        where i = foldr ((+) . snd) 0 envis+++instance Functor (WebLocally x) where+  fmap f (LocalWebInfo co dt id ng sp bn)+       = LocalWebInfo co (f dt) id (map (second . second $ fmap f) ng) sp bn+instance WithField ℝ Manifold x => Comonad (WebLocally x) where+  extract = _thisNodeData+  extend f this@(LocalWebInfo co _ id ng sp bn)+      = LocalWebInfo co (f this) id (map (second . second $ extend f) ng) sp bn+  duplicate this@(LocalWebInfo co _ id ng sp bn)+      = LocalWebInfo co this id (map (second $ second duplicate) ng) sp bn++-- ^ 'fmap' from the co-Kleisli category of 'WebLocally'.+localFmapWeb :: WithField ℝ Manifold x+                => (WebLocally x y -> z) -> PointsWeb x y -> PointsWeb x z+localFmapWeb f = webLocalInfo >>> fmap f+++tweakWebGeometry :: (WithField ℝ Manifold x, SimpleSpace (Needle x))+         => MetricChoice x -> (WebLocally x y -> [WebNodeId])+                        -> PointsWeb x y -> PointsWeb x y+tweakWebGeometry metricf reknit = webLocalInfo >>> fmapNodesInEnvi`id`+         \(NodeInWeb (x₀, (Neighbourhood info _ lm bound)) _)+             -> let lm' = metricf . Shade (inInterior x₀) $ dualNorm lm+                in Neighbourhood (info^.thisNodeData)+                            (UArr.fromList . map (subtract $ info^.thisNodeId)+                                     $ reknit info)+                            lm' bound+++bidirectionaliseWebLinks :: ∀ x y . PointsWeb x y -> PointsWeb x y+bidirectionaliseWebLinks web@(PointsWeb wnrsrc) = fmapNodesInEnvi bdse web+ where bdse :: NodeInWeb x y -> Neighbourhood x y+       bdse (NodeInWeb (x, Neighbourhood y outgn lm bound) envis)+                = Neighbourhood y (UArr.fromList . fastNub $ incmn ++ UArr.toList outgn)+                      lm bound+        where i = foldr ((+) . snd) 0 envis+              incmn = case i `Map.lookup` incoming of+                Just o -> subtract i<$>o+                Nothing -> []+       incoming = Map.fromListWith (++) $ Hask.foldl'+                   (\(i,acc) (Neighbourhood _ outgn _ _)+                        -> (i+1, acc . (((,[i]).(+i)<$>UArr.toList outgn)++)) )+                     (0,id) wnrsrc `snd` []++++pumpHalfspace :: ∀ v . (SimpleSpace v, Scalar v ~ ℝ)+     => Norm v+     -> v                    -- ^ A vector @v@ for which we want @dv<.>^v ≥ 0@.+     -> (DualVector v, [v])  -- ^ A plane @dv₀@ and some vectors @ws@ with @dv₀<.>^w ≥ 0@,+                             --   which should also fulfill @dv<.>^w ≥ 0@.+     -> Maybe (DualVector v) -- ^ The plane @dv@ fulfilling these properties, if possible.+pumpHalfspace rieM v (prevPlane, ws) = case dualSpaceWitness :: DualSpaceWitness v of+ DualSpaceWitness -> +  let    ϑs = fmap (\u -> let x = prevPlane<.>^u+                              y = thisPlane<.>^u+                          in atan2 (x-y) (x+y)) $ v:ws+          -- ϑ = 0 means we are mid-between the planes, ϑ > π/2 means we are past+          -- `thisPlane`, ϑ < -π/2 we are past `prevPlane`. In other words, positive ϑ+          -- mean we should mix in more of `prevPlane`, negative more of `thisPlane`.+         [ϑmin, ϑmax] = [minimum, maximum] <*> [ϑs]+         δϑ = ϑmax - ϑmin+         vNudged = v ^+^ sumV (zipWith (^*) ws smallPseudorandSeq)+                    -- Introduce a tiny contribution from the other vectors to avoid+                    -- a degenerate 1D-situation in which @thisPlane ∝ prevPlane@.+         dv = rieM<$|vNudged+         thisPlane = dv ^/ (dv<.>^vNudged)+         cas ϑ = cos $ ϑ - pi/4+  in if δϑ <= pi then Just $ let ϑbest = ϑmin + δϑ/2+                             in prevPlane^*cas ϑbest ^+^ thisPlane^*cas (-ϑbest)+                 else Nothing++smallPseudorandSeq :: [ℝ]+smallPseudorandSeq = (*2^^(-45)) . fromIntegral <$> lcg 293633+ where lcg x = x : lcg ((a*x)`mod`m)+       m = 2^31 - 1+       a = 963345    :: Int  -- revised Park-Miller++data LinkingBadness r = LinkingBadness+    { gatherDirectionsBadness :: !r -- ^ Prefer picking neighbours at right angles+                                    --   to the currently-explored-boundary. This+                                    --   is needed while we still have to link to+                                    --   points in different spatial directions.+    , closeSystemBadness :: !r      -- ^ Prefer points directly opposed to the+                                    --   current boundary. This is useful when the+                                    --   system of directions is already complete+                                    --   and we want a nicely symmetric “ball” of+                                    --   neighbours around each point.+    } deriving (Functor)++linkingUndesirability :: ℝ -- ^ Absolute-square distance (euclidean norm squared)+                      -> ℝ -- ^ Directional distance (distance from wall containing+                           --   all already known neighbours)+                      -> LinkingBadness ℝ+                           -- ^ “Badness” of this point as the next neighbour to link to.+                           --   In gatherDirections mode this is large if+                           --   the point is far away, but also if it is+                           --   right normal to the wall. The reason we punish this is that+                           --   adding two points directly opposed to each other would lead+                           --   to an ill-defined wall orientation, i.e. wrong normals+                           --   on the web boundary.+linkingUndesirability distSq wallDist+  | wallDist >= 0  = LinkingBadness+   { gatherDirectionsBadness = distSq^2 / max 0 (distSq-wallDist^2)+   , closeSystemBadness = distSq - wallDist^2/2+   }+  | otherwise     = LinkingBadness (1/0) (1/0) +++bestNeighbours :: ∀ i v . (SimpleSpace v, Scalar v ~ ℝ)+                => Norm v -> [(i,v)] -> ([i], Maybe (DualVector v))+bestNeighbours lm' = first (map fst) . bestNeighbours' lm'++bestNeighbours' :: ∀ i v . (SimpleSpace v, Scalar v ~ ℝ)+                => Norm v -> [(i,v)] -> ([(i,v)], Maybe (DualVector v))+bestNeighbours' lm' = extractSmallestOn (\(_,v) -> Just $ lm'|$|v) >>>+    \(Just ((c₀i,c₀δx), candidates)) -> case dualSpaceWitness :: DualSpaceWitness v of+     DualSpaceWitness ->+       let wall₀ = w₀ ^/ (lm|$|w₀) -- sqrt (w₀<.>^c₀δx)+            where w₀ = lm'<$|c₀δx+       in first ((c₀i,c₀δx):)+              $ gatherGoodNeighbours lm' lm wall₀ [c₀δx] [] candidates+ where lm = dualNorm lm' :: Variance v++gatherGoodNeighbours :: ∀ i v . (SimpleSpace v, Scalar v ~ ℝ)+            => Norm v -> Variance v+               -> DualVector v -> [v] -> [(i,v)]+                    -> [(i, v)] -> ([(i,v)], Maybe (DualVector v))+gatherGoodNeighbours lm' lm wall prev preserved cs+ | dimension == 1  = case extractSmallestOn+                       (\(_,δx) -> do+                          let wallDist = - wall<.>^δx+                          guard (wallDist > 0)+                          return wallDist+                       ) cs of+     Just (r, _) -> ([r], Nothing)+     Nothing -> ([], Just wall)+ where dimension = subbasisDimension (entireBasis :: SubBasis v)+gatherGoodNeighbours lm' lm wall prev preserved cs+  = case dualSpaceWitness :: DualSpaceWitness v of+    DualSpaceWitness ->+     case extractSmallestOn+            (\(_,δx) -> do+                let wallDist = - wall<.>^δx+                    dx = lm' <$| δx+                    distSq = dx<.>^δx+                    βmin = minimum [ 1 - (dx<.>^δxo) / sqrt (distSq*distSqo)+                                            -- β behaves basically like ϑ², where ϑ is+                                            -- the angle between two neighbour candidates.+                                   | (δxo, distSqo) <- prevWMag ]+                guard (wallDist >= 0 && βmin > 1e-3)+                return $ gatherDirectionsBadness+                           (linkingUndesirability distSq wallDist) / βmin )+            cs of+         Just ((i,δx), cs')+           | Just wall' <- pumpHalfspace lm' δx (wall,prev)+                          -> first ((i,δx):)+                       $ gatherGoodNeighbours lm' lm (wall'^/(lm|$|wall'))+                               (δx:prev) [] (preserved++cs')+           | (_:_)<-cs'  -> gatherGoodNeighbours lm' lm wall+                               prev ((i,δx):preserved) cs'+         _ -> let closeSys ((i,δx):_)+                    | Nothing <- pumpHalfspace lm' δx (wall,prev)+                        = ([(i,δx)], Nothing)+                  closeSys (_:cs'') = closeSys cs''+                  closeSys []+                   | null closureCandidates  = ([], Just wall)+                   | otherwise  = ([], Nothing)+                  closureCandidates = +                   [ ((i,δx), badness)+                   | (i,δx) <- preserved++cs+                   , let wallDist = - wall<.>^δx+                         distSq = normSq lm' δx+                   , wallDist > 0+                   , wallDist^2 > 1e-3 * distSq+                   , let badness = linkingUndesirability distSq wallDist ]+              in closeSys . map fst $+                   sortBy (comparing $ closeSystemBadness . snd) closureCandidates+ where prevWMag = map (id &&& normSq lm') prev+++extractSmallestOn :: Ord b => (a -> Maybe b) -> [a] -> Maybe (a, [a])+extractSmallestOn f = extract . map (id &&& f)+ where extract [] = Nothing+       extract ((x, Just o):cs) = Just $ go (o,x) cs+       extract ((x, Nothing):cs) = second (x:) <$> extract cs+       go (_,refx) [] = (refx, [])+       go (ref,refx) ((x, Just o):cs)+        | o < ref   = second (refx:) $ go (o,x) cs+       go ref ((x, _):cs) = second (x:) $ go ref cs++type WNIPath = [WebNodeId]+type NodeSet = ℤSet.IntSet++traversePathInIWeb :: ∀ φ x y . (WithField ℝ Manifold x, Monad φ, HasCallStack)+     => [WebNodeId] -> (PathStep x y -> φ y)+              -> PointsWeb x (WebLocally x y) -> φ (PointsWeb x (WebLocally x y))+traversePathInIWeb path f = go path+ where go [] web = pure web+       go [_] web = pure web+       go (i₀:i₁:is) web = do+                   y' <- f $ PathStep p₀ p₁+                   let Right (Identity web')+                         = treeLeaf i₁ (dataAtNode . thisNodeData . const $ pure y')+                              $ webNodeRsc web+                   go (i₁:is) $ PointsWeb web'+        where Just (_, p₀) = indexWeb web i₀+              Just (_, p₁) = indexWeb web i₁++traversePathsTowards :: ∀ f φ x y+         . (WithField ℝ Manifold x, Monad φ, Monad f, HasCallStack)+     => WebNodeId  -- ^ The node towards which the paths should converge.+       -> (PathStep x y -> φ y)+                   -- ^ The action which to traverse along each path.+       -> (∀ υ . WebLocally x y -> φ υ -> f υ)+                   -- ^ Initialisation/evaluation for each path-traversal.+       -> PointsWeb x y -> f (PointsWeb x y)+traversePathsTowards target pathTravF routeInitF web+  | Nothing <- sn  = error $ "Node "++show target++" not in web."+  | otherwise      = fmap (fmap _thisNodeData) . (`execStateT`envied) . forM_ paths+                       $ \path@(p₀:_) -> StateT $ \webState -> do+                 let Just (_, node₀) = indexWeb webState p₀+                 ((),) <$> routeInitF node₀ (traversePathInIWeb path pathTravF webState)+ where envied = webLocalInfo $ bidirectionaliseWebLinks web+       sn@(Just (targetPos,targetNode)) = indexWeb envied target+       paths = go 0 ℤSet.empty False [[target]] []+        where go :: Int -> NodeSet -> Bool -> [WNIPath] -> [WNIPath] -> [WNIPath]+              go targetDist visitedNodes boundaryCreepingInhibitor workers finishedThreads+               = case continue (round (sqrt $ fromIntegral targetDist :: Double))+                        visitedNodes boundaryCreepingInhibitor workers of+                  (_, [], _, newFinished) -> newFinished ++ finishedThreads+                  (visited', continuation, alternatives, newFinished)+                       -> let newThreads = filter (`ℤSet.notMember`visited')+                                                  (ℤSet.toList alternatives)+                          in go (targetDist+1)+                                (ℤSet.union visited' alternatives)+                                True+                                (continuation ++ fmap pure newThreads)+                                (newFinished ++ finishedThreads)+              continue :: Int -> NodeSet -> Bool -> [WNIPath]+                             -> (NodeSet, [WNIPath], NodeSet, [WNIPath])+              continue _ visitedNodes _ [] = (visitedNodes, [], ℤSet.empty, [])+              continue dfsDepth visitedNodes boundaryCreepingInhibitor ((cursor:nds):paths)+                  = case fst <$> sortBy (comparing snd)+                          [ goDfs dfsDepth (ℤSet.insert i visitedNodes) [] [i]+                          | i <- candidates ] of+                       ((preferred, (visited', pAlts)):alts)+                         | Nothing <- guard boundaryCreepingInhibitor+                                       >> cursorNode ^. webBoundingPlane+                          -> case continue dfsDepth visited'+                                         boundaryCreepingInhibitor paths of+                               (visited'', contin'', alts', newFin)+                                 -> ( visited''+                                    , (preferred++cursor:nds):contin''+                                    , ℤSet.union (ℤSet.fromList+                                                   $ pAlts ++ (last . fst <$> alts))+                                                 alts'+                                    , newFin )+                       alts -> case continue dfsDepth visitedNodes+                                           boundaryCreepingInhibitor paths of+                               (visited'', contin'', alts', newFin)+                                 -> ( visited''+                                    , contin''+                                    , ℤSet.union (ℤSet.fromList+                                                   $ last . fst <$> alts)+                                                 alts'+                                    , if null nds then newFin+                                                  else (cursor:nds):newFin )+               where Just (cursorPos,cursorNode) = indexWeb envied cursor+                     tgtOpp = cursorNode^.nodeLocalScalarProduct+                                  <$| targetPos .-~! cursorPos+                     candidates = [ ngb+                                  | (ngb, (_, _)) <- cursorNode^.nodeNeighbours+                                  , ngb`ℤSet.notMember`visitedNodes ]+                     goDfs :: Int -> NodeSet -> [WebNodeId] -> WNIPath+                                -> ((WNIPath, (NodeSet, [WebNodeId])), ℝ)+                     goDfs d2go visited' oldAlts (p:old)+                            = case sortBy (comparing snd) candidates of+                                ((preferred,oppositionQ):alts)+                                  -> let visited'' = ℤSet.insert preferred visited'+                                         alts' = filter (/=preferred) oldAlts+                                                    ++ (fst<$>alts)+                                     in if d2go>1+                                         then goDfs (d2go-1) visited'' alts'+                                                   (preferred:p:old)+                                         else ( (preferred:p:old, (visited'', alts'))+                                              , oppositionQ )+                                [] -> let δn = explorePos .-~! cursorPos+                                          oppositionQ = tgtOpp<.>^δn+                                      in ((p:old, (visited', oldAlts)), oppositionQ)+                      where Just (explorePos,exploreNode) = indexWeb envied p+                            candidates = [ (ngb, tgtOpp<.>^δn)+                                         | (ngb, (_, ngbN)) <- exploreNode^.nodeNeighbours+                                         , ngb`ℤSet.notMember`visited'+                                         , isNothing (ngbN^.webBoundingPlane)+                                         , let δn = ngbN^.thisNodeCoord .-~! cursorPos ]

binary file changed (52444 → 52978 bytes)

binary file changed (247358 → 277311 bytes)

binary file changed (absent → 279648 bytes)

binary file changed (absent → 271740 bytes)

binary file changed (328843 → 626463 bytes)

binary file changed (213891 → 310066 bytes)

@@ -1,5 +1,5 @@ Name:                manifolds-Version:             0.4.1.0+Version:             0.4.4.0 Category:            Math Synopsis:            Coordinate-free hypersurfaces Description:         Manifolds, a generalisation of the notion of &#x201c;smooth curves&#x201d; or surfaces,@@ -28,7 +28,7 @@ Homepage:            https://github.com/leftaroundabout/manifolds Maintainer:          (@) sagemueller $ geo.uni-koeln.de Build-Type:          Simple-Cabal-Version:       >=1.10+Cabal-Version:       >=1.18 Extra-Doc-Files:     images/examples/*.png,                      images/examples/ShadeCombinations/2Dconvolution-skewed.png                      images/examples/TreesAndWebs/*.png@@ -40,14 +40,14 @@  Library   Build-Depends:     base>=4.5 && < 6-                     , manifolds-core == 0.4.1.0+                     , manifolds-core == 0.4.4.0                      , transformers                      , vector-space>=0.8                      , free-vector-spaces>=0.1.1                      , linear                      , MemoTrie                      , vector-                     , linearmap-category >= 0.3.2 && < 0.4+                     , linearmap-category >= 0.3.4 && < 0.4                      , containers                      , comonad                      , free@@ -56,8 +56,11 @@                      , number-show >= 0.1 && < 0.2                      , tagged                      , deepseq+                     , placeholders                      , lens+                     , call-stack                      , constrained-categories >= 0.2.3 && < 0.3.1+                     , pragmatic-show   other-extensions:  FlexibleInstances                      , TypeFamilies                      , FlexibleContexts@@ -74,7 +77,9 @@                      Data.Manifold.TreeCover                      Data.Manifold.Shade                      Data.Manifold.Web+                     Data.Manifold.Web.Internal                      Data.Manifold.DifferentialEquation+                     Data.Manifold.Function.LocalModel                      Data.SimplicialComplex                      Data.Function.Differentiable                      Data.Function.Affine@@ -95,3 +100,26 @@                    Util.Associate                    Util.LtdShow   default-language: Haskell2010++test-suite test+  default-language:+    Haskell2010+  type:+    exitcode-stdio-1.0+  hs-source-dirs:+    test/tasty+  main-is:+    test.hs+  build-depends:+      base >= 4 && < 5+    , tasty >= 0.7+    , tasty-hunit+    , tasty-quickcheck+    , manifolds+    , pragmatic-show+    , containers+    , vector-space+    , constrained-categories+    , linearmap-category+    , lens+

@@ -0,0 +1,476 @@+-- |+-- Module      : Main+-- Copyright   : (c) Justus Sagemüller 2017+-- License     : GPL v3+-- +-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de+-- Stability   : experimental+-- Portability : portable+-- ++{-# LANGUAGE OverloadedLists, TypeFamilies, FlexibleContexts, UndecidableInstances #-}+{-# LANGUAGE TypeOperators #-}++module Main where++import Data.Manifold.Types+import Data.Manifold.PseudoAffine+import Data.Manifold.TreeCover+import Data.Manifold.Web+import Data.Manifold.Web.Internal+import Data.Manifold.Function.LocalModel+import Data.VectorSpace+import Math.LinearMap.Category+import Prelude hiding (id, fst, snd)+import Control.Category.Constrained (id)+import Control.Arrow.Constrained (fst,snd)++import Test.Tasty+import Test.Tasty.HUnit+import qualified Test.Tasty.QuickCheck as QC+import Test.Tasty.QuickCheck ((==>))++import Data.Foldable (toList)+import Data.List (nub)+import qualified Data.Graph as Graph+import qualified Data.Set as Set+import Control.Arrow+import Control.Lens++import qualified Text.Show.Pragmatic as SP+++main = defaultMain tests++tests :: TestTree+tests = testGroup "Tests"+ [ 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”"+        $ 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],[11,5],[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] ]+    , 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+  (≈) :: e -> e -> Bool+instance (SimpleSpace v, Needle v~v, Interior v~v, Floating (Scalar v))+             => AEq (Shade' v) where+  Shade' c₀ σ₀ ≈ Shade' c₁ σ₁+    = (σ₀|$|δ) < ε && (σ₀|$|δ) < ε+     && all (is1 . (σ₀|$|)) (normSpanningSystem' σ₁)+     && all (is1 . (σ₁|$|)) (normSpanningSystem' σ₀)+   where δ = c₁ ^-^ c₀+         ε = 1e-2+         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₁ σ₁+    = (dualNorm σ₀|$|δ) < ε && (dualNorm σ₀|$|δ) < ε+     && all (is1 . (dualNorm σ₀|$|)) (normSpanningSystem σ₁)+     && all (is1 . (dualNorm σ₁|$|)) (normSpanningSystem σ₀)+   where δ = c₁ ^-^ c₀+         ε = 1e-2+         is1 x = abs (x-1) < ε+instance AEq a => AEq (Maybe a) where+  Just x ≈ Just y = x ≈ y+  Nothing ≈ Nothing = True+  _ ≈ _ = 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₁+                                        +infix 1 @?≈       +(@?≈) :: (AEq e, Show e) => e -> e -> Assertion+a@?≈b+ | a≈b        = return ()+ | otherwise  = assertFailure $ "Expected "++show b++", but got "++show a++