diff --git a/Data/Function/Differentiable.hs b/Data/Function/Differentiable.hs --- a/Data/Function/Differentiable.hs +++ b/Data/Function/Differentiable.hs @@ -676,7 +676,7 @@ where npr (LinearManifoldWitness BoundarylessWitness) (ClosedScalarWitness :: ClosedScalarWitness s) = PreRegion $ ppr . ngt - where PreRegion ppr = positivePreRegion' + where PreRegion ppr = positivePreRegion' :: PreRegion s s ngt = actuallyLinearEndo $ negateV id preRegionToInfFrom, preRegionFromMinInfTo :: RealDimension s => s -> PreRegion s s @@ -689,14 +689,14 @@ where prif (LinearManifoldWitness BoundarylessWitness) (ClosedScalarWitness :: ClosedScalarWitness s) xs = PreRegion $ ppr . trl - where PreRegion ppr = positivePreRegion' + where PreRegion ppr = positivePreRegion' :: PreRegion s s trl = actuallyAffineEndo (-xs) id preRegionFromMinInfTo' = prif (linearManifoldWitness :: LinearManifoldWitness s) (closedScalarWitness :: ClosedScalarWitness s) where prif (LinearManifoldWitness BoundarylessWitness) (ClosedScalarWitness :: ClosedScalarWitness s) xe = PreRegion $ ppr . flp - where PreRegion ppr = positivePreRegion' + where PreRegion ppr = positivePreRegion' :: PreRegion s s flp = actuallyAffineEndo xe (negateV id) intervalPreRegion :: ∀ s . RealDimension s => (s,s) -> PreRegion s s diff --git a/Data/Manifold/DifferentialEquation.hs b/Data/Manifold/DifferentialEquation.hs --- a/Data/Manifold/DifferentialEquation.hs +++ b/Data/Manifold/DifferentialEquation.hs @@ -36,7 +36,6 @@ DifferentialEqn, ODE , constLinearDEqn , constLinearODE - , constLinearPDE , iterateFilterDEqn_static -- * Cost functions for error bounds , maxDeviationsGoal @@ -59,6 +58,8 @@ import Data.Manifold.Types import Data.Manifold.PseudoAffine +import Data.Manifold.Shade +import Data.Manifold.Function.LocalModel import Data.Function.Differentiable import Data.Function.Differentiable.Data import Data.Manifold.TreeCover @@ -82,6 +83,7 @@ import Data.Foldable.Constrained import Data.Traversable.Constrained (Traversable, traverse) +import Control.Lens -- | An ordinary differential equation is one that does not need any a-priori -- partial derivatives to compute the derivative for integration in some @@ -90,44 +92,37 @@ -- be an arbitrary one-dimensional space (i.e. basically real intervals or 'S¹'). -- In these cases, there is always only one partial derivative: that which we -- integrate over, in the only possible direction for propagation. -type ODE x y = DifferentialEqn x ℝ⁰ y +type ODE x y = DifferentialEqn AffineModel x y -constLinearDEqn :: ∀ x y ð . ( SimpleSpace x - , SimpleSpace y, AffineManifold y - , SimpleSpace ð, AffineManifold ð - , Scalar x ~ ℝ, Scalar y ~ ℝ, Scalar ð ~ ℝ ) - => ((y,ð) +> (x +> y)) -> ((x +> y) +> (y,ð)) -> DifferentialEqn x ð y +constLinearDEqn :: ∀ x y . ( SimpleSpace x + , SimpleSpace y, AffineManifold y + , Scalar x ~ ℝ, Scalar y ~ ℝ ) + => (y +> (x +> y)) -> ((x +> y) +> y) -> DifferentialEqn QuadraticModel x y constLinearDEqn = case ( linearManifoldWitness :: LinearManifoldWitness x , dualSpaceWitness :: DualSpaceWitness x , linearManifoldWitness :: LinearManifoldWitness y - , dualSpaceWitness :: DualSpaceWitness y - , linearManifoldWitness :: LinearManifoldWitness ð - , dualSpaceWitness :: DualSpaceWitness ð ) of + , dualSpaceWitness :: DualSpaceWitness y ) of ( LinearManifoldWitness BoundarylessWitness, DualSpaceWitness - ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ) -> \bwt'inv bwt' -> \(Shade (_x,y) δxy) -> LocalDifferentialEqn - { _predictDerivatives - = \(Shade' ð δð) -> - let j = bwt'inv $ (y,ð) - δj = bwt' `transformNorm` - sumSubspaceNorms (transformNorm (zeroV&&&id) $ dualNorm δxy) δð - in return $ Shade' j δj - , _rescanDerivatives - = \shy shjApriori _ - -> ( mixShade's $ shy + { _rescanDifferentialEqn + = \(QuadraticModel shy' shj'Apriori _) -> + let shy = dualShade shy' + shjApriori = dualShade shj'Apriori + in ( mixShade's $ shy :| [ projectShade - (Embedding (arr bwt'inv <<< id&&&zeroV) - (arr bwt' >>> fst)) + (Embedding (arr bwt'inv) + (arr bwt')) shjApriori ] , return $ projectShade - (Embedding (arr bwt'inv <<< zeroV&&&id) - (arr bwt' >>> snd)) - shjApriori + (Embedding (arr bwt') + (arr bwt'inv)) + shy ) } -constLinearODE :: ∀ x y . ( SimpleSpace x, Scalar x ~ ℝ, SimpleSpace y, Scalar y ~ ℝ ) +constLinearODE :: ∀ x y . ( SimpleSpace x, Scalar x ~ ℝ + , AffineManifold y, SimpleSpace y, Scalar y ~ ℝ ) => ((x +> y) +> y) -> ODE x y constLinearODE = case ( linearManifoldWitness :: LinearManifoldWitness x , dualSpaceWitness :: DualSpaceWitness x @@ -135,40 +130,15 @@ , dualSpaceWitness :: DualSpaceWitness y ) of ( LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ) -> \bwt' -> - let bwt'inv = (bwt'\$) + let bwt'inv = pseudoInverse bwt' in \(Shade (_x,y) δxy) -> LocalDifferentialEqn - (let j = bwt'inv y - δj = (bwt'>>>zeroV&&&id) `transformNorm` dualNorm δxy - in \_ -> return $ Shade' j δj ) - (\shy _ _ -> (pure shy, Just $ Shade' Origin mempty) ) + (\(AffineModel shy' _) -> + let shy = dualShade shy' + in ( return $ shy & shadeNarrowness %~ scaleNorm 0.01 + , return $ projectShade (Embedding (arr bwt') + (arr bwt'inv)) shy ) + ) -constLinearPDE :: ∀ x y ð . - ( WithField ℝ SimpleSpace x - , WithField ℝ SimpleSpace y - , WithField ℝ SimpleSpace ð, AffineManifold ð ) - => ((x +> y) +> ð) -> (ð +> (x +> y)) -> DifferentialEqn x ð y -constLinearPDE = case ( linearManifoldWitness :: LinearManifoldWitness x - , dualSpaceWitness :: DualSpaceWitness x - , linearManifoldWitness :: LinearManifoldWitness y - , dualSpaceWitness :: DualSpaceWitness y - , linearManifoldWitness :: LinearManifoldWitness ð - , dualSpaceWitness :: DualSpaceWitness ð ) of - ( LinearManifoldWitness BoundarylessWitness, DualSpaceWitness - ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness - ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ) - -> \bwt' bwt'inv (Shade (_x,y) δxy) - -> LocalDifferentialEqn - { _predictDerivatives - = \(Shade' ð δð) -> - let j = bwt'inv $ ð - δj = bwt' `transformNorm` δð - in return $ Shade' j δj - , _rescanDerivatives - = \shy shjApriori _ - -> ( return shy - , return $ projectShade (Embedding (arr bwt'inv) (arr bwt')) shjApriori - ) - } -- | A function that variates, relatively speaking, most strongly -- for arguments around 1. In the zero-limit it approaches a constant diff --git a/Data/Manifold/Function/LocalModel.hs b/Data/Manifold/Function/LocalModel.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/Function/LocalModel.hs @@ -0,0 +1,310 @@ +-- | +-- Module : Data.Manifold.Function.LocalModel +-- Copyright : (c) Justus Sagemüller 2017 +-- License : GPL v3 +-- +-- Maintainer : (@) jsagemue $ uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- + +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UnicodeSyntax #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE ConstraintKinds #-} + +module Data.Manifold.Function.LocalModel ( + -- * The model class + LocalModel (..), ModellableRelation + -- ** Local data fit models + , AffineModel(..), QuadraticModel(..) + , estimateLocalJacobian, estimateLocalHessian + , propagationCenteredModel + , propagationCenteredQuadraticModel + , quadraticModel_derivatives + -- ** Differential equations + , DifferentialEqn, LocalDifferentialEqn(..) + , propagateDEqnSolution_loc, LocalDataPropPlan(..) + -- ** Range interpolation + , rangeWithinVertices + ) where + + +import Data.Manifold.Types +import Data.Manifold.PseudoAffine +import Data.Manifold.Types.Primitive ((^)) +import Data.Manifold.Shade +import Data.Manifold.Riemannian + +import Data.VectorSpace +import Math.LinearMap.Category + +import Data.List.NonEmpty (NonEmpty (..)) +import qualified Data.List.NonEmpty as NE + +import qualified Prelude as Hask + +import Control.Category.Constrained.Prelude +import Control.Arrow.Constrained + +import Control.Lens +import Control.Lens.TH + + +newtype LocalDifferentialEqn ㄇ x y = LocalDifferentialEqn { + _rescanDifferentialEqn :: ㄇ x y + -> (Maybe (Shade' y), Maybe (Shade' (LocalLinear x y))) + } +makeLenses ''LocalDifferentialEqn + +type DifferentialEqn ㄇ x y = Shade (x,y) -> LocalDifferentialEqn ㄇ x y + +data LocalDataPropPlan x y = LocalDataPropPlan + { _sourcePosition :: !(Interior x) + , _targetPosOffset :: !(Needle x) + , _sourceData, _targetAPrioriData :: !y + , _relatedData :: [(Needle x, y)] + } +deriving instance (Show (Interior x), Show y, Show (Needle x)) + => Show (LocalDataPropPlan x y) + +makeLenses ''LocalDataPropPlan + + +{-# DEPRECATED estimateLocalJacobian "Use `fitLocally`" #-} +estimateLocalJacobian :: ∀ x y . ( WithField ℝ Manifold x, Refinable y + , SimpleSpace (Needle x), SimpleSpace (Needle y) ) + => Metric x -> [(Local x, Shade' y)] + -> Maybe (Shade' (LocalLinear x y)) +estimateLocalJacobian = elj ( pseudoAffineWitness :: PseudoAffineWitness x + , pseudoAffineWitness :: PseudoAffineWitness y ) + where elj ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) + , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) ) + mex [(Local x₁, Shade' y₁ ey₁),(Local x₀, Shade' y₀ ey₀)] + = return $ Shade' (dx-+|>δy) + (Norm . LinearFunction $ \δj -> δx ⊗ (σey<$|δj $ δx)) + where Just δx = x₁.-~.x₀ + δx' = (mex<$|δx) + dx = δx'^/(δx'<.>^δx) + Just δy = y₁.-~.y₀ + σey = convolveMetric ([]::[y]) ey₀ ey₁ + elj _ mex (po:ps) + | DualSpaceWitness <- dualSpaceWitness :: DualNeedleWitness y + , length ps > 1 + = mixShade's =<< (:|) <$> estimateLocalJacobian mex ps + <*> sequenceA [estimateLocalJacobian mex [po,pi] | pi<-ps] + elj _ _ _ = return $ Shade' zeroV mempty + + +data AffineModel x y = AffineModel { + _affineModelOffset :: Shade y + , _affineModelLCoeff :: Shade ( Needle x +>Needle y) + } +deriving instance (Show (Shade y), Show (Shade (Needle x+>Needle y))) + => Show (AffineModel x y) +makeLenses ''AffineModel + + +data QuadraticModel x y = QuadraticModel { + _quadraticModelOffset :: Shade y + , _quadraticModelLCoeff :: Shade ( Needle x +>Needle y) + , _quadraticModelQCoeff :: Shade (Needle x⊗〃+>Needle y) + } +deriving instance ( Show (Shade y) + , Show (Shade (Needle x+>Needle y)) + , Show (Shade (Needle x⊗〃+>Needle y)) ) + => Show (QuadraticModel x y) +makeLenses ''QuadraticModel + +type QModelTup s x y = ( Needle y, (Needle x+>Needle y + , SymmetricTensor s (Needle x)+>(Needle y)) ) + + + +quadratic_linearRegression :: ∀ s x y . + ( WithField s PseudoAffine x + , WithField s PseudoAffine y, Geodesic y + , SimpleSpace (Needle x), SimpleSpace (Needle y) ) + => NE.NonEmpty (Needle x, Shade' y) -> QuadraticModel x y +quadratic_linearRegression = case ( dualSpaceWitness :: DualSpaceWitness (Needle x) + , dualSpaceWitness :: DualSpaceWitness (Needle y) ) of + (DualSpaceWitness, DualSpaceWitness) -> gLinearRegression + (\δx -> lfun $ \(c,(b,a)) -> (a $ squareV δx) ^+^ (b $ δx) ^+^ c ) + (\cmy (cBest, (bBest, aBest)) σ + -> let (σc, (σb, σa)) = second summandSpaceNorms $ summandSpaceNorms σ + in QuadraticModel (Shade (cmy⊙+^cBest $ ([]::[y])) σc) + (Shade bBest σb) + (Shade aBest σa) ) + +gLinearRegression :: ∀ s x y ㄇ ψ. + ( WithField s PseudoAffine x + , WithField s PseudoAffine y, Geodesic y + , SimpleSpace (Needle x), SimpleSpace (Needle y) + , SimpleSpace ψ, Scalar ψ ~ s ) + => (Needle x -> ψ -+> Needle y) + -> (Interior y -> ψ -> Variance ψ -> ㄇ x y) + -> NE.NonEmpty (Needle x, Shade' y) -> ㄇ x y +gLinearRegression fwdCalc analyse = qlr (pseudoAffineWitness, geodesicWitness) + where qlr :: (PseudoAffineWitness y, GeodesicWitness y) + -> NE.NonEmpty (Needle x, Shade' y) -> ㄇ x y + qlr (PseudoAffineWitness (SemimanifoldWitness _), GeodesicWitness _) ps + = analyse cmy ψ σψ + where Just cmy = pointsBarycenter $ _shade'Ctr.snd<$>ps + Just vsxy = Hask.mapM (\(x, Shade' y ey) -> (x,).(,ey)<$>y.-~.cmy) ps + ψ = linearFit_bestModel regResult + σψ = dualNorm . (case linearFit_χν² regResult of + χν² | χν² > 0, recip χν² > 0 + -> scaleNorm (recip $ 1 + sqrt χν²) + _ -> {-Dbg.trace ("Fit for regression model requires" + ++" well-defined χν² (which needs positive number of degrees of freedom)." + ++"\n Data: "++show (length ps + * subbasisDimension (entireBasis :: SubBasis (Needle y))) + ++"\n Model parameters: "++show (subbasisDimension + (entireBasis :: SubBasis ψ)) )-} + id) + $ linearFit_modelUncertainty regResult + regResult = linearRegression (arr . fwdCalc) (NE.toList vsxy) + +quadraticModel_derivatives :: ∀ x y . + ( PseudoAffine x, PseudoAffine y + , SimpleSpace (Needle x), SimpleSpace (Needle y) + , Scalar (Needle y) ~ Scalar (Needle x) ) => + QuadraticModel x y -> (Shade' y, (Shade' (LocalLinear x y), Shade' (LocalBilinear x y))) +quadraticModel_derivatives (QuadraticModel sh shð shð²) + | (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) + :: PseudoAffineWitness y <- pseudoAffineWitness + , DualSpaceWitness :: DualSpaceWitness (Needle x) <- dualSpaceWitness + , DualSpaceWitness :: DualSpaceWitness (Needle y) <- dualSpaceWitness + = (dualShade sh, ( dualShade shð + , linIsoTransformShade (2*^id) $ dualShade shð² )) + +{-# DEPRECATED estimateLocalHessian "Use `fitLocally`" #-} +estimateLocalHessian :: ∀ x y . ( WithField ℝ Manifold x, Refinable y, Geodesic y + , FlatSpace (Needle x), FlatSpace (Needle y) ) + => NonEmpty (Local x, Shade' y) -> QuadraticModel x y +estimateLocalHessian pts = quadratic_linearRegression $ first getLocalOffset <$> pts + + +propagationCenteredModel :: ∀ x y ㄇ . + ( ModellableRelation x y, LocalModel ㄇ ) + => LocalDataPropPlan x (Shade' y) -> ㄇ x y +propagationCenteredModel propPlan = case fitLocally (NE.toList ptsFromCenter) of + Just ㄇ->ㄇ + where ctrOffset = propPlan^.targetPosOffset^/2 + ptsFromCenter = (negateV ctrOffset, propPlan^.sourceData) + :| [(δx^-^ctrOffset, shy) + | (δx, shy) + <- (propPlan^.targetPosOffset, propPlan^.targetAPrioriData) + : propPlan^.relatedData + ] + + +propagationCenteredQuadraticModel :: ∀ x y . + ( ModellableRelation x y ) + => LocalDataPropPlan x (Shade' y) -> QuadraticModel x y +propagationCenteredQuadraticModel = propagationCenteredModel + + +propagateDEqnSolution_loc :: ∀ x y ㄇ . (ModellableRelation x y, LocalModel ㄇ) + => DifferentialEqn ㄇ x y + -> LocalDataPropPlan x (Shade' y) + -> Maybe (Shade' y) +propagateDEqnSolution_loc f propPlan + = pdesl (dualSpaceWitness :: DualNeedleWitness x) + (dualSpaceWitness :: DualNeedleWitness y) + (boundarylessWitness :: BoundarylessWitness x) + (pseudoAffineWitness :: PseudoAffineWitness y) + (geodesicWitness :: GeodesicWitness y) + where pdesl DualSpaceWitness DualSpaceWitness BoundarylessWitness + (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) + (GeodesicWitness _) + | Nothing <- jacobian = Nothing + | otherwise = pure result + where (_,jacobian) = f shxy ^. rescanDifferentialEqn + $ propagationCenteredModel propPlan + jacobianSh :: Shade (LocalLinear x y) + Just jacobianSh = dualShade' <$> jacobian + mx = propPlan^.sourcePosition .+~^ propPlan^.targetPosOffset ^/ 2 :: x + (Shade _ expax' :: Shade x) + = coverAllAround (propPlan^.sourcePosition) + [δx | (δx,_) <- propPlan^.relatedData] + shxy = coverAllAround (mx, mυ) + [ (δx ^-^ propPlan^.targetPosOffset ^/ 2, pυ ^+^ v) + | (δx,neυ) <- (zeroV, propPlan^.sourceData) + : (second id + <$> propPlan^.relatedData) + , let Just pυ = neυ^.shadeCtr .-~. mυ + , v <- normSpanningSystem' (neυ^.shadeNarrowness) + ] + where Just mυ = middleBetween (propPlan^.sourceData.shadeCtr) + (propPlan^.targetAPrioriData.shadeCtr) + expax = dualNorm expax' + result :: Shade' y + result = convolveShade' (propPlan^.sourceData) + (dualShade . linearProjectShade (lfun ($ δx)) $ jacobianSh) + δx = propPlan^.targetPosOffset + + +type ModellableRelation x y = ( WithField ℝ Manifold x + , Refinable y, Geodesic y + , FlatSpace (Needle x), FlatSpace (Needle y) ) + +class LocalModel ㄇ where + fitLocally :: ModellableRelation x y + => [(Needle x, Shade' y)] -> Maybe (ㄇ x y) + tweakLocalOffset :: ModellableRelation x y + => Lens' (ㄇ x y) (Shade y) + +modelParametersOverdetMargin :: Int -> Int +modelParametersOverdetMargin n = n + round (sqrt $ fromIntegral n) - 1 + + +-- | Dimension of the space of affine functions on @v@. +p¹Dimension :: ∀ v p . FiniteDimensional v => p v -> Int +p¹Dimension _ = 1 + d + where d = subbasisDimension (entireBasis :: SubBasis v) + +-- | Dimension of the space of quadratic functions on @v@. +p²Dimension :: ∀ v p . FiniteDimensional v => p v -> Int +p²Dimension _ = 1 + d + (d*(d+1))`div`2 + where d = subbasisDimension (entireBasis :: SubBasis v) + +instance LocalModel AffineModel where + fitLocally = aFitL dualSpaceWitness + where aFitL :: ∀ x y . ModellableRelation x y + => DualSpaceWitness (Needle y) + -> [(Needle x, Shade' y)] -> Maybe (AffineModel x y) + aFitL DualSpaceWitness dataPts + | (p₀:ps, pω:_) <- splitAt (modelParametersOverdetMargin + $ p¹Dimension ([]::[Needle x])) dataPts + = Just . gLinearRegression + (\δx -> lfun $ \(b,a) -> (a $ δx) ^+^ b ) + (\cmy (bBest, aBest) σ + -> let (σb, σa) = summandSpaceNorms σ + in AffineModel (Shade (cmy⊙+^bBest $ ([]::[y])) + $ scaleNorm 2 σb) + -- The magic factor 2 seems dubious ↗, but testing indicates + -- that this is necessary to not overrate the accuracy. + -- TODO: check the algorithms in linearmap-category. + (Shade aBest σa) ) + $ (p₀:|ps++[pω]) + | otherwise = Nothing + tweakLocalOffset = affineModelOffset + +instance LocalModel QuadraticModel where + fitLocally = qFitL + where qFitL :: ∀ x y . ModellableRelation x y + => [(Needle x, Shade' y)] -> Maybe (QuadraticModel x y) + qFitL dataPts + | (p₀:ps, pω:_) <- splitAt (modelParametersOverdetMargin + $ p²Dimension ([]::[Needle x])) dataPts + = Just . quadratic_linearRegression + $ (p₀:|ps++[pω]) + | otherwise = Nothing + tweakLocalOffset = quadraticModelOffset diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs --- a/Data/Manifold/PseudoAffine.hs +++ b/Data/Manifold/PseudoAffine.hs @@ -50,12 +50,12 @@ module Data.Manifold.PseudoAffine ( -- * Manifold class - Manifold + Manifold(inInterior) , Semimanifold(..), Needle' , PseudoAffine(..) -- * Type definitions -- ** Needles - , Local(..) + , Local(..), (⊙+^) -- ** Metrics , Metric, Metric' , RieMetric, RieMetric' @@ -97,7 +97,7 @@ import Data.CoNat -import qualified Prelude +import qualified Prelude as Hask import qualified Control.Applicative as Hask import Control.Category.Constrained.Prelude hiding ((^)) @@ -118,6 +118,9 @@ boundarylessWitness :: BoundarylessWitness m default boundarylessWitness :: (m ~ Interior m) => BoundarylessWitness m boundarylessWitness = BoundarylessWitness + inInterior :: m -> Interior m + default inInterior :: (m ~ Interior m) => m -> Interior m + inInterior = id instance (PseudoAffine m, LSpace (Needle m), Interior m ~ m) => Manifold m @@ -135,8 +138,22 @@ -- ≡ locallyTrivialDiffeomorphism p .+~^ 'coerceNeedle' v -- @ locallyTrivialDiffeomorphism :: x -> ξ - coerceNeedle :: Functor p (->) (->) => p (x,ξ) -> (Needle x -+> Needle ξ) - coerceNeedle' :: Functor p (->) (->) => p (x,ξ) -> (Needle' x -+> Needle' ξ) + coerceNeedle :: Hask.Functor p => p (x,ξ) -> (Needle x -+> Needle ξ) + coerceNeedle' :: Hask.Functor p => p (x,ξ) -> (Needle' x -+> Needle' ξ) + coerceNorm :: Hask.Functor p => p (x,ξ) -> Metric x -> Metric ξ + coerceNorm p = case ( oppositeLocalCoercion :: CanonicalDiffeomorphism ξ x + , dualSpaceWitness :: DualSpaceWitness (Needle x) + , dualSpaceWitness :: DualSpaceWitness (Needle ξ) ) of + (CanonicalDiffeomorphism, DualSpaceWitness, DualSpaceWitness) + -> case ( coerceNeedle (swap<$>p), coerceNeedle' p ) of + (f, f') -> \(Norm n) -> Norm $ f' . n . f + coerceVariance :: Hask.Functor p => p (x,ξ) -> Metric' x -> Metric' ξ + coerceVariance p = case ( oppositeLocalCoercion :: CanonicalDiffeomorphism ξ x + , dualSpaceWitness :: DualSpaceWitness (Needle x) + , dualSpaceWitness :: DualSpaceWitness (Needle ξ) ) of + (CanonicalDiffeomorphism, DualSpaceWitness, DualSpaceWitness) + -> case ( coerceNeedle p, coerceNeedle' (swap<$>p) ) of + (f, f') -> \(Norm n) -> Norm $ f . n . f' oppositeLocalCoercion :: CanonicalDiffeomorphism ξ x default oppositeLocalCoercion :: LocallyCoercible ξ x => CanonicalDiffeomorphism ξ x oppositeLocalCoercion = CanonicalDiffeomorphism @@ -443,3 +460,10 @@ type DualNeedleWitness x = DualSpaceWitness (Needle x) + + +infix 6 ⊙+^ +-- | Proxy-version of `translateP`. +(⊙+^) :: ∀ x proxy . Semimanifold x => Interior x -> Needle x -> proxy x -> Interior x +(⊙+^) x v _ = tp x v + where Tagged tp = translateP :: Tagged x (Interior x -> Needle x -> Interior x) diff --git a/Data/Manifold/Shade.hs b/Data/Manifold/Shade.hs --- a/Data/Manifold/Shade.hs +++ b/Data/Manifold/Shade.hs @@ -44,17 +44,10 @@ , factoriseShade, orthoShades, (✠), intersectShade's, linIsoTransformShade , embedShade, projectShade , Refinable, subShade', refineShade', convolveShade', coerceShade - , mixShade's, dualShade - -- * Misc - -- ** Shades + , mixShade's, dualShade, dualShade', wellDefinedShade', linearProjectShade , shadesMerge, pointsShades', pseudoECM, convolveMetric , WithAny(..), shadeWithAny, shadeWithoutAnything - -- ** Local data fit models - , estimateLocalJacobian, estimateLocalHessian, QuadraticModel(..) - -- ** Differential equations - , DifferentialEqn, LocalDifferentialEqn(..) - , propagateDEqnSolution_loc, LocalDataPropPlan(..) - -- ** Range interpolation + -- * Misc , rangeOnGeodesic, rangeWithinVertices ) where @@ -99,6 +92,7 @@ import GHC.Generics (Generic) import Text.Show.Number +import qualified Text.Show.Pragmatic as SP -- | A 'Shade' is a very crude description of a region within a manifold. It @@ -108,8 +102,10 @@ -- -- For a /precise/ description of an arbitrarily-shaped connected subset of a manifold, -- there is 'Region', whose implementation is vastly more complex. -data Shade x = Shade { _shadeCtr :: !(Interior x) - , _shadeExpanse :: !(Metric' x) } +data Shade x where + Shade :: (Semimanifold x, SimpleSpace (Needle x)) + => { _shadeCtr :: !(Interior x) + , _shadeExpanse :: !(Metric' x) } -> Shade x deriving instance (Show (Interior x), Show (Metric' x), WithField ℝ PseudoAffine x) => Show (Shade x) @@ -119,28 +115,7 @@ data Shade' x = Shade' { _shade'Ctr :: !(Interior x) , _shade'Narrowness :: !(Metric x) } -data LocalDifferentialEqn x ð y = LocalDifferentialEqn { - _predictDerivatives :: Shade' ð -> Maybe (Shade' (LocalLinear x y)) - , _rescanDerivatives :: Shade' y -> Shade' (LocalLinear x y) - -> Shade' (LocalBilinear x y) - -> (Maybe (Shade' y), Maybe (Shade' ð)) - } -makeLenses ''LocalDifferentialEqn -type DifferentialEqn x ð y = Shade (x,y) -> LocalDifferentialEqn x ð y - -data LocalDataPropPlan x ym yr = LocalDataPropPlan - { _sourcePosition :: !(Interior x) - , _targetPosOffset :: !(Needle x) - , _sourceData, _targetAPrioriData :: !ym - , _relatedData :: [(Needle x, yr)] - } -deriving instance (Show (Interior x), Show ym, Show yr, Show (Needle x)) - => Show (LocalDataPropPlan x ym yr) - -makeLenses ''LocalDataPropPlan - - class IsShade shade where -- type (*) shade :: *->* -- | Access the center of a 'Shade' or a 'Shade''. @@ -156,7 +131,8 @@ , PseudoAffine y, SimpleSpace (Needle y) , Scalar (Needle x) ~ Scalar (Needle y) ) => shade (x,y) -> (shade x, shade y) - coerceShade :: (Manifold x, Manifold y, LocallyCoercible x y) => shade x -> shade y + coerceShade :: ( Manifold x, Manifold y, LocallyCoercible x y + , SimpleSpace (Needle y) ) => shade x -> shade y -- | ASCII version of '✠'. orthoShades :: ( PseudoAffine x, SimpleSpace (Needle x) , PseudoAffine y, SimpleSpace (Needle y) @@ -168,7 +144,7 @@ -- | Squash a shade down into a lower dimensional space. projectShade :: ( Semimanifold x, Semimanifold y , Object (Affine s) (Interior x), Object (Affine s) (Interior y) - , SemiInner (Needle x), SemiInner (Needle y) ) + , SimpleSpace (Needle x), SemiInner (Needle y) ) => Embedding (Affine s) (Interior x) (Interior y) -> shade y -> shade x -- | Include a shade in a higher-dimensional space. Notice that this behaves @@ -177,13 +153,13 @@ -- pillar” pointing in the projection's orthogonal complement. embedShade :: ( Semimanifold x, Semimanifold y , Object (Affine s) (Interior x), Object (Affine s) (Interior y) - , SemiInner (Needle x), SemiInner (Needle y) ) + , SemiInner (Needle x), SimpleSpace (Needle y) ) => Embedding (Affine s) (Interior x) (Interior y) -> shade x -> shade y linearProjectShade :: ∀ s x y - . (Num' s, LinearSpace x, LinearSpace y, Scalar x ~ s, Scalar y ~ s) + . (Num' s, LinearSpace x, SimpleSpace y, Scalar x ~ s, Scalar y ~ s) => (x+>y) -> Shade x -> Shade y linearProjectShade = case ( linearManifoldWitness :: LinearManifoldWitness x , linearManifoldWitness :: LinearManifoldWitness y @@ -235,7 +211,7 @@ fs DualSpaceWitness DualSpaceWitness (Shade x δx) (Shade y δy) = Shade (x,y) $ sumSubspaceNorms δx δy coerceShade = cS dualSpaceWitness dualSpaceWitness - where cS :: ∀ x y . (LocallyCoercible x y) + where cS :: ∀ x y . (LocallyCoercible x y, SimpleSpace (Needle y)) => DualNeedleWitness x -> DualNeedleWitness y -> Shade x -> Shade y cS DualSpaceWitness DualSpaceWitness = \(Shade x δxym) -> Shade (internCoerce x) (tN δxym) @@ -246,7 +222,7 @@ CanonicalDiffeomorphism -> locallyTrivialDiffeomorphism linIsoTransformShade = lits linearManifoldWitness linearManifoldWitness dualSpaceWitness dualSpaceWitness - where lits :: ∀ x y . ( LinearSpace x, LinearSpace y + where lits :: ∀ x y . ( LinearSpace x, SimpleSpace y , Scalar x ~ Scalar y, Num' (Scalar x) ) => LinearManifoldWitness x -> LinearManifoldWitness y -> DualSpaceWitness x -> DualSpaceWitness y @@ -257,8 +233,9 @@ f (Shade x δx) = Shade (f $ x) (transformNorm (adjoint $ f) δx) embedShade = ps' (semimanifoldWitness, semimanifoldWitness) - where ps' :: ∀ s x y . ( Object (Affine s) (Interior x), Object (Affine s) (Interior y) - , SemiInner (Needle x), SemiInner (Needle y) ) + where ps' :: ∀ s x y . ( Semimanifold y + , Object (Affine s) (Interior x), Object (Affine s) (Interior y) + , SemiInner (Needle x), SimpleSpace (Needle y) ) => (SemimanifoldWitness x, SemimanifoldWitness y) -> Embedding (Affine s) (Interior x) (Interior y) -> Shade x -> Shade y @@ -267,8 +244,9 @@ where y = q $ x (_,j) = evalAffine q x projectShade = ps' (semimanifoldWitness, semimanifoldWitness) - where ps' :: ∀ s x y . ( Object (Affine s) (Interior x), Object (Affine s) (Interior y) - , SemiInner (Needle x), SemiInner (Needle y) ) + where ps' :: ∀ s x y . ( Semimanifold x + , Object (Affine s) (Interior x), Object (Affine s) (Interior y) + , SimpleSpace (Needle x), SemiInner (Needle y) ) => (SemimanifoldWitness x, SemimanifoldWitness y) -> Embedding (Affine s) (Interior x) (Interior y) -> Shade y -> Shade x @@ -412,10 +390,11 @@ Just pinterp = case geodesicWitness :: GeodesicWitness x of GeodesicWitness _ -> geodesicBetween c ζ -fullShade :: WithField ℝ PseudoAffine x => Interior x -> Metric' x -> Shade x +fullShade :: (Semimanifold x, SimpleSpace (Needle x)) + => Interior x -> Metric' x -> Shade x fullShade ctr expa = Shade ctr expa -fullShade' :: WithField ℝ PseudoAffine x => Interior x -> Metric x -> Shade' x +fullShade' :: WithField ℝ SimpleSpace x => Interior x -> Metric x -> Shade' x fullShade' ctr expa = Shade' ctr expa @@ -425,11 +404,11 @@ #if GLASGOW_HASKELL < 800 pattern (:±) :: () #else -pattern (:±) :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) +pattern (:±) :: (Semimanifold x, SimpleSpace (Needle x)) #endif - => (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => (Semimanifold x, SimpleSpace (Needle x)) => Interior x -> [Needle x] -> Shade x -pattern x :± shs <- Shade x (varianceSpanningSystem -> shs) +pattern x :± shs <- (Shade x (varianceSpanningSystem -> shs)) where x :± shs = fullShade x $ spanVariance shs -- | Similar to ':±', but instead of expanding the shade, each vector /restricts/ it. @@ -459,7 +438,7 @@ -- Hence the result type is a list. pointsShades :: (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) => [Interior x] -> [Shade x] -pointsShades = map snd . pointsShades' mempty . map fromInterior +pointsShades = map snd . pointsShades' mempty . map ((,()) . fromInterior) coverAllAround :: ∀ x s . ( Fractional' s, WithField s PseudoAffine x , SimpleSpace (Needle x) ) @@ -488,10 +467,11 @@ pointsCovers = case pseudoAffineWitness :: PseudoAffineWitness x of (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) -> \ps -> map (\(ps', Shade x₀ _) - -> coverAllAround x₀ [v | p<-ps' + -> coverAllAround x₀ [v | (p,())<-ps' , let Just v = p.-~.fromInterior x₀]) - (pointsShades' mempty (fromInterior<$>ps) :: [([x], Shade x)]) + (pointsShades' mempty ((,()).fromInterior<$>ps) + :: [([(x,())], Shade x)]) pointsShade's :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) => [Interior x] -> [Shade' x] @@ -503,20 +483,20 @@ pointsCover's = case dualSpaceWitness :: DualNeedleWitness x of DualSpaceWitness -> map (\(Shade c e :: Shade x) -> Shade' c $ dualNorm e) . pointsCovers -pseudoECM :: ∀ x p . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x), Hask.Functor p) - => p x -> NonEmpty x -> (x, ([x],[x])) +pseudoECM :: ∀ x y p . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x), Hask.Functor p) + => p x -> NonEmpty (x,y) -> (x, ([(x,y)],[(x,y)])) pseudoECM = case semimanifoldWitness :: SemimanifoldWitness x of SemimanifoldWitness _ -> - \_ (p₀ NE.:| psr) -> foldl' ( \(acc, (rb,nr)) (i,p) + \_ ((p₀,y₀) NE.:| psr) -> foldl' ( \(acc, (rb,nr)) (i,(p,y)) -> case (p.-~.acc, toInterior acc) of (Just δ, Just acci) - -> (acci .+~^ δ^/i, (p:rb, nr)) - _ -> (acc, (rb, p:nr)) ) + -> (acci .+~^ δ^/i, ((p,y):rb, nr)) + _ -> (acc, (rb, (p,y):nr)) ) (p₀, mempty) - ( zip [1..] $ p₀:psr ) + ( zip [1..] $ (p₀,y₀):psr ) -pointsShades' :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => Metric' x -> [x] -> [([x], Shade x)] +pointsShades' :: ∀ x y . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) + => Metric' x -> [(x,y)] -> [([(x,y)], Shade x)] pointsShades' _ [] = [] pointsShades' minExt ps = case (expa, toInterior ctr) of (Just e, Just c) @@ -525,7 +505,7 @@ ++ pointsShades' minExt unreachable where (ctr,(inc'd,unreachable)) = pseudoECM ([]::[x]) $ NE.fromList ps expa = ( (<>minExt) . spanVariance . map (^/ fromIntegral (length ps)) ) - <$> mapM (.-~.ctr) ps + <$> mapM ((.-~.ctr) . fst) ps -- | Attempt to reduce the number of shades to fewer (ideally, a single one). @@ -931,143 +911,14 @@ intersectShade's (sh:|shs) = Hask.foldrM refineShade' sh shs -estimateLocalJacobian :: ∀ x y . ( WithField ℝ Manifold x, Refinable y - , SimpleSpace (Needle x), SimpleSpace (Needle y) ) - => Metric x -> [(Local x, Shade' y)] - -> Maybe (Shade' (LocalLinear x y)) -estimateLocalJacobian = elj ( pseudoAffineWitness :: PseudoAffineWitness x - , pseudoAffineWitness :: PseudoAffineWitness y ) - where elj ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) - , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) ) - mex [(Local x₁, Shade' y₁ ey₁),(Local x₀, Shade' y₀ ey₀)] - = return $ Shade' (dx-+|>δy) - (Norm . LinearFunction $ \δj -> δx ⊗ (σey<$|δj $ δx)) - where Just δx = x₁.-~.x₀ - δx' = (mex<$|δx) - dx = δx'^/(δx'<.>^δx) - Just δy = y₁.-~.y₀ - σey = convolveMetric ([]::[y]) ey₀ ey₁ - elj _ mex (po:ps) - | DualSpaceWitness <- dualSpaceWitness :: DualNeedleWitness y - , length ps > 1 - = mixShade's =<< (:|) <$> estimateLocalJacobian mex ps - <*> sequenceA [estimateLocalJacobian mex [po,pi] | pi<-ps] - elj _ _ _ = return $ Shade' zeroV mempty -data QuadraticModel x y = QuadraticModel { - _quadraticModelOffset :: Interior y - , _quadraticModel :: Shade (Needle y, (Needle x+>Needle y, Needle x⊗〃+>Needle y)) - } -quadratic_linearRegression :: ∀ s x y . - ( WithField s PseudoAffine x - , WithField s PseudoAffine y, Geodesic y - , SimpleSpace (Needle x), SimpleSpace (Needle y) ) - => NE.NonEmpty (Needle x, Shade' y) -> QuadraticModel x y -quadratic_linearRegression = qlr - ( dualSpaceWitness, pseudoAffineWitness - , linearManifoldWitness, dualSpaceWitness - , geodesicWitness ) - where qlr :: ( DualSpaceWitness (Needle x) - , PseudoAffineWitness y, LinearManifoldWitness (Needle y) - , DualSpaceWitness (Needle y) - , GeodesicWitness y ) - -> NE.NonEmpty (Needle x, Shade' y) -> QuadraticModel x y - qlr ( DualSpaceWitness - , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) - , LinearManifoldWitness BoundarylessWitness, DualSpaceWitness - , GeodesicWitness _ ) ps - = QuadraticModel cmy - $ coverAllAround mBest (convexPolytopeRepresentatives dm) - where Just cmy = pointsBarycenter $ _shade'Ctr.snd<$>ps - Just vsxy = Hask.mapM (\(x, Shade' y ey) -> (x,).(,ey)<$>y.-~.cmy) ps - (mBest :: ( Needle y, (Needle x+>Needle y - , SymmetricTensor s (Needle x)+>(Needle y)) - ) - , dm) - = linearRegressionWVar - (\δx -> lfun $ \(c,(b,a)) -> (a $ squareV δx) - ^+^ (b $ δx) ^+^ c ) - (NE.toList vsxy) -estimateLocalHessian :: ∀ x y . ( WithField ℝ Manifold x, Refinable y, Geodesic y - , FlatSpace (Needle x), FlatSpace (Needle y) ) - => NonEmpty (Local x, Shade' y) -> QuadraticModel x y -estimateLocalHessian pts = elj ( pseudoAffineWitness :: PseudoAffineWitness x - , pseudoAffineWitness :: PseudoAffineWitness y ) - where elj ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) - , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) ) - = theModel - where localPts :: NonEmpty (Needle x, Shade' y) - localPts = pts >>= \(Local x, Shade' y ey) - -> NE.fromList [ (x, Shade' (y.+~^σ*^δy) ey) - | δy <- normSpanningSystem' ey - , σ <- [-1,1] ] - theModel = quadratic_linearRegression localPts -propagateDEqnSolution_loc :: ∀ x y ð . ( WithField ℝ Manifold x - , Refinable y, Geodesic (Interior y) - , WithField ℝ AffineManifold ð, Geodesic ð - , SimpleSpace (Needle x), SimpleSpace (Needle ð) ) - => DifferentialEqn x ð y - -> LocalDataPropPlan x (Shade' y, Shade' ð) (Shade' y) - -> Maybe (Shade' y) -propagateDEqnSolution_loc f propPlan - = pdesl (dualSpaceWitness :: DualNeedleWitness x) - (dualSpaceWitness :: DualNeedleWitness y) - (boundarylessWitness :: BoundarylessWitness x) - (pseudoAffineWitness :: PseudoAffineWitness y) - where pdesl DualSpaceWitness DualSpaceWitness BoundarylessWitness - (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) - | Nothing <- jacobian = Nothing - | otherwise = pure result - where jacobian = (f shxy ^. predictDerivatives $ shð) - >>= \j -> mixShade's $ j:|[aprioriDirDrv] - Just (Shade' j₀ jExpa) = jacobian - jacobianSh :: Shade (LocalLinear x y) - Just jacobianSh = dualShade' <$> jacobian - aprioriDirDrv :: Shade' (LocalLinear x y) - Just aprioriDirDrv = estimateLocalJacobian expax - [ (Local zeroV :: Local x, propPlan^.sourceData._1) - , (Local δx, propPlan^.targetAPrioriData._1) ] - mx = propPlan^.sourcePosition .+~^ propPlan^.targetPosOffset ^/ 2 :: x - Just shð = middleBetween (propPlan^.sourceData._2) - (propPlan^.targetAPrioriData._2) - shxy = coverAllAround (mx, mυ) - [ (δx ^-^ propPlan^.targetPosOffset ^/ 2, pυ ^+^ v) - | (δx,neυ) <- (zeroV, propPlan^.sourceData._1) - : (second id - <$> propPlan^.relatedData) - , let Just pυ = neυ^.shadeCtr .-~. mυ - , v <- normSpanningSystem' (neυ^.shadeNarrowness) - ] - where Just mυ = middleBetween (propPlan^.sourceData._1.shadeCtr) - (propPlan^.targetAPrioriData._1.shadeCtr) - (Shade _ expax' :: Shade x) - = coverAllAround (propPlan^.sourcePosition) - [δx | (δx,_) <- propPlan^.relatedData] - expax = dualNorm expax' - result :: Shade' y - Just result = wellDefinedShade' $ convolveShade' - (case wellDefinedShade' $ propPlan^.sourceData._1 of {Just s->s}) - (case wellDefinedShade' . dualShade - . linearProjectShade (lfun ($ δx)) - $ jacobianSh - of {Just s->s}) - where δyb = j₀ $ δx - δx = propPlan^.targetPosOffset - - - - - - - - -- | Essentially the same as @(x,y)@, but not considered as a product topology. -- The 'Semimanifold' etc. instances just copy the topology of @x@, ignoring @y@. data x`WithAny`y @@ -1123,7 +974,7 @@ shadeWithAny :: y -> Shade x -> Shade (x`WithAny`y) shadeWithAny y (Shade x xe) = Shade (WithAny y x) xe -shadeWithoutAnything :: Shade (x`WithAny`y) -> Shade x +shadeWithoutAnything :: Semimanifold x => Shade (x`WithAny`y) -> Shade x shadeWithoutAnything (Shade (WithAny _ b) e) = Shade b e @@ -1140,6 +991,10 @@ prettyShowShade' :: LtdErrorShow x => Shade' x -> String prettyShowShade' sh = prettyShowsPrecShade' 0 sh [] +instance LtdErrorShow x => SP.Show (Shade' x) where + showsPrec = prettyShowsPrecShade' +instance LtdErrorShow x => SP.Show (Shade x) where + showsPrec = prettyShowsPrecShade wellDefinedShade' :: LinearSpace (Needle x) => Shade' x -> Maybe (Shade' x) @@ -1157,6 +1012,16 @@ => LtdErrorShowWitness m ltdErrorShowWitness = LtdErrorShowWitness pseudoAffineWitness showsPrecShade'_errorLtdC :: Int -> Shade' m -> ShowS + prettyShowsPrecShade :: Int -> Shade m -> ShowS + prettyShowsPrecShade p sh@(Shade c e') + = showParen (p>6) $ v + . (":±["++) . flip (foldr id) (intersperse (',':) u) . (']':) + where v = showsPrecShade'_errorLtdC 6 (Shade' c e :: Shade' m) + u :: [ShowS] = case ltdErrorShowWitness :: LtdErrorShowWitness m of + LtdErrorShowWitness (PseudoAffineWitness (SemimanifoldWitness _)) -> + [ showsPrecShade'_errorLtdC 6 (Shade' δ e :: Shade' (Needle m)) + | δ <- varianceSpanningSystem e'] + e = dualNorm' e' prettyShowsPrecShade' :: Int -> Shade' m -> ShowS prettyShowsPrecShade' p sh@(Shade' c e) = showParen (p>6) $ v @@ -1213,6 +1078,25 @@ where (shx,shy) = factoriseShade sh shshx = showsPrecShade'_errorLtdC 0 shx shshy = showsPrecShade'_errorLtdC 0 shy + +instance ∀ v . + (HilbertSpace v, SemiInner v, FiniteDimensional v, LtdErrorShow v, Scalar v ~ ℝ) + => LtdErrorShow (LinearMap ℝ v ℝ) where + showsPrecShade'_errorLtdC p sh = showParen (p>7) $ + ("().<"++) . showsPrecShade'_errorLtdC 7 + (linIsoTransformShade (arr fromLinearForm) sh :: Shade' v) +instance ∀ v . + (HilbertSpace v, SemiInner v, FiniteDimensional v, LtdErrorShow v, Scalar v ~ ℝ) + => LtdErrorShow (LinearMap ℝ v (ℝ,ℝ)) where + showsPrecShade'_errorLtdC p sh = showParen (p>7) $ + ( "Left ().<"++) . showsPrecShade'_errorLtdC 7 shx + . ("^+^Right().<"++) . showsPrecShade'_errorLtdC 7 shy + where (shx,shy) = factoriseShade + (linIsoTransformShade (lfun $ \f + -> ( fromLinearForm $ fst . f + , fromLinearForm $ snd . f ) ) sh + :: Shade' (v,v)) + instance LtdErrorShow x => Show (Shade' x) where showsPrec = prettyShowsPrecShade' diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs --- a/Data/Manifold/TreeCover.hs +++ b/Data/Manifold/TreeCover.hs @@ -45,22 +45,21 @@ , Refinable, subShade', refineShade', convolveShade', coerceShade , mixShade's -- * Shade trees - , ShadeTree(..), fromLeafPoints, onlyLeaves, indexShadeTree, positionIndex - -- * View helpers - , onlyNodes + , ShadeTree, fromLeafPoints, fromLeafPoints_, onlyLeaves, onlyLeaves_ + , indexShadeTree, treeLeaf, positionIndex + -- ** View helpers + , entireTree, onlyNodes, trunkBranches, nLeaves, treeDepth -- ** Auxiliary types - , SimpleTree, Trees, NonEmptyTree, GenericTree(..) + , SimpleTree, Trees, NonEmptyTree, GenericTree(..), 朳 -- * Misc - , HasFlatView(..), shadesMerge, smoothInterpolate + , HasFlatView(..), shadesMerge , allTwigs, twigsWithEnvirons, Twig, TwigEnviron, seekPotentialNeighbours - , completeTopShading, flexTwigsShading, coerceShadeTree - , WithAny(..), Shaded, fmapShaded, joinShaded - , constShaded, zipTreeWithList, stripShadedUntopological + , completeTopShading, flexTwigsShading, traverseTrunkBranchChoices + , Shaded(..), fmapShaded + , constShaded, zipTreeWithList , stiAsIntervalMapping, spanShading - , estimateLocalJacobian - , DifferentialEqn, LocalDifferentialEqn(..) - , propagateDEqnSolution_loc, LocalDataPropPlan(..) - , rangeOnGeodesic + , DBranch, DBranch'(..), Hourglass(..) + , unsafeFmapTree -- ** Triangulation-builders , TriangBuild, doTriangBuild , AutoTriang, breakdownAutoTriang @@ -106,6 +105,7 @@ import Data.Functor.Identity import Control.Monad.Trans.State import Control.Monad.Trans.Writer +import Control.Monad.Trans.List import Control.Monad.Trans.OuterMaybe import Control.Monad.Trans.Class import qualified Data.Foldable as Hask @@ -123,6 +123,7 @@ import GHC.Generics (Generic) import Data.Type.Coercion +import Development.Placeholders type Depth = Int @@ -158,7 +159,7 @@ -- | Hourglass as the geometric shape (two opposing ~conical volumes, sharing -- only a single point in the middle); has nothing to do with time. data Hourglass s = Hourglass { upperBulb, lowerBulb :: !s } - deriving (Generic, Hask.Functor, Hask.Foldable, Show) + deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable, Show) instance (NFData s) => NFData (Hourglass s) instance (Semigroup s) => Semigroup (Hourglass s) where Hourglass u l <> Hourglass u' l' = Hourglass (u<>u') (l<>l') @@ -184,24 +185,28 @@ oneBulb LowerBulb f (Hourglass u l) = Hourglass u (f l) +type LeafCount = Int +type LeafIndex = Int -data ShadeTree x = PlainLeaves [x] - | DisjointBranches !Int (NonEmpty (ShadeTree x)) - | OverlappingBranches !Int !(Shade x) (NonEmpty (DBranch x)) - deriving (Generic) +type ShadeTree x = x`Shaded`() + +data Shaded x y = PlainLeaves [(x,y)] + | DisjointBranches !LeafCount (NonEmpty (x`Shaded`y)) + | OverlappingBranches !LeafCount !(Shade x) (NonEmpty (DBranch x y)) + deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable) deriving instance ( WithField ℝ PseudoAffine x, Show x , Show (Interior x), Show (Needle' x), Show (Metric' x) ) => Show (ShadeTree x) data DBranch' x c = DBranch { boughDirection :: !(Needle' x) , boughContents :: !(Hourglass c) } - deriving (Generic, Hask.Functor, Hask.Foldable) -type DBranch x = DBranch' x (ShadeTree x) + deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable) +type DBranch x y = DBranch' x (x`Shaded`y) deriving instance ( WithField ℝ PseudoAffine x, Show (Needle' x), Show c ) => Show (DBranch' x c) newtype DBranches' x c = DBranches (NonEmpty (DBranch' x c)) - deriving (Generic, Hask.Functor, Hask.Foldable) + deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable) deriving instance ( WithField ℝ PseudoAffine x, Show (Needle' x), Show c ) => Show (DBranches' x c) @@ -210,17 +215,31 @@ DBranches b1 <> DBranches b2 = DBranches $ NE.zipWith (\(DBranch d1 c1) (DBranch _ c2) -> DBranch d1 $ c1<>c2 ) b1 b2 + + +trunkBranches :: x`Shaded`y -> NonEmpty (LeafIndex, x`Shaded`y) +trunkBranches (OverlappingBranches _ _ brs) + = (`evalState`0) + . forM (brs >>= \(DBranch _ (Hourglass t b)) -> t:|[b]) $ \st -> do + i₀ <- get + put $ i₀ + nLeaves st + return (i₀, st) +trunkBranches (DisjointBranches _ brs) = (`evalState`0) . forM brs $ \st -> do + i₀ <- get + put $ i₀ + nLeaves st + return (i₀, st) +trunkBranches t = pure (0,t) directionChoices :: WithField ℝ Manifold x - => [DBranch x] - -> [ ( (Needle' x, ShadeTree x) - ,[(Needle' x, ShadeTree x)] ) ] + => [DBranch x y] + -> [ ( (Needle' x, x`Shaded`y) + ,[(Needle' x, x`Shaded`y)] ) ] directionChoices = map (snd *** map snd) . directionIChoices 0 directionIChoices :: (WithField ℝ PseudoAffine x, AdditiveGroup (Needle' x)) - => Int -> [DBranch x] - -> [ ( (Int, (Needle' x, ShadeTree x)) - ,[(Int, (Needle' x, ShadeTree x))] ) ] + => Int -> [DBranch x y] + -> [ ( (Int, (Needle' x, x`Shaded`y)) + ,[(Int, (Needle' x, x`Shaded`y))] ) ] directionIChoices _ [] = [] directionIChoices i₀ (DBranch ѧ (Hourglass t b) : hs) = ( top, bot : map fst uds ) @@ -230,13 +249,12 @@ bot = (i₀+1,(negateV ѧ,b)) uds = directionIChoices (i₀+2) hs -traverseDirectionChoices :: ( WithField ℝ PseudoAffine x, LSpace (Needle x) - , Hask.Applicative f ) - => ( (Int, (Needle' x, ShadeTree x)) - -> [(Int, (Needle' x, ShadeTree x))] - -> f (ShadeTree x) ) - -> [DBranch x] - -> f [DBranch x] +traverseDirectionChoices :: ( AdditiveGroup (Needle' x), Hask.Applicative f ) + => ( (Int, (Needle' x, x`Shaded`y)) + -> [(Int, (Needle' x, x`Shaded`y))] + -> f (x`Shaded`z) ) + -> [DBranch x y] + -> f [DBranch x z] traverseDirectionChoices f dbs = td [] . scanLeafNums 0 $ dbs >>= \(DBranch ѧ (Hourglass τ β)) @@ -252,7 +270,25 @@ scanLeafNums i₀ ((v,t):vts) = (i₀, (v,t)) : scanLeafNums (i₀ + nLeaves t) vts -indexDBranches :: NonEmpty (DBranch x) -> NonEmpty (DBranch' x (Int, ShadeTree x)) + +traverseTrunkBranchChoices :: Hask.Applicative f + => ( (Int, x`Shaded`y) -> x`Shaded`y -> f (x`Shaded`z) ) + -> x`Shaded`y -> f (x`Shaded`z) +traverseTrunkBranchChoices f (OverlappingBranches n sh bs) + = OverlappingBranches n sh . NE.fromList <$> go 0 id (NE.toList bs) + where go _ _ [] = pure [] + go i₀ prbs (tbs@(DBranch v (Hourglass τ β)) : dbs) + = (:) . DBranch v <$> + (Hourglass <$> (f (i₀, τ) . OverlappingBranches (n-nτ) sh + . NE.fromList . prbs $ DBranch v (Hourglass hole β) : dbs) + <*> (f (i₀+nτ, β) . OverlappingBranches (n-nβ) sh + . NE.fromList . prbs $ DBranch v (Hourglass τ hole) : dbs)) + <*> go (i₀+nτ+nβ) (prbs . (tbs:)) dbs + where [nτ, nβ] = nLeaves<$>[τ,β] + hole = PlainLeaves [] + + +indexDBranches :: NonEmpty (DBranch x y) -> NonEmpty (DBranch' x (Int, x`Shaded`y)) indexDBranches (DBranch d (Hourglass t b) :| l) -- this could more concisely be written as a traversal = DBranch d (Hourglass (0,t) (nt,b)) :| ixDBs (nt + nb) l where nt = nLeaves t; nb = nLeaves b @@ -261,31 +297,18 @@ = DBranch δ (Hourglass (i₀,τ) (i₀+nτ,β)) : ixDBs (i₀ + nτ + nβ) l where nτ = nLeaves τ; nβ = nLeaves β -instance (NFData x, NFData (Needle' x)) => NFData (ShadeTree x) where +instance (NFData x, NFData (Needle' x), NFData y) => NFData (x`Shaded`y) where rnf (PlainLeaves xs) = rnf xs rnf (DisjointBranches n bs) = n `seq` rnf (NE.toList bs) rnf (OverlappingBranches n sh bs) = n `seq` sh `seq` rnf (NE.toList bs) -instance (NFData x, NFData (Needle' x)) => NFData (DBranch x) +instance (NFData x, NFData (Needle' x), NFData y) => NFData (DBranch x y) --- | Experimental. There might be a more powerful instance possible. -instance (AffineManifold x) => Semimanifold (ShadeTree x) where - type Needle (ShadeTree x) = Diff x - fromInterior = id - toInterior = pure - translateP = Tagged (.+~^) - PlainLeaves xs .+~^ v = PlainLeaves $ (.+^v)<$>xs - OverlappingBranches n sh br .+~^ v - = OverlappingBranches n (sh.+~^v) - $ fmap (\(DBranch d c) -> DBranch d $ (.+~^v)<$>c) br - DisjointBranches n br .+~^ v = DisjointBranches n $ (.+~^v)<$>br - semimanifoldWitness = case semimanifoldWitness :: SemimanifoldWitness x of - SemimanifoldWitness BoundarylessWitness -> SemimanifoldWitness BoundarylessWitness -- | WRT union. instance (WithField ℝ Manifold x, SimpleSpace (Needle x)) => Semigroup (ShadeTree x) where PlainLeaves [] <> t = t t <> PlainLeaves [] = t - t <> s = fromLeafPoints $ onlyLeaves t ++ onlyLeaves s + t <> s = fromLeafPoints $ onlyLeaves_ t ++ onlyLeaves_ s -- Could probably be done more efficiently sconcat = mconcat . NE.toList instance (WithField ℝ Manifold x, SimpleSpace (Needle x)) => Monoid (ShadeTree x) where @@ -294,7 +317,7 @@ mconcat l = case filter ne l of [] -> mempty [t] -> t - l' -> fromLeafPoints $ onlyLeaves =<< l' + l' -> fromLeafPoints $ onlyLeaves_ =<< l' where ne (PlainLeaves []) = False; ne _ = True @@ -304,14 +327,17 @@ -- -- <> fromLeafPoints :: ∀ x. (WithField ℝ Manifold x, SimpleSpace (Needle x)) - => [x] -> ShadeTree x -fromLeafPoints = fromLeafPoints' sShIdPartition + => [x] -> ShadeTree x +fromLeafPoints = fromLeafPoints_ . map (,()) +fromLeafPoints_ :: ∀ x y. (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => [(x,y)] -> x`Shaded`y +fromLeafPoints_ = fromLeafPoints' sShIdPartition + -- | The leaves of a shade tree are numbered. For a given index, this function -- attempts to find the leaf with that ID, within its immediate environment. -indexShadeTree :: ∀ x . WithField ℝ Manifold x - => ShadeTree x -> Int -> Either Int ([ShadeTree x], x) +indexShadeTree :: ∀ x y . x`Shaded`y -> Int -> Either Int ([x`Shaded`y], (x,y)) indexShadeTree _ i | i<0 = Left i indexShadeTree sh@(PlainLeaves lvs) i = case length lvs of @@ -330,21 +356,42 @@ ) (Left i) (toList brs>>=toList) | otherwise = Left $ i-n +treeLeaf :: ∀ x y f . Hask.Functor f + => Int -> (y -> f y) -> x`Shaded`y -> Either Int (f (x`Shaded`y)) +treeLeaf i _ _ + | i<0 = Left i +treeLeaf i f sh@(PlainLeaves lvs) = case length lvs of + n | i Right . fmap (PlainLeaves . (pre++) . (:post) . (x,)) $ f node + | otherwise -> Left $ i-n +treeLeaf i f (DisjointBranches n brs) + | i (treeLeaf i' f) + result -> return result + ) (Left i) brs + | otherwise = Left $ i-n +treeLeaf i f sh@(OverlappingBranches n _ brs) + | i (treeLeaf i' f) + result -> return result + ) (Left i) (toList brs>>=toList) + | otherwise = Left $ i-n -- | “Inverse indexing” of a tree. This is roughly a nearest-neighbour search, -- but not guaranteed to give the correct result unless evaluated at the -- precise position of a tree leaf. -positionIndex :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x)) +positionIndex :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => Maybe (Metric x) -- ^ For deciding (at the lowest level) what “close” means; -- this is optional for any tree of depth >1. - -> ShadeTree x -- ^ The tree to index into + -> x`Shaded`y -- ^ The tree to index into -> x -- ^ Position to look up - -> Maybe (Int, ([ShadeTree x], x)) + -> Maybe (Int, ([x`Shaded`y], (x,y))) -- ^ Index of the leaf near to the query point, the “path” of -- environment trees leading down to its position (in decreasing - -- order of size), and actual position of the found node. + -- order of size), and actual position+info of the found node. positionIndex (Just m) sh@(PlainLeaves lvs) x - = case catMaybes [ ((i,p),) . normSq m <$> p.-~.x + = case catMaybes [ ((i,p),) . normSq m <$> fst p.-~.x | (i,p) <- zip [0..] lvs] of [] -> empty l | ((i,p),_) <- minimumBy (comparing snd) l @@ -371,27 +418,14 @@ -fromFnGraphPoints :: ∀ x y . ( WithField ℝ Manifold x, WithField ℝ Manifold y - , SimpleSpace (Needle x), SimpleSpace (Needle y) ) - => [(x,y)] -> ShadeTree (x,y) -fromFnGraphPoints = case ( dualSpaceWitness :: DualNeedleWitness x - , boundarylessWitness :: BoundarylessWitness x - , dualSpaceWitness :: DualNeedleWitness y - , boundarylessWitness :: BoundarylessWitness y ) of - (DualSpaceWitness,BoundarylessWitness,DualSpaceWitness,BoundarylessWitness) - -> fromLeafPoints' $ - \(Shade c expa) xs -> case - [ DBranch (v, zeroV) mempty - | v <- normSpanningSystem' (transformNorm (id&&&zeroV) expa :: Metric' x) ] of - (b:bs) -> sShIdPartition' c xs $ b:|bs -fromLeafPoints' :: ∀ x. (WithField ℝ Manifold x, SimpleSpace (Needle x)) => - (Shade x -> [x] -> NonEmpty (DBranch' x [x])) -> [x] -> ShadeTree x +fromLeafPoints' :: ∀ x y. (WithField ℝ Manifold x, SimpleSpace (Needle x)) => + (Shade x -> [(x,y)] -> NonEmpty (DBranch' x [(x,y)])) -> [(x,y)] -> x`Shaded`y fromLeafPoints' sShIdPart = go boundarylessWitness mempty - where go :: BoundarylessWitness x -> Metric' x -> [x] -> ShadeTree x + where go :: BoundarylessWitness x -> Metric' x -> [(x,y)] -> x`Shaded`y go bw@BoundarylessWitness preShExpa = \xs -> case pointsShades' (scaleNorm (1/3) preShExpa) xs of - [] -> mempty + [] -> PlainLeaves [] [(_,rShade)] -> let trials = sShIdPart rShade xs in case reduce rShade trials of Just redBrchs @@ -405,8 +439,8 @@ where branchProc redSh = fmap (fmap $ go bw redSh) - reduce :: Shade x -> NonEmpty (DBranch' x [x]) - -> Maybe (NonEmpty (DBranch' x [x])) + reduce :: Shade x -> NonEmpty (DBranch' x [(x,y)]) + -> Maybe (NonEmpty (DBranch' x [(x,y)])) reduce sh@(Shade c _) brCandidates = case findIndex deficient cards of Just i | (DBranch _ reBr, o:ok) @@ -422,16 +456,17 @@ sShIdPartition' :: (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => Interior x -> [x] -> NonEmpty (DBranch' x [x])->NonEmpty (DBranch' x [x]) + => Interior x -> [(x,y)] -> NonEmpty (DBranch' x [(x,y)]) + -> NonEmpty (DBranch' x [(x,y)]) sShIdPartition' c xs st - = foldr (\p -> let (i,h) = ssi p + = foldr (\(p,y) -> let (i,h) = ssi p in asList $ update_nth (\(DBranch d c) - -> DBranch d (oneBulb h (p:) c)) + -> DBranch d (oneBulb h ((p,y):) c)) i ) st xs where ssi = subshadeId' (fromInterior c) (boughDirection<$>st) sShIdPartition :: (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => Shade x -> [x] -> NonEmpty (DBranch' x [x]) + => Shade x -> [(x,y)] -> NonEmpty (DBranch' x [(x,y)]) sShIdPartition (Shade c expa) xs | b:bs <- [DBranch v mempty | v <- normSpanningSystem' expa] = sShIdPartition' c xs $ b:|bs @@ -470,59 +505,48 @@ sortByKey = map snd . sortBy (comparing fst) -trunks :: ∀ x. (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => ShadeTree x -> [Shade x] +trunks :: ∀ x y . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) + => x`Shaded`y -> [Shade x] trunks t = case (pseudoAffineWitness :: PseudoAffineWitness x, t) of (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness), PlainLeaves lvs) - -> pointsCovers . catMaybes $ toInterior<$>lvs - (_, DisjointBranches _ brs) -> Hask.foldMap trunks brs - (_, OverlappingBranches _ sh _) -> [sh] + -> pointsCovers . catMaybes $ toInterior.fst<$>lvs + (_, DisjointBranches _ brs) -> Hask.foldMap trunks brs + (_, OverlappingBranches _ sh _) -> [sh] -nLeaves :: ShadeTree x -> Int +nLeaves :: x`Shaded`y -> Int nLeaves (PlainLeaves lvs) = length lvs nLeaves (DisjointBranches n _) = n nLeaves (OverlappingBranches n _ _) = n +treeDepth :: x`Shaded`y -> Int +treeDepth (PlainLeaves lvs) = 0 +treeDepth (DisjointBranches _ brs) = 1 + maximum (treeDepth<$>brs) +treeDepth (OverlappingBranches _ _ brs) + = 1 + maximum (maximum . fmap treeDepth<$>brs) -instance ImpliesMetric ShadeTree where - type MetricRequirement ShadeTree x = (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - inferMetric = stInfMet - where stInfMet :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => ShadeTree x -> Metric x - stInfMet (OverlappingBranches _ (Shade _ e) _) = dualNorm' e - stInfMet (PlainLeaves lvs) - = case pointsShades $ Hask.toList . toInterior =<< lvs :: [Shade x] of - (Shade _ sh:_) -> dualNorm' sh - _ -> mempty - stInfMet (DisjointBranches _ (br:|_)) = inferMetric br - inferMetric' = stInfMet - where stInfMet :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => ShadeTree x -> Metric' x - stInfMet (OverlappingBranches _ (Shade _ e) _) = e - stInfMet (PlainLeaves lvs) - = case pointsShades $ Hask.toList . toInterior =<< lvs :: [Shade x] of - (Shade _ sh:_) -> sh - _ -> mempty - stInfMet (DisjointBranches _ (br:|_)) = inferMetric' br -overlappingBranches :: Shade x -> NonEmpty (DBranch x) -> ShadeTree x + +overlappingBranches :: Shade x -> NonEmpty (DBranch x y) -> x`Shaded`y overlappingBranches shx brs = OverlappingBranches n shx brs where n = sum $ fmap (sum . fmap nLeaves) brs -unsafeFmapLeaves :: (x -> x) -> ShadeTree x -> ShadeTree x +unsafeFmapLeaves_ :: (x -> x) -> x`Shaded`y -> x`Shaded`y +unsafeFmapLeaves_ = unsafeFmapLeaves . first + +unsafeFmapLeaves :: ((x,y) -> (x,y')) -> x`Shaded`y -> x`Shaded`y' unsafeFmapLeaves f (PlainLeaves lvs) = PlainLeaves $ fmap f lvs unsafeFmapLeaves f (DisjointBranches n brs) - = DisjointBranches n $ unsafeFmapLeaves f <$> brs + = DisjointBranches n $ unsafeFmapLeaves f <$> brs unsafeFmapLeaves f (OverlappingBranches n sh brs) = OverlappingBranches n sh $ fmap (unsafeFmapLeaves f) <$> brs -unsafeFmapTree :: (NonEmpty x -> NonEmpty y) - -> (Needle' x -> Needle' y) - -> (Shade x -> Shade y) - -> ShadeTree x -> ShadeTree y +unsafeFmapTree :: (NonEmpty (x,y) -> NonEmpty (ξ,υ)) + -> (Needle' x -> Needle' ξ) + -> (Shade x -> Shade ξ) + -> x`Shaded`y -> ξ`Shaded`υ unsafeFmapTree _ _ _ (PlainLeaves []) = PlainLeaves [] unsafeFmapTree f _ _ (PlainLeaves lvs) = PlainLeaves . toList . f $ NE.fromList lvs unsafeFmapTree f fn fs (DisjointBranches n brs) @@ -534,22 +558,13 @@ ) brs in overlappingBranches (fs sh) brs' -coerceShadeTree :: ∀ x y . (LocallyCoercible x y, Manifold x, Manifold y) - => ShadeTree x -> ShadeTree y -coerceShadeTree = case ( dualSpaceWitness :: DualNeedleWitness x - , dualSpaceWitness :: DualNeedleWitness y ) of - (DualSpaceWitness,DualSpaceWitness) - -> unsafeFmapTree (fmap locallyTrivialDiffeomorphism) - (coerceNeedle' ([]::[(x,y)]) $) - coerceShade - -type Twig x = (Int, ShadeTree x) -type TwigEnviron x = [Twig x] +type Twig x y = (Int, x`Shaded`y) +type TwigEnviron x y = [Twig x y] -allTwigs :: ∀ x . WithField ℝ PseudoAffine x => ShadeTree x -> [Twig x] +allTwigs :: ∀ x y . WithField ℝ PseudoAffine x => x`Shaded`y -> [Twig x y] allTwigs tree = go 0 tree [] where go n₀ (DisjointBranches _ dp) = snd (foldl' (\(n₀',prev) br -> (n₀'+nLeaves br, prev . go n₀' br)) (n₀,id) dp) @@ -567,15 +582,15 @@ -- | Example: https://nbviewer.jupyter.org/github/leftaroundabout/manifolds/blob/master/test/Trees-and-Webs.ipynb#pseudorandomCloudTree -- -- <> -twigsWithEnvirons :: ∀ x. (WithField ℝ Manifold x, SimpleSpace (Needle x)) - => ShadeTree x -> [(Twig x, TwigEnviron x)] +twigsWithEnvirons :: ∀ x y. (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => x`Shaded`y -> [(Twig x y, TwigEnviron x y)] twigsWithEnvirons = execWriter . traverseTwigsWithEnvirons (writer . (snd.fst&&&pure)) -traverseTwigsWithEnvirons :: ∀ x f . +traverseTwigsWithEnvirons :: ∀ x y f . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x), Hask.Applicative f) - => ( (Twig x, TwigEnviron x) -> f (ShadeTree x) ) -> ShadeTree x -> f (ShadeTree x) + => ( (Twig x y, TwigEnviron x y) -> f (x`Shaded`y) ) -> x`Shaded`y -> f (x`Shaded`y) traverseTwigsWithEnvirons f = fst . go pseudoAffineWitness [] . (0,) - where go :: PseudoAffineWitness x -> TwigEnviron x -> Twig x -> (f (ShadeTree x), Bool) + where go :: PseudoAffineWitness x -> TwigEnviron x y -> Twig x y -> (f (x`Shaded`y), Bool) go sw _ (i₀, DisjointBranches nlvs djbs) = ( fmap (DisjointBranches nlvs) . Hask.traverse (fst . go sw []) $ NE.zip ioffs djbs @@ -612,7 +627,7 @@ go (PseudoAffineWitness (SemimanifoldWitness _)) envi plvs@(i₀, (PlainLeaves _)) = (f $ purgeRemotes (plvs, envi), True) - twigProximæ :: PseudoAffineWitness x -> Interior x -> ShadeTree x -> TwigEnviron x + twigProximæ :: PseudoAffineWitness x -> Interior x -> x`Shaded`y -> TwigEnviron x y twigProximæ sw x₀ (DisjointBranches _ djbs) = Hask.foldMap (\(i₀,st) -> first (+i₀) <$> twigProximæ sw x₀ st) $ NE.zip ioffs djbs @@ -627,7 +642,7 @@ where overlap = bdir<.>^δxb twigProximæ _ _ plainLeaves = [(0, plainLeaves)] - twigsaveTrim :: (DBranch x -> TwigEnviron x) -> ShadeTree x -> TwigEnviron x + twigsaveTrim :: (DBranch x y -> TwigEnviron x y) -> x`Shaded`y -> TwigEnviron x y twigsaveTrim f ct@(OverlappingBranches _ _ dbs) = case Hask.mapM (\(i₀,dbr) -> noLeaf $ first(+i₀)<$>f dbr) $ NE.zip ioffs dbs of @@ -637,7 +652,7 @@ noLeaf bqs = pure bqs ioffs = NE.scanl (\i -> (+i) . sum . fmap nLeaves . toList) 0 dbs - purgeRemotes :: (Twig x, TwigEnviron x) -> (Twig x, TwigEnviron x) + purgeRemotes :: (Twig x y, TwigEnviron x y) -> (Twig x y, TwigEnviron x y) purgeRemotes = id -- See 7d1f3a4 for the implementation; this didn't work reliable. completeTopShading :: ∀ x y . ( WithField ℝ PseudoAffine x, WithField ℝ PseudoAffine y @@ -647,14 +662,14 @@ , dualSpaceWitness :: DualNeedleWitness y ) of (DualSpaceWitness, DualSpaceWitness) -> pointsShade's . catMaybes - $ toInterior . (_topological &&& _untopological) <$> plvs + $ toInterior <$> plvs completeTopShading (DisjointBranches _ bqs) = take 1 . completeTopShading =<< NE.toList bqs completeTopShading t = case ( dualSpaceWitness :: DualNeedleWitness x , dualSpaceWitness :: DualNeedleWitness y ) of (DualSpaceWitness, DualSpaceWitness) -> pointsCover's . catMaybes - . map (toInterior <<< _topological &&& _untopological) $ onlyLeaves t + . map toInterior $ onlyLeaves t transferAsNormsDo :: ∀ v . LSpace v => Norm v -> Variance v -> v-+>v @@ -682,8 +697,8 @@ fts (xc, (Shade' yc expay, jtg)) = unsafeFmapLeaves applδj t where Just δyc = yc.-~.yc₀ tfm = transferAsNormsDo expay₀ (dualNorm expay) - applδj (WithAny y x) - = WithAny (yc₀ .+~^ ((tfm $ δy) ^+^ (jtg $ δx) ^+^ δyc)) x + applδj (x,y) + = (x, yc₀ .+~^ ((tfm $ δy) ^+^ (jtg $ δx) ^+^ δyc)) where Just δx = x.-~.xc Just δy = y.-~.(yc₀.+~^(j₀ $ δx)) @@ -705,24 +720,26 @@ -seekPotentialNeighbours :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => ShadeTree x -> x`Shaded`[Int] +seekPotentialNeighbours :: ∀ x y . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) + => x`Shaded`y -> x`Shaded`(y,[Int]) seekPotentialNeighbours tree = zipTreeWithList tree - $ snd<$>leavesWithPotentialNeighbours tree + $ case snd<$>leavesWithPotentialNeighbours tree of + (n:ns) -> n:|ns -leavesWithPotentialNeighbours :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => ShadeTree x -> [(x, [Int])] +leavesWithPotentialNeighbours :: ∀ x y + . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) + => x`Shaded`y -> [((x,y), [Int])] leavesWithPotentialNeighbours = map (second snd) . go pseudoAffineWitness 0 0 [] - where go :: PseudoAffineWitness x -> Depth -> Int -> [Wall x] -> ShadeTree x - -> [(x, ([Wall x], [Int]))] + where go :: PseudoAffineWitness x -> Depth -> Int -> [Wall x] -> x`Shaded`y + -> [((x,y), ([Wall x], [Int]))] go (PseudoAffineWitness (SemimanifoldWitness _)) depth n₀ walls (PlainLeaves lvs) - = [ (x, ( [ wall & wallDistance .~ d + = [ ((x,y), ( [ wall & wallDistance .~ d | wall <- walls , Just vw <- [toInterior x>>=(.-~.wall^.wallAnchor)] , let d = (wall^.wallNormal)<.>^vw , d < wall^.wallDistance ] , [] )) - | x <- lvs ] + | (x,y) <- lvs ] go pw depth n₀ walls (DisjointBranches _ dp) = snd (foldl' (\(n₀',prev) br -> ( n₀'+nLeaves br , prev . (go pw depth n₀' walls br++))) @@ -731,9 +748,9 @@ depth n₀ walls (OverlappingBranches _ (Shade brCtr _) dp) = reassemble $ snd (foldl' assignWalls (n₀,id) . directionIChoices 0 $ NE.toList dp) [] - where assignWalls :: (Int, DList (x, ([Wall x],[Int]))) - -> ((Int,(Needle' x, ShadeTree x)), [(Int,(Needle' x, ShadeTree x))]) - -> (Int, DList (x, ([Wall x], [Int]))) + where assignWalls :: (Int, DList ((x,y), ([Wall x],[Int]))) + -> ((Int,(Needle' x, x`Shaded`y)), [(Int,(Needle' x, x`Shaded`y))]) + -> (Int, DList ((x,y), ([Wall x], [Int]))) assignWalls (n₀',prev) ((iDir,(thisDir,br)),otherDirs) = ( n₀'+nLeaves br , prev . (go pw (depth+1) n₀' @@ -747,7 +764,7 @@ updWall wall = wall & wallDistance %~ min bcDist where Just vbw = brCtr.-~.wall^.wallAnchor bcDist = (wall^.wallNormal)<.>^vbw - reassemble :: [(x, ([Wall x],[Int]))] -> [(x, ([Wall x],[Int]))] + reassemble :: [((x,y), ([Wall x],[Int]))] -> [((x,y), ([Wall x],[Int]))] reassemble pts = [ (x, (higherWalls, newGroups++deeperGroups)) | (x, (allWalls, deeperGroups)) <- pts , let (levelWalls,higherWalls) @@ -902,6 +919,8 @@ -- 'NonEmptyTree' x ≅ (x, 'Trees' x) -- @ type NonEmptyTree = GenericTree NonEmpty [] + +type LeafyTree x y = GenericTree [] (ListT (Either y)) x newtype GenericTree c b x = GenericTree { treeBranches :: c (x,GenericTree b b x) } deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable) @@ -912,22 +931,49 @@ instance (Hask.MonadPlus c) => Monoid (GenericTree c b x) where mempty = GenericTree Hask.mzero mappend = (<>) -deriving instance Show (c (x, GenericTree b b x)) => Show (GenericTree c b x) +instance Show (c (x, GenericTree b b x)) => Show (GenericTree c b x) where + showsPrec p (GenericTree t) = showParen (p>9) $ ('朳':) . showsPrec 10 t +deriving instance Eq (c (x, GenericTree b b x)) => Eq (GenericTree c b x) +-- | @U+6733 CJK UNIFIED IDEOGRAPH tree@. +-- The main purpose of this is to give 'GenericTree' a more concise 'Show' instance. +朳 :: c (x, GenericTree b b x) -> GenericTree c b x +朳 = GenericTree + -- | Imitate the specialised 'ShadeTree' structure with a simpler, generic tree. onlyNodes :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) => ShadeTree x -> Trees x onlyNodes (PlainLeaves []) = GenericTree [] onlyNodes (PlainLeaves ps) = let (ctr,_) = pseudoECM ([]::[x]) $ NE.fromList ps - in GenericTree [ (ctr, GenericTree $ (,mempty) <$> ps) ] + in GenericTree [ (ctr, GenericTree $ (,mempty).fst <$> ps) ] onlyNodes (DisjointBranches _ brs) = Hask.foldMap onlyNodes brs onlyNodes (OverlappingBranches _ (Shade ctr _) brs) = GenericTree [ ( fromInterior ctr , Hask.foldMap (Hask.foldMap onlyNodes) brs ) ] +entireTree :: ∀ x y . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) + => x`Shaded`y -> LeafyTree x y +entireTree (PlainLeaves lvs) + = let (ctr,_) = pseudoECM ([]::[x]) $ NE.fromList lvs + in GenericTree [ (ctr, GenericTree . ListT $ Right + [ (x, GenericTree . lift $ Left y) + | (x,y)<-lvs ] ) + ] +entireTree (DisjointBranches _ brs) + = GenericTree [ (x, GenericTree subt) + | GenericTree sub <- NE.toList $ fmap entireTree brs + , (x, GenericTree subt) <- sub ] +entireTree (OverlappingBranches _ (Shade ctr _) brs) + = GenericTree [ ( fromInterior ctr + , GenericTree . ListT . Right + $ Hask.foldMap (Hask.foldMap $ treeBranches . entireTree) brs ) ] + -- | Left (and, typically, also right) inverse of 'fromLeafNodes'. -onlyLeaves :: WithField ℝ PseudoAffine x => ShadeTree x -> [x] +onlyLeaves_ :: WithField ℝ PseudoAffine x => ShadeTree x -> [x] +onlyLeaves_ = map fst . onlyLeaves + +onlyLeaves :: WithField ℝ PseudoAffine x => x`Shaded`y -> [(x,y)] onlyLeaves tree = dismantle tree [] where dismantle (PlainLeaves xs) = (xs++) dismantle (OverlappingBranches _ _ brs) @@ -970,33 +1016,23 @@ -constShaded :: y -> ShadeTree x -> x`Shaded`y -constShaded y = unsafeFmapTree (WithAny y<$>) id (shadeWithAny y) - -stripShadedUntopological :: x`Shaded`y -> ShadeTree x -stripShadedUntopological = unsafeFmapTree (fmap _topological) id shadeWithoutAnything - -fmapShaded :: (y -> υ) -> (x`Shaded`y) -> (x`Shaded`υ) -fmapShaded f = unsafeFmapTree (fmap $ \(WithAny y x) -> WithAny (f y) x) - id - (\(Shade yx shx) -> Shade (fmap f yx) shx) +constShaded :: y -> x`Shaded`y₀ -> x`Shaded`y +constShaded y = unsafeFmapTree (fmap . second $ const y) id id -joinShaded :: (x`WithAny`y)`Shaded`z -> x`Shaded`(y,z) -joinShaded = unsafeFmapTree (fmap $ \(WithAny z (WithAny y x)) -> WithAny (y,z) x) - id - (\(Shade (WithAny z (WithAny y x)) shx) - -> Shade (WithAny (y,z) x) shx ) +fmapShaded :: (Semimanifold x, SimpleSpace (Needle x)) + => (y -> υ) -> (x`Shaded`y) -> (x`Shaded`υ) +fmapShaded f = unsafeFmapTree (fmap $ second f) id id -zipTreeWithList :: ShadeTree x -> [y] -> (x`Shaded`y) -zipTreeWithList tree = go tree . cycle - where go (PlainLeaves lvs) ys = PlainLeaves $ zipWith WithAny ys lvs +zipTreeWithList :: x`Shaded`w -> NonEmpty y -> (x`Shaded`(w,y)) +zipTreeWithList tree = go tree . NE.toList . NE.cycle + where go (PlainLeaves lvs) ys = PlainLeaves $ zipWith (\(x,w) y -> (x,(w,y))) lvs ys go (DisjointBranches n brs) ys = DisjointBranches n . NE.fromList $ snd (foldl (\(ys',prev) br -> (drop (nLeaves br) ys', prev . (go br ys':)) ) (ys,id) $ NE.toList brs) [] - go (OverlappingBranches n (Shade xoc shx) brs) ys - = OverlappingBranches n (Shade (WithAny (head ys) xoc) shx) . NE.fromList + go (OverlappingBranches n shx brs) ys + = OverlappingBranches n shx . NE.fromList $ snd (foldl (\(ys',prev) (DBranch dir (Hourglass top bot)) -> case drop (nLeaves top) ys' of ys'' -> ( drop (nLeaves bot) ys'' @@ -1005,17 +1041,14 @@ ) ) (ys,id) $ NE.toList brs) [] --- | This is to 'ShadeTree' as 'Data.Map.Map' is to 'Data.Set.Set'. -type x`Shaded`y = ShadeTree (x`WithAny`y) - stiWithDensity :: ∀ x y . ( WithField ℝ PseudoAffine x, LinearSpace y, Scalar y ~ ℝ , SimpleSpace (Needle x) ) => x`Shaded`y -> x -> Cℝay y stiWithDensity (PlainLeaves lvs) | [Shade baryc expa :: Shade x] <- pointsShades . catMaybes - $ toInterior . _topological <$> lvs + $ toInterior . fst <$> lvs = let nlvs = fromIntegral $ length lvs :: ℝ - indiShapes = [(Shade pi expa, y) | WithAny y p <- lvs + indiShapes = [(Shade pi expa, y) | (p,y) <- lvs , Just pi <- [toInterior p]] in \x -> let lcCoeffs = [ occlusion psh x | (psh, _) <- indiShapes ] dens = sum lcCoeffs @@ -1025,7 +1058,7 @@ = \x -> foldr1 qGather $ (`stiWithDensity`x)<$>lvs where qGather (Cℝay 0 _) o = o qGather o _ = o -stiWithDensity (OverlappingBranches n (Shade (WithAny _ bc) extend) brs) +stiWithDensity (OverlappingBranches n (Shade bc extend) brs) = ovbSWD (dualSpaceWitness, pseudoAffineWitness) where ovbSWD :: (DualNeedleWitness x, PseudoAffineWitness x) -> x -> Cℝay y ovbSWD (DualSpaceWitness, PseudoAffineWitness (SemimanifoldWitness _)) x @@ -1051,35 +1084,19 @@ -> ( xloc, ( (yloc, recip $ shd|$|(0,1)) , dependence (dualNorm shd) ) ) -smoothInterpolate :: ∀ x y . ( WithField ℝ Manifold x, LinearSpace y, Scalar y ~ ℝ - , SimpleSpace (Needle x) ) - => NonEmpty (x,y) -> x -> y -smoothInterpolate = si boundarylessWitness - where si :: BoundarylessWitness x -> NonEmpty (x,y) -> x -> y - si BoundarylessWitness l = \x -> - case ltr x of - Cℝay 0 _ -> defy - Cℝay _ y -> y - where defy = linearCombo [(y, 1/n) | WithAny y _ <- l'] - n = fromIntegral $ length l' - l' = (uncurry WithAny . swap) <$> NE.toList l - ltr = stiWithDensity $ fromLeafPoints l' - spanShading :: ∀ x y . ( WithField ℝ Manifold x, WithField ℝ Manifold y , SimpleSpace (Needle x), SimpleSpace (Needle y) ) => (Shade x -> Shade y) -> ShadeTree x -> x`Shaded`y -spanShading f = unsafeFmapTree addYs id addYSh - where addYs :: NonEmpty x -> NonEmpty (x`WithAny`y) - addYs l = foldr (NE.<|) (fmap (WithAny $ fromInterior ymid) l ) - (fmap (`WithAny` fromInterior xmid) yexamp) +spanShading f = unsafeFmapTree (addYs . fmap fst) id id + where addYs :: NonEmpty x -> NonEmpty (x,y) + addYs l = foldr (NE.<|) (fmap (,fromInterior ymid) l ) + (fmap (fromInterior xmid,) yexamp) where [xsh@(Shade xmid _)] = pointsCovers . catMaybes . toList $ toInterior<$>l Shade ymid yexpa = f xsh yexamp = [ ymid .+~^ σ*^δy | δy <- varianceSpanningSystem yexpa, σ <- [-1,1] ] - addYSh :: Shade x -> Shade (x`WithAny`y) - addYSh xsh = shadeWithAny (fromInterior . _shadeCtr $ f xsh) xsh diff --git a/Data/Manifold/Types.hs b/Data/Manifold/Types.hs --- a/Data/Manifold/Types.hs +++ b/Data/Manifold/Types.hs @@ -55,7 +55,7 @@ -- ** Lines , Line(..), lineAsPlaneIntersection -- ** Hyperplanes - , Cutplane(..) + , Cutplane(..), normalPlane , fathomCutDistance, sideOfCut, cutPosBetween -- * Linear mappings , LinearMap, LocalLinear @@ -319,7 +319,12 @@ -- behave locally as a plane, globally as an (/n/−1)-dimensional submanifold. data Cutplane x = Cutplane { sawHandle :: x , cutNormal :: Stiefel1 (Needle x) } +deriving instance (Show x, Show (Needle' x)) => Show (Cutplane x) +normalPlane :: x -- ^ Some point lying in the desired plane. + -> Needle' x -- ^ Co-vector perpendicular to the plane. Must be nonzero. + -> Cutplane x +normalPlane x n = Cutplane x $ Stiefel1 n sideOfCut :: (WithField ℝ PseudoAffine x, LinearSpace (Needle x)) diff --git a/Data/Manifold/Types/Stiefel.hs b/Data/Manifold/Types/Stiefel.hs --- a/Data/Manifold/Types/Stiefel.hs +++ b/Data/Manifold/Types/Stiefel.hs @@ -17,6 +17,7 @@ -- scalings, and we prefer that definition since it doesn't require a notion of -- unit length (which is only defined in inner-product spaces). +{-# LANGUAGE StandaloneDeriving, FlexibleContexts, UndecidableInstances #-} @@ -45,3 +46,4 @@ newtype Stiefel1 v = Stiefel1 { getStiefel1N :: DualVector v } +deriving instance (Show (DualVector v)) => Show (Stiefel1 v) diff --git a/Data/Manifold/Web.hs b/Data/Manifold/Web.hs --- a/Data/Manifold/Web.hs +++ b/Data/Manifold/Web.hs @@ -38,7 +38,7 @@ -- ** Construction , fromWebNodes, fromShadeTree_auto, fromShadeTree, fromShaded -- ** Lookup - , nearestNeighbour, indexWeb, webEdges, toGraph + , nearestNeighbour, indexWeb, toGraph, webBoundary -- ** Decomposition , sliceWeb_lin -- , sampleWebAlongLine_lin , sampleWeb_2Dcartesian_lin, sampleEntireWeb_2Dcartesian_lin @@ -46,25 +46,25 @@ , localFocusWeb -- * Uncertain functions , differentiateUncertainWebFunction, differentiate²UncertainWebFunction + , localModels_CGrid -- * Differential equations -- ** Fixed resolution - , iterateFilterDEqn_static + , iterateFilterDEqn_static, iterateFilterDEqn_static_selective -- ** Automatic resolution , filterDEqnSolutions_adaptive, iterateFilterDEqn_adaptive -- ** Configuration , InconsistencyStrategy(..) , InformationMergeStrategy(..) - , naïve, inconsistencyAware, indicateInconsistencies + , naïve, inconsistencyAware, indicateInconsistencies, postponeInconsistencies , PropagationInconsistency(..) -- * Misc , ConvexSet(..), ellipsoid, ellipsoidSet, coerceWebDomain - , rescanPDEOnWeb, rescanPDELocally, webOnions + , rescanPDELocally, webOnions, knitShortcuts ) where import Data.List hiding (filter, all, foldr1) import Data.Maybe -import qualified Data.Set as Set import qualified Data.Map as Map import qualified Data.Vector as Arr import qualified Data.Vector.Mutable as MArr @@ -77,7 +77,7 @@ import Control.DeepSeq import Data.VectorSpace -import Math.LinearMap.Category +import Math.LinearMap.Category hiding (trace) import Data.Tagged import Data.Function (on) @@ -91,8 +91,10 @@ import Data.SetLike.Intersection import Data.Manifold.Riemannian import Data.Manifold.Atlas +import Data.Manifold.Function.LocalModel import Data.Manifold.Function.Quadratic import Data.Function.Affine +import Data.Manifold.Web.Internal import Data.Embedding import qualified Prelude as Hask hiding(foldl, sum, sequence) @@ -103,6 +105,7 @@ import Control.Monad.Trans.State import Control.Monad.Trans.List import Control.Monad.Trans.Except +import Control.Monad.Trans.Writer hiding (censor) import Data.Functor.Identity (Identity(..)) import qualified Data.Foldable as Hask import Data.Foldable (all, toList) @@ -119,100 +122,28 @@ import Control.Comonad (Comonad(..)) import Control.Comonad.Cofree -import Control.Lens ((&), (%~), (^.), (.~), (+~)) +import Control.Lens ((&), (%~), (^.), (.~), (+~), ix) import Control.Lens.TH import GHC.Generics (Generic) - -type WebNodeId = Int - -data Neighbourhood x = Neighbourhood { - _neighbours :: UArr.Vector WebNodeId - , _localScalarProduct :: Metric x - } - deriving (Generic) -makeLenses ''Neighbourhood - -deriving instance ( WithField ℝ PseudoAffine x - , SimpleSpace (Needle x), Show (Needle' x) ) - => Show (Neighbourhood x) - -data WebLocally x y = LocalWebInfo { - _thisNodeCoord :: x - , _thisNodeData :: y - , _thisNodeId :: WebNodeId - , _nodeNeighbours :: [(WebNodeId, (Needle x, WebLocally x y))] - , _nodeLocalScalarProduct :: Metric x - , _nodeIsOnBoundary :: Bool - } deriving (Generic) -makeLenses ''WebLocally - -data NeighbourhoodVector x = NeighbourhoodVector - { _nvectId :: Int - , _theNVect :: Needle x - , _nvectNormal :: Needle' x - , _nvectLength :: Scalar (Needle x) - , _otherNeighboursOverlap :: Scalar (Needle x) - } -makeLenses ''NeighbourhoodVector - -data PropagationInconsistency x υ = PropagationInconsistency { - _inconsistentPropagatedData :: [(x,υ)] - , _inconsistentAPrioriData :: υ } - | PropagationInconsistencies [PropagationInconsistency x υ] - deriving (Show) -makeLenses ''PropagationInconsistency - -instance Monoid (PropagationInconsistency x υ) where - mempty = PropagationInconsistencies [] - mappend p q = mconcat [p,q] - mconcat = PropagationInconsistencies - -instance (NFData x, NFData (Metric x)) => NFData (Neighbourhood x) - --- | A 'PointsWeb' is almost, but not quite a mesh. It is a stongly connected† --- directed graph, backed by a tree for fast nearest-neighbour lookup of points. --- --- †In general, there can be disconnected components, but every connected --- component is strongly connected. -data PointsWeb :: * -> * -> * where - PointsWeb :: { - webNodeRsc :: ShadeTree x - , webNodeAssocData :: Arr.Vector (y, Neighbourhood x) - } -> PointsWeb x y - deriving (Generic, Hask.Functor, Hask.Foldable, Hask.Traversable) - -instance (NFData x, NFData (Metric x), NFData (Needle' x), NFData y) => NFData (PointsWeb x y) - -instance Foldable (PointsWeb x) (->) (->) where - ffoldl = uncurry . Hask.foldl' . curry - foldMap = Hask.foldMap -instance Traversable (PointsWeb x) (PointsWeb x) (->) (->) where - traverse f (PointsWeb rsc asd) - = fmap (PointsWeb rsc . (`Arr.zip`ngss) . Arr.fromList) - . traverse f $ Arr.toList ys - where (ys,ngss) = Arr.unzip asd - - - -type MetricChoice x = Shade x -> Metric x +import Development.Placeholders fromWebNodes :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (MetricChoice x) -> [(x,y)] -> PointsWeb x y fromWebNodes = case boundarylessWitness :: BoundarylessWitness x of BoundarylessWitness -> - \mf -> fromShaded mf . fromLeafPoints . map (uncurry WithAny . swap) + \mf -> fromShaded mf . fromLeafPoints_ fromTopWebNodes :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (MetricChoice x) -> [((x,[Int+Needle x]),y)] -> PointsWeb x y fromTopWebNodes = case boundarylessWitness :: BoundarylessWitness x of BoundarylessWitness -> - \mf -> fromTopShaded mf . fromLeafPoints - . map (uncurry WithAny . swap . regroup') + \mf -> fromTopShaded mf . fromLeafPoints_ . map regroup' -fromShadeTree_auto :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => ShadeTree x -> PointsWeb x () +fromShadeTree_auto :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => ShadeTree x -> PointsWeb x () fromShadeTree_auto = fromShaded (dualNorm' . _shadeExpanse) . constShaded () fromShadeTree :: ∀ x . (WithField ℝ Manifold x, SimpleSpace (Needle x)) @@ -226,192 +157,224 @@ -- Riemannian metric). -> (x`Shaded`y) -- ^ Source tree. -> PointsWeb x y -fromShaded metricf = smoothenWebTopology metricf - . fromTopShaded metricf . fmapShaded (first (map Left) . swap) - . joinShaded . seekPotentialNeighbours +fromShaded metricf = knitShortcuts metricf . autoLinkWeb . unlinkedFromShaded metricf toShaded :: WithField ℝ PseudoAffine x => PointsWeb x y -> (x`Shaded`y) -toShaded (PointsWeb shd asd) = zipTreeWithList shd $ Arr.toList (fst<$>asd) +toShaded (PointsWeb shd) = fmap _dataAtNode shd +unlinkedFromShaded :: ∀ x y . SimpleSpace (Needle x) + => MetricChoice x -> (x`Shaded`y) -> PointsWeb x y +unlinkedFromShaded metricf = PointsWeb<< Neighbourhood y mempty nm (Just dv) + where nm = metricf $notImplemented + dv = head $ normSpanningSystem nm + + + +autoLinkWeb :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => PointsWeb x y -> PointsWeb x y +autoLinkWeb = runIdentity . traverseNodesInEnvi ( pure . fetchNgbs [] + . (id &&& findEnviPts (0,1)) ) + where fetchNgbs :: [(WebNodeIdOffset, Needle x)] + -> (NodeInWeb x y, [[(WebNodeIdOffset, (x, Neighbourhood x y))]]) + -> Neighbourhood x y + fetchNgbs alreadyFound + ( NodeInWeb (x, Neighbourhood y aprNgbs locMetr (Just wall)) + layersAroundThis + , enviLayers ) + | (δi, (v, nh)) : _ <- newNgbCandidates + = fetchNgbs + ((δi, v) : alreadyFound) + ( NodeInWeb (x, Neighbourhood y (UArr.cons δi aprNgbs) locMetr + $ if dimension > 1 + then pumpHalfspace locMetr v + (wall, snd<$>alreadyFound) + else case alreadyFound of + [] -> Just $ locMetr<$|v + [_] -> Nothing + ) + layersAroundThis + , enviLayers ) + where newNgbCandidates + = [ (δi, (v, nh)) + | envi <- enviLayers + , (δi, ((v,_), nh)) <- sortBy (comparing $ snd . fst . snd) + [ (δi, ((v, if dimension > 1 + then gatherDirectionsBadness + $ linkingUndesirability distSq wallDist + else distSq + ), nh)) + | (δi,(xp,nh)) <- envi + , let Just v = xp.-~.x + distSq = normSq locMetr v + wallDist = walln<.>^v + , wallDist >= 0 + , distSq > wallDist^2 + || dimension==1 -- in 1D, we must allow linking + -- to the direct opposite + -- (there IS no other direction) + , not . any (==δi) $ UArr.toList aprNgbs + ++ map fst alreadyFound + ] ] + locMetr' = dualNorm locMetr + walln = wall ^/ (- (locMetr'|$|wall)) + fetchNgbs _ (NodeInWeb (_, d) _, _) = d + findEnviPts (iw,wedgeSize) (NodeInWeb tr ((envi,iSpl):envis)) + = (zip [-iw-iSpl ..] preds ++ zip [wedgeSize-iw ..] succs) + : findEnviPts (iw+iSpl, wedgeSize + iSpl + length succs) + (NodeInWeb tr envis) + where (preds, succs) = splitAt iSpl $ onlyLeaves envi + findEnviPts _ _ = [] + dimension = subbasisDimension (entireBasis :: SubBasis (Needle x)) + fromTopShaded :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (MetricChoice x) -> (x`Shaded`([Int+Needle x], y)) -- ^ Source tree, with topology information -- (IDs of neighbour-candidates, or needles pointing to them) -> PointsWeb x y -fromTopShaded metricf shd = PointsWeb shd' assocData - where shd' = stripShadedUntopological shd - assocData = Hask.foldMap locMesh $ allTwigs shd - - locMesh :: (Int, ShadeTree (x`WithAny`([Int+Needle x], y))) - -> Arr.Vector (y, Neighbourhood x) - locMesh (i₀, locT) = Arr.map findNeighbours $ Arr.fromList locLeaves - where locLeaves :: [ (Int, x`WithAny`([Int+Needle x], y)) ] - locLeaves = map (first (+i₀)) . zip [0..] $ onlyLeaves locT - findNeighbours :: (Int, x`WithAny`([Int+Needle x], y)) -> (y, Neighbourhood x) - findNeighbours (i, WithAny (vns,y) x) - = (y, cullNeighbours locRieM - (i, WithAny([ (i,v) - | (i,WithAny _ xN) <- locLeaves - , Just v <- [xN.-~.x] ] - ++ aprioriNgbs) - x)) - where aprioriNgbs :: [(Int, Needle x)] - aprioriNgbs = catMaybes - [ (second $ const v) <$> - positionIndex (pure locRieM) shd' xN - | Right v <- vns - , let xN = xi.+~^v :: x ] - ++ [ (i,v) | Left i <- vns - , Right (_,xN) <- [indexShadeTree shd' i] - , Just v <- [xN.-~.x] ] - Just xi = toInterior x - - locRieM :: Metric x - locRieM = case pointsCovers . catMaybes . map (toInterior . _topological) - $ onlyLeaves locT of - [sh₀] -> metricf sh₀ +fromTopShaded metricf shd = $notImplemented -cullNeighbours :: ∀ x . (WithField ℝ PseudoAffine x, SimpleSpace (Needle x)) - => Metric x -> (Int, x`WithAny`[(Int,Needle x)]) -> Neighbourhood x -cullNeighbours locRieM (i, WithAny vns x) - = Neighbourhood (UArr.fromList . sort $ _nvectId<$>execState seek mempty) - locRieM - where seek :: State [NeighbourhoodVector x] () - seek = do - Hask.forM_ ( fastNubBy (comparing fst) $ vns ) - $ \(iNgb, v) -> - when (iNgb/=i) `id`do - oldNgbs <- get - let w₀ = locRieM<$|v - l = sqrt $ w₀<.>^v - onOverlap = sum [ o^2 | nw<-oldNgbs - , let o = (nw^.nvectNormal)<.>^v - , o > 0 ] - when (l > onOverlap) `id`do - let w = w₀^/sqrt l^3 - newCandidates - = NeighbourhoodVector iNgb v w l 0 - : [ ongb & otherNeighboursOverlap .~ 0 - | ongb <- oldNgbs - , let o = w<.>^(ongb^.theNVect) - newOverlap = (if o > 0 then (o^2+) else id) - $ ongb^.otherNeighboursOverlap - , newOverlap < ongb^.nvectLength ] - put $ recalcOverlaps newCandidates - recalcOverlaps [] = [] - recalcOverlaps (ngb:ngbs) - = (ngb & otherNeighboursOverlap +~ furtherOvl) - : recalcOverlaps [ ngb' & otherNeighboursOverlap +~ max 0 o ^ 2 - | ngb' <- ngbs - , let o = (ngb^.nvectNormal)<.>^(ngb'^.theNVect) ] - where furtherOvl = sum [ o^2 | nw<-ngbs - , let o = (nw^.nvectNormal)<.>^(ngb^.theNVect) - , o > 0 ] - -- | Re-calculate the links in a web, so as to give each point a satisfyingly -- “complete-spanning” environment. smoothenWebTopology :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => MetricChoice x -> PointsWeb x y -> PointsWeb x y -smoothenWebTopology mc = swt - where swt (PointsWeb shd net) = PointsWeb shd . go allNodes Set.empty - . fst $ makeIndexLinksSymmetric net - where allNodes = Set.fromList . Arr.toList $ fst <$> Arr.indexed net - go activeSet pastLinks asd - | all (isNothing.fst) refined - , Set.null (Set.difference symmetryTouched pastLinks) - = Arr.imap finalise asd' - | otherwise = go (Set.fromList - [ j | (Just i, (_,Neighbourhood ngbs' _)) - <-refined - , j <- i : UArr.toList ngbs' ] - `Set.union` (Set.map fst symmetryTouched)) - updtLinks - asd' - where refined = reseek<$>Set.toList activeSet - where reseek i = ( guard isNews >> pure i - , (y, Neighbourhood newNgbs locRieM) ) - where isNews = newNgbs /= oldNgbs - && or [ not $ Set.member (i,j) pastLinks - | j <- UArr.toList newNgbs ] - (y,Neighbourhood oldNgbs locRieM) = asd Arr.! i - nextNeighbours = fastNub - $ UArr.toList oldNgbs - ++ (UArr.toList._neighbours.snd.(asd Arr.!) - =<< UArr.toList oldNgbs) - x = xLookup Arr.! i - Neighbourhood newNgbs _ - = cullNeighbours locRieM - ( i, WithAny [ (j,v) - | j <- nextNeighbours - , Just v - <- [x .-~. xLookup Arr.! j] ] - x ) - (asd', symmetryTouched) = makeIndexLinksSymmetric - $ asd Arr.// [(i,n) | (Just i,n) <- refined] - updtLinks = Set.unions - [ pastLinks - , Set.fromList - [ (i,j) | (Just i,(_,Neighbourhood n _)) <- refined - , j<-UArr.toList n ] - , symmetryTouched ] - finalise i (y, Neighbourhood n em) - = (y, cullNeighbours em (i, WithAny [ (j,v) - | j<-UArr.toList n - , let xN = xLookup Arr.! j - , Just v <- [xN.-~.x] ] - x )) - where x = xLookup Arr.! i - xLookup = Arr.fromList $ onlyLeaves shd +smoothenWebTopology = knitShortcuts -makeIndexLinksSymmetric - :: Arr.Vector (y, Neighbourhood x) - -> (Arr.Vector (y, Neighbourhood x), Set.Set (WebNodeId,WebNodeId)) -makeIndexLinksSymmetric orig = runST (do - result <- Arr.thaw orig - touched <- newSTRef $ Set.empty - (`Arr.imapM_`orig) $ \i (_,Neighbourhood ngbs _) -> do - UArr.forM_ ngbs $ \j -> do - (yn, Neighbourhood nngbs lsc) <- MArr.read result j - when (not $ i`UArr.elem`nngbs) `id`do - MArr.write result j (yn, Neighbourhood (UArr.snoc nngbs i) lsc) - modifySTRef touched $ Set.insert (j,i) - final <- Arr.freeze result - allTouched <- readSTRef touched - return (final, allTouched) - ) -indexWeb :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) - => PointsWeb x y -> WebNodeId -> Maybe (x,y) -indexWeb (PointsWeb rsc assocD) i - | i>=0, i WebNodeId -> y -unsafeIndexWebData (PointsWeb _ asd) i = fst (asd Arr.! i) +type OSNeedle x = (Needle' x, Needle x) +type OSNode x y = (OSNeedle x, WebLocally x y) +type CPCone x = (Needle' x, ℝ) -webEdges :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) - => PointsWeb x y -> [((x,y), (x,y))] -webEdges web@(PointsWeb rsc assoc) = (lookId***lookId) <$> toList allEdges - where allEdges :: Set.Set (WebNodeId,WebNodeId) - allEdges = Hask.foldMap (\(i,(_, Neighbourhood ngbs _)) - -> Set.fromList [(min i i', max i i') - | i'<-UArr.toList ngbs ] - ) $ Arr.indexed assoc - lookId i | Just xy <- indexWeb web i = xy +-- | Consider at each node not just the connections to already known neighbours, but +-- also the connections to /their/ neighbours. If these next-neighbours turn out +-- to be actually situated closer, link to them directly. +knitShortcuts :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => MetricChoice x -> PointsWeb x y -> PointsWeb x y +knitShortcuts metricf w₀ = tweakWebGeometry metricf closeObtuseAngles + $ pseudoFixMaximise (rateLinkings w₀) w₀ + where pseudoFixMaximise oldBadness oldSt + | newBadness < oldBadness = pseudoFixMaximise newBadness newSt + | otherwise = newSt + where newSt = tweakWebGeometry metricf pickNewNeighbours + $ bidirectionaliseWebLinks oldSt + newBadness = rateLinkings newSt + rateLinkings :: PointsWeb x y -> Double + rateLinkings = geometricMeanOf rateNode . webLocalInfo + rateNode :: WebLocally x y -> Double + rateNode info = geometricMeanOf + (\(_, (δx,_)) -> info^.nodeLocalScalarProduct|$|δx) + $ info^.nodeNeighbours + + pickNewNeighbours :: WebLocally x y -> [WebNodeId] + pickNewNeighbours me = fst <$> go Nothing [] candidates + where go Nothing prevs (cs:ccs) = case bestNeighbours' lm' cs of + (links, Nothing) -> links + (links, Just newWall) + | Just _ <- me^.webBoundingPlane -> links + | otherwise -> + links ++ go (Just newWall) ((snd<$>links) ++ prevs) ccs + go (Just wall) prevs (cs:ccs) = case gatherGoodNeighbours + lm' lm wall prevs [] cs of + (links, Nothing) -> links + (links, Just newWall) + | Nothing <- me^.webBoundingPlane + , (_:_) <-ccs -> + links ++ go (Just newWall) ((snd<$>links) ++ prevs) ccs + | otherwise -> links + go _ _ [] = [] + lm' = me^.nodeLocalScalarProduct :: Metric x + lm = dualNorm lm' + candidates :: [[(WebNodeId, Needle x)]] + candidates = preferred : other + where (preferred, other) = case localOnion me [] of + _l₀:l₁:l₂:ls -> ( first _thisNodeId . swap <$> (l₁++l₂) + , map (first _thisNodeId . swap) <$> ls ) + [_l₀,l₁] -> (first _thisNodeId . swap <$> l₁, []) + + closeObtuseAngles :: WebLocally x y -> [WebNodeId] + closeObtuseAngles me = go [ (dv,v) ^/ sqrt (dv<.>^v) + | (i,(v,_)) <- me^.nodeNeighbours + , let dv = metric<$|v ] + candidates + where go :: [OSNeedle x] -> [OSNode x y] -> [WebNodeId] + go existing fillSrc = case constructUninhabitedCone existing of + Nothing -> fst <$> me^.nodeNeighbours + Just cone -> case findInCone cone fillSrc of + Just ((fv,filler),fillSrc') + -> (filler^.thisNodeId) : go (fv:existing) fillSrc' + Nothing -> fst <$> me^.nodeNeighbours + constructUninhabitedCone :: [OSNeedle x] -> Maybe (CPCone x) + constructUninhabitedCone vs = find (not.(`any`vs).includes) + $ coneBetween <$> choices dimension vs + where coneBetween :: [(Needle' x, a)] -> (Needle' x, ℝ) + coneBetween dvs = (coneDir, (coMetric|$|coneDir)/sqrt 2) + where coneDir = sumV $ fst <$> dvs + findInCone :: CPCone x -> [OSNode x y] + -> Maybe (OSNode x y, [OSNode x y]) + findInCone cone ((po,pn):ps) | cone`includes`po = Just ((po,pn), ps) + findInCone (coneDir, _) ((po,pn):_) + | Just wall <- pn^.webBoundingPlane + , BoundarylessWitness <- boundarylessWitness :: BoundarylessWitness x + , DualSpaceWitness <- dualSpaceWitness :: DualSpaceWitness (Needle x) + , testp <- pn^.thisNodeCoord .+~^ (coMetric<$|wall) + , (metric |$| testp.-~!me^.thisNodeCoord) > (metric|$|snd po) + = Nothing + findInCone cone (p:ps) = second (p:) <$> findInCone cone ps + findInCone _ [] = Nothing + includes :: CPCone x -> OSNeedle x -> Bool + (coneDir, narrowing)`includes`(_, v) = coneDir<.>^v >= narrowing + candidates :: [OSNode x y] + candidates = case localOnion me [] of + _l₀:_l₁:ls -> concat [ snd <$> sortBy (comparing fst) + [ (distSq, ((dv,v) ^/ sqrt distSq, node)) + | (v, node) <- layer + , let dv = metric<$|v + distSq = dv<.>^v ] + | layer <- ls ] + _ -> [] + metric = me^.nodeLocalScalarProduct + coMetric = dualNorm metric + dimension = subbasisDimension (entireBasis :: SubBasis (Needle x)) -coerceWebDomain :: ∀ a b y . (Manifold a, Manifold b, LocallyCoercible a b) +choices :: Int -> [a] -> [[a]] +choices n l = go n l id [] + where go 0 _ f = (f[]:) + go _ [] _ = id + go n (x:xs) f = go n xs f . go (n-1) xs ((x:).f) + +meanOf :: (Hask.Foldable f, Fractional n) => (a -> n) -> f a -> n +meanOf f = renormalise . Hask.foldl' accs (0, 0::Int) + where renormalise (acc,n) = acc/fromIntegral n + accs (acc,n) x = (acc+f x, succ n) + +geometricMeanOf :: (Hask.Foldable f, Floating n) => (a -> n) -> f a -> n +geometricMeanOf f = exp . meanOf (log . f) + + + +webBoundary :: WithField ℝ Manifold x => PointsWeb x y -> [(Cutplane x, y)] +webBoundary = webLocalInfo >>> Hask.toList >>> Hask.concatMap`id` + \info -> [ (Cutplane (info^.thisNodeCoord) (Stiefel1 wall), info^.thisNodeData) + | Just wall <- [info^.webBoundingPlane] ] + + +coerceWebDomain :: ∀ a b y + . (Manifold a, Manifold b, LocallyCoercible a b, SimpleSpace (Needle b)) => PointsWeb a y -> PointsWeb b y -coerceWebDomain (PointsWeb rsc assoc) - = case oppositeLocalCoercion :: CanonicalDiffeomorphism b a of - CanonicalDiffeomorphism - -> PointsWeb ( coerceShadeTree rsc ) - ( fmap (second $ localScalarProduct - %~transformNorm (arr $ coerceNeedle ([]::[(b,a)]))) - assoc ) +coerceWebDomain (PointsWeb web) = PointsWeb + $ unsafeFmapTree ( fmap $ \(x, Neighbourhood y ngbs lscl bndry) + -> ( locallyTrivialDiffeomorphism x + , Neighbourhood y ngbs + (coerceNorm ([]::[(a,b)]) lscl) + (fmap crcNeedle' bndry) ) ) + crcNeedle' coerceShade web + where crcNeedle' = case ( dualSpaceWitness :: DualSpaceWitness (Needle a) + , dualSpaceWitness :: DualSpaceWitness (Needle b) ) of + (DualSpaceWitness, DualSpaceWitness) -> arr $ coerceNeedle' ([]::[(a,b)]) data InterpolationIv y = InterpolationIv { @@ -441,7 +404,13 @@ , Geodesic x, Geodesic y ) => PointsWeb x y -> Cutplane x -> [(x,y)] sliceWeb_lin web = sliceEdgs - where edgs = webEdges web + where edgs :: [((x,y),(x,y))] + edgs = [ (gnodes i₀, gnodes i₁) + | (i₀,i₁) <- fastNub [ (i₀,i₁) + | (il,ir) <- edges graph + , let [i₀,i₁] = sort [il,ir] ] + ] + (graph, gnodes) = toGraph web sliceEdgs cp = [ (xi d, yi d) -- Brute-force search through all edges | ((x₀,y₀), (x₁,y₁)) <- edgs , Just d <- [cutPosBetween cp (x₀,x₁)] @@ -518,151 +487,189 @@ y₁ = maximum (snd<$>pts) pts = fst . fst <$> toList (localFocusWeb web) -webLocalInfo :: ∀ x y . WithField ℝ Manifold x - => PointsWeb x y -> PointsWeb x (WebLocally x y) -webLocalInfo origWeb = result - where result = wli $ localFocusWeb origWeb - wli (PointsWeb rsc asd) = PointsWeb rsc asd' - where asd' = Arr.imap localInfo asd - localInfo i (((x,y), ngbCo), ngbH) - = ( LocalWebInfo { - _thisNodeCoord = x - , _thisNodeData = y - , _thisNodeId = i - , _nodeNeighbours = [ (iNgb, (δx, neighbour)) - | iNgb <- UArr.toList $ ngbH^.neighbours - , let neighbour = unsafeIndexWebData result iNgb - Just δx = _thisNodeCoord neighbour.-~.x - ] - , _nodeLocalScalarProduct = ngbH^.localScalarProduct - , _nodeIsOnBoundary = anyUnopposed (ngbH^.localScalarProduct) ngbCo - }, ngbH ) - anyUnopposed rieM ngbCo = (`any`ngbCo) $ \(v,_) - -> not $ (`any`ngbCo) $ \(v',_) - -> (rieM<$|v) <.>^ v' < 0 + +hardbakeChunk :: WebChunk x y -> PointsWeb x y +hardbakeChunk = _thisChunk + +entireWeb :: PointsWeb x y -> WebChunk x y +entireWeb web = WebChunk web [] + localFocusWeb :: WithField ℝ Manifold x => PointsWeb x y -> PointsWeb x ((x,y), [(Needle x, y)]) -localFocusWeb (PointsWeb rsc asd) = PointsWeb rsc asd'' - where asd' = Arr.imap (\i (y,n) -> case indexShadeTree rsc i of - Right (_,x) -> ((x,y),n) ) asd - asd''= Arr.map (\((x,y),n) -> - (((x,y), [ ( case x'.-~.x of - Just v -> v - , y') - | j<-UArr.toList (n^.neighbours) - , let ((x',y'),_) = asd' Arr.! j - ]), n) - ) asd' +localFocusWeb = webLocalInfo >>> fmap `id`\n + -> ( (n^.thisNodeCoord, n^.thisNodeData) + , [ (δx, ngb^.thisNodeData) + | (_, (δx, ngb)) <- n^.nodeNeighbours ] ) + + +treewiseTraverseLocalWeb :: ∀ f x y . (WithField ℝ Manifold x, Hask.Applicative f) + => (WebLocally x y -> f y) + -> (∀ t i w . (Hask.Traversable t, Ord i) => (w -> f w) -> t (i, w) -> f (t w) ) + -> PointsWeb x y -> f (PointsWeb x y) +treewiseTraverseLocalWeb fl ct = fmap hardbakeChunk . twt . entireWeb + where twt = treewiseTraverseLocalWeb' fl $ ct twt + +treewiseTraverseLocalWeb' :: ∀ f x y . (WithField ℝ Manifold x, Hask.Applicative f) + => (WebLocally x y -> f y) + -> (NonEmpty (Int, WebChunk x y) -> f (NonEmpty (WebChunk x y))) + -> WebChunk x y -> f (WebChunk x y) +treewiseTraverseLocalWeb' fl ct domain + = $notImplemented{- + where rezoomed (PlainLeaves _) _ = localTraverseWebChunk fl domain + rezoomed tree pos + | pos == i₀, nLeaves tree == lDomain + = fmap reassemble $ ct (NE.zipWith + (\jb (i₀b,t') + -> (jb, domain & overrideStart .~ i₀+i₀b + & overriddenData + .~ Arr.slice i₀b (nLeaves t') domainData )) + (0:|[1..]) branches) + | otherwise = go branches + where branches = trunkBranches tree + go (_:|((i₀nb,nb):brs)) + | pos+i₀nb <= i₀ = go $ (i₀nb,nb):|brs + go ((i₀b,t):|_) = rezoomed t (pos+i₀b) + reassemble :: NonEmpty (WebChunk x y) -> WebChunk x y + reassemble brs = domain & overriddenData + .~ Hask.foldMap _overriddenData brs + lDomain = Arr.length domainData + -} + + + localOnion :: ∀ x y . WithField ℝ Manifold x - => WebLocally x y -> [[WebLocally x y]] -localOnion origin = go Map.empty $ Map.singleton (origin^.thisNodeId) (1, origin) - where go previous next + => WebLocally x y -> [WebNodeId] -> [[(Needle x, WebLocally x y)]] +localOnion origin directCandidates = map sortBCDistance . go Map.empty . Map.fromList + $ (origin^.thisNodeId, (1, origin)) : seeds + where seeds :: [(WebNodeId, (Int, WebLocally x y))] + seeds = [ (nid, (1, ninfo)) + | nid <- directCandidates + , (_,(_,ninfo)) <- origin^.nodeNeighbours + , ninfo^.thisNodeId == nid ] + go previous next | Map.null next = [] - | otherwise = ( snd <$> sortBy (comparing $ negate . fst) + | otherwise = ( computeOffset . snd + <$> sortBy (comparing $ negate . fst) (Hask.toList next) ) : go (Map.union previous next) (Map.fromListWith (\(n,ninfo) (n',_) -> (n+n'::Int, ninfo)) - [ (nnid,(1,nneigh)) - | (nid,(_,ninfo))<-Map.toList next - , (nnid,(_,nneigh))<-ninfo^.nodeNeighbours - , Map.notMember nnid previous ]) + [ (nnid,(1,nneigh)) + | (nid,(_,ninfo))<-Map.toList next + , (nnid,(_,nneigh))<-ninfo^.nodeNeighbours + , Map.notMember nnid previous && Map.notMember nnid next ]) + computeOffset p = case p^.thisNodeCoord .-~. origin^.thisNodeCoord of + Just v -> (v,p) + sortBCDistance = map snd . sortBy (comparing fst) . map (bcDist&&&id) + where bcDist (v,_) + = normSq (origin^.nodeLocalScalarProduct) $ v^-^seedBarycenterOffs + seedBarycenterOffs = sumV ngbOffs ^/ fromIntegral (length directCandidates + 1) + where ngbOffs = [ v | (_, (_, n)) <- seeds + , let Just v = n^.thisNodeCoord .-~. origin^.thisNodeCoord ] webOnions :: ∀ x y . WithField ℝ Manifold x => PointsWeb x y -> PointsWeb x [[(x,y)]] -webOnions = localFmapWeb (map (map $ _thisNodeCoord&&&_thisNodeData) . localOnion) +webOnions = localFmapWeb (map (map $ _thisNodeCoord&&&_thisNodeData <<< snd) + . (`localOnion`[])) -nearestNeighbour :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) +nearestNeighbour :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => PointsWeb x y -> x -> Maybe (x,y) -nearestNeighbour (PointsWeb rsc asd) x = fmap lkBest $ positionIndex empty rsc x - where lkBest (iEst, (_, xEst)) = (xProx, yProx) - where (iProx, (xProx, _)) = minimumBy (comparing $ snd . snd) - $ (iEst, (xEst, normSq locMetr vEst)) - : neighbours - (yProx, _) = asd Arr.! iProx - (_, Neighbourhood neighbourIds locMetr) = asd Arr.! iEst - neighbours = [ (i, (xNgb, normSq locMetr v)) - | i <- UArr.toList neighbourIds - , let Right (_, xNgb) = indexShadeTree rsc i - Just v = xNgb.-~.x - ] - Just vEst = xEst.-~.x +nearestNeighbour = webLocalInfo >>> \(PointsWeb rsc) x + -> fmap (fine x) (positionIndex empty rsc x) + where fine :: x -> (Int, ( [Shaded x (Neighbourhood x (WebLocally x y))] + , (x, Neighbourhood x (WebLocally x y)) )) + -> (x,y) + fine x (_, (_, (xc, (Neighbourhood c _ locMetr _)))) + = snd . minimumBy (comparing fst) + . map (first $ (c^.nodeLocalScalarProduct|$|) + . (^-^vc)) + $ (zeroV, (xc, c^.thisNodeData)) + : [ (δx, (ngb^.thisNodeCoord, ngb^.thisNodeData)) + | (_, (δx, ngb)) <- c^.nodeNeighbours ] + where Just vc = x.-~.xc -instance Hask.Functor (WebLocally x) where - fmap f (LocalWebInfo co dt id ng sp bn) - = LocalWebInfo co (f dt) id (map (second . second $ fmap f) ng) sp bn -instance WithField ℝ Manifold x => Comonad (WebLocally x) where - extract = _thisNodeData - extend f this@(LocalWebInfo co _ id ng sp bn) - = LocalWebInfo co (f this) id (map (second . second $ extend f) ng) sp bn - duplicate this@(LocalWebInfo co _ id ng sp bn) - = LocalWebInfo co this id (map (second $ second duplicate) ng) sp bn - -- ^ 'fmap' from the co-Kleisli category of 'WebLocally'. -localFmapWeb :: WithField ℝ Manifold x - => (WebLocally x y -> z) -> PointsWeb x y -> PointsWeb x z -localFmapWeb f = webLocalInfo >>> fmap f +localTraverseWeb :: (WithField ℝ Manifold x, Hask.Applicative m) + => (WebLocally x y -> m z) -> PointsWeb x y -> m (PointsWeb x z) +localTraverseWeb f = webLocalInfo >>> Hask.traverse f -traverseWebWithStrategy :: ( WithField ℝ Manifold x, Hask.Applicative m ) - => InconsistencyStrategy m x y -> (WebLocally x y -> Maybe y) - -> PointsWeb x y -> m (PointsWeb x y) -traverseWebWithStrategy strat f = webLocalInfo - >>> traverse (\info -> handleInconsistency strat - (info^.thisNodeData) (f info)) +-- ^ 'fmap' from the co-Kleisli category of 'WebLocally', restricted to some +-- contiguous part of a web. +localTraverseWebChunk :: (WithField ℝ Manifold x, Hask.Applicative m) + => (WebLocally x y -> m y) -> WebChunk x y -> m (WebChunk x y) +localTraverseWebChunk f (WebChunk this outlayers) + = fmap (\c -> WebChunk c outlayers) $ localTraverseWeb f this differentiateUncertainWebLocally :: ∀ x y - . ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , WithField ℝ Refinable y, SimpleSpace (Needle y) ) + . ( ModellableRelation x y ) => WebLocally x (Shade' y) -> Shade' (LocalLinear x y) -differentiateUncertainWebLocally info - = case estimateLocalJacobian - (info^.nodeLocalScalarProduct) - [ ( Local δx :: Local x, ngb^.thisNodeData ) - | (δx,ngb) <- (zeroV, info) - : (snd<$>info^.nodeNeighbours) - ] of - Just j -> j - _ -> Shade' zeroV mempty +differentiateUncertainWebLocally = duwl + ( dualSpaceWitness :: DualSpaceWitness (Needle x) + , dualSpaceWitness :: DualSpaceWitness (Needle y) ) + where duwl (DualSpaceWitness, DualSpaceWitness) info + = case fitLocally $ + (\(δx,ngb) -> (δx, ngb^.thisNodeData) ) + <$> (zeroV,info) : envi + of + Just (AffineModel _ j :: AffineModel x y) -> dualShade j + where _:directEnvi:remoteEnvi = localOnion info [] + envi = directEnvi ++ concat remoteEnvi + differentiateUncertainWebFunction :: ∀ x y - . ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , WithField ℝ Manifold y, SimpleSpace (Needle y), Refinable y ) + . ( ModellableRelation x y ) => PointsWeb x (Shade' y) -> PointsWeb x (Shade' (LocalLinear x y)) differentiateUncertainWebFunction = localFmapWeb differentiateUncertainWebLocally differentiate²UncertainWebLocally :: ∀ x y - . ( WithField ℝ Manifold x, FlatSpace (Needle x) - , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y) ) + . ( ModellableRelation x y ) => WebLocally x (Shade' y) -> Shade' (Needle x ⊗〃+> Needle y) differentiate²UncertainWebLocally = d²uwl - ( pseudoAffineWitness :: PseudoAffineWitness x - , pseudoAffineWitness :: PseudoAffineWitness y - , dualSpaceWitness :: DualSpaceWitness (Needle x) + ( dualSpaceWitness :: DualSpaceWitness (Needle x) , dualSpaceWitness :: DualSpaceWitness (Needle y) ) - where d²uwl ( PseudoAffineWitness (SemimanifoldWitness _) - , PseudoAffineWitness (SemimanifoldWitness _) - , DualSpaceWitness, DualSpaceWitness ) info - = case estimateLocalHessian $ - (\ngb -> case (ngb^.thisNodeCoord .-~. info^.thisNodeCoord) of - Just δx -> (Local δx :: Local x, ngb^.thisNodeData) ) - <$> info :| envi + where d²uwl (DualSpaceWitness, DualSpaceWitness) info + = case fitLocally $ + (\(δx,ngb) -> (δx, ngb^.thisNodeData) ) + <$> (zeroV,info) : envi of - QuadraticModel _ h -> dualShade $ projectShade - (fromEmbedProject (acoSnd.acoSnd ^/ 2) - (snd.snd ^* 2) ) h - where xVol :: SymmetricTensor ℝ (Needle x) - xVol = squareVs $ fst.snd<$>info^.nodeNeighbours - _:directEnvi:remoteEnvi = localOnion info - envi = directEnvi ++ take (nMinData - length directEnvi) (concat remoteEnvi) - nMinData = 1 + regular_neighboursCount - (subbasisDimension (entireBasis :: SubBasis (Needle x))) + Just (QuadraticModel _ _ h :: QuadraticModel x y) + -> linIsoTransformShade (2*^id) $ dualShade h + where _:directEnvi:remoteEnvi = localOnion info [] + envi = directEnvi ++ concat remoteEnvi + +-- | Calculate a quadratic fit with uncertainty margin centered around the connection +-- between any two adjacent nodes. In case of a regular grid (which we by no means +-- require here!) this corresponds to the vector quantities of an Arakawa type C/D +-- grid (cf. A. Arakawa, V.R. Lamb (1977): +-- Computational design of the basic dynamical processes of the UCLA general circulation model) +localModels_CGrid :: ∀ x y ㄇ . ( ModellableRelation x y, LocalModel ㄇ ) + => PointsWeb x (Shade' y) -> [(x, ㄇ x y)] +localModels_CGrid = Hask.concatMap theCGrid . Hask.toList . webLocalInfo + where theCGrid :: WebLocally x (Shade' y) -> [(x, ㄇ x y)] + theCGrid node = [ ( pn .-~^ δx^/2 + , propagationCenteredModel + ( LocalDataPropPlan + pn + (negateV δx) + (ngbNode^.thisNodeData) + (node^.thisNodeData) + (fmap (second _thisNodeData) + . concat . tail + $ localOnion ngbNode [node^.thisNodeId] ) + ) ) + | (nid, (δx, ngbNode)) <- node^.nodeNeighbours + , nid > node^.thisNodeId + , Just pn <- [toInterior $ ngbNode^.thisNodeCoord] + ] + + acoSnd :: ∀ s v y . ( Object (Affine s) y, Object (Affine s) v , LinearSpace v, Scalar v ~ s ) => Affine s y (v,y) acoSnd = case ( linearManifoldWitness :: LinearManifoldWitness v @@ -671,54 +678,37 @@ (LinearManifoldWitness BoundarylessWitness, DualSpaceWitness, DualSpaceWitness) -> const zeroV &&& id --- | Heuristic formula, matches the number of neighbours each vertex has in a one- --- and two-dimensional count -regular_neighboursCount :: Int -> Int -regular_neighboursCount d - | d>0 = (regular_neighboursCount (d-1) + 1)*2 - | otherwise = 0 - differentiate²UncertainWebFunction :: ∀ x y - . ( WithField ℝ Manifold x, FlatSpace (Needle x) - , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y) ) + . ( ModellableRelation x y ) => PointsWeb x (Shade' y) -> PointsWeb x (Shade' (Needle x ⊗〃+> Needle y)) differentiate²UncertainWebFunction = localFmapWeb differentiate²UncertainWebLocally -rescanPDELocally :: ∀ x y ð . - ( WithField ℝ Manifold x, FlatSpace (Needle x) - , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y) ) - => DifferentialEqn x ð y -> WebLocally x (Shade' y) - -> (Maybe (Shade' y), Maybe (Shade' ð)) +rescanPDELocally :: ∀ x y ㄇ . + ( ModellableRelation x y, LocalModel ㄇ ) + => DifferentialEqn ㄇ x y -> WebLocally x (Shade' y) -> Maybe (Shade' y) rescanPDELocally = case ( dualSpaceWitness :: DualNeedleWitness x , dualSpaceWitness :: DualNeedleWitness y , boundarylessWitness :: BoundarylessWitness x , pseudoAffineWitness :: PseudoAffineWitness y ) of ( DualSpaceWitness,DualSpaceWitness,BoundarylessWitness , PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness) ) - -> \f info -> let xc = info^.thisNodeCoord + -> \f info + -> if isJust $ info^.webBoundingPlane + then return $ info^.thisNodeData + else let xc = info^.thisNodeCoord yc = info^.thisNodeData.shadeCtr in case f $ coverAllAround (xc, yc) [ (δx, (ngb^.thisNodeData.shadeCtr.-~!yc) ^+^ v) | (_,(δx,ngb))<-info^.nodeNeighbours , v <- normSpanningSystem' (ngb^.thisNodeData.shadeNarrowness)] of - LocalDifferentialEqn _ rescan - -> rescan (info^.thisNodeData) - (differentiateUncertainWebLocally info) - (differentiate²UncertainWebLocally info) - -rescanPDEOnWeb :: ( WithField ℝ Manifold x, FlatSpace (Needle x) - , WithField ℝ Refinable y, Geodesic y, FlatSpace (Needle y) - , Hask.Applicative m ) - => InconsistencyStrategy m x (Shade' y, Shade' ð) - -> DifferentialEqn x ð y -> PointsWeb x (Shade' y) - -> m (PointsWeb x (Shade' y, Shade' ð)) -rescanPDEOnWeb strat deq = traverseWebWithStrategy strat - (fzip . rescanPDELocally deq . fmap fst) - . fmap (\shy -> (shy, error - "No default value for inconsistent PDE-rescanning on web")) + LocalDifferentialEqn rescan -> fst + ( rescan $ case fitLocally $ map (id *** _thisNodeData) + =<< (localOnion info []) of + Just ㄇ -> ㄇ) + >>= intersectShade's . (:|[info^.thisNodeData]) toGraph :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => PointsWeb x y -> (Graph, Vertex -> (x, y)) @@ -726,8 +716,8 @@ (graphFromEdges' edgs) where edgs :: [(Int, Int, [Int])] edgs = Arr.toList - . Arr.imap (\i (_, Neighbourhood ngbs _) -> (i, i, UArr.toList ngbs)) - $ webNodeAssocData wb + . Arr.imap (\i (Neighbourhood _ ngbs _ _) -> (i, i, (i+) <$> UArr.toList ngbs)) + . Arr.fromList . Hask.toList $ webNodeRsc wb @@ -798,6 +788,14 @@ Just r -> pure r Nothing -> throwE $ PropagationInconsistency n o ) +postponeInconsistencies :: Hask.Monad m => (NonEmpty υ -> Maybe υ) + -> InformationMergeStrategy [] (WriterT [PropagationInconsistency x υ] m) + (x,υ) υ +postponeInconsistencies merge = InformationMergeStrategy + (\o n -> case merge $ o :| fmap snd n of + Just r -> pure r + Nothing -> writer (o,[PropagationInconsistency n o]) ) + maybeAlt :: Hask.Alternative f => Maybe a -> f a maybeAlt (Just x) = pure x maybeAlt Nothing = Hask.empty @@ -809,14 +807,10 @@ deriving instance Hask.Functor (InconsistencyStrategy m x) -iterateFilterDEqn_static :: ( WithField ℝ Manifold x, FlatSpace (Needle x) - , Refinable y, Geodesic y, FlatSpace (Needle y) - , WithField ℝ AffineManifold ð, Geodesic ð - , SimpleSpace (Needle ð) - , Hask.MonadPlus m ) +iterateFilterDEqn_static :: ( ModellableRelation x y, Hask.MonadPlus m, LocalModel ㄇ ) => InformationMergeStrategy [] m (x,Shade' y) iy -> Embedding (->) (Shade' y) iy - -> DifferentialEqn x ð y + -> DifferentialEqn ㄇ x y -> PointsWeb x (Shade' y) -> Cofree m (PointsWeb x (Shade' y)) iterateFilterDEqn_static strategy shading f = fmap (fmap (shading >-$)) @@ -824,42 +818,154 @@ . fmap (shading $->) -filterDEqnSolutions_static :: ∀ x y iy ð m . - ( WithField ℝ Manifold x, FlatSpace (Needle x) - , Refinable y, Geodesic y, FlatSpace (Needle y) - , WithField ℝ AffineManifold ð, Geodesic ð - , SimpleSpace (Needle ð) - , Hask.MonadPlus m ) +iterateFilterDEqn_static_selective :: ( ModellableRelation x y + , Hask.MonadPlus m, badness ~ ℝ + , LocalModel ㄇ ) + => InformationMergeStrategy [] m (x,Shade' y) iy + -> Embedding (->) (Shade' y) iy + -> (x -> iy -> badness) + -> DifferentialEqn ㄇ x y + -> PointsWeb x (Shade' y) -> Cofree m (PointsWeb x (Shade' y)) +iterateFilterDEqn_static_selective strategy shading badness f + = fmap (fmap (shading >-$)) + . coiter (filterDEqnSolutions_static_selective strategy shading badness f) + . fmap (shading $->) + + +filterDEqnSolutions_static :: ∀ x y ㄇ iy m . + ( ModellableRelation x y, Hask.MonadPlus m, LocalModel ㄇ ) => InformationMergeStrategy [] m (x,Shade' y) iy -> Embedding (->) (Shade' y) iy - -> DifferentialEqn x ð y -> PointsWeb x iy -> m (PointsWeb x iy) + -> DifferentialEqn ㄇ x y -> PointsWeb x iy -> m (PointsWeb x iy) filterDEqnSolutions_static = case geodesicWitness :: GeodesicWitness y of GeodesicWitness _ -> \strategy shading f -> webLocalInfo >>> fmap (id &&& rescanPDELocally f . fmap (shading>-$)) >>> localFocusWeb >>> Hask.traverse ( \((_,(me,updShy)), ngbs) -> let oldValue = me^.thisNodeData :: iy - in case updShy of - (Just shy, Just shð) -> case ngbs of + in if isJust $ me ^. webBoundingPlane + then return oldValue + else case updShy of + Just shy -> case ngbs of [] -> pure oldValue _:_ | BoundarylessWitness <- (boundarylessWitness::BoundarylessWitness x) - -> maybeAlt - ( sequenceA [ fzip sj - >>= \ngbShyð -> (ngbInfo^.thisNodeCoord,)<$> + -> sequenceA [ maybeAlt sj + >>= \ngbShyð -> fmap ((me^.thisNodeCoord .+~^ δx,) + . (shading>-$)) + . mergeInformation strategy oldValue . Hask.toList + $ (ngbInfo^.thisNodeCoord,)<$> propagateDEqnSolution_loc f (LocalDataPropPlan (ngbInfo^.thisNodeCoord) (negateV δx) ngbShyð - (shy, shð) - (fmap (second ((shading>-$) . _thisNodeData) - . snd) $ ngbInfo^.nodeNeighbours) + shy + (fmap (second ((shading>-$) . _thisNodeData)) + . concat . tail $ localOnion ngbInfo + [me^.thisNodeId]) ) | (δx, (ngbInfo,sj)) <- ngbs - ] ) + ] >>= mergeInformation strategy (shading$->shy) _ -> mergeInformation strategy oldValue empty ) + + +data Average a = Average { weight :: Int + , averageAcc :: a + } deriving (Hask.Functor) +instance Num a => Monoid (Average a) where + mempty = Average 0 0 + mappend (Average w₀ a₀) (Average w₁ a₁) = Average (w₀+w₁) (a₀+a₁) +instance Hask.Applicative Average where + pure = Average 1 + Average w₀ a₀ <*> Average w₁ a₁ = Average (w₀*w₁) (a₀ a₁) + +average :: Fractional a => Average a -> a +average (Average w a) = a / fromIntegral w + +averaging :: VectorSpace a => [a] -> Average a +averaging l = Average (length l) (sumV l) + +filterDEqnSolutions_static_selective :: ∀ x y ㄇ iy m badness . + ( ModellableRelation x y + , Hask.MonadPlus m, badness ~ ℝ + , LocalModel ㄇ ) + => InformationMergeStrategy [] m (x,Shade' y) iy -> Embedding (->) (Shade' y) iy + -> (x -> iy -> badness) + -> DifferentialEqn ㄇ x y + -> PointsWeb x iy -> m (PointsWeb x iy) +filterDEqnSolutions_static_selective = case geodesicWitness :: GeodesicWitness y of + GeodesicWitness _ -> \strategy shading badness f + -> -- Integration step: determine at each point from the function values + -- what the derivatives should be, and use them to propagate the solution + -- in all directions. We only spend a single computation step on regions + -- where nothing much changes (indicating the a-priori information is + -- too weak yet to make any predictions anyway), but multiple steps in + -- regions where good progress is noticed. + fmap fst . (runWriterT :: WriterT (Average badness) m (PointsWeb x iy) + -> m (PointsWeb x iy, Average badness)) + . treewiseTraverseLocalWeb ( \me + -> let oldValue = me^.thisNodeData :: iy + badHere = badness $ me^.thisNodeCoord + oldBadness = badHere oldValue + in if isJust $ me ^. webBoundingPlane + then return oldValue + else case me^.nodeNeighbours of + [] -> pure oldValue + _:_ | BoundarylessWitness <- (boundarylessWitness::BoundarylessWitness x) + -> WriterT . fmap (\updated + -> (updated, pure (oldBadness / badHere updated))) + $ sequenceA [ fmap ((me^.thisNodeCoord .+~^ δx,) + . (shading>-$)) + . mergeInformation strategy oldValue . Hask.toList + $ (ngbInfo^.thisNodeCoord,)<$> + propagateDEqnSolution_loc + f (LocalDataPropPlan + (ngbInfo^.thisNodeCoord) + (negateV δx) + (shading >-$ ngbInfo^.thisNodeData) + (shading >-$ oldValue) + (fmap (second ((shading>-$) . _thisNodeData)) + . concat . tail $ localOnion + ngbInfo [me^.thisNodeId] ) + ) + | (_, (δx, ngbInfo)) <- me^.nodeNeighbours + ] + >>= mergeInformation strategy oldValue ) + (\combiner branchData -> WriterT $ do + (branchResults,improvements) + <- runWriterT $ Hask.traverse + (\(i,branch) -> fmap (i,) + . censor (pure . (i,) . average) + $ combiner branch) + branchData + let (best, _) = maximumBy (comparing snd) improvements + (branchResults',improvements') + <- runWriterT $ Hask.traverse + (\(i,branch) -> if i==best + then censor (pure . (i,) . average) + $ combiner branch + else WriterT $ return (branch, pure (i,1)) ) + branchResults + return ( branchResults' + , liftA2 (*) (averaging $ snd<$>improvements) + (averaging $ snd<$>improvements') ) + ) + >=> -- Boundary-condition / differentiation step: update the local values + -- based on a-priori boundary conditions, possibly dependent on + -- numerical derivatives of the current solution estimate. + localTraverseWeb (\me -> fmap (shading$->) + . maybeAlt $ rescanPDELocally f me) + . fmap (shading>-$) + +-- | The +-- is insufficiently polymorphic, requiring @w ~ w'@. +censor :: Functor m (->) (->) => (w -> w') -> WriterT w m a -> WriterT w' m a +censor = mapWriterT . fmap . second + + + handleInconsistency :: InconsistencyStrategy m x a -> a -> Maybe a -> m a handleInconsistency AbortOnInconsistency _ i = i handleInconsistency IgnoreInconsistencies _ (Just v) = Identity v @@ -887,14 +993,13 @@ oldAndNew' (_, l) = (False,) <$> l -filterDEqnSolutions_adaptive :: ∀ x y ð badness m - . ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , WithField ℝ AffineManifold y, Refinable y, Geodesic y - , WithField ℝ AffineManifold ð, Geodesic ð, SimpleSpace (Needle ð) - , badness ~ ℝ, Hask.Monad m ) +filterDEqnSolutions_adaptive :: ∀ x y ㄇ ð badness m + . ( ModellableRelation x y, AffineManifold y + , badness ~ ℝ, Hask.Monad m + , LocalModel ㄇ ) => MetricChoice x -- ^ Scalar product on the domain, for regularising the web. -> InconsistencyStrategy m x (Shade' y) - -> DifferentialEqn x ð y + -> DifferentialEqn ㄇ x y -> (x -> Shade' y -> badness) -> PointsWeb x (SolverNodeState x y) -> m (PointsWeb x (SolverNodeState x y)) @@ -905,7 +1010,7 @@ -> m (PointsWeb x ( (WebLocally x (SolverNodeState x y) , [(Shade' y, badness)]) )) tryPreproc BoundarylessWitness (GeodesicWitness _) - = traverse addPropagation $ webLocalInfo oldState + = Hask.traverse addPropagation $ webLocalInfo oldState where addPropagation wl | null neighbourInfo = pure (wl, []) | otherwise = (wl,) . map (id&&&badness undefined) @@ -917,9 +1022,8 @@ (neigh^.thisNodeCoord) (negateV δx) (convexSetHull $ neigh^.thisNodeData - .solverNodeStatus - , undefined) - (thisShy, undefined) + .solverNodeStatus) + thisShy [ second (convexSetHull . _solverNodeStatus . _thisNodeData) nn | (_,nn)<-neigh^.nodeNeighbours ] ) @@ -1013,9 +1117,8 @@ (n^.thisNodeCoord) (stepV ^-^ δx) (convexSetHull $ - n^.thisNodeData.solverNodeStatus - , undefined) - (aprioriInterpolate, undefined) + n^.thisNodeData.solverNodeStatus) + aprioriInterpolate (second (convexSetHull ._solverNodeStatus ._thisNodeData) @@ -1067,8 +1170,7 @@ , let Just vOld = ngb^.thisNodeCoord .-~. xOld ] -recomputeJacobian :: ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , WithField ℝ Manifold y, SimpleSpace (Needle y), Refinable y ) +recomputeJacobian :: ( ModellableRelation x y ) => PointsWeb x (SolverNodeState x y) -> PointsWeb x (SolverNodeState x y) recomputeJacobian = webLocalInfo @@ -1079,12 +1181,11 @@ iterateFilterDEqn_adaptive - :: ( WithField ℝ Manifold x, SimpleSpace (Needle x) - , WithField ℝ AffineManifold y, Refinable y, Geodesic y, Hask.Monad m - , WithField ℝ AffineManifold ð, Geodesic ð, SimpleSpace (Needle ð) ) + :: ( ModellableRelation x y, AffineManifold y + , LocalModel ㄇ, Hask.Monad m ) => MetricChoice x -- ^ Scalar product on the domain, for regularising the web. -> InconsistencyStrategy m x (Shade' y) - -> DifferentialEqn x ð y + -> DifferentialEqn ㄇ x y -> (x -> Shade' y -> ℝ) -- ^ Badness function for local results. -> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)] iterateFilterDEqn_adaptive mf strategy f badness diff --git a/Data/Manifold/Web/Internal.hs b/Data/Manifold/Web/Internal.hs new file mode 100644 --- /dev/null +++ b/Data/Manifold/Web/Internal.hs @@ -0,0 +1,616 @@ +-- | +-- Module : Data.Manifold.Web.Internal +-- Copyright : (c) Justus Sagemüller 2017 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE DeriveFoldable #-} +{-# LANGUAGE DeriveTraversable #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE UnicodeSyntax #-} + + +module Data.Manifold.Web.Internal where + + +import Prelude hiding ((^)) + +import qualified Data.Vector.Unboxed as UArr + +import Data.Manifold.Types +import Data.Manifold.Types.Primitive +import Data.Manifold.PseudoAffine +import Data.Manifold.Shade +import Data.Manifold.TreeCover +import Data.Function.Affine +import Data.VectorSpace (Scalar, (^+^), (^/), (^*), sumV) +import Math.LinearMap.Category ( SimpleSpace, LSpace, DualVector, Norm, Variance + , (<.>^), dualNorm, (<$|), (|$|), normSq + , dualSpaceWitness, DualSpaceWitness(..) + , FiniteDimensional (..) ) + +import qualified Data.Foldable as Hask +import qualified Data.Traversable as Hask +import Data.List (sortBy) +import qualified Data.List.NonEmpty as NE +import Data.List.NonEmpty (NonEmpty((:|))) +import qualified Data.Map as Map +import qualified Data.Foldable.Constrained as CCt +import Data.Functor.Identity +import Data.Function ((&)) +import Data.Ord (comparing) +import Data.List.FastNub (fastNub) +import qualified Data.IntSet as ℤSet +import Data.IntSet (IntSet) +import Data.Maybe (isNothing) +import Control.Arrow +import Control.Monad (guard, forM_) +import Control.Comonad +import Control.Monad.Trans.State + +import Control.DeepSeq + +import GHC.Generics (Generic) + +import Control.Lens +import Control.Lens.TH + +import Data.CallStack (HasCallStack) + + +type WebNodeId = Int +type WebNodeIdOffset = Int + +data Neighbourhood x y = Neighbourhood { + _dataAtNode :: y + , _neighbours :: UArr.Vector WebNodeIdOffset + , _localScalarProduct :: Metric x + , _webBoundaryAtNode :: Maybe (Needle' x) + } + deriving (Generic, Functor, Foldable, Traversable) +makeLenses ''Neighbourhood + +deriving instance ( WithField ℝ PseudoAffine x + , SimpleSpace (Needle x), Show (Needle' x), Show y ) + => Show (Neighbourhood x y) + +data WebLocally x y = LocalWebInfo { + _thisNodeCoord :: x + , _thisNodeData :: y + , _thisNodeId :: WebNodeId + , _nodeNeighbours :: [(WebNodeId, (Needle x, WebLocally x y))] + , _nodeLocalScalarProduct :: Metric x + , _webBoundingPlane :: Maybe (Needle' x) + } deriving (Generic) +makeLenses ''WebLocally + +data NeighbourhoodVector x = NeighbourhoodVector + { _nvectId :: Int + , _theNVect :: Needle x + , _nvectNormal :: Needle' x + , _nvectLength :: Scalar (Needle x) + , _otherNeighboursOverlap :: Scalar (Needle x) + } +makeLenses ''NeighbourhoodVector + +data PropagationInconsistency x υ = PropagationInconsistency { + _inconsistentPropagatedData :: [(x,υ)] + , _inconsistentAPrioriData :: υ } + | PropagationInconsistencies [PropagationInconsistency x υ] + deriving (Show) +makeLenses ''PropagationInconsistency + +instance Monoid (PropagationInconsistency x υ) where + mempty = PropagationInconsistencies [] + mappend p q = mconcat [p,q] + mconcat = PropagationInconsistencies + +instance (NFData x, NFData (Metric x), NFData (Needle' x), NFData y) + => NFData (Neighbourhood x y) + +-- | A 'PointsWeb' is almost, but not quite a mesh. It is a stongly connected† +-- directed graph, backed by a tree for fast nearest-neighbour lookup of points. +-- +-- †In general, there can be disconnected components, but every connected +-- component is strongly connected. +newtype PointsWeb :: * -> * -> * where + PointsWeb :: { + webNodeRsc :: x`Shaded`Neighbourhood x y + } -> PointsWeb x y + deriving (Generic, Functor, Foldable, Traversable) + +instance (NFData x, NFData (Metric x), NFData (Needle' x), NFData y) => NFData (PointsWeb x y) + +instance CCt.Foldable (PointsWeb x) (->) (->) where + ffoldl = uncurry . Hask.foldl' . curry + foldMap = Hask.foldMap + + +data WebChunk x y = WebChunk { + _thisChunk :: PointsWeb x y + , _layersAroundChunk :: [(x`Shaded`Neighbourhood x y, WebNodeId)] + } + +makeLenses ''WebChunk + +data NodeInWeb x y = NodeInWeb { + _thisNodeOnly :: (x, Neighbourhood x y) + , _layersAroundNode :: [(x`Shaded`Neighbourhood x y, WebNodeId)] + } +makeLenses ''NodeInWeb + +data PathStep x y = PathStep { + _pathStepStart :: WebLocally x y + , _pathStepEnd :: WebLocally x y + } +makeLenses ''PathStep + + +type MetricChoice x = Shade x -> Metric x + + +traverseInnermostChunks :: ∀ f x y z . Applicative f + => (WebChunk x y -> f (PointsWeb x z)) -> PointsWeb x y -> f (PointsWeb x z) +traverseInnermostChunks f = go [] + where go :: [(x`Shaded`Neighbourhood x y, WebNodeId)] -> PointsWeb x y -> f (PointsWeb x z) + go outlayers (w@(PointsWeb (PlainLeaves _))) + = f (WebChunk w outlayers) + go outlayers (PointsWeb w) = PointsWeb <$> traverseTrunkBranchChoices travel w + where travel :: (Int, (Shaded x (Neighbourhood x y))) + -> Shaded x (Neighbourhood x y) + -> f (Shaded x (Neighbourhood x z)) + travel (i₀, br) obrs + = webNodeRsc <$> go ((obrs,i₀) : outlayers) (PointsWeb br) + +traverseNodesInEnvi :: ∀ f x y z . Applicative f + => (NodeInWeb x y -> f (Neighbourhood x z)) + -> PointsWeb x y -> f (PointsWeb x z) +traverseNodesInEnvi f = traverseInnermostChunks fc + where fc :: WebChunk x y -> f (PointsWeb x z) + fc (WebChunk (PointsWeb (PlainLeaves lvs)) outlayers) + = PointsWeb . PlainLeaves <$> Hask.traverse fn (ixedFoci lvs) + where fn ((i, (x, ngbh)), nearbyLeaves) + = (x,) <$> f (NodeInWeb (x,ngbh) + $ (PlainLeaves nearbyLeaves, i) : outlayers) + +fmapNodesInEnvi :: (NodeInWeb x y -> Neighbourhood x z) -> PointsWeb x y -> PointsWeb x z +fmapNodesInEnvi f = runIdentity . traverseNodesInEnvi (Identity . f) + + +ixedFoci :: [a] -> [((Int, a), [a])] +ixedFoci = go 0 + where go _ [] = [] + go i (x:xs) = ((i,x), xs) : map (second (x:)) (go (i+1) xs) + + +indexWeb :: PointsWeb x y -> WebNodeId -> Maybe (x,y) +indexWeb (PointsWeb rsc) i = case indexShadeTree rsc i of + Right (_, (x, Neighbourhood y _ _ _)) -> Just (x, y) + _ -> Nothing + +unsafeIndexWebData :: PointsWeb x y -> WebNodeId -> y +unsafeIndexWebData web i = case indexWeb web i of + Just (x,y) -> y + + +jumpNodeOffset :: WebNodeIdOffset -> NodeInWeb x y -> NodeInWeb x y +jumpNodeOffset 0 node = node +jumpNodeOffset δi (NodeInWeb x environment) + = case zoomoutWebChunk δie $ WebChunk (PointsWeb $ PlainLeaves [x]) environment of + (WebChunk bigChunk envi', δi') + -> case pickNodeInWeb bigChunk δi' of + NodeInWeb x' envi'' -> NodeInWeb x' $ envi'' ++ envi' + where δie | δi < 0 = δi + | otherwise = δi - 1 + +webAroundChunk :: WebChunk x y -> PointsWeb x y +webAroundChunk (WebChunk chunk []) = chunk +webAroundChunk (WebChunk (PointsWeb (PlainLeaves lvs)) + ((PlainLeaves lvsAround, i) : envi)) + = webAroundChunk $ WebChunk (PointsWeb . PlainLeaves $ lvsBefore++lvs++lvsAfter) envi + where (lvsBefore, lvsAfter) = splitAt i lvsAround +webAroundChunk (WebChunk (PointsWeb chunk) + ((OverlappingBranches nw ew (DBranch dir + (Hourglass (PlainLeaves[]) d) :| brs), 0) : envi)) + = webAroundChunk $ WebChunk (PointsWeb $ OverlappingBranches (nw+nLeaves chunk) ew + (DBranch dir (Hourglass chunk d) :| brs)) + envi +webAroundChunk (WebChunk (PointsWeb chunk) + ((OverlappingBranches nw ew (DBranch dir + (Hourglass u (PlainLeaves[])) :| brs), i) : envi)) + | i==nLeaves u + = webAroundChunk $ WebChunk (PointsWeb $ OverlappingBranches (nw+nLeaves chunk) ew + (DBranch dir (Hourglass u chunk) :| brs)) + envi +webAroundChunk (WebChunk chunk + (( OverlappingBranches nw ew (br₀@(DBranch _ (Hourglass u d)) + :|br₁:brs) + , i) : envi)) + = case webAroundChunk (WebChunk chunk [(OverlappingBranches nw ew (br₁:|brs), i')]) + of PointsWeb (OverlappingBranches nw' ew' (br₁':|brs')) + -> webAroundChunk $ WebChunk + (PointsWeb $ OverlappingBranches nw' ew' (br₀:|br₁':brs')) + envi + where i' = i - nLeaves u - nLeaves d +webAroundChunk (WebChunk _ ((OverlappingBranches nw ew branches, i):_)) + = error $ "Environment with branch sizes "++show (fmap nLeaves . Hask.toList<$>(Hask.toList branches)) + ++" does not have a gap at #"++show i +webAroundChunk (WebChunk _ ((PlainLeaves _, _):_)) + = error "Encountered non-PlainLeaves chunk in a PlainLeaves environment." + + +zoomoutWebChunk :: WebNodeIdOffset -> WebChunk x y -> (WebChunk x y, WebNodeId) +zoomoutWebChunk δi (WebChunk chunk ((outlayer, olp) : outlayers)) + | δi < -olp || δi >= nLeaves outlayer - olp + = zoomoutWebChunk δiOut $ WebChunk widerChunk outlayers + | otherwise = (WebChunk widerChunk outlayers, δiIn) + where δiOut | δi < 0 = δi + olp + | otherwise = δi + olp - nLeaves outlayer + δiIn | δi < 0 = δi + olp + | otherwise = δi + olp + nLeaves (webNodeRsc chunk) + widerChunk = webAroundChunk $ WebChunk chunk [(outlayer,olp)] +zoomoutWebChunk δi (WebChunk _ e) + = error $ "Can't zoom out δ"++show δi + ++" from a chunk with "++show (length e)++" environment layers." + +pickNodeInWeb :: PointsWeb x y -> WebNodeId -> NodeInWeb x y +pickNodeInWeb = go [] id + where go _ _ (PointsWeb w) i + | i<0 || i>=n = error + $ "Trying to pick node #"++show i++" in web with "++show n++" nodes." + where n = nLeaves w + go lyrsAcc lMod (PointsWeb (PlainLeaves lvs)) i + | (preds, node:succs)<-splitAt i lvs + = NodeInWeb node $ lMod (PlainLeaves $ preds++succs, i) : lyrsAcc + go lyrsAcc lMod + (PointsWeb (OverlappingBranches nw ew (DBranch dir (Hourglass u d):|brs))) i + | i < nu = go (lMod (OverlappingBranches (nw-nu) ew + (DBranch dir (Hourglass gap d):|brs), 0) : lyrsAcc) + id (PointsWeb u) i + | i < nu+nd = go (lMod (OverlappingBranches (nw-nd) ew + (DBranch dir (Hourglass u gap):|brs), nu) : lyrsAcc) + id (PointsWeb d) (i-nu) + | (b:rs)<-brs = go + lyrsAcc + (lMod . \(OverlappingBranches nwe ewe brse, ne) + -> ( OverlappingBranches (nwe+nu+nd) ewe + $ NE.cons (DBranch dir (Hourglass u d)) brse + , ne+nu+nd ) ) + (PointsWeb $ OverlappingBranches (nw-nu-nd) ew (b:|rs)) + (i-nu-nd) + where gap = PlainLeaves [] + [nu,nd] = nLeaves<$>[u,d] + + +webLocalInfo :: ∀ x y . WithField ℝ Manifold x + => PointsWeb x y -> PointsWeb x (WebLocally x y) +webLocalInfo = fmapNodesInEnvi linkln + where linkln :: NodeInWeb x y -> Neighbourhood x (WebLocally x y) + linkln node@(NodeInWeb (x, locloc@(Neighbourhood y ngbs metric nBoundary)) envis) + = locloc & dataAtNode .~ LocalWebInfo { + _thisNodeCoord = x + , _thisNodeData = y + , _thisNodeId = i + , _nodeNeighbours = [ (i + δi, (δx, ngb)) + | δi <- UArr.toList ngbs + , let ngbNode@(NodeInWeb (xn, _) _) + = jumpNodeOffset δi node + Just δx = xn .-~. x + Neighbourhood ngb _ _ _ = linkln ngbNode ] + , _nodeLocalScalarProduct = metric + , _webBoundingPlane = nBoundary + } + where i = foldr ((+) . snd) 0 envis + + +instance Functor (WebLocally x) where + fmap f (LocalWebInfo co dt id ng sp bn) + = LocalWebInfo co (f dt) id (map (second . second $ fmap f) ng) sp bn +instance WithField ℝ Manifold x => Comonad (WebLocally x) where + extract = _thisNodeData + extend f this@(LocalWebInfo co _ id ng sp bn) + = LocalWebInfo co (f this) id (map (second . second $ extend f) ng) sp bn + duplicate this@(LocalWebInfo co _ id ng sp bn) + = LocalWebInfo co this id (map (second $ second duplicate) ng) sp bn + +-- ^ 'fmap' from the co-Kleisli category of 'WebLocally'. +localFmapWeb :: WithField ℝ Manifold x + => (WebLocally x y -> z) -> PointsWeb x y -> PointsWeb x z +localFmapWeb f = webLocalInfo >>> fmap f + + +tweakWebGeometry :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) + => MetricChoice x -> (WebLocally x y -> [WebNodeId]) + -> PointsWeb x y -> PointsWeb x y +tweakWebGeometry metricf reknit = webLocalInfo >>> fmapNodesInEnvi`id` + \(NodeInWeb (x₀, (Neighbourhood info _ lm bound)) _) + -> let lm' = metricf . Shade (inInterior x₀) $ dualNorm lm + in Neighbourhood (info^.thisNodeData) + (UArr.fromList . map (subtract $ info^.thisNodeId) + $ reknit info) + lm' bound + + +bidirectionaliseWebLinks :: ∀ x y . PointsWeb x y -> PointsWeb x y +bidirectionaliseWebLinks web@(PointsWeb wnrsrc) = fmapNodesInEnvi bdse web + where bdse :: NodeInWeb x y -> Neighbourhood x y + bdse (NodeInWeb (x, Neighbourhood y outgn lm bound) envis) + = Neighbourhood y (UArr.fromList . fastNub $ incmn ++ UArr.toList outgn) + lm bound + where i = foldr ((+) . snd) 0 envis + incmn = case i `Map.lookup` incoming of + Just o -> subtract i<$>o + Nothing -> [] + incoming = Map.fromListWith (++) $ Hask.foldl' + (\(i,acc) (Neighbourhood _ outgn _ _) + -> (i+1, acc . (((,[i]).(+i)<$>UArr.toList outgn)++)) ) + (0,id) wnrsrc `snd` [] + + + +pumpHalfspace :: ∀ v . (SimpleSpace v, Scalar v ~ ℝ) + => Norm v + -> v -- ^ A vector @v@ for which we want @dv<.>^v ≥ 0@. + -> (DualVector v, [v]) -- ^ A plane @dv₀@ and some vectors @ws@ with @dv₀<.>^w ≥ 0@, + -- which should also fulfill @dv<.>^w ≥ 0@. + -> Maybe (DualVector v) -- ^ The plane @dv@ fulfilling these properties, if possible. +pumpHalfspace rieM v (prevPlane, ws) = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> + let ϑs = fmap (\u -> let x = prevPlane<.>^u + y = thisPlane<.>^u + in atan2 (x-y) (x+y)) $ v:ws + -- ϑ = 0 means we are mid-between the planes, ϑ > π/2 means we are past + -- `thisPlane`, ϑ < -π/2 we are past `prevPlane`. In other words, positive ϑ + -- mean we should mix in more of `prevPlane`, negative more of `thisPlane`. + [ϑmin, ϑmax] = [minimum, maximum] <*> [ϑs] + δϑ = ϑmax - ϑmin + vNudged = v ^+^ sumV (zipWith (^*) ws smallPseudorandSeq) + -- Introduce a tiny contribution from the other vectors to avoid + -- a degenerate 1D-situation in which @thisPlane ∝ prevPlane@. + dv = rieM<$|vNudged + thisPlane = dv ^/ (dv<.>^vNudged) + cas ϑ = cos $ ϑ - pi/4 + in if δϑ <= pi then Just $ let ϑbest = ϑmin + δϑ/2 + in prevPlane^*cas ϑbest ^+^ thisPlane^*cas (-ϑbest) + else Nothing + +smallPseudorandSeq :: [ℝ] +smallPseudorandSeq = (*2^^(-45)) . fromIntegral <$> lcg 293633 + where lcg x = x : lcg ((a*x)`mod`m) + m = 2^31 - 1 + a = 963345 :: Int -- revised Park-Miller + +data LinkingBadness r = LinkingBadness + { gatherDirectionsBadness :: !r -- ^ Prefer picking neighbours at right angles + -- to the currently-explored-boundary. This + -- is needed while we still have to link to + -- points in different spatial directions. + , closeSystemBadness :: !r -- ^ Prefer points directly opposed to the + -- current boundary. This is useful when the + -- system of directions is already complete + -- and we want a nicely symmetric “ball” of + -- neighbours around each point. + } deriving (Functor) + +linkingUndesirability :: ℝ -- ^ Absolute-square distance (euclidean norm squared) + -> ℝ -- ^ Directional distance (distance from wall containing + -- all already known neighbours) + -> LinkingBadness ℝ + -- ^ “Badness” of this point as the next neighbour to link to. + -- In gatherDirections mode this is large if + -- the point is far away, but also if it is + -- right normal to the wall. The reason we punish this is that + -- adding two points directly opposed to each other would lead + -- to an ill-defined wall orientation, i.e. wrong normals + -- on the web boundary. +linkingUndesirability distSq wallDist + | wallDist >= 0 = LinkingBadness + { gatherDirectionsBadness = distSq^2 / max 0 (distSq-wallDist^2) + , closeSystemBadness = distSq - wallDist^2/2 + } + | otherwise = LinkingBadness (1/0) (1/0) + + +bestNeighbours :: ∀ i v . (SimpleSpace v, Scalar v ~ ℝ) + => Norm v -> [(i,v)] -> ([i], Maybe (DualVector v)) +bestNeighbours lm' = first (map fst) . bestNeighbours' lm' + +bestNeighbours' :: ∀ i v . (SimpleSpace v, Scalar v ~ ℝ) + => Norm v -> [(i,v)] -> ([(i,v)], Maybe (DualVector v)) +bestNeighbours' lm' = extractSmallestOn (\(_,v) -> Just $ lm'|$|v) >>> + \(Just ((c₀i,c₀δx), candidates)) -> case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> + let wall₀ = w₀ ^/ (lm|$|w₀) -- sqrt (w₀<.>^c₀δx) + where w₀ = lm'<$|c₀δx + in first ((c₀i,c₀δx):) + $ gatherGoodNeighbours lm' lm wall₀ [c₀δx] [] candidates + where lm = dualNorm lm' :: Variance v + +gatherGoodNeighbours :: ∀ i v . (SimpleSpace v, Scalar v ~ ℝ) + => Norm v -> Variance v + -> DualVector v -> [v] -> [(i,v)] + -> [(i, v)] -> ([(i,v)], Maybe (DualVector v)) +gatherGoodNeighbours lm' lm wall prev preserved cs + | dimension == 1 = case extractSmallestOn + (\(_,δx) -> do + let wallDist = - wall<.>^δx + guard (wallDist > 0) + return wallDist + ) cs of + Just (r, _) -> ([r], Nothing) + Nothing -> ([], Just wall) + where dimension = subbasisDimension (entireBasis :: SubBasis v) +gatherGoodNeighbours lm' lm wall prev preserved cs + = case dualSpaceWitness :: DualSpaceWitness v of + DualSpaceWitness -> + case extractSmallestOn + (\(_,δx) -> do + let wallDist = - wall<.>^δx + dx = lm' <$| δx + distSq = dx<.>^δx + βmin = minimum [ 1 - (dx<.>^δxo) / sqrt (distSq*distSqo) + -- β behaves basically like ϑ², where ϑ is + -- the angle between two neighbour candidates. + | (δxo, distSqo) <- prevWMag ] + guard (wallDist >= 0 && βmin > 1e-3) + return $ gatherDirectionsBadness + (linkingUndesirability distSq wallDist) / βmin ) + cs of + Just ((i,δx), cs') + | Just wall' <- pumpHalfspace lm' δx (wall,prev) + -> first ((i,δx):) + $ gatherGoodNeighbours lm' lm (wall'^/(lm|$|wall')) + (δx:prev) [] (preserved++cs') + | (_:_)<-cs' -> gatherGoodNeighbours lm' lm wall + prev ((i,δx):preserved) cs' + _ -> let closeSys ((i,δx):_) + | Nothing <- pumpHalfspace lm' δx (wall,prev) + = ([(i,δx)], Nothing) + closeSys (_:cs'') = closeSys cs'' + closeSys [] + | null closureCandidates = ([], Just wall) + | otherwise = ([], Nothing) + closureCandidates = + [ ((i,δx), badness) + | (i,δx) <- preserved++cs + , let wallDist = - wall<.>^δx + distSq = normSq lm' δx + , wallDist > 0 + , wallDist^2 > 1e-3 * distSq + , let badness = linkingUndesirability distSq wallDist ] + in closeSys . map fst $ + sortBy (comparing $ closeSystemBadness . snd) closureCandidates + where prevWMag = map (id &&& normSq lm') prev + + +extractSmallestOn :: Ord b => (a -> Maybe b) -> [a] -> Maybe (a, [a]) +extractSmallestOn f = extract . map (id &&& f) + where extract [] = Nothing + extract ((x, Just o):cs) = Just $ go (o,x) cs + extract ((x, Nothing):cs) = second (x:) <$> extract cs + go (_,refx) [] = (refx, []) + go (ref,refx) ((x, Just o):cs) + | o < ref = second (refx:) $ go (o,x) cs + go ref ((x, _):cs) = second (x:) $ go ref cs + +type WNIPath = [WebNodeId] +type NodeSet = ℤSet.IntSet + +traversePathInIWeb :: ∀ φ x y . (WithField ℝ Manifold x, Monad φ, HasCallStack) + => [WebNodeId] -> (PathStep x y -> φ y) + -> PointsWeb x (WebLocally x y) -> φ (PointsWeb x (WebLocally x y)) +traversePathInIWeb path f = go path + where go [] web = pure web + go [_] web = pure web + go (i₀:i₁:is) web = do + y' <- f $ PathStep p₀ p₁ + let Right (Identity web') + = treeLeaf i₁ (dataAtNode . thisNodeData . const $ pure y') + $ webNodeRsc web + go (i₁:is) $ PointsWeb web' + where Just (_, p₀) = indexWeb web i₀ + Just (_, p₁) = indexWeb web i₁ + +traversePathsTowards :: ∀ f φ x y + . (WithField ℝ Manifold x, Monad φ, Monad f, HasCallStack) + => WebNodeId -- ^ The node towards which the paths should converge. + -> (PathStep x y -> φ y) + -- ^ The action which to traverse along each path. + -> (∀ υ . WebLocally x y -> φ υ -> f υ) + -- ^ Initialisation/evaluation for each path-traversal. + -> PointsWeb x y -> f (PointsWeb x y) +traversePathsTowards target pathTravF routeInitF web + | Nothing <- sn = error $ "Node "++show target++" not in web." + | otherwise = fmap (fmap _thisNodeData) . (`execStateT`envied) . forM_ paths + $ \path@(p₀:_) -> StateT $ \webState -> do + let Just (_, node₀) = indexWeb webState p₀ + ((),) <$> routeInitF node₀ (traversePathInIWeb path pathTravF webState) + where envied = webLocalInfo $ bidirectionaliseWebLinks web + sn@(Just (targetPos,targetNode)) = indexWeb envied target + paths = go 0 ℤSet.empty False [[target]] [] + where go :: Int -> NodeSet -> Bool -> [WNIPath] -> [WNIPath] -> [WNIPath] + go targetDist visitedNodes boundaryCreepingInhibitor workers finishedThreads + = case continue (round (sqrt $ fromIntegral targetDist :: Double)) + visitedNodes boundaryCreepingInhibitor workers of + (_, [], _, newFinished) -> newFinished ++ finishedThreads + (visited', continuation, alternatives, newFinished) + -> let newThreads = filter (`ℤSet.notMember`visited') + (ℤSet.toList alternatives) + in go (targetDist+1) + (ℤSet.union visited' alternatives) + True + (continuation ++ fmap pure newThreads) + (newFinished ++ finishedThreads) + continue :: Int -> NodeSet -> Bool -> [WNIPath] + -> (NodeSet, [WNIPath], NodeSet, [WNIPath]) + continue _ visitedNodes _ [] = (visitedNodes, [], ℤSet.empty, []) + continue dfsDepth visitedNodes boundaryCreepingInhibitor ((cursor:nds):paths) + = case fst <$> sortBy (comparing snd) + [ goDfs dfsDepth (ℤSet.insert i visitedNodes) [] [i] + | i <- candidates ] of + ((preferred, (visited', pAlts)):alts) + | Nothing <- guard boundaryCreepingInhibitor + >> cursorNode ^. webBoundingPlane + -> case continue dfsDepth visited' + boundaryCreepingInhibitor paths of + (visited'', contin'', alts', newFin) + -> ( visited'' + , (preferred++cursor:nds):contin'' + , ℤSet.union (ℤSet.fromList + $ pAlts ++ (last . fst <$> alts)) + alts' + , newFin ) + alts -> case continue dfsDepth visitedNodes + boundaryCreepingInhibitor paths of + (visited'', contin'', alts', newFin) + -> ( visited'' + , contin'' + , ℤSet.union (ℤSet.fromList + $ last . fst <$> alts) + alts' + , if null nds then newFin + else (cursor:nds):newFin ) + where Just (cursorPos,cursorNode) = indexWeb envied cursor + tgtOpp = cursorNode^.nodeLocalScalarProduct + <$| targetPos .-~! cursorPos + candidates = [ ngb + | (ngb, (_, _)) <- cursorNode^.nodeNeighbours + , ngb`ℤSet.notMember`visitedNodes ] + goDfs :: Int -> NodeSet -> [WebNodeId] -> WNIPath + -> ((WNIPath, (NodeSet, [WebNodeId])), ℝ) + goDfs d2go visited' oldAlts (p:old) + = case sortBy (comparing snd) candidates of + ((preferred,oppositionQ):alts) + -> let visited'' = ℤSet.insert preferred visited' + alts' = filter (/=preferred) oldAlts + ++ (fst<$>alts) + in if d2go>1 + then goDfs (d2go-1) visited'' alts' + (preferred:p:old) + else ( (preferred:p:old, (visited'', alts')) + , oppositionQ ) + [] -> let δn = explorePos .-~! cursorPos + oppositionQ = tgtOpp<.>^δn + in ((p:old, (visited', oldAlts)), oppositionQ) + where Just (explorePos,exploreNode) = indexWeb envied p + candidates = [ (ngb, tgtOpp<.>^δn) + | (ngb, (_, ngbN)) <- exploreNode^.nodeNeighbours + , ngb`ℤSet.notMember`visited' + , isNothing (ngbN^.webBoundingPlane) + , let δn = ngbN^.thisNodeCoord .-~! cursorPos ] diff --git a/images/examples/TreesAndWebs/2D-cartesian-strangeaspect.png b/images/examples/TreesAndWebs/2D-cartesian-strangeaspect.png Binary files a/images/examples/TreesAndWebs/2D-cartesian-strangeaspect.png and b/images/examples/TreesAndWebs/2D-cartesian-strangeaspect.png differ diff --git a/images/examples/TreesAndWebs/2D-cartesiandisk.png b/images/examples/TreesAndWebs/2D-cartesiandisk.png Binary files a/images/examples/TreesAndWebs/2D-cartesiandisk.png and b/images/examples/TreesAndWebs/2D-cartesiandisk.png differ diff --git a/images/examples/TreesAndWebs/2D-hexa-honeycomb.png b/images/examples/TreesAndWebs/2D-hexa-honeycomb.png new file mode 100644 Binary files /dev/null and b/images/examples/TreesAndWebs/2D-hexa-honeycomb.png differ diff --git a/images/examples/TreesAndWebs/2D-hexagonal.png b/images/examples/TreesAndWebs/2D-hexagonal.png new file mode 100644 Binary files /dev/null and b/images/examples/TreesAndWebs/2D-hexagonal.png differ diff --git a/images/examples/TreesAndWebs/2D-normaldistrib.png b/images/examples/TreesAndWebs/2D-normaldistrib.png Binary files a/images/examples/TreesAndWebs/2D-normaldistrib.png and b/images/examples/TreesAndWebs/2D-normaldistrib.png differ diff --git a/images/examples/TreesAndWebs/2D-scatter.png b/images/examples/TreesAndWebs/2D-scatter.png Binary files a/images/examples/TreesAndWebs/2D-scatter.png and b/images/examples/TreesAndWebs/2D-scatter.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.1.0 +Version: 0.4.4.0 Category: Math Synopsis: Coordinate-free hypersurfaces Description: Manifolds, a generalisation of the notion of “smooth curves” or surfaces, @@ -28,7 +28,7 @@ Homepage: https://github.com/leftaroundabout/manifolds Maintainer: (@) sagemueller $ geo.uni-koeln.de Build-Type: Simple -Cabal-Version: >=1.10 +Cabal-Version: >=1.18 Extra-Doc-Files: images/examples/*.png, images/examples/ShadeCombinations/2Dconvolution-skewed.png images/examples/TreesAndWebs/*.png @@ -40,14 +40,14 @@ Library Build-Depends: base>=4.5 && < 6 - , manifolds-core == 0.4.1.0 + , manifolds-core == 0.4.4.0 , transformers , vector-space>=0.8 , free-vector-spaces>=0.1.1 , linear , MemoTrie , vector - , linearmap-category >= 0.3.2 && < 0.4 + , linearmap-category >= 0.3.4 && < 0.4 , containers , comonad , free @@ -56,8 +56,11 @@ , number-show >= 0.1 && < 0.2 , tagged , deepseq + , placeholders , lens + , call-stack , constrained-categories >= 0.2.3 && < 0.3.1 + , pragmatic-show other-extensions: FlexibleInstances , TypeFamilies , FlexibleContexts @@ -74,7 +77,9 @@ Data.Manifold.TreeCover Data.Manifold.Shade Data.Manifold.Web + Data.Manifold.Web.Internal Data.Manifold.DifferentialEquation + Data.Manifold.Function.LocalModel Data.SimplicialComplex Data.Function.Differentiable Data.Function.Affine @@ -95,3 +100,26 @@ Util.Associate Util.LtdShow default-language: Haskell2010 + +test-suite test + default-language: + Haskell2010 + type: + exitcode-stdio-1.0 + hs-source-dirs: + test/tasty + main-is: + test.hs + build-depends: + base >= 4 && < 5 + , tasty >= 0.7 + , tasty-hunit + , tasty-quickcheck + , manifolds + , pragmatic-show + , containers + , vector-space + , constrained-categories + , linearmap-category + , lens + diff --git a/test/tasty/test.hs b/test/tasty/test.hs new file mode 100644 --- /dev/null +++ b/test/tasty/test.hs @@ -0,0 +1,476 @@ +-- | +-- Module : Main +-- Copyright : (c) Justus Sagemüller 2017 +-- License : GPL v3 +-- +-- Maintainer : (@) sagemueller $ geo.uni-koeln.de +-- Stability : experimental +-- Portability : portable +-- + +{-# LANGUAGE OverloadedLists, TypeFamilies, FlexibleContexts, UndecidableInstances #-} +{-# LANGUAGE TypeOperators #-} + +module Main where + +import Data.Manifold.Types +import Data.Manifold.PseudoAffine +import Data.Manifold.TreeCover +import Data.Manifold.Web +import Data.Manifold.Web.Internal +import Data.Manifold.Function.LocalModel +import Data.VectorSpace +import Math.LinearMap.Category +import Prelude hiding (id, fst, snd) +import Control.Category.Constrained (id) +import Control.Arrow.Constrained (fst,snd) + +import Test.Tasty +import Test.Tasty.HUnit +import qualified Test.Tasty.QuickCheck as QC +import Test.Tasty.QuickCheck ((==>)) + +import Data.Foldable (toList) +import Data.List (nub) +import qualified Data.Graph as Graph +import qualified Data.Set as Set +import Control.Arrow +import Control.Lens + +import qualified Text.Show.Pragmatic as SP + + +main = defaultMain tests + +tests :: TestTree +tests = testGroup "Tests" + [ testGroup "Graph structure of webs" + [ testCase "Manually-defined empty web." + $ toList (fst $ toGraph emptyWeb) @?= [] + , testCase "Manually-defined single-point web." + $ toList (fst $ toGraph singletonWeb) @?= [[]] + , testCase "Manually-defined simple triangular web." + $ toList (fst $ toGraph triangularWeb) @?= [[1,2],[0,2],[0,1]] + , testCase "Manually-defined simple quadratic web." + $ toList (fst $ toGraph quadraticWeb) @?= [[1,2],[0,3],[0,3],[1,2]] + , testCase "Envi-aware traversal over simple quadratic web." + $ toList (fst . toGraph $ dummyWebFmap quadraticWeb) @?= [[1,2],[0,3],[0,3],[1,2]] + , testCase "Direct neighbours in empty web." + $ toList (directNeighbours emptyWeb) @?= [] + , testCase "Direct neighbours in single-point web." + $ toList (directNeighbours singletonWeb) @?= [[]] + , testCase "Direct neighbours in simple triangular web." + $ toList (directNeighbours triangularWeb) @?= [[1,2],[0,2],[0,1]] + , testCase "Direct neighbours in simple quadratic web." + $ toList (directNeighbours quadraticWeb) @?= [[1,2],[0,3],[0,3],[1,2]] + , testCase "Direct neighbours in quadratic web with one-direction diagonals." + $ toList (directNeighbours unidirDiagonalLinkedWeb) @?= [[1,2,3],[0,3],[0,1,3],[1,2]] + , testCase "Direct neighbours in 1-dir diag quadratic web after bidirectionalisation." + $ toList (directNeighbours $ bidirectionaliseWebLinks unidirDiagonalLinkedWeb) + @?= [[1,2,3],[0,2,3],[0,1,3],[0,1,2]] + , testCase "Direct neighbours in unsymmetric web." + $ toList (directNeighbours unsymmetricWeb) + @?= [[5],[2,3,0],[4,3],[4,2,5,1],[5],[0,1,6],[5],[4,6]] + , testCase "Next-neighbours in simple quadratic web." + $ toList (nextNeighbours quadraticWeb) @?= + [ [(1,[0,3]),(2,[0,3])] + , [(0,[1,2]),(3,[1,2])] + , [(0,[1,2]),(3,[1,2])] + , [(1,[0,3]),(2,[0,3])] ] + , testCase "Next-neighbours in triangular web (after scrambling)" + $ toList (nextNeighbours $ scrambleKnitting triangularWeb) @?= + [ [(2,[1,0]),(1,[2,0])] + , [(2,[1,0]),(0,[2,1])] + , [(1,[2,0]),(0,[2,1])] ] + , testCase "Layers in a nested web" + $ toList (pointsLocInEnvi nestedWeb) @?= + [ [((1, 朳[(o,朳[ {-LEAF-} (o,朳[]) ])]), 0) + ,((2, 朳[(o,朳[ {- {- -} {- -}-} (o,朳[(o,朳[]),(o,朳[])]) ])]), 0) + ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 0) + ] + , [((1, 朳[(o,朳[ (o,朳[]) {-LEAF-} ])]), 1) + ,((2, 朳[(o,朳[ {- {- -} {- -}-} (o,朳[(o,朳[]),(o,朳[])]) ])]), 0) + ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 0) + ] + , [((1, 朳[(o,朳[ {-LEAF-} (o,朳[]) ])]), 0) + ,((2, 朳[(o,朳[ (o,朳[(o,朳[]),(o,朳[])]) {- {- -} {- -}-} ])]), 2) + ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 0) + ] + , [((1, 朳[(o,朳[ (o,朳[]) {-LEAF-} ])]), 1) + ,((2, 朳[(o,朳[ (o,朳[(o,朳[]),(o,朳[])]) {- {- -} {- -}-} ])]), 2) + ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 0) + ] + , [((1, 朳[(o,朳[ {-LEAF-} (o,朳[]) ])]), 0) + ,((2, 朳[(o,朳[ {- {- -} {- -}-} (o,朳[(o,朳[]),(o,朳[])]) ])]), 0) + ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 4) + ] + , [((1, 朳[(o,朳[ (o,朳[]) {-LEAF-} ])]), 1) + ,((2, 朳[(o,朳[ {- {- -} {- -}-} (o,朳[(o,朳[]),(o,朳[])]) ])]), 0) + ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 4) + ] + , [((1, 朳[(o,朳[ {-LEAF-} (o,朳[]) ])]), 0) + ,((2, 朳[(o,朳[ (o,朳[(o,朳[]),(o,朳[])]) {- {- -} {- -}-} ])]), 2) + ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 4) + ] + , [((1, 朳[(o,朳[ (o,朳[]) {-LEAF-} ])]), 1) + ,((2, 朳[(o,朳[ (o,朳[(o,朳[]),(o,朳[])]) {- {- -} {- -}-} ])]), 2) + ,((4, 朳[(o,朳[(o,朳[(o,朳[(o,朳[]),(o,朳[])]),(o,朳[(o,朳[]),(o,朳[])])])])]), 4) + ] + ] + , testCase "Next-neighbours in nested web." + $ toList (nextNeighbours nestedWeb) @?= + [ [ (1,[0,3,4]), (2,[0,3]) ] + , [ (0,[1,2]) , (3,[1,6]) , (4,[1,5,6]) ] + , [ (0,[1,2]) , (3,[1,6]) ] + , [ (1,[0,3,4]), (6,[3,4,7]) ] + , [ (1,[0,3,4]), (5,[4,7]) , (6,[3,4,7]) ] + , [ (4,[1,5,6]), (7,[5,6]) ] + , [ (3,[1,6]) , (4,[1,5,6]), (7,[5,6]) ] + , [ (5,[4,7]) , (6,[3,4,7]) ] ] + , testCase "Next-neighbours in unsymmetric web." + $ toList (nextNeighbours unsymmetricWeb) @?= + [ [ (5,[0,1,6]) ] + , [ (2,[4,3]) , (3,[4,2,5,1]), (0,[5]) ] + , [ (4,[5]) , (3,[4,2,5,1]) ] + , [ (4,[5]) , (2,[4,3]) , (5,[0,1,6]), (1,[2,3,0]) ] + , [ (5,[0,1,6]) ] + , [ (0,[5]) , (1,[2,3,0]) , (6,[5]) ] + , [ (5,[0,1,6]) ] + , [ (4,[5]) , (6,[5]) ] ] + , testCase "Neighbours in unsymmetric web after scrambling." + $ toList (directNeighbours $ scrambleKnitting unsymmetricWeb) @?= + [ [1,6], [4,3,2,5], [5,4,1], [5,4,0,1,6,2], [0,1,6], [2,3,0], [0,1], [5] ] + ] + , testGroup "Adjacency layers around points in a web" + [ testCase "Onions in nested web" + $ toList (webOnions $ localFmapWeb _thisNodeId nestedWeb) + @?= [ [[(o,0)],[(o,1),(o,2)],[(o,3),(o,4)],[(o,6),(o,5)],[(o,7)]] + , [[(o,1)],[(o,0),(o,3),(o,4)],[(o,6),(o,2),(o,5)],[(o,7)]] + , [[(o,2)],[(o,0),(o,3)],[(o,1),(o,6)],[(o,4),(o,7)],[(o,5)]] + , [[(o,3)],[(o,1),(o,6)],[(o,4),(o,0),(o,7)],[(o,5),(o,2)]] + , [[(o,4)],[(o,1),(o,5),(o,6)],[(o,3),(o,7),(o,0)],[(o,2)]] + , [[(o,5)],[(o,4),(o,7)],[(o,6),(o,1)],[(o,3),(o,0)],[(o,2)]] + , [[(o,6)],[(o,3),(o,4),(o,7)],[(o,1),(o,5)],[(o,0)],[(o,2)]] + , [[(o,7)],[(o,5),(o,6)],[(o,4),(o,3)],[(o,1)],[(o,0)],[(o,2)]] + ] + ] + , testGroup "Neighbour-search for web knitting." + [ testCase "1D line of points" + $ bestNeighbours (euclideanNorm :: Norm ℝ) + (zip [0..] [-1, -0.7 .. 1]) + @?= ([3,4], Nothing) + , testCase "Origin-boundary excluding two points on the x- and y-axes" + $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) + [(0, (1,0)), (1, (0,1))] + @?= ([0,1], Just (sqrt 2/2, sqrt 2/2)) + , testCase "Origin-boundary excluding points in the x≥0 half plane" + $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) + [(0, (1,0)), (1, (0,1)), (2, (0,-1))] + @?= ([0,1,2], Just (1, -1.922877998462862e-16)) + , testCase "Best neighbours in a quadratic surrounding" + $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) + [ (1, (0,-1)), (2, (1,-1)) + , (3, (-1,0)), (4, (1,0)) + , (5, (-1,1)), (6, (0,1)), (7, (1,1)) ] + @?= ([1,3,4,6], Nothing) + , testCase "Best neighbours to the corner of a rectangular grid" + $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) + [ ( 1,(1,0)), ( 2,(2,0)), ( 3,(3,0)) + , (10,(0,1)), (11,(1,1)), (12,(2,1)), (13,(3,1)) + , (20,(0,2)), (21,(1,2)), (22,(2,2)), (23,(3,2)) ] + @?= ([1,10], Just (sqrt 2/2, sqrt 2/2)) + , testCase "Best neighbours in a rectangular grid" + $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) + ((id&&&id) <$> + [ (-2,-1), (-1,-1), ( 0,-1), ( 1,-1), ( 2,-1) + , (-2, 0), (-1, 0),{-ORIGIN-}( 1, 0), ( 2, 0) + , (-2, 1), (-1, 1), ( 0, 1), ( 1, 1), ( 2, 1) ]) + @?= ([(0,-1), (-1,0), (1,0), (0,1)], Nothing) + , testCase "Best neighbours in a big rectangular grid" + $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) + ((id&&&id) <$> + [ (-3,-3), (-2,-3), (-1,-3), ( 0,-3), ( 1,-3), ( 2,-3), ( 3,-3) + , (-3,-2), (-2,-2), (-1,-2), ( 0,-2), ( 1,-2), ( 2,-2), ( 3,-2) + , (-3,-1), (-2,-1), (-1,-1), ( 0,-1), ( 1,-1), ( 2,-1), ( 3,-1) + , (-3, 0), (-2, 0), (-1, 0),{-ORIGIN-}( 1, 0), ( 2, 0), ( 3, 0) + , (-3, 1), (-2, 1), (-1, 1), ( 0, 1), ( 1, 1), ( 2, 1), ( 3, 1) + , (-3, 2), (-2, 2), (-1, 2), ( 0, 2), ( 1, 2), ( 2, 2), ( 3, 2) ]) + @?= ([(0,-1), (-1,0), (1,0), (0,1)], Nothing) + , testCase "Best neighbours in an irregular point-cloud" + $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) + ((id&&&id) <$> + [ (-1,-6) + + , (0,-5), (4,-5),(5,-5) + + , (-1,-4) + + ,(-6,-3), (-4,-3),(2,-3) + + , (-2,-2), (0,-2) + + , (1,-1), (4,-1),(5,-1) + + {-ORIGIN-} + + ,(-5,1), (-3,1),(-2,1), (2,1), (4,1), (5,1) + + , (-3,2),(-2,3), (1,3),(2,3) + + ,(-5,4), (-1,4),(3,4) + + , (-3,5), (3,5) + + , (2,6), (5,6),(6,6) ]) + @?= ([(1,-1), (-2,-2), (2,1), (-6,-3), (-2,1)], Nothing) + , testCase "Best neighbours in degenerate near-boundary constellation" + $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) + ((id &&& (^-^(3.6, 3.0))) <$> reverse + [ (3.15,3.6) + , (3.29,3.4) + , (3.3,3.2), (3.45,3.2), (3.6,3.2) + , (3.15,3.0), (3.3,3.0) {-ORIGIN-} + , (3.15,2.8), (3.3,2.8), (3.45,2.8), (3.6,2.8), (3.75,2.8) + , (3.45,2.6), (3.6,2.6), (3.75,2.6), (3.9,2.6) + , (3.15,2.2) + ]) + @?= ([(3.6,2.8), (3.3,3.0), (3.6,3.2), (3.75,2.8)], Nothing) + , testCase "Best neighbours in point selection from almost-rectangular grid" + $ bestNeighbours (euclideanNorm :: Norm (ℝ,ℝ)) + ([ (235,(0.0,-0.2)) + , (248,(-0.7499999999999996,0.0)) + , (267,(0.0,0.2)) + , (268,(0.15,0.0)) + , (271,(-0.14999,0.0)) + ]) + @?= ([271,267,268,235], Nothing) + , testCase "Best neighbours in point selection of 1D web test" + $ bestNeighbours (euclideanNorm :: Norm ℝ) + ((id &&& (^-^467)) <$> + [ 565.5193483520385, 254.62827644949562 + , 203.3896874080876, 214.87356399193985 ]) + @?= ([565.5193483520385, 254.62827644949562], Nothing) + ] + , testGroup "Automatically building webs" + [ testCase "Minimal, 3-point 1D “web”" + $ let web = fromWebNodes euclideanMetric [(x, ()) | x<-[0,1,2]] + :: PointsWeb ℝ () + in toList (localFmapWeb (\info + -> ( fst <$> info^.nodeNeighbours + , info^.webBoundingPlane ) ) web) + @?= [([1], Just 1), ([0,2], Nothing), ([1], Just $ -1)] + , testCase "Linear 1D “web”" + $ toList (directNeighbours (fromWebNodes euclideanMetric + [(x, ()) | x<-[0, 0.1 .. 2]] :: PointsWeb ℝ () )) + @?= [ [1,9], [0,2], [1,3], [2,4], [3], [6,12], [5,7], [6,8], [7,9], [0,8], [11,15] + , [10,12],[11,5],[14,20],[13,15],[10,14],[17],[16,18],[17,19],[18,20],[13,19] + ] + , testCase "Small linear 1D web with nonuniform spacing" + $ toList (directNeighbours (fromWebNodes euclideanMetric + [ (x, ()) | x<-[ 203.3896874080876 + , 214.87356399193985 + , 254.62827644949562 + , 467.0 + , 565.5193483520385 ] + ] :: PointsWeb ℝ () )) + @?= [ [1], [0,2], [1,3], [4,2], [3] ] + , QC.testProperty "Random 1D web should be strongly connected" + $ \ps -> length ps >= 2 ==> + length (Graph.scc . fst + $ toGraph ( fromWebNodes euclideanMetric + [(x, ()) | x<-Set.toList ps] :: PointsWeb ℝ () ) + ) == 1 + , QC.testProperty "Random 1D web should have only 2 boundary-points" + $ \ps -> length ps >= 2 ==> + length (webBoundary (fromWebNodes euclideanMetric + [(x, ()) | x<-Set.toList ps] :: PointsWeb ℝ () ) + ) == 2 + ] + , testGroup "Shades" + [ testCase "Equality of `Shade`s" + $ (1 :± [1]) @?≈ (1 :± [1] :: Shade ℝ) + , testCase "Equality of `Shade'`s" + $ ((1,0)|±|[(1,-2),(3,4)]) @?≈ ((1,0)|±|[(1,-2),(3,4)] :: Shade' (ℝ,ℝ)) + , testCase "Pragmatically showing" + $ SP.show ((1,0)|±|[(1,-2),(3,4)] :: Shade' (ℝ,ℝ)) + @?= "(1,0)|±|[(5,2),(0,2)]" + , testCase "Pragmatically showing (with orthogonal span)" + $ SP.show ((1,0)|±|[(6,0),(0,2)] :: Shade' (ℝ,ℝ)) + @?= "(1,0)|±|[(6,0),(0,2)]" + ] + , testGroup "Function models for uncertain data" + [ testCase "Fitting a 1D affine model to constant data" + $ fitLocally [ (-1, 5|±|[1]), (0, 5|±|[1]), (1, 5|±|[1]) ] + @?≈ Just ( + AffineModel (5:±[1.15]) (zeroV:±[id^/sqrt 2]) :: AffineModel ℝ ℝ ) + , testCase "Fitting a 2D affine model to constant data" + $ fitLocally [ ((0,1), 5|±|[1]) + , ((-1,0), 5|±|[1]), ((0,0), 5|±|[1]), ((1,0), 5|±|[1]) + , ((0,-1), 5|±|[1]) ] + @?≈ Just ( + AffineModel (5:±[0.9]) (zeroV:±((^/sqrt 2)<$>[fst, snd])) + :: AffineModel (ℝ,ℝ) ℝ ) + , testCase "Fitting a 1D affine model to rising-uncertainty data" + $ fitLocally [ (-1, 3|±|[0.1]), (0, 4|±|[0.5]), (1, 5|±|[1]) ] + @?≈ Just ( + AffineModel (4:±[1/sqrt 2]) (id:±[id^*0.36]) :: AffineModel ℝ ℝ ) + , testCase "Fitting a 1D affine model to quadratic data" + $ fitLocally [ (-1, 3|±|[0.1]), (0, 0|±|[0.1]), (1, 3|±|[0.1]) ] + @?≈ Just ( + AffineModel (2:±[2.94]) (zeroV:±[id^*1.8]) :: AffineModel ℝ ℝ ) + ] + ] + +emptyWeb, singletonWeb, triangularWeb, quadraticWeb, nestedWeb, unsymmetricWeb + , unidirDiagonalLinkedWeb + :: PointsWeb ℝ⁰ () + +emptyWeb = PointsWeb $ PlainLeaves [] + +singletonWeb = PointsWeb $ + PlainLeaves [ (o, Neighbourhood () mempty euclideanNorm Nothing) ] + +triangularWeb = PointsWeb $ + PlainLeaves [ (o, Neighbourhood () [1,2] euclideanNorm Nothing) + , (o, Neighbourhood () [-1,1] euclideanNorm Nothing) + , (o, Neighbourhood () [-2,-1] euclideanNorm Nothing) + ] + +quadraticWeb = PointsWeb $ + OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass + (PlainLeaves [ (o, Neighbourhood () [1,2] euclideanNorm Nothing) + , (o, Neighbourhood () [-1,2] euclideanNorm Nothing) + ]) + (PlainLeaves [ (o, Neighbourhood () [-2,1] euclideanNorm Nothing) + , (o, Neighbourhood () [-2,-1] euclideanNorm Nothing) + ]) + ) + +nestedWeb = PointsWeb $ + OverlappingBranches 8 (Shade o mempty) (pure . DBranch o $ Hourglass + (OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass + (PlainLeaves [ (o, Neighbourhood () [1,2] euclideanNorm Nothing) + , (o, Neighbourhood () [-1,2,3] euclideanNorm Nothing) + ]) + (PlainLeaves [ (o, Neighbourhood () [-2,1] euclideanNorm Nothing) + , (o, Neighbourhood () [-2,3] euclideanNorm Nothing) + ]) + )) + (OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass + (PlainLeaves [ (o, Neighbourhood () [-3,1,2] euclideanNorm Nothing) + , (o, Neighbourhood () [-1,2] euclideanNorm Nothing) + ]) + (PlainLeaves [ (o, Neighbourhood () [-3,-2,1] euclideanNorm Nothing) + , (o, Neighbourhood () [-2,-1] euclideanNorm Nothing) + ]) + )) + ) + +unsymmetricWeb = PointsWeb $ + OverlappingBranches 8 (Shade o mempty) (pure . DBranch o $ Hourglass + (OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass + (PlainLeaves [ (o, Neighbourhood () [5] euclideanNorm Nothing) + , (o, Neighbourhood () [1,2,-1] euclideanNorm Nothing) + ]) + (PlainLeaves [ (o, Neighbourhood () [2,1] euclideanNorm Nothing) + , (o, Neighbourhood () [1,-1,2,-2] euclideanNorm Nothing) + ]) + )) + (OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass + (PlainLeaves [ (o, Neighbourhood () [1] euclideanNorm Nothing) + , (o, Neighbourhood () [-5,-4,1] euclideanNorm Nothing) + ]) + (PlainLeaves [ (o, Neighbourhood () [-1] euclideanNorm Nothing) + , (o, Neighbourhood () [-3,-1] euclideanNorm Nothing) + ]) + )) + ) + +unidirDiagonalLinkedWeb = PointsWeb $ + OverlappingBranches 4 (Shade o mempty) (pure . DBranch o $ Hourglass + (PlainLeaves [ (o, Neighbourhood () [1,2,3] euclideanNorm Nothing) + , (o, Neighbourhood () [-1,2] euclideanNorm Nothing) + ]) + (PlainLeaves [ (o, Neighbourhood () [-2,-1,1] euclideanNorm Nothing) + , (o, Neighbourhood () [-2,-1] euclideanNorm Nothing) + ]) + ) + + + +o = zeroV :: ℝ⁰ + +dummyWebFmap :: PointsWeb ℝ⁰ a -> PointsWeb ℝ⁰ a +dummyWebFmap = localFmapWeb $ \info -> info^.thisNodeData + +directNeighbours :: WithField ℝ Manifold v => PointsWeb v () -> PointsWeb v [WebNodeId] +directNeighbours = localFmapWeb $ + \info -> fst <$> info^.nodeNeighbours + +nextNeighbours :: PointsWeb ℝ⁰ a -> PointsWeb ℝ⁰ [(WebNodeId, [WebNodeId])] +nextNeighbours = webLocalInfo >>> localFmapWeb `id` + \info -> [ ( nId ≡! nId' ≡! (nInfo^.thisNodeId) ≡! (nInfo'^.thisNodeId) + , (fst<$>nInfo^.nodeNeighbours) ≡! (fst<$>nInfo'^.nodeNeighbours) ) + | ((nId,(_,nInfo)),(nId',(_,nInfo'))) + <- zip (info^.nodeNeighbours) + (info^.thisNodeData.nodeNeighbours) + , all (==Origin) [ nInfo''^.thisNodeCoord + | (_,(_,nInfo''))<-nInfo'^.nodeNeighbours ] + ] + +pointsLocInEnvi :: PointsWeb ℝ⁰ a -> PointsWeb ℝ⁰ [((Int, Trees ℝ⁰), WebNodeId)] +pointsLocInEnvi = fmapNodesInEnvi $ + \(NodeInWeb (_, orig) env) + -> fmap (const $ first ((nLeaves&&&onlyNodes) . fmap (const ())) <$> env) orig + + +scrambleKnitting :: PointsWeb ℝ⁰ a -> PointsWeb ℝ⁰ a +scrambleKnitting = tweakWebGeometry euclideanMetric + $ \info -> nub [ i' + | (_, (_, nInfo)) <- info^.nodeNeighbours + , (i',_) <- nInfo^.nodeNeighbours + , i' /= info^.thisNodeId ] + +infixl 4 ≡! +(≡!) :: (Eq a, Show a) => a -> a -> a +x ≡! y | x==y = x + | otherwise = error $ show x++" ≠ "++show y + + +infix 4 ≈ +class AEq e where + (≈) :: e -> e -> Bool +instance (SimpleSpace v, Needle v~v, Interior v~v, Floating (Scalar v)) + => AEq (Shade' v) where + Shade' c₀ σ₀ ≈ Shade' c₁ σ₁ + = (σ₀|$|δ) < ε && (σ₀|$|δ) < ε + && all (is1 . (σ₀|$|)) (normSpanningSystem' σ₁) + && all (is1 . (σ₁|$|)) (normSpanningSystem' σ₀) + where δ = c₁ ^-^ c₀ + ε = 1e-2 + is1 x = abs (x-1) < ε +instance ( SimpleSpace v, DualVector (Needle' v) ~ v, Interior v ~ v + , InnerSpace (Scalar v), Scalar (Needle' v) ~ Scalar v ) + => AEq (Shade v) where + Shade c₀ σ₀ ≈ Shade c₁ σ₁ + = (dualNorm σ₀|$|δ) < ε && (dualNorm σ₀|$|δ) < ε + && all (is1 . (dualNorm σ₀|$|)) (normSpanningSystem σ₁) + && all (is1 . (dualNorm σ₁|$|)) (normSpanningSystem σ₀) + where δ = c₁ ^-^ c₀ + ε = 1e-2 + is1 x = abs (x-1) < ε +instance AEq a => AEq (Maybe a) where + Just x ≈ Just y = x ≈ y + Nothing ≈ Nothing = True + _ ≈ _ = False +instance (AEq (Shade y), AEq (Shade (Needle x +> Needle y))) + => AEq (AffineModel x y) where + AffineModel b₀ a₀ ≈ AffineModel b₁ a₁ = b₀ ≈ b₁ && a₀ ≈ a₁ + +infix 1 @?≈ +(@?≈) :: (AEq e, Show e) => e -> e -> Assertion +a@?≈b + | a≈b = return () + | otherwise = assertFailure $ "Expected "++show b++", but got "++show a + +