manifolds 0.2.2.0 → 0.2.3.0

raw patch · 20 files changed

+1039/−260 lines, 20 filesdep +microlensdep +microlens-thPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: microlens, microlens-th

API changes (from Hackage documentation)

Files

• Control/Monad/Trans/OuterMaybe.hs +22/−0
• Data/CoNat.hs +0/−2
• Data/Embedding.hs +0/−3
• Data/Function/Affine.hs +0/−9
• Data/Function/Differentiable.hs +0/−12
• Data/LinearMap/Category.hs +10/−3
• Data/LinearMap/HerMetric.hs +68/−19
• Data/List/FastNub.hs +9/−0
• Data/Manifold/Cone.hs +0/−11
• Data/Manifold/DifferentialEquation.hs +117/−0
• Data/Manifold/Griddable.hs +0/−24
• Data/Manifold/PseudoAffine.hs +45/−13
• Data/Manifold/Riemannian.hs +1/−22
• Data/Manifold/TreeCover.hs +241/−65
• Data/Manifold/Types.hs +19/−4
• Data/Manifold/Types/Primitive.hs +7/−8
• Data/Manifold/Web.hs +496/−55
• Data/SimplicialComplex.hs +0/−8
• Util/LtdShow.hs +0/−1
• manifolds.cabal +4/−1

@@ -0,0 +1,22 @@+-- |+-- Module      : Control.Monad.Trans.OuterMaybe+-- Copyright   : (c) Justus Sagemüller 2016+-- License     : GPL v3+-- +-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de+-- Stability   : experimental+-- Portability : portable+-- +{-# LANGUAGE DeriveFunctor              #-}++module Control.Monad.Trans.OuterMaybe where++data OuterMaybeT f a = OuterNothing | OuterJust (f a) deriving (Functor)+instance (Applicative f) => Applicative (OuterMaybeT f) where+  pure = OuterJust . pure+  OuterJust fs <*> OuterJust xs = OuterJust $ fs <*> xs+  _ <*> _ = OuterNothing++++

@@ -45,7 +45,6 @@ import Data.VectorSpace import Data.AffineSpace import Data.Basis-import Data.AdditiveGroup import qualified Data.List as List      import qualified Prelude as Hask hiding(foldl)@@ -60,7 +59,6 @@   import qualified Data.Vector as Arr-import qualified Numeric.LinearAlgebra.HMatrix as HMat  import Unsafe.Coerce 

@@ -29,9 +29,6 @@  module Data.Embedding where -import Data.Tagged-import Data.Semigroup- import qualified Prelude as Hask hiding(foldl) import qualified Control.Applicative as Hask import qualified Control.Monad       as Hask

@@ -37,28 +37,19 @@       -import Data.List-import Data.Maybe import Data.Semigroup  import Data.VectorSpace import Data.LinearMap import Data.LinearMap.HerMetric-import Data.MemoTrie (HasTrie(..)) import Data.AffineSpace-import Data.Basis-import Data.Void import Data.Tagged import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine -import Data.CoNat-import Data.VectorSpace.FiniteDimensional- import qualified Prelude import qualified Control.Applicative as Hask -import Data.Constraint.Trivial import Control.Category.Constrained.Prelude hiding ((^)) import Control.Category.Constrained.Reified import Control.Arrow.Constrained

@@ -50,33 +50,21 @@   import Data.List-import qualified Data.Vector.Generic as Arr-import qualified Data.Vector import Data.Maybe import Data.Semigroup-import Data.Function (on) import Data.Embedding-import Data.Fixed  import Data.VectorSpace import Data.LinearMap import Data.LinearMap.Category import Data.LinearMap.HerMetric-import Data.MemoTrie (HasTrie(..)) import Data.AffineSpace import Data.Function.Differentiable.Data import Data.Function.Affine import Data.Basis-import Data.Complex hiding (magnitude)-import Data.Void import Data.Tagged import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine--import Data.CoNat-import Data.VectorSpace.FiniteDimensional--import qualified Numeric.LinearAlgebra.HMatrix as HMat  import qualified Prelude import qualified Control.Applicative as Hask

@@ -31,15 +31,12 @@ module Data.LinearMap.Category where  import Data.Tagged-import Data.Semigroup -import Data.MemoTrie import Data.VectorSpace import Data.LinearMap import Data.VectorSpace.FiniteDimensional import Data.AffineSpace import Data.Basis-import Data.AdditiveGroup      import qualified Prelude as Hask hiding(foldl) import qualified Control.Applicative as Hask@@ -139,6 +136,16 @@  instance (SmoothScalar s) => EnhancedCat (->) (Linear s) where   arr (DenseLinear mat) = fromPackedVector . HMat.app mat . asPackedVector++-- | Inverse function application (for isomorphisms), or+--   least-square solution of a linear equation.+--   Note that least-square is not really well-defined,+--   without reference to a norm / scalar product; the operator uses+--   the implicit norm induced from the 'FiniteDimensional' representation.+(<\$) :: ( SmoothScalar s, FiniteDimensional v, FiniteDimensional w+         , Scalar v ~ s, Scalar w ~ s+         ) => Linear s v w -> w -> v+DenseLinear mat <\$ v = fromPackedVector . (mat HMat.<\>) $ asPackedVector v  type DenseLinearFuncValue s = GenericAgent (Linear s) 

@@ -30,17 +30,19 @@   -- * One-dimensional axes and product spaces   , factoriseMetric, factoriseMetric'   , productMetric, productMetric'-  , metricAsLength, metricFromLength, metric'AsLength+  , tryMetricAsLength, metricAsLength, metricFromLength, metric'AsLength   -- * Utility for metrics   , transformMetric, transformMetric', dualCoCoProduct   , dualiseMetric, dualiseMetric'-  , recipMetric, recipMetric'+  , recipMetric, recipMetric', safeRecipMetric, safeRecipMetric'   -- ** Eigenvectors   , eigenSpan, eigenSpan'   , eigenCoSpan, eigenCoSpan'   , eigenSystem, HasEigenSystem, EigenVector+  -- ** Scaling operations   , metriNormalise, metriNormalise'   , metriScale', metriScale+  , volumeRatio, euclideanRelativeMetricVolume   , adjoint   , extendMetric   , applyLinMapMetric, applyLinMapMetric'@@ -58,6 +60,7 @@   , linMapAsTensProd, linMapFromTensProd   , covariance   , outerProducts+  , orthogonalComplementSpan   ) where      @@ -66,10 +69,8 @@ import Data.VectorSpace import Data.LinearMap import Data.Basis-import Data.MemoTrie import Data.Semigroup import Data.Tagged-import Data.Void import qualified Data.List as List  import qualified Prelude as Hask@@ -297,20 +298,27 @@ --   a space to its dual, the inverse maps from the dual into the --   (double-dual) space &#x2013; i.e., it is a metric on the dual space. --   Deprecated: the singular case isn't properly handled.+recipMetric :: HasMetric v => HerMetric' v -> HerMetric v+recipMetric m' | Option (Just m) <- safeRecipMetric m'  = m+recipMetric _ = singularMetric+ recipMetric' :: HasMetric v => HerMetric v -> HerMetric' v-recipMetric' (HerMetric Nothing) = singularMetric'-recipMetric' (HerMetric (Just (DenseLinear m)))-          | isInfinite' detm  = singularMetric'-          | otherwise         = matrixMetric' minv- where (minv, (detm, _)) = HMat.invlndet m+recipMetric' m | Option (Just m') <- safeRecipMetric' m  = m'+recipMetric' _ = singularMetric' -recipMetric :: HasMetric v => HerMetric' v -> HerMetric v-recipMetric (HerMetric' Nothing) = singularMetric-recipMetric (HerMetric' (Just (DenseLinear m)))-          | isInfinite' detm  = singularMetric-          | otherwise         = matrixMetric minv+safeRecipMetric :: HasMetric v => HerMetric' v -> Option (HerMetric v)+safeRecipMetric (HerMetric' Nothing) = empty+safeRecipMetric (HerMetric' (Just (DenseLinear m)))+          | isInfinite' detm  = empty+          | otherwise         = return $ matrixMetric minv  where (minv, (detm, _)) = HMat.invlndet m +safeRecipMetric' :: HasMetric v => HerMetric v -> Option (HerMetric' v)+safeRecipMetric' (HerMetric Nothing) = empty+safeRecipMetric' (HerMetric (Just (DenseLinear m)))+          | isInfinite' detm  = empty+          | otherwise         = return $ matrixMetric' minv+ where (minv, (detm, _)) = HMat.invlndet m  isInfinite' :: (Eq a, Num a) => a -> Bool isInfinite' 0 = False@@ -449,7 +457,7 @@          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, [])+  eigenSystem (_, _) = (fmap Stiefel1 completeBasisValues, [])   -- | Constraint that a space's scalars need to fulfill so it can be used for 'HerMetric'.@@ -728,11 +736,29 @@            = HMat.invlndet . getDenseMatrix $ fst . m . (id&&&zeroV)  +volumeRatio :: HasMetric v => HerMetric v -> HerMetric v -> Scalar v+volumeRatio (HerMetric Nothing) (HerMetric Nothing) = 1+volumeRatio (HerMetric _) (HerMetric Nothing) = 0+volumeRatio (HerMetric (Just (DenseLinear m₁)))+            (HerMetric (Just (DenseLinear m₂)))+    = HMat.det m₂ / HMat.det m₁+volumeRatio (HerMetric Nothing) (HerMetric _) = 1/0++euclideanRelativeMetricVolume :: (HasMetric v, InnerSpace v) => HerMetric v -> Scalar v+euclideanRelativeMetricVolume (HerMetric Nothing) = 1/0+euclideanRelativeMetricVolume (HerMetric (Just (DenseLinear m))) = recip $ HMat.det m++tryMetricAsLength :: HerMetric ℝ -> Option ℝ+tryMetricAsLength m = case metricSq m 1 of+   o | o > 0      -> pure . sqrt $ recip o+     | otherwise  -> empty++-- | Unsafe version of 'tryMetricAsLength', only works reliable if the metric+--   is strictly positive definite. metricAsLength :: HerMetric ℝ -> ℝ metricAsLength m = case metricSq m 1 of-   o | o > 0      -> sqrt $ recip o+   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 ℝ@@ -765,13 +791,27 @@  -- | Same as 'spanHilbertSubspace', but with the standard 'euclideanMetric' (i.e., the --   basis vectors will be orthonormal in the usual sense, in both @w@ and @v@).-spanSubHilbertSpace :: forall s v w+spanSubHilbertSpace :: ∀ s v w         . (HasMetric v, InnerSpace v, Scalar v ~ s, IsFreeSpace w, Scalar w ~ s)       => [v]           -> Option (Embedding (Linear s) w v) spanSubHilbertSpace = spanHilbertSubspace euclideanMetric'  +orthogonalComplementSpan :: ∀ v . (HasMetric v, Scalar v ~ ℝ)+                            => [Stiefel1 (DualSpace v)] -> [Stiefel1 v]+orthogonalComplementSpan avoidSpace+           = fst ( iterate nextOVect ( [], ( cycle completeBasisValues+                                           , pseudoRieszPair <$> avoidSpace ) )+                    !! (d - lav) )+ where Tagged d = dimension :: Tagged v Int+       lav = length avoidSpace+       nextOVect (result, (v:src, avoid))+           | Option (Just newAvoid@(vfin', _)) <- mkPseudoRieszPair vPurged+                          = (Stiefel1 vfin':result, (src, newAvoid : avoid))+        where vPurged = foldl (\vp (av', av) -> vp ^-^ av ^* (vp^<.>av')) v avoid++ -- | The /n/-th Stiefel manifold is the space of all possible configurations of --   /n/ orthonormal vectors. In the case /n/ = 1, simply the subspace of normalised --   vectors, i.e. equivalent to the 'UnitSphere'. Even so, it strictly speaking@@ -781,8 +821,17 @@ --   vectors modulo scaling by positive factors. newtype Stiefel1 v = Stiefel1 { getStiefel1N :: DualSpace v } -+pseudoRieszPair :: (HasMetric v, Scalar v ~ ℝ) => Stiefel1 v -> (v, DualSpace v)+pseudoRieszPair (Stiefel1 v')+              = (fromPackedVector $ HMat.scale (1/HMat.norm_2 vp) vp, v')+ where vp = asPackedVector v' +mkPseudoRieszPair :: (HasMetric v, Scalar v ~ ℝ) => DualSpace v -> Option (v, DualSpace v)+mkPseudoRieszPair v'+   | nv' > 0    = pure (fromPackedVector $ HMat.scale (1/nv') vp, v')+   | otherwise  = empty+ where vp = asPackedVector v'+       nv' = HMat.norm_2 vp   

