manifolds 0.4.0.0 → 0.4.1.0

raw patch · 15 files changed

+1969/−1056 lines, 15 filesdep +freedep +number-showdep ~linearmap-categorydep ~manifolds-corebinary-addedPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: free, number-show

Dependency ranges changed: linearmap-category, manifolds-core

API changes (from Hackage documentation)

Files

• Data/Function/Affine.hs +53/−0
• Data/Function/Differentiable.hs +131/−64
• Data/Function/Differentiable/Data.hs +1/−1
• Data/Manifold/Atlas.hs +29/−0
• Data/Manifold/Cone.hs +1/−1
• Data/Manifold/DifferentialEquation.hs +86/−29
• Data/Manifold/Function/Quadratic.hs +124/−0
• Data/Manifold/PseudoAffine.hs +24/−49
• Data/Manifold/Riemannian.hs +74/−13
• Data/Manifold/Shade.hs +1218/−0
• Data/Manifold/TreeCover.hs +11/−798
• Data/Manifold/Types.hs +17/−8
• Data/Manifold/Web.hs +192/−89
• images/examples/TreesAndWebs/2D-cartesian_onion.png binary
• manifolds.cabal +8/−4

@@ -33,6 +33,8 @@               Affine(..)             , evalAffine             , fromOffsetSlope+            -- * Misc+            , lensEmbedding, correspondingDirections             ) where      @@ -46,9 +48,11 @@ import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine import Data.Manifold.Atlas+import Data.Embedding  import qualified Prelude import qualified Control.Applicative as Hask+import Data.Functor (($>))  import Control.Category.Constrained.Prelude hiding ((^)) import Control.Category.Constrained.Reified@@ -58,8 +62,10 @@  import Math.LinearMap.Category +import Control.Lens  + data Affine s d c where     Affine :: (ChartIndex d :->: (c, LinearMap s (Needle d) (Needle c)))                -> Affine s d c@@ -198,3 +204,50 @@        -> \y0 ðx'y -> Affine . trie $ chartReferencePoint                     >>> \x₀ -> let δy = ðx'y $ x₀                                in (y0.+~^δy, ðx'y)+++instance EnhancedCat (Embedding (Affine s)) (Embedding (LinearMap s)) where+  arr (Embedding e p) = Embedding (arr e) (arr p)+++lensEmbedding :: ∀ k s x c .+                 ( Num' s+                 , LinearSpace x, LinearSpace c, Object k x, Object k c+                 , Scalar x ~ s, Scalar c ~ s+                 , EnhancedCat k (LinearMap s) )+                  => Lens' x c -> Embedding k c x+lensEmbedding l = Embedding (arr $ (arr $ LinearFunction (\c -> zeroV & l .~ c)+                                     :: LinearMap s c x) )+                            (arr $ (arr $ LinearFunction (^.l)+                                     :: LinearMap s x c) )+++correspondingDirections :: ∀ s x c t+                        . ( WithField s AffineManifold c+                          , WithField s AffineManifold x+                          , SemiInner (Needle c), SemiInner (Needle x)+                          , RealFrac' s+                          , Traversable t )+              => (Interior c, Interior x)+                  -> t (Needle c, Needle x) -> Maybe (Embedding (Affine s) c x)+correspondingDirections (c₀, x₀) dirMap+   = freeEmbeddings $> Embedding (Affine . trie $ c2x boundarylessWitness)+                                 (Affine . trie $ x2c boundarylessWitness)+ where freeEmbeddings = fzip ( embedFreeSubspace $ fst<$>dirMap+                             , embedFreeSubspace $ snd<$>dirMap )+       c2t :: Lens' (Needle c) (t s)+       c2t = case freeEmbeddings of Just (Lens ct, _) -> ct+       x2t :: Lens' (Needle x) (t s)+       x2t = case freeEmbeddings of Just (_, Lens xt) -> xt+       c2x :: BoundarylessWitness c -> ChartIndex c+                            -> (x, LinearMap s (Needle c) (Needle x))+       c2x BoundarylessWitness ιc+              = ( x₀ .+~^ (zeroV & x2t .~ δc^.c2t)+                , arr . LinearFunction $ \dc -> zeroV & x2t .~ dc^.c2t )+        where Just δc = chartReferencePoint ιc .-~. c₀+       x2c :: BoundarylessWitness x -> ChartIndex x+                            -> (c, LinearMap s (Needle x) (Needle c))+       x2c BoundarylessWitness ιx+              = ( c₀ .+~^ (zeroV & c2t .~ δx^.x2t)+                , arr . LinearFunction $ \dx -> zeroV & c2t .~ dx^.x2t )+        where Just δx = chartReferencePoint ιx .-~. x₀

