manifolds 0.2.0.1 → 0.2.2.0

raw patch · 14 files changed

+723/−160 lines, 14 filesbinary-addedPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Data/LinearMap/Category.hs +6/−0
• Data/LinearMap/HerMetric.hs +160/−22
• Data/Manifold/PseudoAffine.hs +57/−3
• Data/Manifold/TreeCover.hs +232/−132
• Data/Manifold/Types.hs +1/−1
• Data/Manifold/Web.hs +219/−0
• Data/SetLike/Intersection.hs +37/−0
• Data/VectorSpace/FiniteDimensional.hs +6/−0
• images/examples/ShadeCombinations/2Dconvolution-skewed.png binary
• images/examples/cartesiandisk-2d-PointsWeb.png binary
• images/examples/normaldistrib-2d-PointsWeb.png binary
• images/examples/normaldistrib-2d-ShadeTree.png binary
• images/examples/simple-2d-PointsWeb.png binary
• manifolds.cabal +5/−2

@@ -70,6 +70,12 @@ identMat = DenseLinear $ HMat.ident n  where (Tagged n) = dimension :: Tagged v Int +-- | Coerce the matrix representations of two linear mappings; the result makes+--   sense iff the spaces are canonically isomorphic (certainly if they+--   are good instances of 'Data.Manifold.PseudoAffine.LocallyCoercible').+unsafeCoerceLinear :: Linear s a b -> Linear s c d+unsafeCoerceLinear (DenseLinear m) = DenseLinear m+ convertLinear :: ∀ v w s . ( FiniteDimensional v, FiniteDimensional w                            , Scalar v ~ s, Scalar w ~ s )                    => Isomorphism (->) (v:-*w) (Linear s v w)