@@ -6,6 +6,8 @@  import Data.List import Data.Function+import Data.Ord+import Control.Arrow ((&&&))   type FastNub a = (Eq a, Ord a) -- S̶h̶o̶u̶l̶d̶ ̶r̶e̶a̶l̶l̶y̶ ̶b̶e̶ ̶(̶E̶q̶ ̶a̶,̶̶ ̶H̶a̶s̶h̶a̶b̶l̶e̶ ̶a̶)̶@@ -63,3 +65,10 @@        fis (x:xs) (y:ys) | x<y  = fis xs (y:ys)                          | x>y  = fis (x:xs) ys                          | otherwise  = x : fis xs ys+++-- | This function is also defined in "GHC.Exts", but only in a version that requires+--   𝓞(𝑛⋅log 𝑛) function applications, as opposed to 𝑛 here.+sortWith :: Ord b => (a -> b) -> [a] -> [a]+sortWith f = map snd . sortBy (comparing fst) . map (f &&& id)+

@@ -31,22 +31,12 @@       -import Data.List import qualified Data.Vector.Generic as Arr-import qualified Data.Vector import Data.Maybe import Data.Semigroup-import Data.Function (on)-import Data.Fixed  import Data.VectorSpace-import Data.LinearMap import Data.LinearMap.HerMetric-import Data.MemoTrie (HasTrie(..))-import Data.AffineSpace-import Data.Basis-import Data.Complex hiding (magnitude)-import Data.Void import Data.Tagged import Data.Manifold.Types.Primitive @@ -64,7 +54,6 @@ import Data.Foldable.Constrained  import Data.Manifold.PseudoAffine-import Data.Embedding   

@@ -0,0 +1,117 @@+-- |+-- Module      : Data.Manifold.DifferentialEquation+-- 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 MultiParamTypeClasses      #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE GADTs                      #-}+{-# LANGUAGE RankNTypes                 #-}+{-# LANGUAGE TupleSections              #-}+{-# LANGUAGE ParallelListComp           #-}+{-# LANGUAGE UnicodeSyntax              #-}+{-# LANGUAGE ConstraintKinds            #-}+{-# LANGUAGE PatternGuards              #-}+{-# LANGUAGE LambdaCase                 #-}+{-# LANGUAGE TypeOperators              #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE LiberalTypeSynonyms        #-}+++module Data.Manifold.DifferentialEquation (+            -- * Formulating simple differential eqns.+              DifferentialEqn+            , constLinearDEqn+            , filterDEqnSolution_static, iterateFilterDEqn_static+            -- * Cost functions for error bounds+            , maxDeviationsGoal+            , uncertaintyGoal+            , uncertaintyGoal'+            , euclideanVolGoal+            ) where+++import Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NE+import Data.Semigroup++import Data.VectorSpace+import Data.LinearMap.HerMetric+import Data.LinearMap.Category+import Data.AffineSpace+import Data.Basis++import Data.Manifold.Types+import Data.Manifold.PseudoAffine+import Data.Function.Differentiable+import Data.Function.Differentiable.Data+import Data.Manifold.TreeCover+import Data.Manifold.Web++import qualified Numeric.LinearAlgebra.HMatrix as HMat+import qualified Data.List as List++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 qualified Data.Foldable       as Hask+import qualified Data.Traversable as Hask++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)+++constLinearDEqn :: (WithField ℝ LinearManifold x, WithField ℝ LinearManifold y)+              => Linear ℝ (DualSpace y) (Linear ℝ y x) -> DifferentialEqn x y+constLinearDEqn bwt = factoriseShade+    >>> \(_x, Shade y δy) -> let j = bwt'm HMat.<\> (asPackedVector y)+                                 δj = bwt' `transformMetric` recipMetric δy+                             in Shade' (fromPackedVector j) δj+ where bwt'@(DenseLinear bwt'm) = adjoint bwt+++-- | A function that variates, relatively speaking, most strongly+--   for arguments around 1. In the zero-limit it approaches a constant+--   (but with arbitrarily large derivative); for η → ∞ the derivative+--   approaches 0.+--   +--   The idea is that if you consider the ratio of two function values,+--   it will be close to 1 if both arguments on the same side of 1,+--   even if their ratio is large.+--   Only if both arguments are close to 1, or lie on opposite sides+--   of it, will the ratio of the function values will be significant.+goalSensitive :: ℝ -> ℝ+goalSensitive η =  0.3 + sqrt (η * (1 + η/(1+η)) / (3 + η))++euclideanVolGoal :: WithField ℝ EuclidSpace y => ℝ -> x -> Shade' y -> ℝ+euclideanVolGoal vTgt _ (Shade' _ shy) = goalSensitive η+ where η = euclideanRelativeMetricVolume shy / vTgt++maxDeviationsGoal :: WithField ℝ EuclidSpace y => [Needle y] -> x -> Shade' y -> ℝ+maxDeviationsGoal = uncertaintyGoal . projector's++uncertaintyGoal :: WithField ℝ EuclidSpace y => Metric' y -> x -> Shade' y -> ℝ+uncertaintyGoal = uncertaintyGoal' . const++uncertaintyGoal' :: WithField ℝ EuclidSpace y => (x -> Metric' y) -> x -> Shade' y -> ℝ+uncertaintyGoal' f x (Shade' _ shy)+         = List.sum [goalSensitive $ 1 / metricSq' m q | q <- shySpan]+ where shySpan = eigenSpan' shy+       m = f x

@@ -37,47 +37,24 @@  import Data.List hiding (filter, all, elem, sum) import Data.Maybe-import qualified Data.Map as Map-import qualified Data.Vector as Arr-import Data.List.NonEmpty (NonEmpty(..))-import Data.List.FastNub-import qualified Data.List.NonEmpty as NE-import Data.Semigroup -import Data.VectorSpace-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 ((^), (^.)) import Data.Manifold.PseudoAffine import Data.Manifold.TreeCover (Shade(..), fullShade, shadeCtr, shadeExpanse)      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.Class import qualified Data.Foldable       as Hask import Data.Foldable (all, elem, toList, sum) 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(..), Traversable) import Control.Arrow.Constrained@@ -85,7 +62,6 @@ import Data.Foldable.Constrained  import Text.Printf-import GHC.Generics (Generic)   data GridAxis m g = GridAxInterval (Shade m)