@@ -77,6 +77,9 @@   +type RealDimension s+       = ( RealFloat' s, SimpleSpace s, Show s, Atlas s, HasTrie (ChartIndex s)+         , s ~ Needle s, s ~ Interior s, s ~ Scalar s, s ~ DualVector s )   discretisePathIn :: (WithField ℝ Manifold y, SimpleSpace (Needle y))@@ -265,10 +268,12 @@   -unsafe_dev_ε_δ :: RealDimension a+unsafe_dev_ε_δ :: ∀ a . RealDimension a                 => String -> (a -> a) -> LinDevPropag a a-unsafe_dev_ε_δ errHint f d-            = let ε'² = normSq d 1+unsafe_dev_ε_δ = case ( linearManifoldWitness :: LinearManifoldWitness a+                      , closedScalarWitness :: ClosedScalarWitness a ) of+ (LinearManifoldWitness _, ClosedScalarWitness) -> \errHint f d+           -> let ε'² = normSq d 1               in if ε'²>0                   then let δ = f . sqrt $ recip ε'²                        in if δ > 0@@ -278,9 +283,12 @@                                     ++show(sqrt $ recip ε'²)                                     ++ " gives non-positive δ="++show δ++"."                   else mempty-dev_ε_δ :: RealDimension a+dev_ε_δ :: ∀ a . RealDimension a          => (a -> a) -> Metric a -> Maybe (Metric a)-dev_ε_δ f d = let ε'² = normSq d 1+dev_ε_δ = case ( linearManifoldWitness :: LinearManifoldWitness a+                      , closedScalarWitness :: ClosedScalarWitness a ) of+ (LinearManifoldWitness _, ClosedScalarWitness) -> \f d+           -> let ε'² = normSq d 1               in if ε'²>0                   then let δ = f . sqrt $ recip ε'²                        in if δ > 0@@ -288,8 +296,11 @@                            else empty                   else pure mempty -as_devεδ :: RealDimension a => LinDevPropag a a -> a -> a-as_devεδ ldp ε | ε>0+as_devεδ :: ∀ a . RealDimension a => LinDevPropag a a -> a -> a+as_devεδ = asdevεδ linearManifoldWitness closedScalarWitness where+ asdevεδ :: LinearManifoldWitness a -> ClosedScalarWitness a -> LinDevPropag a a -> a -> a+ asdevεδ (LinearManifoldWitness _) ClosedScalarWitness+         ldp ε | ε>0                , δ'² <- normSq (ldp $ spanNorm [recip ε]) 1                , δ'² > 0                     = sqrt $ recip δ'²@@ -466,40 +477,61 @@   -instance (LocallyScalable s v, LinearManifold v, LocallyScalable s a, RealFloat' s)+instance ∀ v s a . (LinearSpace v, Scalar v ~ s, LocallyScalable s a, RealFloat' s)     => AdditiveGroup (DfblFuncValue s a v) where-  zeroV = point zeroV-  GenericAgent (AffinDiffable ef f) ^+^ GenericAgent (AffinDiffable eg g)-         = GenericAgent $ AffinDiffable (ef<>eg) (f^+^g)-  α^+^β = dfblFnValsCombine (\a b -> (a^+^b, arr addV, const mempty)) α β-  negateV (GenericAgent (AffinDiffable ef f))-         = GenericAgent $ AffinDiffable ef (negateV f)-  negateV α = dfblFnValsFunc (\a -> (negateV a, negateV id, const mempty)) α+  zeroV = case ( linearManifoldWitness :: LinearManifoldWitness v+               , dualSpaceWitness :: DualSpaceWitness v ) of+     (LinearManifoldWitness _, DualSpaceWitness) -> point zeroV+  (^+^) = case ( linearManifoldWitness :: LinearManifoldWitness v+               , dualSpaceWitness :: DualSpaceWitness v ) of+     (LinearManifoldWitness _, DualSpaceWitness)+         -> curry $ \case+        (GenericAgent (AffinDiffable ef f), GenericAgent (AffinDiffable eg g))+              -> GenericAgent $ AffinDiffable (ef<>eg) (f^+^g)+        (α,β) -> dfblFnValsCombine (\a b -> (a^+^b, arr addV, const mempty)) α β+  negateV = case ( linearManifoldWitness :: LinearManifoldWitness v+                 , dualSpaceWitness :: DualSpaceWitness v ) of+      (LinearManifoldWitness _, DualSpaceWitness) -> \case+         (GenericAgent (AffinDiffable ef f))+           -> GenericAgent $ AffinDiffable ef (negateV f)+         α -> dfblFnValsFunc (\a -> (negateV a, negateV id, const mempty)) α   -instance (RealDimension n, LocallyScalable n a)+instance ∀ n a . (RealDimension n, LocallyScalable n a)             => Num (DfblFuncValue n a n) where-  fromInteger i = point $ fromInteger i-  (+) = (^+^)-  (*) = dfblFnValsCombine $+  fromInteger = case ( linearManifoldWitness :: LinearManifoldWitness n+                     , closedScalarWitness :: ClosedScalarWitness n ) of+      (LinearManifoldWitness _, ClosedScalarWitness) -> point . fromInteger+  (+) = case closedScalarWitness :: ClosedScalarWitness n of+      ClosedScalarWitness -> (^+^)+  (*) = case ( linearManifoldWitness :: LinearManifoldWitness n+             , closedScalarWitness :: ClosedScalarWitness n ) of+      (LinearManifoldWitness _, ClosedScalarWitness) -> dfblFnValsCombine $           \a b -> ( a*b                   , arr $ addV <<< (scale $ a)***(scale $ b)-                  , unsafe_dev_ε_δ(show a++"*"++show b) sqrt+                  , unsafe_dev_ε_δ(show a++"*"++show b) (sqrt :: n->n)                        >>> \d¹₂ -> (d¹₂,d¹₂)                            -- ε δa δb = (a+δa)·(b+δb) - (a·b + (a·δa + b·δb))                             --         = δa·δb                            --   so choose δa = δb = √ε                   )-  negate = negateV-  abs = dfblFnValsFunc dfblAbs-   where dfblAbs a-          | a>0        = (a, id, unsafe_dev_ε_δ("abs "++show a) $ \ε -> a + ε/2) -          | a<0        = (-a, negateV id, unsafe_dev_ε_δ("abs "++show a) $ \ε -> ε/2 - a)-          | otherwise  = (0, zeroV, scaleNorm (sqrt 0.5))-  signum = dfblFnValsFunc dfblSgn-   where dfblSgn a-          | a>0        = (1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ const a)-          | a<0        = (-1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ \_ -> -a)-          | otherwise  = (0, zeroV, const $ spanNorm [1])+  negate = case closedScalarWitness :: ClosedScalarWitness n of+     ClosedScalarWitness -> negateV+  abs = mkabs linearManifoldWitness closedScalarWitness+   where mkabs :: LinearManifoldWitness n -> ClosedScalarWitness n+                     -> DfblFuncValue n a n -> DfblFuncValue n a n+         mkabs (LinearManifoldWitness _) ClosedScalarWitness = dfblFnValsFunc dfblAbs+          where dfblAbs a+                 | a>0        = (a, id, unsafe_dev_ε_δ("abs "++show a) $ \ε -> a + ε/2) +                 | a<0        = (-a, negateV id, unsafe_dev_ε_δ("abs "++show a) $ \ε -> ε/2 - a)+                 | otherwise  = (0, zeroV, scaleNorm (sqrt 0.5))+  signum = mksgn linearManifoldWitness closedScalarWitness+   where mksgn :: LinearManifoldWitness n -> ClosedScalarWitness n+                     -> DfblFuncValue n a n -> DfblFuncValue n a n+         mksgn (LinearManifoldWitness _) ClosedScalarWitness = dfblFnValsFunc dfblSgn+          where dfblSgn a+                 | a>0        = (1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ const a)+                 | a<0        = (-1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ \_ -> -a)+                 | otherwise  = (0, zeroV, const $ spanNorm [1])   @@ -516,10 +548,13 @@   -- | Important special operator needed to compute intersection of 'Region's.-minDblfuncs :: (LocallyScalable s m, RealDimension s)+minDblfuncs :: ∀ s m . (LocallyScalable s m, RealDimension s)      => Differentiable s m s -> Differentiable s m s -> Differentiable s m s-minDblfuncs (Differentiable f) (Differentiable g) = Differentiable h- where h x+minDblfuncs (Differentiable f) (Differentiable g)+             = Differentiable $ h linearManifoldWitness closedScalarWitness+ where h :: LinearManifoldWitness s -> ClosedScalarWitness s+             -> m -> (s, Needle m+>Needle s, LinDevPropag m s)+       h (LinearManifoldWitness _) ClosedScalarWitness x          | fx < gx   = ( fx, jf                        , \d -> devf d <> devg d                                <> transformNorm δj@@ -542,9 +577,11 @@   -genericisePreRegion :: (RealDimension s, LocallyScalable s m)+genericisePreRegion :: ∀ s m . (RealDimension s, LocallyScalable s m)                           => PreRegion s m -> PreRegion s m-genericisePreRegion GlobalRegion = PreRegion $ const 1+genericisePreRegion GlobalRegion = case ( linearManifoldWitness :: LinearManifoldWitness s+                                        , closedScalarWitness :: ClosedScalarWitness s ) of+    (LinearManifoldWitness _, ClosedScalarWitness) -> PreRegion $ const 1 genericisePreRegion (RealSubray PositiveHalfSphere xl) = preRegionToInfFrom' xl genericisePreRegion (RealSubray NegativeHalfSphere xr) = preRegionFromMinInfTo' xr genericisePreRegion r = r@@ -571,13 +608,19 @@ regionProd (Region a₀ ra) (Region b₀ rb) = Region (a₀,b₀) (preRegionProd ra rb)  -- | Cartesian product of two pre-regions.-preRegionProd :: (RealDimension s, LocallyScalable s a, LocallyScalable s b)+preRegionProd :: ∀ s a b . (RealDimension s, LocallyScalable s a, LocallyScalable s b)                   => PreRegion s a -> PreRegion s b -> PreRegion s (a,b)-preRegionProd GlobalRegion GlobalRegion = GlobalRegion-preRegionProd GlobalRegion (PreRegion rb) = PreRegion $ rb . snd-preRegionProd (PreRegion ra) GlobalRegion = PreRegion $ ra . fst-preRegionProd (PreRegion ra) (PreRegion rb) = PreRegion $ minDblfuncs (ra.fst) (rb.snd)-preRegionProd ra rb = preRegionProd (genericisePreRegion ra) (genericisePreRegion rb)+preRegionProd = prp linearManifoldWitness closedScalarWitness+ where prp :: LinearManifoldWitness s -> ClosedScalarWitness s+                 -> PreRegion s a -> PreRegion s b -> PreRegion s (a,b)+       prp _ _ GlobalRegion GlobalRegion = GlobalRegion+       prp (LinearManifoldWitness _) ClosedScalarWitness GlobalRegion (PreRegion rb)+                    = PreRegion $ rb . snd+       prp (LinearManifoldWitness _) ClosedScalarWitness (PreRegion ra) GlobalRegion+                    = PreRegion $ ra . fst+       prp (LinearManifoldWitness _) ClosedScalarWitness (PreRegion ra) (PreRegion rb)+                    = PreRegion $ minDblfuncs (ra.fst) (rb.snd)+       prp _ _ ra rb = preRegionProd (genericisePreRegion ra) (genericisePreRegion rb)   positivePreRegion, negativePreRegion :: (RealDimension s) => PreRegion s s@@ -585,9 +628,13 @@ negativePreRegion = RealSubray NegativeHalfSphere 0  -positivePreRegion', negativePreRegion' :: (RealDimension s) => PreRegion s s-positivePreRegion' = PreRegion $ Differentiable prr- where prr x = ( 1 - 1/xp1+positivePreRegion', negativePreRegion' :: ∀ s . (RealDimension s) => PreRegion s s+positivePreRegion' = PreRegion . Differentiable+                       $ prr linearManifoldWitness closedScalarWitness+ where prr :: LinearManifoldWitness s -> ClosedScalarWitness s+           -> s -> (s, Needle s+>Needle s, LinDevPropag s s)+       prr (LinearManifoldWitness _) ClosedScalarWitness+           x = ( 1 - 1/xp1                , (1/xp1²) *^ id                , unsafe_dev_ε_δ("positivePreRegion@"++show x) δ )                  -- ε = (1 − 1/(1+x)) + (-δ · 1/(x+1)²) − (1 − 1/(1+x−δ))@@ -624,26 +671,42 @@                   | otherwise  = ε * x / ((1+ε)/x + ε)               xp1 = (x+1)               xp1² = xp1 ^ 2-negativePreRegion' = PreRegion $ ppr . ngt- where PreRegion ppr = positivePreRegion'-       ngt = actuallyLinearEndo $ negateV id+negativePreRegion' = npr (linearManifoldWitness :: LinearManifoldWitness s)+                         (closedScalarWitness :: ClosedScalarWitness s)+ where npr (LinearManifoldWitness BoundarylessWitness)+           (ClosedScalarWitness :: ClosedScalarWitness s)+                  = PreRegion $ ppr . ngt+        where PreRegion ppr = positivePreRegion'+              ngt = actuallyLinearEndo $ negateV id  preRegionToInfFrom, preRegionFromMinInfTo :: RealDimension s => s -> PreRegion s s preRegionToInfFrom = RealSubray PositiveHalfSphere preRegionFromMinInfTo = RealSubray NegativeHalfSphere -preRegionToInfFrom', preRegionFromMinInfTo' :: RealDimension s => s -> PreRegion s s-preRegionToInfFrom' xs = PreRegion $ ppr . trl- where PreRegion ppr = positivePreRegion'-       trl = actuallyAffineEndo (-xs) id-preRegionFromMinInfTo' xe = PreRegion $ ppr . flp- where PreRegion ppr = positivePreRegion'-       flp = actuallyAffineEndo xe (negateV id)+preRegionToInfFrom', preRegionFromMinInfTo' :: ∀ s . RealDimension s => s -> PreRegion s s+preRegionToInfFrom' = prif (linearManifoldWitness :: LinearManifoldWitness s)+                           (closedScalarWitness :: ClosedScalarWitness s)+ where prif (LinearManifoldWitness BoundarylessWitness)+            (ClosedScalarWitness :: ClosedScalarWitness s)+            xs = PreRegion $ ppr . trl+        where PreRegion ppr = positivePreRegion'+              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'+              flp = actuallyAffineEndo xe (negateV id) -intervalPreRegion :: RealDimension s => (s,s) -> PreRegion s s-intervalPreRegion (lb,rb) = PreRegion $ Differentiable prr+intervalPreRegion :: ∀ s . RealDimension s => (s,s) -> PreRegion s s+intervalPreRegion (lb,rb) = PreRegion . Differentiable+                             $ prr linearManifoldWitness closedScalarWitness  where m = lb + radius; radius = (rb - lb)/2-       prr x = ( 1 - ((x-m)/radius)^2+       prr :: LinearManifoldWitness s -> ClosedScalarWitness s+                -> s -> (s, Needle s+>Needle s, LinDevPropag s s)+       prr (LinearManifoldWitness _) ClosedScalarWitness+           x = ( 1 - ((x-m)/radius)^2                , (2*(m-x)/radius^2) *^ id                , unsafe_dev_ε_δ("intervalPreRegion@"++show x) $ (*radius) . sqrt ) @@ -1049,8 +1112,8 @@   -instance ( RealDimension n, WithField n Manifold a-         , LocallyScalable n a, SimpleSpace (Needle a) )+instance ∀ n a . ( RealDimension n, WithField n Manifold a+                 , LocallyScalable n a, SimpleSpace (Needle a) )             => Floating (RWDfblFuncValue n a n) where   pi = point pi   @@ -1058,7 +1121,8 @@     $ \x -> let ex = exp x             in if ex*2 == ex  -- numerical trouble...                 then if x<0 then ( 0, zeroV, unsafe_dev_ε_δ("exp "++show x) $ \ε -> log ε - x )-                            else ( ex, ex*^id, unsafe_dev_ε_δ("exp "++show x) $ \_ -> 1e-300 )+                            else ( ex, ex*^id+                                 , unsafe_dev_ε_δ("exp "++show x) $ \_ -> 1e-300 :: n )                 else ( ex, ex *^ id, unsafe_dev_ε_δ("exp "++show x)                           $ \ε -> case acosh(ε/(2*ex) + 1) of                                     δ | δ==δ      -> δ@@ -1315,10 +1379,13 @@  -- | Like 'Data.VectorSpace.lerp', but gives a differentiable function --   instead of a Hask one.-lerp_diffable :: ( WithField s LinearManifold m, Atlas m-                 , HasTrie (ChartIndex m), RealDimension s )+lerp_diffable :: ∀ m s . ( LinearSpace m, Scalar m ~ s, Atlas m+                         , HasTrie (ChartIndex m), RealDimension s )       => m -> m -> Differentiable s s m-lerp_diffable a b = actuallyAffine a . arr $ flipBilin scale $ b.-.a+lerp_diffable = case ( linearManifoldWitness :: LinearManifoldWitness m+                     , dualSpaceWitness :: DualSpaceWitness m ) of+     (LinearManifoldWitness BoundarylessWitness, DualSpaceWitness)+         -> \a b -> actuallyAffine a . arr $ flipBilin scale $ b.-.a   

@@ -91,7 +91,7 @@ --   includes that point) to define a connected subset of a manifold. data PreRegion s m where   GlobalRegion :: PreRegion s m-  RealSubray :: RealDimension s => S⁰ -> s -> PreRegion s s+  RealSubray :: Num' s => S⁰ -> s -> PreRegion s s   PreRegion :: (Differentiable s m s) -- A function that is positive at reference point /p/,                                       -- decreases and crosses zero at the region's                                       -- boundaries. (If it goes positive again somewhere

@@ -9,6 +9,8 @@ --   {-# LANGUAGE TypeFamilies              #-}+{-# LANGUAGE ConstraintKinds           #-}+{-# LANGUAGE FlexibleContexts          #-} {-# LANGUAGE FlexibleInstances         #-} {-# LANGUAGE EmptyDataDecls, EmptyCase #-} {-# LANGUAGE CPP                       #-}@@ -25,9 +27,14 @@ import Data.Void  import Data.VectorSpace.Free+import Math.LinearMap.Category  import Control.Arrow +import Data.MemoTrie (HasTrie)++import qualified Linear.Affine as LinAff+ class Semimanifold m => Atlas m where   type ChartIndex m :: *   chartReferencePoint :: ChartIndex m -> m@@ -49,6 +56,8 @@ VectorSpaceAtlas(Num s, V2 s) VectorSpaceAtlas(Num s, V3 s) VectorSpaceAtlas(Num s, V4 s)+VectorSpaceAtlas((LinearSpace v, Scalar v ~ s, TensorSpace w, Scalar w ~ s), LinearMap s v w)+VectorSpaceAtlas((TensorSpace v, Scalar v ~ s, TensorSpace w, Scalar w ~ s), Tensor s v w)  instance (Atlas x, Atlas y) => Atlas (x,y) where   type ChartIndex (x,y) = (ChartIndex x, ChartIndex y)@@ -78,3 +87,23 @@   interiorChartReferencePoint _ NegativeHalfSphere = S² pi 0   lookupAtlas (S² ϑ _) | ϑ<pi/2     = PositiveHalfSphere                        | otherwise  = NegativeHalfSphere++instance (LinearSpace (a n), Needle (a n) ~ a n, Interior (a n) ~ a n)+              => Atlas (LinAff.Point a n) where+  type ChartIndex (LinAff.Point a n) = ()+  interiorChartReferencePoint _ () = LinAff.P zeroV+  lookupAtlas _ = ()++++-- | The 'AffineSpace' class plus manifold constraints.+type AffineManifold m = ( Atlas m, Manifold m, AffineSpace m+                        , Needle m ~ Diff m, HasTrie (ChartIndex m) )++-- | An euclidean space is a real affine space whose tangent space is a Hilbert space.+type EuclidSpace x = ( AffineManifold x, InnerSpace (Diff x)+                     , DualVector (Diff x) ~ Diff x, Floating (Scalar (Diff x)) )++euclideanMetric :: EuclidSpace x => proxy x -> Metric x+euclideanMetric _ = euclideanNorm+

@@ -143,7 +143,7 @@ projCD¹ToCℝay (CD¹ h m) = Cℝay (bijectIntvtoℝplus h) m  -stiefel1Project :: LinearManifold v =>+stiefel1Project :: LinearSpace v =>              DualVector v       -- ^ Must be nonzero.                  -> Stiefel1 v stiefel1Project = Stiefel1

@@ -33,10 +33,11 @@  module Data.Manifold.DifferentialEquation (             -- * Formulating simple differential eqns.-              DifferentialEqn+              DifferentialEqn, ODE+            , constLinearDEqn             , constLinearODE             , constLinearPDE-            , filterDEqnSolution_static, iterateFilterDEqn_static+            , iterateFilterDEqn_static             -- * Cost functions for error bounds             , maxDeviationsGoal             , uncertaintyGoal@@ -62,7 +63,10 @@ import Data.Function.Differentiable.Data import Data.Manifold.TreeCover import Data.Manifold.Web+import Data.Manifold.Atlas +import Data.Embedding+ import qualified Data.List as List  import qualified Prelude as Hask hiding(foldl, sum, sequence)@@ -79,39 +83,92 @@ import Data.Traversable.Constrained (Traversable, traverse)  +-- | An ordinary differential equation is one that does not need any a-priori+--   partial derivatives to compute the derivative for integration in some+--   propagation direction. Classically, ODEs are usually understood as+--   @DifferentialEquation ℝ ℝ⁰ y@, but actually @x@ can at least+--   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 -constLinearODE :: ∀ x y . ( WithField ℝ LinearManifold x, SimpleSpace x-                          , WithField ℝ LinearManifold y, SimpleSpace y )-              => ((x +> y) +> y) -> DifferentialEqn x y-constLinearODE = case ( dualSpaceWitness :: DualNeedleWitness x-                      , dualSpaceWitness :: DualNeedleWitness y ) of-   (DualSpaceWitness, DualSpaceWitness) -> \bwt' ->+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 = 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'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+                             :| [ projectShade+                                   (Embedding (arr bwt'inv <<< id&&&zeroV)+                                              (arr bwt'    >>> fst))+                                   shjApriori ]+                   , return $ projectShade+                                   (Embedding (arr bwt'inv <<< zeroV&&&id)+                                              (arr bwt'    >>> snd))+                                   shjApriori+                   )+         }++constLinearODE :: ∀ x y . ( SimpleSpace x, Scalar x ~ ℝ, SimpleSpace y, Scalar y ~ ℝ )+              => ((x +> y) +> y) -> ODE x y+constLinearODE = case ( linearManifoldWitness :: LinearManifoldWitness x+                      , dualSpaceWitness :: DualSpaceWitness x+                      , linearManifoldWitness :: LinearManifoldWitness y+                      , dualSpaceWitness :: DualSpaceWitness y ) of+   ( LinearManifoldWitness BoundarylessWitness, DualSpaceWitness+    ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ) -> \bwt' ->     let bwt'inv = (bwt'\$)     in \(Shade (_x,y) δxy) -> LocalDifferentialEqn             (let j = bwt'inv y                  δj = (bwt'>>>zeroV&&&id) `transformNorm` dualNorm δxy-             in return $ Shade' j δj )-            (\_ -> pure )+             in \_ -> return $ Shade' j δj )+            (\shy _ _ -> (pure shy, Just $ Shade' Origin mempty) ) -constLinearPDE :: ∀ x y y' .-                  ( WithField ℝ LinearManifold x, SimpleSpace x-                  , WithField ℝ LinearManifold y, SimpleSpace y, FiniteFreeSpace y-                  , WithField ℝ LinearManifold y', SimpleSpace y' )-              => ((x +> (y,y')) +> (y, y')) -> DifferentialEqn x (y,y')-constLinearPDE = undefined{-case ( dualSpaceWitness :: DualNeedleWitness x-                      , dualSpaceWitness :: DualNeedleWitness y-                      , dualSpaceWitness :: DualSpaceWitness y' ) of-   (DualSpaceWitness, DualSpaceWitness, DualSpaceWitness) -> \bwt' ->-    let bwt'inv = (bwt'\$)-    in  \(Shade (_x,(y,y')) δxy) (Shade' jApriori σjApriori)-                            -> let j = bwt'inv $ (zeroV,y')-                                   δj = (bwt'>>>zeroV&&&id)-                                         `transformNorm` dualNorm δxy-                                   (_,y'Apriori) = bwt' $ jApriori-                                   Norm δy' = (arr $ LinearFunction bwt'inv . (zeroV&&&id))-                                         `transformNorm` σjApriori-                             in (Shade' (y,y'Apriori) . Norm $ zeroV *** δy' , )-                              <$> mixShade's (Shade' jApriori σjApriori :| [Shade' j δj])-}+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,124 @@+-- |+-- Module      : Data.Manifold.Function.Quadratic+-- Copyright   : (c) Justus Sagemüller 2016+-- License     : GPL v3+-- +-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de+-- Stability   : experimental+-- Portability : portable+-- ++{-# LANGUAGE TypeFamilies             #-}+{-# LANGUAGE FlexibleContexts         #-}+{-# LANGUAGE LiberalTypeSynonyms      #-}+{-# LANGUAGE GADTs                    #-}+{-# LANGUAGE TypeOperators            #-}+{-# LANGUAGE UnicodeSyntax            #-}+{-# LANGUAGE ScopedTypeVariables      #-}+++module Data.Manifold.Function.Quadratic (+              Quadratic(..), evalQuadratic+            ) where++++import Data.Semigroup+import qualified Data.List.NonEmpty as NE++import Data.MemoTrie+import Data.VectorSpace+import Data.AffineSpace+import Data.Tagged+import Data.Manifold.PseudoAffine+import Data.Manifold.Atlas+import Data.Manifold.Riemannian+import Data.Function.Affine++import Prelude                       hiding (id, ($), fmap, fst)+import Control.Category.Constrained.Prelude (id, ($), fmap, fst)+import Control.Arrow.Constrained ((>>>), (&&&), (***), second)++import Math.LinearMap.Category++++data Quadratic s d c where+    Quadratic :: (ChartIndex d :->: ( c, ( LinearMap s (Needle d) (Needle c)+                                         , LinearMap s (SymmetricTensor s (Needle d))+                                                     (Needle c) ) )+                 )  -> Quadratic s d c++affineQuadratic :: (WithField s AffineManifold d, WithField s AffineManifold c)+        => Affine s d c -> Quadratic s d c+affineQuadratic (Affine f) = Quadratic . trie+                  $ untrie f >>> second (id &&& const zeroV)++instance ( Atlas x, HasTrie (ChartIndex x), LinearSpace (Needle x), Scalar (Needle x) ~ s+         , Manifold y, Scalar (Needle y) ~ s )+              => Semimanifold (Quadratic s x y) where+  type Needle (Quadratic s x y) = Quadratic s x (Needle y)+  toInterior = pure+  fromInterior = id+  (.+~^) = case ( semimanifoldWitness :: SemimanifoldWitness y+                , boundarylessWitness :: BoundarylessWitness y ) of+    (SemimanifoldWitness _, BoundarylessWitness) -> \(Quadratic f) (Quadratic g)+      -> Quadratic . trie $ \ix -> case (untrie f ix, untrie g ix) of+          ((fx₀,f'), (gx₀,g')) -> (fx₀.+~^gx₀, f'^+^g')+  translateP = Tagged (.+~^)+  semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness y of+    SemimanifoldWitness _ -> SemimanifoldWitness BoundarylessWitness+instance ( Atlas x, HasTrie (ChartIndex x), LinearSpace (Needle x), Scalar (Needle x) ~ s+         , Manifold y, Scalar (Needle y) ~ s )+              => PseudoAffine (Quadratic s x y) where+  (.-~!) = case ( semimanifoldWitness :: SemimanifoldWitness y+                , boundarylessWitness :: BoundarylessWitness y ) of+    (SemimanifoldWitness _, BoundarylessWitness) -> \(Quadratic f) (Quadratic g)+      -> Quadratic . trie $ \ix -> case (untrie f ix, untrie g ix) of+          ((fx₀,f'), (gx₀,g')) -> (fx₀.-~!gx₀, f'^-^g')+  pseudoAffineWitness = case semimanifoldWitness :: SemimanifoldWitness y of+    SemimanifoldWitness _ -> PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+instance ( Atlas x, HasTrie (ChartIndex x), LinearSpace (Needle x), Scalar (Needle x) ~ s+         , Manifold y, Scalar (Needle y) ~ s )+              => AffineSpace (Quadratic s x y) where+  type Diff (Quadratic s x y) = Quadratic s x (Needle y)+  (.+^) = (.+~^); (.-.) = (.-~!)+instance ( Atlas x, HasTrie (ChartIndex x), LinearSpace (Needle x), Scalar (Needle x) ~ s+         , LinearSpace y, Scalar y ~ s, Num' s )+            => AdditiveGroup (Quadratic s x y) where+  zeroV = case linearManifoldWitness :: LinearManifoldWitness y of+       LinearManifoldWitness _ -> Quadratic . trie $ const (zeroV, zeroV)+  (^+^) = case ( linearManifoldWitness :: LinearManifoldWitness y+               , dualSpaceWitness :: DualSpaceWitness y ) of+      (LinearManifoldWitness BoundarylessWitness, DualSpaceWitness) -> (.+~^)+  negateV = case linearManifoldWitness :: LinearManifoldWitness y of+       LinearManifoldWitness _ -> \(Quadratic f) -> Quadratic . trie $+             untrie f >>> negateV***negateV+instance ( Atlas x, HasTrie (ChartIndex x), LinearSpace (Needle x), Scalar (Needle x) ~ s+         , LinearSpace y, Scalar y ~ s, Num' s )+            => VectorSpace (Quadratic s x y) where+  type Scalar (Quadratic s x y) = s+  (*^) = case linearManifoldWitness :: LinearManifoldWitness y of+       LinearManifoldWitness _ -> \μ (Quadratic f) -> Quadratic . trie $+             untrie f >>> (μ*^)***(μ*^)++evalQuadratic :: ∀ s x y . ( Manifold x, Atlas x, HasTrie (ChartIndex x)+                           , Manifold y+                           , s ~ Scalar (Needle x), s ~ Scalar (Needle y) )+               => Quadratic s x y -> x+                    -> (y, ( LinearMap s (Needle x) (Needle y)+                           , LinearMap s (SymmetricTensor s (Needle x)) (Needle y) ))+evalQuadratic = ea (boundarylessWitness, boundarylessWitness)+ where ea :: (BoundarylessWitness x, BoundarylessWitness y)+             -> Quadratic s x y -> x -> (y, ( LinearMap s (Needle x) (Needle y)+                                            , LinearMap s (SymmetricTensor s (Needle x)) (Needle y) ))+       ea (BoundarylessWitness, BoundarylessWitness)+          (Quadratic f) x = ( fx₀.+~^(ðx'f₀ $ v).+~^(ð²x'f $ squareV v)+                            , ( ðx'f₀ ^+^ 2*^((currySymBilin $ ð²x'f) $ v)+                              , ð²x'f+                              ) )+        where Just v = x .-~. chartReferencePoint chIx+              chIx = lookupAtlas x+              (fx₀, (ðx'f₀, ð²x'f)) = untrie f chIx++

@@ -57,7 +57,7 @@             -- ** Needles             , Local(..)             -- ** Metrics-            , Metric, Metric', euclideanMetric+            , Metric, Metric'             , RieMetric, RieMetric'             -- ** Constraints             , SemimanifoldWitness(..)@@ -65,14 +65,10 @@             , BoundarylessWitness(..)             , boundarylessWitness             , DualNeedleWitness -            , RealDimension, AffineManifold-            , LinearManifold             , WithField-            , HilbertManifold-            , EuclidSpace             , LocallyScalable             -- ** Local functions-            , LocalLinear, LocalAffine+            , LocalLinear, LocalBilinear, LocalAffine             -- * Misc             , alerpB, palerp, palerpB, LocallyCoercible(..), CanonicalDiffeomorphism(..)             , ImpliesMetric(..), coerceMetric, coerceMetric'@@ -93,6 +89,7 @@ import qualified Linear.Affine as LinAff import Data.Embedding import Data.LinearMap+import Data.VectorSpace.Free import Math.LinearMap.Category import Data.AffineSpace import Data.Tagged@@ -108,6 +105,8 @@ import Control.Monad.Constrained import Data.Foldable.Constrained +import Control.Lens (Lens', lens, (^.), (&), (%~), (.~))+ import GHC.Exts (Constraint)  @@ -147,6 +146,8 @@                   => p (x,ξ) -> CanonicalDiffeomorphism (Interior x) (Interior ξ)   interiorLocalCoercion _ = CanonicalDiffeomorphism +type NumPrime n = (Num' n, Eq n)+ #define identityCoercion(c,t)                   \ instance (c) => LocallyCoercible (t) (t) where { \   locallyTrivialDiffeomorphism = id;              \@@ -154,16 +155,16 @@   coerceNeedle' _ = id;                             \   oppositeLocalCoercion = CanonicalDiffeomorphism;   \   interiorLocalCoercion _ = CanonicalDiffeomorphism }-identityCoercion(NumberManifold s, ZeroDim s)-identityCoercion(NumberManifold s, V0 s)+identityCoercion(NumPrime s, ZeroDim s)+identityCoercion(NumPrime s, V0 s) identityCoercion((), ℝ)-identityCoercion(NumberManifold s, V1 s)+identityCoercion(NumPrime s, V1 s) identityCoercion((), (ℝ,ℝ))-identityCoercion(NumberManifold s, V2 s)+identityCoercion(NumPrime s, V2 s) identityCoercion((), (ℝ,(ℝ,ℝ))) identityCoercion((), ((ℝ,ℝ),ℝ))-identityCoercion(NumberManifold s, V3 s)-identityCoercion(NumberManifold s, V4 s)+identityCoercion(NumPrime s, V3 s)+identityCoercion(NumPrime s, V4 s)   data CanonicalDiffeomorphism a b where@@ -180,15 +181,14 @@                            , Num' s )  type LocalLinear x y = LinearMap (Scalar (Needle x)) (Needle x) (Needle y)-type LocalAffine x y = (Needle y, LocalLinear x y)+type LocalBilinear x y = LinearMap (Scalar (Needle x))+                                   (SymmetricTensor (Scalar (Needle x)) (Needle x))+                                   (Needle y) --- | Basically just an “updated” version of the 'VectorSpace' class.---   Every vector space is a manifold, this constraint makes it explicit.-type LinearManifold x = ( AffineManifold x, Needle x ~ x, LSpace x ) -type LinearManifold' x = ( PseudoAffine x, AffineSpace x, Diff x ~ x-                         , Interior x ~ x, Needle x ~ x, LSpace x ) +type LocalAffine x y = (Needle y, LocalLinear x y)+ -- | Require some constraint on a manifold, and also fix the type of the manifold's --   underlying field. For example, @WithField ℝ 'HilbertManifold' v@ constrains --   @v@ to be a real (i.e., 'Double'-) Hilbert space.@@ -198,34 +198,7 @@ --   applied, for @type@ constraints this is by default not allowed). type WithField s c x = ( c x, s ~ Scalar (Needle x), s ~ Scalar (Needle' x) ) --- | The 'RealFloat' class plus manifold constraints.-type RealDimension r = ( PseudoAffine r, Interior r ~ r, Needle r ~ r, r ~ ℝ) --- | The 'AffineSpace' class plus manifold constraints.-type AffineManifold m = ( PseudoAffine m, Interior m ~ m, AffineSpace m-                        , Needle m ~ Diff m, LinearManifold' (Diff m) )---- | A Hilbert space is a /complete/ inner product space. Being a vector space, it is---   also a manifold.--- ---   (Stricly speaking, that doesn't have much to do with the completeness criterion;---   but since 'Manifold's are at the moment confined to finite dimension, they are in---   fact (trivially) complete.)-type HilbertManifold x = ( LinearManifold x, InnerSpace x-                         , Interior x ~ x, Needle x ~ x, DualVector x ~ x-                         , Floating (Scalar x) )---- | An euclidean space is a real affine space whose tangent space is a Hilbert space.-type EuclidSpace x = ( AffineManifold x, InnerSpace (Diff x)-                     , DualVector (Diff x) ~ Diff x, Floating (Scalar (Diff x)) )--type NumberManifold n = ( Num' n, Manifold n, Interior n ~ n, Needle n ~ n-                        , LSpace n, DualVector n ~ n, Scalar n ~ n )--euclideanMetric :: EuclidSpace x => proxy x -> Metric x-euclideanMetric _ = euclideanNorm-- -- | A co-needle can be understood as a “paper stack”, with which you can measure --   the length that a needle reaches in a given direction by counting the number --   of holes punched through them.@@ -289,11 +262,11 @@  deriveAffine(KnownNat n, FreeVect n ℝ) -instance (NumberManifold s) => LocallyCoercible (ZeroDim s) (V0 s) where+instance (NumPrime s) => LocallyCoercible (ZeroDim s) (V0 s) where   locallyTrivialDiffeomorphism Origin = V0   coerceNeedle _ = LinearFunction $ \Origin -> V0   coerceNeedle' _ = LinearFunction $ \Origin -> V0-instance (NumberManifold s) => LocallyCoercible (V0 s) (ZeroDim s) where+instance (NumPrime s) => LocallyCoercible (V0 s) (ZeroDim s) where   locallyTrivialDiffeomorphism V0 = Origin   coerceNeedle _ = LinearFunction $ \V0 -> Origin   coerceNeedle' _ = LinearFunction $ \V0 -> Origin@@ -376,13 +349,15 @@             -> CanonicalDiffeomorphism  -instance LinearManifold (a n) => Semimanifold (LinAff.Point a n) where+instance (LinearSpace (a n), Needle (a n) ~ a n, Interior (a n) ~ a n)+            => Semimanifold (LinAff.Point a n) where   type Needle (LinAff.Point a n) = a n   fromInterior = id   toInterior = pure   LinAff.P v .+~^ w = LinAff.P $ v ^+^ w   translateP = Tagged $ \(LinAff.P v) w -> LinAff.P $ v ^+^ w-instance LinearManifold (a n) => PseudoAffine (LinAff.Point a n) where+instance (LinearSpace (a n), Needle (a n) ~ a n, Interior (a n) ~ a n)+            => PseudoAffine (LinAff.Point a n) where   LinAff.P v .-~. LinAff.P w = return $ v ^-^ w  

@@ -41,12 +41,15 @@ {-# LANGUAGE LiberalTypeSynonyms        #-} {-# LANGUAGE CPP                        #-} {-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE DefaultSignatures          #-}   module Data.Manifold.Riemannian  where   import Data.Maybe+import Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NE import qualified Data.Vector as Arr  import Data.VectorSpace@@ -59,6 +62,7 @@ import Data.Manifold.Types.Primitive ((^), empty, embed, coEmbed) import Data.Manifold.Types.Stiefel import Data.Manifold.PseudoAffine+import Data.Manifold.Atlas (AffineManifold)      import Data.CoNat @@ -76,6 +80,9 @@ import Data.Foldable.Constrained  +data GeodesicWitness x where+  GeodesicWitness :: Geodesic (Interior x)+       => SemimanifoldWitness x -> GeodesicWitness x  class Semimanifold x => Geodesic x where   geodesicBetween ::@@ -85,7 +92,13 @@             --   If the two points are actually connected by a path...        -> Maybe (D¹ -> x) -- ^ ...then this is the interpolation function. Attention:                            --   the type will change to 'Differentiable' in the future.+  geodesicWitness :: GeodesicWitness x+  default geodesicWitness :: Geodesic (Interior x) => GeodesicWitness x+  geodesicWitness = GeodesicWitness semimanifoldWitness+  middleBetween :: x -> x -> Maybe x+  middleBetween p₀ p₁ = ($ D¹ 0) <$> geodesicBetween p₀ p₁ + interpolate :: (Geodesic x, IntervalLike i) => x -> x -> Maybe (i -> x) interpolate a b = (. toClosedInterval) <$> geodesicBetween a b @@ -94,36 +107,55 @@  #define deriveAffineGD(x)                                         \ instance Geodesic x where {                                        \-  geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹ \+  geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹; \+  middleBetween a b = return $ alerp a b (1/2) \  }  deriveAffineGD (ℝ) -instance Geodesic (ZeroDim ℝ) where+instance Geodesic (ZeroDim s) where   geodesicBetween Origin Origin = return $ \_ -> Origin+  middleBetween Origin Origin = return Origin -instance (Geodesic a, Geodesic b) => Geodesic (a,b) where+instance ∀ a b . (Geodesic a, Geodesic b) => Geodesic (a,b) where   geodesicBetween (a,b) (α,β) = liftA2 (&&&) (geodesicBetween a α) (geodesicBetween b β)+  geodesicWitness = case ( geodesicWitness :: GeodesicWitness a+                         , geodesicWitness :: GeodesicWitness b ) of+     (GeodesicWitness _, GeodesicWitness _) -> GeodesicWitness semimanifoldWitness+  middleBetween (a,b) (α,β) = fzip (middleBetween a α, middleBetween b β) -instance (Geodesic a, Geodesic b, Geodesic c) => Geodesic (a,b,c) where+instance ∀ a b c . (Geodesic a, Geodesic b, Geodesic c) => Geodesic (a,b,c) where   geodesicBetween (a,b,c) (α,β,γ)       = liftA3 (\ia ib ic t -> (ia t, ib t, ic t))            (geodesicBetween a α) (geodesicBetween b β) (geodesicBetween c γ)+  geodesicWitness = case ( geodesicWitness :: GeodesicWitness a+                         , geodesicWitness :: GeodesicWitness b+                         , geodesicWitness :: GeodesicWitness c ) of+     (GeodesicWitness _, GeodesicWitness _, GeodesicWitness _)+         -> GeodesicWitness semimanifoldWitness  -- instance (KnownNat n) => Geodesic (FreeVect n ℝ) where --   geodesicBetween (FreeVect v) (FreeVect w) --       = return $ \(D¹ t) -> let μv = (1-t)/2; μw = (t+1)/2 --                             in FreeVect $ Arr.zipWith (\vi wi -> μv*vi + μw*wi) v w -instance (Geodesic v, FiniteFreeSpace v, WithField ℝ HilbertManifold v)+instance ∀ v . ( Geodesic v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v)+               , LinearSpace v, Scalar v ~ ℝ, Geodesic (DualVector v)+               , InnerSpace (DualVector v) )              => Geodesic (Stiefel1 v) where-  geodesicBetween (Stiefel1 p') (Stiefel1 q')-      = (\f -> \(D¹ t) -> Stiefel1 . f . D¹ $ g * tan (ϑ*t))+  geodesicBetween = gb dualSpaceWitness+   where gb :: DualSpaceWitness v -> Stiefel1 v -> Stiefel1 v -> Maybe (D¹ -> Stiefel1 v)+         gb DualSpaceWitness (Stiefel1 p') (Stiefel1 q')+           = (\f -> \(D¹ t) -> Stiefel1 . f . D¹ $ g * tan (ϑ*t))             <$> geodesicBetween p q-   where p = normalized p'; q = normalized q'-         l = magnitude $ p^-^q-         ϑ = asin $ l/2-         g = sqrt $ 4/l^2 - 1+          where p = normalized p'; q = normalized q'+                l = magnitude $ p^-^q+                ϑ = asin $ l/2+                g = sqrt $ 4/l^2 - 1+  middleBetween = gb dualSpaceWitness+   where gb :: DualSpaceWitness v -> Stiefel1 v -> Stiefel1 v -> Maybe (Stiefel1 v)+         gb DualSpaceWitness  (Stiefel1 p) (Stiefel1 q)+             = Stiefel1 <$> middleBetween (normalized p) (normalized q)   instance Geodesic S⁰ where@@ -138,6 +170,10 @@                         <$> geodesicBetween (pi-φ) (-ϕ-pi)     | otherwise       = (>>> S¹ . \ψ -> signum ψ*pi - ψ)                         <$> geodesicBetween (-pi-φ) (pi-ϕ)+  middleBetween (S¹ φ) (S¹ ϕ)+    | abs (φ-ϕ) < pi  = S¹ <$> middleBetween φ ϕ+    | φ > 0           = S¹ <$> middleBetween (pi-φ) (-ϕ-pi)+    | otherwise       = S¹ <$> middleBetween (-pi-φ) (pi-ϕ)   -- instance Geodesic (Cℝay S⁰) where@@ -198,6 +234,15 @@ deriveAffineGD (ℝ³) deriveAffineGD (ℝ⁴) +instance (TensorSpace v, Scalar v ~ ℝ, TensorSpace w, Scalar w ~ ℝ)+             => Geodesic (Tensor ℝ v w) where+  geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹+instance (LinearSpace v, Scalar v ~ ℝ, TensorSpace w, Scalar w ~ ℝ)+             => Geodesic (LinearMap ℝ v w) where+  geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹+instance (TensorSpace v, Scalar v ~ ℝ, TensorSpace w, Scalar w ~ ℝ)+             => Geodesic (LinearFunction ℝ v w) where+  geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹   -- | One-dimensional manifolds, whose closure is homeomorpic to the unit interval.@@ -232,5 +277,21 @@   -middleBetween :: Geodesic m => m -> m -> Maybe m-middleBetween p₀ p₁ = ($ D¹ 0) <$> geodesicBetween p₀ p₁++pointsBarycenter :: Geodesic m => NonEmpty m -> Maybe m+pointsBarycenter (p:|[]) = Just p+pointsBarycenter ps = case ( pointsBarycenter (NE.fromList group₀)+                           , pointsBarycenter (NE.fromList group₁) ) of+            (Just bc₀, Just bc₁)+                | δn == 0      -> middleBetween bc₀ bc₁+                | otherwise    -> ($ D¹ (fromIntegral δn/fromIntegral ntot))+                                    <$> geodesicBetween bc₀ bc₁+            _                  -> Nothing+ where psl = Hask.toList ps+       (group₀, group₁) = splitAt nl psl+       ntot = length psl+       nl = ntot`quot`2+       δn = ntot  - nl*2+++type FlatSpace x = (AffineManifold x, Geodesic x, SimpleSpace x)

@@ -0,0 +1,1218 @@+-- |+-- Module      : Data.Manifold.Shade+-- Copyright   : (c) Justus Sagemüller 2016+-- License     : GPL v3+-- +-- Maintainer  : (@) jsagemue $ uni-koeln.de+-- Stability   : experimental+-- Portability : portable+-- +{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE UndecidableInstances       #-}+{-# LANGUAGE StandaloneDeriving         #-}+{-# LANGUAGE DeriveGeneric              #-}+{-# LANGUAGE DeriveFunctor              #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE GADTs                      #-}+{-# LANGUAGE RankNTypes                 #-}+{-# LANGUAGE ParallelListComp           #-}+{-# LANGUAGE UnicodeSyntax              #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ViewPatterns               #-}+{-# LANGUAGE LambdaCase                 #-}+{-# LANGUAGE TypeOperators              #-}+{-# LANGUAGE CPP                        #-}+{-# LANGUAGE TupleSections              #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE LiberalTypeSynonyms        #-}+{-# LANGUAGE DefaultSignatures          #-}+{-# LANGUAGE TemplateHaskell            #-}+++module Data.Manifold.Shade (+       -- * Shades +         Shade(..), pattern(:±), Shade'(..), (|±|), IsShade+       -- ** Lenses+       , shadeCtr, shadeExpanse, shadeNarrowness+       -- ** Construction+       , fullShade, fullShade', pointsShades, pointsShade's+       , pointsCovers, pointsCover's, coverAllAround+       -- ** Evaluation+       , occlusion, prettyShowsPrecShade', prettyShowShade', LtdErrorShow+       -- ** Misc+       , factoriseShade, orthoShades, (✠), intersectShade's, linIsoTransformShade+       , embedShade, projectShade+       , Refinable, subShade', refineShade', convolveShade', coerceShade+       , mixShade's, dualShade+       -- * Misc+       -- ** Shades+       , shadesMerge, pointsShades', pseudoECM, convolveMetric+       , WithAny(..), shadeWithAny, shadeWithoutAnything+       -- ** Local data fit models+       , estimateLocalJacobian, estimateLocalHessian, QuadraticModel(..)+       -- ** Differential equations+       , DifferentialEqn, LocalDifferentialEqn(..)+       , propagateDEqnSolution_loc, LocalDataPropPlan(..)+       -- ** Range interpolation+       , rangeOnGeodesic, rangeWithinVertices+    ) where+++import Data.List hiding (filter, all, elem, sum, foldr1)+import Data.Maybe+import Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NE+import Data.Semigroup+import Control.DeepSeq+import Data.MemoTrie++import Data.VectorSpace+import Data.AffineSpace+import Math.LinearMap.Category+import Data.Tagged+import Linear (_x,_y,_z,_w)++import Data.Manifold.Types+import Data.Manifold.Types.Primitive ((^))+import Data.Manifold.PseudoAffine+import Data.Manifold.Riemannian+import Data.Manifold.Atlas+import Data.Function.Affine+import Data.Manifold.Function.Quadratic++import Data.Embedding++import Control.Lens (Lens', (^.), view, _1, _2, mapping, (&))+import Control.Lens.TH++import qualified Prelude as Hask hiding(foldl, sum, sequence)+import qualified Control.Applicative as Hask+import qualified Data.Foldable       as Hask+import Data.Foldable (all, elem, toList, sum, foldr1)++import Control.Category.Constrained.Prelude hiding+     ((^), all, elem, sum, forM, Foldable(..), foldr1, Traversable, traverse)+import Control.Arrow.Constrained+import Control.Monad.Constrained hiding (forM)++import GHC.Generics (Generic)++import Text.Show.Number+++-- | A 'Shade' is a very crude description of a region within a manifold. It+--   can be interpreted as either an ellipsoid shape, or as the Gaussian peak+--   of a normal distribution (use <http://hackage.haskell.org/package/manifold-random>+--   for actually sampling from that distribution).+-- +--   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) }+deriving instance (Show (Interior x), Show (Metric' x), WithField ℝ PseudoAffine x)+                => Show (Shade x)++-- | A &#x201c;co-shade&#x201d; can describe ellipsoid regions as well, but unlike+--   'Shade' it can be unlimited / infinitely wide in some directions.+--   It does OTOH need to have nonzero thickness, which 'Shade' needs not.+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''.+  shadeCtr :: Lens' (shade x) (Interior x)+--  -- | Convert between 'Shade' and 'Shade' (which must be neither singular nor infinite).+--  unsafeDualShade :: WithField ℝ Manifold x => shade x -> shade* x+  -- | Check the statistical likelihood-density of a point being within a shade.+  --   This is taken as a normal distribution.+  occlusion :: ( PseudoAffine x, SimpleSpace (Needle x)+               , s ~ (Scalar (Needle x)), RealFloat' s )+                => shade x -> x -> s+  factoriseShade :: ( PseudoAffine x, SimpleSpace (Needle x)+                    , 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+  -- | ASCII version of '✠'.+  orthoShades :: ( PseudoAffine x, SimpleSpace (Needle x)+           , PseudoAffine y, SimpleSpace (Needle y)+           , Scalar (Needle x) ~ Scalar (Needle y) )+                => shade x -> shade y -> shade (x,y)+  linIsoTransformShade :: ( SimpleSpace x, SimpleSpace y, Scalar x ~ Scalar y+                          , Num' (Scalar x) )+                          => (x+>y) -> shade x -> shade y+  -- | 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) )+                        => Embedding (Affine s) (Interior x) (Interior y)+                              -> shade y -> shade x+  -- | Include a shade in a higher-dimensional space. Notice that this behaves+  --   fundamentally different for 'Shade' and 'Shade''. For 'Shade', it gives+  --   a “flat image” of the region, whereas for 'Shade'' it gives an “extrusion+  --   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) )+                        => 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)+                  => (x+>y) -> Shade x -> Shade y+linearProjectShade = case ( linearManifoldWitness :: LinearManifoldWitness x+                          , linearManifoldWitness :: LinearManifoldWitness y+                          , dualSpaceWitness :: DualSpaceWitness x+                          , dualSpaceWitness :: DualSpaceWitness y ) of+   ( LinearManifoldWitness BoundarylessWitness+    ,LinearManifoldWitness BoundarylessWitness+    ,DualSpaceWitness, DualSpaceWitness )+       -> \f (Shade x ex) -> Shade (f $ x) (transformVariance f ex)+++infixl 5 ✠+-- | Combine two shades on independent subspaces to a shade with the same+--   properties on the subspaces (see 'factoriseShade') and no covariance.+(✠) :: ( IsShade shade, PseudoAffine x, SimpleSpace (Needle x)+       , PseudoAffine y, SimpleSpace (Needle y)+       , Scalar (Needle x) ~ Scalar (Needle y) )+                => shade x -> shade y -> shade (x,y)+(✠) = orthoShades++instance IsShade Shade where+  shadeCtr f (Shade c e) = fmap (`Shade`e) $ f c+  occlusion = occ pseudoAffineWitness dualSpaceWitness+   where occ :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x)+                        , Scalar (Needle x) ~ s, RealFloat' s )+                    => PseudoAffineWitness x -> DualNeedleWitness x -> Shade x -> x -> s+         occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness (Shade p₀ δ)+                 = \p -> case toInterior p >>= (.-~.p₀) of+           (Just vd) | mSq <- normSq δinv vd+                     , mSq == mSq  -- avoid NaN+                     -> exp (negate mSq)+           _         -> zeroV+          where δinv = dualNorm δ+  factoriseShade = fs dualSpaceWitness dualSpaceWitness+   where fs :: ∀ x y . ( PseudoAffine x, SimpleSpace (Needle x)+                       , PseudoAffine y, SimpleSpace (Needle y)+                       , Scalar (Needle x) ~ Scalar (Needle y) )+               => DualNeedleWitness x -> DualNeedleWitness y+                       -> Shade (x,y) -> (Shade x, Shade y)+         fs DualSpaceWitness DualSpaceWitness (Shade (x₀,y₀) δxy)+                   = (Shade x₀ δx, Shade y₀ δy)+          where (δx,δy) = summandSpaceNorms δxy+  orthoShades = fs dualSpaceWitness dualSpaceWitness+   where fs :: ∀ x y . ( PseudoAffine x, SimpleSpace (Needle x)+                       , PseudoAffine y, SimpleSpace (Needle y)+                       , Scalar (Needle x) ~ Scalar (Needle y) )+               => DualNeedleWitness x -> DualNeedleWitness y+                       -> Shade x -> Shade y ->  Shade (x,y)+         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)+                => DualNeedleWitness x -> DualNeedleWitness y -> Shade x -> Shade y+         cS DualSpaceWitness DualSpaceWitness+                    = \(Shade x δxym) -> Shade (internCoerce x) (tN δxym)+          where tN = case oppositeLocalCoercion :: CanonicalDiffeomorphism y x of+                      CanonicalDiffeomorphism ->+                       transformNorm . arr $ coerceNeedle' ([]::[(y,x)])+                internCoerce = case interiorLocalCoercion ([]::[(x,y)]) of+                      CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism+  linIsoTransformShade = lits linearManifoldWitness linearManifoldWitness+                              dualSpaceWitness dualSpaceWitness+   where lits :: ∀ x y . ( LinearSpace x, LinearSpace y+                         , Scalar x ~ Scalar y, Num' (Scalar x) )+               => LinearManifoldWitness x -> LinearManifoldWitness y+                   -> DualSpaceWitness x -> DualSpaceWitness y+                       -> (x+>y) -> Shade x -> Shade y+         lits (LinearManifoldWitness BoundarylessWitness)+              (LinearManifoldWitness BoundarylessWitness)+              DualSpaceWitness DualSpaceWitness+              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) )+                        => (SemimanifoldWitness x, SemimanifoldWitness y)+               -> Embedding (Affine s) (Interior x) (Interior y)+                              -> Shade x -> Shade y+         ps' (SemimanifoldWitness _, SemimanifoldWitness _)+              (Embedding q _) (Shade x e) = Shade y (transformVariance j e)+          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) )+                        => (SemimanifoldWitness x, SemimanifoldWitness y)+               -> Embedding (Affine s) (Interior x) (Interior y)+                              -> Shade y -> Shade x+         ps' (SemimanifoldWitness _, SemimanifoldWitness _)+              (Embedding _ q) (Shade x e) = Shade y (transformVariance j e)+          where y = q $ x+                (_,j) = evalAffine q x+++dualShade :: ∀ x . (PseudoAffine x, SimpleSpace (Needle x))+                => Shade x -> Shade' x+dualShade = case dualSpaceWitness :: DualSpaceWitness (Needle x) of+    DualSpaceWitness -> \(Shade c e) -> Shade' c $ dualNorm e++dualShade' :: ∀ x . (PseudoAffine x, SimpleSpace (Needle x))+                => Shade' x -> Shade x+dualShade' = case dualSpaceWitness :: DualSpaceWitness (Needle x) of+    DualSpaceWitness -> \(Shade' c e) -> Shade c $ dualNorm' e++instance ImpliesMetric Shade where+  type MetricRequirement Shade x = (Manifold x, SimpleSpace (Needle x))+  inferMetric' (Shade _ e) = e+  inferMetric = im dualSpaceWitness+   where im :: (Manifold x, SimpleSpace (Needle x))+                   => DualNeedleWitness x -> Shade x -> Metric x+         im DualSpaceWitness (Shade _ e) = dualNorm e++instance ImpliesMetric Shade' where+  type MetricRequirement Shade' x = (Manifold x, SimpleSpace (Needle x))+  inferMetric (Shade' _ e) = e+  inferMetric' (Shade' _ e) = dualNorm e++shadeExpanse :: Lens' (Shade x) (Metric' x)+shadeExpanse f (Shade c e) = fmap (Shade c) $ f e++instance IsShade Shade' where+  shadeCtr f (Shade' c e) = fmap (`Shade'`e) $ f c+  occlusion = occ pseudoAffineWitness+   where occ :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x)+                        , Scalar (Needle x) ~ s, RealFloat' s )+                    => PseudoAffineWitness x -> Shade' x -> x -> s+         occ (PseudoAffineWitness (SemimanifoldWitness _)) (Shade' p₀ δinv) p+               = case toInterior p >>= (.-~.p₀) of+           (Just vd) | mSq <- normSq δinv vd+                     , mSq == mSq  -- avoid NaN+                     -> exp (negate mSq)+           _         -> zeroV+  factoriseShade (Shade' (x₀,y₀) δxy) = (Shade' x₀ δx, Shade' y₀ δy)+   where (δx,δy) = summandSpaceNorms δxy+  orthoShades (Shade' x δx) (Shade' y δy) = Shade' (x,y) $ sumSubspaceNorms δx δy+  coerceShade = cS+   where cS :: ∀ x y . (LocallyCoercible x y) => Shade' x -> Shade' y+         cS = \(Shade' x δxym) -> Shade' (internCoerce x) (tN δxym)+          where tN = case oppositeLocalCoercion :: CanonicalDiffeomorphism y x of+                      CanonicalDiffeomorphism ->+                       transformNorm . arr $ coerceNeedle ([]::[(y,x)])+                internCoerce = case interiorLocalCoercion ([]::[(x,y)]) of+                      CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism+  linIsoTransformShade = lits linearManifoldWitness linearManifoldWitness+                              dualSpaceWitness dualSpaceWitness+   where lits :: ∀ x y . ( SimpleSpace x, SimpleSpace y+                         , Scalar x ~ Scalar y, RealFloat' (Scalar x) )+               => LinearManifoldWitness x -> LinearManifoldWitness y+                   -> DualSpaceWitness x -> DualSpaceWitness y+                       -> (x+>y) -> Shade' x -> Shade' y+         lits (LinearManifoldWitness BoundarylessWitness)+              (LinearManifoldWitness BoundarylessWitness)+              DualSpaceWitness DualSpaceWitness+               f (Shade' x δx)+          = Shade' (f $ x) (transformNorm (pseudoInverse 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) )+                        => (SemimanifoldWitness x, SemimanifoldWitness y)+               -> Embedding (Affine s) (Interior x) (Interior y)+                              -> Shade' x -> Shade' y+         ps (SemimanifoldWitness _, SemimanifoldWitness _)+             (Embedding q p) (Shade' x e) = Shade' y (transformNorm j e)+          where y = q $ x+                (_,j) = evalAffine p y+  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) )+                        => (SemimanifoldWitness x, SemimanifoldWitness y)+               -> Embedding (Affine s) (Interior x) (Interior y)+                              -> Shade' y -> Shade' x+         ps (SemimanifoldWitness _, SemimanifoldWitness _)+             (Embedding p q) (Shade' x e) = Shade' y (transformNorm j e)+          where y = q $ x+                (_,j) = evalAffine p y+++shadeNarrowness :: Lens' (Shade' x) (Metric x)+shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e++instance ∀ x . (PseudoAffine x) => Semimanifold (Shade x) where+  type Needle (Shade x) = Needle x+  fromInterior = id+  toInterior = pure+  translateP = Tagged (.+~^)+  (.+~^) = case semimanifoldWitness :: SemimanifoldWitness x of+             SemimanifoldWitness BoundarylessWitness+                   -> \(Shade c e) v -> Shade (c.+~^v) e+  (.-~^) = case semimanifoldWitness :: SemimanifoldWitness x of+             SemimanifoldWitness BoundarylessWitness+                   -> \(Shade c e) v -> Shade (c.-~^v) e+  semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of+                         (SemimanifoldWitness BoundarylessWitness)+                          -> SemimanifoldWitness BoundarylessWitness++instance (WithField ℝ PseudoAffine x, Geodesic (Interior x), SimpleSpace (Needle x))+             => Geodesic (Shade x) where+  geodesicBetween = gb dualSpaceWitness+   where gb :: DualNeedleWitness x -> Shade x -> Shade x -> Maybe (D¹ -> Shade x)+         gb DualSpaceWitness (Shade c (Norm e)) (Shade ζ (Norm η)) = pure interp+          where interp t@(D¹ q) = Shade (pinterp t)+                                 (Norm . arr . lerp ed ηd $ (q+1)/2)+                ed@(LinearMap _) = arr e+                ηd@(LinearMap _) = arr η+                Just pinterp = geodesicBetween c ζ++instance (AffineManifold x) => Semimanifold (Shade' x) where+  type Needle (Shade' x) = Needle x+  fromInterior = id+  toInterior = pure+  translateP = Tagged (.+~^)+  (.+~^) = case boundarylessWitness :: BoundarylessWitness x of+      BoundarylessWitness -> \(Shade' c e) v -> Shade' (c.+~^v) e+  (.-~^) = case boundarylessWitness :: BoundarylessWitness x of+      BoundarylessWitness -> \(Shade' c e) v -> Shade' (c.-~^v) e+  semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of+     SemimanifoldWitness BoundarylessWitness -> SemimanifoldWitness BoundarylessWitness++instance ∀ x . (WithField ℝ AffineManifold x, Geodesic x, SimpleSpace (Needle x))+            => Geodesic (Shade' x) where+  geodesicBetween (Shade' c e) (Shade' ζ η) = pure interp+   where sharedSpan = sharedNormSpanningSystem e η+         interp t = Shade' (pinterp t)+                           (spanNorm [ v ^/ (alerpB 1 (recip qη) t)+                                     | (v,qη) <- sharedSpan ])+         Just pinterp = case geodesicWitness :: GeodesicWitness x of+            GeodesicWitness _ -> geodesicBetween c ζ++fullShade :: WithField ℝ PseudoAffine x => Interior x -> Metric' x -> Shade x+fullShade ctr expa = Shade ctr expa++fullShade' :: WithField ℝ PseudoAffine x => Interior x -> Metric x -> Shade' x+fullShade' ctr expa = Shade' ctr expa+++infixl 6 :±, |±|++-- | Span a 'Shade' from a center point and multiple deviation-vectors.+#if GLASGOW_HASKELL < 800+pattern (:±) :: ()+#else+pattern (:±) :: (WithField ℝ Manifold x, SimpleSpace (Needle x))+#endif+             => (WithField ℝ Manifold x, SimpleSpace (Needle x))+                         => Interior x -> [Needle x] -> Shade x+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.+--   Iff these form a orthogonal basis (in whatever sense applicable), then both+--   methods will be equivalent.+-- +--   Note that '|±|' is only possible, as such, in an inner-product space; in+--   general you need reciprocal vectors ('Needle'') to define a 'Shade''.+(|±|) :: ∀ x . WithField ℝ EuclidSpace x => x -> [Needle x] -> Shade' x+(|±|) = case boundarylessWitness :: BoundarylessWitness x of+   BoundarylessWitness -> \x shs -> Shade' x $ spanNorm [v^/(v<.>v) | v<-shs]++++                 +++-- | Attempt to find a 'Shade' that describes the distribution of given points.+--   At least in an affine space (and thus locally in any manifold), this can be used to+--   estimate the parameters of a normal distribution from which some points were+--   sampled. Note that some points will be &#x201c;outside&#x201d; of the shade,+--   as happens for a normal distribution with some statistical likelyhood.+--   (Use 'pointsCovers' if you need to prevent that.)+-- +--   For /nonconnected/ manifolds it will be necessary to yield separate shades+--   for each connected component. And for an empty input list, there is no shade!+--   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++coverAllAround :: ∀ x s . ( Fractional' s, WithField s PseudoAffine x+                          , SimpleSpace (Needle x) )+                  => Interior x -> [Needle x] -> Shade x+coverAllAround x₀ offs = Shade x₀+         $ guaranteeIn dualSpaceWitness offs+               (scaleNorm (1/fromIntegral (length offs)) $ spanVariance offs)+ where guaranteeIn :: DualNeedleWitness x -> [Needle x] -> Metric' x -> Metric' x+       guaranteeIn w@DualSpaceWitness offs ex+          = case offs >>= \v -> guard ((ex'|$|v) > 1) >> [(v, spanVariance [v])] of+             []   -> ex+             outs -> guaranteeIn w (fst<$>outs)+                                 ( densifyNorm $+                                    ex <> scaleNorm+                                                (sqrt . recip . fromIntegral+                                                            $ 2 * length outs)+                                                (mconcat $ snd<$>outs)+                                 )+        where ex' = dualNorm ex++-- | Like 'pointsShades', but ensure that all points are actually in+--   the shade, i.e. if @['Shade' x₀ ex]@ is the result then+--   @'metric' (recipMetric ex) (p-x₀) ≤ 1@ for all @p@ in the list.+pointsCovers :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))+                          => [Interior x] -> [Shade x]+pointsCovers = case pseudoAffineWitness :: PseudoAffineWitness x of+                 (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) ->+                  \ps -> map (\(ps', Shade x₀ _)+                                -> coverAllAround x₀ [v | p<-ps'+                                                        , let Just v+                                                                 = p.-~.fromInterior x₀])+                             (pointsShades' mempty (fromInterior<$>ps) :: [([x], Shade x)])++pointsShade's :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))+                     => [Interior x] -> [Shade' x]+pointsShade's = case dualSpaceWitness :: DualNeedleWitness x of+ DualSpaceWitness -> map (\(Shade c e :: Shade x) -> Shade' c $ dualNorm e) . pointsShades++pointsCover's :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))+                     => [Interior x] -> [Shade' x]+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 = case semimanifoldWitness :: SemimanifoldWitness x of+ SemimanifoldWitness _ ->+   \_ (p₀ NE.:| psr) -> foldl' ( \(acc, (rb,nr)) (i,p)+                                -> case (p.-~.acc, toInterior acc) of +                                      (Just δ, Just acci)+                                        -> (acci .+~^ δ^/i, (p:rb, nr))+                                      _ -> (acc, (rb, p:nr)) )+                             (p₀, mempty)+                             ( zip [1..] $ p₀:psr )++pointsShades' :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))+                                => Metric' x -> [x] -> [([x], Shade x)]+pointsShades' _ [] = []+pointsShades' minExt ps = case (expa, toInterior ctr) of +                           (Just e, Just c)+                             -> (ps, fullShade c e) : pointsShades' minExt unreachable+                           _ -> pointsShades' minExt inc'd+                                  ++ pointsShades' minExt unreachable+ where (ctr,(inc'd,unreachable)) = pseudoECM ([]::[x]) $ NE.fromList ps+       expa = ( (<>minExt) . spanVariance . map (^/ fromIntegral (length ps)) )+              <$> mapM (.-~.ctr) ps+       ++-- | Attempt to reduce the number of shades to fewer (ideally, a single one).+--   In the simplest cases these should guaranteed cover the same area;+--   for non-flat manifolds it only works in a heuristic sense.+shadesMerge :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x))+                 => ℝ -- ^ How near (inverse normalised distance, relative to shade expanse)+                      --   two shades must be to be merged. If this is zero, any shades+                      --   in the same connected region of a manifold are merged.+                 -> [Shade x] -- ^ A list of /n/ shades.+                 -> [Shade x] -- ^ /m/ &#x2264; /n/ shades which cover at least the same area.+shadesMerge fuzz (sh₁@(Shade c₁ e₁) : shs)+    = case extractJust (tryMerge pseudoAffineWitness dualSpaceWitness)+                 shs of+          (Just mg₁, shs') -> shadesMerge fuzz+                                $ shs'++[mg₁] -- Append to end to prevent undue weighting+                                              -- of first shade and its mergers.+          (_, shs') -> sh₁ : shadesMerge fuzz shs' + where tryMerge :: PseudoAffineWitness x -> DualNeedleWitness x+                         -> Shade x -> Maybe (Shade x)+       tryMerge (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) DualSpaceWitness+                    (Shade c₂ e₂)+           | Just v <- c₁.-~.c₂+           , [e₁',e₂'] <- dualNorm<$>[e₁, e₂] +           , b₁ <- e₂'|$|v+           , b₂ <- e₁'|$|v+           , fuzz*b₁*b₂ <= b₁ + b₂+                  = Just $ let cc = c₂ .+~^ v ^/ 2+                               Just cv₁ = c₁.-~.cc+                               Just cv₂ = c₂.-~.cc+                           in Shade cc $ e₁ <> e₂ <> spanVariance [cv₁, cv₂]+           | otherwise  = Nothing+shadesMerge _ shs = shs++-- | Weakened version of 'intersectShade's'. What this function calculates is+--   rather the /weighted mean/ of ellipsoid regions. If you interpret the+--   shades as uncertain physical measurements with normal distribution,+--   it gives the maximum-likelyhood result for multiple measurements of the+--   same quantity.+mixShade's :: ∀ y . (WithField ℝ Manifold y, SimpleSpace (Needle y))+                 => NonEmpty (Shade' y) -> Maybe (Shade' y)+mixShade's = ms pseudoAffineWitness dualSpaceWitness+ where ms :: PseudoAffineWitness y -> DualNeedleWitness y+                  -> NonEmpty (Shade' y) -> Maybe (Shade' y)+       ms (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness+                 (Shade' c₀ (Norm e₁):|shs) = sequenceA ciso >> pure mixed+        where ciso = [ci.-~.c₀ | Shade' ci shi <- shs]+              cis = [v | Just v <- ciso]+              σe = arr . sumV $ e₁ : (applyNorm . _shade'Narrowness<$>shs)+              cc = σe \$ sumV [ei $ ci | ci <- cis+                                       | Shade' _ (Norm ei) <- shs]+              mixed = Shade' (c₀+^cc) $ densifyNorm ( mconcat+                             [ Norm $ ei ^/ (1+(normSq ni $ ci^-^cc))+                             | ni@(Norm ei) <- Norm e₁ : (_shade'Narrowness<$>shs)+                             | ci <- zeroV : cis+                             ] )+              Tagged (+^) = translateP :: Tagged y (Interior y->Needle y->Interior y)+  -- cc should minimise the quadratic form+  -- β(cc) = ∑ᵢ ⟨cc−cᵢ|eᵢ|cc−cᵢ⟩+  -- = ⟨cc|e₁|cc⟩ + ∑ᵢ₌₁… ⟨cc−c₂|e₂|cc−c₂⟩+  -- = ⟨cc|e₁|cc⟩ + ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ − 2⋅⟨cᵢ|eᵢ|cc⟩ + ⟨cᵢ|eᵢ|cᵢ⟩ )+  -- It is thus+  -- β(cc + δ⋅v) − β cc+  -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩+  --     + ∑ᵢ₌₁…( ⟨cc + δ⋅v|eᵢ|cc + δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|cc + δ⋅v⟩ + ⟨cᵢ|eᵢ|cᵢ⟩ )+  --     − ⟨cc|e₁|cc⟩+  --     − ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ + 2⋅⟨cᵢ|eᵢ|cc⟩ − ⟨cᵢ|eᵢ|cᵢ⟩ )+  -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩+  --     + ∑ᵢ₌₁…( ⟨cc + δ⋅v|eᵢ|cc + δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|δ⋅v⟩ )+  --     − ⟨cc|e₁|cc⟩+  --     − ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ )+  -- = 2⋅⟨δ⋅v|e₁|cc⟩ + ⟨δ⋅v|e₁|δ⋅v⟩+  --     + ∑ᵢ₌₁…( 2⋅⟨δ⋅v|eᵢ|cc⟩ + ⟨δ⋅v|eᵢ|δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|δ⋅v⟩ )+  -- = 2⋅⟨δ⋅v|∑ᵢeᵢ|cc⟩ − 2⋅∑ᵢ₌₁… ⟨cᵢ|eᵢ|δ⋅v⟩ + 𝓞(δ²)+  -- This should vanish for all v, which is fulfilled by+  -- (∑ᵢeᵢ)|cc⟩ = ∑ᵢ₌₁… eᵢ|cᵢ⟩.++-- | Evaluate the shade as a quadratic form; essentially+-- @+-- minusLogOcclusion sh x = x <.>^ (sh^.shadeExpanse $ x - sh^.shadeCtr)+-- @+-- where 'shadeExpanse' gives a metric (matrix) that characterises the+-- width of the shade.+minusLogOcclusion' :: ∀ x s . ( PseudoAffine x, LinearSpace (Needle x)+                              , s ~ (Scalar (Needle x)), RealFloat' s )+              => Shade' x -> x -> s+minusLogOcclusion' (Shade' p₀ δinv)+        = occ (pseudoAffineWitness :: PseudoAffineWitness x)+              (dualSpaceWitness :: DualNeedleWitness x)+ where occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness+           p = case toInterior p >>= (.-~.p₀) of+         (Just vd) | mSq <- normSq δinv vd+                   , mSq == mSq  -- avoid NaN+                   -> mSq+         _         -> 1/0+minusLogOcclusion :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x)+                             , s ~ (Scalar (Needle x)), RealFloat' s )+              => Shade x -> x -> s+minusLogOcclusion (Shade p₀ δ)+        = occ (pseudoAffineWitness :: PseudoAffineWitness x)+              (dualSpaceWitness :: DualNeedleWitness x)+ where occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness+            = \p -> case toInterior p >>= (.-~.p₀) of+         (Just vd) | mSq <- normSq δinv vd+                   , mSq == mSq  -- avoid NaN+                   -> mSq+         _         -> 1/0+        where δinv = dualNorm δ+++++{-# WARNING rangeOnGeodesic "This function never worked properly. Use 'rangeWithinVertices'." #-}+rangeOnGeodesic :: ∀ i m . +      ( WithField ℝ PseudoAffine m, Geodesic m, SimpleSpace (Needle m)+      , WithField ℝ IntervalLike i, SimpleSpace (Needle i) )+                     => m -> m -> Maybe (Shade i -> Shade m)+rangeOnGeodesic = case ( semimanifoldWitness :: SemimanifoldWitness i+                       , dualSpaceWitness :: DualNeedleWitness i+                       , dualSpaceWitness :: DualNeedleWitness m ) of+ (SemimanifoldWitness _, DualSpaceWitness, DualSpaceWitness) ->+  \p₀ p₁ -> geodesicBetween p₀ p₁ >>=+      \interp -> case pointsShades =<<+                       [ mapMaybe (toInterior . interp . D¹) [-(1-ε), 1-ε]+                       | ε <- [0.0001, 0.001, 0.01, 0.1] ] of+                      defaultSh:_ -> Just $+                       \(Shade t₀ et) -> case pointsShades+                         . mapMaybe (toInterior+                               . interp . (toClosedInterval :: i -> D¹))+                         $ fromInterior <$> t₀ : [ t₀+^v+                                                 | v<-normSpanningSystem et ] of+                       [sh] -> sh+                       _ -> defaultSh+                      _ -> Nothing+ where Tagged (+^) = translateP :: Tagged i (Interior i->Needle i->Interior i)+++rangeWithinVertices :: ∀ s i m t+        . ( RealFrac' s+          , WithField s PseudoAffine i, WithField s PseudoAffine m+          , Geodesic i, Geodesic m+          , SimpleSpace (Needle i), SimpleSpace (Needle m)+          , AffineManifold (Interior i), AffineManifold (Interior m)+          , Object (Affine s) (Interior i), Object (Affine s) (Interior m)+          , Hask.Traversable t )+          => (Interior i,Interior m) -> t (i,m) -> Maybe (Shade i -> Shade m)+rangeWithinVertices+      = case ( semimanifoldWitness :: SemimanifoldWitness i+             , semimanifoldWitness :: SemimanifoldWitness m ) of+  (SemimanifoldWitness BoundarylessWitness, SemimanifoldWitness BoundarylessWitness)+      -> \(cii,cmi) verts ->+       let ci = fromInterior cii+           cm = fromInterior cmi+       in do+           vs <- sequenceA [ fzip ( middleBetween pi ci >>= (.-~.ci)+                                  , middleBetween pm cm >>= (.-~.cm) )+                           | (pi, pm) <- Hask.toList verts ]+           affinSys <- (correspondingDirections (cii,cmi) vs+                                 :: Maybe (Embedding (Affine (Scalar (Needle i)))+                                                     (Interior i) (Interior m)))+           return $ embedShade affinSys+          +++++data DebugView x where+  DebugView :: ( Show x, Show (Needle x+>Needle' x), LinearShowable (Needle x)+               , Needle' x ~ Needle x ) => DebugView x++-- | Class of manifolds which can use 'Shade'' as a basic set type.+--   This is easily possible for vector spaces with the default implementations.+class (WithField ℝ PseudoAffine y, SimpleSpace (Needle y)) => Refinable y where+  debugView :: Maybe (DebugView y)+  default debugView :: ( Show y, Show (Needle y+>Needle' y)+                       , Needle' y~Needle y, LinearShowable (Needle y) )+                         => Maybe (DebugView y)+  debugView = Just DebugView+  +  -- | @a `subShade'` b ≡ True@ means @a@ is fully contained in @b@, i.e. from+  --   @'minusLogOcclusion'' a p < 1@ follows also @minusLogOcclusion' b p < 1@.+  subShade' :: Shade' y -> Shade' y -> Bool+  subShade' (Shade' ac ae) (Shade' tc te)+        = case pseudoAffineWitness :: PseudoAffineWitness y of+   PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+    | Just v <- tc.-~.ac+    , v² <- normSq te v+    , v² <= 1+     -> all (\(y',μ) -> case μ of+            Nothing -> True  -- 'te' has infinite extension in this direction+            Just ξ+              | ξ<1 -> False -- 'ae' would be vaster than 'te' in this direction+              | ω <- abs $ y'<.>^v+                    -> (ω + 1/ξ)^2 <= 1 - v² + ω^2+                 -- See @images/constructions/subellipse-check-heuristic.svg@+         ) $ sharedSeminormSpanningSystem te ae+   _ -> False+  +  -- | Intersection between two shades.+  refineShade' :: Shade' y -> Shade' y -> Maybe (Shade' y)+  refineShade' (Shade' c₀ (Norm e₁)) (Shade' c₀₂ (Norm e₂))+      = case ( dualSpaceWitness :: DualNeedleWitness y+             , pseudoAffineWitness :: PseudoAffineWitness y ) of+          (DualSpaceWitness, PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))+               -> do+           c₂ <- c₀₂.-~.c₀+           let σe = arr $ e₁^+^e₂+               e₁c₂ = e₁ $ c₂+               e₂c₂ = e₂ $ c₂+               cc = σe \$ e₂c₂+               cc₂ = cc ^-^ c₂+               e₁cc = e₁ $ cc+               e₂cc = e₂ $ cc+               α = 2 + e₂c₂<.>^cc₂+           guard (α > 0)+           let ee = σe ^/ α+               c₂e₁c₂ = e₁c₂<.>^c₂+               c₂e₂c₂ = e₂c₂<.>^c₂+               c₂eec₂ = (c₂e₁c₂ + c₂e₂c₂) / α+           return $ case middle . sort+                $ quadraticEqnSol c₂e₁c₂+                                  (2 * (e₁cc<.>^c₂))+                                  (e₁cc<.>^cc - 1)+                ++quadraticEqnSol c₂e₂c₂+                                  (2 * (e₂cc<.>^c₂ - c₂e₂c₂))+                                  (e₂cc<.>^cc - 2 * (e₂c₂<.>^cc) + c₂e₂c₂ - 1) of+            [γ₁,γ₂] | abs (γ₁+γ₂) < 2 -> let+               cc' = cc ^+^ ((γ₁+γ₂)/2)*^c₂+               rγ = abs (γ₁ - γ₂) / 2+               η = if rγ * c₂eec₂ /= 0 && 1 - rγ^2 * c₂eec₂ > 0+                   then sqrt (1 - rγ^2 * c₂eec₂) / (rγ * c₂eec₂)+                   else 0+             in Shade' (c₀.+~^cc')+                       (Norm (arr ee) <> spanNorm [ee $ c₂^*η])+            _ -> Shade' (c₀.+~^cc) (Norm $ arr ee)+   where quadraticEqnSol a b c+             | a == 0, b /= 0       = [-c/b]+             | a /= 0 && disc == 0  = [- b / (2*a)]+             | a /= 0 && disc > 0   = [ (σ * sqrt disc - b) / (2*a)+                                      | σ <- [-1, 1] ]+             | otherwise            = []+          where disc = b^2 - 4*a*c+         middle (_:x:y:_) = [x,y]+         middle l = l+  -- ⟨x−c₁|e₁|x−c₁⟩ < 1  ∧  ⟨x−c₂|e₂|x−c₂⟩ < 1+  -- We search (cc,ee) such that this implies+  -- ⟨x−cc|ee|x−cc⟩ < 1.+  -- Let WLOG c₁ = 0, so+  -- ⟨x|e₁|x⟩ < 1.+  -- cc should minimise the quadratic form+  -- β(cc) = ⟨cc−c₁|e₁|cc−c₁⟩ + ⟨cc−c₂|e₂|cc−c₂⟩+  -- = ⟨cc|e₁|cc⟩ + ⟨cc−c₂|e₂|cc−c₂⟩+  -- = ⟨cc|e₁|cc⟩ + ⟨cc|e₂|cc⟩ − 2⋅⟨c₂|e₂|cc⟩ + ⟨c₂|e₂|c₂⟩+  -- It is thus+  -- β(cc + δ⋅v) − β cc+  -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + ⟨cc + δ⋅v|e₂|cc + δ⋅v⟩ − 2⋅⟨c₂|e₂|cc + δ⋅v⟩ + ⟨c₂|e₂|c₂⟩+  --     − ⟨cc|e₁|cc⟩ − ⟨cc|e₂|cc⟩ + 2⋅⟨c₂|e₂|cc⟩ − ⟨c₂|e₂|c₂⟩+  -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + ⟨cc + δ⋅v|e₂|cc + δ⋅v⟩ − 2⋅⟨c₂|e₂|δ⋅v⟩+  --     − ⟨cc|e₁|cc⟩ − ⟨cc|e₂|cc⟩+  -- = 2⋅⟨δ⋅v|e₁|cc⟩ + ⟨δ⋅v|e₁|δ⋅v⟩ + 2⋅⟨δ⋅v|e₂|cc⟩ + ⟨δ⋅v|e₂|δ⋅v⟩ − 2⋅⟨c₂|e₂|δ⋅v⟩+  -- = 2⋅δ⋅⟨v|e₁+e₂|cc⟩ − 2⋅δ⋅⟨v|e₂|c₂⟩ + 𝓞(δ²)+  -- This should vanish for all v, which is fulfilled by+  -- (e₁+e₂)|cc⟩ = e₂|c₂⟩.+  -- +  -- If we now choose+  -- ee = (e₁+e₂) / α+  -- then+  -- ⟨x−cc|ee|x−cc⟩ ⋅ α+  --  = ⟨x−cc|ee|x⟩ ⋅ α − ⟨x−cc|ee|cc⟩ ⋅ α+  --  = ⟨x|ee|x−cc⟩ ⋅ α − ⟨x−cc|e₂|c₂⟩+  --  = ⟨x|ee|x⟩ ⋅ α − ⟨x|ee|cc⟩ ⋅ α − ⟨x−cc|e₂|c₂⟩+  --  = ⟨x|e₁+e₂|x⟩ − ⟨x|e₂|c₂⟩ − ⟨x−cc|e₂|c₂⟩+  --  = ⟨x|e₁|x⟩ + ⟨x|e₂|x⟩ − ⟨x|e₂|c₂⟩ − ⟨x−cc|e₂|c₂⟩+  --  < 1 + ⟨x|e₂|x−c₂⟩ − ⟨x−cc|e₂|c₂⟩+  --  = 1 + ⟨x−c₂|e₂|x−c₂⟩ + ⟨c₂|e₂|x−c₂⟩ − ⟨x−cc|e₂|c₂⟩+  --  < 2 + ⟨x−c₂−x+cc|e₂|c₂⟩+  --  = 2 + ⟨cc−c₂|e₂|c₂⟩+  -- Really we want+  -- ⟨x−cc|ee|x−cc⟩ ⋅ α < α+  -- So choose α = 2 + ⟨cc−c₂|e₂|c₂⟩.+  -- +  -- The ellipsoid "cc±√ee" captures perfectly the intersection+  -- of the boundary of the shades, but it tends to significantly+  -- overshoot the interior intersection in perpendicular direction,+  -- i.e. in direction of c₂−c₁. E.g.+  -- https://github.com/leftaroundabout/manifolds/blob/bc0460b9/manifolds/images/examples/ShadeCombinations/EllipseIntersections.png+  -- 1. Really, the relevant points are those where either of the+  --    intersector badnesses becomes 1. The intersection shade should+  --    be centered between those points. We perform according corrections,+  --    but only in c₂ direction, so this can be handled efficiently+  --    as a 1D quadratic equation.+  --    Consider+  --       dⱼ c := ⟨c−cⱼ|eⱼ|c−cⱼ⟩ =! 1+  --       dⱼ (cc + γ⋅c₂)+  --           = ⟨cc+γ⋅c₂−cⱼ|eⱼ|cc+γ⋅c₂−cⱼ⟩+  --           = ⟨cc−cⱼ|eⱼ|cc−cⱼ⟩ + 2⋅γ⋅⟨c₂|eⱼ|cc−cⱼ⟩ + γ²⋅⟨c₂|eⱼ|c₂⟩+  --           =! 1+  --    So+  --    γⱼ = (- b ± √(b²−4⋅a⋅c)) / 2⋅a+  --     where a = ⟨c₂|eⱼ|c₂⟩+  --           b = 2 ⋅ (⟨c₂|eⱼ|cc⟩ − ⟨c₂|eⱼ|cⱼ⟩)+  --           c = ⟨cc|eⱼ|cc⟩ − 2⋅⟨cc|eⱼ|cⱼ⟩ + ⟨cⱼ|eⱼ|cⱼ⟩ − 1+  --    The ± sign should be chosen to get the smaller |γ| (otherwise+  --    we end up on the wrong side of the shade), i.e.+  --    γⱼ = (sgn bⱼ ⋅ √(bⱼ²−4⋅aⱼ⋅cⱼ) − bⱼ) / 2⋅aⱼ+  -- 2. Trim the result in that direction to the actual+  --    thickness of the lens-shaped intersection: we want+  --    ⟨rγ⋅c₂|ee'|rγ⋅c₂⟩ = 1+  --    for a squeezed version of ee,+  --    ee' = ee + ee|η⋅c₂⟩⟨η⋅c₂|ee+  --    ee' = ee + η² ⋅ ee|c₂⟩⟨c₂|ee+  --    ⟨rγ⋅c₂|ee'|rγ⋅c₂⟩+  --        = rγ² ⋅ (⟨c₂|ee|c₂⟩ + η² ⋅ ⟨c₂|ee|c₂⟩²)+  --        = rγ² ⋅ ⟨c₂|ee|c₂⟩ + η² ⋅ rγ² ⋅ ⟨c₂|ee|c₂⟩²+  --    η² = (1 − rγ²⋅⟨c₂|ee|c₂⟩) / (rγ² ⋅ ⟨c₂|ee|c₂⟩²)+  --    η = √(1 − rγ²⋅⟨c₂|ee|c₂⟩) / (rγ ⋅ ⟨c₂|ee|c₂⟩)+  --    With ⟨c₂|ee|c₂⟩ = (⟨c₂|e₁|c₂⟩ + ⟨c₂|e₂|c₂⟩)/α.++  +  -- | If @p@ is in @a@ (red) and @δ@ is in @b@ (green),+  --   then @p.+~^δ@ is in @convolveShade' a b@ (blue).+  -- +--   Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/ShadeCombinations.ipynb#shadeConvolutions+-- +-- <<images/examples/ShadeCombinations/2Dconvolution-skewed.png>>+  convolveMetric :: Hask.Functor p => p y -> Metric y -> Metric y -> Metric y+  convolveMetric _ ey eδ = case wellDefinedNorm result of+          Just r  -> r+          Nothing -> case debugView :: Maybe (DebugView y) of+            Just DebugView -> error $ "Can not convolve norms "+                               ++show (arr (applyNorm ey) :: Needle y+>Needle' y)+                               ++" and "++show (arr (applyNorm eδ) :: Needle y+>Needle' y)+   where eδsp = sharedSeminormSpanningSystem ey eδ+         result = spanNorm [ f ^* ζ crl | (f,crl) <- eδsp ]+         ζ = case filter (>0) . catMaybes $ snd<$>eδsp of+            [] -> const 0+            nzrelap+               -> let cre₁ = 1/minimum nzrelap+                      cre₂ =  maximum nzrelap+                      edgeFactor = sqrt ( (1 + cre₁)^2 + (1 + cre₂)^2 )+                                / (sqrt (1 + cre₁^2) + sqrt (1 + cre₂^2))+                  in \case+                        Nothing -> 0+                        Just 0  -> 0+                        Just sq -> edgeFactor / (recip sq + 1)+  +  convolveShade' :: Shade' y -> Shade' (Needle y) -> Shade' y+  convolveShade' = defaultConvolveShade'+  +defaultConvolveShade' :: ∀ y . Refinable y => Shade' y -> Shade' (Needle y) -> Shade' y+defaultConvolveShade' = case (pseudoAffineWitness :: PseudoAffineWitness y) of+  PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+    -> \(Shade' y₀ ey) (Shade' δ₀ eδ) -> Shade' (y₀.+~^δ₀)+                                          $ convolveMetric ([]::[y]) ey eδ++instance Refinable ℝ where+  refineShade' (Shade' cl el) (Shade' cr er)+         = case (normSq el 1, normSq er 1) of+             (0, _) -> return $ Shade' cr er+             (_, 0) -> return $ Shade' cl el+             (ql,qr) | ql>0, qr>0+                    -> let [rl,rr] = sqrt . recip <$> [ql,qr]+                           b = maximum $ zipWith (-) [cl,cr] [rl,rr]+                           t = minimum $ zipWith (+) [cl,cr] [rl,rr]+                       in guard (b<t) >>+                           let cm = (b+t)/2+                               rm = (t-b)/2+                           in return $ Shade' cm (spanNorm [recip rm])+--   convolveShade' (Shade' y₀ ey) (Shade' δ₀ eδ)+--          = case (metricSq ey 1, metricSq eδ 1) of+--              (wy,wδ) | wy>0, wδ>0+--                  -> Shade' (y₀.+~^δ₀)+--                            ( projector . recip+--                                   $ recip (sqrt wy) + recip (sqrt wδ) )+--              (_ , _) -> Shade' y₀ zeroV++instance ∀ a b . ( Refinable a, Refinable b+                 , Scalar (DualVector (DualVector (Needle b)))+                      ~ Scalar (DualVector (DualVector (Needle a))) )+    => Refinable (a,b) where+  debugView = case ( debugView :: Maybe (DebugView a)+                   , debugView :: Maybe (DebugView b)+                   , dualSpaceWitness :: DualSpaceWitness (Needle a)+                   , dualSpaceWitness :: DualSpaceWitness (Needle b) ) of+      (Just DebugView, Just DebugView, DualSpaceWitness, DualSpaceWitness)+              -> Just DebugView+  +instance Refinable ℝ⁰+instance Refinable ℝ¹+instance Refinable ℝ²+instance Refinable ℝ³+instance Refinable ℝ⁴+                            +instance ( SimpleSpace a, SimpleSpace b+         , Refinable a, Refinable b+         , Scalar a ~ ℝ, Scalar b ~ ℝ+         , Scalar (DualVector a) ~ ℝ, Scalar (DualVector b) ~ ℝ+         , Scalar (DualVector (DualVector a)) ~ ℝ, Scalar (DualVector (DualVector b)) ~ ℝ )+            => Refinable (LinearMap ℝ a b) where+  debugView = Nothing++intersectShade's :: ∀ y . Refinable y => NonEmpty (Shade' y) -> Maybe (Shade' y)+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+      = WithAny { _untopological :: y+                , _topological :: !x  }+ deriving (Hask.Functor, Show, Generic)++instance (NFData x, NFData y) => NFData (WithAny x y)++instance ∀ x y . (Semimanifold x) => Semimanifold (x`WithAny`y) where+  type Needle (WithAny x y) = Needle x+  type Interior (WithAny x y) = Interior x `WithAny` y+  WithAny y x .+~^ δx = WithAny y $ x.+~^δx+  fromInterior (WithAny y x) = WithAny y $ fromInterior x+  toInterior (WithAny y x) = fmap (WithAny y) $ toInterior x+  translateP = tpWD+   where tpWD :: ∀ x y . Semimanifold x => Tagged (WithAny x y)+                            (Interior x`WithAny`y -> Needle x -> Interior x`WithAny`y)+         tpWD = Tagged `id` \(WithAny y x) δx -> WithAny y $ tpx x δx+          where Tagged tpx = translateP :: Tagged x (Interior x -> Needle x -> Interior x)+  semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of+      SemimanifoldWitness BoundarylessWitness -> SemimanifoldWitness BoundarylessWitness+            +instance (PseudoAffine x) => PseudoAffine (x`WithAny`y) where+  WithAny _ x .-~. WithAny _ ξ = x.-~.ξ+  pseudoAffineWitness = case pseudoAffineWitness :: PseudoAffineWitness x of+      PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)+       -> PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)++instance (AffineSpace x) => AffineSpace (x`WithAny`y) where+  type Diff (WithAny x y) = Diff x+  WithAny _ x .-. WithAny _ ξ = x.-.ξ+  WithAny y x .+^ δx = WithAny y $ x.+^δx ++instance (VectorSpace x, Monoid y) => VectorSpace (x`WithAny`y) where+  type Scalar (WithAny x y) = Scalar x+  μ *^ WithAny y x = WithAny y $ μ*^x ++instance (AdditiveGroup x, Monoid y) => AdditiveGroup (x`WithAny`y) where+  zeroV = WithAny mempty zeroV+  negateV (WithAny y x) = WithAny y $ negateV x+  WithAny y x ^+^ WithAny υ ξ = WithAny (mappend y υ) (x^+^ξ)++instance (AdditiveGroup x) => Hask.Applicative (WithAny x) where+  pure x = WithAny x zeroV+  WithAny f x <*> WithAny t ξ = WithAny (f t) (x^+^ξ)+  +instance (AdditiveGroup x) => Hask.Monad (WithAny x) where+  return x = WithAny x zeroV+  WithAny y x >>= f = WithAny r $ x^+^q+   where WithAny r q = f y++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 (Shade (WithAny _ b) e) = Shade b e++                      +++++extractJust :: (a->Maybe b) -> [a] -> (Maybe b, [a])+extractJust f [] = (Nothing,[])+extractJust f (x:xs) | Just r <- f x  = (Just r, xs)+                     | otherwise      = second (x:) $ extractJust f xs+++prettyShowShade' :: LtdErrorShow x => Shade' x -> String+prettyShowShade' sh = prettyShowsPrecShade' 0 sh []++++wellDefinedShade' :: LinearSpace (Needle x) => Shade' x -> Maybe (Shade' x)+wellDefinedShade' (Shade' c e) = Shade' c <$> wellDefinedNorm e++++data LtdErrorShowWitness m where+   LtdErrorShowWitness :: (LtdErrorShow (Interior m), LtdErrorShow (Needle m))+                  => PseudoAffineWitness m -> LtdErrorShowWitness m++class Refinable m => LtdErrorShow m where+  ltdErrorShowWitness :: LtdErrorShowWitness m+  default ltdErrorShowWitness :: (LtdErrorShow (Interior m), LtdErrorShow (Needle m))+                         => 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 sh+         u :: [ShowS] = case ltdErrorShowWitness :: LtdErrorShowWitness m of+           LtdErrorShowWitness (PseudoAffineWitness (SemimanifoldWitness _)) ->+             [ showsPrecShade'_errorLtdC 6 (Shade' δ e :: Shade' (Needle m))+             | δ <- varianceSpanningSystem e']+         e' = dualNorm e++instance LtdErrorShow ℝ⁰ where+  showsPrecShade'_errorLtdC _ _ = ("zeroV"++)+instance LtdErrorShow ℝ where+  showsPrecShade'_errorLtdC _ (Shade' v u) = errorLtdShow (δ/2) v+   where δ = case u<$|1 of+          σ | σ>0 -> sqrt $ 1/σ+          _       -> v*10+instance LtdErrorShow ℝ² where+  showsPrecShade'_errorLtdC _ sh = ("V2 "++) . shshx . (' ':) . shshy+   where shx = projectShade (lensEmbedding _x) sh :: Shade' ℝ+         shy = projectShade (lensEmbedding _y) sh :: Shade' ℝ+         shshx = showsPrecShade'_errorLtdC 0 shx +         shshy = showsPrecShade'_errorLtdC 0 shy +instance LtdErrorShow ℝ³ where+  showsPrecShade'_errorLtdC _ sh = ("V3 "++) . shshx . (' ':) . shshy . (' ':) . shshz+   where shx = projectShade (lensEmbedding _x) sh :: Shade' ℝ+         shy = projectShade (lensEmbedding _y) sh :: Shade' ℝ+         shz = projectShade (lensEmbedding _z) sh :: Shade' ℝ+         shshx = showsPrecShade'_errorLtdC 0 shx +         shshy = showsPrecShade'_errorLtdC 0 shy +         shshz = showsPrecShade'_errorLtdC 0 shz +instance LtdErrorShow ℝ⁴ where+  showsPrecShade'_errorLtdC _ sh+           = ("V4 "++) . shshx . (' ':) . shshy . (' ':) . shshz . (' ':) . shshw+   where shx = projectShade (lensEmbedding _x) sh :: Shade' ℝ+         shy = projectShade (lensEmbedding _y) sh :: Shade' ℝ+         shz = projectShade (lensEmbedding _z) sh :: Shade' ℝ+         shw = projectShade (lensEmbedding _w) sh :: Shade' ℝ+         shshx = showsPrecShade'_errorLtdC 0 shx +         shshy = showsPrecShade'_errorLtdC 0 shy +         shshz = showsPrecShade'_errorLtdC 0 shz +         shshw = showsPrecShade'_errorLtdC 0 shw +instance ∀ x y .+         ( LtdErrorShow x, LtdErrorShow y+         , Scalar (DualVector (Needle' x)) ~ Scalar (DualVector (Needle' y)) )+              => LtdErrorShow (x,y) where+  ltdErrorShowWitness = case ( ltdErrorShowWitness :: LtdErrorShowWitness x+                             , ltdErrorShowWitness :: LtdErrorShowWitness y ) of+   (  LtdErrorShowWitness(PseudoAffineWitness(SemimanifoldWitness BoundarylessWitness))+    , LtdErrorShowWitness(PseudoAffineWitness(SemimanifoldWitness BoundarylessWitness)) )+    ->LtdErrorShowWitness(PseudoAffineWitness(SemimanifoldWitness BoundarylessWitness))+  showsPrecShade'_errorLtdC _ sh = ('(':) . shshx . (',':) . shshy . (')':)+   where (shx,shy) = factoriseShade sh+         shshx = showsPrecShade'_errorLtdC 0 shx +         shshy = showsPrecShade'_errorLtdC 0 shy +                       +instance LtdErrorShow x => Show (Shade' x) where+  showsPrec = prettyShowsPrecShade'

@@ -14,26 +14,17 @@ {-# LANGUAGE DeriveFunctor              #-} {-# LANGUAGE DeriveFoldable             #-} {-# LANGUAGE DeriveTraversable          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE TypeFamilies               #-} {-# LANGUAGE FunctionalDependencies     #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE GADTs                      #-} {-# LANGUAGE RankNTypes                 #-} {-# LANGUAGE TupleSections              #-}-{-# LANGUAGE ParallelListComp           #-}-{-# LANGUAGE MonadComprehensions        #-} {-# LANGUAGE UnicodeSyntax              #-}-{-# LANGUAGE ConstraintKinds            #-}-{-# LANGUAGE PatternGuards              #-} {-# LANGUAGE PatternSynonyms            #-}-{-# LANGUAGE ViewPatterns               #-} {-# LANGUAGE LambdaCase                 #-} {-# LANGUAGE TypeOperators              #-}-{-# LANGUAGE CPP                        #-} {-# LANGUAGE ScopedTypeVariables        #-}-{-# LANGUAGE LiberalTypeSynonyms        #-}-{-# LANGUAGE RecordWildCards            #-} {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE TemplateHaskell            #-} @@ -47,9 +38,10 @@        , fullShade, fullShade', pointsShades, pointsShade's        , pointsCovers, pointsCover's, coverAllAround        -- ** Evaluation-       , occlusion+       , occlusion, prettyShowsPrecShade', prettyShowShade'        -- ** Misc        , factoriseShade, intersectShade's, linIsoTransformShade+       , embedShade, projectShade        , Refinable, subShade', refineShade', convolveShade', coerceShade        , mixShade's        -- * Shade trees@@ -72,6 +64,8 @@        -- ** Triangulation-builders        , TriangBuild, doTriangBuild        , AutoTriang, breakdownAutoTriang+       -- ** External+       , AffineManifold, euclideanMetric     ) where  @@ -92,10 +86,13 @@ import Data.Tagged  import Data.SimplicialComplex+import Data.Manifold.Shade import Data.Manifold.Types import Data.Manifold.Types.Primitive ((^), empty) import Data.Manifold.PseudoAffine import Data.Manifold.Riemannian+import Data.Manifold.Atlas+import Data.Function.Affine      import Data.Embedding import Data.CoNat@@ -127,51 +124,7 @@ import Data.Type.Coercion  --- | Possibly / Partially / asymPtotically singular metric.-data PSM x = PSM {-       psmExpanse :: !(Metric' x)-     , relevantEigenspan :: ![Needle' x]-     }-        --- | A 'Shade' is a very crude description of a region within a manifold. It---   can be interpreted as either an ellipsoid shape, or as the Gaussian peak---   of a normal distribution (use <http://hackage.haskell.org/package/manifold-random>---   for actually sampling from that distribution).--- ---   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) }-deriving instance (Show (Interior x), Show (Metric' x), WithField ℝ PseudoAffine x)-                => Show (Shade x)---- | A &#x201c;co-shade&#x201d; can describe ellipsoid regions as well, but unlike---   'Shade' it can be unlimited / infinitely wide in some directions.---   It does OTOH need to have nonzero thickness, which 'Shade' needs not.-data Shade' x = Shade' { _shade'Ctr :: !(Interior x)-                       , _shade'Narrowness :: !(Metric x) }-deriving instance (Show (Interior x), Show (Metric x), WithField ℝ PseudoAffine x)-                => Show (Shade' x)--data LocalDifferentialEqn x y = LocalDifferentialEqn {-      _predictDerivatives :: Maybe (Shade' (LocalLinear x y))-    , _rescanDerivatives :: Shade' (LocalLinear x y) -> Shade' y -> Maybe (Shade' 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- type Depth = Int data Wall x = Wall { _wallID :: (Depth,(Int,Int))                    , _wallAnchor :: Interior x@@ -181,182 +134,6 @@ makeLenses ''Wall  -class IsShade shade where---  type (*) shade :: *->*-  -- | Access the center of a 'Shade' or a 'Shade''.-  shadeCtr :: Lens' (shade x) (Interior x)---  -- | Convert between 'Shade' and 'Shade' (which must be neither singular nor infinite).---  unsafeDualShade :: WithField ℝ Manifold x => shade x -> shade* x-  -- | Check the statistical likelihood-density of a point being within a shade.-  --   This is taken as a normal distribution.-  occlusion :: ( PseudoAffine x, SimpleSpace (Needle x)-               , s ~ (Scalar (Needle x)), RealDimension s )-                => shade x -> x -> s-  factoriseShade :: ( Manifold x, SimpleSpace (Needle x)-                    , Manifold 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-  linIsoTransformShade :: ( LinearManifold x, LinearManifold y-                          , SimpleSpace x, SimpleSpace y, Scalar x ~ Scalar y )-                          => (x+>y) -> shade x -> shade y--instance IsShade Shade where-  shadeCtr f (Shade c e) = fmap (`Shade`e) $ f c-  occlusion = occ pseudoAffineWitness dualSpaceWitness-   where occ :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x)-                        , Scalar (Needle x) ~ s, RealDimension s )-                    => PseudoAffineWitness x -> DualNeedleWitness x -> Shade x -> x -> s-         occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness (Shade p₀ δ)-                 = \p -> case toInterior p >>= (.-~.p₀) of-           (Just vd) | mSq <- normSq δinv vd-                     , mSq == mSq  -- avoid NaN-                     -> exp (negate mSq)-           _         -> zeroV-          where δinv = dualNorm δ-  factoriseShade = fs dualSpaceWitness dualSpaceWitness-   where fs :: ∀ x y . ( Manifold x, SimpleSpace (Needle x)-                       , Manifold y, SimpleSpace (Needle y)-                       , Scalar (Needle x) ~ Scalar (Needle y) )-               => DualNeedleWitness x -> DualNeedleWitness y-                       -> Shade (x,y) -> (Shade x, Shade y)-         fs DualSpaceWitness DualSpaceWitness (Shade (x₀,y₀) δxy)-                   = (Shade x₀ δx, Shade y₀ δy)-          where (δx,δy) = summandSpaceNorms δxy-  coerceShade = cS dualSpaceWitness dualSpaceWitness-   where cS :: ∀ x y . (LocallyCoercible x y)-                => DualNeedleWitness x -> DualNeedleWitness y -> Shade x -> Shade y-         cS DualSpaceWitness DualSpaceWitness-                    = \(Shade x δxym) -> Shade (internCoerce x) (tN δxym)-          where tN = case oppositeLocalCoercion :: CanonicalDiffeomorphism y x of-                      CanonicalDiffeomorphism ->-                       transformNorm . arr $ coerceNeedle' ([]::[(y,x)])-                internCoerce = case interiorLocalCoercion ([]::[(x,y)]) of-                      CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism-  linIsoTransformShade = lits dualSpaceWitness dualSpaceWitness-   where lits :: ∀ x y . ( LinearManifold x, LinearManifold y-                         , Scalar (Needle x) ~ Scalar (Needle y) )-               => DualSpaceWitness x -> DualSpaceWitness y-                       -> (x+>y) -> Shade x -> Shade y-         lits DualSpaceWitness DualSpaceWitness f (Shade x δx)-                  = Shade (f $ x) (transformNorm (adjoint $ f) δx)--instance ImpliesMetric Shade where-  type MetricRequirement Shade x = (Manifold x, SimpleSpace (Needle x))-  inferMetric' (Shade _ e) = e-  inferMetric = im dualSpaceWitness-   where im :: (Manifold x, SimpleSpace (Needle x))-                   => DualNeedleWitness x -> Shade x -> Metric x-         im DualSpaceWitness (Shade _ e) = dualNorm e--instance ImpliesMetric Shade' where-  type MetricRequirement Shade' x = (Manifold x, SimpleSpace (Needle x))-  inferMetric (Shade' _ e) = e-  inferMetric' (Shade' _ e) = dualNorm e--shadeExpanse :: Lens' (Shade x) (Metric' x)-shadeExpanse f (Shade c e) = fmap (Shade c) $ f e--instance IsShade Shade' where-  shadeCtr f (Shade' c e) = fmap (`Shade'`e) $ f c-  occlusion = occ pseudoAffineWitness-   where occ :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x)-                        , Scalar (Needle x) ~ s, RealDimension s )-                    => PseudoAffineWitness x -> Shade' x -> x -> s-         occ (PseudoAffineWitness (SemimanifoldWitness _)) (Shade' p₀ δinv) p-               = case toInterior p >>= (.-~.p₀) of-           (Just vd) | mSq <- normSq δinv vd-                     , mSq == mSq  -- avoid NaN-                     -> exp (negate mSq)-           _         -> zeroV-  factoriseShade (Shade' (x₀,y₀) δxy) = (Shade' x₀ δx, Shade' y₀ δy)-   where (δx,δy) = summandSpaceNorms δxy-  coerceShade = cS-   where cS :: ∀ x y . (LocallyCoercible x y) => Shade' x -> Shade' y-         cS = \(Shade' x δxym) -> Shade' (internCoerce x) (tN δxym)-          where tN = case oppositeLocalCoercion :: CanonicalDiffeomorphism y x of-                      CanonicalDiffeomorphism ->-                       transformNorm . arr $ coerceNeedle ([]::[(y,x)])-                internCoerce = case interiorLocalCoercion ([]::[(x,y)]) of-                      CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism-  linIsoTransformShade f (Shade' x δx)-          = Shade' (f $ x) (transformNorm (pseudoInverse f) δx)--shadeNarrowness :: Lens' (Shade' x) (Metric x)-shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e--instance ∀ x . (PseudoAffine x) => Semimanifold (Shade x) where-  type Needle (Shade x) = Needle x-  fromInterior = id-  toInterior = pure-  translateP = Tagged (.+~^)-  (.+~^) = case semimanifoldWitness :: SemimanifoldWitness x of-             SemimanifoldWitness BoundarylessWitness-                   -> \(Shade c e) v -> Shade (c.+~^v) e-  (.-~^) = case semimanifoldWitness :: SemimanifoldWitness x of-             SemimanifoldWitness BoundarylessWitness-                   -> \(Shade c e) v -> Shade (c.-~^v) e-  semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of-                         (SemimanifoldWitness BoundarylessWitness)-                          -> SemimanifoldWitness BoundarylessWitness--instance (WithField ℝ PseudoAffine x, Geodesic (Interior x), SimpleSpace (Needle x))-             => Geodesic (Shade x) where-  geodesicBetween = gb dualSpaceWitness-   where gb :: DualNeedleWitness x -> Shade x -> Shade x -> Maybe (D¹ -> Shade x)-         gb DualSpaceWitness (Shade c (Norm e)) (Shade ζ (Norm η)) = pure interp-          where interp t@(D¹ q) = Shade (pinterp t)-                                 (Norm . arr . lerp ed ηd $ (q+1)/2)-                ed@(LinearMap _) = arr e-                ηd@(LinearMap _) = arr η-                Just pinterp = geodesicBetween c ζ--instance (AffineManifold x) => Semimanifold (Shade' x) where-  type Needle (Shade' x) = Diff x-  fromInterior = id-  toInterior = pure-  translateP = Tagged (.+~^)-  Shade' c e .+~^ v = Shade' (c.+^v) e-  Shade' c e .-~^ v = Shade' (c.-^v) e--instance (WithField ℝ AffineManifold x, Geodesic x, SimpleSpace (Needle x))-            => Geodesic (Shade' x) where-  geodesicBetween (Shade' c e) (Shade' ζ η) = pure interp-   where sharedSpan = sharedNormSpanningSystem e η-         interp t = Shade' (pinterp t)-                           (spanNorm [ v ^/ (alerpB 1 (recip qη) t)-                                     | (v,qη) <- sharedSpan ])-         Just pinterp = geodesicBetween c ζ--fullShade :: WithField ℝ PseudoAffine x => Interior x -> Metric' x -> Shade x-fullShade ctr expa = Shade ctr expa--fullShade' :: WithField ℝ PseudoAffine x => Interior x -> Metric x -> Shade' x-fullShade' ctr expa = Shade' ctr expa----- | Span a 'Shade' from a center point and multiple deviation-vectors.-#if GLASGOW_HASKELL < 800-pattern (:±) :: ()-#else-pattern (:±) :: (WithField ℝ Manifold x, SimpleSpace (Needle x))-#endif-             => (WithField ℝ Manifold x, SimpleSpace (Needle x))-                         => Interior x -> [Needle x] -> Shade x-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.---   Iff these form a orthogonal basis (in whatever sense applicable), then both---   methods will be equivalent.--- ---   Note that '|±|' is only possible, as such, in an inner-product space; in---   general you need reciprocal vectors ('Needle'') to define a 'Shade''.-(|±|) :: WithField ℝ EuclidSpace x => x -> [Needle x] -> Shade' x-x |±| shs = Shade' x $ spanNorm [v^/(v<.>v) | v<-shs]--- subshadeId' :: ∀ x . (WithField ℝ PseudoAffine x, LinearSpace (Needle x))                    => x -> NonEmpty (Needle' x) -> x -> (Int, HourglassBulb) subshadeId' c expvs x = case ( dualSpaceWitness :: DualNeedleWitness x@@ -375,220 +152,9 @@                    --- | Attempt to find a 'Shade' that describes the distribution of given points.---   At least in an affine space (and thus locally in any manifold), this can be used to---   estimate the parameters of a normal distribution from which some points were---   sampled. Note that some points will be &#x201c;outside&#x201d; of the shade,---   as happens for a normal distribution with some statistical likelyhood.---   (Use 'pointsCovers' if you need to prevent that.)--- ---   For /nonconnected/ manifolds it will be necessary to yield separate shades---   for each connected component. And for an empty input list, there is no shade!---   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 -coverAllAround :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                  => Interior x -> [Needle x] -> Shade x-coverAllAround x₀ offs = Shade x₀-         $ guaranteeIn dualSpaceWitness offs-               (scaleNorm (1/fromIntegral (length offs)) $ spanVariance offs)- where guaranteeIn :: DualNeedleWitness x -> [Needle x] -> Metric' x -> Metric' x-       guaranteeIn w@DualSpaceWitness offs ex-          = case offs >>= \v -> guard ((ex'|$|v) > 1) >> [(v, spanVariance [v])] of-             []   -> ex-             outs -> guaranteeIn w (fst<$>outs)-                                 ( densifyNorm $-                                    ex <> scaleNorm-                                                (sqrt . recip . fromIntegral-                                                            $ 2 * length outs)-                                                (mconcat $ snd<$>outs)-                                 )-        where ex' = dualNorm ex --- | Like 'pointsShades', but ensure that all points are actually in---   the shade, i.e. if @['Shade' x₀ ex]@ is the result then---   @'metric' (recipMetric ex) (p-x₀) ≤ 1@ for all @p@ in the list.-pointsCovers :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                          => [Interior x] -> [Shade x]-pointsCovers = case pseudoAffineWitness :: PseudoAffineWitness x of-                 (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) ->-                  \ps -> map (\(ps', Shade x₀ _)-                                -> coverAllAround x₀ [v | p<-ps'-                                                        , let Just v-                                                                 = p.-~.fromInterior x₀])-                             (pointsShades' mempty (fromInterior<$>ps) :: [([x], Shade x)]) -pointsShade's :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                     => [Interior x] -> [Shade' x]-pointsShade's = case dualSpaceWitness :: DualNeedleWitness x of- DualSpaceWitness -> map (\(Shade c e :: Shade x) -> Shade' c $ dualNorm e) . pointsShades--pointsCover's :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                     => [Interior x] -> [Shade' x]-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 = case semimanifoldWitness :: SemimanifoldWitness x of- SemimanifoldWitness _ ->-   \_ (p₀ NE.:| psr) -> foldl' ( \(acc, (rb,nr)) (i,p)-                                -> case (p.-~.acc, toInterior acc) of -                                      (Just δ, Just acci)-                                        -> (acci .+~^ δ^/i, (p:rb, nr))-                                      _ -> (acc, (rb, p:nr)) )-                             (p₀, mempty)-                             ( zip [1..] $ p₀:psr )--pointsShades' :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x))-                                => Metric' x -> [x] -> [([x], Shade x)]-pointsShades' _ [] = []-pointsShades' minExt ps = case (expa, toInterior ctr) of -                           (Just e, Just c)-                             -> (ps, fullShade c e) : pointsShades' minExt unreachable-                           _ -> pointsShades' minExt inc'd-                                  ++ pointsShades' minExt unreachable- where (ctr,(inc'd,unreachable)) = pseudoECM ([]::[x]) $ NE.fromList ps-       expa = ( (<>minExt) . spanVariance . map (^/ fromIntegral (length ps)) )-              <$> mapM (.-~.ctr) ps-       ---- | Attempt to reduce the number of shades to fewer (ideally, a single one).---   In the simplest cases these should guaranteed cover the same area;---   for non-flat manifolds it only works in a heuristic sense.-shadesMerge :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x))-                 => ℝ -- ^ How near (inverse normalised distance, relative to shade expanse)-                      --   two shades must be to be merged. If this is zero, any shades-                      --   in the same connected region of a manifold are merged.-                 -> [Shade x] -- ^ A list of /n/ shades.-                 -> [Shade x] -- ^ /m/ &#x2264; /n/ shades which cover at least the same area.-shadesMerge fuzz (sh₁@(Shade c₁ e₁) : shs)-    = case extractJust (tryMerge pseudoAffineWitness dualSpaceWitness)-                 shs of-          (Just mg₁, shs') -> shadesMerge fuzz-                                $ shs'++[mg₁] -- Append to end to prevent undue weighting-                                              -- of first shade and its mergers.-          (_, shs') -> sh₁ : shadesMerge fuzz shs' - where tryMerge :: PseudoAffineWitness x -> DualNeedleWitness x-                         -> Shade x -> Maybe (Shade x)-       tryMerge (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) DualSpaceWitness-                    (Shade c₂ e₂)-           | Just v <- c₁.-~.c₂-           , [e₁',e₂'] <- dualNorm<$>[e₁, e₂] -           , b₁ <- e₂'|$|v-           , b₂ <- e₁'|$|v-           , fuzz*b₁*b₂ <= b₁ + b₂-                  = Just $ let cc = c₂ .+~^ v ^/ 2-                               Just cv₁ = c₁.-~.cc-                               Just cv₂ = c₂.-~.cc-                           in Shade cc $ e₁ <> e₂ <> spanVariance [cv₁, cv₂]-           | otherwise  = Nothing-shadesMerge _ shs = shs---- | Weakened version of 'intersectShade's'. What this function calculates is---   rather the /weighted mean/ of ellipsoid regions. If you interpret the---   shades as uncertain physical measurements with normal distribution,---   it gives the maximum-likelyhood result for multiple measurements of the---   same quantity.-mixShade's :: ∀ y . (WithField ℝ Manifold y, SimpleSpace (Needle y))-                 => NonEmpty (Shade' y) -> Maybe (Shade' y)-mixShade's = ms pseudoAffineWitness dualSpaceWitness- where ms :: PseudoAffineWitness y -> DualNeedleWitness y-                  -> NonEmpty (Shade' y) -> Maybe (Shade' y)-       ms (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness-                 (Shade' c₀ (Norm e₁):|shs) = sequenceA ciso >> pure mixed-        where ciso = [ci.-~.c₀ | Shade' ci shi <- shs]-              cis = [v | Just v <- ciso]-              σe = arr . sumV $ e₁ : (applyNorm . _shade'Narrowness<$>shs)-              cc = σe \$ sumV [ei $ ci | ci <- cis-                                       | Shade' _ (Norm ei) <- shs]-              mixed = Shade' (c₀+^cc) $ densifyNorm ( mconcat-                             [ Norm $ ei ^/ (1+(normSq ni $ ci^-^cc))-                             | ni@(Norm ei) <- Norm e₁ : (_shade'Narrowness<$>shs)-                             | ci <- zeroV : cis-                             ] )-              Tagged (+^) = translateP :: Tagged y (Interior y->Needle y->Interior y)-  -- cc should minimise the quadratic form-  -- β(cc) = ∑ᵢ ⟨cc−cᵢ|eᵢ|cc−cᵢ⟩-  -- = ⟨cc|e₁|cc⟩ + ∑ᵢ₌₁… ⟨cc−c₂|e₂|cc−c₂⟩-  -- = ⟨cc|e₁|cc⟩ + ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ − 2⋅⟨cᵢ|eᵢ|cc⟩ + ⟨cᵢ|eᵢ|cᵢ⟩ )-  -- It is thus-  -- β(cc + δ⋅v) − β cc-  -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩-  --     + ∑ᵢ₌₁…( ⟨cc + δ⋅v|eᵢ|cc + δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|cc + δ⋅v⟩ + ⟨cᵢ|eᵢ|cᵢ⟩ )-  --     − ⟨cc|e₁|cc⟩-  --     − ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ + 2⋅⟨cᵢ|eᵢ|cc⟩ − ⟨cᵢ|eᵢ|cᵢ⟩ )-  -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩-  --     + ∑ᵢ₌₁…( ⟨cc + δ⋅v|eᵢ|cc + δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|δ⋅v⟩ )-  --     − ⟨cc|e₁|cc⟩-  --     − ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ )-  -- = 2⋅⟨δ⋅v|e₁|cc⟩ + ⟨δ⋅v|e₁|δ⋅v⟩-  --     + ∑ᵢ₌₁…( 2⋅⟨δ⋅v|eᵢ|cc⟩ + ⟨δ⋅v|eᵢ|δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|δ⋅v⟩ )-  -- = 2⋅⟨δ⋅v|∑ᵢeᵢ|cc⟩ − 2⋅∑ᵢ₌₁… ⟨cᵢ|eᵢ|δ⋅v⟩ + 𝓞(δ²)-  -- This should vanish for all v, which is fulfilled by-  -- (∑ᵢeᵢ)|cc⟩ = ∑ᵢ₌₁… eᵢ|cᵢ⟩.---- | Evaluate the shade as a quadratic form; essentially--- @--- minusLogOcclusion sh x = x <.>^ (sh^.shadeExpanse $ x - sh^.shadeCtr)--- @--- where 'shadeExpanse' gives a metric (matrix) that characterises the--- width of the shade.-minusLogOcclusion' :: ∀ x s . ( PseudoAffine x, LinearSpace (Needle x)-                              , s ~ (Scalar (Needle x)), RealDimension s )-              => Shade' x -> x -> s-minusLogOcclusion' (Shade' p₀ δinv)-        = occ (pseudoAffineWitness :: PseudoAffineWitness x)-              (dualSpaceWitness :: DualNeedleWitness x)- where occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness-           p = case toInterior p >>= (.-~.p₀) of-         (Just vd) | mSq <- normSq δinv vd-                   , mSq == mSq  -- avoid NaN-                   -> mSq-         _         -> 1/0-minusLogOcclusion :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x)-                             , s ~ (Scalar (Needle x)), RealDimension s )-              => Shade x -> x -> s-minusLogOcclusion (Shade p₀ δ)-        = occ (pseudoAffineWitness :: PseudoAffineWitness x)-              (dualSpaceWitness :: DualNeedleWitness x)- where occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness-            = \p -> case toInterior p >>= (.-~.p₀) of-         (Just vd) | mSq <- normSq δinv vd-                   , mSq == mSq  -- avoid NaN-                   -> mSq-         _         -> 1/0-        where δinv = dualNorm δ-----rangeOnGeodesic :: ∀ i m . -      ( WithField ℝ PseudoAffine m, Geodesic m, SimpleSpace (Needle m)-      , WithField ℝ IntervalLike i, SimpleSpace (Needle i) )-                     => m -> m -> Maybe (Shade i -> Shade m)-rangeOnGeodesic = case ( semimanifoldWitness :: SemimanifoldWitness i-                       , dualSpaceWitness :: DualNeedleWitness i-                       , dualSpaceWitness :: DualNeedleWitness m ) of- (SemimanifoldWitness _, DualSpaceWitness, DualSpaceWitness) ->-  \p₀ p₁ -> (`fmap`(geodesicBetween p₀ p₁))-    $ \interp -> \(Shade t₀ et)-                -> case pointsShades-                         . mapMaybe (toInterior-                               . interp . (toClosedInterval :: i -> D¹))-                         $ fromInterior <$> t₀ : [ t₀+^v-                                                 | v<-normSpanningSystem et ] of-             [sh] -> sh-             _ -> case pointsShades $ mapMaybe (toInterior . interp . D¹)-                        [-0.999, 0.999] of-                [sh] -> sh- where Tagged (+^) = translateP :: Tagged i (Interior i->Needle i->Interior i)---- -- | 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 }@@ -712,6 +278,8 @@         = 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@@ -976,299 +544,8 @@                                  coerceShade  --- | Class of manifolds which can use 'Shade'' as a basic set type.---   This is easily possible for vector spaces with the default implementations.-class (WithField ℝ PseudoAffine y, SimpleSpace (Needle y)) => Refinable y where-  -- | @a `subShade'` b ≡ True@ means @a@ is fully contained in @b@, i.e. from-  --   @'minusLogOcclusion'' a p < 1@ follows also @minusLogOcclusion' b p < 1@.-  subShade' :: Shade' y -> Shade' y -> Bool-  subShade' (Shade' ac ae) (Shade' tc te)-        = case pseudoAffineWitness :: PseudoAffineWitness y of-   PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)-    | Just v <- tc.-~.ac-    , v² <- normSq te v-    , v² <= 1-     -> all (\(y',μ) -> case μ of-            Nothing -> True  -- 'te' has infinite extension in this direction-            Just ξ-              | ξ<1 -> False -- 'ae' would be vaster than 'te' in this direction-              | ω <- abs $ y'<.>^v-                    -> (ω + 1/ξ)^2 <= 1 - v² + ω^2-                 -- See @images/constructions/subellipse-check-heuristic.svg@-         ) $ sharedSeminormSpanningSystem te ae-   _ -> False-  -  -- | Intersection between two shades.-  refineShade' :: Shade' y -> Shade' y -> Maybe (Shade' y)-  refineShade' (Shade' c₀ (Norm e₁)) (Shade' c₀₂ (Norm e₂))-      = case ( dualSpaceWitness :: DualNeedleWitness y-             , pseudoAffineWitness :: PseudoAffineWitness y ) of-          (DualSpaceWitness, PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))-               -> do-           c₂ <- c₀₂.-~.c₀-           let σe = arr $ e₁^+^e₂-               e₁c₂ = e₁ $ c₂-               e₂c₂ = e₂ $ c₂-               cc = σe \$ e₂c₂-               cc₂ = cc ^-^ c₂-               e₁cc = e₁ $ cc-               e₂cc = e₂ $ cc-               α = 2 + e₂c₂<.>^cc₂-           guard (α > 0)-           let ee = σe ^/ α-               c₂e₁c₂ = e₁c₂<.>^c₂-               c₂e₂c₂ = e₂c₂<.>^c₂-               c₂eec₂ = (c₂e₁c₂ + c₂e₂c₂) / α-           return $ case middle . sort-                $ quadraticEqnSol c₂e₁c₂-                                  (2 * (e₁cc<.>^c₂))-                                  (e₁cc<.>^cc - 1)-                ++quadraticEqnSol c₂e₂c₂-                                  (2 * (e₂cc<.>^c₂ - c₂e₂c₂))-                                  (e₂cc<.>^cc - 2 * (e₂c₂<.>^cc) + c₂e₂c₂ - 1) of-            [γ₁,γ₂] | abs (γ₁+γ₂) < 2 -> let-               cc' = cc ^+^ ((γ₁+γ₂)/2)*^c₂-               rγ = abs (γ₁ - γ₂) / 2-               η = if rγ * c₂eec₂ /= 0 && 1 - rγ^2 * c₂eec₂ > 0-                   then sqrt (1 - rγ^2 * c₂eec₂) / (rγ * c₂eec₂)-                   else 0-             in Shade' (c₀.+~^cc')-                       (Norm (arr ee) <> spanNorm [ee $ c₂^*η])-            _ -> Shade' (c₀.+~^cc) (Norm $ arr ee)-   where quadraticEqnSol a b c-             | a == 0, b /= 0       = [-c/b]-             | a /= 0 && disc == 0  = [- b / (2*a)]-             | a /= 0 && disc > 0   = [ (σ * sqrt disc - b) / (2*a)-                                      | σ <- [-1, 1] ]-             | otherwise            = []-          where disc = b^2 - 4*a*c-         middle (_:x:y:_) = [x,y]-         middle l = l-  -- ⟨x−c₁|e₁|x−c₁⟩ < 1  ∧  ⟨x−c₂|e₂|x−c₂⟩ < 1-  -- We search (cc,ee) such that this implies-  -- ⟨x−cc|ee|x−cc⟩ < 1.-  -- Let WLOG c₁ = 0, so-  -- ⟨x|e₁|x⟩ < 1.-  -- cc should minimise the quadratic form-  -- β(cc) = ⟨cc−c₁|e₁|cc−c₁⟩ + ⟨cc−c₂|e₂|cc−c₂⟩-  -- = ⟨cc|e₁|cc⟩ + ⟨cc−c₂|e₂|cc−c₂⟩-  -- = ⟨cc|e₁|cc⟩ + ⟨cc|e₂|cc⟩ − 2⋅⟨c₂|e₂|cc⟩ + ⟨c₂|e₂|c₂⟩-  -- It is thus-  -- β(cc + δ⋅v) − β cc-  -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + ⟨cc + δ⋅v|e₂|cc + δ⋅v⟩ − 2⋅⟨c₂|e₂|cc + δ⋅v⟩ + ⟨c₂|e₂|c₂⟩-  --     − ⟨cc|e₁|cc⟩ − ⟨cc|e₂|cc⟩ + 2⋅⟨c₂|e₂|cc⟩ − ⟨c₂|e₂|c₂⟩-  -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + ⟨cc + δ⋅v|e₂|cc + δ⋅v⟩ − 2⋅⟨c₂|e₂|δ⋅v⟩-  --     − ⟨cc|e₁|cc⟩ − ⟨cc|e₂|cc⟩-  -- = 2⋅⟨δ⋅v|e₁|cc⟩ + ⟨δ⋅v|e₁|δ⋅v⟩ + 2⋅⟨δ⋅v|e₂|cc⟩ + ⟨δ⋅v|e₂|δ⋅v⟩ − 2⋅⟨c₂|e₂|δ⋅v⟩-  -- = 2⋅δ⋅⟨v|e₁+e₂|cc⟩ − 2⋅δ⋅⟨v|e₂|c₂⟩ + 𝓞(δ²)-  -- This should vanish for all v, which is fulfilled by-  -- (e₁+e₂)|cc⟩ = e₂|c₂⟩.-  -- -  -- If we now choose-  -- ee = (e₁+e₂) / α-  -- then-  -- ⟨x−cc|ee|x−cc⟩ ⋅ α-  --  = ⟨x−cc|ee|x⟩ ⋅ α − ⟨x−cc|ee|cc⟩ ⋅ α-  --  = ⟨x|ee|x−cc⟩ ⋅ α − ⟨x−cc|e₂|c₂⟩-  --  = ⟨x|ee|x⟩ ⋅ α − ⟨x|ee|cc⟩ ⋅ α − ⟨x−cc|e₂|c₂⟩-  --  = ⟨x|e₁+e₂|x⟩ − ⟨x|e₂|c₂⟩ − ⟨x−cc|e₂|c₂⟩-  --  = ⟨x|e₁|x⟩ + ⟨x|e₂|x⟩ − ⟨x|e₂|c₂⟩ − ⟨x−cc|e₂|c₂⟩-  --  < 1 + ⟨x|e₂|x−c₂⟩ − ⟨x−cc|e₂|c₂⟩-  --  = 1 + ⟨x−c₂|e₂|x−c₂⟩ + ⟨c₂|e₂|x−c₂⟩ − ⟨x−cc|e₂|c₂⟩-  --  < 2 + ⟨x−c₂−x+cc|e₂|c₂⟩-  --  = 2 + ⟨cc−c₂|e₂|c₂⟩-  -- Really we want-  -- ⟨x−cc|ee|x−cc⟩ ⋅ α < α-  -- So choose α = 2 + ⟨cc−c₂|e₂|c₂⟩.-  -- -  -- The ellipsoid "cc±√ee" captures perfectly the intersection-  -- of the boundary of the shades, but it tends to significantly-  -- overshoot the interior intersection in perpendicular direction,-  -- i.e. in direction of c₂−c₁. E.g.-  -- https://github.com/leftaroundabout/manifolds/blob/bc0460b9/manifolds/images/examples/ShadeCombinations/EllipseIntersections.png-  -- 1. Really, the relevant points are those where either of the-  --    intersector badnesses becomes 1. The intersection shade should-  --    be centered between those points. We perform according corrections,-  --    but only in c₂ direction, so this can be handled efficiently-  --    as a 1D quadratic equation.-  --    Consider-  --       dⱼ c := ⟨c−cⱼ|eⱼ|c−cⱼ⟩ =! 1-  --       dⱼ (cc + γ⋅c₂)-  --           = ⟨cc+γ⋅c₂−cⱼ|eⱼ|cc+γ⋅c₂−cⱼ⟩-  --           = ⟨cc−cⱼ|eⱼ|cc−cⱼ⟩ + 2⋅γ⋅⟨c₂|eⱼ|cc−cⱼ⟩ + γ²⋅⟨c₂|eⱼ|c₂⟩-  --           =! 1-  --    So-  --    γⱼ = (- b ± √(b²−4⋅a⋅c)) / 2⋅a-  --     where a = ⟨c₂|eⱼ|c₂⟩-  --           b = 2 ⋅ (⟨c₂|eⱼ|cc⟩ − ⟨c₂|eⱼ|cⱼ⟩)-  --           c = ⟨cc|eⱼ|cc⟩ − 2⋅⟨cc|eⱼ|cⱼ⟩ + ⟨cⱼ|eⱼ|cⱼ⟩ − 1-  --    The ± sign should be chosen to get the smaller |γ| (otherwise-  --    we end up on the wrong side of the shade), i.e.-  --    γⱼ = (sgn bⱼ ⋅ √(bⱼ²−4⋅aⱼ⋅cⱼ) − bⱼ) / 2⋅aⱼ-  -- 2. Trim the result in that direction to the actual-  --    thickness of the lens-shaped intersection: we want-  --    ⟨rγ⋅c₂|ee'|rγ⋅c₂⟩ = 1-  --    for a squeezed version of ee,-  --    ee' = ee + ee|η⋅c₂⟩⟨η⋅c₂|ee-  --    ee' = ee + η² ⋅ ee|c₂⟩⟨c₂|ee-  --    ⟨rγ⋅c₂|ee'|rγ⋅c₂⟩-  --        = rγ² ⋅ (⟨c₂|ee|c₂⟩ + η² ⋅ ⟨c₂|ee|c₂⟩²)-  --        = rγ² ⋅ ⟨c₂|ee|c₂⟩ + η² ⋅ rγ² ⋅ ⟨c₂|ee|c₂⟩²-  --    η² = (1 − rγ²⋅⟨c₂|ee|c₂⟩) / (rγ² ⋅ ⟨c₂|ee|c₂⟩²)-  --    η = √(1 − rγ²⋅⟨c₂|ee|c₂⟩) / (rγ ⋅ ⟨c₂|ee|c₂⟩)-  --    With ⟨c₂|ee|c₂⟩ = (⟨c₂|e₁|c₂⟩ + ⟨c₂|e₂|c₂⟩)/α. -  -  -- | If @p@ is in @a@ (red) and @δ@ is in @b@ (green),-  --   then @p.+~^δ@ is in @convolveShade' a b@ (blue).-  -- ---   Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/ShadeCombinations.ipynb#shadeConvolutions--- --- <<images/examples/ShadeCombinations/2Dconvolution-skewed.png>>-  convolveMetric :: Hask.Functor p => p y -> Metric y -> Metric y -> Metric y-  convolveMetric _ ey eδ = spanNorm [ f ^* ζ crl-                                    | (f,crl) <- eδsp ]-   where eδsp = sharedSeminormSpanningSystem ey eδ-         ζ = case filter (>0) . catMaybes $ snd<$>eδsp of-            [] -> const 0-            nzrelap-               -> let cre₁ = 1/minimum nzrelap-                      cre₂ =  maximum nzrelap-                      edgeFactor = sqrt ( (1 + cre₁)^2 + (1 + cre₂)^2 )-                                / (sqrt (1 + cre₁^2) + sqrt (1 + cre₂^2))-                  in \case-                        Nothing -> 0-                        Just 0  -> 0-                        Just sq -> edgeFactor / (recip sq + 1)-  -  convolveShade' :: Shade' y -> Shade' (Needle y) -> Shade' y-  convolveShade' = defaultConvolveShade'-  -defaultConvolveShade' :: ∀ y . Refinable y => Shade' y -> Shade' (Needle y) -> Shade' y-defaultConvolveShade' = case (pseudoAffineWitness :: PseudoAffineWitness y) of-  PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)-    -> \(Shade' y₀ ey) (Shade' δ₀ eδ) -> Shade' (y₀.+~^δ₀)-                                          $ convolveMetric ([]::[y]) ey eδ -instance Refinable ℝ where-  refineShade' (Shade' cl el) (Shade' cr er)-         = case (normSq el 1, normSq er 1) of-             (0, _) -> return $ Shade' cr er-             (_, 0) -> return $ Shade' cl el-             (ql,qr) | ql>0, qr>0-                    -> let [rl,rr] = sqrt . recip <$> [ql,qr]-                           b = maximum $ zipWith (-) [cl,cr] [rl,rr]-                           t = minimum $ zipWith (+) [cl,cr] [rl,rr]-                       in guard (b<t) >>-                           let cm = (b+t)/2-                               rm = (t-b)/2-                           in return $ Shade' cm (spanNorm [recip rm])---   convolveShade' (Shade' y₀ ey) (Shade' δ₀ eδ)---          = case (metricSq ey 1, metricSq eδ 1) of---              (wy,wδ) | wy>0, wδ>0---                  -> Shade' (y₀.+~^δ₀)---                            ( projector . recip---                                   $ recip (sqrt wy) + recip (sqrt wδ) )---              (_ , _) -> Shade' y₀ zeroV--instance ( Refinable a, Interior a ~ a, Refinable b, Interior b ~ b-         , Scalar (DualVector (DualVector (Needle b)))-                      ~ Scalar (DualVector (DualVector (Needle a))) )-    => Refinable (a,b)-  -instance Refinable ℝ⁰-instance Refinable ℝ¹-instance Refinable ℝ²-instance Refinable ℝ³-instance Refinable ℝ⁴-                            -instance ( SimpleSpace a, SimpleSpace b-         , Scalar a ~ ℝ, Scalar b ~ ℝ-         , Scalar (DualVector a) ~ ℝ, Scalar (DualVector b) ~ ℝ-         , Scalar (DualVector (DualVector a)) ~ ℝ, Scalar (DualVector (DualVector b)) ~ ℝ )-            => Refinable (LinearMap ℝ a b)--intersectShade's :: ∀ y . Refinable y => NonEmpty (Shade' y) -> Maybe (Shade' y)-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----propagateDEqnSolution_loc :: ∀ x y . ( WithField ℝ Manifold x-                                     , Refinable y, Geodesic (Interior y)-                                     , SimpleSpace (Needle x) )-           => 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)- where pdesl DualSpaceWitness DualSpaceWitness BoundarylessWitness-             (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))-          | Nothing <- jacobian  = Nothing-          | otherwise            = pure result-         where jacobian = f shxy ^. predictDerivatives-               Just (Shade' j₀ jExpa) = jacobian--               mx = propPlan^.sourcePosition .+~^ propPlan^.targetPosOffset ^/ 2-               Just my = middleBetween (propPlan^.sourceData.shadeCtr)-                                       (propPlan^.targetAPrioriData.shadeCtr)-               shxy = coverAllAround (mx, my)-                                     [ (δx ^-^ propPlan^.targetPosOffset ^/ 2, py ^+^ v)-                                     | (δx,ney) <- (zeroV, propPlan^.sourceData)-                                                  : (propPlan^.relatedData)-                                     , let Just py = ney^.shadeCtr .-~. my-                                     , v <- normSpanningSystem' (ney^.shadeNarrowness)-                                     ]-               (Shade _ expax' :: Shade x)-                    = coverAllAround (propPlan^.sourcePosition)-                                     [δx | (δx,_) <- propPlan^.relatedData]-               expax = dualNorm expax'-               result :: Shade' y-               result = convolveShade'-                        (propPlan^.sourceData)-                        (Shade' δyb $ applyLinMapNorm jExpa dx)-                where δyb = j₀ $ δx-               δx = propPlan^.targetPosOffset-               dx = δx'^/(δx'<.>^δx)-                where δx' = expax<$|δx--applyLinMapNorm :: ∀ x y . (LSpace x, LSpace y, Scalar x ~ Scalar y)-           => Norm (x+>y) -> DualVector x -> Norm y-applyLinMapNorm = case dualSpaceWitness :: DualSpaceWitness y of-  DualSpaceWitness -> \n dx -> transformNorm (arr $ LinearFunction (dx-+|>)) n--ignoreDirectionalDependence :: ∀ x y . (LSpace x, LSpace y, Scalar x ~ Scalar y)-           => (x, DualVector x) -> Norm (x+>y) -> Norm (x+>y)-ignoreDirectionalDependence = case dualSpaceWitness :: DualSpaceWitness y of-  DualSpaceWitness -> \(v,v') -> transformNorm . arr . LinearFunction $-         \j -> j . arr (LinearFunction $ \x -> x ^-^ v^*(v'<.>^x))- type Twig x = (Int, ShadeTree x) type TwigEnviron x = [Twig x] @@ -1693,64 +970,6 @@   --- | 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-      = WithAny { _untopological :: y-                , _topological :: !x  }- deriving (Hask.Functor, Show, Generic)--instance (NFData x, NFData y) => NFData (WithAny x y)--instance ∀ x y . (Semimanifold x) => Semimanifold (x`WithAny`y) where-  type Needle (WithAny x y) = Needle x-  type Interior (WithAny x y) = Interior x `WithAny` y-  WithAny y x .+~^ δx = WithAny y $ x.+~^δx-  fromInterior (WithAny y x) = WithAny y $ fromInterior x-  toInterior (WithAny y x) = fmap (WithAny y) $ toInterior x-  translateP = tpWD-   where tpWD :: ∀ x y . Semimanifold x => Tagged (WithAny x y)-                            (Interior x`WithAny`y -> Needle x -> Interior x`WithAny`y)-         tpWD = Tagged `id` \(WithAny y x) δx -> WithAny y $ tpx x δx-          where Tagged tpx = translateP :: Tagged x (Interior x -> Needle x -> Interior x)-  semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of-      SemimanifoldWitness BoundarylessWitness -> SemimanifoldWitness BoundarylessWitness-            -instance (PseudoAffine x) => PseudoAffine (x`WithAny`y) where-  WithAny _ x .-~. WithAny _ ξ = x.-~.ξ-  pseudoAffineWitness = case pseudoAffineWitness :: PseudoAffineWitness x of-      PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)-       -> PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)--instance (AffineSpace x) => AffineSpace (x`WithAny`y) where-  type Diff (WithAny x y) = Diff x-  WithAny _ x .-. WithAny _ ξ = x.-.ξ-  WithAny y x .+^ δx = WithAny y $ x.+^δx --instance (VectorSpace x, Monoid y) => VectorSpace (x`WithAny`y) where-  type Scalar (WithAny x y) = Scalar x-  μ *^ WithAny y x = WithAny y $ μ*^x --instance (AdditiveGroup x, Monoid y) => AdditiveGroup (x`WithAny`y) where-  zeroV = WithAny mempty zeroV-  negateV (WithAny y x) = WithAny y $ negateV x-  WithAny y x ^+^ WithAny υ ξ = WithAny (mappend y υ) (x^+^ξ)--instance (AdditiveGroup x) => Hask.Applicative (WithAny x) where-  pure x = WithAny x zeroV-  WithAny f x <*> WithAny t ξ = WithAny (f t) (x^+^ξ)-  -instance (AdditiveGroup x) => Hask.Monad (WithAny x) where-  return x = WithAny x zeroV-  WithAny y x >>= f = WithAny r $ x^+^q-   where WithAny r q = f y--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 (Shade (WithAny _ b) e) = Shade b e- constShaded :: y -> ShadeTree x -> x`Shaded`y constShaded y = unsafeFmapTree (WithAny y<$>) id (shadeWithAny y) @@ -1789,7 +1008,7 @@ -- | 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, WithField ℝ LinearManifold y+stiWithDensity :: ∀ x y . ( WithField ℝ PseudoAffine x, LinearSpace y, Scalar y ~ ℝ                           , SimpleSpace (Needle x) )          => x`Shaded`y -> x -> Cℝay y stiWithDensity (PlainLeaves lvs)@@ -1832,7 +1051,7 @@                  -> ( xloc, ( (yloc, recip $ shd|$|(0,1))                             , dependence (dualNorm shd) ) ) -smoothInterpolate :: ∀ x y . ( WithField ℝ Manifold x, WithField ℝ LinearManifold y+smoothInterpolate :: ∀ x y . ( WithField ℝ Manifold x, LinearSpace y, Scalar y ~ ℝ                              , SimpleSpace (Needle x) )              => NonEmpty (x,y) -> x -> y smoothInterpolate = si boundarylessWitness@@ -1911,10 +1130,4 @@   ---extractJust :: (a->Maybe b) -> [a] -> (Maybe b, [a])-extractJust f [] = (Nothing,[])-extractJust f (x:xs) | Just r <- f x  = (Just r, xs)-                     | otherwise      = second (x:) $ extractJust f xs 

@@ -18,6 +18,7 @@ {-# LANGUAGE FunctionalDependencies   #-} {-# LANGUAGE FlexibleContexts         #-} {-# LANGUAGE LiberalTypeSynonyms      #-}+{-# LANGUAGE StandaloneDeriving       #-} {-# LANGUAGE GADTs                    #-} {-# LANGUAGE RankNTypes               #-} {-# LANGUAGE TupleSections            #-}@@ -58,6 +59,8 @@         , fathomCutDistance, sideOfCut, cutPosBetween         -- * Linear mappings         , LinearMap, LocalLinear+        -- * Misc+        , StiefelScalar    ) where  @@ -81,6 +84,7 @@ import Math.LinearMap.Category  import qualified Prelude+import qualified Data.Traversable as Hask  import Control.Category.Constrained.Prelude hiding ((^)) import Control.Arrow.Constrained@@ -89,6 +93,8 @@  import Data.Type.Coercion +type StiefelScalar s = (RealFloat s, UArr.Unbox s)+ #define deriveAffine(c,t)                \ instance (c) => Semimanifold (t) where {  \   type Needle (t) = Diff (t);              \@@ -101,6 +107,7 @@   newtype Stiefel1Needle v = Stiefel1Needle { getStiefel1Tangent :: UArr.Vector (Scalar v) }+deriving instance (Eq (Scalar v), UArr.Unbox (Scalar v)) => Eq (Stiefel1Needle v) newtype Stiefel1Basis v = Stiefel1Basis { getStiefel1Basis :: Int } s1bTrie :: ∀ v b. FiniteFreeSpace v => (Stiefel1Basis v->b) -> Stiefel1Basis v:->:b s1bTrie = \f -> St1BTrie $ fmap (f . Stiefel1Basis) allIs@@ -163,11 +170,12 @@  deriveAffine((FiniteFreeSpace v, UArr.Unbox (Scalar v)), Stiefel1Needle v) -instance ∀ v . (LSpace v, FiniteFreeSpace v, UArr.Unbox (Scalar v))+instance ∀ v . (LSpace v, FiniteFreeSpace v, Eq (Scalar v), UArr.Unbox (Scalar v))               => TensorSpace (Stiefel1Needle v) where   type TensorProduct (Stiefel1Needle v) w = Array w   scalarSpaceWitness = case scalarSpaceWitness :: ScalarSpaceWitness v of          ScalarSpaceWitness -> ScalarSpaceWitness+  linearManifoldWitness = LinearManifoldWitness BoundarylessWitness   zeroTensor = Tensor $ Arr.replicate (freeDimension ([]::[v]) - 1) zeroV   toFlatTensor = LinearFunction $ Tensor . Arr.convert . getStiefel1Tangent   fromFlatTensor = LinearFunction $ Stiefel1Needle . Arr.convert . getTensorProduct@@ -184,6 +192,7 @@   fzipTensorWith = bilinearFunction $ \f (Tensor a, Tensor b)                      -> Tensor $ Arr.zipWith (curry $ arr f) a b   coerceFmapTensorProduct _ Coercion = Coercion+  wellDefinedTensor (Tensor a) = Tensor <$> Hask.traverse wellDefinedVector a  asTensor :: Coercion (LinearMap s a b) (Tensor s (DualVector a) b) asTensor = Coercion@@ -194,7 +203,7 @@             => LinearMap s a b -> a -> b (+$>) = getLinearFunction . getLinearFunction applyLinear   -instance ∀ v . (LSpace v, FiniteFreeSpace v, UArr.Unbox (Scalar v))+instance ∀ v . (LSpace v, FiniteFreeSpace v, Eq (Scalar v), UArr.Unbox (Scalar v))               => LinearSpace (Stiefel1Needle v) where   type DualVector (Stiefel1Needle v) = Stiefel1Needle v   linearId = LinearMap . Arr.generate d $ \i -> Stiefel1Needle . Arr.generate d $@@ -228,9 +237,9 @@   composeLinear = bilinearFunction $ \f (LinearMap g)                      -> LinearMap $ Arr.map (getLinearFunction applyLinear f$) g -instance ∀ k v .-   ( WithField k LinearManifold v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v)-   , RealFloat k, UArr.Unbox k ) => Semimanifold (Stiefel1 v) where +instance ∀ v .+   ( LinearSpace v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v)+   , StiefelScalar (Scalar v) ) => Semimanifold (Stiefel1 v) where   type Needle (Stiefel1 v) = Stiefel1Needle v   fromInterior = id   toInterior = pure@@ -261,9 +270,9 @@                 insi ti v = Arr.generate d $ \i -> if | i<im      -> v Arr.! i                                                       | i>im      -> v Arr.! (i-1)                                                        | otherwise -> ti-instance ∀ k v .-   ( WithField k LinearManifold v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v)-   , RealFloat k, UArr.Unbox k ) => PseudoAffine (Stiefel1 v) where +instance ∀ v .+   ( LinearSpace v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v)+   , StiefelScalar (Scalar v) ) => PseudoAffine (Stiefel1 v) where   (.-~.) = dpst dualSpaceWitness    where dpst :: DualSpaceWitness v -> Stiefel1 v -> Stiefel1 v -> Maybe (Stiefel1Needle v)          dpst DualSpaceWitness (Stiefel1 s) (Stiefel1 t)

@@ -45,22 +45,27 @@               -- ** Local environments             , localFocusWeb               -- * Uncertain functions-            , differentiateUncertainWebFunction+            , differentiateUncertainWebFunction, differentiate²UncertainWebFunction               -- * Differential equations               -- ** Fixed resolution-            , filterDEqnSolution_static, iterateFilterDEqn_static+            , iterateFilterDEqn_static               -- ** Automatic resolution             , filterDEqnSolutions_adaptive, iterateFilterDEqn_adaptive               -- ** Configuration             , InconsistencyStrategy(..)+            , InformationMergeStrategy(..)+            , naïve, inconsistencyAware, indicateInconsistencies+            , PropagationInconsistency(..)               -- * Misc-            , ConvexSet(..), ellipsoid, coerceWebDomain+            , ConvexSet(..), ellipsoid, ellipsoidSet, coerceWebDomain+            , rescanPDEOnWeb, rescanPDELocally, webOnions             ) 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 import qualified Data.Vector.Unboxed as UArr@@ -81,9 +86,14 @@ import Data.Manifold.Types import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine+import Data.Manifold.Shade import Data.Manifold.TreeCover import Data.SetLike.Intersection import Data.Manifold.Riemannian+import Data.Manifold.Atlas+import Data.Manifold.Function.Quadratic+import Data.Function.Affine+import Data.Embedding      import qualified Prelude as Hask hiding(foldl, sum, sequence) import qualified Control.Applicative as Hask@@ -92,6 +102,7 @@ import Data.STRef (newSTRef, modifySTRef, readSTRef) import Control.Monad.Trans.State import Control.Monad.Trans.List+import Control.Monad.Trans.Except import Data.Functor.Identity (Identity(..)) import qualified Data.Foldable       as Hask import Data.Foldable (all, toList)@@ -107,6 +118,7 @@ import Data.Traversable.Constrained (Traversable, traverse)  import Control.Comonad (Comonad(..))+import Control.Comonad.Cofree import Control.Lens ((&), (%~), (^.), (.~), (+~)) import Control.Lens.TH @@ -145,7 +157,18 @@           } 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†@@ -532,7 +555,24 @@                                 ]), n)                  ) asd' +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+        | Map.null next = []+        | otherwise  = ( 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 ]) +webOnions :: ∀ x y . WithField ℝ Manifold x+            => PointsWeb x y -> PointsWeb x [[(x,y)]]+webOnions = localFmapWeb (map (map $ _thisNodeCoord&&&_thisNodeData) . localOnion)+ nearestNeighbour :: (WithField ℝ Manifold x, SimpleSpace (Needle x))                       => PointsWeb x y -> x -> Maybe (x,y) nearestNeighbour (PointsWeb rsc asd) x = fmap lkBest $ positionIndex empty rsc x@@ -595,11 +635,62 @@              -> PointsWeb x (Shade' (LocalLinear x y)) differentiateUncertainWebFunction = localFmapWeb differentiateUncertainWebLocally -rescanPDELocally :: ∀ x y .-     ( WithField ℝ Manifold x, SimpleSpace (Needle x)-     , WithField ℝ Refinable y, SimpleSpace (Needle y) )-         => DifferentialEqn x y -> WebLocally x (Shade' y)-                                -> Maybe (Shade' y)+differentiate²UncertainWebLocally :: ∀ x y+   . ( WithField ℝ Manifold x, FlatSpace (Needle x)+     , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle 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 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+                          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)))++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+              , dualSpaceWitness :: DualSpaceWitness (Needle v)+              , dualSpaceWitness :: DualSpaceWitness (Needle y) ) of+   (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) )+         => 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 = case ( dualSpaceWitness :: DualNeedleWitness x                         , dualSpaceWitness :: DualNeedleWitness y                         , boundarylessWitness :: BoundarylessWitness x@@ -614,16 +705,20 @@                                      , v <- normSpanningSystem'                                               (ngb^.thisNodeData.shadeNarrowness)] of                         LocalDifferentialEqn _ rescan-                            -> rescan (differentiateUncertainWebLocally info)-                                      (info^.thisNodeData)+                            -> rescan (info^.thisNodeData)+                                      (differentiateUncertainWebLocally info)+                                      (differentiate²UncertainWebLocally info) -rescanPDEOnWeb :: ( WithField ℝ Manifold x, SimpleSpace (Needle x)-                  , WithField ℝ Refinable y, SimpleSpace (Needle y)+rescanPDEOnWeb :: ( WithField ℝ Manifold x, FlatSpace (Needle x)+                  , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y)                   , Hask.Applicative m )-                => InconsistencyStrategy m x (Shade' y)-                  -> DifferentialEqn x y -> PointsWeb x (Shade' y)-                                   -> m (PointsWeb x (Shade' y))-rescanPDEOnWeb strat = traverseWebWithStrategy strat . rescanPDELocally+                => 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"))  toGraph :: (WithField ℝ Manifold x, SimpleSpace (Needle x))               => PointsWeb x y -> (Graph, Vertex -> (x, y))@@ -644,12 +739,15 @@       -- ^ If @p@ is in all intersectors, it must also be in the hull.     , convexSetIntersectors :: [Shade' x]     }-deriving instance ( WithField ℝ Manifold x, SimpleSpace (Needle x)-                  , Show (Interior x), Show (Needle' x) ) => Show (ConvexSet x)+deriving instance LtdErrorShow x => Show (ConvexSet x)  ellipsoid :: Shade' x -> ConvexSet x ellipsoid s = ConvexSet s [s] +ellipsoidSet :: Embedding (->) (Maybe (Shade' x)) (ConvexSet x)+ellipsoidSet = Embedding (\case {Just s -> ConvexSet s [s]; Nothing -> EmptyConvex})+                         (\case {ConvexSet h _ -> Just h; EmptyConvex -> Nothing})+ intersectors :: ConvexSet x -> Maybe (NonEmpty (Shade' x)) intersectors (ConvexSet h []) = pure (h:|[]) intersectors (ConvexSet _ (i:sts)) = pure (i:|sts)@@ -684,7 +782,26 @@ dupHead (x:|xs) = x:|x:xs  +newtype InformationMergeStrategy n m y' y = InformationMergeStrategy+    { mergeInformation :: y -> n y' -> m y } +naïve :: (NonEmpty y -> y) -> InformationMergeStrategy [] Identity (x,y) y+naïve merge = InformationMergeStrategy (\o n -> Identity . merge $ o :| fmap snd n)++inconsistencyAware :: (NonEmpty y -> m y) -> InformationMergeStrategy [] m (x,y) y+inconsistencyAware merge = InformationMergeStrategy (\o n -> merge $ o :| fmap snd n)++indicateInconsistencies :: (NonEmpty υ -> Maybe υ)+         -> InformationMergeStrategy [] (Except (PropagationInconsistency x υ)) (x,υ) υ+indicateInconsistencies merge = InformationMergeStrategy+           (\o n -> case merge $ o :| fmap snd n of+               Just r  -> pure r+               Nothing -> throwE $ PropagationInconsistency n o )++maybeAlt :: Hask.Alternative f => Maybe a -> f a+maybeAlt (Just x) = pure x+maybeAlt Nothing = Hask.empty+ data InconsistencyStrategy m x y where     AbortOnInconsistency :: InconsistencyStrategy Maybe x y     IgnoreInconsistencies :: InconsistencyStrategy Identity x y@@ -692,76 +809,56 @@ deriving instance Hask.Functor (InconsistencyStrategy m x)  -iterateFilterDEqn_static :: ( WithField ℝ Manifold x, SimpleSpace (Needle x)-                            , Refinable y, Geodesic (Interior y)-                            , Hask.Applicative m )-       => InconsistencyStrategy m x (Shade' y) -> DifferentialEqn x y-                 -> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)]-iterateFilterDEqn_static strategy f-                           = map (fmap convexSetHull)-                           . itWhileJust strategy-                                (filterDEqnSolutions_static (ellipsoid<$>strategy) f)-                           . fmap (`ConvexSet`[])+iterateFilterDEqn_static :: ( WithField ℝ Manifold x, FlatSpace (Needle x)+                            , Refinable y, Geodesic y, FlatSpace (Needle y)+                            , WithField ℝ AffineManifold ð, Geodesic ð+                            , SimpleSpace (Needle ð)+                            , Hask.MonadPlus m )+       => InformationMergeStrategy [] m (x,Shade' y) iy+           -> Embedding (->) (Shade' y) iy+           -> DifferentialEqn x ð y+                 -> PointsWeb x (Shade' y) -> Cofree m (PointsWeb x (Shade' y))+iterateFilterDEqn_static strategy shading f+                           = fmap (fmap (shading >-$))+                           . coiter (filterDEqnSolutions_static strategy shading f)+                           . fmap (shading $->) -filterDEqnSolution_static :: ∀ x y m . ( WithField ℝ Manifold x, SimpleSpace (Needle x)-                                       , Refinable y, Geodesic (Interior y) )-       => InconsistencyStrategy m x (Shade' y) -> DifferentialEqn x y-            -> PointsWeb x (Shade' y) -> m (PointsWeb x (Shade' y))-filterDEqnSolution_static strat@AbortOnInconsistency f-    = case boundarylessWitness :: BoundarylessWitness x of-     BoundarylessWitness ->-        rescanPDEOnWeb strat f >=> webLocalInfo-           >>> Hask.traverse `id`\me -> case me^.nodeNeighbours of-                  []   -> return $ me^.thisNodeData-                  ngbs -> refineShade' (me^.thisNodeData)-                            =<< intersectShade's-                            =<< ( sequenceA $ NE.fromList-                                  [ propagateDEqnSolution_loc-                                       f (LocalDataPropPlan-                                             (ngbInfo^.thisNodeCoord)-                                             (negateV δx)-                                             (ngbInfo^.thisNodeData)-                                             (me^.thisNodeData)-                                             (fmap (second _thisNodeData . snd)-                                                       $ ngbInfo^.nodeNeighbours)-                                          )-                                  | (_, (δx, ngbInfo)) <- ngbs-                                  ] ) -filterDEqnSolutions_static :: ∀ x y m .-                              ( WithField ℝ Manifold x, SimpleSpace (Needle x)-                              , Refinable y, Geodesic (Interior y)-                              , Hask.Applicative m )-       => InconsistencyStrategy m x (ConvexSet y) -> DifferentialEqn x y-            -> PointsWeb x (ConvexSet y) -> m (PointsWeb x (ConvexSet y))-filterDEqnSolutions_static strategy f-       = webLocalInfo-           >>> fmap (id &&& rescanPDELocally f . fmap convexSetHull)-           >>> localFocusWeb >>> Hask.traverse `id`\((_,(me,updShy)), ngbs)-          -> let oldValue = me^.thisNodeData :: ConvexSet y+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 )+       => InformationMergeStrategy [] m  (x,Shade' y) iy -> Embedding (->) (Shade' y) 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 -> case ngbs of+              (Just shy, Just shð) -> case ngbs of                   []  -> pure oldValue                   _:_ | BoundarylessWitness <- (boundarylessWitness::BoundarylessWitness x)-                    -> handleInconsistency strategy oldValue-                          $ ( sequenceA $ NE.fromList-                                  [ sj >>= \ngbShy ->+                    -> maybeAlt+                          ( sequenceA [ fzip sj+                                >>= \ngbShyð -> (ngbInfo^.thisNodeCoord,)<$>                                      propagateDEqnSolution_loc                                        f (LocalDataPropPlan                                              (ngbInfo^.thisNodeCoord)                                              (negateV δx)-                                             ngbShy-                                             shy-                                             (fmap (second (convexSetHull . _thisNodeData)+                                             ngbShyð+                                             (shy, shð)+                                             (fmap (second ((shading>-$) . _thisNodeData)                                                     . snd) $ ngbInfo^.nodeNeighbours)                                           )                                   | (δx, (ngbInfo,sj)) <- ngbs                                   ] )-                            >>= intersectShade's-                            >>= pure . ((oldValue<>) . ellipsoid)-                            >>= \case EmptyConvex -> empty-                                      c           -> pure c-              _ -> handleInconsistency strategy oldValue empty+                            >>= mergeInformation strategy (shading$->shy)+              _ -> mergeInformation strategy oldValue empty+        )  handleInconsistency :: InconsistencyStrategy m x a -> a -> Maybe a -> m a handleInconsistency AbortOnInconsistency _ i = i@@ -790,23 +887,25 @@ oldAndNew' (_, l) = (False,) <$> l  -filterDEqnSolutions_adaptive :: ∀ x y badness m+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 )        => 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)) filterDEqnSolutions_adaptive mf strategy f badness' oldState-            = fmap recomputeJacobian $ filterGo boundarylessWitness-                                         =<< tryPreproc boundarylessWitness- where tryPreproc :: BoundarylessWitness x+            = fmap recomputeJacobian $ filterGo boundarylessWitness geodesicWitness+                                         =<< tryPreproc boundarylessWitness geodesicWitness+ where tryPreproc :: BoundarylessWitness x -> GeodesicWitness y                       -> m (PointsWeb x ( (WebLocally x (SolverNodeState x y)                                         , [(Shade' y, badness)]) ))-       tryPreproc BoundarylessWitness = traverse addPropagation $ webLocalInfo oldState+       tryPreproc BoundarylessWitness (GeodesicWitness _)+               = traverse addPropagation $ webLocalInfo oldState         where addPropagation wl                  | null neighbourInfo = pure (wl, [])                  | otherwise           = (wl,) . map (id&&&badness undefined)@@ -818,8 +917,9 @@                                            (neigh^.thisNodeCoord)                                            (negateV δx)                                            (convexSetHull $ neigh^.thisNodeData-                                                                  .solverNodeStatus)-                                           (thisShy)+                                                                  .solverNodeStatus+                                           , undefined)+                                           (thisShy, undefined)                                            [ second (convexSetHull                                                      . _solverNodeStatus . _thisNodeData) nn                                            | (_,nn)<-neigh^.nodeNeighbours ] )@@ -832,11 +932,12 @@        errTgtModulation = (1-) . (`mod'`1) . negate . sqrt $ fromIntegral totalAge        badness x = badness' x . (shadeNarrowness %~ (scaleNorm errTgtModulation))               -       filterGo :: BoundarylessWitness x+       filterGo :: BoundarylessWitness x -> GeodesicWitness y                    -> (PointsWeb x ( (WebLocally x (SolverNodeState x y)                                    , [(Shade' y, badness)]) ))                    -> m (PointsWeb x (SolverNodeState x y))-       filterGo BoundarylessWitness preproc'd   = fmap (smoothenWebTopology mf+       filterGo BoundarylessWitness (GeodesicWitness _) preproc'd+             = fmap (smoothenWebTopology mf                                      . fromTopWebNodes mf . concat . fmap retraceBonds                                         . Hask.toList . webLocalInfo . webLocalInfo)              $ Hask.traverse (uncurry localChange) preproc'd@@ -912,8 +1013,9 @@                                                       (n^.thisNodeCoord)                                                       (stepV ^-^ δx)                                                       (convexSetHull $-                                                        n^.thisNodeData.solverNodeStatus)-                                                      (aprioriInterpolate)+                                                        n^.thisNodeData.solverNodeStatus+                                                      , undefined)+                                                      (aprioriInterpolate, undefined)                                                       (second (convexSetHull                                                                ._solverNodeStatus                                                                ._thisNodeData)@@ -978,10 +1080,11 @@  iterateFilterDEqn_adaptive      :: ( WithField ℝ Manifold x, SimpleSpace (Needle x)-        , WithField ℝ AffineManifold y, Refinable y, Geodesic y, Hask.Monad m )+        , WithField ℝ AffineManifold y, Refinable y, Geodesic y, Hask.Monad m+        , WithField ℝ AffineManifold ð, Geodesic ð, SimpleSpace (Needle ð) )        => 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