@@ -21,7 +21,7 @@   , toDualWith, fromDualWith   , metricSq, metricSq', metric, metric', metrics, metrics'   -- * Defining metrics-  , projector, projector'+  , projector, projector', projectors, projector's   , euclideanMetric'   -- * Metrics induce inner products   , spanHilbertSubspace@@ -35,8 +35,10 @@   , transformMetric, transformMetric', dualCoCoProduct   , dualiseMetric, dualiseMetric'   , recipMetric, recipMetric'+  -- ** Eigenvectors   , eigenSpan, eigenSpan'   , eigenCoSpan, eigenCoSpan'+  , eigenSystem, HasEigenSystem, EigenVector   , metriNormalise, metriNormalise'   , metriScale', metriScale   , adjoint@@ -55,6 +57,7 @@   , Stiefel1(..)   , linMapAsTensProd, linMapFromTensProd   , covariance+  , outerProducts   ) where      @@ -166,15 +169,25 @@ --   describing the ellipsoid span of the vectors /e/&#x2080; and 2&#x22c5;/e/&#x2081;. --   Metrics generated this way are positive definite if no negative coefficients have --   been introduced with the '*^' scaling operator or with '^-^'.+--   +--   Note: @projector a ^+^ projector b ^+^ ...@ is more efficiently written as+--   @projectors [a, b, ...]@ projector :: HasMetric v => DualSpace v -> HerMetric v-projector u = matrixMetric $ HMat.outer uDecomp uDecomp- where uDecomp = asPackedVector u+projector u = HerMetric . pure $ u ⊗ u  projector' :: HasMetric v => v -> HerMetric' v-projector' v = matrixMetric' $ HMat.outer vDecomp vDecomp- where vDecomp = asPackedVector v+projector' v = HerMetric' . pure $ v ⊗ v +-- | Efficient shortcut for the 'sumV' of multiple 'projector's.+projectors :: HasMetric v => [DualSpace v] -> HerMetric v+projectors [] = zeroV+projectors us = HerMetric . pure . outerProducts $ zip us us +projector's :: HasMetric v => [v] -> HerMetric' v+projector's [] = zeroV+projector's vs = HerMetric' . pure . outerProducts $ zip vs vs++ singularMetric :: forall v . HasMetric v => HerMetric v singularMetric = matrixMetric $ HMat.scale (1/0) (HMat.ident dim)  where (Tagged dim) = dimension :: Tagged v Int@@ -305,14 +318,17 @@   --- | The eigenbasis of a /positive definite/ metric, with each eigenvector scaled---   to the square root of the eigenvalue.+-- | The eigenbasis of a metric, with each eigenvector scaled to the+--   square root of the eigenvalue. If the metric is not positive+--   definite (i.e. if it has zero eigenvalues), then the 'eigenSpan'+--   will contain zero vectors. --    --   This constitutes, in a sense, --   a decomposition of a metric into a set of 'projector'' vectors. If those---   are 'sumV'ed again, the original metric is obtained. (This holds even for---   non-Hilbert/Banach spaces, even though the concept of eigenbasis and---   &#x201c;scaled length&#x201d; doesn't really makes sense then in the usual way!)+--   are 'sumV'ed again (use 'projectors's' for this), then the original metric+--   is obtained. (This holds even for non-Hilbert/Banach spaces,+--   although the concept of eigenbasis and+--   &#x201c;scaled length&#x201d; doesn't really make sense there.) eigenSpan :: (HasMetric v, Scalar v ~ ℝ) => HerMetric' v -> [v] eigenSpan (HerMetric' Nothing) = [] eigenSpan (HerMetric' (Just (DenseLinear m))) = map fromPackedVector eigSpan@@ -325,18 +341,117 @@  where (μs,vsm) = HMat.eigSH' m        eigSpan = zipWith (HMat.scale . sqrt) (HMat.toList μs) (HMat.toColumns vsm) +-- | The reciprocal-space counterparts of the nonzero-EV eigenvectors, as can+--   be obtained from 'eigenSpan'. The systems of vectors/dual vectors+--   behave as orthonormal groups WRT each other, i.e. for each @f@+--   in @'eigenCoSpan' m@ there will be exactly one @v@ in @'eigenSpan' m@+--   such that @f<.>^v ≡ 1@; the other @f<.>^v@ pairings are zero.+-- +--   Furthermore, @'metric' m f ≡ 1@ for each @f@ in the co-span, which might+--   be seen as the actual defining characteristic of these span/co-span systems. eigenCoSpan :: (HasMetric v, Scalar v ~ ℝ) => HerMetric' v -> [DualSpace v] eigenCoSpan (HerMetric' Nothing) = [] eigenCoSpan (HerMetric' (Just (DenseLinear m))) = map fromPackedVector eigSpan  where (μs,vsm) = HMat.eigSH' m-       eigSpan = zipWith (HMat.scale . recip . sqrt) (HMat.toList μs) (HMat.toColumns vsm)+       eigSpan = map (uncurry $ HMat.scale . recip . sqrt)+                 . filter ((>0) . fst)+                 $ zip (HMat.toList μs) (HMat.toColumns vsm) eigenCoSpan' :: (HasMetric v, Scalar v ~ ℝ) => HerMetric v -> [v] eigenCoSpan' (HerMetric Nothing) = [] eigenCoSpan' (HerMetric (Just (DenseLinear m))) = map fromPackedVector eigSpan  where (μs,vsm) = HMat.eigSH' m-       eigSpan = zipWith (HMat.scale . recip . sqrt) (HMat.toList μs) (HMat.toColumns vsm)+       eigSpan = map (uncurry $ HMat.scale . recip . sqrt)+                 . filter ((>0) . fst)+                 $ zip (HMat.toList μs) (HMat.toColumns vsm)  +class HasEigenSystem m where+  type EigenVector m :: *+  -- | Generalised combination of 'eigenSpan' and 'eigenCoSpan'; this will give a+  --   maximum spanning set of vector-covector pairs @(f,v)@ such that @f<.>^v ≡ 1@+  --   and @metric m f ≡ 1@, whereas all @f@ and @v'@ from different tuples+  --   are orthogonal.+  --   It also yields the /kernel/ of singular metric, spanned by a set of stiefel-manifold+  --   points, i.e. vectors of unspecified length that correspond to the eigenvalue 0.+  -- +  --   You may also consider this as a /factorisation/ of a linear operator+  --   @𝐴 : 𝑉 → 𝑉'@ into mappings @𝑅 : 𝑉 → ℝⁿ@ and @𝐿 : ℝⁿ → 𝑉'@ (or, equivalently+  --   because ℝⁿ is a Hilbert space, @𝑅' : ℝⁿ → V'@ and @𝐿' : V → ℝⁿ@, which+  --   gives you an SVD-style inverse).+  eigenSystem :: m -> ( [Stiefel1 (EigenVector m)]+                      , [(EigenVector m, DualSpace (EigenVector m))] )++instance (HasMetric v, Scalar v ~ ℝ) => HasEigenSystem (HerMetric' v) where+  type EigenVector (HerMetric' v) = v+  eigenSystem (HerMetric' Nothing) = (fmap Stiefel1 completeBasisValues, [])+  eigenSystem (HerMetric' (Just (DenseLinear m))) = concat***concat $ unzip eigSpan+   where (μs,vsm) = HMat.eigSH' m+         eigSpan = zipWith (\μ v+                    -> if μ>0+                        then let sμ = sqrt μ+                             in ([], [( fromPackedVector $ HMat.scale sμ v+                                      , fromPackedVector $ HMat.scale (recip sμ) v )])+                        else ([Stiefel1 $ fromPackedVector v], [])+                   ) (HMat.toList μs) (HMat.toColumns vsm)++instance (HasMetric v, Scalar v ~ ℝ) => HasEigenSystem (HerMetric v) where+  type EigenVector (HerMetric v) = DualSpace v+  eigenSystem (HerMetric Nothing) = (fmap Stiefel1 completeBasisValues, [])+  eigenSystem (HerMetric (Just (DenseLinear m))) = concat***concat $ unzip eigSpan+   where (μs,vsm) = HMat.eigSH' m+         eigSpan = zipWith (\μ v+                    -> if μ>0+                        then let sμ = sqrt μ+                             in ([], [( fromPackedVector $ HMat.scale sμ v+                                      , fromPackedVector $ HMat.scale (recip sμ) v )])+                        else ([Stiefel1 $ fromPackedVector v], [])+                   ) (HMat.toList μs) (HMat.toColumns vsm)++instance (HasMetric v, Scalar v ~ ℝ) => HasEigenSystem (HerMetric' v, HerMetric' v) where+  type EigenVector (HerMetric' v, HerMetric' v) = v+  eigenSystem (n, HerMetric' (Just (DenseLinear m))) | not $ null nSpan+                                      = (++nKernel).concat***concat $ unzip eigSpan+   where (μs,vsm) = HMat.eigSH' $ fromv2ℝn HMat.<> m HMat.<> fromℝn2v'+                    -- m :: v' -> v+         eigSpan = zipWith (\μ v+                    -> if μ>0+                        then let sμ = sqrt μ+                             in ([], [( fromPackedVector $+                                        fromℝn2v HMat.#> HMat.scale sμ v+                                      , fromPackedVector $+                                        fromℝn2v' HMat.#> HMat.scale (recip sμ) v )+                                      ])+                        else ([Stiefel1 $ fromPackedVector v], [])+                   ) (HMat.toList μs) (HMat.toColumns vsm)+         fromv2ℝn = HMat.fromRows $ map (asPackedVector . snd) nSpan+         fromℝn2v' = HMat.tr fromv2ℝn+         fromℝn2v = HMat.fromColumns $ map (asPackedVector . fst) nSpan+         (nKernel, nSpan) = eigenSystem n+  eigenSystem (_, HerMetric' Nothing) = (fmap Stiefel1 completeBasisValues, [])++instance (HasMetric v, Scalar v ~ ℝ) => HasEigenSystem (HerMetric v, HerMetric v) where+  type EigenVector (HerMetric v, HerMetric v) = DualSpace v+  eigenSystem (n, HerMetric (Just (DenseLinear m))) | not $ null nSpan+                                      = (++nKernel).concat***concat $ unzip eigSpan+   where (μs,vsm) = HMat.eigSH' $ fromv'2ℝn HMat.<> m HMat.<> fromℝn2v+                    -- m :: v -> v'+         eigSpan = zipWith (\μ v+                    -> if μ>0+                        then let sμ = sqrt μ+                             in ([], [( fromPackedVector $+                                        fromℝn2v' HMat.#> HMat.scale sμ v+                                      , fromPackedVector $+                                        fromℝn2v HMat.#> HMat.scale (recip sμ) v )+                                      ])+                        else ([Stiefel1 $ fromPackedVector v], [])+                   ) (HMat.toList μs) (HMat.toColumns vsm)+         fromv'2ℝn = HMat.fromRows $ map (asPackedVector . snd) nSpan+         fromℝn2v = HMat.tr fromv'2ℝn+         fromℝn2v' = HMat.fromColumns $ map (asPackedVector . fst) nSpan+         (nKernel, nSpan) = eigenSystem n+  eigenSystem (_, HerMetric Nothing) = (fmap Stiefel1 completeBasisValues, [])++ -- | Constraint that a space's scalars need to fulfill so it can be used for 'HerMetric'. type MetricScalar s = ( SmoothScalar s                       , Ord s  -- We really rather wouldn't require this...@@ -562,13 +677,11 @@ factoriseMetric :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)                => HerMetric (v,w) -> (HerMetric v, HerMetric w) factoriseMetric (HerMetric Nothing) = (HerMetric Nothing, HerMetric Nothing)-factoriseMetric met = (sumV *** sumV) . unzip-                   $ (projector.fst &&& projector.snd) <$> eigenSpan' met+factoriseMetric met = (projectors *** projectors) . unzip $ eigenSpan' met  factoriseMetric' :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)                => HerMetric' (v,w) -> (HerMetric' v, HerMetric' w)-factoriseMetric' met = (sumV *** sumV) . unzip-                   $ (projector'.fst &&& projector'.snd) <$> eigenSpan met+factoriseMetric' met = (projector's *** projector's) . unzip $ eigenSpan met  productMetric :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)                => HerMetric v -> HerMetric w -> HerMetric (v,w)@@ -617,15 +730,16 @@  metricAsLength :: HerMetric ℝ -> ℝ metricAsLength m = case metricSq m 1 of-   o | o > 0    -> recip o-     | o < 0    -> error "Metric fails to be positive definite!"-     | o == 0   -> error "Trying to use zero metric as length."+   o | o > 0      -> sqrt $ recip o+     | o < 0      -> error "Metric fails to be positive definite!"+     | o == 0     -> error "Trying to use zero metric as length."+     | otherwise  -> error "Metric yields NaN."  metricFromLength :: ℝ -> HerMetric ℝ metricFromLength = projector . recip  metric'AsLength :: HerMetric' ℝ -> ℝ-metric'AsLength = recip . (`metric'`1)  -- do we really want `recip` here?+metric'AsLength = sqrt . (`metric'`1)   spanHilbertSubspace :: ∀ s v w@@ -678,7 +792,7 @@     | null eigSp  = showString "zeroV"     | otherwise   = showParen (p>5)                       . foldr1 ((.) . (.(" ^+^ "++)))-                      $ ((("projector "++).).showsPrec 6)<$>eigSp+                      $ ((("projector "++).).showsPrec 10)<$>eigSp    where eigSp = eigenSpan' m  instance (HasMetric v, Scalar v ~ Double, Show v) => Show (HerMetric' v) where@@ -686,7 +800,7 @@     | null eigSp  = showString "zeroV"     | otherwise   = showParen (p>5)                       . foldr1 ((.) . (.(" ^+^ "++)))-                      $ ((("projector' "++).).showsPrec 6)<$>eigSp+                      $ ((("projector' "++).).showsPrec 10)<$>eigSp    where eigSp = eigenSpan m  @@ -705,3 +819,27 @@                     => DualSpace v ⊗ w -> v:-*w linMapFromTensProd (DensTensProd m) = linear $                          asPackedVector >>> HMat.app m >>> fromPackedVector++++(⊗) :: (HasMetric v, FiniteDimensional w, Scalar v ~ s, Scalar w ~ s)+                    => w -> DualSpace v -> Linear s v w+w ⊗ v' = DenseLinear $ HMat.outer wDecomp v'Decomp+ where wDecomp = asPackedVector w+       v'Decomp = asPackedVector v'++outerProducts :: (HasMetric v, FiniteDimensional w, Scalar v ~ s, Scalar w ~ s)+                    => [(w, DualSpace v)] -> Linear s v w+outerProducts [] = zeroV+outerProducts pds = DenseLinear $ HMat.fromColumns (asPackedVector.fst<$>pds)+                          HMat.<> HMat.fromRows    (asPackedVector.snd<$>pds)++instance ∀ v w s . ( HasMetric v, FiniteDimensional w+                   , Show (DualSpace v), Show w, Scalar v ~ s, Scalar w ~ s )+    => Show (Linear s v w) where+  showsPrec p f = showParen (p>9) $ ("outerProducts "++)+        . shows [ (w, v' :: DualSpace v)+                | (v,v') <- zip completeBasisValues completeBasisValues+                , let w = f $ v ]+  +