@@ -65,17 +65,14 @@             -- ** Local functions             , LocalLinear, LocalAffine             -- * Misc-            , palerp, LocallyCoercible(..)+            , alerpB, palerp, palerpB, LocallyCoercible(..)+            , ImpliesMetric(..)             ) where       -import Data.List-import qualified Data.Vector.Generic as Arr-import qualified Data.Vector import Data.Maybe import Data.Semigroup-import Data.Function (on) import Data.Fixed  import Data.VectorSpace@@ -83,19 +80,13 @@ import Data.LinearMap import Data.LinearMap.HerMetric import Data.LinearMap.Category-import Data.MemoTrie (HasTrie(..)) import Data.AffineSpace-import Data.Basis-import Data.Complex hiding (magnitude)-import Data.Void import Data.Tagged import Data.Manifold.Types.Primitive  import Data.CoNat import Data.VectorSpace.FiniteDimensional -import qualified Numeric.LinearAlgebra.HMatrix as HMat- import qualified Prelude import qualified Control.Applicative as Hask @@ -104,6 +95,7 @@ import Control.Monad.Constrained import Data.Foldable.Constrained +import GHC.Exts (Constraint)   @@ -191,6 +183,7 @@ class ( Semimanifold x, Semimanifold (Interior x)       , Needle (Interior x) ~ Needle x, Interior (Interior x) ~ Interior x)         => PseudoAffine x where+  {-# MINIMAL (.-~.) | (.-~!) #-}   -- | The path reaching from one point to another.   --   Should only yield 'Nothing' if   -- @@ -212,6 +205,14 @@   --   manifold”. To adress this problem, these functions basically consider only the   --   /interior/ of the space.   (.-~.) :: x -> Interior x -> Option (Needle x)+  p.-~.q = return $ p.-~!q+  +  -- | Unsafe version of '.-~.'. If the two points lie in disjoint regions,+  --   the behaviour is undefined.+  (.-~!) :: x -> Interior x -> Needle x+  p.-~!q = case p.-~.q of+      Option (Just v) -> v+         @@ -320,17 +321,29 @@ --   its end points. --  --   A proper, really well-defined (on global scales) interpolation---   only makes sense on a Riemannian manifold, as geodesics.---   This is a task to be tackled in the future.+--   only makes sense on a Riemannian manifold, as 'Data.Manifold.Riemannian.Geodesic'. palerp :: ∀ x. Manifold x     => Interior x -> Interior x -> Option (Scalar (Needle x) -> x) palerp p1 p2 = case (fromInterior p2 :: x) .-~. p1 of   Option (Just v) -> return $ \t -> p1 .+~^ t *^ v   _ -> empty +-- | Like 'palerp', but actually restricted to the interval between the points,+--   with a signature like 'Data.Manifold.Riemannian.geodesicBetween'+--   rather than 'Data.AffineSpace.alerp'.+palerpB :: ∀ x. WithField ℝ Manifold x => Interior x -> Interior x -> Option (D¹ -> x)+palerpB p1 p2 = case (fromInterior p2 :: x) .-~. p1 of+  Option (Just v) -> return $ \(D¹ t) -> p1 .+~^ ((t+1)/2) *^ v+  _ -> empty +-- | Like 'alerp', but actually restricted to the interval between the points.+alerpB :: ∀ x. (AffineSpace x, VectorSpace (Diff x), Scalar (Diff x) ~ ℝ)+                   => x -> x -> D¹ -> x+alerpB p1 p2 = case p2 .-. p1 of+  v -> \(D¹ t) -> p1 .+^ ((t+1)/2) *^ v  + hugeℝVal :: ℝ hugeℝVal = 1e+100 @@ -582,5 +595,24 @@   ++class ImpliesMetric s where+  {-# MINIMAL inferMetric | inferMetric' #-}+  type MetricRequirement s x :: Constraint+  type MetricRequirement s x = Semimanifold x+  inferMetric :: (MetricRequirement s x, HasMetric (Needle x))+                     => s x -> Option (Metric x)+  inferMetric = safeRecipMetric <=< inferMetric'+  inferMetric' :: (MetricRequirement s x, HasMetric (Needle x))+                     => s x -> Option (Metric' x)+  inferMetric' = safeRecipMetric' <=< inferMetric++instance ImpliesMetric HerMetric where+  type MetricRequirement HerMetric x = x ~ Needle x+  inferMetric = pure++instance ImpliesMetric HerMetric' where+  type MetricRequirement HerMetric' x = x ~ Needle x+  inferMetric' = pure  

@@ -46,47 +46,27 @@ module Data.Manifold.Riemannian  where  -import Data.List hiding (filter, all, elem, sum) import Data.Maybe-import qualified Data.Map as Map import qualified Data.Vector as Arr-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.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.Manifold.Types import Data.Manifold.Types.Primitive ((^), empty, embed, coEmbed) import Data.Manifold.PseudoAffine import Data.VectorSpace.FiniteDimensional     -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.Class import qualified Data.Foldable       as Hask-import Data.Foldable (all, elem, toList, sum) import qualified Data.Traversable as Hask-import Data.Traversable (forM)  import qualified Numeric.LinearAlgebra.HMatrix as HMat @@ -96,10 +76,9 @@ import Control.Monad.Constrained hiding (forM) import Data.Foldable.Constrained -import GHC.Generics (Generic)  -class PseudoAffine x => Geodesic x where+class Semimanifold x => Geodesic x where   geodesicBetween ::           x -- ^ Starting point; the interpolation will yield this at -1.        -> x -- ^ End point, for +1.

@@ -46,9 +46,10 @@        -- ** Evaluation        , occlusion        -- ** Misc-       , factoriseShade, intersectShade's, Refinable, refineShade', convolveShade', coerceShade+       , factoriseShade, intersectShade's+       , Refinable, subShade', refineShade', convolveShade', coerceShade        -- * Shade trees-       , ShadeTree(..), fromLeafPoints, onlyLeaves, indexShadeTree+       , ShadeTree(..), fromLeafPoints, onlyLeaves, indexShadeTree, positionIndex        -- * View helpers        , onlyNodes        -- ** Auxiliary types@@ -56,9 +57,9 @@        -- * Misc        , sShSaw, chainsaw, HasFlatView(..), shadesMerge, smoothInterpolate        , twigsWithEnvirons, completeTopShading, flexTwigsShading-       , WithAny(..), Shaded, stiAsIntervalMapping, spanShading+       , WithAny(..), Shaded, fmapShaded, stiAsIntervalMapping, spanShading        , constShaded, stripShadedUntopological-       , DifferentialEqn, filterDEqnSolution_loc+       , DifferentialEqn, propagateDEqnSolution_loc        -- ** Triangulation-builders        , TriangBuild, doTriangBuild, singleFullSimplex, autoglueTriangulation        , AutoTriang, elementaryTriang, breakdownAutoTriang@@ -78,34 +79,28 @@  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.SetLike.Intersection+import Data.Manifold.Riemannian      import Data.Embedding import Data.CoNat +import Lens.Micro (Lens')+ 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.OuterMaybe import Control.Monad.Trans.Class import qualified Data.Foldable       as Hask import Data.Foldable (all, elem, toList, sum, foldr1)@@ -153,7 +148,7 @@ class IsShade shade where --  type (*) shade :: *->*   -- | Access the center of a 'Shade' or a 'Shade''.-  shadeCtr :: Functor f (->) (->) => (Interior x->f (Interior x)) -> shade x -> f (shade x)+  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.@@ -179,7 +174,15 @@   coerceShade (Shade x (HerMetric' δxym))           = Shade (locallyTrivialDiffeomorphism x) (HerMetric' $ unsafeCoerceLinear<$>δxym) -shadeExpanse :: Functor f (->) (->) => (Metric' x -> f (Metric' x)) -> Shade x -> f (Shade x)+instance ImpliesMetric Shade where+  type MetricRequirement Shade x = Manifold x+  inferMetric' (Shade _ e) = pure e++instance ImpliesMetric Shade' where+  type MetricRequirement Shade' x = Manifold x+  inferMetric (Shade' _ e) = pure e++shadeExpanse :: Lens' (Shade x) (Metric' x) shadeExpanse f (Shade c e) = fmap (Shade c) $ f e  instance IsShade Shade' where@@ -195,7 +198,7 @@   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 :: Lens' (Shade' x) (Metric x) shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e  instance (AffineManifold x) => Semimanifold (Shade x) where@@ -206,6 +209,32 @@   Shade c e .+~^ v = Shade (c.+^v) e   Shade c e .-~^ v = Shade (c.-^v) e +instance (WithField ℝ AffineManifold x, Geodesic x) => Geodesic (Shade x) where+  geodesicBetween (Shade c e) (Shade ζ η) = pure interp+   where ([], sharedSpan) = eigenSystem (e,η)+         interp t = Shade (pinterp t)+                          (projector's [ v ^* (alerpB qe qη t)+                                       | ([qe,qη], (v,_)) <- zip coeffs sharedSpan ])+         coeffs = [ [metric' m v' | m <- [e,η]] | (_,v') <- sharedSpan ]+         Option (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) => Geodesic (Shade' x) where+  geodesicBetween (Shade' c e) (Shade' ζ η) = pure interp+   where ([], sharedSpan) = eigenSystem (e,η)+         interp t = Shade' (pinterp t)+                           (projectors [ v' ^/ (alerpB qe qη t)+                                       | ([qe,qη], (v',_)) <- zip coeffs sharedSpan ])+         coeffs = [ [recip $ metric m v | m <- [e,η]] | (_,v) <- sharedSpan ]+         Option (Just pinterp) = geodesicBetween c ζ+ fullShade :: WithField ℝ Manifold x => x -> Metric' x -> Shade x fullShade ctr expa = Shade ctr expa @@ -433,6 +462,15 @@        scanLeafNums i₀ ((v,t):vts) = (i₀, (v,t)) : scanLeafNums (i₀ + nLeaves t) vts  +indexDBranches :: NonEmpty (DBranch x) -> NonEmpty (DBranch' x (Int, ShadeTree x))+indexDBranches (DBranch d (Hourglass t b) :| l) -- this could more concisely be written as a traversal+              = DBranch d (Hourglass (0,t) (nt,b)) :| ixDBs (nt + nb) l+ where nt = nLeaves t; nb = nLeaves b+       ixDBs _ [] = []+       ixDBs i₀ (DBranch δ (Hourglass τ β) : l)+               = DBranch δ (Hourglass (i₀,τ) (i₀+nτ,β)) : ixDBs (i₀ + nτ + nβ) l+        where nτ = nLeaves τ; nβ = nLeaves β+ instance (NFData x, NFData (Needle' x)) => NFData (ShadeTree x) where   rnf (PlainLeaves xs) = rnf xs   rnf (DisjointBranches n bs) = n `seq` rnf (NE.toList bs)@@ -477,6 +515,8 @@ fromLeafPoints = fromLeafPoints' sShIdPartition  +-- | The leaves of a shade tree are numbered. For a given index, this function+--   attempts to find the leaf with that ID, within its immediate environment. indexShadeTree :: ∀ x . WithField ℝ Manifold x        => ShadeTree x -> Int -> Either Int ([ShadeTree x], x) indexShadeTree _ i@@ -498,6 +538,44 @@     | otherwise  = Left $ i-n  +-- | “Inverse indexing” of a tree. This is roughly a nearest-neighbour search,+--   but not guaranteed to give the correct result unless evaluated at the+--   precise position of a tree leaf.+positionIndex :: ∀ x . WithField ℝ Manifold x+       => Option (Metric x)  -- ^ For deciding (at the lowest level) what “close” means;+                             --   this is optional for any tree of depth >1.+        -> ShadeTree x       -- ^ The tree to index into+        -> x                 -- ^ Position to look up+        -> Option (Int, ([ShadeTree x], x))+                   -- ^ Index of the leaf near to the query point, the “path” of+                   --   environment trees leading down to its position (in decreasing+                   --   order of size), and actual position of the found node.+positionIndex (Option (Just m)) sh@(PlainLeaves lvs) x+        = case catMaybes [ ((i,p),) . metricSq m <$> getOption (p.-~.x)+                            | (i,p) <- zip [0..] lvs] of+           [] -> empty+           l | ((i,p),_) <- minimumBy (comparing snd) l+              -> pure (i, ([sh], p))+positionIndex m (DisjointBranches _ brs) x+        = fst . foldl' (\case+                          (q@(Option (Just _)), i₀) -> const (q, i₀)+                          (_, i₀) -> \t' -> ( first (+i₀) <$> positionIndex m t' x+                                            , i₀+nLeaves t' ) )+                       (empty, 0)+              $        brs+positionIndex _ sh@(OverlappingBranches n (Shade c ce) brs) x+   | Option (Just vx) <- x.-~.c+        = let (_,(i₀,t')) = maximumBy (comparing fst)+                       [ (σ*ω, t')+                       | DBranch d (Hourglass t'u t'd) <- NE.toList $ indexDBranches brs+                       , let ω = d<.>^vx+                       , (t',σ) <- [(t'u, 1), (t'd, -1)] ]+          in ((+i₀) *** first (sh:))+                 <$> positionIndex (return $ recipMetric ce) t' x+positionIndex _ _ _ = empty+++ fromFnGraphPoints :: ∀ x y . (WithField ℝ Manifold x, WithField ℝ Manifold y)                      => [(x,y)] -> ShadeTree (x,y) fromFnGraphPoints = fromLeafPoints' fg_sShIdPart@@ -602,6 +680,17 @@ nLeaves (DisjointBranches n _) = n nLeaves (OverlappingBranches n _ _) = n ++instance ImpliesMetric ShadeTree where+  type MetricRequirement ShadeTree x = WithField ℝ Manifold x+  inferMetric' (OverlappingBranches _ (Shade _ e) _) = pure e+  inferMetric' (PlainLeaves lvs) = case pointsShades lvs of+        (Shade _ sh:_) -> pure sh+        _ -> empty+  inferMetric' (DisjointBranches _ (br:|_)) = inferMetric' br+++ overlappingBranches :: Shade x -> NonEmpty (DBranch x) -> ShadeTree x overlappingBranches shx brs = OverlappingBranches n shx brs  where n = sum $ fmap (sum . fmap nLeaves) brs@@ -624,7 +713,7 @@       in DisjointBranches (sum $ nLeaves<$>brs') brs' unsafeFmapTree f fn fs (OverlappingBranches n sh brs)     = let brs' = fmap (\(DBranch dir br)-                        -> DBranch (fn dir) (unsafeFmapTree f fn fs<$>br)+                      -> DBranch (fn dir) (unsafeFmapTree f fn fs<$>br)                       ) brs       in overlappingBranches (fs sh) brs' @@ -636,37 +725,125 @@   --   @'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 ]+                                  [ ac.+~^σ*^v | σ<-[-1,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-                                     [ (v, 1 / (1 + metricSq oExpa w))-                                     | v <- (*^) <$> [-1,1] <*> span-                                     , let p = ctr .+~^ v  :: y-                                           Option (Just w) = p.-~.oCtr-                                     ]-                         | ctr                  <- [c,     ζ    ]-                         | span <- eigenCoSpan'<$> [e,     η    ]-                         | (oCtr,oExpa)         <- [(ζ,η), (c,e)]-                         ]-              Option (Just c'2ζ') = ζ'.-~.c'-              Option (Just c2ζ') = ζ'.-~.c-              Option (Just ζ2c') = c'.-~.ζ-              μc = metricSq e c2ζ'-              μζ = metricSq η ζ2c'-              iCtr = c' .+~^ c'2ζ' ^* (μζ/(μc + μζ)) -- weighted mean between c' and ζ'.-              Option (Just rc) = c.-~.iCtr-              Option (Just rζ) = ζ.-~.iCtr-              rcⰰ = toDualWith e rc-              rζⰰ = toDualWith η rζ-              μe = rcⰰ<.>^rc-              μη = rζⰰ<.>^rζ-              iExpa = (e^+^η)^/2 ^+^ projector rcⰰ^/(1-μe) ^+^ projector rζⰰ^/(1-μη)+  refineShade' (Shade' c₀ (HerMetric (Just e₁))) +               (Shade' c₀₂ (HerMetric (Just e₂)))+           | Option (Just c₂) <- c₀₂.-~.c₀+           , e₁c₂ <- e₁ $ c₂+           , e₂c₂ <- e₂ $ c₂+           , cc <- σe <\$ e₂c₂+           , cc₂ <- cc ^-^ c₂+           , e₁cc <- e₁ $ cc+           , e₂cc <- e₂ $ cc+           , α <- 2 + cc₂<.>^e₂c₂+           , α > 0+           , ee <- σe ^/ α+           , c₂e₁c₂ <- c₂^<.>e₁c₂+           , c₂e₂c₂ <- c₂^<.>e₂c₂+           , c₂eec₂ <- (c₂e₁c₂ + c₂e₂c₂) / α+           , [γ₁,γ₂] <- middle . sort+                $ quadraticEqnSol c₂e₁c₂+                                  (2 * (c₂^<.>e₁cc))+                                  (cc^<.>e₁cc - 1)+               ++ quadraticEqnSol c₂e₂c₂+                                  (2 * (c₂^<.>e₂cc - c₂e₂c₂))+                                  (cc^<.>e₂cc - 2 * (cc^<.>e₂c₂) + c₂e₂c₂ - 1)+           , 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+                  = return $+                 Shade' (c₀.+~^cc')+                        (HerMetric (Just ee) ^+^ projector (ee $ c₂^*η))+           +           | otherwise          = empty+   where σe = e₁^+^e₂+         quadraticEqnSol a b c+             | a /= 0 && disc > 0  = [ (σ * sqrt disc - b) / (2*a)+                                     | σ <- [-1, 1] ]+             | otherwise           = [0]+          where disc = b^2 - 4*a*c+         middle (_:x:y:_) = [x,y]+         middle l = l+  refineShade' (Shade' _ (HerMetric Nothing)) s₂ = pure s₂+  refineShade' s₁ (Shade' _ (HerMetric Nothing)) = pure s₁+  -- ⟨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).@@ -728,31 +905,31 @@ 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+propagateDEqnSolution_loc :: ∀ x y . (WithField ℝ Manifold x, Refinable y)+           => DifferentialEqn x y -> ((x, Shade' y), NonEmpty (Needle x, Shade' y))+                   -> NonEmpty (Shade' y)+propagateDEqnSolution_loc f ((x, shy@(Shade' y _)), neighbours) = ycs  where jShade@(Shade' j₀ jExpa) = f shxy        [shxy] = pointsCovers [ (xs, ys')-                             | (xs, Shade' ys yse) <- (x,shy):NE.toList neighbours+                             | (xs, Shade' ys yse)+                                 <- (x,shy):(first (x.+~^)<$>NE.toList neighbours)                              , δy <- eigenCoSpan' yse                              , ys' <- [ys.+~^δy, ys.-~^δy] ]-       [Shade' _ expax] = pointsCover's $ x : (fst<$>NE.toList neighbours)+       [Shade' _ expax] = pointsCover's $ x : ((x.+~^).fst<$>NE.toList neighbours)        marginδs :: NonEmpty (Needle x, (Needle y, Metric y))        marginδs = [ (δxm, (δym, expany))-                  | (xn, Shade' yn expany) <- neighbours-                  , let (Option (Just δxm)) = xn.-~.x-                        (Option (Just δym)) = yn.-~.y+                  | (δxm, Shade' yn expany) <- neighbours+                  , let (Option (Just δym)) = yn.-~.y                   ]        back2Centre :: (Needle x, (Needle y, Metric y)) -> Shade' y        back2Centre (δx, (δym, expany))             = convolveShade'                 (Shade' y expany)-                (Shade' δyb $ applyLinMapMetric jExpa δx')+                (Shade' δyb $ applyLinMapMetric jExpa (δx'^/(δx'<.>^δx)))         where δyb = δym ^-^ (j₀ $ δx)               δx' = toDualWith expax δx-       yc :: Option (Shade' y)-       yc = intersectShade's $ back2Centre <$> marginδs+       ycs :: NonEmpty (Shade' y)+       ycs = back2Centre <$> marginδs        xSpan = eigenCoSpan' expax  @@ -763,12 +940,6 @@     => 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-  pure = OuterJust . pure-  OuterJust fs <*> OuterJust xs = OuterJust $ fs <*> xs-  _ <*> _ = OuterNothing- traverseTwigsWithEnvirons :: ∀ x f .             (WithField ℝ Manifold x, Hask.Applicative f)     => ( ((Int, ShadeTree x), [(Int, ShadeTree x)]) -> f (ShadeTree x))@@ -1444,6 +1615,11 @@  stripShadedUntopological :: x`Shaded`y -> ShadeTree x stripShadedUntopological = unsafeFmapTree (fmap _topological) id shadeWithoutAnything++fmapShaded :: (y -> υ) -> (x`Shaded`y) -> (x`Shaded`υ)+fmapShaded f = unsafeFmapTree (fmap $ \(WithAny y x) -> WithAny (f y) x)+                              id+                              (\(Shade yx shx) -> Shade (fmap f yx) shx)  -- | This is to 'ShadeTree' as 'Data.Map.Map' is to 'Data.Set.Set'. type x`Shaded`y = ShadeTree (x`WithAny`y)

@@ -50,9 +50,12 @@         , D¹(..), D²(..)         , ℝay         , CD¹(..), Cℝay(..)-        -- * Cut-planes+        -- * Affine subspaces+        -- ** Lines+        , Line(..), lineAsPlaneIntersection+        -- ** Hyperplanes         , Cutplane(..)-        , fathomCutDistance, sideOfCut+        , fathomCutDistance, sideOfCut, cutPosBetween         -- * Linear mappings         , Linear, LocalLinear, denseLinear    ) where@@ -63,9 +66,7 @@ import Data.MemoTrie (HasTrie(..)) import Data.Basis import Data.Fixed-import Data.Void import Data.Tagged-import Data.Monoid import Data.Semigroup import qualified Numeric.LinearAlgebra.HMatrix as HMat import qualified Data.Vector.Generic as Arr@@ -221,6 +222,8 @@   +data Line x = Line { lineHandle :: x+                   , lineDirection :: Stiefel1 (Needle' x) }   @@ -258,4 +261,16 @@  where fathom v = (cn <.>^ v) / scaleDist        scaleDist = metric' met cn           ++cutPosBetween :: WithField ℝ Manifold x => Cutplane x -> (x,x) -> Option D¹+cutPosBetween (Cutplane h (Stiefel1 cn)) (x₀,x₁)+    | Option (Just [d₀,d₁]) <- map (cn<.>^) <$> sequenceA [x₀.-~.h, x₁.-~.h]+    , d₀*d₁ < 0+                  = pure . D¹ $ d₁ / (d₁ - d₀)+    | otherwise   = empty+++lineAsPlaneIntersection :: WithField ℝ Manifold x => Line x -> [Cutplane x]+lineAsPlaneIntersection (Line h dir)+      = [Cutplane h nrml | nrml <- orthogonalComplementSpan [dir]] 

@@ -44,7 +44,7 @@         -- * Projective spaces         , ℝP¹,  ℝP²(..)         -- * Intervals\/disks\/cones-        , D¹(..), D²(..)+        , D¹(..), fromIntv0to1, D²(..)         , ℝay         , CD¹(..), Cℝay(..)         -- * Tensor products@@ -59,7 +59,6 @@ import Data.VectorSpace import Data.AffineSpace import Data.Basis-import Data.Complex hiding (magnitude) import Data.Void import Data.Monoid @@ -67,12 +66,12 @@  import Control.Applicative (Const(..), Alternative(..)) +import Lens.Micro ((^.))+ import qualified Prelude  import Control.Category.Constrained.Prelude hiding ((^)) import Control.Arrow.Constrained-import Control.Monad.Constrained-import Data.Foldable.Constrained  import Data.Embedding @@ -155,6 +154,10 @@ --   the two points -1 and 1 of 'S⁰', i.e. this is simply a closed interval. newtype D¹ = D¹ { xParamD¹ :: Double -- ^ Range @[-1, 1]@.                 }+fromIntv0to1 :: ℝ -> D¹+fromIntv0to1 x | x<0        = D¹ (-1)+               | x>1        = D¹ 1+               | otherwise  = D¹ $ (x+1)/2  -- | The standard, closed unit disk. Homeomorphic to the cone over 'S¹', but not in the --   the obvious, &#x201c;flat&#x201d; way. (And not at all, despite@@ -286,8 +289,4 @@ (^) = (Prelude.^)  -infixl 8 ^.-{-# INLINE (^.) #-}-(^.) :: s -> (forall f . Prelude.Functor f => (a->f a) -> s->f s) -> a-o ^. g = getConst (g Const o) 

@@ -16,7 +16,7 @@ {-# LANGUAGE DeriveTraversable          #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE TypeFamilies               #-}-{-# LANGUAGE FunctionalDependencies     #-}+{-# LANGUAGE MultiParamTypeClasses      #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE GADTs                      #-} {-# LANGUAGE RankNTypes                 #-}@@ -25,17 +25,32 @@ {-# LANGUAGE UnicodeSyntax              #-} {-# LANGUAGE ConstraintKinds            #-} {-# LANGUAGE PatternGuards              #-}-{-# LANGUAGE PatternSynonyms            #-}-{-# LANGUAGE ViewPatterns               #-} {-# LANGUAGE LambdaCase                 #-} {-# LANGUAGE TypeOperators              #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE LiberalTypeSynonyms        #-}-{-# LANGUAGE RecordWildCards            #-}-{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE TemplateHaskell            #-}  -module Data.Manifold.Web where+module Data.Manifold.Web (+              -- * The web data type+              PointsWeb+              -- ** Construction+            , fromWebNodes, fromShadeTree_auto, fromShadeTree, fromShaded+              -- ** Lookup+            , nearestNeighbour, indexWeb, webEdges, toGraph+              -- ** Decomposition+            , sliceWeb_lin -- , sampleWebAlongLine_lin+              -- ** Local environments+            , localFocusWeb+              -- * Differential equations+              -- ** Fixed resolution+            , filterDEqnSolution_static, iterateFilterDEqn_static+              -- ** Automatic resolution+            , filterDEqnSolutions_adaptive, iterateFilterDEqn_adaptive+              -- * Misc+            , ConvexSet(..), ellipsoid+            ) where   import Data.List hiding (filter, all, elem, sum, foldr1)@@ -44,49 +59,35 @@ 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.List.FastNub (fastNubBy) import Data.Ord (comparing)+import Data.Semigroup 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.Function (on)+import Data.Fixed (mod') -import Data.SimplicialComplex import Data.Manifold.Types-import Data.Manifold.Types.Primitive ((^), empty)+import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine-import Data.Function.Differentiable-import Data.Function.Differentiable.Data import Data.Manifold.TreeCover+import Data.SetLike.Intersection+import Data.Manifold.Riemannian     -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 Control.Monad.Trans.List import qualified Data.Foldable       as Hask-import Data.Foldable (all, elem, toList, sum, foldr1)+import Data.Foldable (all, toList) import qualified Data.Traversable as Hask import Data.Traversable (forM)--import qualified Numeric.LinearAlgebra.HMatrix as HMat+import Data.Graph  import Control.Category.Constrained.Prelude hiding      ((^), all, elem, sum, forM, Foldable(..), foldr1, Traversable, traverse)@@ -95,20 +96,36 @@ import Data.Foldable.Constrained import Data.Traversable.Constrained (Traversable, traverse) +import Control.Comonad (Comonad(..))+import Lens.Micro ((&), (%~), (^.), (.~))+import Lens.Micro.TH+ import GHC.Generics (Generic)   type WebNodeId = Int-type NeighbourRefs = UArr.Vector WebNodeId +data Neighbourhood x = Neighbourhood {+     neighbours :: UArr.Vector WebNodeId+   , localScalarProduct :: Metric x+   }+  deriving (Generic)++instance (NFData x, NFData (HerMetric (Needle x))) => NFData (Neighbourhood x)++-- | A 'PointsWeb' is almost, but not quite a mesh. It is a stongly connected†+--   directed graph, backed by a tree for fast nearest-neighbour lookup of points.+-- +--   †In general, there can be disconnected components, but every connected+--   component is strongly connected. data PointsWeb :: * -> * -> * where    PointsWeb :: {        webNodeRsc :: ShadeTree x-     , webNodeAssocData :: Arr.Vector (y, NeighbourRefs)+     , webNodeAssocData :: Arr.Vector (y, Neighbourhood x)      } -> PointsWeb x y   deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable) -instance (NFData x, NFData (Needle' x), NFData y) => NFData (PointsWeb x y)+instance (NFData x, NFData (HerMetric (Needle x)), NFData (Needle' x), NFData y) => NFData (PointsWeb x y)  instance Foldable (PointsWeb x) (->) (->) where   ffoldl = uncurry . Hask.foldl' . curry@@ -121,10 +138,18 @@   +type MetricChoice x = Shade x -> Metric x++ fromWebNodes :: ∀ x y . WithField ℝ Manifold x-                    => (Shade x->Metric x) -> [(x,y)] -> PointsWeb x y+                    => (MetricChoice x) -> [(x,y)] -> PointsWeb x y fromWebNodes mf = fromShaded mf . fromLeafPoints . map (uncurry WithAny . swap) +fromTopWebNodes :: ∀ x y . WithField ℝ Manifold x+                    => (MetricChoice x) -> [((x,[Needle x]),y)] -> PointsWeb x y+fromTopWebNodes mf = fromTopShaded mf . fromLeafPoints+                   . map (uncurry WithAny . swap . regroup')+ fromShadeTree_auto :: ∀ x . WithField ℝ Manifold x => ShadeTree x -> PointsWeb x () fromShadeTree_auto = fromShaded (recipMetric . _shadeExpanse) . constShaded () @@ -133,33 +158,46 @@ fromShadeTree mf = fromShaded mf . constShaded ()  fromShaded :: ∀ x y . WithField ℝ Manifold x-     => (Shade x -> Metric x) -- ^ Local scalar-product generator. You can always+     => (MetricChoice 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 +fromShaded metricf = fromTopShaded metricf . fmapShaded ([],)++fromTopShaded :: ∀ x y . WithField ℝ Manifold x+     => (MetricChoice x)+     -> (x`Shaded`([Needle x], y))  -- ^ Source tree, with a priori topology information+                                    --   (needles pointing to already-known neighbour candidates)+     -> PointsWeb x y+fromTopShaded 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+       locMesh :: ( (Int, ShadeTree (x`WithAny`([Needle x], y)))+                  , [(Int, ShadeTree (x`WithAny`([Needle x], y)))])+                   -> Arr.Vector (y, Neighbourhood x)+       locMesh ((i₀, locT), neighRegions) = Arr.map findNeighbours $ Arr.fromList locLeaves+        where locLeaves :: [ (Int, x`WithAny`([Needle x], y)) ]+              locLeaves = map (first (+i₀)) . zip [0..] $ onlyLeaves locT+              vicinityLeaves :: [(Int, x)]               vicinityLeaves = Hask.foldMap-                                (\(i₀n, ngbR) -> Arr.map (first (+i₀n))-                                               . Arr.indexed-                                               . Arr.fromList+                                (\(i₀n, ngbR) -> map ((+i₀n) *** _topological)+                                               . zip [0..]                                                $ 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) ->+              findNeighbours :: (Int, x`WithAny`([Needle x], y)) -> (y, Neighbourhood x)+              findNeighbours (i, WithAny (vns,y) x)+                         = (y, Neighbourhood+                                 (UArr.fromList $ fst<$>execState seek mempty)+                                 locRieM )+               where seek :: State [(Int, (Needle x, Needle' x))] ()+                     seek = do+                        Hask.forM_ ( fastNubBy (comparing fst)+                                      $ map (second _topological) locLeaves+                                           ++ vicinityLeaves ++ aprioriNgbs )+                                  $ \(iNgb, xNgb) ->                            when (iNgb/=i) `id`do                               let (Option (Just v)) = xNgb.-~.x                               oldNgbs <- get@@ -171,6 +209,12 @@                                          | neighbour@(_,(nv,_))<-oldNgbs                                          , visibleOverlap w nv                                          ]+                     aprioriNgbs :: [(Int, x)]+                     aprioriNgbs = catMaybes+                                    [ getOption $ (second $ const xN) <$>+                                          positionIndex (pure locRieM) shd' xN+                                    | v <- vns+                                    , let xN = x.+~^v :: x ]                              visibleOverlap :: Needle' x -> Needle x -> Bool               visibleOverlap w v = o < 1@@ -188,32 +232,429 @@   , Right (_,x) <- indexShadeTree rsc i  = pure (x, fst (assocD Arr.! i))   | otherwise                            = empty +unsafeIndexWebData :: PointsWeb x y -> WebNodeId -> y+unsafeIndexWebData (PointsWeb _ asd) i = fst (asd Arr.! i)+ 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))+       allEdges = Hask.foldMap (\(i,(_, Neighbourhood 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)])+data InterpolationIv y = InterpolationIv {+          _interpolationSegRange :: (ℝ,ℝ)+        , _interpolationFunction :: ℝ -> y+        }++type InterpolationSeq y = [InterpolationIv y]++mkInterpolationSeq_lin :: (x~ℝ, Geodesic y)+           => [(x,y)] -> InterpolationSeq y+mkInterpolationSeq_lin [(xψ,yψ), (xω,yω)]+       = return $ InterpolationIv+           (xψ,xω)+           (\x -> let drel = fromIntv0to1 $ (x-xψ)/(xω-xψ)+                  in yio drel )+ where Option (Just yio) = geodesicBetween yψ yω+mkInterpolationSeq_lin (p₀:p₁:ps)+    = mkInterpolationSeq_lin [p₀,p₁] <> mkInterpolationSeq_lin (p₁:ps)+mkInterpolationSeq_lin _ = []+++-- | Fetch a point between any two neighbouring web nodes on opposite+--   sides of the plane, and linearly interpolate the values onto the+--   cut plane.+sliceWeb_lin :: ∀ x y . (WithField ℝ Manifold x, Geodesic x, Geodesic y)+               => PointsWeb x y -> Cutplane x -> [(x,y)]+sliceWeb_lin web = sliceEdgs+ where edgs = webEdges web+       sliceEdgs cp = [ (xi d, yi d)  -- Brute-force search through all edges+                      | ((x₀,y₀), (x₁,y₁)) <- edgs+                      , Option (Just d) <- [cutPosBetween cp (x₀,x₁)]+                      , Option (Just xi) <- [geodesicBetween x₀ x₁]+                      , Option (Just yi) <- [geodesicBetween y₀ y₁]+                      ]++-- sampleWebAlongLine_lin :: ∀ x y . (WithField ℝ Manifold x, Geodesic x, Geodesic y)+--                => PointsWeb x y -> x -> Needle x -> [(x,y)]+-- sampleWebAlongLine_lin web x₀ dir = sampleWebAlongLines_lin web x₀ [(dir, maxBound)]+++data GridPlanes x = GridPlanes {+        _gridPlaneNormal :: Needle' x+      , _gridPlaneSpacing :: Needle x+      , _gridPlanesCount :: Int+      }+data GridSetup x = GridSetup {+        _gridStartCorner :: x+      , _gridSplitDirs :: [GridPlanes x]+      }++cartesianGrid2D :: (x~ℝ, y~ℝ) => ((x,x), Int) -> ((y,y), Int) -> GridSetup (x,y)+cartesianGrid2D ((x₀,x₁), nx) ((y₀,y₁), ny)+    = GridSetup (x₀,y₀) [ GridPlanes (0,1) (0, (y₁-y₀)/fromIntegral ny) ny+                        , GridPlanes (1,0) ((x₁-x₀)/fromIntegral nx, 0) ny ]++splitToGridLines :: (WithField ℝ Manifold x, Geodesic x, Geodesic y)+          => PointsWeb x y -> GridSetup x -> [((x, GridPlanes x), [(x,y)])]+splitToGridLines web (GridSetup x₀ [GridPlanes dirΩ spcΩ nΩ, linePln])+    = [ ((x₀', linePln), sliceWeb_lin web $ Cutplane x₀' (Stiefel1 dirΩ))+      | k <- [0 .. nΩ-1]+      , let x₀' = x₀.+~^(fromIntegral k *^ spcΩ) ]++sampleWebAlongGrid_lin :: ∀ x y . (WithField ℝ Manifold x, Geodesic x, Geodesic y)+               => PointsWeb x y -> GridSetup x -> [(x,Option y)]+sampleWebAlongGrid_lin web grid = finalLine =<< splitToGridLines web grid+ where finalLine :: ((x, GridPlanes x), [(x,y)]) -> [(x,Option y)]+       finalLine ((x₀, GridPlanes _ dir nSpl), verts)+          | length verts < 2  = take nSpl $ (,empty)<$>iterate (.+~^dir) x₀+       finalLine ((x₀, GridPlanes _ dir nSpl), verts)  = take nSpl $ go (x₀,0) intpseq +        where intpseq = mkInterpolationSeq_lin+                         [ (metric metr $ x.-~!x₀, y) | (x,y) <- verts ]+              go (x,_) [] = (,empty)<$>iterate (.+~^dir) x+              go xt (InterpolationIv (_,te) f:fs)+                        = case break ((<te) . snd) $ iterate ((.+~^dir)***(+1)) xt of+                             (thisRange, xtn:_)+                                 -> ((id***pure.f)<$>thisRange) ++ go xtn fs+       Option (Just metr) = inferMetric $ webNodeRsc web+       +sampleWeb_2Dcartesian_lin :: (x~ℝ, y~ℝ, Geodesic z)+             => PointsWeb (x,y) z -> ((x,x),Int) -> ((y,y),Int) -> [(y,[(x,Option z)])]+sampleWeb_2Dcartesian_lin web (xspec@(_,nx)) yspec+       = go . sampleWebAlongGrid_lin web $ cartesianGrid2D xspec yspec+ where go [] = []+       go l@(((_,y),_):_) = let (ln,l') = splitAt nx l+                             in (y, map (\((x,_),z) -> (x,z)) ln) : go l'+       +sampleEntireWeb_2Dcartesian_lin :: (x~ℝ, y~ℝ, Geodesic z)+             => PointsWeb (x,y) z -> Int -> Int -> [(y,[(x,Option z)])]+sampleEntireWeb_2Dcartesian_lin web nx ny+       = sampleWeb_2Dcartesian_lin web ((x₀,x₁),nx) ((y₀,y₁),ny)+ where x₀ = minimum (fst<$>pts)+       x₁ = maximum (fst<$>pts)+       y₀ = minimum (snd<$>pts)+       y₁ = maximum (snd<$>pts)+       pts = fst . fst <$> toList (localFocusWeb web)++webLocalInfo :: ∀ x y . WithField ℝ Manifold x+            => PointsWeb x y -> PointsWeb x (WebLocally x y)+webLocalInfo origWeb = result+ where result = wli $ localFocusWeb origWeb+       wli (PointsWeb rsc asd) = PointsWeb rsc asd'+        where asd' = Arr.imap localInfo asd+       localInfo i (((x,y), ngbCo), ngbH)+            = ( LocalWebInfo {+                  _thisNodeCoord = x+                , _thisNodeData = y+                , _containingWeb = result+                , _thisNodeId = i+                , _nodeNeighbours = zip (UArr.toList $ neighbours ngbH) ngbCo+                , _nodeLocalScalarProduct = localScalarProduct ngbH+                , _nodeIsOnBoundary = anyUnopposed (localScalarProduct ngbH) ngbCo+                }, ngbH )+       anyUnopposed rieM ngbCo = (`any`ngbCo) $ \(v,_)+                         -> not $ (`any`ngbCo) $ \(v',_)+                              -> toDualWith rieM v <.>^ v' < 0++localFocusWeb :: WithField ℝ Manifold x+                   => PointsWeb x y -> PointsWeb x ((x,y), [(Needle x, y)]) localFocusWeb (PointsWeb rsc asd) = PointsWeb rsc asd''  where asd' = Arr.imap (\i (y,n) -> case indexShadeTree rsc i of                                          Right (_,x) -> ((x,y),n) ) asd-       asd''= Arr.map (\(xy,n) ->-                       ((xy, [fst (asd' Arr.! j) | j<-UArr.toList n]), n)+       asd''= Arr.map (\((x,y),n) ->+                       (((x,y), [ ( case x'.-~.x of+                                     Option (Just v) -> v+                                  , y')+                                | j<-UArr.toList (neighbours n)+                                , let ((x',y'),_) = asd' Arr.! j+                                ]), n)                  ) asd'  +nearestNeighbour :: WithField ℝ Manifold x+                      => PointsWeb x y -> x -> Option (x,y)+nearestNeighbour (PointsWeb rsc asd) x = fmap lkBest $ positionIndex empty rsc x+ where lkBest (iEst, (_, xEst)) = (xProx, yProx)+        where (iProx, (xProx, _)) = minimumBy (comparing $ snd . snd)+                                     $ (iEst, (xEst, metricSq locMetr vEst))+                                         : neighbours+              (yProx, _) = asd Arr.! iProx+              (_, Neighbourhood neighbourIds locMetr) = asd Arr.! iEst+              neighbours = [ (i, (xNgb, metricSq locMetr v))+                           | i <- UArr.toList neighbourIds+                           , let Right (_, xNgb) = indexShadeTree rsc i+                                 Option (Just v) = xNgb.-~.x+                           ]+              Option (Just vEst) = xEst.-~.x++++data WebLocally x y = LocalWebInfo {+      _thisNodeCoord :: x+    , _thisNodeData :: y+    , _containingWeb :: PointsWeb x (WebLocally x y)+    , _thisNodeId :: WebNodeId+    , _nodeNeighbours :: [(WebNodeId, (Needle x, y))]+    , _nodeLocalScalarProduct :: Metric x+    , _nodeIsOnBoundary :: Bool+    } deriving (Generic)+makeLenses ''WebLocally++instance Hask.Functor (WebLocally x) where+  fmap f (LocalWebInfo co dt wb id ng sp bn)+       = LocalWebInfo co (f dt) (fmap (fmap f) wb) id (map (second $ second f) ng) sp bn+instance WithField ℝ Manifold x => Comonad (WebLocally x) where+  extract = _thisNodeData+  duplicate lweb = unsafeIndexWebData deepened $ _thisNodeId lweb+   where deepened = webLocalInfo $ _containingWeb lweb++++++toGraph :: WithField ℝ Manifold x => PointsWeb x y -> (Graph, Vertex -> (x, y))+toGraph wb = second (>>> \(i,_,_) -> case indexWeb wb i of {Option (Just xy) -> xy})+                (graphFromEdges' edgs)+ where edgs :: [(Int, Int, [Int])]+       edgs = Arr.toList+            . Arr.imap (\i (_, Neighbourhood ngbs _) -> (i, i, UArr.toList ngbs))+                    $ webNodeAssocData wb+++++data ConvexSet x+    = EmptyConvex+    | ConvexSet {+      convexSetHull :: Shade' x+      -- ^ If @p@ is in all intersectors, it must also be in the hull.+    , convexSetIntersectors :: [Shade' x]+    }++ellipsoid :: Shade' x -> ConvexSet x+ellipsoid s = ConvexSet s [s]++intersectors :: ConvexSet x -> Option (NonEmpty (Shade' x))+intersectors (ConvexSet h []) = pure (h:|[])+intersectors (ConvexSet _ (i:sts)) = pure (i:|sts)+intersectors _ = empty++-- | Under intersection.+instance Refinable x => Semigroup (ConvexSet x) where+  a<>b = sconcat (a:|[b])+  sconcat csets+    | Option (Just allIntersectors) <- sconcat <$> Hask.traverse intersectors csets+    , IntersectT ists <- rmTautologyIntersect perfectRefine $ IntersectT allIntersectors+    , Option (Just hull') <- intersectShade's ists+                 = ConvexSet hull' (NE.toList ists)+    | otherwise  = EmptyConvex+   where perfectRefine sh₁ sh₂+           | sh₁`subShade'`sh₂   = pure sh₁+           | sh₂`subShade'`sh₁   = pure sh₂+           | otherwise           = empty++++itWhileJust :: (a -> Option a) -> a -> [a]+itWhileJust f x | Option (Just y) <- f x  = x : itWhileJust f y+itWhileJust _ x = [x]++dupHead :: NonEmpty a -> NonEmpty a+dupHead (x:|xs) = x:|x:xs+++iterateFilterDEqn_static :: (WithField ℝ Manifold x, Refinable y)+       => DifferentialEqn x y -> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)]+iterateFilterDEqn_static f = map (fmap convexSetHull)+                           . itWhileJust (filterDEqnSolutions_static f)+                           . fmap (`ConvexSet`[])+ 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)+                            =<< intersectShade's+                                  ( propagateDEqnSolution_loc f ((x,shy), NE.fromList ngbs) )++filterDEqnSolutions_static :: (WithField ℝ Manifold x, Refinable y)+       => DifferentialEqn x y -> PointsWeb x (ConvexSet y) -> Option (PointsWeb x (ConvexSet y))+filterDEqnSolutions_static f = localFocusWeb >>> Hask.traverse `id`+            \((x, shy@(ConvexSet hull _)), ngbs) -> if null ngbs+              then pure shy+              else ((shy<>) . ellipsoid)+                      <$> intersectShade's +                            ( propagateDEqnSolution_loc f+                               ((x,hull), second convexSetHull<$>NE.fromList ngbs) )+                     >>= \case EmptyConvex -> empty+                               c           -> pure c+++data SolverNodeState y = SolverNodeInfo {+      _solverNodeStatus :: ConvexSet y+    , _solverNodeBadness :: ℝ+    , _solverNodeAge :: Int+    }+makeLenses ''SolverNodeState+++type OldAndNew d = (Option d, [d])++oldAndNew :: OldAndNew d -> [d]+oldAndNew (Option (Just x), l) = x : l+oldAndNew (_, l) = l++oldAndNew' :: OldAndNew d -> [(Bool, d)]+oldAndNew' (Option (Just x), l) = (True, x) : fmap (False,) l+oldAndNew' (_, l) = (False,) <$> l+++filterDEqnSolutions_adaptive :: ∀ x y badness+        . (WithField ℝ Manifold x, Refinable y, badness ~ ℝ)+       => MetricChoice x      -- ^ Scalar product on the domain, for regularising the web.+       -> DifferentialEqn x y +       -> (x -> Shade' y -> badness)+             -> PointsWeb x (SolverNodeState y)+                        -> Option (PointsWeb x (SolverNodeState y))+filterDEqnSolutions_adaptive mf f badness' oldState+         = fmap (fromTopWebNodes mf . concat . fmap retraceBonds+                                        . Hask.toList . webLocalInfo . webLocalInfo)+             $ Hask.traverse (uncurry localChange) preproc'd+ where preproc'd :: PointsWeb x ((WebLocally x (SolverNodeState y), [(Shade' y, badness)]))+       preproc'd = fmap addPropagation $ webLocalInfo oldState+        where addPropagation wl+                 | null neighbourHulls = (wl, [])+                 | otherwise           = (wl, map (id&&&badness undefined) propFromNgbs)+               where propFromNgbs = NE.toList $ propagateDEqnSolution_loc f+                                     ( (thisPos, thisShy), NE.fromList neighbourHulls )+                     thisPos = _thisNodeCoord wl :: x+                     thisShy = convexSetHull . _solverNodeStatus $ _thisNodeData wl+                     neighbourHulls = second (convexSetHull . _solverNodeStatus) . snd+                                        <$> _nodeNeighbours wl+       smallBadnessGradient, largeBadnessGradient :: ℝ+       (smallBadnessGradient, largeBadnessGradient)+           = ( badnessGradRated!!(n`div`4), badnessGradRated!!(n*3`div`4) )+        where n = length badnessGradRated+              badnessGradRated = sort [ ngBad / bad+                                      | ( LocalWebInfo {+                                            _thisNodeData+                                              = SolverNodeInfo _ bad _+                                          , _nodeNeighbours=ngbs        }+                                        , ngbProps) <- Hask.toList preproc'd+                                      , (_, ngBad) <- ngbProps+                                      , ngBad>bad ]+       localChange :: WebLocally x (SolverNodeState y) -> [(Shade' y, badness)]+                             -> Option (OldAndNew (x, SolverNodeState y))+       localChange localInfo@LocalWebInfo{+                         _thisNodeCoord = x+                       , _thisNodeData = SolverNodeInfo+                                            shy@(ConvexSet hull _) prevBadness age+                       , _nodeNeighbours = ngbs+                       }+                   ngbProps+        | null ngbs  = return (pure (x, SolverNodeInfo shy prevBadness (age+1)), [])+        | otherwise  = do+               let neighbourHulls = second (convexSetHull . _solverNodeStatus) . snd+                                       <$> NE.fromList ngbs+                   (environAge, unfreshness)+                      = maximum&&&minimum $ age : (_solverNodeAge . snd . snd <$> ngbs)+               case find (\(_, badnessN)+                               -> badnessN / prevBadness > smallBadnessGradient)+                              $ ngbProps of+                 Nothing | age < environAge   -- point is an obsolete step-stone;+                   -> return (empty,empty)    -- do not further use it.+                 _otherwise -> do+                   shy' <- ((shy<>) . ellipsoid)+                            <$> intersectShade's (fst <$> NE.fromList ngbProps)+                   newBadness <- case shy' of+                      EmptyConvex        -> empty+                      ConvexSet hull' _  -> return $ badness x hull'+                   let updatedNode = SolverNodeInfo shy' newBadness (age+1)+                   stepStones <-+                     if unfreshness < 3+                      then return []+                      else fmap concat . forM (zip (snd<$>ngbs) ngbProps)+                                   $ \( (vN, SolverNodeInfo (ConvexSet hullN _)+                                                          _ ageN)+                                        , (_, nBadnessProp'd) ) -> do+                       case ageN of+                        _  | ageN > 0+                           , badnessGrad <- nBadnessProp'd / prevBadness+                           , badnessGrad > largeBadnessGradient -> do+                                 let stepV = vN^/2+                                     xStep = x .+~^ stepV+                                 shyStep <- intersectShade's $+                                            propagateDEqnSolution_loc f+                                            ( (xStep, hull)+                                            , NE.cons (negateV stepV, hull)+                                                $ fmap (\(vN',hullN')+                                                         -> (vN'^-^stepV, hullN') )+                                                    neighbourHulls )+                                 return [( xStep+                                         , SolverNodeInfo (ellipsoid shyStep)+                                                 (badness xStep shyStep) 1+                                         )]+                        _otherwise -> return []+                   let updated = (x, updatedNode)+                   return $ (pure updated, stepStones)+       +       totalAge = maximum $ _solverNodeAge . _thisNodeData . fst <$> preproc'd+       errTgtModulation = (1-) . (`mod'`1) . negate . sqrt $ fromIntegral totalAge+       badness x = badness' x . (shadeNarrowness %~ (^* errTgtModulation))+       +       retraceBonds :: WebLocally x (WebLocally x (OldAndNew (x, SolverNodeState y)))+                       -> [((x, [Needle x]), SolverNodeState y)]+       retraceBonds locWeb@LocalWebInfo{ _thisNodeId = myId+                                       , _thisNodeCoord = xOld+                                       , _nodeLocalScalarProduct = locMetr }+            = [ ( (x, fst<$>neighbourCandidates), snsy)+              | (isOld, (x, snsy)) <- focused+              , let neighbourCandidates+                     = [ (v,nnId)+                       | (_,ngb) <- knownNgbs+                       , (Option (Just v), nnId)+                          <- case oldAndNew $ ngb^.thisNodeData of+                                   [] -> [ (xN.-~.x, nnId)+                                         | (nnId, (_,nnWeb)) <- ngb^.nodeNeighbours+                                         , nnId /= myId+                                         , (xN,_) <- oldAndNew nnWeb ]+                                   l -> [(xN.-~.x, ngb^.thisNodeId) | (xN,_) <- l]+                       ]+                    possibleConflicts = [ metricSq locMetr v+                                        | (v,nnId)<-neighbourCandidates+                                        , nnId > myId ]+              , isOld || null possibleConflicts+                  || minimum possibleConflicts > oldMinDistSq / 4+              ]+        where focused = oldAndNew' $ locWeb^.thisNodeData^.thisNodeData+              knownNgbs = snd <$> locWeb^.nodeNeighbours+              oldMinDistSq = minimum [ metricSq locMetr vOld+                                     | (_,ngb) <- knownNgbs+                                     , let Option (Just vOld) = ngb^.thisNodeCoord .-~. xOld+                                     ]+                              +++iterateFilterDEqn_adaptive :: (WithField ℝ Manifold x, Refinable y)+       => MetricChoice x      -- ^ Scalar product on the domain, for regularising the web.+       -> DifferentialEqn x y+       -> (x -> Shade' y -> ℝ) -- ^ Badness function for local results.+             -> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)]+iterateFilterDEqn_adaptive mf f badness+    = map (fmap (convexSetHull . _solverNodeStatus))+    . itWhileJust (filterDEqnSolutions_adaptive mf f badness)+    . fmap (\((x,shy),_) -> SolverNodeInfo (ellipsoid shy)+                                           (badness x shy)+                                           1+           )+    . localFocusWeb++  

@@ -64,22 +64,15 @@  import Data.List hiding (filter, all, elem) import Data.Maybe-import qualified Data.Map as Map import qualified Data.Vector as Arr-import Data.List.NonEmpty (NonEmpty(..)) import Data.List.FastNub import qualified Data.List.NonEmpty as NE import Data.Semigroup import Data.Ord (comparing) -import Data.VectorSpace-import Data.LinearMap import Data.LinearMap.Category-import Data.Void import Data.Tagged-import Data.Proxy -import Data.Manifold.Types import Data.Manifold.Types.Primitive ((^), empty) import Data.Manifold.PseudoAffine     @@ -101,7 +94,6 @@ import Control.Monad.Constrained import Data.Foldable.Constrained -import GHC.Generics (Generic)  infixr 5 :<|, .<. 

@@ -6,7 +6,6 @@ module Util.LtdShow (LtdShow(..)) where  import qualified Data.Vector as V-import Data.Vector(fromList, toList, (!), singleton)  type Array = V.Vector 

@@ -1,5 +1,5 @@ Name:                manifolds-Version:             0.2.2.0+Version:             0.2.3.0 Category:            Math Synopsis:            Coordinate-free hypersurfaces Description:         Manifolds, a generalisation of the notion of &#x201c;smooth curves&#x201d; or surfaces,@@ -49,6 +49,7 @@                      , void                      , tagged                      , deepseq+                     , microlens >= 0.4 && <= 0.5, microlens-th                      , trivial-constraint >= 0.4                      , constrained-categories >= 0.2.3 && < 0.3   other-extensions:  FlexibleInstances@@ -66,6 +67,7 @@                      Data.Manifold.PseudoAffine                      Data.Manifold.TreeCover                      Data.Manifold.Web+                     Data.Manifold.DifferentialEquation                      Data.SimplicialComplex                      Data.LinearMap.HerMetric                      Data.Function.Differentiable@@ -82,6 +84,7 @@                    Data.Function.Differentiable.Data                    Data.Function.Affine                    Data.VectorSpace.FiniteDimensional+                   Control.Monad.Trans.OuterMaybe                    Util.Associate                    Util.LtdShow   default-language: Haskell2010