binary file changed (absent → 41642 bytes)

@@ -1,5 +1,5 @@ Name:                manifolds-Version:             0.4.0.0+Version:             0.4.1.0 Category:            Math Synopsis:            Coordinate-free hypersurfaces Description:         Manifolds, a generalisation of the notion of &#x201c;smooth curves&#x201d; or surfaces,@@ -40,18 +40,20 @@  Library   Build-Depends:     base>=4.5 && < 6-                     , manifolds-core == 0.4.0.0+                     , manifolds-core == 0.4.1.0                      , transformers                      , vector-space>=0.8                      , free-vector-spaces>=0.1.1                      , linear                      , MemoTrie                      , vector-                     , linearmap-category > 0.3 && < 0.4+                     , linearmap-category >= 0.3.2 && < 0.4                      , containers                      , comonad+                     , free                      , semigroups                      , void+                     , number-show >= 0.1 && < 0.2                      , tagged                      , deepseq                      , lens@@ -70,10 +72,12 @@   Exposed-modules:   Data.Manifold                      Data.Manifold.PseudoAffine                      Data.Manifold.TreeCover+                     Data.Manifold.Shade                      Data.Manifold.Web                      Data.Manifold.DifferentialEquation                      Data.SimplicialComplex                      Data.Function.Differentiable+                     Data.Function.Affine                      Data.Manifold.Types                      Data.Manifold.Types.Stiefel                      Data.Manifold.Griddable@@ -85,8 +89,8 @@                    Data.Manifold.Cone                    Data.CoNat                    Data.Embedding+                   Data.Manifold.Function.Quadratic                    Data.Function.Differentiable.Data-                   Data.Function.Affine                    Control.Monad.Trans.OuterMaybe                    Util.Associate                    Util.LtdShow