@@ -32,6 +32,7 @@ {-# LANGUAGE FunctionalDependencies   #-} {-# LANGUAGE FlexibleContexts         #-} {-# LANGUAGE LiberalTypeSynonyms      #-}+{-# LANGUAGE DataKinds                #-} {-# LANGUAGE GADTs                    #-} {-# LANGUAGE RankNTypes               #-} {-# LANGUAGE TupleSections            #-}@@ -64,7 +65,7 @@             -- ** Local functions             , LocalLinear, LocalAffine             -- * Misc-            , palerp+            , palerp, LocallyCoercible(..)             ) where      @@ -78,6 +79,7 @@ import Data.Fixed  import Data.VectorSpace+import Data.Embedding import Data.LinearMap import Data.LinearMap.HerMetric import Data.LinearMap.Category@@ -219,6 +221,29 @@ class (PseudoAffine m, LinearManifold (Needle m), Interior m ~ m) => Manifold m instance (PseudoAffine m, LinearManifold (Needle m), Interior m ~ m) => Manifold m +++-- | Instances of this class must be diffeomorphic manifolds, and even have+--   /canonically isomorphic/ tangent spaces, so that+--   @'fromPackedVector' . 'asPackedVector' :: 'Needle' x -> 'Needle' ξ@+--   defines a meaningful “representational identity“ between these spaces.+class (PseudoAffine x, PseudoAffine ξ, Scalar (Needle x) ~ Scalar (Needle ξ))+         => LocallyCoercible x ξ where+  -- | Must be compatible with the canonical isomorphism on the tangent spaces,+  --   i.e.+  -- @+  -- locallyTrivialDiffeomorphism (p .+~^ 'fromPackedVector' v)+  --   ≡ locallyTrivialDiffeomorphism p .+~^ 'fromPackedVector' v+  -- @+  locallyTrivialDiffeomorphism :: x -> ξ+  +instance LocallyCoercible ℝ ℝ where locallyTrivialDiffeomorphism = id+instance LocallyCoercible (ℝ,ℝ) (ℝ,ℝ) where locallyTrivialDiffeomorphism = id+instance LocallyCoercible (ℝ,(ℝ,ℝ)) (ℝ,(ℝ,ℝ)) where locallyTrivialDiffeomorphism = id+instance LocallyCoercible ((ℝ,ℝ),ℝ) ((ℝ,ℝ),ℝ) where locallyTrivialDiffeomorphism = id+++ type LocallyScalable s x = ( PseudoAffine x                            , HasMetric (Needle x)                            , s ~ Scalar (Needle x) )@@ -268,8 +293,8 @@ type EuclidSpace x = ( AffineManifold x, InnerSpace (Diff x)                      , DualSpace (Diff x) ~ Diff x, Floating (Scalar (Diff x)) ) -euclideanMetric :: EuclidSpace x => Tagged x (Metric x)-euclideanMetric = Tagged euclideanMetric'+euclideanMetric :: EuclidSpace x => proxy x -> Metric x+euclideanMetric _ = euclideanMetric'   -- | A co-needle can be understood as a “paper stack”, with which you can measure@@ -331,6 +356,14 @@   (.+~^) = (.+^) instance SmoothScalar s => PseudoAffine (FinVecArrRep t b s) where   a.-~.b = pure (a.-.b)+instance SmoothScalar s => LocallyCoercible (FinVecArrRep t b s) (FinVecArrRep t b s) where+  locallyTrivialDiffeomorphism = id+instance (SmoothScalar s, LinearManifold b, Scalar b ~ s)+           => LocallyCoercible (FinVecArrRep t b s) b where+  locallyTrivialDiffeomorphism = (concreteArrRep$<-$)+instance (SmoothScalar s, LinearManifold b, Scalar b ~ s)+           => LocallyCoercible b (FinVecArrRep t b s) where+  locallyTrivialDiffeomorphism = (concreteArrRep$->$)     instance Semimanifold (ZeroDim k) where@@ -360,6 +393,10 @@                 Tagged tb = translateP :: Tagged b (Interior b -> Needle b -> Interior b) instance (PseudoAffine a, PseudoAffine b) => PseudoAffine (a,b) where   (a,b).-~.(c,d) = liftA2 (,) (a.-~.c) (b.-~.d)+instance (PseudoAffine a, PseudoAffine b, PseudoAffine c)+     => LocallyCoercible (a,(b,c)) ((a,b),c) where locallyTrivialDiffeomorphism = regroup+instance (PseudoAffine a, PseudoAffine b, PseudoAffine c)+     => LocallyCoercible ((a,b),c) (a,(b,c)) where locallyTrivialDiffeomorphism = regroup'  instance (Semimanifold a, Semimanifold b, Semimanifold c) => Semimanifold (a,b,c) where   type Needle (a,b,c) = (Needle a, Needle b, Needle c)@@ -379,6 +416,18 @@                 Tagged tc = translateP :: Tagged c (Interior c -> Needle c -> Interior c) instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) => PseudoAffine (a,b,c) where   (a,b,c).-~.(d,e,f) = liftA3 (,,) (a.-~.d) (b.-~.e) (c.-~.f)+instance (PseudoAffine a, PseudoAffine b, PseudoAffine c)+     => LocallyCoercible (a,b,c) ((a,b),c) where+  locallyTrivialDiffeomorphism (a,b,c) = ((a,b),c)+instance (PseudoAffine a, PseudoAffine b, PseudoAffine c)+     => LocallyCoercible (a,b,c) (a,(b,c)) where+  locallyTrivialDiffeomorphism (a,b,c) = (a,(b,c))+instance (PseudoAffine a, PseudoAffine b, PseudoAffine c)+     => LocallyCoercible ((a,b),c) (a,b,c) where+  locallyTrivialDiffeomorphism ((a,b),c) = (a,b,c)+instance (PseudoAffine a, PseudoAffine b, PseudoAffine c)+     => LocallyCoercible (a,(b,c)) (a,b,c) where+  locallyTrivialDiffeomorphism (a,(b,c)) = (a,b,c)  instance (MetricScalar a, KnownNat n) => Semimanifold (FreeVect n a) where   type Needle (FreeVect n a) = FreeVect n a@@ -388,6 +437,11 @@   (.+~^) = (.+^) instance (MetricScalar a, KnownNat n) => PseudoAffine (FreeVect n a) where   a.-~.b = pure (a.-.b)+instance LocallyCoercible ℝ (ℝ ^ S Z) where+  locallyTrivialDiffeomorphism = replicVector+instance LocallyCoercible (ℝ ^ S Z) ℝ where+  locallyTrivialDiffeomorphism = (<.>^replicVector 1)+  instance (HasMetric a, FiniteDimensional b, Scalar a~Scalar b) => Semimanifold (a⊗b) where   type Needle (a⊗b) = a ⊗ b

@@ -22,9 +22,12 @@ {-# LANGUAGE RankNTypes                 #-} {-# LANGUAGE TupleSections              #-} {-# LANGUAGE ParallelListComp           #-}+{-# LANGUAGE MonadComprehensions        #-} {-# LANGUAGE UnicodeSyntax              #-} {-# LANGUAGE ConstraintKinds            #-} {-# LANGUAGE PatternGuards              #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ViewPatterns               #-} {-# LANGUAGE LambdaCase                 #-} {-# LANGUAGE TypeOperators              #-} {-# LANGUAGE ScopedTypeVariables        #-}@@ -35,26 +38,27 @@  module Data.Manifold.TreeCover (        -- * Shades -         Shade(..), Shade'(..), IsShade+         Shade(..), pattern(:±), Shade'(..), (|±|), IsShade        -- ** Lenses        , shadeCtr, shadeExpanse, shadeNarrowness        -- ** Construction-       , fullShade, fullShade', pointsShades+       , fullShade, fullShade', pointsShades, pointsCovers        -- ** Evaluation        , occlusion        -- ** Misc-       , factoriseShade, intersectShade's+       , factoriseShade, intersectShade's, Refinable, refineShade', convolveShade', coerceShade        -- * Shade trees-       , ShadeTree(..), fromLeafPoints-       -- * Simple view helpers-       , onlyNodes, onlyLeaves+       , ShadeTree(..), fromLeafPoints, onlyLeaves, indexShadeTree+       -- * View helpers+       , onlyNodes        -- ** Auxiliary types        , SimpleTree, Trees, NonEmptyTree, GenericTree(..)        -- * Misc        , sShSaw, chainsaw, HasFlatView(..), shadesMerge, smoothInterpolate        , twigsWithEnvirons, completeTopShading, flexTwigsShading        , WithAny(..), Shaded, stiAsIntervalMapping, spanShading-       , DifferentialEqn, filterDEqnSolution_static+       , constShaded, stripShadedUntopological+       , DifferentialEqn, filterDEqnSolution_loc        -- ** Triangulation-builders        , TriangBuild, doTriangBuild, singleFullSimplex, autoglueTriangulation        , AutoTriang, elementaryTriang, breakdownAutoTriang@@ -90,6 +94,7 @@ import Data.Manifold.PseudoAffine import Data.Function.Differentiable import Data.Function.Differentiable.Data+import Data.SetLike.Intersection      import Data.Embedding import Data.CoNat@@ -158,6 +163,7 @@   factoriseShade :: ( Manifold x, RealDimension (Scalar (Needle x))                     , Manifold y, RealDimension (Scalar (Needle y)) )                 => shade (x,y) -> (shade x, shade y)+  coerceShade :: (Manifold x, Manifold y, LocallyCoercible x y) => shade x -> shade y  instance IsShade Shade where   shadeCtr f (Shade c e) = fmap (`Shade`e) $ f c@@ -170,6 +176,8 @@          δinv = recipMetric δ   factoriseShade (Shade (x₀,y₀) δxy) = (Shade x₀ δx, Shade y₀ δy)    where (δx,δy) = factoriseMetric' δxy+  coerceShade (Shade x (HerMetric' δxym))+          = Shade (locallyTrivialDiffeomorphism x) (HerMetric' $ unsafeCoerceLinear<$>δxym)  shadeExpanse :: Functor f (->) (->) => (Metric' x -> f (Metric' x)) -> Shade x -> f (Shade x) shadeExpanse f (Shade c e) = fmap (Shade c) $ f e@@ -184,6 +192,8 @@            _               -> zeroV   factoriseShade (Shade' (x₀,y₀) δxy) = (Shade' x₀ δx, Shade' y₀ δy)    where (δx,δy) = factoriseMetric δxy+  coerceShade (Shade' x (HerMetric δxym))+          = Shade' (locallyTrivialDiffeomorphism x) (HerMetric $ unsafeCoerceLinear<$>δxym)  shadeNarrowness :: Functor f (->) (->) => (Metric x -> f (Metric x)) -> Shade' x -> f (Shade' x) shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e@@ -202,6 +212,24 @@ fullShade' :: WithField ℝ Manifold x => x -> Metric x -> Shade' x fullShade' ctr expa = Shade' ctr expa ++-- | Span a 'Shade' from a center point and multiple deviation-vectors.+pattern (:±) :: () => WithField ℝ Manifold x => x -> [Needle x] -> Shade x+pattern x :± shs <- Shade x (eigenSpan -> shs)+ where x :± shs = fullShade x $ projector's 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 $ projectors [v^/(v<.>v) | v<-shs]+++ subshadeId' :: WithField ℝ Manifold x                    => x -> NonEmpty (Needle' x) -> x -> (Int, HourglassBulb) subshadeId' c expvs x = case x .-~. c of@@ -215,20 +243,40 @@                    --- | Attempt to find a 'Shade' that &#x201c;covers&#x201d; the given points.+-- | 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.+--   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 list result.+--   Hence the result type is a list. pointsShades :: WithField ℝ Manifold x => [x] -> [Shade x] pointsShades = map snd . pointsShades' zeroV +-- | 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 ℝ Manifold x => [x] -> [Shade x]+pointsCovers = map guaranteeIn . pointsShades' zeroV+ where guaranteeIn (ps, Shade x₀ ex) +          = case ps >>= \p -> let Option (Just v) = p.-~.x₀+                              in guard (metric ex' v > 1) >> [(p,projector' v)]+             of []   -> Shade x₀ ex+                outs -> guaranteeIn ( fst<$>outs+                                    , Shade x₀+                                         $ ex ^+^ sumV (snd<$>outs)+                                                    ^/ fromIntegral (2 * length outs) )+        where ex' = recipMetric ex+ pointsShade's :: WithField ℝ Manifold x => [x] -> [Shade' x] pointsShade's = map (\(Shade c e) -> Shade' c $ recipMetric e) . pointsShades +pointsCover's :: WithField ℝ Manifold x => [x] -> [Shade' x]+pointsCover's = map (\(Shade c e) -> Shade' c $ recipMetric e) . pointsCovers+ pseudoECM :: WithField ℝ Manifold x => NonEmpty x -> (x, ([x],[x])) pseudoECM (p₀ NE.:| psr) = foldl' ( \(acc, (rb,nr)) (i,p)                                   -> case p.-~.acc of @@ -245,7 +293,7 @@                            _ -> pointsShades' minExt inc'd                                   ++ pointsShades' minExt unreachable  where (ctr,(inc'd,unreachable)) = pseudoECM $ NE.fromList ps-       expa = ( (^+^minExt) . (^/ fromIntegral(length ps)) . sumV . map projector' )+       expa = ( (^+^minExt) . (^/ fromIntegral(length ps)) . projector's )               <$> mapM (.-~.ctr) ps         @@ -273,7 +321,7 @@                   = Just $ let cc = c₂ .+~^ v ^/ 2                                Option (Just cv₁) = c₁.-~.cc                                Option (Just cv₂) = c₂.-~.cc-                           in Shade cc . sumV $ [e₁, e₂] ++ (projector'<$>[cv₁, cv₂])+                           in Shade cc $ e₁ ^+^ e₂ ^+^ projector's [cv₁, cv₂]            | otherwise  = Nothing shadesMerge _ shs = shs @@ -365,21 +413,24 @@        uds = directionChoices hs  traverseDirectionChoices :: (WithField ℝ Manifold x, Hask.Applicative f)-               => (     (Needle' x, ShadeTree x)-                    -> [(Needle' x, ShadeTree x)]+               => (    (Int, (Needle' x, ShadeTree x))+                    -> [(Int, (Needle' x, ShadeTree x))]                     -> f (ShadeTree x) )                  -> [DBranch x]                  -> f [DBranch x]-traverseDirectionChoices f dbs = td [] (dbs >>=-                            \(DBranch ѧ (Hourglass τ β))-                              -> [(ѧ,τ), (negateV ѧ,β)])- where td pds ((ѧ,t):(v,b):vds)+traverseDirectionChoices f dbs+           = td [] . scanLeafNums 0+               $ dbs >>= \(DBranch ѧ (Hourglass τ β))+                              -> [(ѧ,τ), (negateV ѧ,β)]+ where td pds (ѧt@(_,(ѧ,_)):vb:vds)          = liftA3 (\t' b' -> (DBranch ѧ (Hourglass t' b') :))-             (f (ѧ,t) $ pds++(v,b):uds)-             (f (v,b) $ pds++(ѧ,t):uds)-             $ td ((ѧ,t):(v,b):pds) vds+             (f ѧt $ vb:uds)+             (f vb $ ѧt:uds)+             $ td (ѧt:vb:pds) vds         where uds = pds ++ vds        td _ _ = pure []+       scanLeafNums _ [] = []+       scanLeafNums i₀ ((v,t):vts) = (i₀, (v,t)) : scanLeafNums (i₀ + nLeaves t) vts   instance (NFData x, NFData (Needle' x)) => NFData (ShadeTree x) where@@ -419,14 +470,34 @@  -- | Build a quite nicely balanced tree from a cloud of points, on any real manifold. -- ---   Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/generate-ShadeTrees.ipynb#pseudorandomCloudTree+--   Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/Trees-and-Webs.ipynb#pseudorandomCloudTree --  -- <<images/examples/simple-2d-ShadeTree.png>> fromLeafPoints :: ∀ x. WithField ℝ Manifold x => [x] -> ShadeTree x fromLeafPoints = fromLeafPoints' sShIdPartition  +indexShadeTree :: ∀ x . WithField ℝ Manifold x+       => ShadeTree x -> Int -> Either Int ([ShadeTree x], x)+indexShadeTree _ i+    | i<0        = Left i+indexShadeTree sh@(PlainLeaves lvs) i = case length lvs of+  n | i<n       -> Right ([sh], lvs!!i)+    | otherwise -> Left $ i-n+indexShadeTree (DisjointBranches n brs) i+    | i<n        = foldl (\case +                             Left i' -> (`indexShadeTree`i')+                             result  -> return result+                         ) (Left i) brs+    | otherwise  = Left $ i-n+indexShadeTree sh@(OverlappingBranches n _ brs) i+    | i<n        = first (sh:) <$> foldl (\case +                             Left i' -> (`indexShadeTree`i')+                             result  -> return result+                         ) (Left i) (toList brs>>=toList)+    | otherwise  = Left $ i-n + fromFnGraphPoints :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y)                      => [(x,y)] -> ShadeTree (x,y) fromFnGraphPoints = fromLeafPoints' fg_sShIdPart@@ -521,7 +592,7 @@   trunks :: ∀ x. WithField ℝ Manifold x => ShadeTree x -> [Shade x]-trunks (PlainLeaves lvs) = pointsShades lvs+trunks (PlainLeaves lvs) = pointsCovers lvs trunks (DisjointBranches _ brs) = Hask.foldMap trunks brs trunks (OverlappingBranches _ sh _) = [sh] @@ -558,11 +629,19 @@       in overlappingBranches (fs sh) brs'  -intersectShade's :: ∀ y . WithField ℝ Manifold y => [Shade' y] -> Option (Shade' y)-intersectShade's [] = error "Global `Shade'` not implemented, so can't do intersection of zero co-shades."-intersectShade's (sh:shs) = Hask.foldrM inter2 sh shs- where inter2 :: Shade' y -> Shade' y -> Option (Shade' y)-       inter2 (Shade' c e) (Shade' ζ η)+-- | 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 ℝ Manifold 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) tsh = all ((<1) . minusLogOcclusion' tsh)+                                  [ ac.+~^σ*^v | σ<-[0,1], v<-eigenCoSpan' ae ]+  +  -- | Specialised intersection operation. If @p@ is in @a@ and @b@, then it is+  --   also in @refineShade' a b@. (The converse may not hold.)+  refineShade' :: Shade' y -> Shade' y -> Option (Shade' y)+  refineShade' (Shade' c e) (Shade' ζ η)            | μe < 1 && μη < 1  = return $ Shade' iCtr iExpa            | otherwise         = empty         where [c', ζ'] = [ ctr.+~^linearCombo@@ -588,42 +667,92 @@               μe = rcⰰ<.>^rc               μη = rζⰰ<.>^rζ               iExpa = (e^+^η)^/2 ^+^ projector rcⰰ^/(1-μe) ^+^ projector rζⰰ^/(1-μη)+  +  -- | 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>>+  convolveShade' :: Shade' y -> Shade' (Needle y) -> Shade' y+  convolveShade' (Shade' y₀ ey) (Shade' δ₀ eδ)+          = Shade' (y₀.+~^δ₀)+                   ( projectors [ f ^* ζ crl+                                | (f,_) <- eδsp+                                | crl <- corelap ] )+   where (_,eδsp) = eigenSystem (ey,eδ)+         corelap = map (metric ey . snd) eδsp+         ζ = case filter (>0) corelap 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+                        0  -> 0+                        sq -> edgeFactor / (recip sq + 1)+   +instance Refinable ℝ where+  refineShade' (Shade' cl el) (Shade' cr er)+         = case (metricSq el 1, metricSq 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 (projector $ 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, Refinable b) => Refinable (a,b)+  +                             +intersectShade's :: ∀ y . Refinable y => NonEmpty (Shade' y) -> Option (Shade' y)+intersectShade's (sh:|shs) = Hask.foldrM refineShade' sh shs -type DifferentialEqn x y = Shade' (x,y) -> Shade' (LocalLinear x y)  -filterDEqnSolution_loc :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y)-           => DifferentialEqn x y -> (Shade' (x,y), [Shade' (x,y)])-                   -> Option (Shade' y, LocalLinear x y)-filterDEqnSolution_loc f (shxy@(Shade' (x,y) expa), neighbours) = (,j₀) <$> yc++type DifferentialEqn x y = Shade (x,y) -> Shade' (LocalLinear x y)+++filterDEqnSolution_loc :: ∀ x y . (WithField ℝ Manifold x, Refinable y)+           => DifferentialEqn x y -> ((x, Shade' y), NonEmpty (x, Shade' y))+                   -> Option (Shade' y)+filterDEqnSolution_loc f ((x, shy@(Shade' y expay)), neighbours) = yc  where jShade@(Shade' j₀ jExpa) = f shxy-       marginδs :: [(Needle x, (Needle y, Metric y))]+       [shxy] = pointsCovers [ (xs, ys')+                             | (xs, Shade' ys yse) <- (x,shy):NE.toList neighbours+                             , δy <- eigenCoSpan' yse+                             , ys' <- [ys.+~^δy, ys.-~^δy] ]+       [Shade' _ expax] = pointsCover's $ x : (fst<$>NE.toList neighbours)+       marginδs :: NonEmpty (Needle x, (Needle y, Metric y))        marginδs = [ (δxm, (δym, expany))-                  | Shade' (xn, yn) expan <- neighbours-                  , let (Option (Just δx)) = x.-~.xn-                        (expanx, expany) = factoriseMetric expan-                        (Option (Just yc'n))-                               = covariance $ recipMetric' expan-                        xntoMarg = metriNormalise expanx δx-                        (Option (Just δxm))-                           = (xn .+~^ xntoMarg :: x) .-~. x-                        (Option (Just δym))-                           = (yn .+~^ (yc'n $ xntoMarg) :: y-                               ) .-~. y+                  | (xn, Shade' yn expany) <- neighbours+                  , let (Option (Just δxm)) = xn.-~.x+                        (Option (Just δym)) = yn.-~.y                   ]        back2Centre :: (Needle x, (Needle y, Metric y)) -> Shade' y        back2Centre (δx, (δym, expany))-            = Shade' (y.+~^δyb) . recipMetric-                $ recipMetric' expany-                  ^+^ recipMetric' (applyLinMapMetric jExpa δx')+            = convolveShade'+                (Shade' y expany)+                (Shade' δyb $ applyLinMapMetric jExpa δx')         where δyb = δym ^-^ (j₀ $ δx)               δx' = toDualWith expax δx        yc :: Option (Shade' y)        yc = intersectShade's $ back2Centre <$> marginδs-       (expax, expay) = factoriseMetric expa        xSpan = eigenCoSpan' expax  @@ -631,8 +760,8 @@ -- The simple list-yielding version (see rev. b4a427d59ec82889bab2fde39225b14a57b694df -- may well be more efficient than this version via a traversal. twigsWithEnvirons :: ∀ x. WithField ℝ Manifold x-    => ShadeTree x -> [(ShadeTree x, [ShadeTree x])]-twigsWithEnvirons = execWriter . traverseTwigsWithEnvirons (writer . (fst&&&pure))+    => ShadeTree x -> [((Int, ShadeTree x), [(Int, ShadeTree x)])]+twigsWithEnvirons = execWriter . traverseTwigsWithEnvirons (writer . (snd.fst&&&pure))  data OuterMaybeT f a = OuterNothing | OuterJust (f a) deriving (Hask.Functor) instance (Hask.Applicative f) => Hask.Applicative (OuterMaybeT f) where@@ -642,91 +771,73 @@  traverseTwigsWithEnvirons :: ∀ x f .             (WithField ℝ Manifold x, Hask.Applicative f)-    => ((ShadeTree x, [ShadeTree x]) -> f (ShadeTree x))+    => ( ((Int, ShadeTree x), [(Int, ShadeTree x)]) -> f (ShadeTree x))          -> ShadeTree x -> f (ShadeTree x)-traverseTwigsWithEnvirons f = fst . go []- where go :: [ShadeTree x] -> ShadeTree x -> (f (ShadeTree x), Bool)-       go _ (DisjointBranches nlvs djbs) = ( fmap (DisjointBranches nlvs)-                                               $ Hask.traverse (fst . go []) djbs-                                           , False )-       go envi ct@(OverlappingBranches nlvs rob@(Shade robc _) brs)+traverseTwigsWithEnvirons f = fst . go [] . (0,)+ where go :: [(Int, ShadeTree x)] -> (Int, ShadeTree x)+                          -> (f (ShadeTree x), Bool)+       go _ (i₀, DisjointBranches nlvs djbs) = ( fmap (DisjointBranches nlvs)+                                                   . Hask.traverse (fst . go [])+                                                   $ NE.zip ioffs djbs+                                               , False )+        where ioffs = NE.scanl (\i -> (+i) . nLeaves) i₀ djbs+       go envi ct@(i₀, (OverlappingBranches nlvs rob@(Shade robc _) brs))                 = ( case descentResult of                      OuterNothing -> f-                         $ purgeRemotes (ct, Hask.foldMap (twigProximæ robc) envi)+                         $ purgeRemotes+                            (ct, Hask.foldMap (\(io,te)+                                            -> first (+io) <$> twigProximæ robc te) envi)                      OuterJust dR -> fmap (OverlappingBranches nlvs rob . NE.fromList) dR                   , False )         where descentResult = traverseDirectionChoices tdc $ NE.toList brs-              tdc (vy, ty) alts = case go envi'' ty of+              tdc (io, (vy, ty)) alts = case go envi'' (i₀+io, ty) of                                    (_, True) -> OuterNothing                                    (down, _) -> OuterJust down-               where envi'' = filter (trunks >>> \(Shade ce _:_)+               where envi'' = filter (snd >>> trunks >>> \(Shade ce _:_)                                          -> let Option (Just δyenv) = ce.-~.robc                                                 qq = vy<.>^δyenv                                             in qq > -1 && qq < 5                                        ) envi'-                              ++ map snd alts+                              ++ map ((+i₀)***snd) alts               envi' = approach =<< envi-              approach apt@(OverlappingBranches _ (Shade envc _) _)-                  = twigsaveTrim hither apt+              approach (i₀e, apt@(OverlappingBranches _ (Shade envc _) _))+                  = first (+i₀e) <$> twigsaveTrim hither apt                where Option (Just δxenv) = robc .-~. envc                      hither (DBranch bdir (Hourglass bdc₁ bdc₂))-                       | bdir<.>^δxenv > 0  = [bdc₁]-                       | otherwise          = [bdc₂]+                       | bdir<.>^δxenv > 0  = [(0           , bdc₁)]+                       | otherwise          = [(nLeaves bdc₁, bdc₂)]               approach q = [q]-       go envi plvs@(PlainLeaves _) = (f $ purgeRemotes (plvs, envi), True)+       go envi plvs@(i₀, (PlainLeaves _))+                         = (f $ purgeRemotes (plvs, envi), True)        -       twigProximæ :: x -> ShadeTree x -> [ShadeTree x]-       twigProximæ x₀ (DisjointBranches _ djbs) = Hask.foldMap (twigProximæ x₀) djbs+       twigProximæ :: x -> ShadeTree x -> [(Int, ShadeTree x)]+       twigProximæ x₀ (DisjointBranches _ djbs)+               = Hask.foldMap (\(i₀,st) -> first (+i₀) <$> twigProximæ x₀ st)+                    $ NE.zip ioffs djbs+        where ioffs = NE.scanl (\i -> (+i) . nLeaves) 0 djbs        twigProximæ x₀ ct@(OverlappingBranches _ (Shade xb qb) brs)                    = twigsaveTrim hither ct         where Option (Just δxb) = x₀ .-~. xb               hither (DBranch bdir (Hourglass bdc₁ bdc₂))                  | bdir<.>^δxb > 0  = twigProximæ x₀ bdc₁-                 | otherwise        = twigProximæ x₀ bdc₂-       twigProximæ _ plainLeaves = [plainLeaves]+                 | otherwise        = first (+nLeaves bdc₁)+                                     <$> twigProximæ x₀ bdc₂+       twigProximæ _ plainLeaves = [(0, plainLeaves)]        -       twigsaveTrim :: (DBranch x -> [ShadeTree x])-                       -> ShadeTree x -> [ShadeTree x]+       twigsaveTrim :: (DBranch x -> [(Int,ShadeTree x)])+                       -> ShadeTree x -> [(Int,ShadeTree x)]        twigsaveTrim f ct@(OverlappingBranches _ _ dbs)-                 = case Hask.mapM (f >>> noLeaf) dbs of+                 = case Hask.mapM (\(i₀,dbr) -> noLeaf $ first(+i₀)<$>f dbr)+                                 $ NE.zip ioffs dbs of                       Just pqe -> Hask.fold pqe-                      _        -> [ct]-        where noLeaf [PlainLeaves _] = empty+                      _        -> [(0,ct)]+        where noLeaf [(_,PlainLeaves _)] = empty               noLeaf bqs = pure bqs+              ioffs = NE.scanl (\i -> (+i) . sum . fmap nLeaves . toList) 0 dbs        -       purgeRemotes :: (ShadeTree x, [ShadeTree x]) -> (ShadeTree x, [ShadeTree x])-       purgeRemotes (ctm@(OverlappingBranches _ sm@(Shade xm _) _), candidates)-                                       = (ctm, filter unobscured closeby)-        where closeby = filter proximate candidates-              proximate (OverlappingBranches _ sh@(Shade xh _) _)-                    = minusLogOcclusion sh xm * minusLogOcclusion sm xh-                       < 1024  -- = (2⋅4²)².  The four-radius occlusion occurs-                               -- if two 𝑟-sized shades have just enough space-                               -- to fit another 𝑟-shade between them; then-                               -- we don't consider the shades neighbours-                               -- anymore. A factor √2 for the discrepancy-                               -- between standard deviation and max distance.-              proximate _ = True-              unobscured ht@(OverlappingBranches _ (Shade xh _) _)-                     = all (don'tObscure (xh, onlyLeaves ht)) closeby-              don'tObscure (xh,lvsh) (OverlappingBranches _ sb@(Shade xb eb) _)-                          = vmc⋅vhc >= 0 || vm⋅vh >= 0-               where Option (Just vm) = pbm .-~. xb-                     Option (Just vh) = pbh .-~. xb-                     Option (Just vmc) = xm .-~. xb-                     Option (Just vhc) = xh .-~. xb-                     [pbm, pbh] = [ maximumBy (comparing $ \l ->-                                               let Option (Just w) = l.-~.xb-                                               in v⋅w ) lvs-                                  | lvs <- [lvsm, lvsh]-                                  | v <- [vhc, vmc] ]-                     (⋅) :: Needle x -> Needle x -> ℝ-                     v⋅w = toDualWith mb v <.>^ w-                     mb = recipMetric eb-              don'tObscure _ _ = True-              lvsm = onlyLeaves ctm-       purgeRemotes xyz = xyz-    +       purgeRemotes :: ((Int,ShadeTree x), [(Int,ShadeTree x)])+                    -> ((Int,ShadeTree x), [(Int,ShadeTree x)])+       purgeRemotes = id -- See 7d1f3a4 for the implementation; this didn't work reliable.       completeTopShading :: (WithField ℝ Manifold x, WithField ℝ Manifold y)                    => x`Shaded`y -> [Shade' (x,y)]@@ -734,7 +845,7 @@                      = pointsShade's $ (_topological &&& _untopological) <$> plvs completeTopShading (DisjointBranches _ bqs)                      = take 1 . completeTopShading =<< NE.toList bqs-completeTopShading t = pointsShade's . map (_topological &&& _untopological) $ onlyLeaves t+completeTopShading t = pointsCover's . map (_topological &&& _untopological) $ onlyLeaves t  flexTopShading :: ∀ x y f . ( WithField ℝ Manifold x, WithField ℝ Manifold y                             , Applicative f (->) (->) )@@ -769,25 +880,8 @@                   => (Shade' (x,y) -> f (x, (Shade' y, LocalLinear x y)))                       -> x`Shaded`y -> f (x`Shaded`y) flexTwigsShading f = traverseTwigsWithEnvirons locFlex- where locFlex :: ∀ μ . (x`Shaded`y, μ) -> f (x`Shaded`y)-       locFlex (lsh, _) = flexTopShading f lsh--filterDEqnSolution_static :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y)-           => DifferentialEqn x y-               -> x`Shaded`y -> Option (x`Shaded`y)-filterDEqnSolution_static deq tr = traverseTwigsWithEnvirons locSoltn tr- where locSoltn :: (x`Shaded`y, [x`Shaded`y]) -> Option (x`Shaded`y)-       locSoltn (local, environs) = do-            let enviShades = completeTopShading =<< environs-            flexed <- flexTopShading-                           (\oSh@(Shade' (ox,_) _) -> -                              (ox,) <$> filterDEqnSolution_loc deq (oSh, enviShades)-                           ) local-            top'@(Shade' (top'x,_) top'exp)-                     <- intersectShade's $ completeTopShading =<< [local, flexed]-            let (_, top'ySh) = factoriseShade top'-            j' <- covariance $ recipMetric' top'exp-            flexTopShading (const $ pure (top'x, (top'ySh, j'))) flexed+ where locFlex :: ∀ μ . ((Int, x`Shaded`y), μ) -> f (x`Shaded`y)+       locFlex ((_,lsh), _) = flexTopShading f lsh                   @@ -1085,7 +1179,7 @@ elementaryTriang :: ∀ n n' x . (KnownNat n', n~S n', WithField ℝ EuclidSpace x)                       => Simplex n x -> AutoTriang n x elementaryTriang t = AutoTriang (fullOpenSimplex m t >> return ())- where (Tagged m) = euclideanMetric :: Tagged x (Metric x)+ where m = euclideanMetric t  breakdownAutoTriang :: ∀ n n' x . (KnownNat n', n ~ S n') => AutoTriang n x -> [Simplex n x] breakdownAutoTriang (AutoTriang t) = doTriangBuild t@@ -1345,6 +1439,12 @@ 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)++stripShadedUntopological :: x`Shaded`y -> ShadeTree x+stripShadedUntopological = unsafeFmapTree (fmap _topological) id shadeWithoutAnything+ -- | This is to 'ShadeTree' as 'Data.Map.Map' is to 'Data.Set.Set'. type x`Shaded`y = ShadeTree (x`WithAny`y) @@ -1380,7 +1480,7 @@  stiAsIntervalMapping :: (x ~ ℝ, y ~ ℝ)             => x`Shaded`y -> [(x, ((y, Diff y), Linear ℝ x y))]-stiAsIntervalMapping = twigsWithEnvirons >=> pure.fst >=> completeTopShading >=> pure.+stiAsIntervalMapping = twigsWithEnvirons >=> pure.snd.fst >=> completeTopShading >=> pure.              \(Shade' (xloc, yloc) shd)                  -> ( xloc, ( (yloc, recip $ metric shd (0,1))                             , case covariance (recipMetric' shd) of@@ -1404,7 +1504,7 @@  where addYs :: NonEmpty x -> NonEmpty (x`WithAny`y)        addYs l = foldr (NE.<|) (fmap ( WithAny ymid) l     )                                (fmap (`WithAny`xmid) yexamp)-          where [xsh@(Shade xmid _)] = pointsShades $ toList l+          where [xsh@(Shade xmid _)] = pointsCovers $ toList l                 Shade ymid yexpa = f xsh                 yexamp = [ ymid .+~^ σ*^δy                          | δy <- eigenSpan yexpa, σ <- [-1,1] ]

@@ -54,7 +54,7 @@         , Cutplane(..)         , fathomCutDistance, sideOfCut         -- * Linear mappings-        , Linear, denseLinear+        , Linear, LocalLinear, denseLinear    ) where  

@@ -0,0 +1,219 @@+-- |+-- Module      : Data.Manifold.Web+-- Copyright   : (c) Justus Sagemüller 2016+-- License     : GPL v3+-- +-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de+-- Stability   : experimental+-- Portability : portable+-- +{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE UndecidableInstances       #-}+{-# LANGUAGE StandaloneDeriving         #-}+{-# LANGUAGE DeriveGeneric              #-}+{-# LANGUAGE DeriveFunctor              #-}+{-# LANGUAGE DeriveFoldable             #-}+{-# LANGUAGE DeriveTraversable          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE FunctionalDependencies     #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE GADTs                      #-}+{-# LANGUAGE RankNTypes                 #-}+{-# LANGUAGE TupleSections              #-}+{-# LANGUAGE ParallelListComp           #-}+{-# LANGUAGE UnicodeSyntax              #-}+{-# LANGUAGE ConstraintKinds            #-}+{-# LANGUAGE PatternGuards              #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ViewPatterns               #-}+{-# LANGUAGE LambdaCase                 #-}+{-# LANGUAGE TypeOperators              #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE LiberalTypeSynonyms        #-}+{-# LANGUAGE RecordWildCards            #-}+{-# LANGUAGE DataKinds                  #-}+++module Data.Manifold.Web where+++import Data.List hiding (filter, all, elem, sum, foldr1)+import Data.Maybe+import qualified Data.Set as Set+import qualified Data.Vector as Arr+import qualified Data.Vector.Unboxed as UArr+import Data.List.NonEmpty (NonEmpty(..))+import Data.List.FastNub+import qualified Data.List.NonEmpty as NE+import Data.Semigroup+import Data.Ord (comparing)+import Control.DeepSeq++import Data.VectorSpace+import Data.AffineSpace+import Data.LinearMap+import Data.LinearMap.HerMetric+import Data.LinearMap.Category+import Data.AffineSpace+import Data.Basis+import Data.Complex hiding (magnitude)+import Data.Void+import Data.Tagged+import Data.Proxy++import Data.SimplicialComplex+import Data.Manifold.Types+import Data.Manifold.Types.Primitive ((^), empty)+import Data.Manifold.PseudoAffine+import Data.Function.Differentiable+import Data.Function.Differentiable.Data+import Data.Manifold.TreeCover+    +import Data.Embedding+import Data.CoNat++import qualified Prelude as Hask hiding(foldl, sum, sequence)+import qualified Control.Applicative as Hask+import qualified Control.Monad       as Hask hiding(forM_, sequence)+import Data.Functor.Identity+import Control.Monad.Trans.State+import Control.Monad.Trans.Writer+import Control.Monad.Trans.Maybe+import Control.Monad.Trans.Class+import qualified Data.Foldable       as Hask+import Data.Foldable (all, elem, toList, sum, foldr1)+import qualified Data.Traversable as Hask+import Data.Traversable (forM)++import qualified Numeric.LinearAlgebra.HMatrix as HMat++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 Data.Foldable.Constrained+import Data.Traversable.Constrained (Traversable, traverse)++import GHC.Generics (Generic)+++type WebNodeId = Int+type NeighbourRefs = UArr.Vector WebNodeId++data PointsWeb :: * -> * -> * where+   PointsWeb :: {+       webNodeRsc :: ShadeTree x+     , webNodeAssocData :: Arr.Vector (y, NeighbourRefs)+     } -> PointsWeb x y+  deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable)++instance (NFData x, NFData (Needle' x), NFData y) => NFData (PointsWeb x y)++instance Foldable (PointsWeb x) (->) (->) where+  ffoldl = uncurry . Hask.foldl' . curry+  foldMap = Hask.foldMap+instance Traversable (PointsWeb x) (PointsWeb x) (->) (->) where+  traverse f (PointsWeb rsc asd)+           = fmap (PointsWeb rsc . (`Arr.zip`ngss) . Arr.fromList)+              . traverse f $ Arr.toList ys+   where (ys,ngss) = Arr.unzip asd++++fromWebNodes :: ∀ x y . WithField ℝ Manifold x+                    => (Shade x->Metric x) -> [(x,y)] -> PointsWeb x y+fromWebNodes mf = fromShaded mf . fromLeafPoints . map (uncurry WithAny . swap)++fromShadeTree_auto :: ∀ x . WithField ℝ Manifold x => ShadeTree x -> PointsWeb x ()+fromShadeTree_auto = fromShaded (recipMetric . _shadeExpanse) . constShaded ()++fromShadeTree :: ∀ x . WithField ℝ Manifold x+     => (Shade x -> Metric x) -> ShadeTree x -> PointsWeb x ()+fromShadeTree mf = fromShaded mf . constShaded ()++fromShaded :: ∀ x y . WithField ℝ Manifold x+     => (Shade x -> Metric x) -- ^ Local scalar-product generator. You can always+                              --   use @'recipMetric' . '_shadeExpanse'@ (but this+                              --   may give distortions compared to an actual+                              --   Riemannian metric).+     -> (x`Shaded`y)          -- ^ Source tree.+     -> PointsWeb x y+fromShaded metricf shd = PointsWeb shd' assocData + where shd' = stripShadedUntopological shd+       assocData = Hask.foldMap locMesh $ twigsWithEnvirons shd+       +       locMesh :: ((Int, ShadeTree (x`WithAny`y)), [(Int, ShadeTree (x`WithAny`y))])+                   -> Arr.Vector (y, NeighbourRefs)+       locMesh ((i₀, locT), neighRegions) = Arr.map findNeighbours locLeaves+        where locLeaves = Arr.map (first (+i₀)) . Arr.indexed . Arr.fromList+                                          $ onlyLeaves locT+              vicinityLeaves = Hask.foldMap+                                (\(i₀n, ngbR) -> Arr.map (first (+i₀n))+                                               . Arr.indexed+                                               . Arr.fromList+                                               $ onlyLeaves ngbR+                                ) neighRegions+              findNeighbours :: (Int, x`WithAny`y) -> (y, NeighbourRefs)+              findNeighbours (i, WithAny y x)+                         = (y, UArr.fromList $ fst<$>execState seek mempty)+               where seek = do+                        Hask.forM_ (locLeaves Arr.++ vicinityLeaves)+                                  $ \(iNgb, WithAny _ xNgb) ->+                           when (iNgb/=i) `id`do+                              let (Option (Just v)) = xNgb.-~.x+                              oldNgbs <- get+                              when (all (\(_,(_,nw)) -> visibleOverlap nw v) oldNgbs) `id`do+                                 let w = w₀ ^/ (w₀<.>^v)+                                      where w₀ = toDualWith locRieM v+                                 put $ (iNgb, (v,w))+                                       : [ neighbour+                                         | neighbour@(_,(nv,_))<-oldNgbs+                                         , visibleOverlap w nv+                                         ]+              +              visibleOverlap :: Needle' x -> Needle x -> Bool+              visibleOverlap w v = o < 1+               where o = w<.>^v+              +              locRieM :: Metric x+              locRieM = case pointsCovers . map _topological+                                  $ onlyLeaves locT+                                   ++ Hask.foldMap (onlyLeaves . snd) neighRegions of+                          [sh₀] -> metricf sh₀++indexWeb :: WithField ℝ Manifold x => PointsWeb x y -> WebNodeId -> Option (x,y)+indexWeb (PointsWeb rsc assocD) i+  | i>=0, i<Arr.length assocD+  , Right (_,x) <- indexShadeTree rsc i  = pure (x, fst (assocD Arr.! i))+  | otherwise                            = empty++webEdges :: ∀ x y . WithField ℝ Manifold x+            => PointsWeb x y -> [((x,y), (x,y))]+webEdges web@(PointsWeb rsc assoc) = (lookId***lookId) <$> toList allEdges+ where allEdges :: Set.Set (WebNodeId,WebNodeId)+       allEdges = Hask.foldMap (\(i,(_,ngbs))+                    -> Set.fromList [(min i i', max i i')+                                    | i'<-UArr.toList ngbs ]+                               ) $ Arr.indexed assoc+       lookId i | Option (Just xy) <- indexWeb web i  = xy+++localFocusWeb :: WithField ℝ Manifold x => PointsWeb x y -> PointsWeb x ((x,y), [(x,y)])+localFocusWeb (PointsWeb rsc asd) = PointsWeb rsc asd''+ where asd' = Arr.imap (\i (y,n) -> case indexShadeTree rsc i of+                                         Right (_,x) -> ((x,y),n) ) asd+       asd''= Arr.map (\(xy,n) ->+                       ((xy, [fst (asd' Arr.! j) | j<-UArr.toList n]), n)+                 ) asd'+++filterDEqnSolution_static :: (WithField ℝ Manifold x, Refinable y)+       => DifferentialEqn x y -> PointsWeb x (Shade' y) -> Option (PointsWeb x (Shade' y))+filterDEqnSolution_static f = localFocusWeb >>> Hask.traverse `id`+                   \((x,shy), ngbs) -> if null ngbs+                     then pure shy+                     else refineShade' shy+                            =<< filterDEqnSolution_loc f ((x,shy), NE.fromList ngbs)++

@@ -0,0 +1,37 @@++-- |+-- Module      : Data.SetLike.Intersection+-- Copyright   : (c) Justus Sagemüller 2016+-- License     : GPL v3+-- +-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de+-- Stability   : experimental+-- Portability : portable+-- +++module Data.SetLike.Intersection where++import Data.Semigroup+import qualified Data.List.NonEmpty as NE+import Data.List.NonEmpty (NonEmpty(..))+++newtype IntersectT s x = IntersectT { getIntersectors :: NonEmpty (s x) }+++singleIntersect :: s x -> IntersectT s x+singleIntersect = IntersectT . pure++rmTautologyIntersect ::+         (s x -> s x -> Option (s x)) -- ^ Subset-finder+      -> IntersectT s x -> IntersectT s x+rmTautologyIntersect smaller (IntersectT isoa) = IntersectT $ rti isoa+ where rti (s₀:|ss) = reduce [] ss+        where reduce [] [] = s₀:|[]+              reduce (sp₀:sp) [] = NE.cons s₀ $ rti (sp₀:|sp)+              reduce sp (s₁:sr) = case smaller s₀ s₁ of+               Option (Just si) -> rti $ si :| (sp ++ sr)+               Option Nothing   -> reduce (s₁:sp) sr+            +

@@ -86,6 +86,12 @@   completeBasis :: Tagged v [Basis v]   completeBasis = liftA2 (\dim f -> f <$> [0 .. dim - 1]) dimension indexBasis   +  completeBasisValues :: [v]+  completeBasisValues = defCBVs+   where defCBVs :: ∀ v . FiniteDimensional v => [v]+         defCBVs = basisValue <$> cb+          where Tagged cb = completeBasis :: Tagged v [Basis v]+     asPackedVector :: v -> HMat.Vector (Scalar v)   asPackedVector v = HMat.fromList $ snd <$> decompose v   

binary file changed (absent → 84463 bytes)

binary file changed (absent → 36607 bytes)

binary file changed (absent → 203879 bytes)

binary file changed (absent → 231921 bytes)

binary file changed (absent → 101994 bytes)

@@ -1,5 +1,5 @@ Name:                manifolds-Version:             0.2.0.1+Version:             0.2.2.0 Category:            Math Synopsis:            Coordinate-free hypersurfaces Description:         Manifolds, a generalisation of the notion of “smooth curves” or surfaces,@@ -29,7 +29,8 @@ Maintainer:          (@) sagemueller $ geo.uni-koeln.de Build-Type:          Simple Cabal-Version:       >=1.10-Extra-Doc-Files:     images/examples/*.png+Extra-Doc-Files:     images/examples/*.png,+                     images/examples/ShadeCombinations/2Dconvolution-skewed.png  Source-Repository head     type: git@@ -64,6 +65,7 @@   Exposed-modules:   Data.Manifold                      Data.Manifold.PseudoAffine                      Data.Manifold.TreeCover+                     Data.Manifold.Web                      Data.SimplicialComplex                      Data.LinearMap.HerMetric                      Data.Function.Differentiable@@ -72,6 +74,7 @@                      Data.Manifold.Riemannian   Other-modules:   Data.List.FastNub                    Data.Manifold.Types.Primitive+                   Data.SetLike.Intersection                    Data.Manifold.Cone                    Data.CoNat                    Data.Embedding