-- | -- Module : Data.Manifold.Web -- Copyright : (c) Justus Sagemüller 2016 -- License : GPL v3 -- -- Maintainer : (@) sagemueller $ geo.uni-koeln.de -- Stability : experimental -- Portability : portable -- {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE DeriveFunctor #-} {-# LANGUAGE DeriveFoldable #-} {-# LANGUAGE DeriveTraversable #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE TupleSections #-} {-# LANGUAGE ParallelListComp #-} {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE PatternGuards #-} {-# LANGUAGE LambdaCase #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE LiberalTypeSynonyms #-} {-# LANGUAGE TemplateHaskell #-} module Data.Manifold.Web ( -- * The web data type PointsWeb -- ** Construction , fromWebNodes, fromShadeTree_auto, fromShadeTree, fromShaded -- ** Lookup , nearestNeighbour, indexWeb, toGraph, webBoundary -- ** Decomposition , sliceWeb_lin -- , sampleWebAlongLine_lin , sampleWeb_2Dcartesian_lin, sampleEntireWeb_2Dcartesian_lin -- ** Local environments , localFocusWeb -- * Uncertain functions , differentiateUncertainWebFunction, differentiate²UncertainWebFunction , localModels_CGrid -- * Differential equations -- ** Fixed resolution , iterateFilterDEqn_static, iterateFilterDEqn_pathwise , iterateFilterDEqn_static_selective -- ** Automatic resolution , filterDEqnSolutions_adaptive, iterateFilterDEqn_adaptive -- ** Configuration , InconsistencyStrategy(..) , InformationMergeStrategy(..) , naïve, inconsistencyAware, indicateInconsistencies, postponeInconsistencies , PropagationInconsistency(..) -- * Misc , ConvexSet(..), ellipsoid, ellipsoidSet, coerceWebDomain , rescanPDELocally, localOnion, webOnions, knitShortcuts ) where import Data.List hiding (filter, all, foldr1) import Data.Maybe import qualified Data.Map as Map import qualified Data.Vector as Arr import qualified Data.Vector.Mutable as MArr import qualified Data.Vector.Unboxed as UArr import qualified Data.Array as PArr import Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NE import Data.List.FastNub (fastNub,fastNubBy) import Data.Ord (comparing) import Data.Semigroup import Control.DeepSeq import Data.VectorSpace import Math.LinearMap.Category hiding (trace) import Data.Tagged import Data.Function (on) import Data.Fixed (mod') import Data.Manifold.Types import Data.Manifold.Types.Primitive import Data.Manifold.PseudoAffine import Data.Manifold.Shade import Data.Manifold.TreeCover import Data.SetLike.Intersection import Data.Manifold.Riemannian import Data.Manifold.Atlas import Data.Manifold.Function.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) import qualified Control.Applicative as Hask import qualified Control.Monad as Hask hiding(forM_, sequence) import Control.Monad.ST (runST) import Data.STRef (newSTRef, modifySTRef, readSTRef) 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) import qualified Data.Traversable as Hask import Data.Traversable (forM) import Data.Graph import Control.Category.Constrained.Prelude hiding ((^), all, elem, sum, forM, Foldable(..), foldr1, Traversable, traverse) import Control.Arrow.Constrained import Control.Monad.Constrained hiding (forM) import Data.Foldable.Constrained import Data.Traversable.Constrained (Traversable, traverse) import Control.Comonad (Comonad(..)) import Control.Comonad.Cofree import Control.Lens ((&), (%~), (^.), (.~), (+~), ix, iover, indexing) import Control.Lens.TH import GHC.Generics (Generic) import Development.Placeholders unlinkedFromWebNodes :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (MetricChoice x) -> [(x,y)] -> PointsWeb x y unlinkedFromWebNodes = case boundarylessWitness :: BoundarylessWitness x of BoundarylessWitness -> \mf -> unlinkedFromShaded mf . fromLeafPoints_ 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_ 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 regroup' 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)) => (Shade x -> Metric x) -> ShadeTree x -> PointsWeb x () fromShadeTree mf = fromShaded mf . constShaded () fromShaded :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => (MetricChoice x) -- ^ Local scalar-product generator. You can always -- use @'recipMetric' . '_shadeExpanse'@ (but this -- may give distortions compared to an actual -- Riemannian metric). -> (x`Shaded`y) -- ^ Source tree. -> PointsWeb x y fromShaded metricf = knitShortcuts metricf . autoLinkWeb . unlinkedFromShaded metricf toShaded :: WithField ℝ PseudoAffine x => PointsWeb x y -> (x`Shaded`y) 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 = $notImplemented -- | 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 = knitShortcuts type OSNeedle x = (Needle' x, Needle x) type OSNode x y = (OSNeedle x, WebLocally x y) type CPCone x = (Needle' x, ℝ) -- | 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)) 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 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 { _interpolationSegRange :: (ℝ,ℝ) , _interpolationFunction :: ℝ -> y } type InterpolationSeq y = [InterpolationIv y] mkInterpolationSeq_lin :: (x~ℝ, Geodesic y) => [(x,y)] -> InterpolationSeq y mkInterpolationSeq_lin [(xψ,yψ), (xω,yω)] = return $ InterpolationIv (xψ,xω) (\x -> let drel = fromIntv0to1 $ (x-xψ)/(xω-xψ) in yio drel ) where Just yio = geodesicBetween yψ yω mkInterpolationSeq_lin (p₀:p₁:ps) = mkInterpolationSeq_lin [p₀,p₁] <> mkInterpolationSeq_lin (p₁:ps) mkInterpolationSeq_lin _ = [] -- | Fetch a point between any two neighbouring web nodes on opposite -- sides of the plane, and linearly interpolate the values onto the -- cut plane. sliceWeb_lin :: ∀ x y . ( WithField ℝ Manifold x, SimpleSpace (Needle x) , Geodesic x, Geodesic y ) => PointsWeb x y -> Cutplane x -> [(x,y)] sliceWeb_lin web = sliceEdgs 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₁)] , Just xi <- [geodesicBetween x₀ x₁] , Just yi <- [geodesicBetween y₀ y₁] ] data GridPlanes x = GridPlanes { _gridPlaneNormal :: Needle' x , _gridPlaneSpacing :: Needle x , _gridPlanesCount :: Int } deriving instance (Show x, Show (Needle x), Show (Needle' x)) => Show (GridPlanes x) data GridSetup x = GridSetup { _gridStartCorner :: x , _gridSplitDirs :: [GridPlanes x] } deriving instance (Show x, Show (Needle x), Show (Needle' x)) => Show (GridSetup x) cartesianGrid2D :: (x~ℝ, y~ℝ) => ((x,x), Int) -> ((y,y), Int) -> GridSetup (x,y) cartesianGrid2D ((x₀,x₁), nx) ((y₀,y₁), ny) = GridSetup (x₀+dx/2, y₀+dy/2) [ GridPlanes (0,1) (0, dy) ny, GridPlanes (1,0) (dx, 0) nx ] where dx = (x₁-x₀)/fromIntegral nx dy = (y₁-y₀)/fromIntegral ny splitToGridLines :: ( WithField ℝ Manifold x, SimpleSpace (Needle x) , Geodesic x, Geodesic y ) => PointsWeb x y -> GridSetup x -> [((x, GridPlanes x), [(x,y)])] splitToGridLines web (GridSetup x₀ [GridPlanes dirΩ spcΩ nΩ, linePln]) = [ ((x₀', linePln), sliceWeb_lin web $ Cutplane x₀' (Stiefel1 dirΩ)) | k <- [0 .. nΩ-1] , let x₀' = x₀i.+~^(fromIntegral k *^ spcΩ) ] where Just x₀i = toInterior x₀ sampleWebAlongGrid_lin :: ∀ x y . ( WithField ℝ Manifold x, SimpleSpace (Needle x) , Geodesic x, Geodesic y ) => PointsWeb x y -> GridSetup x -> [(x,Maybe y)] sampleWebAlongGrid_lin web grid = finalLine boundarylessWitness =<< splitToGridLines web grid where finalLine :: BoundarylessWitness x -> ((x, GridPlanes x), [(x,y)]) -> [(x,Maybe y)] finalLine BoundarylessWitness ((x₀, GridPlanes _ dir nSpl), verts) | length verts < 2 = take nSpl $ (,empty)<$>iterate (.+~^dir) x₀ finalLine BoundarylessWitness ((x₀, GridPlanes dx dir nSpl), verts) = take nSpl $ go (x₀,0) intpseq where intpseq = mkInterpolationSeq_lin $ sortBy (comparing fst) [ (dx <.>^ (x.-~!x₀), y) | (x,y) <- verts ] go (x,_) [] = (,empty)<$>iterate (.+~^dir) x go xt (InterpolationIv (tb,te) f:fs) = case span (( [ (x, if t^dir sampleWeb_2Dcartesian_lin :: (x~ℝ, y~ℝ, Geodesic z) => PointsWeb (x,y) z -> ((x,x),Int) -> ((y,y),Int) -> [(y,[(x,Maybe z)])] sampleWeb_2Dcartesian_lin web (xspec@(_,nx)) yspec = go . sampleWebAlongGrid_lin web $ cartesianGrid2D xspec yspec where go [] = [] go l@(((_,y),_):_) = let (ln,l') = splitAt nx l in (y, map (\((x,_),z) -> (x,z)) ln) : go l' sampleEntireWeb_2Dcartesian_lin :: (x~ℝ, y~ℝ, Geodesic z) => PointsWeb (x,y) z -> Int -> Int -> [(y,[(x,Maybe z)])] sampleEntireWeb_2Dcartesian_lin web nx ny = sampleWeb_2Dcartesian_lin web ((x₀,x₁),nx) ((y₀,y₁),ny) where x₀ = minimum (fst<$>pts) x₁ = maximum (fst<$>pts) y₀ = minimum (snd<$>pts) y₁ = maximum (snd<$>pts) pts = fst . fst <$> toList (localFocusWeb web) 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 = 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 -> [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 = ( 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 && 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 <<< snd) . (`localOnion`[])) nearestNeighbour :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => PointsWeb x y -> x -> Maybe (x,y) 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 -- ^ 'fmap' from the co-Kleisli category of 'WebLocally'. 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 -- ^ '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 . ( ModellableRelation x y ) => WebLocally x (Shade' y) -> Shade' (LocalLinear x y) 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 . ( ModellableRelation x y ) => PointsWeb x (Shade' y) -> PointsWeb x (Shade' (LocalLinear x y)) differentiateUncertainWebFunction = localFmapWeb differentiateUncertainWebLocally differentiate²UncertainWebLocally :: ∀ x y . ( ModellableRelation x y ) => WebLocally x (Shade' y) -> Shade' (Needle x ⊗〃+> Needle y) differentiate²UncertainWebLocally = d²uwl ( dualSpaceWitness :: DualSpaceWitness (Needle x) , dualSpaceWitness :: DualSpaceWitness (Needle y) ) where d²uwl (DualSpaceWitness, DualSpaceWitness) info = case fitLocally $ (\(δx,ngb) -> (δx, ngb^.thisNodeData) ) <$> (zeroV,info) : envi of 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 , dualSpaceWitness :: DualSpaceWitness (Needle v) , dualSpaceWitness :: DualSpaceWitness (Needle y) ) of (LinearManifoldWitness BoundarylessWitness, DualSpaceWitness, DualSpaceWitness) -> const zeroV &&& id differentiate²UncertainWebFunction :: ∀ x y . ( ModellableRelation x y ) => PointsWeb x (Shade' y) -> PointsWeb x (Shade' (Needle x ⊗〃+> Needle y)) differentiate²UncertainWebFunction = localFmapWeb differentiate²UncertainWebLocally 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 -> 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 -> fst ( rescan $ case fitLocally $ map (id *** _thisNodeData) =<< (localOnion info []) of Just ㄇ -> ㄇ) >>= intersectShade's . (:|[info^.thisNodeData]) fromGraph :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x)) => MetricChoice x -> Graph -> (Vertex -> (x, y)) -> PointsWeb x y fromGraph metricf gr dataLookup = introduceLinks $ unlinkedFromWebNodes metricf [(fst (dataLookup v), v) | v <- vertices gr] where introduceLinks :: PointsWeb x Vertex -> PointsWeb x y introduceLinks (PointsWeb w) = PointsWeb $ iover (indexing Hask.traverse) (\wi (Neighbourhood vert _ sclPr bound) -> let neighbours = gr PArr.! wi neighbourwis = (vertToWebNode Map.!) <$> neighbours (x, y) = dataLookup vert in Neighbourhood y (UArr.fromList $ subtract wi<$>neighbourwis) sclPr (snd (bestNeighbours sclPr [ ((), fst (dataLookup ni).-~!x) | ni<-neighbours ])) ) w where webNodeToVert = Map.fromList assocs vertToWebNode = Map.fromList $ swap<$>assocs assocs = zip [0..] [vert | Neighbourhood vert _ _ _ <- toList w] toGraph :: (WithField ℝ Manifold x, SimpleSpace (Needle x)) => PointsWeb x y -> (Graph, Vertex -> (x, y)) toGraph wb = second (>>> \(i,_,_) -> case indexWeb wb i of {Just xy -> xy}) (graphFromEdges' edgs) where edgs :: [(Int, Int, [Int])] edgs = Arr.toList . Arr.imap (\i (Neighbourhood _ ngbs _ _) -> (i, i, (i+) <$> UArr.toList ngbs)) . Arr.fromList . Hask.toList $ webNodeRsc wb data ConvexSet x = EmptyConvex | ConvexSet { convexSetHull :: Shade' x -- ^ If @p@ is in all intersectors, it must also be in the hull. , convexSetIntersectors :: [Shade' x] } deriving instance LtdErrorShow x => Show (ConvexSet x) ellipsoid :: Shade' x -> ConvexSet x ellipsoid s = ConvexSet s [s] ellipsoidSet :: Embedding (->) (Maybe (Shade' x)) (ConvexSet x) ellipsoidSet = Embedding (\case {Just s -> ConvexSet s [s]; Nothing -> EmptyConvex}) (\case {ConvexSet h _ -> Just h; EmptyConvex -> Nothing}) intersectors :: ConvexSet x -> Maybe (NonEmpty (Shade' x)) intersectors (ConvexSet h []) = pure (h:|[]) intersectors (ConvexSet _ (i:sts)) = pure (i:|sts) intersectors _ = empty -- | Under intersection. instance Refinable x => Semigroup (ConvexSet x) where a<>b = sconcat (a:|[b]) sconcat csets | Just allIntersectors <- sconcat <$> Hask.traverse intersectors csets , IntersectT ists <- rmTautologyIntersect perfectRefine $ IntersectT allIntersectors , Just hull' <- intersectShade's ists = ConvexSet hull' (NE.toList ists) | otherwise = EmptyConvex where perfectRefine sh₁ sh₂ | sh₁`subShade'`sh₂ = pure sh₁ | sh₂`subShade'`sh₁ = pure sh₂ | otherwise = empty itWhileJust :: InconsistencyStrategy m x y -> (a -> m a) -> a -> [a] itWhileJust AbortOnInconsistency f x | Just y <- f x = x : itWhileJust AbortOnInconsistency f y itWhileJust IgnoreInconsistencies f x | Identity y <- f x = x : itWhileJust IgnoreInconsistencies f y itWhileJust (HighlightInconsistencies yh) f x | Identity y <- f x = x : itWhileJust (HighlightInconsistencies yh) f y itWhileJust _ _ x = [x] dupHead :: NonEmpty a -> NonEmpty a dupHead (x:|xs) = x:|x:xs newtype InformationMergeStrategy n m y' y = InformationMergeStrategy { mergeInformation :: y -> n y' -> m y } naïve :: (NonEmpty y -> y) -> InformationMergeStrategy [] Identity (x,y) y naïve merge = InformationMergeStrategy (\o n -> Identity . merge $ o :| fmap snd n) inconsistencyAware :: (NonEmpty y -> m y) -> InformationMergeStrategy [] m (x,y) y inconsistencyAware merge = InformationMergeStrategy (\o n -> merge $ o :| fmap snd n) indicateInconsistencies :: (NonEmpty υ -> Maybe υ) -> InformationMergeStrategy [] (Except (PropagationInconsistency x υ)) (x,υ) υ indicateInconsistencies merge = InformationMergeStrategy (\o n -> case merge $ o :| fmap snd n of Just r -> pure r Nothing -> throwE $ PropagationInconsistency n o ) 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 data InconsistencyStrategy m x y where AbortOnInconsistency :: InconsistencyStrategy Maybe x y IgnoreInconsistencies :: InconsistencyStrategy Identity x y HighlightInconsistencies :: y -> InconsistencyStrategy Identity x y deriving instance Hask.Functor (InconsistencyStrategy m x) iterateFilterDEqn_static :: ( ModellableRelation x y, Hask.MonadPlus m, LocalModel ㄇ ) => InformationMergeStrategy [] m (x,Shade' y) iy -> Embedding (->) (Shade' y) iy -> DifferentialEqn ㄇ x y -> PointsWeb x (Shade' y) -> Cofree m (PointsWeb x (Shade' y)) iterateFilterDEqn_static strategy shading f = fmap (fmap (shading >-$)) . coiter (filterDEqnSolutions_static strategy shading f) . fmap (shading $->) iterateFilterDEqn_pathwise :: ( ModellableRelation x y, Hask.MonadPlus m, Hask.Traversable m, LocalModel ㄇ ) => InformationMergeStrategy [] m (x,Shade' y) iy -> Embedding (->) (Shade' y) iy -> DifferentialEqn ㄇ x y -> PointsWeb x (Shade' y) -> Cofree m (PointsWeb x (Shade' y)) iterateFilterDEqn_pathwise strategy shading f = fmap (fmap (shading >-$)) . (`evalState`7438) . unfoldM (\oldWeb -> do r <- get let i = r `mod` nLeaves (webNodeRsc oldWeb) m = 2^31 - 1 a = 963345 :: Int -- revised Park-Miller put $ (a*r)`mod`m return ( oldWeb , filterDEqnSolutions_static strategy shading f =<) 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) 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 if isJust $ me ^. webBoundingPlane then return oldValue else case updShy of Just shy -> case ngbs of [] -> pure oldValue _:_ | BoundarylessWitness <- (boundarylessWitness::BoundarylessWitness x) -> 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 (fmap (second ((shading>-$) . _thisNodeData)) . concat . tail $ localOnion ngbInfo [me^.thisNodeId]) ) | (δx, (ngbInfo,sj)) <- ngbs ] >>= mergeInformation strategy (shading$->shy) _ -> mergeInformation strategy oldValue empty ) filterDEqnSolutions_pathsTowards :: ∀ x y ㄇ iy m . ( ModellableRelation x y, Hask.MonadPlus m, LocalModel ㄇ ) => WebNodeId -> InformationMergeStrategy [] m (x,Shade' y) iy -> Embedding (->) (Shade' y) iy -> DifferentialEqn ㄇ x y -> PointsWeb x iy -> m (PointsWeb x iy) filterDEqnSolutions_pathsTowards = case ( geodesicWitness :: GeodesicWitness y , boundarylessWitness :: BoundarylessWitness x ) of (GeodesicWitness _, BoundarylessWitness) -> \targetNode strategy shading f -> traversePathsTowards targetNode (\(PathStep stepStart stepEnd) -> StateT $ \odeState -> let apriori = shading >-$ stepEnd^.thisNodeData in case propagateDEqnSolution_loc f (LocalDataPropPlan{ _sourcePosition = stepStart^.thisNodeCoord , _targetPosOffset = (stepEnd^.thisNodeCoord) .-~! (stepStart^.thisNodeCoord) , _sourceData = odeState , _targetAPrioriData = apriori , _relatedData = (fmap (second ((shading>-$) . _thisNodeData)) . concat . tail $ localOnion stepStart [stepEnd^.thisNodeId]) }) of Nothing -> undefined <$> mergeInformation strategy (stepEnd^.thisNodeData) [] Just propd -> (, propd) <$> mergeInformation strategy (stepEnd^.thisNodeData) [ ( stepEnd^.thisNodeCoord, apriori ) , ( stepStart^.thisNodeCoord, propd ) ] ) (\startPoint pathTrav -> evalStateT pathTrav $ shading >-$ startPoint^.thisNodeData) 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 handleInconsistency IgnoreInconsistencies b _ = Identity b handleInconsistency (HighlightInconsistencies _) _ (Just v) = Identity v handleInconsistency (HighlightInconsistencies b) _ _ = Identity b data SolverNodeState x y = SolverNodeInfo { _solverNodeStatus :: ConvexSet y , _solverNodeJacobian :: Shade' (LocalLinear x y) , _solverNodeBadness :: ℝ , _solverNodeAge :: Int } makeLenses ''SolverNodeState type OldAndNew d = (Maybe d, [d]) oldAndNew :: OldAndNew d -> [d] oldAndNew (Just x, l) = x : l oldAndNew (_, l) = l oldAndNew' :: OldAndNew d -> [(Bool, d)] oldAndNew' (Just x, l) = (True, x) : fmap (False,) l oldAndNew' (_, l) = (False,) <$> l 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 -> (x -> Shade' y -> badness) -> PointsWeb x (SolverNodeState x y) -> m (PointsWeb x (SolverNodeState x y)) filterDEqnSolutions_adaptive mf strategy f badness' oldState = fmap recomputeJacobian $ filterGo boundarylessWitness geodesicWitness =<< tryPreproc boundarylessWitness geodesicWitness where tryPreproc :: BoundarylessWitness x -> GeodesicWitness y -> m (PointsWeb x ( (WebLocally x (SolverNodeState x y) , [(Shade' y, badness)]) )) tryPreproc BoundarylessWitness (GeodesicWitness _) = Hask.traverse addPropagation $ webLocalInfo oldState where addPropagation wl | null neighbourInfo = pure (wl, []) | otherwise = (wl,) . map (id&&&badness undefined) <$> propFromNgbs where propFromNgbs :: m [Shade' y] propFromNgbs = mapM (handleInconsistency strategy thisShy) [ propagateDEqnSolution_loc f (LocalDataPropPlan (neigh^.thisNodeCoord) (negateV δx) (convexSetHull $ neigh^.thisNodeData .solverNodeStatus) thisShy [ second (convexSetHull . _solverNodeStatus . _thisNodeData) nn | (_,nn)<-neigh^.nodeNeighbours ] ) | (δx, neigh) <- neighbourInfo ] -- ( (thisPos, thisShy), NE.fromList neighbourHulls ) thisPos = _thisNodeCoord wl :: x thisShy = convexSetHull . _solverNodeStatus $ _thisNodeData wl neighbourInfo = snd <$> _nodeNeighbours wl totalAge = maximum $ _solverNodeAge <$> oldState errTgtModulation = (1-) . (`mod'`1) . negate . sqrt $ fromIntegral totalAge badness x = badness' x . (shadeNarrowness %~ (scaleNorm errTgtModulation)) filterGo :: BoundarylessWitness x -> GeodesicWitness y -> (PointsWeb x ( (WebLocally x (SolverNodeState x y) , [(Shade' y, badness)]) )) -> m (PointsWeb x (SolverNodeState x y)) filterGo BoundarylessWitness (GeodesicWitness _) preproc'd = fmap (smoothenWebTopology mf . fromTopWebNodes mf . concat . fmap retraceBonds . Hask.toList . webLocalInfo . webLocalInfo) $ Hask.traverse (uncurry localChange) preproc'd where smallBadnessGradient, largeBadnessGradient :: ℝ (smallBadnessGradient, largeBadnessGradient) = ( badnessGradRated!!(n`div`4), badnessGradRated!!(n*3`div`4) ) where n = case length badnessGradRated of 0 -> error "No statistics available for badness-grading." l -> l badnessGradRated :: [badness] badnessGradRated = sort [ ngBad / bad | ( LocalWebInfo { _thisNodeData = SolverNodeInfo _ _ bad _ , _nodeNeighbours=ngbs } , ngbProps) <- Hask.toList preproc'd , (_, ngBad) <- ngbProps , ngBad>bad ] localChange :: WebLocally x (SolverNodeState x y) -> [(Shade' y, badness)] -> m (OldAndNew (x, SolverNodeState x y)) localChange localInfo@LocalWebInfo{ _thisNodeCoord = x , _thisNodeData = SolverNodeInfo shy@(ConvexSet hull _) prevJacobi prevBadness age , _nodeNeighbours = ngbs } ngbProps | null ngbs = return ( pure (x, SolverNodeInfo shy prevJacobi prevBadness (age+1)) , [] ) | otherwise = do let (environAge, unfreshness) = maximum&&&minimum $ age : (_solverNodeAge . _thisNodeData . snd . snd <$> ngbs) case find (\(_, badnessN) -> badnessN / prevBadness > smallBadnessGradient) $ ngbProps of Nothing | age < environAge -- point is an obsolete step-stone; -> return (empty,empty) -- do not further use it. _otherwise -> do shy' <- handleInconsistency (ellipsoid<$>strategy) shy $ ((shy<>) . ellipsoid) <$> intersectShade's (fst <$> NE.fromList ngbProps) newBadness <- handleInconsistency (badness x<$>strategy) prevBadness $ case shy' of EmptyConvex -> empty ConvexSet hull' _ -> return $ badness x hull' let updatedNode = SolverNodeInfo shy' prevJacobi newBadness (age+1) stepStones <- if unfreshness < 3 then return [] else fmap concat . forM (zip (second _thisNodeData.snd<$>ngbs) ngbProps) $ \( (vN, SolverNodeInfo (ConvexSet hullN _) _ _ ageN) , (_, nBadnessProp'd) ) -> do case ageN of _ | ageN > 0 , badnessGrad <- nBadnessProp'd / prevBadness , badnessGrad > largeBadnessGradient -> do let stepV = vN^/2 xStep = x .+~^ stepV aprioriInterpolate :: Shade' y Just aprioriInterpolate = middleBetween hull hullN case intersectShade's =<< (sequenceA $ NE.fromList [ propagateDEqnSolution_loc f (LocalDataPropPlan (n^.thisNodeCoord) (stepV ^-^ δx) (convexSetHull $ n^.thisNodeData.solverNodeStatus) aprioriInterpolate (second (convexSetHull ._solverNodeStatus ._thisNodeData) . snd <$> n^.nodeNeighbours) ) -- ( (xStep, hull) -- , NE.cons (negateV stepV, hull) -- $ fmap (\(vN',hullN') -- -> (vN'^-^stepV, hullN') ) ) | (_, (δx, n)) <- ngbs ]) of Just shyStep -> return [( xStep , SolverNodeInfo (ellipsoid shyStep) prevJacobi (badness xStep shyStep) 1 )] _ -> return [] _otherwise -> return [] let updated = (x, updatedNode) return $ (pure updated, stepStones) retraceBonds :: WebLocally x (WebLocally x (OldAndNew (x, SolverNodeState x y))) -> [((x, [Int+Needle x]), SolverNodeState x y)] retraceBonds locWeb@LocalWebInfo{ _thisNodeId = myId , _thisNodeCoord = xOld , _nodeLocalScalarProduct = locMetr } = [ ( (x, Right . fst<$>neighbourCandidates), snsy) | (isOld, (x, snsy)) <- focused , let neighbourCandidates = [ (v,nnId) | (_,ngb) <- knownNgbs , (Just v, nnId) <- case oldAndNew $ ngb^.thisNodeData of [] -> [ (xN.-~.x, nnId) | (nnId, (_,nnWeb)) <- ngb^.nodeNeighbours , nnId /= myId , (xN,_) <- oldAndNew $ nnWeb^.thisNodeData ] l -> [(xN.-~.x, ngb^.thisNodeId) | (xN,_) <- l] ] possibleConflicts = [ normSq locMetr v | (v,nnId)<-neighbourCandidates , nnId > myId ] , isOld || null possibleConflicts || minimum possibleConflicts > oldMinDistSq / 4 ] where focused = oldAndNew' $ locWeb^.thisNodeData^.thisNodeData knownNgbs = second _thisNodeData . snd <$> locWeb^.nodeNeighbours oldMinDistSq = minimum [ normSq locMetr vOld | (_,ngb) <- knownNgbs , let Just vOld = ngb^.thisNodeCoord .-~. xOld ] recomputeJacobian :: ( ModellableRelation x y ) => PointsWeb x (SolverNodeState x y) -> PointsWeb x (SolverNodeState x y) recomputeJacobian = webLocalInfo >>> fmap ((_thisNodeData &&& differentiateUncertainWebLocally . fmap (convexSetHull . _solverNodeStatus)) >>> \(nst, shj) -> nst & solverNodeJacobian .~ shj ) iterateFilterDEqn_adaptive :: ( 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 -> (x -> Shade' y -> ℝ) -- ^ Badness function for local results. -> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)] iterateFilterDEqn_adaptive mf strategy f badness = map (fmap (convexSetHull . _solverNodeStatus)) . itWhileJust strategy (filterDEqnSolutions_adaptive mf strategy f badness) . recomputeJacobian . fmap (\((x,shy),_) -> SolverNodeInfo (ellipsoid shy) (Shade' zeroV mempty) (badness x shy) 1 ) . localFocusWeb