module Data.Manifold.Web (
PointsWeb
, fromWebNodes, fromShadeTree_auto, fromShadeTree, fromShaded
, nearestNeighbour, indexWeb, webEdges, toGraph
, sliceWeb_lin
, localFocusWeb
, filterDEqnSolution_static, iterateFilterDEqn_static
, filterDEqnSolutions_adaptive, iterateFilterDEqn_adaptive
, ConvexSet(..), ellipsoid
) where
import Data.List hiding (filter, all, elem, sum, foldr1)
import Data.Maybe
import qualified Data.Set as Set
import qualified Data.Vector as Arr
import qualified Data.Vector.Unboxed as UArr
import Data.List.NonEmpty (NonEmpty(..))
import qualified Data.List.NonEmpty as NE
import Data.List.FastNub (fastNubBy)
import Data.Ord (comparing)
import Data.Semigroup
import Control.DeepSeq
import Data.VectorSpace
import Data.LinearMap.HerMetric
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.TreeCover
import Data.SetLike.Intersection
import Data.Manifold.Riemannian
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.Trans.State
import Control.Monad.Trans.List
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 Lens.Micro ((&), (%~), (^.), (.~))
import Lens.Micro.TH
import GHC.Generics (Generic)
type WebNodeId = Int
data Neighbourhood x = Neighbourhood {
neighbours :: UArr.Vector WebNodeId
, localScalarProduct :: Metric x
}
deriving (Generic)
instance (NFData x, NFData (HerMetric (Needle x))) => NFData (Neighbourhood x)
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 (HerMetric (Needle 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
fromWebNodes :: ∀ x y . WithField ℝ Manifold x
=> (MetricChoice x) -> [(x,y)] -> PointsWeb x y
fromWebNodes mf = fromShaded mf . fromLeafPoints . map (uncurry WithAny . swap)
fromTopWebNodes :: ∀ x y . WithField ℝ Manifold x
=> (MetricChoice x) -> [((x,[Needle x]),y)] -> PointsWeb x y
fromTopWebNodes mf = fromTopShaded mf . fromLeafPoints
. map (uncurry WithAny . swap . regroup')
fromShadeTree_auto :: ∀ x . WithField ℝ Manifold x => ShadeTree x -> PointsWeb x ()
fromShadeTree_auto = fromShaded (recipMetric . _shadeExpanse) . constShaded ()
fromShadeTree :: ∀ x . WithField ℝ Manifold x
=> (Shade x -> Metric x) -> ShadeTree x -> PointsWeb x ()
fromShadeTree mf = fromShaded mf . constShaded ()
fromShaded :: ∀ x y . WithField ℝ Manifold x
=> (MetricChoice x)
-> (x`Shaded`y)
-> PointsWeb x y
fromShaded metricf = fromTopShaded metricf . fmapShaded ([],)
fromTopShaded :: ∀ x y . WithField ℝ Manifold x
=> (MetricChoice x)
-> (x`Shaded`([Needle x], y))
-> PointsWeb x y
fromTopShaded metricf shd = PointsWeb shd' assocData
where shd' = stripShadedUntopological shd
assocData = Hask.foldMap locMesh $ twigsWithEnvirons shd
locMesh :: ( (Int, ShadeTree (x`WithAny`([Needle x], y)))
, [(Int, ShadeTree (x`WithAny`([Needle x], y)))])
-> Arr.Vector (y, Neighbourhood x)
locMesh ((i₀, locT), neighRegions) = Arr.map findNeighbours $ Arr.fromList locLeaves
where locLeaves :: [ (Int, x`WithAny`([Needle x], y)) ]
locLeaves = map (first (+i₀)) . zip [0..] $ onlyLeaves locT
vicinityLeaves :: [(Int, x)]
vicinityLeaves = Hask.foldMap
(\(i₀n, ngbR) -> map ((+i₀n) *** _topological)
. zip [0..]
$ onlyLeaves ngbR
) neighRegions
findNeighbours :: (Int, x`WithAny`([Needle x], y)) -> (y, Neighbourhood x)
findNeighbours (i, WithAny (vns,y) x)
= (y, Neighbourhood
(UArr.fromList $ fst<$>execState seek mempty)
locRieM )
where seek :: State [(Int, (Needle x, Needle' x))] ()
seek = do
Hask.forM_ ( fastNubBy (comparing fst)
$ map (second _topological) locLeaves
++ vicinityLeaves ++ aprioriNgbs )
$ \(iNgb, xNgb) ->
when (iNgb/=i) `id`do
let (Option (Just v)) = xNgb.-~.x
oldNgbs <- get
when (all (\(_,(_,nw)) -> visibleOverlap nw v) oldNgbs) `id`do
let w = w₀ ^/ (w₀<.>^v)
where w₀ = toDualWith locRieM v
put $ (iNgb, (v,w))
: [ neighbour
| neighbour@(_,(nv,_))<-oldNgbs
, visibleOverlap w nv
]
aprioriNgbs :: [(Int, x)]
aprioriNgbs = catMaybes
[ getOption $ (second $ const xN) <$>
positionIndex (pure locRieM) shd' xN
| v <- vns
, let xN = x.+~^v :: x ]
visibleOverlap :: Needle' x -> Needle x -> Bool
visibleOverlap w v = o < 1
where o = w<.>^v
locRieM :: Metric x
locRieM = case pointsCovers . map _topological
$ onlyLeaves locT
++ Hask.foldMap (onlyLeaves . snd) neighRegions of
[sh₀] -> metricf sh₀
indexWeb :: WithField ℝ Manifold x => PointsWeb x y -> WebNodeId -> Option (x,y)
indexWeb (PointsWeb rsc assocD) i
| i>=0, i<Arr.length assocD
, Right (_,x) <- indexShadeTree rsc i = pure (x, fst (assocD Arr.! i))
| otherwise = empty
unsafeIndexWebData :: PointsWeb x y -> WebNodeId -> y
unsafeIndexWebData (PointsWeb _ asd) i = fst (asd Arr.! i)
webEdges :: ∀ x y . WithField ℝ Manifold x
=> PointsWeb x y -> [((x,y), (x,y))]
webEdges web@(PointsWeb rsc assoc) = (lookId***lookId) <$> toList allEdges
where allEdges :: Set.Set (WebNodeId,WebNodeId)
allEdges = Hask.foldMap (\(i,(_, Neighbourhood ngbs _))
-> Set.fromList [(min i i', max i i')
| i'<-UArr.toList ngbs ]
) $ Arr.indexed assoc
lookId i | Option (Just xy) <- indexWeb web i = xy
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 $ (xxψ)/(xωxψ)
in yio drel )
where Option (Just yio) = geodesicBetween yψ yω
mkInterpolationSeq_lin (p₀:p₁:ps)
= mkInterpolationSeq_lin [p₀,p₁] <> mkInterpolationSeq_lin (p₁:ps)
mkInterpolationSeq_lin _ = []
sliceWeb_lin :: ∀ x y . (WithField ℝ Manifold x, Geodesic x, Geodesic y)
=> PointsWeb x y -> Cutplane x -> [(x,y)]
sliceWeb_lin web = sliceEdgs
where edgs = webEdges web
sliceEdgs cp = [ (xi d, yi d)
| ((x₀,y₀), (x₁,y₁)) <- edgs
, Option (Just d) <- [cutPosBetween cp (x₀,x₁)]
, Option (Just xi) <- [geodesicBetween x₀ x₁]
, Option (Just yi) <- [geodesicBetween y₀ y₁]
]
data GridPlanes x = GridPlanes {
_gridPlaneNormal :: Needle' x
, _gridPlaneSpacing :: Needle x
, _gridPlanesCount :: Int
}
data GridSetup x = GridSetup {
_gridStartCorner :: x
, _gridSplitDirs :: [GridPlanes x]
}
cartesianGrid2D :: (x~ℝ, y~ℝ) => ((x,x), Int) -> ((y,y), Int) -> GridSetup (x,y)
cartesianGrid2D ((x₀,x₁), nx) ((y₀,y₁), ny)
= GridSetup (x₀,y₀) [ GridPlanes (0,1) (0, (y₁y₀)/fromIntegral ny) ny
, GridPlanes (1,0) ((x₁x₀)/fromIntegral nx, 0) ny ]
splitToGridLines :: (WithField ℝ Manifold x, Geodesic x, Geodesic y)
=> PointsWeb x y -> GridSetup x -> [((x, GridPlanes x), [(x,y)])]
splitToGridLines web (GridSetup x₀ [GridPlanes dirΩ spcΩ nΩ, linePln])
= [ ((x₀', linePln), sliceWeb_lin web $ Cutplane x₀' (Stiefel1 dirΩ))
| k <- [0 .. nΩ1]
, let x₀' = x₀.+~^(fromIntegral k *^ spcΩ) ]
sampleWebAlongGrid_lin :: ∀ x y . (WithField ℝ Manifold x, Geodesic x, Geodesic y)
=> PointsWeb x y -> GridSetup x -> [(x,Option y)]
sampleWebAlongGrid_lin web grid = finalLine =<< splitToGridLines web grid
where finalLine :: ((x, GridPlanes x), [(x,y)]) -> [(x,Option y)]
finalLine ((x₀, GridPlanes _ dir nSpl), verts)
| length verts < 2 = take nSpl $ (,empty)<$>iterate (.+~^dir) x₀
finalLine ((x₀, GridPlanes _ dir nSpl), verts) = take nSpl $ go (x₀,0) intpseq
where intpseq = mkInterpolationSeq_lin
[ (metric metr $ x.-~!x₀, y) | (x,y) <- verts ]
go (x,_) [] = (,empty)<$>iterate (.+~^dir) x
go xt (InterpolationIv (_,te) f:fs)
= case break ((<te) . snd) $ iterate ((.+~^dir)***(+1)) xt of
(thisRange, xtn:_)
-> ((id***pure.f)<$>thisRange) ++ go xtn fs
Option (Just metr) = inferMetric $ webNodeRsc web
sampleWeb_2Dcartesian_lin :: (x~ℝ, y~ℝ, Geodesic z)
=> PointsWeb (x,y) z -> ((x,x),Int) -> ((y,y),Int) -> [(y,[(x,Option z)])]
sampleWeb_2Dcartesian_lin web (xspec@(_,nx)) yspec
= go . sampleWebAlongGrid_lin web $ cartesianGrid2D xspec yspec
where go [] = []
go l@(((_,y),_):_) = let (ln,l') = splitAt nx l
in (y, map (\((x,_),z) -> (x,z)) ln) : go l'
sampleEntireWeb_2Dcartesian_lin :: (x~ℝ, y~ℝ, Geodesic z)
=> PointsWeb (x,y) z -> Int -> Int -> [(y,[(x,Option z)])]
sampleEntireWeb_2Dcartesian_lin web nx ny
= sampleWeb_2Dcartesian_lin web ((x₀,x₁),nx) ((y₀,y₁),ny)
where x₀ = minimum (fst<$>pts)
x₁ = maximum (fst<$>pts)
y₀ = minimum (snd<$>pts)
y₁ = maximum (snd<$>pts)
pts = fst . fst <$> toList (localFocusWeb web)
webLocalInfo :: ∀ x y . WithField ℝ Manifold x
=> PointsWeb x y -> PointsWeb x (WebLocally x y)
webLocalInfo origWeb = result
where result = wli $ localFocusWeb origWeb
wli (PointsWeb rsc asd) = PointsWeb rsc asd'
where asd' = Arr.imap localInfo asd
localInfo i (((x,y), ngbCo), ngbH)
= ( LocalWebInfo {
_thisNodeCoord = x
, _thisNodeData = y
, _containingWeb = result
, _thisNodeId = i
, _nodeNeighbours = zip (UArr.toList $ neighbours ngbH) ngbCo
, _nodeLocalScalarProduct = localScalarProduct ngbH
, _nodeIsOnBoundary = anyUnopposed (localScalarProduct ngbH) ngbCo
}, ngbH )
anyUnopposed rieM ngbCo = (`any`ngbCo) $ \(v,_)
-> not $ (`any`ngbCo) $ \(v',_)
-> toDualWith rieM v <.>^ v' < 0
localFocusWeb :: WithField ℝ Manifold x
=> PointsWeb x y -> PointsWeb x ((x,y), [(Needle x, y)])
localFocusWeb (PointsWeb rsc asd) = PointsWeb rsc asd''
where asd' = Arr.imap (\i (y,n) -> case indexShadeTree rsc i of
Right (_,x) -> ((x,y),n) ) asd
asd''= Arr.map (\((x,y),n) ->
(((x,y), [ ( case x'.-~.x of
Option (Just v) -> v
, y')
| j<-UArr.toList (neighbours n)
, let ((x',y'),_) = asd' Arr.! j
]), n)
) asd'
nearestNeighbour :: WithField ℝ Manifold x
=> PointsWeb x y -> x -> Option (x,y)
nearestNeighbour (PointsWeb rsc asd) x = fmap lkBest $ positionIndex empty rsc x
where lkBest (iEst, (_, xEst)) = (xProx, yProx)
where (iProx, (xProx, _)) = minimumBy (comparing $ snd . snd)
$ (iEst, (xEst, metricSq locMetr vEst))
: neighbours
(yProx, _) = asd Arr.! iProx
(_, Neighbourhood neighbourIds locMetr) = asd Arr.! iEst
neighbours = [ (i, (xNgb, metricSq locMetr v))
| i <- UArr.toList neighbourIds
, let Right (_, xNgb) = indexShadeTree rsc i
Option (Just v) = xNgb.-~.x
]
Option (Just vEst) = xEst.-~.x
data WebLocally x y = LocalWebInfo {
_thisNodeCoord :: x
, _thisNodeData :: y
, _containingWeb :: PointsWeb x (WebLocally x y)
, _thisNodeId :: WebNodeId
, _nodeNeighbours :: [(WebNodeId, (Needle x, y))]
, _nodeLocalScalarProduct :: Metric x
, _nodeIsOnBoundary :: Bool
} deriving (Generic)
makeLenses ''WebLocally
instance Hask.Functor (WebLocally x) where
fmap f (LocalWebInfo co dt wb id ng sp bn)
= LocalWebInfo co (f dt) (fmap (fmap f) wb) id (map (second $ second f) ng) sp bn
instance WithField ℝ Manifold x => Comonad (WebLocally x) where
extract = _thisNodeData
duplicate lweb = unsafeIndexWebData deepened $ _thisNodeId lweb
where deepened = webLocalInfo $ _containingWeb lweb
toGraph :: WithField ℝ Manifold x => PointsWeb x y -> (Graph, Vertex -> (x, y))
toGraph wb = second (>>> \(i,_,_) -> case indexWeb wb i of {Option (Just xy) -> xy})
(graphFromEdges' edgs)
where edgs :: [(Int, Int, [Int])]
edgs = Arr.toList
. Arr.imap (\i (_, Neighbourhood ngbs _) -> (i, i, UArr.toList ngbs))
$ webNodeAssocData wb
data ConvexSet x
= EmptyConvex
| ConvexSet {
convexSetHull :: Shade' x
, convexSetIntersectors :: [Shade' x]
}
ellipsoid :: Shade' x -> ConvexSet x
ellipsoid s = ConvexSet s [s]
intersectors :: ConvexSet x -> Option (NonEmpty (Shade' x))
intersectors (ConvexSet h []) = pure (h:|[])
intersectors (ConvexSet _ (i:sts)) = pure (i:|sts)
intersectors _ = empty
instance Refinable x => Semigroup (ConvexSet x) where
a<>b = sconcat (a:|[b])
sconcat csets
| Option (Just allIntersectors) <- sconcat <$> Hask.traverse intersectors csets
, IntersectT ists <- rmTautologyIntersect perfectRefine $ IntersectT allIntersectors
, Option (Just hull') <- intersectShade's ists
= ConvexSet hull' (NE.toList ists)
| otherwise = EmptyConvex
where perfectRefine sh₁ sh₂
| sh₁`subShade'`sh₂ = pure sh₁
| sh₂`subShade'`sh₁ = pure sh₂
| otherwise = empty
itWhileJust :: (a -> Option a) -> a -> [a]
itWhileJust f x | Option (Just y) <- f x = x : itWhileJust f y
itWhileJust _ x = [x]
dupHead :: NonEmpty a -> NonEmpty a
dupHead (x:|xs) = x:|x:xs
iterateFilterDEqn_static :: (WithField ℝ Manifold x, Refinable y)
=> DifferentialEqn x y -> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)]
iterateFilterDEqn_static f = map (fmap convexSetHull)
. itWhileJust (filterDEqnSolutions_static f)
. fmap (`ConvexSet`[])
filterDEqnSolution_static :: (WithField ℝ Manifold x, Refinable y)
=> DifferentialEqn x y -> PointsWeb x (Shade' y) -> Option (PointsWeb x (Shade' y))
filterDEqnSolution_static f = localFocusWeb >>> Hask.traverse `id`
\((x,shy), ngbs) -> if null ngbs
then pure shy
else refineShade' shy
=<< intersectShade's
( propagateDEqnSolution_loc f ((x,shy), NE.fromList ngbs) )
filterDEqnSolutions_static :: (WithField ℝ Manifold x, Refinable y)
=> DifferentialEqn x y -> PointsWeb x (ConvexSet y) -> Option (PointsWeb x (ConvexSet y))
filterDEqnSolutions_static f = localFocusWeb >>> Hask.traverse `id`
\((x, shy@(ConvexSet hull _)), ngbs) -> if null ngbs
then pure shy
else ((shy<>) . ellipsoid)
<$> intersectShade's
( propagateDEqnSolution_loc f
((x,hull), second convexSetHull<$>NE.fromList ngbs) )
>>= \case EmptyConvex -> empty
c -> pure c
data SolverNodeState y = SolverNodeInfo {
_solverNodeStatus :: ConvexSet y
, _solverNodeBadness :: ℝ
, _solverNodeAge :: Int
}
makeLenses ''SolverNodeState
type OldAndNew d = (Option d, [d])
oldAndNew :: OldAndNew d -> [d]
oldAndNew (Option (Just x), l) = x : l
oldAndNew (_, l) = l
oldAndNew' :: OldAndNew d -> [(Bool, d)]
oldAndNew' (Option (Just x), l) = (True, x) : fmap (False,) l
oldAndNew' (_, l) = (False,) <$> l
filterDEqnSolutions_adaptive :: ∀ x y badness
. (WithField ℝ Manifold x, Refinable y, badness ~ ℝ)
=> MetricChoice x
-> DifferentialEqn x y
-> (x -> Shade' y -> badness)
-> PointsWeb x (SolverNodeState y)
-> Option (PointsWeb x (SolverNodeState y))
filterDEqnSolutions_adaptive mf f badness' oldState
= fmap (fromTopWebNodes mf . concat . fmap retraceBonds
. Hask.toList . webLocalInfo . webLocalInfo)
$ Hask.traverse (uncurry localChange) preproc'd
where preproc'd :: PointsWeb x ((WebLocally x (SolverNodeState y), [(Shade' y, badness)]))
preproc'd = fmap addPropagation $ webLocalInfo oldState
where addPropagation wl
| null neighbourHulls = (wl, [])
| otherwise = (wl, map (id&&&badness undefined) propFromNgbs)
where propFromNgbs = NE.toList $ propagateDEqnSolution_loc f
( (thisPos, thisShy), NE.fromList neighbourHulls )
thisPos = _thisNodeCoord wl :: x
thisShy = convexSetHull . _solverNodeStatus $ _thisNodeData wl
neighbourHulls = second (convexSetHull . _solverNodeStatus) . snd
<$> _nodeNeighbours wl
smallBadnessGradient, largeBadnessGradient :: ℝ
(smallBadnessGradient, largeBadnessGradient)
= ( badnessGradRated!!(n`div`4), badnessGradRated!!(n*3`div`4) )
where n = length badnessGradRated
badnessGradRated = sort [ ngBad / bad
| ( LocalWebInfo {
_thisNodeData
= SolverNodeInfo _ bad _
, _nodeNeighbours=ngbs }
, ngbProps) <- Hask.toList preproc'd
, (_, ngBad) <- ngbProps
, ngBad>bad ]
localChange :: WebLocally x (SolverNodeState y) -> [(Shade' y, badness)]
-> Option (OldAndNew (x, SolverNodeState y))
localChange localInfo@LocalWebInfo{
_thisNodeCoord = x
, _thisNodeData = SolverNodeInfo
shy@(ConvexSet hull _) prevBadness age
, _nodeNeighbours = ngbs
}
ngbProps
| null ngbs = return (pure (x, SolverNodeInfo shy prevBadness (age+1)), [])
| otherwise = do
let neighbourHulls = second (convexSetHull . _solverNodeStatus) . snd
<$> NE.fromList ngbs
(environAge, unfreshness)
= maximum&&&minimum $ age : (_solverNodeAge . snd . snd <$> ngbs)
case find (\(_, badnessN)
-> badnessN / prevBadness > smallBadnessGradient)
$ ngbProps of
Nothing | age < environAge
-> return (empty,empty)
_otherwise -> do
shy' <- ((shy<>) . ellipsoid)
<$> intersectShade's (fst <$> NE.fromList ngbProps)
newBadness <- case shy' of
EmptyConvex -> empty
ConvexSet hull' _ -> return $ badness x hull'
let updatedNode = SolverNodeInfo shy' newBadness (age+1)
stepStones <-
if unfreshness < 3
then return []
else fmap concat . forM (zip (snd<$>ngbs) ngbProps)
$ \( (vN, SolverNodeInfo (ConvexSet hullN _)
_ ageN)
, (_, nBadnessProp'd) ) -> do
case ageN of
_ | ageN > 0
, badnessGrad <- nBadnessProp'd / prevBadness
, badnessGrad > largeBadnessGradient -> do
let stepV = vN^/2
xStep = x .+~^ stepV
shyStep <- intersectShade's $
propagateDEqnSolution_loc f
( (xStep, hull)
, NE.cons (negateV stepV, hull)
$ fmap (\(vN',hullN')
-> (vN'^-^stepV, hullN') )
neighbourHulls )
return [( xStep
, SolverNodeInfo (ellipsoid shyStep)
(badness xStep shyStep) 1
)]
_otherwise -> return []
let updated = (x, updatedNode)
return $ (pure updated, stepStones)
totalAge = maximum $ _solverNodeAge . _thisNodeData . fst <$> preproc'd
errTgtModulation = (1) . (`mod'`1) . negate . sqrt $ fromIntegral totalAge
badness x = badness' x . (shadeNarrowness %~ (^* errTgtModulation))
retraceBonds :: WebLocally x (WebLocally x (OldAndNew (x, SolverNodeState y)))
-> [((x, [Needle x]), SolverNodeState y)]
retraceBonds locWeb@LocalWebInfo{ _thisNodeId = myId
, _thisNodeCoord = xOld
, _nodeLocalScalarProduct = locMetr }
= [ ( (x, fst<$>neighbourCandidates), snsy)
| (isOld, (x, snsy)) <- focused
, let neighbourCandidates
= [ (v,nnId)
| (_,ngb) <- knownNgbs
, (Option (Just v), nnId)
<- case oldAndNew $ ngb^.thisNodeData of
[] -> [ (xN.-~.x, nnId)
| (nnId, (_,nnWeb)) <- ngb^.nodeNeighbours
, nnId /= myId
, (xN,_) <- oldAndNew nnWeb ]
l -> [(xN.-~.x, ngb^.thisNodeId) | (xN,_) <- l]
]
possibleConflicts = [ metricSq locMetr v
| (v,nnId)<-neighbourCandidates
, nnId > myId ]
, isOld || null possibleConflicts
|| minimum possibleConflicts > oldMinDistSq / 4
]
where focused = oldAndNew' $ locWeb^.thisNodeData^.thisNodeData
knownNgbs = snd <$> locWeb^.nodeNeighbours
oldMinDistSq = minimum [ metricSq locMetr vOld
| (_,ngb) <- knownNgbs
, let Option (Just vOld) = ngb^.thisNodeCoord .-~. xOld
]
iterateFilterDEqn_adaptive :: (WithField ℝ Manifold x, Refinable y)
=> MetricChoice x
-> DifferentialEqn x y
-> (x -> Shade' y -> ℝ)
-> PointsWeb x (Shade' y) -> [PointsWeb x (Shade' y)]
iterateFilterDEqn_adaptive mf f badness
= map (fmap (convexSetHull . _solverNodeStatus))
. itWhileJust (filterDEqnSolutions_adaptive mf f badness)
. fmap (\((x,shy),_) -> SolverNodeInfo (ellipsoid shy)
(badness x shy)
1
)
. localFocusWeb