diff --git a/Data/Function/Affine.hs b/Data/Function/Affine.hs --- a/Data/Function/Affine.hs +++ b/Data/Function/Affine.hs @@ -33,6 +33,8 @@ Affine(..) , evalAffine , fromOffsetSlope + -- * Misc + , lensEmbedding, correspondingDirections ) where @@ -46,9 +48,11 @@ import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine import Data.Manifold.Atlas +import Data.Embedding import qualified Prelude import qualified Control.Applicative as Hask +import Data.Functor (($>)) import Control.Category.Constrained.Prelude hiding ((^)) import Control.Category.Constrained.Reified @@ -58,8 +62,10 @@ import Math.LinearMap.Category +import Control.Lens + data Affine s d c where Affine :: (ChartIndex d :->: (c, LinearMap s (Needle d) (Needle c))) -> Affine s d c @@ -198,3 +204,50 @@ -> \y0 ðx'y -> Affine . trie $ chartReferencePoint >>> \x₀ -> let δy = ðx'y $ x₀ in (y0.+~^δy, ðx'y) + + +instance EnhancedCat (Embedding (Affine s)) (Embedding (LinearMap s)) where + arr (Embedding e p) = Embedding (arr e) (arr p) + + +lensEmbedding :: ∀ k s x c . + ( Num' s + , LinearSpace x, LinearSpace c, Object k x, Object k c + , Scalar x ~ s, Scalar c ~ s + , EnhancedCat k (LinearMap s) ) + => Lens' x c -> Embedding k c x +lensEmbedding l = Embedding (arr $ (arr $ LinearFunction (\c -> zeroV & l .~ c) + :: LinearMap s c x) ) + (arr $ (arr $ LinearFunction (^.l) + :: LinearMap s x c) ) + + +correspondingDirections :: ∀ s x c t + . ( WithField s AffineManifold c + , WithField s AffineManifold x + , SemiInner (Needle c), SemiInner (Needle x) + , RealFrac' s + , Traversable t ) + => (Interior c, Interior x) + -> t (Needle c, Needle x) -> Maybe (Embedding (Affine s) c x) +correspondingDirections (c₀, x₀) dirMap + = freeEmbeddings $> Embedding (Affine . trie $ c2x boundarylessWitness) + (Affine . trie $ x2c boundarylessWitness) + where freeEmbeddings = fzip ( embedFreeSubspace $ fst<$>dirMap + , embedFreeSubspace $ snd<$>dirMap ) + c2t :: Lens' (Needle c) (t s) + c2t = case freeEmbeddings of Just (Lens ct, _) -> ct + x2t :: Lens' (Needle x) (t s) + x2t = case freeEmbeddings of Just (_, Lens xt) -> xt + c2x :: BoundarylessWitness c -> ChartIndex c + -> (x, LinearMap s (Needle c) (Needle x)) + c2x BoundarylessWitness ιc + = ( x₀ .+~^ (zeroV & x2t .~ δc^.c2t) + , arr . LinearFunction $ \dc -> zeroV & x2t .~ dc^.c2t ) + where Just δc = chartReferencePoint ιc .-~. c₀ + x2c :: BoundarylessWitness x -> ChartIndex x + -> (c, LinearMap s (Needle x) (Needle c)) + x2c BoundarylessWitness ιx + = ( c₀ .+~^ (zeroV & c2t .~ δx^.x2t) + , arr . LinearFunction $ \dx -> zeroV & c2t .~ dx^.x2t ) + where Just δx = chartReferencePoint ιx .-~. x₀ diff --git a/Data/Function/Differentiable.hs b/Data/Function/Differentiable.hs --- a/Data/Function/Differentiable.hs +++ b/Data/Function/Differentiable.hs @@ -77,6 +77,9 @@ +type RealDimension s + = ( RealFloat' s, SimpleSpace s, Show s, Atlas s, HasTrie (ChartIndex s) + , s ~ Needle s, s ~ Interior s, s ~ Scalar s, s ~ DualVector s ) discretisePathIn :: (WithField ℝ Manifold y, SimpleSpace (Needle y)) @@ -265,10 +268,12 @@ -unsafe_dev_ε_δ :: RealDimension a +unsafe_dev_ε_δ :: ∀ a . RealDimension a => String -> (a -> a) -> LinDevPropag a a -unsafe_dev_ε_δ errHint f d - = let ε'² = normSq d 1 +unsafe_dev_ε_δ = case ( linearManifoldWitness :: LinearManifoldWitness a + , closedScalarWitness :: ClosedScalarWitness a ) of + (LinearManifoldWitness _, ClosedScalarWitness) -> \errHint f d + -> let ε'² = normSq d 1 in if ε'²>0 then let δ = f . sqrt $ recip ε'² in if δ > 0 @@ -278,9 +283,12 @@ ++show(sqrt $ recip ε'²) ++ " gives non-positive δ="++show δ++"." else mempty -dev_ε_δ :: RealDimension a +dev_ε_δ :: ∀ a . RealDimension a => (a -> a) -> Metric a -> Maybe (Metric a) -dev_ε_δ f d = let ε'² = normSq d 1 +dev_ε_δ = case ( linearManifoldWitness :: LinearManifoldWitness a + , closedScalarWitness :: ClosedScalarWitness a ) of + (LinearManifoldWitness _, ClosedScalarWitness) -> \f d + -> let ε'² = normSq d 1 in if ε'²>0 then let δ = f . sqrt $ recip ε'² in if δ > 0 @@ -288,8 +296,11 @@ else empty else pure mempty -as_devεδ :: RealDimension a => LinDevPropag a a -> a -> a -as_devεδ ldp ε | ε>0 +as_devεδ :: ∀ a . RealDimension a => LinDevPropag a a -> a -> a +as_devεδ = asdevεδ linearManifoldWitness closedScalarWitness where + asdevεδ :: LinearManifoldWitness a -> ClosedScalarWitness a -> LinDevPropag a a -> a -> a + asdevεδ (LinearManifoldWitness _) ClosedScalarWitness + ldp ε | ε>0 , δ'² <- normSq (ldp $ spanNorm [recip ε]) 1 , δ'² > 0 = sqrt $ recip δ'² @@ -466,40 +477,61 @@ -instance (LocallyScalable s v, LinearManifold v, LocallyScalable s a, RealFloat' s) +instance ∀ v s a . (LinearSpace v, Scalar v ~ s, LocallyScalable s a, RealFloat' s) => AdditiveGroup (DfblFuncValue s a v) where - zeroV = point zeroV - GenericAgent (AffinDiffable ef f) ^+^ GenericAgent (AffinDiffable eg g) - = GenericAgent $ AffinDiffable (ef<>eg) (f^+^g) - α^+^β = dfblFnValsCombine (\a b -> (a^+^b, arr addV, const mempty)) α β - negateV (GenericAgent (AffinDiffable ef f)) - = GenericAgent $ AffinDiffable ef (negateV f) - negateV α = dfblFnValsFunc (\a -> (negateV a, negateV id, const mempty)) α + zeroV = case ( linearManifoldWitness :: LinearManifoldWitness v + , dualSpaceWitness :: DualSpaceWitness v ) of + (LinearManifoldWitness _, DualSpaceWitness) -> point zeroV + (^+^) = case ( linearManifoldWitness :: LinearManifoldWitness v + , dualSpaceWitness :: DualSpaceWitness v ) of + (LinearManifoldWitness _, DualSpaceWitness) + -> curry $ \case + (GenericAgent (AffinDiffable ef f), GenericAgent (AffinDiffable eg g)) + -> GenericAgent $ AffinDiffable (ef<>eg) (f^+^g) + (α,β) -> dfblFnValsCombine (\a b -> (a^+^b, arr addV, const mempty)) α β + negateV = case ( linearManifoldWitness :: LinearManifoldWitness v + , dualSpaceWitness :: DualSpaceWitness v ) of + (LinearManifoldWitness _, DualSpaceWitness) -> \case + (GenericAgent (AffinDiffable ef f)) + -> GenericAgent $ AffinDiffable ef (negateV f) + α -> dfblFnValsFunc (\a -> (negateV a, negateV id, const mempty)) α -instance (RealDimension n, LocallyScalable n a) +instance ∀ n a . (RealDimension n, LocallyScalable n a) => Num (DfblFuncValue n a n) where - fromInteger i = point $ fromInteger i - (+) = (^+^) - (*) = dfblFnValsCombine $ + fromInteger = case ( linearManifoldWitness :: LinearManifoldWitness n + , closedScalarWitness :: ClosedScalarWitness n ) of + (LinearManifoldWitness _, ClosedScalarWitness) -> point . fromInteger + (+) = case closedScalarWitness :: ClosedScalarWitness n of + ClosedScalarWitness -> (^+^) + (*) = case ( linearManifoldWitness :: LinearManifoldWitness n + , closedScalarWitness :: ClosedScalarWitness n ) of + (LinearManifoldWitness _, ClosedScalarWitness) -> dfblFnValsCombine $ \a b -> ( a*b , arr $ addV <<< (scale $ a)***(scale $ b) - , unsafe_dev_ε_δ(show a++"*"++show b) sqrt + , unsafe_dev_ε_δ(show a++"*"++show b) (sqrt :: n->n) >>> \d¹₂ -> (d¹₂,d¹₂) -- ε δa δb = (a+δa)·(b+δb) - (a·b + (a·δa + b·δb)) -- = δa·δb -- so choose δa = δb = √ε ) - negate = negateV - abs = dfblFnValsFunc dfblAbs - where dfblAbs a - | a>0 = (a, id, unsafe_dev_ε_δ("abs "++show a) $ \ε -> a + ε/2) - | a<0 = (-a, negateV id, unsafe_dev_ε_δ("abs "++show a) $ \ε -> ε/2 - a) - | otherwise = (0, zeroV, scaleNorm (sqrt 0.5)) - signum = dfblFnValsFunc dfblSgn - where dfblSgn a - | a>0 = (1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ const a) - | a<0 = (-1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ \_ -> -a) - | otherwise = (0, zeroV, const $ spanNorm [1]) + negate = case closedScalarWitness :: ClosedScalarWitness n of + ClosedScalarWitness -> negateV + abs = mkabs linearManifoldWitness closedScalarWitness + where mkabs :: LinearManifoldWitness n -> ClosedScalarWitness n + -> DfblFuncValue n a n -> DfblFuncValue n a n + mkabs (LinearManifoldWitness _) ClosedScalarWitness = dfblFnValsFunc dfblAbs + where dfblAbs a + | a>0 = (a, id, unsafe_dev_ε_δ("abs "++show a) $ \ε -> a + ε/2) + | a<0 = (-a, negateV id, unsafe_dev_ε_δ("abs "++show a) $ \ε -> ε/2 - a) + | otherwise = (0, zeroV, scaleNorm (sqrt 0.5)) + signum = mksgn linearManifoldWitness closedScalarWitness + where mksgn :: LinearManifoldWitness n -> ClosedScalarWitness n + -> DfblFuncValue n a n -> DfblFuncValue n a n + mksgn (LinearManifoldWitness _) ClosedScalarWitness = dfblFnValsFunc dfblSgn + where dfblSgn a + | a>0 = (1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ const a) + | a<0 = (-1, zeroV, unsafe_dev_ε_δ("signum "++show a) $ \_ -> -a) + | otherwise = (0, zeroV, const $ spanNorm [1]) @@ -516,10 +548,13 @@ -- | Important special operator needed to compute intersection of 'Region's. -minDblfuncs :: (LocallyScalable s m, RealDimension s) +minDblfuncs :: ∀ s m . (LocallyScalable s m, RealDimension s) => Differentiable s m s -> Differentiable s m s -> Differentiable s m s -minDblfuncs (Differentiable f) (Differentiable g) = Differentiable h - where h x +minDblfuncs (Differentiable f) (Differentiable g) + = Differentiable $ h linearManifoldWitness closedScalarWitness + where h :: LinearManifoldWitness s -> ClosedScalarWitness s + -> m -> (s, Needle m+>Needle s, LinDevPropag m s) + h (LinearManifoldWitness _) ClosedScalarWitness x | fx < gx = ( fx, jf , \d -> devf d <> devg d <> transformNorm δj @@ -542,9 +577,11 @@ -genericisePreRegion :: (RealDimension s, LocallyScalable s m) +genericisePreRegion :: ∀ s m . (RealDimension s, LocallyScalable s m) => PreRegion s m -> PreRegion s m -genericisePreRegion GlobalRegion = PreRegion $ const 1 +genericisePreRegion GlobalRegion = case ( linearManifoldWitness :: LinearManifoldWitness s + , closedScalarWitness :: ClosedScalarWitness s ) of + (LinearManifoldWitness _, ClosedScalarWitness) -> PreRegion $ const 1 genericisePreRegion (RealSubray PositiveHalfSphere xl) = preRegionToInfFrom' xl genericisePreRegion (RealSubray NegativeHalfSphere xr) = preRegionFromMinInfTo' xr genericisePreRegion r = r @@ -571,13 +608,19 @@ regionProd (Region a₀ ra) (Region b₀ rb) = Region (a₀,b₀) (preRegionProd ra rb) -- | Cartesian product of two pre-regions. -preRegionProd :: (RealDimension s, LocallyScalable s a, LocallyScalable s b) +preRegionProd :: ∀ s a b . (RealDimension s, LocallyScalable s a, LocallyScalable s b) => PreRegion s a -> PreRegion s b -> PreRegion s (a,b) -preRegionProd GlobalRegion GlobalRegion = GlobalRegion -preRegionProd GlobalRegion (PreRegion rb) = PreRegion $ rb . snd -preRegionProd (PreRegion ra) GlobalRegion = PreRegion $ ra . fst -preRegionProd (PreRegion ra) (PreRegion rb) = PreRegion $ minDblfuncs (ra.fst) (rb.snd) -preRegionProd ra rb = preRegionProd (genericisePreRegion ra) (genericisePreRegion rb) +preRegionProd = prp linearManifoldWitness closedScalarWitness + where prp :: LinearManifoldWitness s -> ClosedScalarWitness s + -> PreRegion s a -> PreRegion s b -> PreRegion s (a,b) + prp _ _ GlobalRegion GlobalRegion = GlobalRegion + prp (LinearManifoldWitness _) ClosedScalarWitness GlobalRegion (PreRegion rb) + = PreRegion $ rb . snd + prp (LinearManifoldWitness _) ClosedScalarWitness (PreRegion ra) GlobalRegion + = PreRegion $ ra . fst + prp (LinearManifoldWitness _) ClosedScalarWitness (PreRegion ra) (PreRegion rb) + = PreRegion $ minDblfuncs (ra.fst) (rb.snd) + prp _ _ ra rb = preRegionProd (genericisePreRegion ra) (genericisePreRegion rb) positivePreRegion, negativePreRegion :: (RealDimension s) => PreRegion s s @@ -585,9 +628,13 @@ negativePreRegion = RealSubray NegativeHalfSphere 0 -positivePreRegion', negativePreRegion' :: (RealDimension s) => PreRegion s s -positivePreRegion' = PreRegion $ Differentiable prr - where prr x = ( 1 - 1/xp1 +positivePreRegion', negativePreRegion' :: ∀ s . (RealDimension s) => PreRegion s s +positivePreRegion' = PreRegion . Differentiable + $ prr linearManifoldWitness closedScalarWitness + where prr :: LinearManifoldWitness s -> ClosedScalarWitness s + -> s -> (s, Needle s+>Needle s, LinDevPropag s s) + prr (LinearManifoldWitness _) ClosedScalarWitness + x = ( 1 - 1/xp1 , (1/xp1²) *^ id , unsafe_dev_ε_δ("positivePreRegion@"++show x) δ ) -- ε = (1 − 1/(1+x)) + (-δ · 1/(x+1)²) − (1 − 1/(1+x−δ)) @@ -624,26 +671,42 @@ | otherwise = ε * x / ((1+ε)/x + ε) xp1 = (x+1) xp1² = xp1 ^ 2 -negativePreRegion' = PreRegion $ ppr . ngt - where PreRegion ppr = positivePreRegion' - ngt = actuallyLinearEndo $ negateV id +negativePreRegion' = npr (linearManifoldWitness :: LinearManifoldWitness s) + (closedScalarWitness :: ClosedScalarWitness s) + where npr (LinearManifoldWitness BoundarylessWitness) + (ClosedScalarWitness :: ClosedScalarWitness s) + = PreRegion $ ppr . ngt + where PreRegion ppr = positivePreRegion' + ngt = actuallyLinearEndo $ negateV id preRegionToInfFrom, preRegionFromMinInfTo :: RealDimension s => s -> PreRegion s s preRegionToInfFrom = RealSubray PositiveHalfSphere preRegionFromMinInfTo = RealSubray NegativeHalfSphere -preRegionToInfFrom', preRegionFromMinInfTo' :: RealDimension s => s -> PreRegion s s -preRegionToInfFrom' xs = PreRegion $ ppr . trl - where PreRegion ppr = positivePreRegion' - trl = actuallyAffineEndo (-xs) id -preRegionFromMinInfTo' xe = PreRegion $ ppr . flp - where PreRegion ppr = positivePreRegion' - flp = actuallyAffineEndo xe (negateV id) +preRegionToInfFrom', preRegionFromMinInfTo' :: ∀ s . RealDimension s => s -> PreRegion s s +preRegionToInfFrom' = prif (linearManifoldWitness :: LinearManifoldWitness s) + (closedScalarWitness :: ClosedScalarWitness s) + where prif (LinearManifoldWitness BoundarylessWitness) + (ClosedScalarWitness :: ClosedScalarWitness s) + xs = PreRegion $ ppr . trl + where PreRegion ppr = positivePreRegion' + trl = actuallyAffineEndo (-xs) id +preRegionFromMinInfTo' = prif (linearManifoldWitness :: LinearManifoldWitness s) + (closedScalarWitness :: ClosedScalarWitness s) + where prif (LinearManifoldWitness BoundarylessWitness) + (ClosedScalarWitness :: ClosedScalarWitness s) + xe = PreRegion $ ppr . flp + where PreRegion ppr = positivePreRegion' + flp = actuallyAffineEndo xe (negateV id) -intervalPreRegion :: RealDimension s => (s,s) -> PreRegion s s -intervalPreRegion (lb,rb) = PreRegion $ Differentiable prr +intervalPreRegion :: ∀ s . RealDimension s => (s,s) -> PreRegion s s +intervalPreRegion (lb,rb) = PreRegion . Differentiable + $ prr linearManifoldWitness closedScalarWitness where m = lb + radius; radius = (rb - lb)/2 - prr x = ( 1 - ((x-m)/radius)^2 + prr :: LinearManifoldWitness s -> ClosedScalarWitness s + -> s -> (s, Needle s+>Needle s, LinDevPropag s s) + prr (LinearManifoldWitness _) ClosedScalarWitness + x = ( 1 - ((x-m)/radius)^2 , (2*(m-x)/radius^2) *^ id , unsafe_dev_ε_δ("intervalPreRegion@"++show x) $ (*radius) . sqrt ) @@ -1049,8 +1112,8 @@ -instance ( RealDimension n, WithField n Manifold a - , LocallyScalable n a, SimpleSpace (Needle a) ) +instance ∀ n a . ( RealDimension n, WithField n Manifold a + , LocallyScalable n a, SimpleSpace (Needle a) ) => Floating (RWDfblFuncValue n a n) where pi = point pi @@ -1058,7 +1121,8 @@ $ \x -> let ex = exp x in if ex*2 == ex -- numerical trouble... then if x<0 then ( 0, zeroV, unsafe_dev_ε_δ("exp "++show x) $ \ε -> log ε - x ) - else ( ex, ex*^id, unsafe_dev_ε_δ("exp "++show x) $ \_ -> 1e-300 ) + else ( ex, ex*^id + , unsafe_dev_ε_δ("exp "++show x) $ \_ -> 1e-300 :: n ) else ( ex, ex *^ id, unsafe_dev_ε_δ("exp "++show x) $ \ε -> case acosh(ε/(2*ex) + 1) of δ | δ==δ -> δ @@ -1315,10 +1379,13 @@ -- | Like 'Data.VectorSpace.lerp', but gives a differentiable function -- instead of a Hask one. -lerp_diffable :: ( WithField s LinearManifold m, Atlas m - , HasTrie (ChartIndex m), RealDimension s ) +lerp_diffable :: ∀ m s . ( LinearSpace m, Scalar m ~ s, Atlas m + , HasTrie (ChartIndex m), RealDimension s ) => m -> m -> Differentiable s s m -lerp_diffable a b = actuallyAffine a . arr $ flipBilin scale $ b.-.a +lerp_diffable = case ( linearManifoldWitness :: LinearManifoldWitness m + , dualSpaceWitness :: DualSpaceWitness m ) of + (LinearManifoldWitness BoundarylessWitness, DualSpaceWitness) + -> \a b -> actuallyAffine a . arr $ flipBilin scale $ b.-.a diff --git a/Data/Function/Differentiable/Data.hs b/Data/Function/Differentiable/Data.hs --- a/Data/Function/Differentiable/Data.hs +++ b/Data/Function/Differentiable/Data.hs @@ -91,7 +91,7 @@ -- includes that point) to define a connected subset of a manifold. data PreRegion s m where GlobalRegion :: PreRegion s m - RealSubray :: RealDimension s => S⁰ -> s -> PreRegion s s + RealSubray :: Num' s => S⁰ -> s -> PreRegion s s PreRegion :: (Differentiable s m s) -- A function that is positive at reference point /p/, -- decreases and crosses zero at the region's -- boundaries. (If it goes positive again somewhere diff --git a/Data/Manifold/Atlas.hs b/Data/Manifold/Atlas.hs --- a/Data/Manifold/Atlas.hs +++ b/Data/Manifold/Atlas.hs @@ -9,6 +9,8 @@ -- {-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE EmptyDataDecls, EmptyCase #-} {-# LANGUAGE CPP #-} @@ -25,9 +27,14 @@ import Data.Void import Data.VectorSpace.Free +import Math.LinearMap.Category import Control.Arrow +import Data.MemoTrie (HasTrie) + +import qualified Linear.Affine as LinAff + class Semimanifold m => Atlas m where type ChartIndex m :: * chartReferencePoint :: ChartIndex m -> m @@ -49,6 +56,8 @@ VectorSpaceAtlas(Num s, V2 s) VectorSpaceAtlas(Num s, V3 s) VectorSpaceAtlas(Num s, V4 s) +VectorSpaceAtlas((LinearSpace v, Scalar v ~ s, TensorSpace w, Scalar w ~ s), LinearMap s v w) +VectorSpaceAtlas((TensorSpace v, Scalar v ~ s, TensorSpace w, Scalar w ~ s), Tensor s v w) instance (Atlas x, Atlas y) => Atlas (x,y) where type ChartIndex (x,y) = (ChartIndex x, ChartIndex y) @@ -78,3 +87,23 @@ interiorChartReferencePoint _ NegativeHalfSphere = S² pi 0 lookupAtlas (S² ϑ _) | ϑ Atlas (LinAff.Point a n) where + type ChartIndex (LinAff.Point a n) = () + interiorChartReferencePoint _ () = LinAff.P zeroV + lookupAtlas _ = () + + + +-- | The 'AffineSpace' class plus manifold constraints. +type AffineManifold m = ( Atlas m, Manifold m, AffineSpace m + , Needle m ~ Diff m, HasTrie (ChartIndex m) ) + +-- | An euclidean space is a real affine space whose tangent space is a Hilbert space. +type EuclidSpace x = ( AffineManifold x, InnerSpace (Diff x) + , DualVector (Diff x) ~ Diff x, Floating (Scalar (Diff x)) ) + +euclideanMetric :: EuclidSpace x => proxy x -> Metric x +euclideanMetric _ = euclideanNorm + diff --git a/Data/Manifold/Cone.hs b/Data/Manifold/Cone.hs --- a/Data/Manifold/Cone.hs +++ b/Data/Manifold/Cone.hs @@ -143,7 +143,7 @@ projCD¹ToCℝay (CD¹ h m) = Cℝay (bijectIntvtoℝplus h) m -stiefel1Project :: LinearManifold v => +stiefel1Project :: LinearSpace v => DualVector v -- ^ Must be nonzero. -> Stiefel1 v stiefel1Project = Stiefel1 diff --git a/Data/Manifold/DifferentialEquation.hs b/Data/Manifold/DifferentialEquation.hs --- a/Data/Manifold/DifferentialEquation.hs +++ b/Data/Manifold/DifferentialEquation.hs @@ -33,10 +33,11 @@ module Data.Manifold.DifferentialEquation ( -- * Formulating simple differential eqns. - DifferentialEqn + DifferentialEqn, ODE + , constLinearDEqn , constLinearODE , constLinearPDE - , filterDEqnSolution_static, iterateFilterDEqn_static + , iterateFilterDEqn_static -- * Cost functions for error bounds , maxDeviationsGoal , uncertaintyGoal @@ -62,7 +63,10 @@ import Data.Function.Differentiable.Data import Data.Manifold.TreeCover import Data.Manifold.Web +import Data.Manifold.Atlas +import Data.Embedding + import qualified Data.List as List import qualified Prelude as Hask hiding(foldl, sum, sequence) @@ -79,39 +83,92 @@ import Data.Traversable.Constrained (Traversable, traverse) +-- | An ordinary differential equation is one that does not need any a-priori +-- partial derivatives to compute the derivative for integration in some +-- propagation direction. Classically, ODEs are usually understood as +-- @DifferentialEquation ℝ ℝ⁰ y@, but actually @x@ can at least +-- be an arbitrary one-dimensional space (i.e. basically real intervals or 'S¹'). +-- In these cases, there is always only one partial derivative: that which we +-- integrate over, in the only possible direction for propagation. +type ODE x y = DifferentialEqn x ℝ⁰ y -constLinearODE :: ∀ x y . ( WithField ℝ LinearManifold x, SimpleSpace x - , WithField ℝ LinearManifold y, SimpleSpace y ) - => ((x +> y) +> y) -> DifferentialEqn x y -constLinearODE = case ( dualSpaceWitness :: DualNeedleWitness x - , dualSpaceWitness :: DualNeedleWitness y ) of - (DualSpaceWitness, DualSpaceWitness) -> \bwt' -> +constLinearDEqn :: ∀ x y ð . ( SimpleSpace x + , SimpleSpace y, AffineManifold y + , SimpleSpace ð, AffineManifold ð + , Scalar x ~ ℝ, Scalar y ~ ℝ, Scalar ð ~ ℝ ) + => ((y,ð) +> (x +> y)) -> ((x +> y) +> (y,ð)) -> DifferentialEqn x ð y +constLinearDEqn = case ( linearManifoldWitness :: LinearManifoldWitness x + , dualSpaceWitness :: DualSpaceWitness x + , linearManifoldWitness :: LinearManifoldWitness y + , dualSpaceWitness :: DualSpaceWitness y + , linearManifoldWitness :: LinearManifoldWitness ð + , dualSpaceWitness :: DualSpaceWitness ð ) of + ( LinearManifoldWitness BoundarylessWitness, DualSpaceWitness + ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness + ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ) -> \bwt'inv bwt' -> + \(Shade (_x,y) δxy) -> LocalDifferentialEqn + { _predictDerivatives + = \(Shade' ð δð) -> + let j = bwt'inv $ (y,ð) + δj = bwt' `transformNorm` + sumSubspaceNorms (transformNorm (zeroV&&&id) $ dualNorm δxy) δð + in return $ Shade' j δj + , _rescanDerivatives + = \shy shjApriori _ + -> ( mixShade's $ shy + :| [ projectShade + (Embedding (arr bwt'inv <<< id&&&zeroV) + (arr bwt' >>> fst)) + shjApriori ] + , return $ projectShade + (Embedding (arr bwt'inv <<< zeroV&&&id) + (arr bwt' >>> snd)) + shjApriori + ) + } + +constLinearODE :: ∀ x y . ( SimpleSpace x, Scalar x ~ ℝ, SimpleSpace y, Scalar y ~ ℝ ) + => ((x +> y) +> y) -> ODE x y +constLinearODE = case ( linearManifoldWitness :: LinearManifoldWitness x + , dualSpaceWitness :: DualSpaceWitness x + , linearManifoldWitness :: LinearManifoldWitness y + , dualSpaceWitness :: DualSpaceWitness y ) of + ( LinearManifoldWitness BoundarylessWitness, DualSpaceWitness + ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ) -> \bwt' -> let bwt'inv = (bwt'\$) in \(Shade (_x,y) δxy) -> LocalDifferentialEqn (let j = bwt'inv y δj = (bwt'>>>zeroV&&&id) `transformNorm` dualNorm δxy - in return $ Shade' j δj ) - (\_ -> pure ) + in \_ -> return $ Shade' j δj ) + (\shy _ _ -> (pure shy, Just $ Shade' Origin mempty) ) -constLinearPDE :: ∀ x y y' . - ( WithField ℝ LinearManifold x, SimpleSpace x - , WithField ℝ LinearManifold y, SimpleSpace y, FiniteFreeSpace y - , WithField ℝ LinearManifold y', SimpleSpace y' ) - => ((x +> (y,y')) +> (y, y')) -> DifferentialEqn x (y,y') -constLinearPDE = undefined{-case ( dualSpaceWitness :: DualNeedleWitness x - , dualSpaceWitness :: DualNeedleWitness y - , dualSpaceWitness :: DualSpaceWitness y' ) of - (DualSpaceWitness, DualSpaceWitness, DualSpaceWitness) -> \bwt' -> - let bwt'inv = (bwt'\$) - in \(Shade (_x,(y,y')) δxy) (Shade' jApriori σjApriori) - -> let j = bwt'inv $ (zeroV,y') - δj = (bwt'>>>zeroV&&&id) - `transformNorm` dualNorm δxy - (_,y'Apriori) = bwt' $ jApriori - Norm δy' = (arr $ LinearFunction bwt'inv . (zeroV&&&id)) - `transformNorm` σjApriori - in (Shade' (y,y'Apriori) . Norm $ zeroV *** δy' , ) - <$> mixShade's (Shade' jApriori σjApriori :| [Shade' j δj])-} +constLinearPDE :: ∀ x y ð . + ( WithField ℝ SimpleSpace x + , WithField ℝ SimpleSpace y + , WithField ℝ SimpleSpace ð, AffineManifold ð ) + => ((x +> y) +> ð) -> (ð +> (x +> y)) -> DifferentialEqn x ð y +constLinearPDE = case ( linearManifoldWitness :: LinearManifoldWitness x + , dualSpaceWitness :: DualSpaceWitness x + , linearManifoldWitness :: LinearManifoldWitness y + , dualSpaceWitness :: DualSpaceWitness y + , linearManifoldWitness :: LinearManifoldWitness ð + , dualSpaceWitness :: DualSpaceWitness ð ) of + ( LinearManifoldWitness BoundarylessWitness, DualSpaceWitness + ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness + ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ) + -> \bwt' bwt'inv (Shade (_x,y) δxy) + -> LocalDifferentialEqn + { _predictDerivatives + = \(Shade' ð δð) -> + let j = bwt'inv $ ð + δj = bwt' `transformNorm` δð + in return $ Shade' j δj + , _rescanDerivatives + = \shy shjApriori _ + -> ( return shy + , return $ projectShade (Embedding (arr bwt'inv) (arr bwt')) shjApriori + ) + } -- | A function that variates, relatively speaking, most strongly -- for arguments around 1. In the zero-limit it approaches a constant diff --git a/Data/Manifold/Function/Quadratic.hs b/Data/Manifold/Function/Quadratic.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/Function/Quadratic.hs @@ -0,0 +1,124 @@ +-- | +-- Module : Data.Manifold.Function.Quadratic +-- Copyright : (c) Justus Sagemüller 2016 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- + +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE LiberalTypeSynonyms #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE ScopedTypeVariables #-} + + +module Data.Manifold.Function.Quadratic ( + Quadratic(..), evalQuadratic + ) where + + + +import Data.Semigroup +import qualified Data.List.NonEmpty as NE + +import Data.MemoTrie +import Data.VectorSpace +import Data.AffineSpace +import Data.Tagged +import Data.Manifold.PseudoAffine +import Data.Manifold.Atlas +import Data.Manifold.Riemannian +import Data.Function.Affine + +import Prelude hiding (id, ($), fmap, fst) +import Control.Category.Constrained.Prelude (id, ($), fmap, fst) +import Control.Arrow.Constrained ((>>>), (&&&), (***), second) + +import Math.LinearMap.Category + + + +data Quadratic s d c where + Quadratic :: (ChartIndex d :->: ( c, ( LinearMap s (Needle d) (Needle c) + , LinearMap s (SymmetricTensor s (Needle d)) + (Needle c) ) ) + ) -> Quadratic s d c + +affineQuadratic :: (WithField s AffineManifold d, WithField s AffineManifold c) + => Affine s d c -> Quadratic s d c +affineQuadratic (Affine f) = Quadratic . trie + $ untrie f >>> second (id &&& const zeroV) + +instance ( Atlas x, HasTrie (ChartIndex x), LinearSpace (Needle x), Scalar (Needle x) ~ s + , Manifold y, Scalar (Needle y) ~ s ) + => Semimanifold (Quadratic s x y) where + type Needle (Quadratic s x y) = Quadratic s x (Needle y) + toInterior = pure + fromInterior = id + (.+~^) = case ( semimanifoldWitness :: SemimanifoldWitness y + , boundarylessWitness :: BoundarylessWitness y ) of + (SemimanifoldWitness _, BoundarylessWitness) -> \(Quadratic f) (Quadratic g) + -> Quadratic . trie $ \ix -> case (untrie f ix, untrie g ix) of + ((fx₀,f'), (gx₀,g')) -> (fx₀.+~^gx₀, f'^+^g') + translateP = Tagged (.+~^) + semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness y of + SemimanifoldWitness _ -> SemimanifoldWitness BoundarylessWitness +instance ( Atlas x, HasTrie (ChartIndex x), LinearSpace (Needle x), Scalar (Needle x) ~ s + , Manifold y, Scalar (Needle y) ~ s ) + => PseudoAffine (Quadratic s x y) where + (.-~!) = case ( semimanifoldWitness :: SemimanifoldWitness y + , boundarylessWitness :: BoundarylessWitness y ) of + (SemimanifoldWitness _, BoundarylessWitness) -> \(Quadratic f) (Quadratic g) + -> Quadratic . trie $ \ix -> case (untrie f ix, untrie g ix) of + ((fx₀,f'), (gx₀,g')) -> (fx₀.-~!gx₀, f'^-^g') + pseudoAffineWitness = case semimanifoldWitness :: SemimanifoldWitness y of + SemimanifoldWitness _ -> PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) +instance ( Atlas x, HasTrie (ChartIndex x), LinearSpace (Needle x), Scalar (Needle x) ~ s + , Manifold y, Scalar (Needle y) ~ s ) + => AffineSpace (Quadratic s x y) where + type Diff (Quadratic s x y) = Quadratic s x (Needle y) + (.+^) = (.+~^); (.-.) = (.-~!) +instance ( Atlas x, HasTrie (ChartIndex x), LinearSpace (Needle x), Scalar (Needle x) ~ s + , LinearSpace y, Scalar y ~ s, Num' s ) + => AdditiveGroup (Quadratic s x y) where + zeroV = case linearManifoldWitness :: LinearManifoldWitness y of + LinearManifoldWitness _ -> Quadratic . trie $ const (zeroV, zeroV) + (^+^) = case ( linearManifoldWitness :: LinearManifoldWitness y + , dualSpaceWitness :: DualSpaceWitness y ) of + (LinearManifoldWitness BoundarylessWitness, DualSpaceWitness) -> (.+~^) + negateV = case linearManifoldWitness :: LinearManifoldWitness y of + LinearManifoldWitness _ -> \(Quadratic f) -> Quadratic . trie $ + untrie f >>> negateV***negateV +instance ( Atlas x, HasTrie (ChartIndex x), LinearSpace (Needle x), Scalar (Needle x) ~ s + , LinearSpace y, Scalar y ~ s, Num' s ) + => VectorSpace (Quadratic s x y) where + type Scalar (Quadratic s x y) = s + (*^) = case linearManifoldWitness :: LinearManifoldWitness y of + LinearManifoldWitness _ -> \μ (Quadratic f) -> Quadratic . trie $ + untrie f >>> (μ*^)***(μ*^) + +evalQuadratic :: ∀ s x y . ( Manifold x, Atlas x, HasTrie (ChartIndex x) + , Manifold y + , s ~ Scalar (Needle x), s ~ Scalar (Needle y) ) + => Quadratic s x y -> x + -> (y, ( LinearMap s (Needle x) (Needle y) + , LinearMap s (SymmetricTensor s (Needle x)) (Needle y) )) +evalQuadratic = ea (boundarylessWitness, boundarylessWitness) + where ea :: (BoundarylessWitness x, BoundarylessWitness y) + -> Quadratic s x y -> x -> (y, ( LinearMap s (Needle x) (Needle y) + , LinearMap s (SymmetricTensor s (Needle x)) (Needle y) )) + ea (BoundarylessWitness, BoundarylessWitness) + (Quadratic f) x = ( fx₀.+~^(ðx'f₀ $ v).+~^(ð²x'f $ squareV v) + , ( ðx'f₀ ^+^ 2*^((currySymBilin $ ð²x'f) $ v) + , ð²x'f + ) ) + where Just v = x .-~. chartReferencePoint chIx + chIx = lookupAtlas x + (fx₀, (ðx'f₀, ð²x'f)) = untrie f chIx + + diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs --- a/Data/Manifold/PseudoAffine.hs +++ b/Data/Manifold/PseudoAffine.hs @@ -57,7 +57,7 @@ -- ** Needles , Local(..) -- ** Metrics - , Metric, Metric', euclideanMetric + , Metric, Metric' , RieMetric, RieMetric' -- ** Constraints , SemimanifoldWitness(..) @@ -65,14 +65,10 @@ , BoundarylessWitness(..) , boundarylessWitness , DualNeedleWitness - , RealDimension, AffineManifold - , LinearManifold , WithField - , HilbertManifold - , EuclidSpace , LocallyScalable -- ** Local functions - , LocalLinear, LocalAffine + , LocalLinear, LocalBilinear, LocalAffine -- * Misc , alerpB, palerp, palerpB, LocallyCoercible(..), CanonicalDiffeomorphism(..) , ImpliesMetric(..), coerceMetric, coerceMetric' @@ -93,6 +89,7 @@ import qualified Linear.Affine as LinAff import Data.Embedding import Data.LinearMap +import Data.VectorSpace.Free import Math.LinearMap.Category import Data.AffineSpace import Data.Tagged @@ -108,6 +105,8 @@ import Control.Monad.Constrained import Data.Foldable.Constrained +import Control.Lens (Lens', lens, (^.), (&), (%~), (.~)) + import GHC.Exts (Constraint) @@ -147,6 +146,8 @@ => p (x,ξ) -> CanonicalDiffeomorphism (Interior x) (Interior ξ) interiorLocalCoercion _ = CanonicalDiffeomorphism +type NumPrime n = (Num' n, Eq n) + #define identityCoercion(c,t) \ instance (c) => LocallyCoercible (t) (t) where { \ locallyTrivialDiffeomorphism = id; \ @@ -154,16 +155,16 @@ coerceNeedle' _ = id; \ oppositeLocalCoercion = CanonicalDiffeomorphism; \ interiorLocalCoercion _ = CanonicalDiffeomorphism } -identityCoercion(NumberManifold s, ZeroDim s) -identityCoercion(NumberManifold s, V0 s) +identityCoercion(NumPrime s, ZeroDim s) +identityCoercion(NumPrime s, V0 s) identityCoercion((), ℝ) -identityCoercion(NumberManifold s, V1 s) +identityCoercion(NumPrime s, V1 s) identityCoercion((), (ℝ,ℝ)) -identityCoercion(NumberManifold s, V2 s) +identityCoercion(NumPrime s, V2 s) identityCoercion((), (ℝ,(ℝ,ℝ))) identityCoercion((), ((ℝ,ℝ),ℝ)) -identityCoercion(NumberManifold s, V3 s) -identityCoercion(NumberManifold s, V4 s) +identityCoercion(NumPrime s, V3 s) +identityCoercion(NumPrime s, V4 s) data CanonicalDiffeomorphism a b where @@ -180,15 +181,14 @@ , Num' s ) type LocalLinear x y = LinearMap (Scalar (Needle x)) (Needle x) (Needle y) -type LocalAffine x y = (Needle y, LocalLinear x y) +type LocalBilinear x y = LinearMap (Scalar (Needle x)) + (SymmetricTensor (Scalar (Needle x)) (Needle x)) + (Needle y) --- | Basically just an “updated” version of the 'VectorSpace' class. --- Every vector space is a manifold, this constraint makes it explicit. -type LinearManifold x = ( AffineManifold x, Needle x ~ x, LSpace x ) -type LinearManifold' x = ( PseudoAffine x, AffineSpace x, Diff x ~ x - , Interior x ~ x, Needle x ~ x, LSpace x ) +type LocalAffine x y = (Needle y, LocalLinear x y) + -- | Require some constraint on a manifold, and also fix the type of the manifold's -- underlying field. For example, @WithField ℝ 'HilbertManifold' v@ constrains -- @v@ to be a real (i.e., 'Double'-) Hilbert space. @@ -198,34 +198,7 @@ -- applied, for @type@ constraints this is by default not allowed). type WithField s c x = ( c x, s ~ Scalar (Needle x), s ~ Scalar (Needle' x) ) --- | The 'RealFloat' class plus manifold constraints. -type RealDimension r = ( PseudoAffine r, Interior r ~ r, Needle r ~ r, r ~ ℝ) --- | The 'AffineSpace' class plus manifold constraints. -type AffineManifold m = ( PseudoAffine m, Interior m ~ m, AffineSpace m - , Needle m ~ Diff m, LinearManifold' (Diff m) ) - --- | A Hilbert space is a /complete/ inner product space. Being a vector space, it is --- also a manifold. --- --- (Stricly speaking, that doesn't have much to do with the completeness criterion; --- but since 'Manifold's are at the moment confined to finite dimension, they are in --- fact (trivially) complete.) -type HilbertManifold x = ( LinearManifold x, InnerSpace x - , Interior x ~ x, Needle x ~ x, DualVector x ~ x - , Floating (Scalar x) ) - --- | An euclidean space is a real affine space whose tangent space is a Hilbert space. -type EuclidSpace x = ( AffineManifold x, InnerSpace (Diff x) - , DualVector (Diff x) ~ Diff x, Floating (Scalar (Diff x)) ) - -type NumberManifold n = ( Num' n, Manifold n, Interior n ~ n, Needle n ~ n - , LSpace n, DualVector n ~ n, Scalar n ~ n ) - -euclideanMetric :: EuclidSpace x => proxy x -> Metric x -euclideanMetric _ = euclideanNorm - - -- | A co-needle can be understood as a “paper stack”, with which you can measure -- the length that a needle reaches in a given direction by counting the number -- of holes punched through them. @@ -289,11 +262,11 @@ deriveAffine(KnownNat n, FreeVect n ℝ) -instance (NumberManifold s) => LocallyCoercible (ZeroDim s) (V0 s) where +instance (NumPrime s) => LocallyCoercible (ZeroDim s) (V0 s) where locallyTrivialDiffeomorphism Origin = V0 coerceNeedle _ = LinearFunction $ \Origin -> V0 coerceNeedle' _ = LinearFunction $ \Origin -> V0 -instance (NumberManifold s) => LocallyCoercible (V0 s) (ZeroDim s) where +instance (NumPrime s) => LocallyCoercible (V0 s) (ZeroDim s) where locallyTrivialDiffeomorphism V0 = Origin coerceNeedle _ = LinearFunction $ \V0 -> Origin coerceNeedle' _ = LinearFunction $ \V0 -> Origin @@ -376,13 +349,15 @@ -> CanonicalDiffeomorphism -instance LinearManifold (a n) => Semimanifold (LinAff.Point a n) where +instance (LinearSpace (a n), Needle (a n) ~ a n, Interior (a n) ~ a n) + => Semimanifold (LinAff.Point a n) where type Needle (LinAff.Point a n) = a n fromInterior = id toInterior = pure LinAff.P v .+~^ w = LinAff.P $ v ^+^ w translateP = Tagged $ \(LinAff.P v) w -> LinAff.P $ v ^+^ w -instance LinearManifold (a n) => PseudoAffine (LinAff.Point a n) where +instance (LinearSpace (a n), Needle (a n) ~ a n, Interior (a n) ~ a n) + => PseudoAffine (LinAff.Point a n) where LinAff.P v .-~. LinAff.P w = return $ v ^-^ w diff --git a/Data/Manifold/Riemannian.hs b/Data/Manifold/Riemannian.hs --- a/Data/Manifold/Riemannian.hs +++ b/Data/Manifold/Riemannian.hs @@ -41,12 +41,15 @@ {-# LANGUAGE LiberalTypeSynonyms #-} {-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} +{-# LANGUAGE DefaultSignatures #-} module Data.Manifold.Riemannian where import Data.Maybe +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NE import qualified Data.Vector as Arr import Data.VectorSpace @@ -59,6 +62,7 @@ import Data.Manifold.Types.Primitive ((^), empty, embed, coEmbed) import Data.Manifold.Types.Stiefel import Data.Manifold.PseudoAffine +import Data.Manifold.Atlas (AffineManifold) import Data.CoNat @@ -76,6 +80,9 @@ import Data.Foldable.Constrained +data GeodesicWitness x where + GeodesicWitness :: Geodesic (Interior x) + => SemimanifoldWitness x -> GeodesicWitness x class Semimanifold x => Geodesic x where geodesicBetween :: @@ -85,7 +92,13 @@ -- If the two points are actually connected by a path... -> Maybe (D¹ -> x) -- ^ ...then this is the interpolation function. Attention: -- the type will change to 'Differentiable' in the future. + geodesicWitness :: GeodesicWitness x + default geodesicWitness :: Geodesic (Interior x) => GeodesicWitness x + geodesicWitness = GeodesicWitness semimanifoldWitness + middleBetween :: x -> x -> Maybe x + middleBetween p₀ p₁ = ($ D¹ 0) <$> geodesicBetween p₀ p₁ + interpolate :: (Geodesic x, IntervalLike i) => x -> x -> Maybe (i -> x) interpolate a b = (. toClosedInterval) <$> geodesicBetween a b @@ -94,36 +107,55 @@ #define deriveAffineGD(x) \ instance Geodesic x where { \ - geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹ \ + geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹; \ + middleBetween a b = return $ alerp a b (1/2) \ } deriveAffineGD (ℝ) -instance Geodesic (ZeroDim ℝ) where +instance Geodesic (ZeroDim s) where geodesicBetween Origin Origin = return $ \_ -> Origin + middleBetween Origin Origin = return Origin -instance (Geodesic a, Geodesic b) => Geodesic (a,b) where +instance ∀ a b . (Geodesic a, Geodesic b) => Geodesic (a,b) where geodesicBetween (a,b) (α,β) = liftA2 (&&&) (geodesicBetween a α) (geodesicBetween b β) + geodesicWitness = case ( geodesicWitness :: GeodesicWitness a + , geodesicWitness :: GeodesicWitness b ) of + (GeodesicWitness _, GeodesicWitness _) -> GeodesicWitness semimanifoldWitness + middleBetween (a,b) (α,β) = fzip (middleBetween a α, middleBetween b β) -instance (Geodesic a, Geodesic b, Geodesic c) => Geodesic (a,b,c) where +instance ∀ a b c . (Geodesic a, Geodesic b, Geodesic c) => Geodesic (a,b,c) where geodesicBetween (a,b,c) (α,β,γ) = liftA3 (\ia ib ic t -> (ia t, ib t, ic t)) (geodesicBetween a α) (geodesicBetween b β) (geodesicBetween c γ) + geodesicWitness = case ( geodesicWitness :: GeodesicWitness a + , geodesicWitness :: GeodesicWitness b + , geodesicWitness :: GeodesicWitness c ) of + (GeodesicWitness _, GeodesicWitness _, GeodesicWitness _) + -> GeodesicWitness semimanifoldWitness -- instance (KnownNat n) => Geodesic (FreeVect n ℝ) where -- geodesicBetween (FreeVect v) (FreeVect w) -- = return $ \(D¹ t) -> let μv = (1-t)/2; μw = (t+1)/2 -- in FreeVect $ Arr.zipWith (\vi wi -> μv*vi + μw*wi) v w -instance (Geodesic v, FiniteFreeSpace v, WithField ℝ HilbertManifold v) +instance ∀ v . ( Geodesic v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v) + , LinearSpace v, Scalar v ~ ℝ, Geodesic (DualVector v) + , InnerSpace (DualVector v) ) => Geodesic (Stiefel1 v) where - geodesicBetween (Stiefel1 p') (Stiefel1 q') - = (\f -> \(D¹ t) -> Stiefel1 . f . D¹ $ g * tan (ϑ*t)) + geodesicBetween = gb dualSpaceWitness + where gb :: DualSpaceWitness v -> Stiefel1 v -> Stiefel1 v -> Maybe (D¹ -> Stiefel1 v) + gb DualSpaceWitness (Stiefel1 p') (Stiefel1 q') + = (\f -> \(D¹ t) -> Stiefel1 . f . D¹ $ g * tan (ϑ*t)) <$> geodesicBetween p q - where p = normalized p'; q = normalized q' - l = magnitude $ p^-^q - ϑ = asin $ l/2 - g = sqrt $ 4/l^2 - 1 + where p = normalized p'; q = normalized q' + l = magnitude $ p^-^q + ϑ = asin $ l/2 + g = sqrt $ 4/l^2 - 1 + middleBetween = gb dualSpaceWitness + where gb :: DualSpaceWitness v -> Stiefel1 v -> Stiefel1 v -> Maybe (Stiefel1 v) + gb DualSpaceWitness (Stiefel1 p) (Stiefel1 q) + = Stiefel1 <$> middleBetween (normalized p) (normalized q) instance Geodesic S⁰ where @@ -138,6 +170,10 @@ <$> geodesicBetween (pi-φ) (-ϕ-pi) | otherwise = (>>> S¹ . \ψ -> signum ψ*pi - ψ) <$> geodesicBetween (-pi-φ) (pi-ϕ) + middleBetween (S¹ φ) (S¹ ϕ) + | abs (φ-ϕ) < pi = S¹ <$> middleBetween φ ϕ + | φ > 0 = S¹ <$> middleBetween (pi-φ) (-ϕ-pi) + | otherwise = S¹ <$> middleBetween (-pi-φ) (pi-ϕ) -- instance Geodesic (Cℝay S⁰) where @@ -198,6 +234,15 @@ deriveAffineGD (ℝ³) deriveAffineGD (ℝ⁴) +instance (TensorSpace v, Scalar v ~ ℝ, TensorSpace w, Scalar w ~ ℝ) + => Geodesic (Tensor ℝ v w) where + geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹ +instance (LinearSpace v, Scalar v ~ ℝ, TensorSpace w, Scalar w ~ ℝ) + => Geodesic (LinearMap ℝ v w) where + geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹ +instance (TensorSpace v, Scalar v ~ ℝ, TensorSpace w, Scalar w ~ ℝ) + => Geodesic (LinearFunction ℝ v w) where + geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹ -- | One-dimensional manifolds, whose closure is homeomorpic to the unit interval. @@ -232,5 +277,21 @@ -middleBetween :: Geodesic m => m -> m -> Maybe m -middleBetween p₀ p₁ = ($ D¹ 0) <$> geodesicBetween p₀ p₁ + +pointsBarycenter :: Geodesic m => NonEmpty m -> Maybe m +pointsBarycenter (p:|[]) = Just p +pointsBarycenter ps = case ( pointsBarycenter (NE.fromList group₀) + , pointsBarycenter (NE.fromList group₁) ) of + (Just bc₀, Just bc₁) + | δn == 0 -> middleBetween bc₀ bc₁ + | otherwise -> ($ D¹ (fromIntegral δn/fromIntegral ntot)) + <$> geodesicBetween bc₀ bc₁ + _ -> Nothing + where psl = Hask.toList ps + (group₀, group₁) = splitAt nl psl + ntot = length psl + nl = ntot`quot`2 + δn = ntot - nl*2 + + +type FlatSpace x = (AffineManifold x, Geodesic x, SimpleSpace x) diff --git a/Data/Manifold/Shade.hs b/Data/Manifold/Shade.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/Shade.hs @@ -0,0 +1,1218 @@ +-- | +-- Module : Data.Manifold.Shade +-- Copyright : (c) Justus Sagemüller 2016 +-- License : GPL v3 +-- +-- Maintainer : (@) jsagemue $ uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ParallelListComp #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE PatternSynonyms #-} +{-# LANGUAGE ViewPatterns #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE LiberalTypeSynonyms #-} +{-# LANGUAGE DefaultSignatures #-} +{-# LANGUAGE TemplateHaskell #-} + + +module Data.Manifold.Shade ( + -- * Shades + Shade(..), pattern(:±), Shade'(..), (|±|), IsShade + -- ** Lenses + , shadeCtr, shadeExpanse, shadeNarrowness + -- ** Construction + , fullShade, fullShade', pointsShades, pointsShade's + , pointsCovers, pointsCover's, coverAllAround + -- ** Evaluation + , occlusion, prettyShowsPrecShade', prettyShowShade', LtdErrorShow + -- ** Misc + , factoriseShade, orthoShades, (✠), intersectShade's, linIsoTransformShade + , embedShade, projectShade + , Refinable, subShade', refineShade', convolveShade', coerceShade + , mixShade's, dualShade + -- * Misc + -- ** Shades + , shadesMerge, pointsShades', pseudoECM, convolveMetric + , WithAny(..), shadeWithAny, shadeWithoutAnything + -- ** Local data fit models + , estimateLocalJacobian, estimateLocalHessian, QuadraticModel(..) + -- ** Differential equations + , DifferentialEqn, LocalDifferentialEqn(..) + , propagateDEqnSolution_loc, LocalDataPropPlan(..) + -- ** Range interpolation + , rangeOnGeodesic, rangeWithinVertices + ) where + + +import Data.List hiding (filter, all, elem, sum, foldr1) +import Data.Maybe +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NE +import Data.Semigroup +import Control.DeepSeq +import Data.MemoTrie + +import Data.VectorSpace +import Data.AffineSpace +import Math.LinearMap.Category +import Data.Tagged +import Linear (_x,_y,_z,_w) + +import Data.Manifold.Types +import Data.Manifold.Types.Primitive ((^)) +import Data.Manifold.PseudoAffine +import Data.Manifold.Riemannian +import Data.Manifold.Atlas +import Data.Function.Affine +import Data.Manifold.Function.Quadratic + +import Data.Embedding + +import Control.Lens (Lens', (^.), view, _1, _2, mapping, (&)) +import Control.Lens.TH + +import qualified Prelude as Hask hiding(foldl, sum, sequence) +import qualified Control.Applicative as Hask +import qualified Data.Foldable as Hask +import Data.Foldable (all, elem, toList, sum, foldr1) + +import Control.Category.Constrained.Prelude hiding + ((^), all, elem, sum, forM, Foldable(..), foldr1, Traversable, traverse) +import Control.Arrow.Constrained +import Control.Monad.Constrained hiding (forM) + +import GHC.Generics (Generic) + +import Text.Show.Number + + +-- | A 'Shade' is a very crude description of a region within a manifold. It +-- can be interpreted as either an ellipsoid shape, or as the Gaussian peak +-- of a normal distribution (use +-- for actually sampling from that distribution). +-- +-- For a /precise/ description of an arbitrarily-shaped connected subset of a manifold, +-- there is 'Region', whose implementation is vastly more complex. +data Shade x = Shade { _shadeCtr :: !(Interior x) + , _shadeExpanse :: !(Metric' x) } +deriving instance (Show (Interior x), Show (Metric' x), WithField ℝ PseudoAffine x) + => Show (Shade x) + +-- | A “co-shade” can describe ellipsoid regions as well, but unlike +-- 'Shade' it can be unlimited / infinitely wide in some directions. +-- It does OTOH need to have nonzero thickness, which 'Shade' needs not. +data Shade' x = Shade' { _shade'Ctr :: !(Interior x) + , _shade'Narrowness :: !(Metric x) } + +data LocalDifferentialEqn x ð y = LocalDifferentialEqn { + _predictDerivatives :: Shade' ð -> Maybe (Shade' (LocalLinear x y)) + , _rescanDerivatives :: Shade' y -> Shade' (LocalLinear x y) + -> Shade' (LocalBilinear x y) + -> (Maybe (Shade' y), Maybe (Shade' ð)) + } +makeLenses ''LocalDifferentialEqn + +type DifferentialEqn x ð y = Shade (x,y) -> LocalDifferentialEqn x ð y + +data LocalDataPropPlan x ym yr = LocalDataPropPlan + { _sourcePosition :: !(Interior x) + , _targetPosOffset :: !(Needle x) + , _sourceData, _targetAPrioriData :: !ym + , _relatedData :: [(Needle x, yr)] + } +deriving instance (Show (Interior x), Show ym, Show yr, Show (Needle x)) + => Show (LocalDataPropPlan x ym yr) + +makeLenses ''LocalDataPropPlan + + +class IsShade shade where +-- type (*) shade :: *->* + -- | Access the center of a 'Shade' or a 'Shade''. + shadeCtr :: Lens' (shade x) (Interior x) +-- -- | Convert between 'Shade' and 'Shade' (which must be neither singular nor infinite). +-- unsafeDualShade :: WithField ℝ Manifold x => shade x -> shade* x + -- | Check the statistical likelihood-density of a point being within a shade. + -- This is taken as a normal distribution. + occlusion :: ( PseudoAffine x, SimpleSpace (Needle x) + , s ~ (Scalar (Needle x)), RealFloat' s ) + => shade x -> x -> s + factoriseShade :: ( PseudoAffine x, SimpleSpace (Needle x) + , PseudoAffine y, SimpleSpace (Needle y) + , Scalar (Needle x) ~ Scalar (Needle y) ) + => shade (x,y) -> (shade x, shade y) + coerceShade :: (Manifold x, Manifold y, LocallyCoercible x y) => shade x -> shade y + -- | ASCII version of '✠'. + orthoShades :: ( PseudoAffine x, SimpleSpace (Needle x) + , PseudoAffine y, SimpleSpace (Needle y) + , Scalar (Needle x) ~ Scalar (Needle y) ) + => shade x -> shade y -> shade (x,y) + linIsoTransformShade :: ( SimpleSpace x, SimpleSpace y, Scalar x ~ Scalar y + , Num' (Scalar x) ) + => (x+>y) -> shade x -> shade y + -- | Squash a shade down into a lower dimensional space. + projectShade :: ( Semimanifold x, Semimanifold y + , Object (Affine s) (Interior x), Object (Affine s) (Interior y) + , SemiInner (Needle x), SemiInner (Needle y) ) + => Embedding (Affine s) (Interior x) (Interior y) + -> shade y -> shade x + -- | Include a shade in a higher-dimensional space. Notice that this behaves + -- fundamentally different for 'Shade' and 'Shade''. For 'Shade', it gives + -- a “flat image” of the region, whereas for 'Shade'' it gives an “extrusion + -- pillar” pointing in the projection's orthogonal complement. + embedShade :: ( Semimanifold x, Semimanifold y + , Object (Affine s) (Interior x), Object (Affine s) (Interior y) + , SemiInner (Needle x), SemiInner (Needle y) ) + => Embedding (Affine s) (Interior x) (Interior y) + -> shade x -> shade y + + +linearProjectShade :: ∀ s x y + . (Num' s, LinearSpace x, LinearSpace y, Scalar x ~ s, Scalar y ~ s) + => (x+>y) -> Shade x -> Shade y +linearProjectShade = case ( linearManifoldWitness :: LinearManifoldWitness x + , linearManifoldWitness :: LinearManifoldWitness y + , dualSpaceWitness :: DualSpaceWitness x + , dualSpaceWitness :: DualSpaceWitness y ) of + ( LinearManifoldWitness BoundarylessWitness + ,LinearManifoldWitness BoundarylessWitness + ,DualSpaceWitness, DualSpaceWitness ) + -> \f (Shade x ex) -> Shade (f $ x) (transformVariance f ex) + + +infixl 5 ✠ +-- | Combine two shades on independent subspaces to a shade with the same +-- properties on the subspaces (see 'factoriseShade') and no covariance. +(✠) :: ( IsShade shade, PseudoAffine x, SimpleSpace (Needle x) + , PseudoAffine y, SimpleSpace (Needle y) + , Scalar (Needle x) ~ Scalar (Needle y) ) + => shade x -> shade y -> shade (x,y) +(✠) = orthoShades + +instance IsShade Shade where + shadeCtr f (Shade c e) = fmap (`Shade`e) $ f c + occlusion = occ pseudoAffineWitness dualSpaceWitness + where occ :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x) + , Scalar (Needle x) ~ s, RealFloat' s ) + => PseudoAffineWitness x -> DualNeedleWitness x -> Shade x -> x -> s + occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness (Shade p₀ δ) + = \p -> case toInterior p >>= (.-~.p₀) of + (Just vd) | mSq <- normSq δinv vd + , mSq == mSq -- avoid NaN + -> exp (negate mSq) + _ -> zeroV + where δinv = dualNorm δ + factoriseShade = fs dualSpaceWitness dualSpaceWitness + where fs :: ∀ x y . ( PseudoAffine x, SimpleSpace (Needle x) + , PseudoAffine y, SimpleSpace (Needle y) + , Scalar (Needle x) ~ Scalar (Needle y) ) + => DualNeedleWitness x -> DualNeedleWitness y + -> Shade (x,y) -> (Shade x, Shade y) + fs DualSpaceWitness DualSpaceWitness (Shade (x₀,y₀) δxy) + = (Shade x₀ δx, Shade y₀ δy) + where (δx,δy) = summandSpaceNorms δxy + orthoShades = fs dualSpaceWitness dualSpaceWitness + where fs :: ∀ x y . ( PseudoAffine x, SimpleSpace (Needle x) + , PseudoAffine y, SimpleSpace (Needle y) + , Scalar (Needle x) ~ Scalar (Needle y) ) + => DualNeedleWitness x -> DualNeedleWitness y + -> Shade x -> Shade y -> Shade (x,y) + fs DualSpaceWitness DualSpaceWitness (Shade x δx) (Shade y δy) + = Shade (x,y) $ sumSubspaceNorms δx δy + coerceShade = cS dualSpaceWitness dualSpaceWitness + where cS :: ∀ x y . (LocallyCoercible x y) + => DualNeedleWitness x -> DualNeedleWitness y -> Shade x -> Shade y + cS DualSpaceWitness DualSpaceWitness + = \(Shade x δxym) -> Shade (internCoerce x) (tN δxym) + where tN = case oppositeLocalCoercion :: CanonicalDiffeomorphism y x of + CanonicalDiffeomorphism -> + transformNorm . arr $ coerceNeedle' ([]::[(y,x)]) + internCoerce = case interiorLocalCoercion ([]::[(x,y)]) of + CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism + linIsoTransformShade = lits linearManifoldWitness linearManifoldWitness + dualSpaceWitness dualSpaceWitness + where lits :: ∀ x y . ( LinearSpace x, LinearSpace y + , Scalar x ~ Scalar y, Num' (Scalar x) ) + => LinearManifoldWitness x -> LinearManifoldWitness y + -> DualSpaceWitness x -> DualSpaceWitness y + -> (x+>y) -> Shade x -> Shade y + lits (LinearManifoldWitness BoundarylessWitness) + (LinearManifoldWitness BoundarylessWitness) + DualSpaceWitness DualSpaceWitness + f (Shade x δx) + = Shade (f $ x) (transformNorm (adjoint $ f) δx) + embedShade = ps' (semimanifoldWitness, semimanifoldWitness) + where ps' :: ∀ s x y . ( Object (Affine s) (Interior x), Object (Affine s) (Interior y) + , SemiInner (Needle x), SemiInner (Needle y) ) + => (SemimanifoldWitness x, SemimanifoldWitness y) + -> Embedding (Affine s) (Interior x) (Interior y) + -> Shade x -> Shade y + ps' (SemimanifoldWitness _, SemimanifoldWitness _) + (Embedding q _) (Shade x e) = Shade y (transformVariance j e) + where y = q $ x + (_,j) = evalAffine q x + projectShade = ps' (semimanifoldWitness, semimanifoldWitness) + where ps' :: ∀ s x y . ( Object (Affine s) (Interior x), Object (Affine s) (Interior y) + , SemiInner (Needle x), SemiInner (Needle y) ) + => (SemimanifoldWitness x, SemimanifoldWitness y) + -> Embedding (Affine s) (Interior x) (Interior y) + -> Shade y -> Shade x + ps' (SemimanifoldWitness _, SemimanifoldWitness _) + (Embedding _ q) (Shade x e) = Shade y (transformVariance j e) + where y = q $ x + (_,j) = evalAffine q x + + +dualShade :: ∀ x . (PseudoAffine x, SimpleSpace (Needle x)) + => Shade x -> Shade' x +dualShade = case dualSpaceWitness :: DualSpaceWitness (Needle x) of + DualSpaceWitness -> \(Shade c e) -> Shade' c $ dualNorm e + +dualShade' :: ∀ x . (PseudoAffine x, SimpleSpace (Needle x)) + => Shade' x -> Shade x +dualShade' = case dualSpaceWitness :: DualSpaceWitness (Needle x) of + DualSpaceWitness -> \(Shade' c e) -> Shade c $ dualNorm' e + +instance ImpliesMetric Shade where + type MetricRequirement Shade x = (Manifold x, SimpleSpace (Needle x)) + inferMetric' (Shade _ e) = e + inferMetric = im dualSpaceWitness + where im :: (Manifold x, SimpleSpace (Needle x)) + => DualNeedleWitness x -> Shade x -> Metric x + im DualSpaceWitness (Shade _ e) = dualNorm e + +instance ImpliesMetric Shade' where + type MetricRequirement Shade' x = (Manifold x, SimpleSpace (Needle x)) + inferMetric (Shade' _ e) = e + inferMetric' (Shade' _ e) = dualNorm e + +shadeExpanse :: Lens' (Shade x) (Metric' x) +shadeExpanse f (Shade c e) = fmap (Shade c) $ f e + +instance IsShade Shade' where + shadeCtr f (Shade' c e) = fmap (`Shade'`e) $ f c + occlusion = occ pseudoAffineWitness + where occ :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x) + , Scalar (Needle x) ~ s, RealFloat' s ) + => PseudoAffineWitness x -> Shade' x -> x -> s + occ (PseudoAffineWitness (SemimanifoldWitness _)) (Shade' p₀ δinv) p + = case toInterior p >>= (.-~.p₀) of + (Just vd) | mSq <- normSq δinv vd + , mSq == mSq -- avoid NaN + -> exp (negate mSq) + _ -> zeroV + factoriseShade (Shade' (x₀,y₀) δxy) = (Shade' x₀ δx, Shade' y₀ δy) + where (δx,δy) = summandSpaceNorms δxy + orthoShades (Shade' x δx) (Shade' y δy) = Shade' (x,y) $ sumSubspaceNorms δx δy + coerceShade = cS + where cS :: ∀ x y . (LocallyCoercible x y) => Shade' x -> Shade' y + cS = \(Shade' x δxym) -> Shade' (internCoerce x) (tN δxym) + where tN = case oppositeLocalCoercion :: CanonicalDiffeomorphism y x of + CanonicalDiffeomorphism -> + transformNorm . arr $ coerceNeedle ([]::[(y,x)]) + internCoerce = case interiorLocalCoercion ([]::[(x,y)]) of + CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism + linIsoTransformShade = lits linearManifoldWitness linearManifoldWitness + dualSpaceWitness dualSpaceWitness + where lits :: ∀ x y . ( SimpleSpace x, SimpleSpace y + , Scalar x ~ Scalar y, RealFloat' (Scalar x) ) + => LinearManifoldWitness x -> LinearManifoldWitness y + -> DualSpaceWitness x -> DualSpaceWitness y + -> (x+>y) -> Shade' x -> Shade' y + lits (LinearManifoldWitness BoundarylessWitness) + (LinearManifoldWitness BoundarylessWitness) + DualSpaceWitness DualSpaceWitness + f (Shade' x δx) + = Shade' (f $ x) (transformNorm (pseudoInverse f) δx) + embedShade = ps (semimanifoldWitness, semimanifoldWitness) + where ps :: ∀ s x y . ( Object (Affine s) (Interior x), Object (Affine s) (Interior y) + , SemiInner (Needle x), SemiInner (Needle y) ) + => (SemimanifoldWitness x, SemimanifoldWitness y) + -> Embedding (Affine s) (Interior x) (Interior y) + -> Shade' x -> Shade' y + ps (SemimanifoldWitness _, SemimanifoldWitness _) + (Embedding q p) (Shade' x e) = Shade' y (transformNorm j e) + where y = q $ x + (_,j) = evalAffine p y + projectShade = ps (semimanifoldWitness, semimanifoldWitness) + where ps :: ∀ s x y . ( Object (Affine s) (Interior x), Object (Affine s) (Interior y) + , SemiInner (Needle x), SemiInner (Needle y) ) + => (SemimanifoldWitness x, SemimanifoldWitness y) + -> Embedding (Affine s) (Interior x) (Interior y) + -> Shade' y -> Shade' x + ps (SemimanifoldWitness _, SemimanifoldWitness _) + (Embedding p q) (Shade' x e) = Shade' y (transformNorm j e) + where y = q $ x + (_,j) = evalAffine p y + + +shadeNarrowness :: Lens' (Shade' x) (Metric x) +shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e + +instance ∀ x . (PseudoAffine x) => Semimanifold (Shade x) where + type Needle (Shade x) = Needle x + fromInterior = id + toInterior = pure + translateP = Tagged (.+~^) + (.+~^) = case semimanifoldWitness :: SemimanifoldWitness x of + SemimanifoldWitness BoundarylessWitness + -> \(Shade c e) v -> Shade (c.+~^v) e + (.-~^) = case semimanifoldWitness :: SemimanifoldWitness x of + SemimanifoldWitness BoundarylessWitness + -> \(Shade c e) v -> Shade (c.-~^v) e + semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of + (SemimanifoldWitness BoundarylessWitness) + -> SemimanifoldWitness BoundarylessWitness + +instance (WithField ℝ PseudoAffine x, Geodesic (Interior x), SimpleSpace (Needle x)) + => Geodesic (Shade x) where + geodesicBetween = gb dualSpaceWitness + where gb :: DualNeedleWitness x -> Shade x -> Shade x -> Maybe (D¹ -> Shade x) + gb DualSpaceWitness (Shade c (Norm e)) (Shade ζ (Norm η)) = pure interp + where interp t@(D¹ q) = Shade (pinterp t) + (Norm . arr . lerp ed ηd $ (q+1)/2) + ed@(LinearMap _) = arr e + ηd@(LinearMap _) = arr η + Just pinterp = geodesicBetween c ζ + +instance (AffineManifold x) => Semimanifold (Shade' x) where + type Needle (Shade' x) = Needle x + fromInterior = id + toInterior = pure + translateP = Tagged (.+~^) + (.+~^) = case boundarylessWitness :: BoundarylessWitness x of + BoundarylessWitness -> \(Shade' c e) v -> Shade' (c.+~^v) e + (.-~^) = case boundarylessWitness :: BoundarylessWitness x of + BoundarylessWitness -> \(Shade' c e) v -> Shade' (c.-~^v) e + semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of + SemimanifoldWitness BoundarylessWitness -> SemimanifoldWitness BoundarylessWitness + +instance ∀ x . (WithField ℝ AffineManifold x, Geodesic x, SimpleSpace (Needle x)) + => Geodesic (Shade' x) where + geodesicBetween (Shade' c e) (Shade' ζ η) = pure interp + where sharedSpan = sharedNormSpanningSystem e η + interp t = Shade' (pinterp t) + (spanNorm [ v ^/ (alerpB 1 (recip qη) t) + | (v,qη) <- sharedSpan ]) + Just pinterp = case geodesicWitness :: GeodesicWitness x of + GeodesicWitness _ -> geodesicBetween c ζ + +fullShade :: WithField ℝ PseudoAffine x => Interior x -> Metric' x -> Shade x +fullShade ctr expa = Shade ctr expa + +fullShade' :: WithField ℝ PseudoAffine x => Interior x -> Metric x -> Shade' x +fullShade' ctr expa = Shade' ctr expa + + +infixl 6 :±, |±| + +-- | Span a 'Shade' from a center point and multiple deviation-vectors. +#if GLASGOW_HASKELL < 800 +pattern (:±) :: () +#else +pattern (:±) :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) +#endif + => (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => Interior x -> [Needle x] -> Shade x +pattern x :± shs <- Shade x (varianceSpanningSystem -> shs) + where x :± shs = fullShade x $ spanVariance shs + +-- | Similar to ':±', but instead of expanding the shade, each vector /restricts/ it. +-- Iff these form a orthogonal basis (in whatever sense applicable), then both +-- methods will be equivalent. +-- +-- Note that '|±|' is only possible, as such, in an inner-product space; in +-- general you need reciprocal vectors ('Needle'') to define a 'Shade''. +(|±|) :: ∀ x . WithField ℝ EuclidSpace x => x -> [Needle x] -> Shade' x +(|±|) = case boundarylessWitness :: BoundarylessWitness x of + BoundarylessWitness -> \x shs -> Shade' x $ spanNorm [v^/(v<.>v) | v<-shs] + + + + + + +-- | Attempt to find a 'Shade' that describes the distribution of given points. +-- At least in an affine space (and thus locally in any manifold), this can be used to +-- estimate the parameters of a normal distribution from which some points were +-- sampled. Note that some points will be “outside” of the shade, +-- as happens for a normal distribution with some statistical likelyhood. +-- (Use 'pointsCovers' if you need to prevent that.) +-- +-- For /nonconnected/ manifolds it will be necessary to yield separate shades +-- for each connected component. And for an empty input list, there is no shade! +-- Hence the result type is a list. +pointsShades :: (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) + => [Interior x] -> [Shade x] +pointsShades = map snd . pointsShades' mempty . map fromInterior + +coverAllAround :: ∀ x s . ( Fractional' s, WithField s PseudoAffine x + , SimpleSpace (Needle x) ) + => Interior x -> [Needle x] -> Shade x +coverAllAround x₀ offs = Shade x₀ + $ guaranteeIn dualSpaceWitness offs + (scaleNorm (1/fromIntegral (length offs)) $ spanVariance offs) + where guaranteeIn :: DualNeedleWitness x -> [Needle x] -> Metric' x -> Metric' x + guaranteeIn w@DualSpaceWitness offs ex + = case offs >>= \v -> guard ((ex'|$|v) > 1) >> [(v, spanVariance [v])] of + [] -> ex + outs -> guaranteeIn w (fst<$>outs) + ( densifyNorm $ + ex <> scaleNorm + (sqrt . recip . fromIntegral + $ 2 * length outs) + (mconcat $ snd<$>outs) + ) + where ex' = dualNorm ex + +-- | Like 'pointsShades', but ensure that all points are actually in +-- the shade, i.e. if @['Shade' x₀ ex]@ is the result then +-- @'metric' (recipMetric ex) (p-x₀) ≤ 1@ for all @p@ in the list. +pointsCovers :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) + => [Interior x] -> [Shade x] +pointsCovers = case pseudoAffineWitness :: PseudoAffineWitness x of + (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) -> + \ps -> map (\(ps', Shade x₀ _) + -> coverAllAround x₀ [v | p<-ps' + , let Just v + = p.-~.fromInterior x₀]) + (pointsShades' mempty (fromInterior<$>ps) :: [([x], Shade x)]) + +pointsShade's :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) + => [Interior x] -> [Shade' x] +pointsShade's = case dualSpaceWitness :: DualNeedleWitness x of + DualSpaceWitness -> map (\(Shade c e :: Shade x) -> Shade' c $ dualNorm e) . pointsShades + +pointsCover's :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) + => [Interior x] -> [Shade' x] +pointsCover's = case dualSpaceWitness :: DualNeedleWitness x of + DualSpaceWitness -> map (\(Shade c e :: Shade x) -> Shade' c $ dualNorm e) . pointsCovers + +pseudoECM :: ∀ x p . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x), Hask.Functor p) + => p x -> NonEmpty x -> (x, ([x],[x])) +pseudoECM = case semimanifoldWitness :: SemimanifoldWitness x of + SemimanifoldWitness _ -> + \_ (p₀ NE.:| psr) -> foldl' ( \(acc, (rb,nr)) (i,p) + -> case (p.-~.acc, toInterior acc) of + (Just δ, Just acci) + -> (acci .+~^ δ^/i, (p:rb, nr)) + _ -> (acc, (rb, p:nr)) ) + (p₀, mempty) + ( zip [1..] $ p₀:psr ) + +pointsShades' :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) + => Metric' x -> [x] -> [([x], Shade x)] +pointsShades' _ [] = [] +pointsShades' minExt ps = case (expa, toInterior ctr) of + (Just e, Just c) + -> (ps, fullShade c e) : pointsShades' minExt unreachable + _ -> pointsShades' minExt inc'd + ++ pointsShades' minExt unreachable + where (ctr,(inc'd,unreachable)) = pseudoECM ([]::[x]) $ NE.fromList ps + expa = ( (<>minExt) . spanVariance . map (^/ fromIntegral (length ps)) ) + <$> mapM (.-~.ctr) ps + + +-- | Attempt to reduce the number of shades to fewer (ideally, a single one). +-- In the simplest cases these should guaranteed cover the same area; +-- for non-flat manifolds it only works in a heuristic sense. +shadesMerge :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => ℝ -- ^ How near (inverse normalised distance, relative to shade expanse) + -- two shades must be to be merged. If this is zero, any shades + -- in the same connected region of a manifold are merged. + -> [Shade x] -- ^ A list of /n/ shades. + -> [Shade x] -- ^ /m/ ≤ /n/ shades which cover at least the same area. +shadesMerge fuzz (sh₁@(Shade c₁ e₁) : shs) + = case extractJust (tryMerge pseudoAffineWitness dualSpaceWitness) + shs of + (Just mg₁, shs') -> shadesMerge fuzz + $ shs'++[mg₁] -- Append to end to prevent undue weighting + -- of first shade and its mergers. + (_, shs') -> sh₁ : shadesMerge fuzz shs' + where tryMerge :: PseudoAffineWitness x -> DualNeedleWitness x + -> Shade x -> Maybe (Shade x) + tryMerge (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) DualSpaceWitness + (Shade c₂ e₂) + | Just v <- c₁.-~.c₂ + , [e₁',e₂'] <- dualNorm<$>[e₁, e₂] + , b₁ <- e₂'|$|v + , b₂ <- e₁'|$|v + , fuzz*b₁*b₂ <= b₁ + b₂ + = Just $ let cc = c₂ .+~^ v ^/ 2 + Just cv₁ = c₁.-~.cc + Just cv₂ = c₂.-~.cc + in Shade cc $ e₁ <> e₂ <> spanVariance [cv₁, cv₂] + | otherwise = Nothing +shadesMerge _ shs = shs + +-- | Weakened version of 'intersectShade's'. What this function calculates is +-- rather the /weighted mean/ of ellipsoid regions. If you interpret the +-- shades as uncertain physical measurements with normal distribution, +-- it gives the maximum-likelyhood result for multiple measurements of the +-- same quantity. +mixShade's :: ∀ y . (WithField ℝ Manifold y, SimpleSpace (Needle y)) + => NonEmpty (Shade' y) -> Maybe (Shade' y) +mixShade's = ms pseudoAffineWitness dualSpaceWitness + where ms :: PseudoAffineWitness y -> DualNeedleWitness y + -> NonEmpty (Shade' y) -> Maybe (Shade' y) + ms (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness + (Shade' c₀ (Norm e₁):|shs) = sequenceA ciso >> pure mixed + where ciso = [ci.-~.c₀ | Shade' ci shi <- shs] + cis = [v | Just v <- ciso] + σe = arr . sumV $ e₁ : (applyNorm . _shade'Narrowness<$>shs) + cc = σe \$ sumV [ei $ ci | ci <- cis + | Shade' _ (Norm ei) <- shs] + mixed = Shade' (c₀+^cc) $ densifyNorm ( mconcat + [ Norm $ ei ^/ (1+(normSq ni $ ci^-^cc)) + | ni@(Norm ei) <- Norm e₁ : (_shade'Narrowness<$>shs) + | ci <- zeroV : cis + ] ) + Tagged (+^) = translateP :: Tagged y (Interior y->Needle y->Interior y) + -- cc should minimise the quadratic form + -- β(cc) = ∑ᵢ ⟨cc−cᵢ|eᵢ|cc−cᵢ⟩ + -- = ⟨cc|e₁|cc⟩ + ∑ᵢ₌₁… ⟨cc−c₂|e₂|cc−c₂⟩ + -- = ⟨cc|e₁|cc⟩ + ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ − 2⋅⟨cᵢ|eᵢ|cc⟩ + ⟨cᵢ|eᵢ|cᵢ⟩ ) + -- It is thus + -- β(cc + δ⋅v) − β cc + -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + -- + ∑ᵢ₌₁…( ⟨cc + δ⋅v|eᵢ|cc + δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|cc + δ⋅v⟩ + ⟨cᵢ|eᵢ|cᵢ⟩ ) + -- − ⟨cc|e₁|cc⟩ + -- − ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ + 2⋅⟨cᵢ|eᵢ|cc⟩ − ⟨cᵢ|eᵢ|cᵢ⟩ ) + -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + -- + ∑ᵢ₌₁…( ⟨cc + δ⋅v|eᵢ|cc + δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|δ⋅v⟩ ) + -- − ⟨cc|e₁|cc⟩ + -- − ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ ) + -- = 2⋅⟨δ⋅v|e₁|cc⟩ + ⟨δ⋅v|e₁|δ⋅v⟩ + -- + ∑ᵢ₌₁…( 2⋅⟨δ⋅v|eᵢ|cc⟩ + ⟨δ⋅v|eᵢ|δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|δ⋅v⟩ ) + -- = 2⋅⟨δ⋅v|∑ᵢeᵢ|cc⟩ − 2⋅∑ᵢ₌₁… ⟨cᵢ|eᵢ|δ⋅v⟩ + 𝓞(δ²) + -- This should vanish for all v, which is fulfilled by + -- (∑ᵢeᵢ)|cc⟩ = ∑ᵢ₌₁… eᵢ|cᵢ⟩. + +-- | Evaluate the shade as a quadratic form; essentially +-- @ +-- minusLogOcclusion sh x = x <.>^ (sh^.shadeExpanse $ x - sh^.shadeCtr) +-- @ +-- where 'shadeExpanse' gives a metric (matrix) that characterises the +-- width of the shade. +minusLogOcclusion' :: ∀ x s . ( PseudoAffine x, LinearSpace (Needle x) + , s ~ (Scalar (Needle x)), RealFloat' s ) + => Shade' x -> x -> s +minusLogOcclusion' (Shade' p₀ δinv) + = occ (pseudoAffineWitness :: PseudoAffineWitness x) + (dualSpaceWitness :: DualNeedleWitness x) + where occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness + p = case toInterior p >>= (.-~.p₀) of + (Just vd) | mSq <- normSq δinv vd + , mSq == mSq -- avoid NaN + -> mSq + _ -> 1/0 +minusLogOcclusion :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x) + , s ~ (Scalar (Needle x)), RealFloat' s ) + => Shade x -> x -> s +minusLogOcclusion (Shade p₀ δ) + = occ (pseudoAffineWitness :: PseudoAffineWitness x) + (dualSpaceWitness :: DualNeedleWitness x) + where occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness + = \p -> case toInterior p >>= (.-~.p₀) of + (Just vd) | mSq <- normSq δinv vd + , mSq == mSq -- avoid NaN + -> mSq + _ -> 1/0 + where δinv = dualNorm δ + + + + +{-# WARNING rangeOnGeodesic "This function never worked properly. Use 'rangeWithinVertices'." #-} +rangeOnGeodesic :: ∀ i m . + ( WithField ℝ PseudoAffine m, Geodesic m, SimpleSpace (Needle m) + , WithField ℝ IntervalLike i, SimpleSpace (Needle i) ) + => m -> m -> Maybe (Shade i -> Shade m) +rangeOnGeodesic = case ( semimanifoldWitness :: SemimanifoldWitness i + , dualSpaceWitness :: DualNeedleWitness i + , dualSpaceWitness :: DualNeedleWitness m ) of + (SemimanifoldWitness _, DualSpaceWitness, DualSpaceWitness) -> + \p₀ p₁ -> geodesicBetween p₀ p₁ >>= + \interp -> case pointsShades =<< + [ mapMaybe (toInterior . interp . D¹) [-(1-ε), 1-ε] + | ε <- [0.0001, 0.001, 0.01, 0.1] ] of + defaultSh:_ -> Just $ + \(Shade t₀ et) -> case pointsShades + . mapMaybe (toInterior + . interp . (toClosedInterval :: i -> D¹)) + $ fromInterior <$> t₀ : [ t₀+^v + | v<-normSpanningSystem et ] of + [sh] -> sh + _ -> defaultSh + _ -> Nothing + where Tagged (+^) = translateP :: Tagged i (Interior i->Needle i->Interior i) + + +rangeWithinVertices :: ∀ s i m t + . ( RealFrac' s + , WithField s PseudoAffine i, WithField s PseudoAffine m + , Geodesic i, Geodesic m + , SimpleSpace (Needle i), SimpleSpace (Needle m) + , AffineManifold (Interior i), AffineManifold (Interior m) + , Object (Affine s) (Interior i), Object (Affine s) (Interior m) + , Hask.Traversable t ) + => (Interior i,Interior m) -> t (i,m) -> Maybe (Shade i -> Shade m) +rangeWithinVertices + = case ( semimanifoldWitness :: SemimanifoldWitness i + , semimanifoldWitness :: SemimanifoldWitness m ) of + (SemimanifoldWitness BoundarylessWitness, SemimanifoldWitness BoundarylessWitness) + -> \(cii,cmi) verts -> + let ci = fromInterior cii + cm = fromInterior cmi + in do + vs <- sequenceA [ fzip ( middleBetween pi ci >>= (.-~.ci) + , middleBetween pm cm >>= (.-~.cm) ) + | (pi, pm) <- Hask.toList verts ] + affinSys <- (correspondingDirections (cii,cmi) vs + :: Maybe (Embedding (Affine (Scalar (Needle i))) + (Interior i) (Interior m))) + return $ embedShade affinSys + + + + + +data DebugView x where + DebugView :: ( Show x, Show (Needle x+>Needle' x), LinearShowable (Needle x) + , Needle' x ~ Needle x ) => DebugView x + +-- | Class of manifolds which can use 'Shade'' as a basic set type. +-- This is easily possible for vector spaces with the default implementations. +class (WithField ℝ PseudoAffine y, SimpleSpace (Needle y)) => Refinable y where + debugView :: Maybe (DebugView y) + default debugView :: ( Show y, Show (Needle y+>Needle' y) + , Needle' y~Needle y, LinearShowable (Needle y) ) + => Maybe (DebugView y) + debugView = Just DebugView + + -- | @a `subShade'` b ≡ True@ means @a@ is fully contained in @b@, i.e. from + -- @'minusLogOcclusion'' a p < 1@ follows also @minusLogOcclusion' b p < 1@. + subShade' :: Shade' y -> Shade' y -> Bool + subShade' (Shade' ac ae) (Shade' tc te) + = case pseudoAffineWitness :: PseudoAffineWitness y of + PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) + | Just v <- tc.-~.ac + , v² <- normSq te v + , v² <= 1 + -> all (\(y',μ) -> case μ of + Nothing -> True -- 'te' has infinite extension in this direction + Just ξ + | ξ<1 -> False -- 'ae' would be vaster than 'te' in this direction + | ω <- abs $ y'<.>^v + -> (ω + 1/ξ)^2 <= 1 - v² + ω^2 + -- See @images/constructions/subellipse-check-heuristic.svg@ + ) $ sharedSeminormSpanningSystem te ae + _ -> False + + -- | Intersection between two shades. + refineShade' :: Shade' y -> Shade' y -> Maybe (Shade' y) + refineShade' (Shade' c₀ (Norm e₁)) (Shade' c₀₂ (Norm e₂)) + = case ( dualSpaceWitness :: DualNeedleWitness y + , pseudoAffineWitness :: PseudoAffineWitness y ) of + (DualSpaceWitness, PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) + -> do + c₂ <- c₀₂.-~.c₀ + let σe = arr $ e₁^+^e₂ + e₁c₂ = e₁ $ c₂ + e₂c₂ = e₂ $ c₂ + cc = σe \$ e₂c₂ + cc₂ = cc ^-^ c₂ + e₁cc = e₁ $ cc + e₂cc = e₂ $ cc + α = 2 + e₂c₂<.>^cc₂ + guard (α > 0) + let ee = σe ^/ α + c₂e₁c₂ = e₁c₂<.>^c₂ + c₂e₂c₂ = e₂c₂<.>^c₂ + c₂eec₂ = (c₂e₁c₂ + c₂e₂c₂) / α + return $ case middle . sort + $ quadraticEqnSol c₂e₁c₂ + (2 * (e₁cc<.>^c₂)) + (e₁cc<.>^cc - 1) + ++quadraticEqnSol c₂e₂c₂ + (2 * (e₂cc<.>^c₂ - c₂e₂c₂)) + (e₂cc<.>^cc - 2 * (e₂c₂<.>^cc) + c₂e₂c₂ - 1) of + [γ₁,γ₂] | abs (γ₁+γ₂) < 2 -> let + cc' = cc ^+^ ((γ₁+γ₂)/2)*^c₂ + rγ = abs (γ₁ - γ₂) / 2 + η = if rγ * c₂eec₂ /= 0 && 1 - rγ^2 * c₂eec₂ > 0 + then sqrt (1 - rγ^2 * c₂eec₂) / (rγ * c₂eec₂) + else 0 + in Shade' (c₀.+~^cc') + (Norm (arr ee) <> spanNorm [ee $ c₂^*η]) + _ -> Shade' (c₀.+~^cc) (Norm $ arr ee) + where quadraticEqnSol a b c + | a == 0, b /= 0 = [-c/b] + | a /= 0 && disc == 0 = [- b / (2*a)] + | a /= 0 && disc > 0 = [ (σ * sqrt disc - b) / (2*a) + | σ <- [-1, 1] ] + | otherwise = [] + where disc = b^2 - 4*a*c + middle (_:x:y:_) = [x,y] + middle l = l + -- ⟨x−c₁|e₁|x−c₁⟩ < 1 ∧ ⟨x−c₂|e₂|x−c₂⟩ < 1 + -- We search (cc,ee) such that this implies + -- ⟨x−cc|ee|x−cc⟩ < 1. + -- Let WLOG c₁ = 0, so + -- ⟨x|e₁|x⟩ < 1. + -- cc should minimise the quadratic form + -- β(cc) = ⟨cc−c₁|e₁|cc−c₁⟩ + ⟨cc−c₂|e₂|cc−c₂⟩ + -- = ⟨cc|e₁|cc⟩ + ⟨cc−c₂|e₂|cc−c₂⟩ + -- = ⟨cc|e₁|cc⟩ + ⟨cc|e₂|cc⟩ − 2⋅⟨c₂|e₂|cc⟩ + ⟨c₂|e₂|c₂⟩ + -- It is thus + -- β(cc + δ⋅v) − β cc + -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + ⟨cc + δ⋅v|e₂|cc + δ⋅v⟩ − 2⋅⟨c₂|e₂|cc + δ⋅v⟩ + ⟨c₂|e₂|c₂⟩ + -- − ⟨cc|e₁|cc⟩ − ⟨cc|e₂|cc⟩ + 2⋅⟨c₂|e₂|cc⟩ − ⟨c₂|e₂|c₂⟩ + -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + ⟨cc + δ⋅v|e₂|cc + δ⋅v⟩ − 2⋅⟨c₂|e₂|δ⋅v⟩ + -- − ⟨cc|e₁|cc⟩ − ⟨cc|e₂|cc⟩ + -- = 2⋅⟨δ⋅v|e₁|cc⟩ + ⟨δ⋅v|e₁|δ⋅v⟩ + 2⋅⟨δ⋅v|e₂|cc⟩ + ⟨δ⋅v|e₂|δ⋅v⟩ − 2⋅⟨c₂|e₂|δ⋅v⟩ + -- = 2⋅δ⋅⟨v|e₁+e₂|cc⟩ − 2⋅δ⋅⟨v|e₂|c₂⟩ + 𝓞(δ²) + -- This should vanish for all v, which is fulfilled by + -- (e₁+e₂)|cc⟩ = e₂|c₂⟩. + -- + -- If we now choose + -- ee = (e₁+e₂) / α + -- then + -- ⟨x−cc|ee|x−cc⟩ ⋅ α + -- = ⟨x−cc|ee|x⟩ ⋅ α − ⟨x−cc|ee|cc⟩ ⋅ α + -- = ⟨x|ee|x−cc⟩ ⋅ α − ⟨x−cc|e₂|c₂⟩ + -- = ⟨x|ee|x⟩ ⋅ α − ⟨x|ee|cc⟩ ⋅ α − ⟨x−cc|e₂|c₂⟩ + -- = ⟨x|e₁+e₂|x⟩ − ⟨x|e₂|c₂⟩ − ⟨x−cc|e₂|c₂⟩ + -- = ⟨x|e₁|x⟩ + ⟨x|e₂|x⟩ − ⟨x|e₂|c₂⟩ − ⟨x−cc|e₂|c₂⟩ + -- < 1 + ⟨x|e₂|x−c₂⟩ − ⟨x−cc|e₂|c₂⟩ + -- = 1 + ⟨x−c₂|e₂|x−c₂⟩ + ⟨c₂|e₂|x−c₂⟩ − ⟨x−cc|e₂|c₂⟩ + -- < 2 + ⟨x−c₂−x+cc|e₂|c₂⟩ + -- = 2 + ⟨cc−c₂|e₂|c₂⟩ + -- Really we want + -- ⟨x−cc|ee|x−cc⟩ ⋅ α < α + -- So choose α = 2 + ⟨cc−c₂|e₂|c₂⟩. + -- + -- The ellipsoid "cc±√ee" captures perfectly the intersection + -- of the boundary of the shades, but it tends to significantly + -- overshoot the interior intersection in perpendicular direction, + -- i.e. in direction of c₂−c₁. E.g. + -- https://github.com/leftaroundabout/manifolds/blob/bc0460b9/manifolds/images/examples/ShadeCombinations/EllipseIntersections.png + -- 1. Really, the relevant points are those where either of the + -- intersector badnesses becomes 1. The intersection shade should + -- be centered between those points. We perform according corrections, + -- but only in c₂ direction, so this can be handled efficiently + -- as a 1D quadratic equation. + -- Consider + -- dⱼ c := ⟨c−cⱼ|eⱼ|c−cⱼ⟩ =! 1 + -- dⱼ (cc + γ⋅c₂) + -- = ⟨cc+γ⋅c₂−cⱼ|eⱼ|cc+γ⋅c₂−cⱼ⟩ + -- = ⟨cc−cⱼ|eⱼ|cc−cⱼ⟩ + 2⋅γ⋅⟨c₂|eⱼ|cc−cⱼ⟩ + γ²⋅⟨c₂|eⱼ|c₂⟩ + -- =! 1 + -- So + -- γⱼ = (- b ± √(b²−4⋅a⋅c)) / 2⋅a + -- where a = ⟨c₂|eⱼ|c₂⟩ + -- b = 2 ⋅ (⟨c₂|eⱼ|cc⟩ − ⟨c₂|eⱼ|cⱼ⟩) + -- c = ⟨cc|eⱼ|cc⟩ − 2⋅⟨cc|eⱼ|cⱼ⟩ + ⟨cⱼ|eⱼ|cⱼ⟩ − 1 + -- The ± sign should be chosen to get the smaller |γ| (otherwise + -- we end up on the wrong side of the shade), i.e. + -- γⱼ = (sgn bⱼ ⋅ √(bⱼ²−4⋅aⱼ⋅cⱼ) − bⱼ) / 2⋅aⱼ + -- 2. Trim the result in that direction to the actual + -- thickness of the lens-shaped intersection: we want + -- ⟨rγ⋅c₂|ee'|rγ⋅c₂⟩ = 1 + -- for a squeezed version of ee, + -- ee' = ee + ee|η⋅c₂⟩⟨η⋅c₂|ee + -- ee' = ee + η² ⋅ ee|c₂⟩⟨c₂|ee + -- ⟨rγ⋅c₂|ee'|rγ⋅c₂⟩ + -- = rγ² ⋅ (⟨c₂|ee|c₂⟩ + η² ⋅ ⟨c₂|ee|c₂⟩²) + -- = rγ² ⋅ ⟨c₂|ee|c₂⟩ + η² ⋅ rγ² ⋅ ⟨c₂|ee|c₂⟩² + -- η² = (1 − rγ²⋅⟨c₂|ee|c₂⟩) / (rγ² ⋅ ⟨c₂|ee|c₂⟩²) + -- η = √(1 − rγ²⋅⟨c₂|ee|c₂⟩) / (rγ ⋅ ⟨c₂|ee|c₂⟩) + -- With ⟨c₂|ee|c₂⟩ = (⟨c₂|e₁|c₂⟩ + ⟨c₂|e₂|c₂⟩)/α. + + + -- | If @p@ is in @a@ (red) and @δ@ is in @b@ (green), + -- then @p.+~^δ@ is in @convolveShade' a b@ (blue). + -- +-- Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/ShadeCombinations.ipynb#shadeConvolutions +-- +-- <> + convolveMetric :: Hask.Functor p => p y -> Metric y -> Metric y -> Metric y + convolveMetric _ ey eδ = case wellDefinedNorm result of + Just r -> r + Nothing -> case debugView :: Maybe (DebugView y) of + Just DebugView -> error $ "Can not convolve norms " + ++show (arr (applyNorm ey) :: Needle y+>Needle' y) + ++" and "++show (arr (applyNorm eδ) :: Needle y+>Needle' y) + where eδsp = sharedSeminormSpanningSystem ey eδ + result = spanNorm [ f ^* ζ crl | (f,crl) <- eδsp ] + ζ = case filter (>0) . catMaybes $ snd<$>eδsp of + [] -> const 0 + nzrelap + -> let cre₁ = 1/minimum nzrelap + cre₂ = maximum nzrelap + edgeFactor = sqrt ( (1 + cre₁)^2 + (1 + cre₂)^2 ) + / (sqrt (1 + cre₁^2) + sqrt (1 + cre₂^2)) + in \case + Nothing -> 0 + Just 0 -> 0 + Just sq -> edgeFactor / (recip sq + 1) + + convolveShade' :: Shade' y -> Shade' (Needle y) -> Shade' y + convolveShade' = defaultConvolveShade' + +defaultConvolveShade' :: ∀ y . Refinable y => Shade' y -> Shade' (Needle y) -> Shade' y +defaultConvolveShade' = case (pseudoAffineWitness :: PseudoAffineWitness y) of + PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) + -> \(Shade' y₀ ey) (Shade' δ₀ eδ) -> Shade' (y₀.+~^δ₀) + $ convolveMetric ([]::[y]) ey eδ + +instance Refinable ℝ where + refineShade' (Shade' cl el) (Shade' cr er) + = case (normSq el 1, normSq er 1) of + (0, _) -> return $ Shade' cr er + (_, 0) -> return $ Shade' cl el + (ql,qr) | ql>0, qr>0 + -> let [rl,rr] = sqrt . recip <$> [ql,qr] + b = maximum $ zipWith (-) [cl,cr] [rl,rr] + t = minimum $ zipWith (+) [cl,cr] [rl,rr] + in guard (b> + let cm = (b+t)/2 + rm = (t-b)/2 + in return $ Shade' cm (spanNorm [recip rm]) +-- convolveShade' (Shade' y₀ ey) (Shade' δ₀ eδ) +-- = case (metricSq ey 1, metricSq eδ 1) of +-- (wy,wδ) | wy>0, wδ>0 +-- -> Shade' (y₀.+~^δ₀) +-- ( projector . recip +-- $ recip (sqrt wy) + recip (sqrt wδ) ) +-- (_ , _) -> Shade' y₀ zeroV + +instance ∀ a b . ( Refinable a, Refinable b + , Scalar (DualVector (DualVector (Needle b))) + ~ Scalar (DualVector (DualVector (Needle a))) ) + => Refinable (a,b) where + debugView = case ( debugView :: Maybe (DebugView a) + , debugView :: Maybe (DebugView b) + , dualSpaceWitness :: DualSpaceWitness (Needle a) + , dualSpaceWitness :: DualSpaceWitness (Needle b) ) of + (Just DebugView, Just DebugView, DualSpaceWitness, DualSpaceWitness) + -> Just DebugView + +instance Refinable ℝ⁰ +instance Refinable ℝ¹ +instance Refinable ℝ² +instance Refinable ℝ³ +instance Refinable ℝ⁴ + +instance ( SimpleSpace a, SimpleSpace b + , Refinable a, Refinable b + , Scalar a ~ ℝ, Scalar b ~ ℝ + , Scalar (DualVector a) ~ ℝ, Scalar (DualVector b) ~ ℝ + , Scalar (DualVector (DualVector a)) ~ ℝ, Scalar (DualVector (DualVector b)) ~ ℝ ) + => Refinable (LinearMap ℝ a b) where + debugView = Nothing + +intersectShade's :: ∀ y . Refinable y => NonEmpty (Shade' y) -> Maybe (Shade' y) +intersectShade's (sh:|shs) = Hask.foldrM refineShade' sh shs + + +estimateLocalJacobian :: ∀ x y . ( WithField ℝ Manifold x, Refinable y + , SimpleSpace (Needle x), SimpleSpace (Needle y) ) + => Metric x -> [(Local x, Shade' y)] + -> Maybe (Shade' (LocalLinear x y)) +estimateLocalJacobian = elj ( pseudoAffineWitness :: PseudoAffineWitness x + , pseudoAffineWitness :: PseudoAffineWitness y ) + where elj ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) + , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) ) + mex [(Local x₁, Shade' y₁ ey₁),(Local x₀, Shade' y₀ ey₀)] + = return $ Shade' (dx-+|>δy) + (Norm . LinearFunction $ \δj -> δx ⊗ (σey<$|δj $ δx)) + where Just δx = x₁.-~.x₀ + δx' = (mex<$|δx) + dx = δx'^/(δx'<.>^δx) + Just δy = y₁.-~.y₀ + σey = convolveMetric ([]::[y]) ey₀ ey₁ + elj _ mex (po:ps) + | DualSpaceWitness <- dualSpaceWitness :: DualNeedleWitness y + , length ps > 1 + = mixShade's =<< (:|) <$> estimateLocalJacobian mex ps + <*> sequenceA [estimateLocalJacobian mex [po,pi] | pi<-ps] + elj _ _ _ = return $ Shade' zeroV mempty + + + +data QuadraticModel x y = QuadraticModel { + _quadraticModelOffset :: Interior y + , _quadraticModel :: Shade (Needle y, (Needle x+>Needle y, Needle x⊗〃+>Needle y)) + } + +quadratic_linearRegression :: ∀ s x y . + ( WithField s PseudoAffine x + , WithField s PseudoAffine y, Geodesic y + , SimpleSpace (Needle x), SimpleSpace (Needle y) ) + => NE.NonEmpty (Needle x, Shade' y) -> QuadraticModel x y +quadratic_linearRegression = qlr + ( dualSpaceWitness, pseudoAffineWitness + , linearManifoldWitness, dualSpaceWitness + , geodesicWitness ) + where qlr :: ( DualSpaceWitness (Needle x) + , PseudoAffineWitness y, LinearManifoldWitness (Needle y) + , DualSpaceWitness (Needle y) + , GeodesicWitness y ) + -> NE.NonEmpty (Needle x, Shade' y) -> QuadraticModel x y + qlr ( DualSpaceWitness + , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) + , LinearManifoldWitness BoundarylessWitness, DualSpaceWitness + , GeodesicWitness _ ) ps + = QuadraticModel cmy + $ coverAllAround mBest (convexPolytopeRepresentatives dm) + where Just cmy = pointsBarycenter $ _shade'Ctr.snd<$>ps + Just vsxy = Hask.mapM (\(x, Shade' y ey) -> (x,).(,ey)<$>y.-~.cmy) ps + (mBest :: ( Needle y, (Needle x+>Needle y + , SymmetricTensor s (Needle x)+>(Needle y)) + ) + , dm) + = linearRegressionWVar + (\δx -> lfun $ \(c,(b,a)) -> (a $ squareV δx) + ^+^ (b $ δx) ^+^ c ) + (NE.toList vsxy) + +estimateLocalHessian :: ∀ x y . ( WithField ℝ Manifold x, Refinable y, Geodesic y + , FlatSpace (Needle x), FlatSpace (Needle y) ) + => NonEmpty (Local x, Shade' y) -> QuadraticModel x y +estimateLocalHessian pts = elj ( pseudoAffineWitness :: PseudoAffineWitness x + , pseudoAffineWitness :: PseudoAffineWitness y ) + where elj ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) + , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) ) + = theModel + where localPts :: NonEmpty (Needle x, Shade' y) + localPts = pts >>= \(Local x, Shade' y ey) + -> NE.fromList [ (x, Shade' (y.+~^σ*^δy) ey) + | δy <- normSpanningSystem' ey + , σ <- [-1,1] ] + theModel = quadratic_linearRegression localPts + + + +propagateDEqnSolution_loc :: ∀ x y ð . ( WithField ℝ Manifold x + , Refinable y, Geodesic (Interior y) + , WithField ℝ AffineManifold ð, Geodesic ð + , SimpleSpace (Needle x), SimpleSpace (Needle ð) ) + => DifferentialEqn x ð y + -> LocalDataPropPlan x (Shade' y, Shade' ð) (Shade' y) + -> Maybe (Shade' y) +propagateDEqnSolution_loc f propPlan + = pdesl (dualSpaceWitness :: DualNeedleWitness x) + (dualSpaceWitness :: DualNeedleWitness y) + (boundarylessWitness :: BoundarylessWitness x) + (pseudoAffineWitness :: PseudoAffineWitness y) + where pdesl DualSpaceWitness DualSpaceWitness BoundarylessWitness + (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) + | Nothing <- jacobian = Nothing + | otherwise = pure result + where jacobian = (f shxy ^. predictDerivatives $ shð) + >>= \j -> mixShade's $ j:|[aprioriDirDrv] + Just (Shade' j₀ jExpa) = jacobian + jacobianSh :: Shade (LocalLinear x y) + Just jacobianSh = dualShade' <$> jacobian + aprioriDirDrv :: Shade' (LocalLinear x y) + Just aprioriDirDrv = estimateLocalJacobian expax + [ (Local zeroV :: Local x, propPlan^.sourceData._1) + , (Local δx, propPlan^.targetAPrioriData._1) ] + mx = propPlan^.sourcePosition .+~^ propPlan^.targetPosOffset ^/ 2 :: x + Just shð = middleBetween (propPlan^.sourceData._2) + (propPlan^.targetAPrioriData._2) + shxy = coverAllAround (mx, mυ) + [ (δx ^-^ propPlan^.targetPosOffset ^/ 2, pυ ^+^ v) + | (δx,neυ) <- (zeroV, propPlan^.sourceData._1) + : (second id + <$> propPlan^.relatedData) + , let Just pυ = neυ^.shadeCtr .-~. mυ + , v <- normSpanningSystem' (neυ^.shadeNarrowness) + ] + where Just mυ = middleBetween (propPlan^.sourceData._1.shadeCtr) + (propPlan^.targetAPrioriData._1.shadeCtr) + (Shade _ expax' :: Shade x) + = coverAllAround (propPlan^.sourcePosition) + [δx | (δx,_) <- propPlan^.relatedData] + expax = dualNorm expax' + result :: Shade' y + Just result = wellDefinedShade' $ convolveShade' + (case wellDefinedShade' $ propPlan^.sourceData._1 of {Just s->s}) + (case wellDefinedShade' . dualShade + . linearProjectShade (lfun ($ δx)) + $ jacobianSh + of {Just s->s}) + where δyb = j₀ $ δx + δx = propPlan^.targetPosOffset + + + + + + + + +-- | Essentially the same as @(x,y)@, but not considered as a product topology. +-- The 'Semimanifold' etc. instances just copy the topology of @x@, ignoring @y@. +data x`WithAny`y + = WithAny { _untopological :: y + , _topological :: !x } + deriving (Hask.Functor, Show, Generic) + +instance (NFData x, NFData y) => NFData (WithAny x y) + +instance ∀ x y . (Semimanifold x) => Semimanifold (x`WithAny`y) where + type Needle (WithAny x y) = Needle x + type Interior (WithAny x y) = Interior x `WithAny` y + WithAny y x .+~^ δx = WithAny y $ x.+~^δx + fromInterior (WithAny y x) = WithAny y $ fromInterior x + toInterior (WithAny y x) = fmap (WithAny y) $ toInterior x + translateP = tpWD + where tpWD :: ∀ x y . Semimanifold x => Tagged (WithAny x y) + (Interior x`WithAny`y -> Needle x -> Interior x`WithAny`y) + tpWD = Tagged `id` \(WithAny y x) δx -> WithAny y $ tpx x δx + where Tagged tpx = translateP :: Tagged x (Interior x -> Needle x -> Interior x) + semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of + SemimanifoldWitness BoundarylessWitness -> SemimanifoldWitness BoundarylessWitness + +instance (PseudoAffine x) => PseudoAffine (x`WithAny`y) where + WithAny _ x .-~. WithAny _ ξ = x.-~.ξ + pseudoAffineWitness = case pseudoAffineWitness :: PseudoAffineWitness x of + PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) + -> PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) + +instance (AffineSpace x) => AffineSpace (x`WithAny`y) where + type Diff (WithAny x y) = Diff x + WithAny _ x .-. WithAny _ ξ = x.-.ξ + WithAny y x .+^ δx = WithAny y $ x.+^δx + +instance (VectorSpace x, Monoid y) => VectorSpace (x`WithAny`y) where + type Scalar (WithAny x y) = Scalar x + μ *^ WithAny y x = WithAny y $ μ*^x + +instance (AdditiveGroup x, Monoid y) => AdditiveGroup (x`WithAny`y) where + zeroV = WithAny mempty zeroV + negateV (WithAny y x) = WithAny y $ negateV x + WithAny y x ^+^ WithAny υ ξ = WithAny (mappend y υ) (x^+^ξ) + +instance (AdditiveGroup x) => Hask.Applicative (WithAny x) where + pure x = WithAny x zeroV + WithAny f x <*> WithAny t ξ = WithAny (f t) (x^+^ξ) + +instance (AdditiveGroup x) => Hask.Monad (WithAny x) where + return x = WithAny x zeroV + WithAny y x >>= f = WithAny r $ x^+^q + where WithAny r q = f y + +shadeWithAny :: y -> Shade x -> Shade (x`WithAny`y) +shadeWithAny y (Shade x xe) = Shade (WithAny y x) xe + +shadeWithoutAnything :: Shade (x`WithAny`y) -> Shade x +shadeWithoutAnything (Shade (WithAny _ b) e) = Shade b e + + + + + + +extractJust :: (a->Maybe b) -> [a] -> (Maybe b, [a]) +extractJust f [] = (Nothing,[]) +extractJust f (x:xs) | Just r <- f x = (Just r, xs) + | otherwise = second (x:) $ extractJust f xs + + +prettyShowShade' :: LtdErrorShow x => Shade' x -> String +prettyShowShade' sh = prettyShowsPrecShade' 0 sh [] + + + +wellDefinedShade' :: LinearSpace (Needle x) => Shade' x -> Maybe (Shade' x) +wellDefinedShade' (Shade' c e) = Shade' c <$> wellDefinedNorm e + + + +data LtdErrorShowWitness m where + LtdErrorShowWitness :: (LtdErrorShow (Interior m), LtdErrorShow (Needle m)) + => PseudoAffineWitness m -> LtdErrorShowWitness m + +class Refinable m => LtdErrorShow m where + ltdErrorShowWitness :: LtdErrorShowWitness m + default ltdErrorShowWitness :: (LtdErrorShow (Interior m), LtdErrorShow (Needle m)) + => LtdErrorShowWitness m + ltdErrorShowWitness = LtdErrorShowWitness pseudoAffineWitness + showsPrecShade'_errorLtdC :: Int -> Shade' m -> ShowS + prettyShowsPrecShade' :: Int -> Shade' m -> ShowS + prettyShowsPrecShade' p sh@(Shade' c e) + = showParen (p>6) $ v + . ("|±|["++) . flip (foldr id) (intersperse (',':) u) . (']':) + where v = showsPrecShade'_errorLtdC 6 sh + u :: [ShowS] = case ltdErrorShowWitness :: LtdErrorShowWitness m of + LtdErrorShowWitness (PseudoAffineWitness (SemimanifoldWitness _)) -> + [ showsPrecShade'_errorLtdC 6 (Shade' δ e :: Shade' (Needle m)) + | δ <- varianceSpanningSystem e'] + e' = dualNorm e + +instance LtdErrorShow ℝ⁰ where + showsPrecShade'_errorLtdC _ _ = ("zeroV"++) +instance LtdErrorShow ℝ where + showsPrecShade'_errorLtdC _ (Shade' v u) = errorLtdShow (δ/2) v + where δ = case u<$|1 of + σ | σ>0 -> sqrt $ 1/σ + _ -> v*10 +instance LtdErrorShow ℝ² where + showsPrecShade'_errorLtdC _ sh = ("V2 "++) . shshx . (' ':) . shshy + where shx = projectShade (lensEmbedding _x) sh :: Shade' ℝ + shy = projectShade (lensEmbedding _y) sh :: Shade' ℝ + shshx = showsPrecShade'_errorLtdC 0 shx + shshy = showsPrecShade'_errorLtdC 0 shy +instance LtdErrorShow ℝ³ where + showsPrecShade'_errorLtdC _ sh = ("V3 "++) . shshx . (' ':) . shshy . (' ':) . shshz + where shx = projectShade (lensEmbedding _x) sh :: Shade' ℝ + shy = projectShade (lensEmbedding _y) sh :: Shade' ℝ + shz = projectShade (lensEmbedding _z) sh :: Shade' ℝ + shshx = showsPrecShade'_errorLtdC 0 shx + shshy = showsPrecShade'_errorLtdC 0 shy + shshz = showsPrecShade'_errorLtdC 0 shz +instance LtdErrorShow ℝ⁴ where + showsPrecShade'_errorLtdC _ sh + = ("V4 "++) . shshx . (' ':) . shshy . (' ':) . shshz . (' ':) . shshw + where shx = projectShade (lensEmbedding _x) sh :: Shade' ℝ + shy = projectShade (lensEmbedding _y) sh :: Shade' ℝ + shz = projectShade (lensEmbedding _z) sh :: Shade' ℝ + shw = projectShade (lensEmbedding _w) sh :: Shade' ℝ + shshx = showsPrecShade'_errorLtdC 0 shx + shshy = showsPrecShade'_errorLtdC 0 shy + shshz = showsPrecShade'_errorLtdC 0 shz + shshw = showsPrecShade'_errorLtdC 0 shw +instance ∀ x y . + ( LtdErrorShow x, LtdErrorShow y + , Scalar (DualVector (Needle' x)) ~ Scalar (DualVector (Needle' y)) ) + => LtdErrorShow (x,y) where + ltdErrorShowWitness = case ( ltdErrorShowWitness :: LtdErrorShowWitness x + , ltdErrorShowWitness :: LtdErrorShowWitness y ) of + ( LtdErrorShowWitness(PseudoAffineWitness(SemimanifoldWitness BoundarylessWitness)) + , LtdErrorShowWitness(PseudoAffineWitness(SemimanifoldWitness BoundarylessWitness)) ) + ->LtdErrorShowWitness(PseudoAffineWitness(SemimanifoldWitness BoundarylessWitness)) + showsPrecShade'_errorLtdC _ sh = ('(':) . shshx . (',':) . shshy . (')':) + where (shx,shy) = factoriseShade sh + shshx = showsPrecShade'_errorLtdC 0 shx + shshy = showsPrecShade'_errorLtdC 0 shy + +instance LtdErrorShow x => Show (Shade' x) where + showsPrec = prettyShowsPrecShade' diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs --- a/Data/Manifold/TreeCover.hs +++ b/Data/Manifold/TreeCover.hs @@ -14,26 +14,17 @@ {-# LANGUAGE DeriveFunctor #-} {-# LANGUAGE DeriveFoldable #-} {-# LANGUAGE DeriveTraversable #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE TupleSections #-} -{-# LANGUAGE ParallelListComp #-} -{-# LANGUAGE MonadComprehensions #-} {-# LANGUAGE UnicodeSyntax #-} -{-# LANGUAGE ConstraintKinds #-} -{-# LANGUAGE PatternGuards #-} {-# LANGUAGE PatternSynonyms #-} -{-# LANGUAGE ViewPatterns #-} {-# LANGUAGE LambdaCase #-} {-# LANGUAGE TypeOperators #-} -{-# LANGUAGE CPP #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE LiberalTypeSynonyms #-} -{-# LANGUAGE RecordWildCards #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE TemplateHaskell #-} @@ -47,9 +38,10 @@ , fullShade, fullShade', pointsShades, pointsShade's , pointsCovers, pointsCover's, coverAllAround -- ** Evaluation - , occlusion + , occlusion, prettyShowsPrecShade', prettyShowShade' -- ** Misc , factoriseShade, intersectShade's, linIsoTransformShade + , embedShade, projectShade , Refinable, subShade', refineShade', convolveShade', coerceShade , mixShade's -- * Shade trees @@ -72,6 +64,8 @@ -- ** Triangulation-builders , TriangBuild, doTriangBuild , AutoTriang, breakdownAutoTriang + -- ** External + , AffineManifold, euclideanMetric ) where @@ -92,10 +86,13 @@ import Data.Tagged import Data.SimplicialComplex +import Data.Manifold.Shade import Data.Manifold.Types import Data.Manifold.Types.Primitive ((^), empty) import Data.Manifold.PseudoAffine import Data.Manifold.Riemannian +import Data.Manifold.Atlas +import Data.Function.Affine import Data.Embedding import Data.CoNat @@ -127,51 +124,7 @@ import Data.Type.Coercion --- | Possibly / Partially / asymPtotically singular metric. -data PSM x = PSM { - psmExpanse :: !(Metric' x) - , relevantEigenspan :: ![Needle' x] - } - --- | A 'Shade' is a very crude description of a region within a manifold. It --- can be interpreted as either an ellipsoid shape, or as the Gaussian peak --- of a normal distribution (use --- for actually sampling from that distribution). --- --- For a /precise/ description of an arbitrarily-shaped connected subset of a manifold, --- there is 'Region', whose implementation is vastly more complex. -data Shade x = Shade { _shadeCtr :: !(Interior x) - , _shadeExpanse :: !(Metric' x) } -deriving instance (Show (Interior x), Show (Metric' x), WithField ℝ PseudoAffine x) - => Show (Shade x) - --- | A “co-shade” can describe ellipsoid regions as well, but unlike --- 'Shade' it can be unlimited / infinitely wide in some directions. --- It does OTOH need to have nonzero thickness, which 'Shade' needs not. -data Shade' x = Shade' { _shade'Ctr :: !(Interior x) - , _shade'Narrowness :: !(Metric x) } -deriving instance (Show (Interior x), Show (Metric x), WithField ℝ PseudoAffine x) - => Show (Shade' x) - -data LocalDifferentialEqn x y = LocalDifferentialEqn { - _predictDerivatives :: Maybe (Shade' (LocalLinear x y)) - , _rescanDerivatives :: Shade' (LocalLinear x y) -> Shade' y -> Maybe (Shade' y) - } -makeLenses ''LocalDifferentialEqn - -type DifferentialEqn x y = Shade (x,y) -> LocalDifferentialEqn x y - -data LocalDataPropPlan x y = LocalDataPropPlan - { _sourcePosition :: !(Interior x) - , _targetPosOffset :: !(Needle x) - , _sourceData, _targetAPrioriData :: !y - , _relatedData :: [(Needle x, y)] - } -deriving instance (Show (Interior x), Show y, Show (Needle x)) => Show (LocalDataPropPlan x y) - -makeLenses ''LocalDataPropPlan - type Depth = Int data Wall x = Wall { _wallID :: (Depth,(Int,Int)) , _wallAnchor :: Interior x @@ -181,182 +134,6 @@ makeLenses ''Wall -class IsShade shade where --- type (*) shade :: *->* - -- | Access the center of a 'Shade' or a 'Shade''. - shadeCtr :: Lens' (shade x) (Interior x) --- -- | Convert between 'Shade' and 'Shade' (which must be neither singular nor infinite). --- unsafeDualShade :: WithField ℝ Manifold x => shade x -> shade* x - -- | Check the statistical likelihood-density of a point being within a shade. - -- This is taken as a normal distribution. - occlusion :: ( PseudoAffine x, SimpleSpace (Needle x) - , s ~ (Scalar (Needle x)), RealDimension s ) - => shade x -> x -> s - factoriseShade :: ( Manifold x, SimpleSpace (Needle x) - , Manifold y, SimpleSpace (Needle y) - , Scalar (Needle x) ~ Scalar (Needle y) ) - => shade (x,y) -> (shade x, shade y) - coerceShade :: (Manifold x, Manifold y, LocallyCoercible x y) => shade x -> shade y - linIsoTransformShade :: ( LinearManifold x, LinearManifold y - , SimpleSpace x, SimpleSpace y, Scalar x ~ Scalar y ) - => (x+>y) -> shade x -> shade y - -instance IsShade Shade where - shadeCtr f (Shade c e) = fmap (`Shade`e) $ f c - occlusion = occ pseudoAffineWitness dualSpaceWitness - where occ :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x) - , Scalar (Needle x) ~ s, RealDimension s ) - => PseudoAffineWitness x -> DualNeedleWitness x -> Shade x -> x -> s - occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness (Shade p₀ δ) - = \p -> case toInterior p >>= (.-~.p₀) of - (Just vd) | mSq <- normSq δinv vd - , mSq == mSq -- avoid NaN - -> exp (negate mSq) - _ -> zeroV - where δinv = dualNorm δ - factoriseShade = fs dualSpaceWitness dualSpaceWitness - where fs :: ∀ x y . ( Manifold x, SimpleSpace (Needle x) - , Manifold y, SimpleSpace (Needle y) - , Scalar (Needle x) ~ Scalar (Needle y) ) - => DualNeedleWitness x -> DualNeedleWitness y - -> Shade (x,y) -> (Shade x, Shade y) - fs DualSpaceWitness DualSpaceWitness (Shade (x₀,y₀) δxy) - = (Shade x₀ δx, Shade y₀ δy) - where (δx,δy) = summandSpaceNorms δxy - coerceShade = cS dualSpaceWitness dualSpaceWitness - where cS :: ∀ x y . (LocallyCoercible x y) - => DualNeedleWitness x -> DualNeedleWitness y -> Shade x -> Shade y - cS DualSpaceWitness DualSpaceWitness - = \(Shade x δxym) -> Shade (internCoerce x) (tN δxym) - where tN = case oppositeLocalCoercion :: CanonicalDiffeomorphism y x of - CanonicalDiffeomorphism -> - transformNorm . arr $ coerceNeedle' ([]::[(y,x)]) - internCoerce = case interiorLocalCoercion ([]::[(x,y)]) of - CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism - linIsoTransformShade = lits dualSpaceWitness dualSpaceWitness - where lits :: ∀ x y . ( LinearManifold x, LinearManifold y - , Scalar (Needle x) ~ Scalar (Needle y) ) - => DualSpaceWitness x -> DualSpaceWitness y - -> (x+>y) -> Shade x -> Shade y - lits DualSpaceWitness DualSpaceWitness f (Shade x δx) - = Shade (f $ x) (transformNorm (adjoint $ f) δx) - -instance ImpliesMetric Shade where - type MetricRequirement Shade x = (Manifold x, SimpleSpace (Needle x)) - inferMetric' (Shade _ e) = e - inferMetric = im dualSpaceWitness - where im :: (Manifold x, SimpleSpace (Needle x)) - => DualNeedleWitness x -> Shade x -> Metric x - im DualSpaceWitness (Shade _ e) = dualNorm e - -instance ImpliesMetric Shade' where - type MetricRequirement Shade' x = (Manifold x, SimpleSpace (Needle x)) - inferMetric (Shade' _ e) = e - inferMetric' (Shade' _ e) = dualNorm e - -shadeExpanse :: Lens' (Shade x) (Metric' x) -shadeExpanse f (Shade c e) = fmap (Shade c) $ f e - -instance IsShade Shade' where - shadeCtr f (Shade' c e) = fmap (`Shade'`e) $ f c - occlusion = occ pseudoAffineWitness - where occ :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x) - , Scalar (Needle x) ~ s, RealDimension s ) - => PseudoAffineWitness x -> Shade' x -> x -> s - occ (PseudoAffineWitness (SemimanifoldWitness _)) (Shade' p₀ δinv) p - = case toInterior p >>= (.-~.p₀) of - (Just vd) | mSq <- normSq δinv vd - , mSq == mSq -- avoid NaN - -> exp (negate mSq) - _ -> zeroV - factoriseShade (Shade' (x₀,y₀) δxy) = (Shade' x₀ δx, Shade' y₀ δy) - where (δx,δy) = summandSpaceNorms δxy - coerceShade = cS - where cS :: ∀ x y . (LocallyCoercible x y) => Shade' x -> Shade' y - cS = \(Shade' x δxym) -> Shade' (internCoerce x) (tN δxym) - where tN = case oppositeLocalCoercion :: CanonicalDiffeomorphism y x of - CanonicalDiffeomorphism -> - transformNorm . arr $ coerceNeedle ([]::[(y,x)]) - internCoerce = case interiorLocalCoercion ([]::[(x,y)]) of - CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism - linIsoTransformShade f (Shade' x δx) - = Shade' (f $ x) (transformNorm (pseudoInverse f) δx) - -shadeNarrowness :: Lens' (Shade' x) (Metric x) -shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e - -instance ∀ x . (PseudoAffine x) => Semimanifold (Shade x) where - type Needle (Shade x) = Needle x - fromInterior = id - toInterior = pure - translateP = Tagged (.+~^) - (.+~^) = case semimanifoldWitness :: SemimanifoldWitness x of - SemimanifoldWitness BoundarylessWitness - -> \(Shade c e) v -> Shade (c.+~^v) e - (.-~^) = case semimanifoldWitness :: SemimanifoldWitness x of - SemimanifoldWitness BoundarylessWitness - -> \(Shade c e) v -> Shade (c.-~^v) e - semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of - (SemimanifoldWitness BoundarylessWitness) - -> SemimanifoldWitness BoundarylessWitness - -instance (WithField ℝ PseudoAffine x, Geodesic (Interior x), SimpleSpace (Needle x)) - => Geodesic (Shade x) where - geodesicBetween = gb dualSpaceWitness - where gb :: DualNeedleWitness x -> Shade x -> Shade x -> Maybe (D¹ -> Shade x) - gb DualSpaceWitness (Shade c (Norm e)) (Shade ζ (Norm η)) = pure interp - where interp t@(D¹ q) = Shade (pinterp t) - (Norm . arr . lerp ed ηd $ (q+1)/2) - ed@(LinearMap _) = arr e - ηd@(LinearMap _) = arr η - Just pinterp = geodesicBetween c ζ - -instance (AffineManifold x) => Semimanifold (Shade' x) where - type Needle (Shade' x) = Diff x - fromInterior = id - toInterior = pure - translateP = Tagged (.+~^) - Shade' c e .+~^ v = Shade' (c.+^v) e - Shade' c e .-~^ v = Shade' (c.-^v) e - -instance (WithField ℝ AffineManifold x, Geodesic x, SimpleSpace (Needle x)) - => Geodesic (Shade' x) where - geodesicBetween (Shade' c e) (Shade' ζ η) = pure interp - where sharedSpan = sharedNormSpanningSystem e η - interp t = Shade' (pinterp t) - (spanNorm [ v ^/ (alerpB 1 (recip qη) t) - | (v,qη) <- sharedSpan ]) - Just pinterp = geodesicBetween c ζ - -fullShade :: WithField ℝ PseudoAffine x => Interior x -> Metric' x -> Shade x -fullShade ctr expa = Shade ctr expa - -fullShade' :: WithField ℝ PseudoAffine x => Interior x -> Metric x -> Shade' x -fullShade' ctr expa = Shade' ctr expa - - --- | Span a 'Shade' from a center point and multiple deviation-vectors. -#if GLASGOW_HASKELL < 800 -pattern (:±) :: () -#else -pattern (:±) :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) -#endif - => (WithField ℝ Manifold x, SimpleSpace (Needle x)) - => Interior x -> [Needle x] -> Shade x -pattern x :± shs <- Shade x (varianceSpanningSystem -> shs) - where x :± shs = fullShade x $ spanVariance shs - --- | Similar to ':±', but instead of expanding the shade, each vector /restricts/ it. --- Iff these form a orthogonal basis (in whatever sense applicable), then both --- methods will be equivalent. --- --- Note that '|±|' is only possible, as such, in an inner-product space; in --- general you need reciprocal vectors ('Needle'') to define a 'Shade''. -(|±|) :: WithField ℝ EuclidSpace x => x -> [Needle x] -> Shade' x -x |±| shs = Shade' x $ spanNorm [v^/(v<.>v) | v<-shs] - - - subshadeId' :: ∀ x . (WithField ℝ PseudoAffine x, LinearSpace (Needle x)) => x -> NonEmpty (Needle' x) -> x -> (Int, HourglassBulb) subshadeId' c expvs x = case ( dualSpaceWitness :: DualNeedleWitness x @@ -375,220 +152,9 @@ --- | Attempt to find a 'Shade' that describes the distribution of given points. --- At least in an affine space (and thus locally in any manifold), this can be used to --- estimate the parameters of a normal distribution from which some points were --- sampled. Note that some points will be “outside” of the shade, --- as happens for a normal distribution with some statistical likelyhood. --- (Use 'pointsCovers' if you need to prevent that.) --- --- For /nonconnected/ manifolds it will be necessary to yield separate shades --- for each connected component. And for an empty input list, there is no shade! --- Hence the result type is a list. -pointsShades :: (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => [Interior x] -> [Shade x] -pointsShades = map snd . pointsShades' mempty . map fromInterior -coverAllAround :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => Interior x -> [Needle x] -> Shade x -coverAllAround x₀ offs = Shade x₀ - $ guaranteeIn dualSpaceWitness offs - (scaleNorm (1/fromIntegral (length offs)) $ spanVariance offs) - where guaranteeIn :: DualNeedleWitness x -> [Needle x] -> Metric' x -> Metric' x - guaranteeIn w@DualSpaceWitness offs ex - = case offs >>= \v -> guard ((ex'|$|v) > 1) >> [(v, spanVariance [v])] of - [] -> ex - outs -> guaranteeIn w (fst<$>outs) - ( densifyNorm $ - ex <> scaleNorm - (sqrt . recip . fromIntegral - $ 2 * length outs) - (mconcat $ snd<$>outs) - ) - where ex' = dualNorm ex --- | Like 'pointsShades', but ensure that all points are actually in --- the shade, i.e. if @['Shade' x₀ ex]@ is the result then --- @'metric' (recipMetric ex) (p-x₀) ≤ 1@ for all @p@ in the list. -pointsCovers :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => [Interior x] -> [Shade x] -pointsCovers = case pseudoAffineWitness :: PseudoAffineWitness x of - (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) -> - \ps -> map (\(ps', Shade x₀ _) - -> coverAllAround x₀ [v | p<-ps' - , let Just v - = p.-~.fromInterior x₀]) - (pointsShades' mempty (fromInterior<$>ps) :: [([x], Shade x)]) -pointsShade's :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => [Interior x] -> [Shade' x] -pointsShade's = case dualSpaceWitness :: DualNeedleWitness x of - DualSpaceWitness -> map (\(Shade c e :: Shade x) -> Shade' c $ dualNorm e) . pointsShades - -pointsCover's :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => [Interior x] -> [Shade' x] -pointsCover's = case dualSpaceWitness :: DualNeedleWitness x of - DualSpaceWitness -> map (\(Shade c e :: Shade x) -> Shade' c $ dualNorm e) . pointsCovers - -pseudoECM :: ∀ x p . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x), Hask.Functor p) - => p x -> NonEmpty x -> (x, ([x],[x])) -pseudoECM = case semimanifoldWitness :: SemimanifoldWitness x of - SemimanifoldWitness _ -> - \_ (p₀ NE.:| psr) -> foldl' ( \(acc, (rb,nr)) (i,p) - -> case (p.-~.acc, toInterior acc) of - (Just δ, Just acci) - -> (acci .+~^ δ^/i, (p:rb, nr)) - _ -> (acc, (rb, p:nr)) ) - (p₀, mempty) - ( zip [1..] $ p₀:psr ) - -pointsShades' :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => Metric' x -> [x] -> [([x], Shade x)] -pointsShades' _ [] = [] -pointsShades' minExt ps = case (expa, toInterior ctr) of - (Just e, Just c) - -> (ps, fullShade c e) : pointsShades' minExt unreachable - _ -> pointsShades' minExt inc'd - ++ pointsShades' minExt unreachable - where (ctr,(inc'd,unreachable)) = pseudoECM ([]::[x]) $ NE.fromList ps - expa = ( (<>minExt) . spanVariance . map (^/ fromIntegral (length ps)) ) - <$> mapM (.-~.ctr) ps - - --- | Attempt to reduce the number of shades to fewer (ideally, a single one). --- In the simplest cases these should guaranteed cover the same area; --- for non-flat manifolds it only works in a heuristic sense. -shadesMerge :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x)) - => ℝ -- ^ How near (inverse normalised distance, relative to shade expanse) - -- two shades must be to be merged. If this is zero, any shades - -- in the same connected region of a manifold are merged. - -> [Shade x] -- ^ A list of /n/ shades. - -> [Shade x] -- ^ /m/ ≤ /n/ shades which cover at least the same area. -shadesMerge fuzz (sh₁@(Shade c₁ e₁) : shs) - = case extractJust (tryMerge pseudoAffineWitness dualSpaceWitness) - shs of - (Just mg₁, shs') -> shadesMerge fuzz - $ shs'++[mg₁] -- Append to end to prevent undue weighting - -- of first shade and its mergers. - (_, shs') -> sh₁ : shadesMerge fuzz shs' - where tryMerge :: PseudoAffineWitness x -> DualNeedleWitness x - -> Shade x -> Maybe (Shade x) - tryMerge (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) DualSpaceWitness - (Shade c₂ e₂) - | Just v <- c₁.-~.c₂ - , [e₁',e₂'] <- dualNorm<$>[e₁, e₂] - , b₁ <- e₂'|$|v - , b₂ <- e₁'|$|v - , fuzz*b₁*b₂ <= b₁ + b₂ - = Just $ let cc = c₂ .+~^ v ^/ 2 - Just cv₁ = c₁.-~.cc - Just cv₂ = c₂.-~.cc - in Shade cc $ e₁ <> e₂ <> spanVariance [cv₁, cv₂] - | otherwise = Nothing -shadesMerge _ shs = shs - --- | Weakened version of 'intersectShade's'. What this function calculates is --- rather the /weighted mean/ of ellipsoid regions. If you interpret the --- shades as uncertain physical measurements with normal distribution, --- it gives the maximum-likelyhood result for multiple measurements of the --- same quantity. -mixShade's :: ∀ y . (WithField ℝ Manifold y, SimpleSpace (Needle y)) - => NonEmpty (Shade' y) -> Maybe (Shade' y) -mixShade's = ms pseudoAffineWitness dualSpaceWitness - where ms :: PseudoAffineWitness y -> DualNeedleWitness y - -> NonEmpty (Shade' y) -> Maybe (Shade' y) - ms (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness - (Shade' c₀ (Norm e₁):|shs) = sequenceA ciso >> pure mixed - where ciso = [ci.-~.c₀ | Shade' ci shi <- shs] - cis = [v | Just v <- ciso] - σe = arr . sumV $ e₁ : (applyNorm . _shade'Narrowness<$>shs) - cc = σe \$ sumV [ei $ ci | ci <- cis - | Shade' _ (Norm ei) <- shs] - mixed = Shade' (c₀+^cc) $ densifyNorm ( mconcat - [ Norm $ ei ^/ (1+(normSq ni $ ci^-^cc)) - | ni@(Norm ei) <- Norm e₁ : (_shade'Narrowness<$>shs) - | ci <- zeroV : cis - ] ) - Tagged (+^) = translateP :: Tagged y (Interior y->Needle y->Interior y) - -- cc should minimise the quadratic form - -- β(cc) = ∑ᵢ ⟨cc−cᵢ|eᵢ|cc−cᵢ⟩ - -- = ⟨cc|e₁|cc⟩ + ∑ᵢ₌₁… ⟨cc−c₂|e₂|cc−c₂⟩ - -- = ⟨cc|e₁|cc⟩ + ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ − 2⋅⟨cᵢ|eᵢ|cc⟩ + ⟨cᵢ|eᵢ|cᵢ⟩ ) - -- It is thus - -- β(cc + δ⋅v) − β cc - -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ - -- + ∑ᵢ₌₁…( ⟨cc + δ⋅v|eᵢ|cc + δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|cc + δ⋅v⟩ + ⟨cᵢ|eᵢ|cᵢ⟩ ) - -- − ⟨cc|e₁|cc⟩ - -- − ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ + 2⋅⟨cᵢ|eᵢ|cc⟩ − ⟨cᵢ|eᵢ|cᵢ⟩ ) - -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ - -- + ∑ᵢ₌₁…( ⟨cc + δ⋅v|eᵢ|cc + δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|δ⋅v⟩ ) - -- − ⟨cc|e₁|cc⟩ - -- − ∑ᵢ₌₁…( ⟨cc|eᵢ|cc⟩ ) - -- = 2⋅⟨δ⋅v|e₁|cc⟩ + ⟨δ⋅v|e₁|δ⋅v⟩ - -- + ∑ᵢ₌₁…( 2⋅⟨δ⋅v|eᵢ|cc⟩ + ⟨δ⋅v|eᵢ|δ⋅v⟩ − 2⋅⟨cᵢ|eᵢ|δ⋅v⟩ ) - -- = 2⋅⟨δ⋅v|∑ᵢeᵢ|cc⟩ − 2⋅∑ᵢ₌₁… ⟨cᵢ|eᵢ|δ⋅v⟩ + 𝓞(δ²) - -- This should vanish for all v, which is fulfilled by - -- (∑ᵢeᵢ)|cc⟩ = ∑ᵢ₌₁… eᵢ|cᵢ⟩. - --- | Evaluate the shade as a quadratic form; essentially --- @ --- minusLogOcclusion sh x = x <.>^ (sh^.shadeExpanse $ x - sh^.shadeCtr) --- @ --- where 'shadeExpanse' gives a metric (matrix) that characterises the --- width of the shade. -minusLogOcclusion' :: ∀ x s . ( PseudoAffine x, LinearSpace (Needle x) - , s ~ (Scalar (Needle x)), RealDimension s ) - => Shade' x -> x -> s -minusLogOcclusion' (Shade' p₀ δinv) - = occ (pseudoAffineWitness :: PseudoAffineWitness x) - (dualSpaceWitness :: DualNeedleWitness x) - where occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness - p = case toInterior p >>= (.-~.p₀) of - (Just vd) | mSq <- normSq δinv vd - , mSq == mSq -- avoid NaN - -> mSq - _ -> 1/0 -minusLogOcclusion :: ∀ x s . ( PseudoAffine x, SimpleSpace (Needle x) - , s ~ (Scalar (Needle x)), RealDimension s ) - => Shade x -> x -> s -minusLogOcclusion (Shade p₀ δ) - = occ (pseudoAffineWitness :: PseudoAffineWitness x) - (dualSpaceWitness :: DualNeedleWitness x) - where occ (PseudoAffineWitness (SemimanifoldWitness _)) DualSpaceWitness - = \p -> case toInterior p >>= (.-~.p₀) of - (Just vd) | mSq <- normSq δinv vd - , mSq == mSq -- avoid NaN - -> mSq - _ -> 1/0 - where δinv = dualNorm δ - - - - -rangeOnGeodesic :: ∀ i m . - ( WithField ℝ PseudoAffine m, Geodesic m, SimpleSpace (Needle m) - , WithField ℝ IntervalLike i, SimpleSpace (Needle i) ) - => m -> m -> Maybe (Shade i -> Shade m) -rangeOnGeodesic = case ( semimanifoldWitness :: SemimanifoldWitness i - , dualSpaceWitness :: DualNeedleWitness i - , dualSpaceWitness :: DualNeedleWitness m ) of - (SemimanifoldWitness _, DualSpaceWitness, DualSpaceWitness) -> - \p₀ p₁ -> (`fmap`(geodesicBetween p₀ p₁)) - $ \interp -> \(Shade t₀ et) - -> case pointsShades - . mapMaybe (toInterior - . interp . (toClosedInterval :: i -> D¹)) - $ fromInterior <$> t₀ : [ t₀+^v - | v<-normSpanningSystem et ] of - [sh] -> sh - _ -> case pointsShades $ mapMaybe (toInterior . interp . D¹) - [-0.999, 0.999] of - [sh] -> sh - where Tagged (+^) = translateP :: Tagged i (Interior i->Needle i->Interior i) - - - - -- | Hourglass as the geometric shape (two opposing ~conical volumes, sharing -- only a single point in the middle); has nothing to do with time. data Hourglass s = Hourglass { upperBulb, lowerBulb :: !s } @@ -712,6 +278,8 @@ = OverlappingBranches n (sh.+~^v) $ fmap (\(DBranch d c) -> DBranch d $ (.+~^v)<$>c) br DisjointBranches n br .+~^ v = DisjointBranches n $ (.+~^v)<$>br + semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of + SemimanifoldWitness BoundarylessWitness -> SemimanifoldWitness BoundarylessWitness -- | WRT union. instance (WithField ℝ Manifold x, SimpleSpace (Needle x)) => Semigroup (ShadeTree x) where @@ -976,299 +544,8 @@ coerceShade --- | Class of manifolds which can use 'Shade'' as a basic set type. --- This is easily possible for vector spaces with the default implementations. -class (WithField ℝ PseudoAffine y, SimpleSpace (Needle y)) => Refinable y where - -- | @a `subShade'` b ≡ True@ means @a@ is fully contained in @b@, i.e. from - -- @'minusLogOcclusion'' a p < 1@ follows also @minusLogOcclusion' b p < 1@. - subShade' :: Shade' y -> Shade' y -> Bool - subShade' (Shade' ac ae) (Shade' tc te) - = case pseudoAffineWitness :: PseudoAffineWitness y of - PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) - | Just v <- tc.-~.ac - , v² <- normSq te v - , v² <= 1 - -> all (\(y',μ) -> case μ of - Nothing -> True -- 'te' has infinite extension in this direction - Just ξ - | ξ<1 -> False -- 'ae' would be vaster than 'te' in this direction - | ω <- abs $ y'<.>^v - -> (ω + 1/ξ)^2 <= 1 - v² + ω^2 - -- See @images/constructions/subellipse-check-heuristic.svg@ - ) $ sharedSeminormSpanningSystem te ae - _ -> False - - -- | Intersection between two shades. - refineShade' :: Shade' y -> Shade' y -> Maybe (Shade' y) - refineShade' (Shade' c₀ (Norm e₁)) (Shade' c₀₂ (Norm e₂)) - = case ( dualSpaceWitness :: DualNeedleWitness y - , pseudoAffineWitness :: PseudoAffineWitness y ) of - (DualSpaceWitness, PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) - -> do - c₂ <- c₀₂.-~.c₀ - let σe = arr $ e₁^+^e₂ - e₁c₂ = e₁ $ c₂ - e₂c₂ = e₂ $ c₂ - cc = σe \$ e₂c₂ - cc₂ = cc ^-^ c₂ - e₁cc = e₁ $ cc - e₂cc = e₂ $ cc - α = 2 + e₂c₂<.>^cc₂ - guard (α > 0) - let ee = σe ^/ α - c₂e₁c₂ = e₁c₂<.>^c₂ - c₂e₂c₂ = e₂c₂<.>^c₂ - c₂eec₂ = (c₂e₁c₂ + c₂e₂c₂) / α - return $ case middle . sort - $ quadraticEqnSol c₂e₁c₂ - (2 * (e₁cc<.>^c₂)) - (e₁cc<.>^cc - 1) - ++quadraticEqnSol c₂e₂c₂ - (2 * (e₂cc<.>^c₂ - c₂e₂c₂)) - (e₂cc<.>^cc - 2 * (e₂c₂<.>^cc) + c₂e₂c₂ - 1) of - [γ₁,γ₂] | abs (γ₁+γ₂) < 2 -> let - cc' = cc ^+^ ((γ₁+γ₂)/2)*^c₂ - rγ = abs (γ₁ - γ₂) / 2 - η = if rγ * c₂eec₂ /= 0 && 1 - rγ^2 * c₂eec₂ > 0 - then sqrt (1 - rγ^2 * c₂eec₂) / (rγ * c₂eec₂) - else 0 - in Shade' (c₀.+~^cc') - (Norm (arr ee) <> spanNorm [ee $ c₂^*η]) - _ -> Shade' (c₀.+~^cc) (Norm $ arr ee) - where quadraticEqnSol a b c - | a == 0, b /= 0 = [-c/b] - | a /= 0 && disc == 0 = [- b / (2*a)] - | a /= 0 && disc > 0 = [ (σ * sqrt disc - b) / (2*a) - | σ <- [-1, 1] ] - | otherwise = [] - where disc = b^2 - 4*a*c - middle (_:x:y:_) = [x,y] - middle l = l - -- ⟨x−c₁|e₁|x−c₁⟩ < 1 ∧ ⟨x−c₂|e₂|x−c₂⟩ < 1 - -- We search (cc,ee) such that this implies - -- ⟨x−cc|ee|x−cc⟩ < 1. - -- Let WLOG c₁ = 0, so - -- ⟨x|e₁|x⟩ < 1. - -- cc should minimise the quadratic form - -- β(cc) = ⟨cc−c₁|e₁|cc−c₁⟩ + ⟨cc−c₂|e₂|cc−c₂⟩ - -- = ⟨cc|e₁|cc⟩ + ⟨cc−c₂|e₂|cc−c₂⟩ - -- = ⟨cc|e₁|cc⟩ + ⟨cc|e₂|cc⟩ − 2⋅⟨c₂|e₂|cc⟩ + ⟨c₂|e₂|c₂⟩ - -- It is thus - -- β(cc + δ⋅v) − β cc - -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + ⟨cc + δ⋅v|e₂|cc + δ⋅v⟩ − 2⋅⟨c₂|e₂|cc + δ⋅v⟩ + ⟨c₂|e₂|c₂⟩ - -- − ⟨cc|e₁|cc⟩ − ⟨cc|e₂|cc⟩ + 2⋅⟨c₂|e₂|cc⟩ − ⟨c₂|e₂|c₂⟩ - -- = ⟨cc + δ⋅v|e₁|cc + δ⋅v⟩ + ⟨cc + δ⋅v|e₂|cc + δ⋅v⟩ − 2⋅⟨c₂|e₂|δ⋅v⟩ - -- − ⟨cc|e₁|cc⟩ − ⟨cc|e₂|cc⟩ - -- = 2⋅⟨δ⋅v|e₁|cc⟩ + ⟨δ⋅v|e₁|δ⋅v⟩ + 2⋅⟨δ⋅v|e₂|cc⟩ + ⟨δ⋅v|e₂|δ⋅v⟩ − 2⋅⟨c₂|e₂|δ⋅v⟩ - -- = 2⋅δ⋅⟨v|e₁+e₂|cc⟩ − 2⋅δ⋅⟨v|e₂|c₂⟩ + 𝓞(δ²) - -- This should vanish for all v, which is fulfilled by - -- (e₁+e₂)|cc⟩ = e₂|c₂⟩. - -- - -- If we now choose - -- ee = (e₁+e₂) / α - -- then - -- ⟨x−cc|ee|x−cc⟩ ⋅ α - -- = ⟨x−cc|ee|x⟩ ⋅ α − ⟨x−cc|ee|cc⟩ ⋅ α - -- = ⟨x|ee|x−cc⟩ ⋅ α − ⟨x−cc|e₂|c₂⟩ - -- = ⟨x|ee|x⟩ ⋅ α − ⟨x|ee|cc⟩ ⋅ α − ⟨x−cc|e₂|c₂⟩ - -- = ⟨x|e₁+e₂|x⟩ − ⟨x|e₂|c₂⟩ − ⟨x−cc|e₂|c₂⟩ - -- = ⟨x|e₁|x⟩ + ⟨x|e₂|x⟩ − ⟨x|e₂|c₂⟩ − ⟨x−cc|e₂|c₂⟩ - -- < 1 + ⟨x|e₂|x−c₂⟩ − ⟨x−cc|e₂|c₂⟩ - -- = 1 + ⟨x−c₂|e₂|x−c₂⟩ + ⟨c₂|e₂|x−c₂⟩ − ⟨x−cc|e₂|c₂⟩ - -- < 2 + ⟨x−c₂−x+cc|e₂|c₂⟩ - -- = 2 + ⟨cc−c₂|e₂|c₂⟩ - -- Really we want - -- ⟨x−cc|ee|x−cc⟩ ⋅ α < α - -- So choose α = 2 + ⟨cc−c₂|e₂|c₂⟩. - -- - -- The ellipsoid "cc±√ee" captures perfectly the intersection - -- of the boundary of the shades, but it tends to significantly - -- overshoot the interior intersection in perpendicular direction, - -- i.e. in direction of c₂−c₁. E.g. - -- https://github.com/leftaroundabout/manifolds/blob/bc0460b9/manifolds/images/examples/ShadeCombinations/EllipseIntersections.png - -- 1. Really, the relevant points are those where either of the - -- intersector badnesses becomes 1. The intersection shade should - -- be centered between those points. We perform according corrections, - -- but only in c₂ direction, so this can be handled efficiently - -- as a 1D quadratic equation. - -- Consider - -- dⱼ c := ⟨c−cⱼ|eⱼ|c−cⱼ⟩ =! 1 - -- dⱼ (cc + γ⋅c₂) - -- = ⟨cc+γ⋅c₂−cⱼ|eⱼ|cc+γ⋅c₂−cⱼ⟩ - -- = ⟨cc−cⱼ|eⱼ|cc−cⱼ⟩ + 2⋅γ⋅⟨c₂|eⱼ|cc−cⱼ⟩ + γ²⋅⟨c₂|eⱼ|c₂⟩ - -- =! 1 - -- So - -- γⱼ = (- b ± √(b²−4⋅a⋅c)) / 2⋅a - -- where a = ⟨c₂|eⱼ|c₂⟩ - -- b = 2 ⋅ (⟨c₂|eⱼ|cc⟩ − ⟨c₂|eⱼ|cⱼ⟩) - -- c = ⟨cc|eⱼ|cc⟩ − 2⋅⟨cc|eⱼ|cⱼ⟩ + ⟨cⱼ|eⱼ|cⱼ⟩ − 1 - -- The ± sign should be chosen to get the smaller |γ| (otherwise - -- we end up on the wrong side of the shade), i.e. - -- γⱼ = (sgn bⱼ ⋅ √(bⱼ²−4⋅aⱼ⋅cⱼ) − bⱼ) / 2⋅aⱼ - -- 2. Trim the result in that direction to the actual - -- thickness of the lens-shaped intersection: we want - -- ⟨rγ⋅c₂|ee'|rγ⋅c₂⟩ = 1 - -- for a squeezed version of ee, - -- ee' = ee + ee|η⋅c₂⟩⟨η⋅c₂|ee - -- ee' = ee + η² ⋅ ee|c₂⟩⟨c₂|ee - -- ⟨rγ⋅c₂|ee'|rγ⋅c₂⟩ - -- = rγ² ⋅ (⟨c₂|ee|c₂⟩ + η² ⋅ ⟨c₂|ee|c₂⟩²) - -- = rγ² ⋅ ⟨c₂|ee|c₂⟩ + η² ⋅ rγ² ⋅ ⟨c₂|ee|c₂⟩² - -- η² = (1 − rγ²⋅⟨c₂|ee|c₂⟩) / (rγ² ⋅ ⟨c₂|ee|c₂⟩²) - -- η = √(1 − rγ²⋅⟨c₂|ee|c₂⟩) / (rγ ⋅ ⟨c₂|ee|c₂⟩) - -- With ⟨c₂|ee|c₂⟩ = (⟨c₂|e₁|c₂⟩ + ⟨c₂|e₂|c₂⟩)/α. - - -- | If @p@ is in @a@ (red) and @δ@ is in @b@ (green), - -- then @p.+~^δ@ is in @convolveShade' a b@ (blue). - -- --- Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/ShadeCombinations.ipynb#shadeConvolutions --- --- <> - convolveMetric :: Hask.Functor p => p y -> Metric y -> Metric y -> Metric y - convolveMetric _ ey eδ = spanNorm [ f ^* ζ crl - | (f,crl) <- eδsp ] - where eδsp = sharedSeminormSpanningSystem ey eδ - ζ = case filter (>0) . catMaybes $ snd<$>eδsp of - [] -> const 0 - nzrelap - -> let cre₁ = 1/minimum nzrelap - cre₂ = maximum nzrelap - edgeFactor = sqrt ( (1 + cre₁)^2 + (1 + cre₂)^2 ) - / (sqrt (1 + cre₁^2) + sqrt (1 + cre₂^2)) - in \case - Nothing -> 0 - Just 0 -> 0 - Just sq -> edgeFactor / (recip sq + 1) - - convolveShade' :: Shade' y -> Shade' (Needle y) -> Shade' y - convolveShade' = defaultConvolveShade' - -defaultConvolveShade' :: ∀ y . Refinable y => Shade' y -> Shade' (Needle y) -> Shade' y -defaultConvolveShade' = case (pseudoAffineWitness :: PseudoAffineWitness y) of - PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) - -> \(Shade' y₀ ey) (Shade' δ₀ eδ) -> Shade' (y₀.+~^δ₀) - $ convolveMetric ([]::[y]) ey eδ -instance Refinable ℝ where - refineShade' (Shade' cl el) (Shade' cr er) - = case (normSq el 1, normSq er 1) of - (0, _) -> return $ Shade' cr er - (_, 0) -> return $ Shade' cl el - (ql,qr) | ql>0, qr>0 - -> let [rl,rr] = sqrt . recip <$> [ql,qr] - b = maximum $ zipWith (-) [cl,cr] [rl,rr] - t = minimum $ zipWith (+) [cl,cr] [rl,rr] - in guard (b> - let cm = (b+t)/2 - rm = (t-b)/2 - in return $ Shade' cm (spanNorm [recip rm]) --- convolveShade' (Shade' y₀ ey) (Shade' δ₀ eδ) --- = case (metricSq ey 1, metricSq eδ 1) of --- (wy,wδ) | wy>0, wδ>0 --- -> Shade' (y₀.+~^δ₀) --- ( projector . recip --- $ recip (sqrt wy) + recip (sqrt wδ) ) --- (_ , _) -> Shade' y₀ zeroV - -instance ( Refinable a, Interior a ~ a, Refinable b, Interior b ~ b - , Scalar (DualVector (DualVector (Needle b))) - ~ Scalar (DualVector (DualVector (Needle a))) ) - => Refinable (a,b) - -instance Refinable ℝ⁰ -instance Refinable ℝ¹ -instance Refinable ℝ² -instance Refinable ℝ³ -instance Refinable ℝ⁴ - -instance ( SimpleSpace a, SimpleSpace b - , Scalar a ~ ℝ, Scalar b ~ ℝ - , Scalar (DualVector a) ~ ℝ, Scalar (DualVector b) ~ ℝ - , Scalar (DualVector (DualVector a)) ~ ℝ, Scalar (DualVector (DualVector b)) ~ ℝ ) - => Refinable (LinearMap ℝ a b) - -intersectShade's :: ∀ y . Refinable y => NonEmpty (Shade' y) -> Maybe (Shade' y) -intersectShade's (sh:|shs) = Hask.foldrM refineShade' sh shs - - -estimateLocalJacobian :: ∀ x y . ( WithField ℝ Manifold x, Refinable y - , SimpleSpace (Needle x), SimpleSpace (Needle y) ) - => Metric x -> [(Local x, Shade' y)] - -> Maybe (Shade' (LocalLinear x y)) -estimateLocalJacobian = elj ( pseudoAffineWitness :: PseudoAffineWitness x - , pseudoAffineWitness :: PseudoAffineWitness y ) - where elj ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) - , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) ) - mex [(Local x₁, Shade' y₁ ey₁),(Local x₀, Shade' y₀ ey₀)] - = return $ Shade' (dx-+|>δy) - (Norm . LinearFunction $ \δj -> δx ⊗ (σey<$|δj $ δx)) - where Just δx = x₁.-~.x₀ - δx' = (mex<$|δx) - dx = δx'^/(δx'<.>^δx) - Just δy = y₁.-~.y₀ - σey = convolveMetric ([]::[y]) ey₀ ey₁ - elj _ mex (po:ps) - | DualSpaceWitness <- dualSpaceWitness :: DualNeedleWitness y - , length ps > 1 - = mixShade's =<< (:|) <$> estimateLocalJacobian mex ps - <*> sequenceA [estimateLocalJacobian mex [po,pi] | pi<-ps] - elj _ _ _ = return $ Shade' zeroV mempty - - - -propagateDEqnSolution_loc :: ∀ x y . ( WithField ℝ Manifold x - , Refinable y, Geodesic (Interior y) - , SimpleSpace (Needle x) ) - => DifferentialEqn x y - -> LocalDataPropPlan x (Shade' y) - -> Maybe (Shade' y) -propagateDEqnSolution_loc f propPlan - = pdesl (dualSpaceWitness :: DualNeedleWitness x) - (dualSpaceWitness :: DualNeedleWitness y) - (boundarylessWitness :: BoundarylessWitness x) - (pseudoAffineWitness :: PseudoAffineWitness y) - where pdesl DualSpaceWitness DualSpaceWitness BoundarylessWitness - (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) - | Nothing <- jacobian = Nothing - | otherwise = pure result - where jacobian = f shxy ^. predictDerivatives - Just (Shade' j₀ jExpa) = jacobian - - mx = propPlan^.sourcePosition .+~^ propPlan^.targetPosOffset ^/ 2 - Just my = middleBetween (propPlan^.sourceData.shadeCtr) - (propPlan^.targetAPrioriData.shadeCtr) - shxy = coverAllAround (mx, my) - [ (δx ^-^ propPlan^.targetPosOffset ^/ 2, py ^+^ v) - | (δx,ney) <- (zeroV, propPlan^.sourceData) - : (propPlan^.relatedData) - , let Just py = ney^.shadeCtr .-~. my - , v <- normSpanningSystem' (ney^.shadeNarrowness) - ] - (Shade _ expax' :: Shade x) - = coverAllAround (propPlan^.sourcePosition) - [δx | (δx,_) <- propPlan^.relatedData] - expax = dualNorm expax' - result :: Shade' y - result = convolveShade' - (propPlan^.sourceData) - (Shade' δyb $ applyLinMapNorm jExpa dx) - where δyb = j₀ $ δx - δx = propPlan^.targetPosOffset - dx = δx'^/(δx'<.>^δx) - where δx' = expax<$|δx - -applyLinMapNorm :: ∀ x y . (LSpace x, LSpace y, Scalar x ~ Scalar y) - => Norm (x+>y) -> DualVector x -> Norm y -applyLinMapNorm = case dualSpaceWitness :: DualSpaceWitness y of - DualSpaceWitness -> \n dx -> transformNorm (arr $ LinearFunction (dx-+|>)) n - -ignoreDirectionalDependence :: ∀ x y . (LSpace x, LSpace y, Scalar x ~ Scalar y) - => (x, DualVector x) -> Norm (x+>y) -> Norm (x+>y) -ignoreDirectionalDependence = case dualSpaceWitness :: DualSpaceWitness y of - DualSpaceWitness -> \(v,v') -> transformNorm . arr . LinearFunction $ - \j -> j . arr (LinearFunction $ \x -> x ^-^ v^*(v'<.>^x)) - type Twig x = (Int, ShadeTree x) type TwigEnviron x = [Twig x] @@ -1693,64 +970,6 @@ --- | Essentially the same as @(x,y)@, but not considered as a product topology. --- The 'Semimanifold' etc. instances just copy the topology of @x@, ignoring @y@. -data x`WithAny`y - = WithAny { _untopological :: y - , _topological :: !x } - deriving (Hask.Functor, Show, Generic) - -instance (NFData x, NFData y) => NFData (WithAny x y) - -instance ∀ x y . (Semimanifold x) => Semimanifold (x`WithAny`y) where - type Needle (WithAny x y) = Needle x - type Interior (WithAny x y) = Interior x `WithAny` y - WithAny y x .+~^ δx = WithAny y $ x.+~^δx - fromInterior (WithAny y x) = WithAny y $ fromInterior x - toInterior (WithAny y x) = fmap (WithAny y) $ toInterior x - translateP = tpWD - where tpWD :: ∀ x y . Semimanifold x => Tagged (WithAny x y) - (Interior x`WithAny`y -> Needle x -> Interior x`WithAny`y) - tpWD = Tagged `id` \(WithAny y x) δx -> WithAny y $ tpx x δx - where Tagged tpx = translateP :: Tagged x (Interior x -> Needle x -> Interior x) - semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of - SemimanifoldWitness BoundarylessWitness -> SemimanifoldWitness BoundarylessWitness - -instance (PseudoAffine x) => PseudoAffine (x`WithAny`y) where - WithAny _ x .-~. WithAny _ ξ = x.-~.ξ - pseudoAffineWitness = case pseudoAffineWitness :: PseudoAffineWitness x of - PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) - -> PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) - -instance (AffineSpace x) => AffineSpace (x`WithAny`y) where - type Diff (WithAny x y) = Diff x - WithAny _ x .-. WithAny _ ξ = x.-.ξ - WithAny y x .+^ δx = WithAny y $ x.+^δx - -instance (VectorSpace x, Monoid y) => VectorSpace (x`WithAny`y) where - type Scalar (WithAny x y) = Scalar x - μ *^ WithAny y x = WithAny y $ μ*^x - -instance (AdditiveGroup x, Monoid y) => AdditiveGroup (x`WithAny`y) where - zeroV = WithAny mempty zeroV - negateV (WithAny y x) = WithAny y $ negateV x - WithAny y x ^+^ WithAny υ ξ = WithAny (mappend y υ) (x^+^ξ) - -instance (AdditiveGroup x) => Hask.Applicative (WithAny x) where - pure x = WithAny x zeroV - WithAny f x <*> WithAny t ξ = WithAny (f t) (x^+^ξ) - -instance (AdditiveGroup x) => Hask.Monad (WithAny x) where - return x = WithAny x zeroV - WithAny y x >>= f = WithAny r $ x^+^q - where WithAny r q = f y - -shadeWithAny :: y -> Shade x -> Shade (x`WithAny`y) -shadeWithAny y (Shade x xe) = Shade (WithAny y x) xe - -shadeWithoutAnything :: Shade (x`WithAny`y) -> Shade x -shadeWithoutAnything (Shade (WithAny _ b) e) = Shade b e - constShaded :: y -> ShadeTree x -> x`Shaded`y constShaded y = unsafeFmapTree (WithAny y<$>) id (shadeWithAny y) @@ -1789,7 +1008,7 @@ -- | This is to 'ShadeTree' as 'Data.Map.Map' is to 'Data.Set.Set'. type x`Shaded`y = ShadeTree (x`WithAny`y) -stiWithDensity :: ∀ x y . ( WithField ℝ PseudoAffine x, WithField ℝ LinearManifold y +stiWithDensity :: ∀ x y . ( WithField ℝ PseudoAffine x, LinearSpace y, Scalar y ~ ℝ , SimpleSpace (Needle x) ) => x`Shaded`y -> x -> Cℝay y stiWithDensity (PlainLeaves lvs) @@ -1832,7 +1051,7 @@ -> ( xloc, ( (yloc, recip $ shd|$|(0,1)) , dependence (dualNorm shd) ) ) -smoothInterpolate :: ∀ x y . ( WithField ℝ Manifold x, WithField ℝ LinearManifold y +smoothInterpolate :: ∀ x y . ( WithField ℝ Manifold x, LinearSpace y, Scalar y ~ ℝ , SimpleSpace (Needle x) ) => NonEmpty (x,y) -> x -> y smoothInterpolate = si boundarylessWitness @@ -1911,10 +1130,4 @@ - - -extractJust :: (a->Maybe b) -> [a] -> (Maybe b, [a]) -extractJust f [] = (Nothing,[]) -extractJust f (x:xs) | Just r <- f x = (Just r, xs) - | otherwise = second (x:) $ extractJust f xs diff --git a/Data/Manifold/Types.hs b/Data/Manifold/Types.hs --- a/Data/Manifold/Types.hs +++ b/Data/Manifold/Types.hs @@ -18,6 +18,7 @@ {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE LiberalTypeSynonyms #-} +{-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE TupleSections #-} @@ -58,6 +59,8 @@ , fathomCutDistance, sideOfCut, cutPosBetween -- * Linear mappings , LinearMap, LocalLinear + -- * Misc + , StiefelScalar ) where @@ -81,6 +84,7 @@ import Math.LinearMap.Category import qualified Prelude +import qualified Data.Traversable as Hask import Control.Category.Constrained.Prelude hiding ((^)) import Control.Arrow.Constrained @@ -89,6 +93,8 @@ import Data.Type.Coercion +type StiefelScalar s = (RealFloat s, UArr.Unbox s) + #define deriveAffine(c,t) \ instance (c) => Semimanifold (t) where { \ type Needle (t) = Diff (t); \ @@ -101,6 +107,7 @@ newtype Stiefel1Needle v = Stiefel1Needle { getStiefel1Tangent :: UArr.Vector (Scalar v) } +deriving instance (Eq (Scalar v), UArr.Unbox (Scalar v)) => Eq (Stiefel1Needle v) newtype Stiefel1Basis v = Stiefel1Basis { getStiefel1Basis :: Int } s1bTrie :: ∀ v b. FiniteFreeSpace v => (Stiefel1Basis v->b) -> Stiefel1Basis v:->:b s1bTrie = \f -> St1BTrie $ fmap (f . Stiefel1Basis) allIs @@ -163,11 +170,12 @@ deriveAffine((FiniteFreeSpace v, UArr.Unbox (Scalar v)), Stiefel1Needle v) -instance ∀ v . (LSpace v, FiniteFreeSpace v, UArr.Unbox (Scalar v)) +instance ∀ v . (LSpace v, FiniteFreeSpace v, Eq (Scalar v), UArr.Unbox (Scalar v)) => TensorSpace (Stiefel1Needle v) where type TensorProduct (Stiefel1Needle v) w = Array w scalarSpaceWitness = case scalarSpaceWitness :: ScalarSpaceWitness v of ScalarSpaceWitness -> ScalarSpaceWitness + linearManifoldWitness = LinearManifoldWitness BoundarylessWitness zeroTensor = Tensor $ Arr.replicate (freeDimension ([]::[v]) - 1) zeroV toFlatTensor = LinearFunction $ Tensor . Arr.convert . getStiefel1Tangent fromFlatTensor = LinearFunction $ Stiefel1Needle . Arr.convert . getTensorProduct @@ -184,6 +192,7 @@ fzipTensorWith = bilinearFunction $ \f (Tensor a, Tensor b) -> Tensor $ Arr.zipWith (curry $ arr f) a b coerceFmapTensorProduct _ Coercion = Coercion + wellDefinedTensor (Tensor a) = Tensor <$> Hask.traverse wellDefinedVector a asTensor :: Coercion (LinearMap s a b) (Tensor s (DualVector a) b) asTensor = Coercion @@ -194,7 +203,7 @@ => LinearMap s a b -> a -> b (+$>) = getLinearFunction . getLinearFunction applyLinear -instance ∀ v . (LSpace v, FiniteFreeSpace v, UArr.Unbox (Scalar v)) +instance ∀ v . (LSpace v, FiniteFreeSpace v, Eq (Scalar v), UArr.Unbox (Scalar v)) => LinearSpace (Stiefel1Needle v) where type DualVector (Stiefel1Needle v) = Stiefel1Needle v linearId = LinearMap . Arr.generate d $ \i -> Stiefel1Needle . Arr.generate d $ @@ -228,9 +237,9 @@ composeLinear = bilinearFunction $ \f (LinearMap g) -> LinearMap $ Arr.map (getLinearFunction applyLinear f$) g -instance ∀ k v . - ( WithField k LinearManifold v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v) - , RealFloat k, UArr.Unbox k ) => Semimanifold (Stiefel1 v) where +instance ∀ v . + ( LinearSpace v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v) + , StiefelScalar (Scalar v) ) => Semimanifold (Stiefel1 v) where type Needle (Stiefel1 v) = Stiefel1Needle v fromInterior = id toInterior = pure @@ -261,9 +270,9 @@ insi ti v = Arr.generate d $ \i -> if | i v Arr.! i | i>im -> v Arr.! (i-1) | otherwise -> ti -instance ∀ k v . - ( WithField k LinearManifold v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v) - , RealFloat k, UArr.Unbox k ) => PseudoAffine (Stiefel1 v) where +instance ∀ v . + ( LinearSpace v, FiniteFreeSpace v, FiniteFreeSpace (DualVector v) + , StiefelScalar (Scalar v) ) => PseudoAffine (Stiefel1 v) where (.-~.) = dpst dualSpaceWitness where dpst :: DualSpaceWitness v -> Stiefel1 v -> Stiefel1 v -> Maybe (Stiefel1Needle v) dpst DualSpaceWitness (Stiefel1 s) (Stiefel1 t) diff --git a/Data/Manifold/Web.hs b/Data/Manifold/Web.hs --- a/Data/Manifold/Web.hs +++ b/Data/Manifold/Web.hs @@ -45,22 +45,27 @@ -- ** Local environments , localFocusWeb -- * Uncertain functions - , differentiateUncertainWebFunction + , differentiateUncertainWebFunction, differentiate²UncertainWebFunction -- * Differential equations -- ** Fixed resolution - , filterDEqnSolution_static, iterateFilterDEqn_static + , iterateFilterDEqn_static -- ** Automatic resolution , filterDEqnSolutions_adaptive, iterateFilterDEqn_adaptive -- ** Configuration , InconsistencyStrategy(..) + , InformationMergeStrategy(..) + , naïve, inconsistencyAware, indicateInconsistencies + , PropagationInconsistency(..) -- * Misc - , ConvexSet(..), ellipsoid, coerceWebDomain + , ConvexSet(..), ellipsoid, ellipsoidSet, coerceWebDomain + , rescanPDEOnWeb, rescanPDELocally, webOnions ) where import Data.List hiding (filter, all, foldr1) import Data.Maybe import qualified Data.Set as Set +import qualified Data.Map as Map import qualified Data.Vector as Arr import qualified Data.Vector.Mutable as MArr import qualified Data.Vector.Unboxed as UArr @@ -81,9 +86,14 @@ import Data.Manifold.Types import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine +import Data.Manifold.Shade import Data.Manifold.TreeCover import Data.SetLike.Intersection import Data.Manifold.Riemannian +import Data.Manifold.Atlas +import Data.Manifold.Function.Quadratic +import Data.Function.Affine +import Data.Embedding import qualified Prelude as Hask hiding(foldl, sum, sequence) import qualified Control.Applicative as Hask @@ -92,6 +102,7 @@ import Data.STRef (newSTRef, modifySTRef, readSTRef) import Control.Monad.Trans.State import Control.Monad.Trans.List +import Control.Monad.Trans.Except import Data.Functor.Identity (Identity(..)) import qualified Data.Foldable as Hask import Data.Foldable (all, toList) @@ -107,6 +118,7 @@ import Data.Traversable.Constrained (Traversable, traverse) import Control.Comonad (Comonad(..)) +import Control.Comonad.Cofree import Control.Lens ((&), (%~), (^.), (.~), (+~)) import Control.Lens.TH @@ -145,7 +157,18 @@ } makeLenses ''NeighbourhoodVector +data PropagationInconsistency x υ = PropagationInconsistency { + _inconsistentPropagatedData :: [(x,υ)] + , _inconsistentAPrioriData :: υ } + | PropagationInconsistencies [PropagationInconsistency x υ] + deriving (Show) +makeLenses ''PropagationInconsistency +instance Monoid (PropagationInconsistency x υ) where + mempty = PropagationInconsistencies [] + mappend p q = mconcat [p,q] + mconcat = PropagationInconsistencies + instance (NFData x, NFData (Metric x)) => NFData (Neighbourhood x) -- | A 'PointsWeb' is almost, but not quite a mesh. It is a stongly connected† @@ -532,7 +555,24 @@ ]), n) ) asd' +localOnion :: ∀ x y . WithField ℝ Manifold x + => WebLocally x y -> [[WebLocally x y]] +localOnion origin = go Map.empty $ Map.singleton (origin^.thisNodeId) (1, origin) + where go previous next + | Map.null next = [] + | otherwise = ( snd <$> sortBy (comparing $ negate . fst) + (Hask.toList next) ) + : go (Map.union previous next) + (Map.fromListWith (\(n,ninfo) (n',_) -> (n+n'::Int, ninfo)) + [ (nnid,(1,nneigh)) + | (nid,(_,ninfo))<-Map.toList next + , (nnid,(_,nneigh))<-ninfo^.nodeNeighbours + , Map.notMember nnid previous ]) +webOnions :: ∀ x y . WithField ℝ Manifold x + => PointsWeb x y -> PointsWeb x [[(x,y)]] +webOnions = localFmapWeb (map (map $ _thisNodeCoord&&&_thisNodeData) . localOnion) + nearestNeighbour :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => PointsWeb x y -> x -> Maybe (x,y) nearestNeighbour (PointsWeb rsc asd) x = fmap lkBest $ positionIndex empty rsc x @@ -595,11 +635,62 @@ -> PointsWeb x (Shade' (LocalLinear x y)) differentiateUncertainWebFunction = localFmapWeb differentiateUncertainWebLocally -rescanPDELocally :: ∀ x y . - ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , WithField ℝ Refinable y, SimpleSpace (Needle y) ) - => DifferentialEqn x y -> WebLocally x (Shade' y) - -> Maybe (Shade' y) +differentiate²UncertainWebLocally :: ∀ x y + . ( WithField ℝ Manifold x, FlatSpace (Needle x) + , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y) ) + => WebLocally x (Shade' y) + -> Shade' (Needle x ⊗〃+> Needle y) +differentiate²UncertainWebLocally = d²uwl + ( pseudoAffineWitness :: PseudoAffineWitness x + , pseudoAffineWitness :: PseudoAffineWitness y + , dualSpaceWitness :: DualSpaceWitness (Needle x) + , dualSpaceWitness :: DualSpaceWitness (Needle y) ) + where d²uwl ( PseudoAffineWitness (SemimanifoldWitness _) + , PseudoAffineWitness (SemimanifoldWitness _) + , DualSpaceWitness, DualSpaceWitness ) info + = case estimateLocalHessian $ + (\ngb -> case (ngb^.thisNodeCoord .-~. info^.thisNodeCoord) of + Just δx -> (Local δx :: Local x, ngb^.thisNodeData) ) + <$> info :| envi + of + QuadraticModel _ h -> dualShade $ projectShade + (fromEmbedProject (acoSnd.acoSnd ^/ 2) + (snd.snd ^* 2) ) h + where xVol :: SymmetricTensor ℝ (Needle x) + xVol = squareVs $ fst.snd<$>info^.nodeNeighbours + _:directEnvi:remoteEnvi = localOnion info + envi = directEnvi ++ take (nMinData - length directEnvi) (concat remoteEnvi) + nMinData = 1 + regular_neighboursCount + (subbasisDimension (entireBasis :: SubBasis (Needle x))) + +acoSnd :: ∀ s v y . ( Object (Affine s) y, Object (Affine s) v + , LinearSpace v, Scalar v ~ s ) => Affine s y (v,y) +acoSnd = case ( linearManifoldWitness :: LinearManifoldWitness v + , dualSpaceWitness :: DualSpaceWitness (Needle v) + , dualSpaceWitness :: DualSpaceWitness (Needle y) ) of + (LinearManifoldWitness BoundarylessWitness, DualSpaceWitness, DualSpaceWitness) + -> const zeroV &&& id + +-- | Heuristic formula, matches the number of neighbours each vertex has in a one- +-- and two-dimensional count +regular_neighboursCount :: Int -> Int +regular_neighboursCount d + | d>0 = (regular_neighboursCount (d-1) + 1)*2 + | otherwise = 0 + + +differentiate²UncertainWebFunction :: ∀ x y + . ( WithField ℝ Manifold x, FlatSpace (Needle x) + , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y) ) + => PointsWeb x (Shade' y) + -> PointsWeb x (Shade' (Needle x ⊗〃+> Needle y)) +differentiate²UncertainWebFunction = localFmapWeb differentiate²UncertainWebLocally + +rescanPDELocally :: ∀ x y ð . + ( WithField ℝ Manifold x, FlatSpace (Needle x) + , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y) ) + => DifferentialEqn x ð y -> WebLocally x (Shade' y) + -> (Maybe (Shade' y), Maybe (Shade' ð)) rescanPDELocally = case ( dualSpaceWitness :: DualNeedleWitness x , dualSpaceWitness :: DualNeedleWitness y , boundarylessWitness :: BoundarylessWitness x @@ -614,16 +705,20 @@ , v <- normSpanningSystem' (ngb^.thisNodeData.shadeNarrowness)] of LocalDifferentialEqn _ rescan - -> rescan (differentiateUncertainWebLocally info) - (info^.thisNodeData) + -> rescan (info^.thisNodeData) + (differentiateUncertainWebLocally info) + (differentiate²UncertainWebLocally info) -rescanPDEOnWeb :: ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , WithField ℝ Refinable y, SimpleSpace (Needle y) +rescanPDEOnWeb :: ( WithField ℝ Manifold x, FlatSpace (Needle x) + , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y) , Hask.Applicative m ) - => InconsistencyStrategy m x (Shade' y) - -> DifferentialEqn x y -> PointsWeb x (Shade' y) - -> m (PointsWeb x (Shade' y)) -rescanPDEOnWeb strat = traverseWebWithStrategy strat . rescanPDELocally + => InconsistencyStrategy m x (Shade' y, Shade' ð) + -> DifferentialEqn x ð y -> PointsWeb x (Shade' y) + -> m (PointsWeb x (Shade' y, Shade' ð)) +rescanPDEOnWeb strat deq = traverseWebWithStrategy strat + (fzip . rescanPDELocally deq . fmap fst) + . fmap (\shy -> (shy, error + "No default value for inconsistent PDE-rescanning on web")) toGraph :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => PointsWeb x y -> (Graph, Vertex -> (x, y)) @@ -644,12 +739,15 @@ -- ^ If @p@ is in all intersectors, it must also be in the hull. , convexSetIntersectors :: [Shade' x] } -deriving instance ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , Show (Interior x), Show (Needle' x) ) => Show (ConvexSet x) +deriving instance LtdErrorShow x => Show (ConvexSet x) ellipsoid :: Shade' x -> ConvexSet x ellipsoid s = ConvexSet s [s] +ellipsoidSet :: Embedding (->) (Maybe (Shade' x)) (ConvexSet x) +ellipsoidSet = Embedding (\case {Just s -> ConvexSet s [s]; Nothing -> EmptyConvex}) + (\case {ConvexSet h _ -> Just h; EmptyConvex -> Nothing}) + intersectors :: ConvexSet x -> Maybe (NonEmpty (Shade' x)) intersectors (ConvexSet h []) = pure (h:|[]) intersectors (ConvexSet _ (i:sts)) = pure (i:|sts) @@ -684,7 +782,26 @@ dupHead (x:|xs) = x:|x:xs +newtype InformationMergeStrategy n m y' y = InformationMergeStrategy + { mergeInformation :: y -> n y' -> m y } +naïve :: (NonEmpty y -> y) -> InformationMergeStrategy [] Identity (x,y) y +naïve merge = InformationMergeStrategy (\o n -> Identity . merge $ o :| fmap snd n) + +inconsistencyAware :: (NonEmpty y -> m y) -> InformationMergeStrategy [] m (x,y) y +inconsistencyAware merge = InformationMergeStrategy (\o n -> merge $ o :| fmap snd n) + +indicateInconsistencies :: (NonEmpty υ -> Maybe υ) + -> InformationMergeStrategy [] (Except (PropagationInconsistency x υ)) (x,υ) υ +indicateInconsistencies merge = InformationMergeStrategy + (\o n -> case merge $ o :| fmap snd n of + Just r -> pure r + Nothing -> throwE $ PropagationInconsistency n o ) + +maybeAlt :: Hask.Alternative f => Maybe a -> f a +maybeAlt (Just x) = pure x +maybeAlt Nothing = Hask.empty + data InconsistencyStrategy m x y where AbortOnInconsistency :: InconsistencyStrategy Maybe x y IgnoreInconsistencies :: InconsistencyStrategy Identity x y @@ -692,76 +809,56 @@ deriving instance Hask.Functor (InconsistencyStrategy m x) -iterateFilterDEqn_static :: ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , Refinable y, Geodesic (Interior y) - , Hask.Applicative m ) - => InconsistencyStrategy m x (Shade' y) -> DifferentialEqn x y - -> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)] -iterateFilterDEqn_static strategy f - = map (fmap convexSetHull) - . itWhileJust strategy - (filterDEqnSolutions_static (ellipsoid<$>strategy) f) - . fmap (`ConvexSet`[]) +iterateFilterDEqn_static :: ( WithField ℝ Manifold x, FlatSpace (Needle x) + , Refinable y, Geodesic y, FlatSpace (Needle y) + , WithField ℝ AffineManifold ð, Geodesic ð + , SimpleSpace (Needle ð) + , Hask.MonadPlus m ) + => InformationMergeStrategy [] m (x,Shade' y) iy + -> Embedding (->) (Shade' y) iy + -> DifferentialEqn x ð y + -> PointsWeb x (Shade' y) -> Cofree m (PointsWeb x (Shade' y)) +iterateFilterDEqn_static strategy shading f + = fmap (fmap (shading >-$)) + . coiter (filterDEqnSolutions_static strategy shading f) + . fmap (shading $->) -filterDEqnSolution_static :: ∀ x y m . ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , Refinable y, Geodesic (Interior y) ) - => InconsistencyStrategy m x (Shade' y) -> DifferentialEqn x y - -> PointsWeb x (Shade' y) -> m (PointsWeb x (Shade' y)) -filterDEqnSolution_static strat@AbortOnInconsistency f - = case boundarylessWitness :: BoundarylessWitness x of - BoundarylessWitness -> - rescanPDEOnWeb strat f >=> webLocalInfo - >>> Hask.traverse `id`\me -> case me^.nodeNeighbours of - [] -> return $ me^.thisNodeData - ngbs -> refineShade' (me^.thisNodeData) - =<< intersectShade's - =<< ( sequenceA $ NE.fromList - [ propagateDEqnSolution_loc - f (LocalDataPropPlan - (ngbInfo^.thisNodeCoord) - (negateV δx) - (ngbInfo^.thisNodeData) - (me^.thisNodeData) - (fmap (second _thisNodeData . snd) - $ ngbInfo^.nodeNeighbours) - ) - | (_, (δx, ngbInfo)) <- ngbs - ] ) -filterDEqnSolutions_static :: ∀ x y m . - ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , Refinable y, Geodesic (Interior y) - , Hask.Applicative m ) - => InconsistencyStrategy m x (ConvexSet y) -> DifferentialEqn x y - -> PointsWeb x (ConvexSet y) -> m (PointsWeb x (ConvexSet y)) -filterDEqnSolutions_static strategy f - = webLocalInfo - >>> fmap (id &&& rescanPDELocally f . fmap convexSetHull) - >>> localFocusWeb >>> Hask.traverse `id`\((_,(me,updShy)), ngbs) - -> let oldValue = me^.thisNodeData :: ConvexSet y +filterDEqnSolutions_static :: ∀ x y iy ð m . + ( WithField ℝ Manifold x, FlatSpace (Needle x) + , Refinable y, Geodesic y, FlatSpace (Needle y) + , WithField ℝ AffineManifold ð, Geodesic ð + , SimpleSpace (Needle ð) + , Hask.MonadPlus m ) + => InformationMergeStrategy [] m (x,Shade' y) iy -> Embedding (->) (Shade' y) iy + -> DifferentialEqn x ð y -> PointsWeb x iy -> m (PointsWeb x iy) +filterDEqnSolutions_static = case geodesicWitness :: GeodesicWitness y of + GeodesicWitness _ -> \strategy shading f + -> webLocalInfo + >>> fmap (id &&& rescanPDELocally f . fmap (shading>-$)) + >>> localFocusWeb >>> Hask.traverse ( \((_,(me,updShy)), ngbs) + -> let oldValue = me^.thisNodeData :: iy in case updShy of - Just shy -> case ngbs of + (Just shy, Just shð) -> case ngbs of [] -> pure oldValue _:_ | BoundarylessWitness <- (boundarylessWitness::BoundarylessWitness x) - -> handleInconsistency strategy oldValue - $ ( sequenceA $ NE.fromList - [ sj >>= \ngbShy -> + -> maybeAlt + ( sequenceA [ fzip sj + >>= \ngbShyð -> (ngbInfo^.thisNodeCoord,)<$> propagateDEqnSolution_loc f (LocalDataPropPlan (ngbInfo^.thisNodeCoord) (negateV δx) - ngbShy - shy - (fmap (second (convexSetHull . _thisNodeData) + ngbShyð + (shy, shð) + (fmap (second ((shading>-$) . _thisNodeData) . snd) $ ngbInfo^.nodeNeighbours) ) | (δx, (ngbInfo,sj)) <- ngbs ] ) - >>= intersectShade's - >>= pure . ((oldValue<>) . ellipsoid) - >>= \case EmptyConvex -> empty - c -> pure c - _ -> handleInconsistency strategy oldValue empty + >>= mergeInformation strategy (shading$->shy) + _ -> mergeInformation strategy oldValue empty + ) handleInconsistency :: InconsistencyStrategy m x a -> a -> Maybe a -> m a handleInconsistency AbortOnInconsistency _ i = i @@ -790,23 +887,25 @@ oldAndNew' (_, l) = (False,) <$> l -filterDEqnSolutions_adaptive :: ∀ x y badness m +filterDEqnSolutions_adaptive :: ∀ x y ð badness m . ( WithField ℝ Manifold x, SimpleSpace (Needle x) , WithField ℝ AffineManifold y, Refinable y, Geodesic y + , WithField ℝ AffineManifold ð, Geodesic ð, SimpleSpace (Needle ð) , badness ~ ℝ, Hask.Monad m ) => MetricChoice x -- ^ Scalar product on the domain, for regularising the web. -> InconsistencyStrategy m x (Shade' y) - -> DifferentialEqn x y + -> DifferentialEqn x ð y -> (x -> Shade' y -> badness) -> PointsWeb x (SolverNodeState x y) -> m (PointsWeb x (SolverNodeState x y)) filterDEqnSolutions_adaptive mf strategy f badness' oldState - = fmap recomputeJacobian $ filterGo boundarylessWitness - =<< tryPreproc boundarylessWitness - where tryPreproc :: BoundarylessWitness x + = fmap recomputeJacobian $ filterGo boundarylessWitness geodesicWitness + =<< tryPreproc boundarylessWitness geodesicWitness + where tryPreproc :: BoundarylessWitness x -> GeodesicWitness y -> m (PointsWeb x ( (WebLocally x (SolverNodeState x y) , [(Shade' y, badness)]) )) - tryPreproc BoundarylessWitness = traverse addPropagation $ webLocalInfo oldState + tryPreproc BoundarylessWitness (GeodesicWitness _) + = traverse addPropagation $ webLocalInfo oldState where addPropagation wl | null neighbourInfo = pure (wl, []) | otherwise = (wl,) . map (id&&&badness undefined) @@ -818,8 +917,9 @@ (neigh^.thisNodeCoord) (negateV δx) (convexSetHull $ neigh^.thisNodeData - .solverNodeStatus) - (thisShy) + .solverNodeStatus + , undefined) + (thisShy, undefined) [ second (convexSetHull . _solverNodeStatus . _thisNodeData) nn | (_,nn)<-neigh^.nodeNeighbours ] ) @@ -832,11 +932,12 @@ errTgtModulation = (1-) . (`mod'`1) . negate . sqrt $ fromIntegral totalAge badness x = badness' x . (shadeNarrowness %~ (scaleNorm errTgtModulation)) - filterGo :: BoundarylessWitness x + filterGo :: BoundarylessWitness x -> GeodesicWitness y -> (PointsWeb x ( (WebLocally x (SolverNodeState x y) , [(Shade' y, badness)]) )) -> m (PointsWeb x (SolverNodeState x y)) - filterGo BoundarylessWitness preproc'd = fmap (smoothenWebTopology mf + filterGo BoundarylessWitness (GeodesicWitness _) preproc'd + = fmap (smoothenWebTopology mf . fromTopWebNodes mf . concat . fmap retraceBonds . Hask.toList . webLocalInfo . webLocalInfo) $ Hask.traverse (uncurry localChange) preproc'd @@ -912,8 +1013,9 @@ (n^.thisNodeCoord) (stepV ^-^ δx) (convexSetHull $ - n^.thisNodeData.solverNodeStatus) - (aprioriInterpolate) + n^.thisNodeData.solverNodeStatus + , undefined) + (aprioriInterpolate, undefined) (second (convexSetHull ._solverNodeStatus ._thisNodeData) @@ -978,10 +1080,11 @@ iterateFilterDEqn_adaptive :: ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , WithField ℝ AffineManifold y, Refinable y, Geodesic y, Hask.Monad m ) + , WithField ℝ AffineManifold y, Refinable y, Geodesic y, Hask.Monad m + , WithField ℝ AffineManifold ð, Geodesic ð, SimpleSpace (Needle ð) ) => MetricChoice x -- ^ Scalar product on the domain, for regularising the web. -> InconsistencyStrategy m x (Shade' y) - -> DifferentialEqn x y + -> DifferentialEqn x ð y -> (x -> Shade' y -> ℝ) -- ^ Badness function for local results. -> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)] iterateFilterDEqn_adaptive mf strategy f badness diff --git a/images/examples/TreesAndWebs/2D-cartesian_onion.png b/images/examples/TreesAndWebs/2D-cartesian_onion.png new file mode 100644 Binary files /dev/null and b/images/examples/TreesAndWebs/2D-cartesian_onion.png differ diff --git a/manifolds.cabal b/manifolds.cabal --- a/manifolds.cabal +++ b/manifolds.cabal @@ -1,5 +1,5 @@ Name: manifolds -Version: 0.4.0.0 +Version: 0.4.1.0 Category: Math Synopsis: Coordinate-free hypersurfaces Description: Manifolds, a generalisation of the notion of “smooth curves” or surfaces, @@ -40,18 +40,20 @@ Library Build-Depends: base>=4.5 && < 6 - , manifolds-core == 0.4.0.0 + , manifolds-core == 0.4.1.0 , transformers , vector-space>=0.8 , free-vector-spaces>=0.1.1 , linear , MemoTrie , vector - , linearmap-category > 0.3 && < 0.4 + , linearmap-category >= 0.3.2 && < 0.4 , containers , comonad + , free , semigroups , void + , number-show >= 0.1 && < 0.2 , tagged , deepseq , lens @@ -70,10 +72,12 @@ Exposed-modules: Data.Manifold Data.Manifold.PseudoAffine Data.Manifold.TreeCover + Data.Manifold.Shade Data.Manifold.Web Data.Manifold.DifferentialEquation Data.SimplicialComplex Data.Function.Differentiable + Data.Function.Affine Data.Manifold.Types Data.Manifold.Types.Stiefel Data.Manifold.Griddable @@ -85,8 +89,8 @@ Data.Manifold.Cone Data.CoNat Data.Embedding + Data.Manifold.Function.Quadratic Data.Function.Differentiable.Data - Data.Function.Affine Control.Monad.Trans.OuterMaybe Util.Associate Util.LtdShow