moonlight-triangulation-1.4.0.2: bench/spade-compare/hs/Moonlight/Triangulation/Bench/SpadeCompare/Removal.hs
{-# LANGUAGE DataKinds #-}
-- | Removal lanes, their gates, and the removal diagnostics.
module Moonlight.Triangulation.Bench.SpadeCompare.Removal where
import Control.Monad (when)
import Data.Foldable (traverse_)
import Data.List (sortOn)
import System.Exit (exitWith, ExitCode (ExitFailure))
import System.FilePath ((</>))
import System.IO (hPutStrLn, stderr)
import qualified Data.Vector as V
import Moonlight.Triangulation
import Moonlight.Triangulation.Dcel
import Moonlight.Triangulation.Handles.HandleDefs
import Moonlight.Triangulation.Math
import Moonlight.Triangulation.PointLocation
import Moonlight.Triangulation.Removal
import Moonlight.Triangulation.Session
import Moonlight.Triangulation.Types
import Moonlight.Triangulation.Bench.SpadeCompare.Support
import Moonlight.Triangulation.Bench.SpadeCompare.Input (randomPoints)
-- | The removal gate: a deterministic prefix of the input, removed by
-- coordinate, and the surviving canonical edge set. Delaunay triangulations
-- of point sets in general position are unique, so the survivor is forced
-- and both implementations must agree bit-for-bit.
writeRemovalGate :: FilePath -> Int -> Int -> IO ()
writeRemovalGate directory pointCount removalCount = do
let points = randomPoints 0x9e3779b97f4a7c15 pointCount
triangulation <- buildTriangulation <$> require (delaunay unitElementDefaults (V.fromList points))
surviving <- removeAllByCoordinate triangulation (V.fromList (take removalCount points))
writeFile
(directory </> ("removal-" <> show pointCount <> "-" <> show removalCount <> "-edges.txt"))
(unlines (canonicalEdges surviving))
-- | Remove vertices by coordinate inside one session: one O(n) dense thaw and
-- publication for the whole program, then O(degree) topology work per removal.
-- The snapshot 'locateAndRemove' fold pays that O(n) boundary per removal.
removeAllByCoordinate
:: DelaunayTriangulation Point
-> V.Vector Point
-> IO (DelaunayTriangulation Point)
removeAllByCoordinate triangulation removals =
case withSession triangulation 0 (removeManyAt removals) of
Left refusal -> fail (show refusal)
Right (outcomes, surviving, _) -> do
V.imapM_ checkRemoved outcomes
pure surviving
where
checkRemoved index outcome = case outcome of
Nothing ->
error
( "removal target is not a vertex: index="
<> show index
<> " point="
<> ( case removals V.! index of
Point x y -> hex64 x <> "," <> hex64 y
)
)
Just _ -> pure ()
-- | Removal by coordinate through the persistent entry point: a thaw and a
-- publication per removal, against the session fold's one of each for the whole
-- program. Same removals, same order, same survivor — only the price differs.
removeAllPersistently
:: DelaunayTriangulation Point
-> V.Vector Point
-> IO (DelaunayTriangulation Point)
removeAllPersistently =
V.foldM'
(\current point -> do
outcome <- require (locateAndRemove current point)
case outcome of
Nothing -> fail ("removal target is not a vertex: " <> pointHex point)
Just result -> pure (removalTriangulation result)
)
writePersistentRemovalGate :: FilePath -> Int -> Int -> IO ()
writePersistentRemovalGate directory pointCount removalCount = do
let points = randomPoints 0x9e3779b97f4a7c15 pointCount
base <- buildTriangulation <$> require (delaunay unitElementDefaults (V.fromList points))
surviving <- removeAllPersistently base (V.fromList (take removalCount points))
writeFile
( directory
</> ("persistent-removal-" <> show pointCount <> "-" <> show removalCount <> "-edges.txt")
)
(unlines (canonicalEdges surviving))
-- | Removal audited against the one oracle that needs no second implementation:
-- the Delaunay triangulation of a point set in general position is unique, so
-- removing a prefix must leave exactly the triangulation a fresh bulk load of
-- the survivors produces. Reports the first removal at which that identity
-- breaks, so a divergence is dated rather than merely observed.
--
-- Every removal also reports the handle swap-compaction relocated into the
-- freed slot, and that handle is the caller's only re-anchor after a removal
-- invalidates the rest. An out-of-range or misnamed handle is invisible to the
-- edge-set oracle, so it is checked here against the vertex the mesh moved.
reportRemovalAudit :: Int -> Int -> IO ()
reportRemovalAudit pointCount removalCount = do
let allPoints = V.fromList (randomPoints 0x9e3779b97f4a7c15 pointCount)
base <- buildTriangulation <$> require (delaunay unitElementDefaults allPoints)
let step current index
| index >= removalCount = putStrLn ("clean-through " <> show removalCount)
| otherwise =
case allPoints V.!? index of
Nothing -> putStrLn ("input-exhausted-at-removal " <> show index)
Just point -> do
outcome <- require (locateAndRemove current point)
case outcome of
Nothing -> putStrLn ("located-nothing-at-removal " <> show index)
Just result -> do
let nextTriangulation = removalTriangulation result
movedPoint =
vertexPoint current (VertexId (fromIntegral (numVertices current - 1)))
case removalOutcomeSwap (removalOutcome result) of
Nothing
| vertexPoint current (VertexId (fromIntegral (numVertices current - 1))) /= point ->
putStrLn ("swap-handle-missing-at-removal " <> show index)
| otherwise -> pure ()
Just (swappedIn@(VertexId rawSwappedIn), _)
| toInteger rawSwappedIn >= toInteger (numVertices nextTriangulation) ->
putStrLn
( "swap-handle-out-of-range-at-removal "
<> show index
<> " handle="
<> show rawSwappedIn
<> " vertices="
<> show (numVertices nextTriangulation)
)
| vertexPoint nextTriangulation swappedIn /= movedPoint ->
putStrLn
( "swap-handle-misnamed-at-removal "
<> show index
<> " handle="
<> show rawSwappedIn
)
| otherwise -> pure ()
fresh <-
buildTriangulation
<$> require (delaunay unitElementDefaults (V.drop (index + 1) allPoints))
let removalEdges = canonicalEdges nextTriangulation
freshEdges = canonicalEdges fresh
if removalEdges == freshEdges
then step nextTriangulation (index + 1)
else do
let (onlyRemoval, onlyFresh) = sortedDifference removalEdges freshEdges
putStrLn ("diverged-at-removal " <> show index)
putStrLn
( "removed-point "
<> case point of
Point x y -> hex64 x <> "," <> hex64 y
)
putStrLn ("removal-vertices " <> show (numVertices nextTriangulation))
putStrLn ("fresh-vertices " <> show (numVertices fresh))
putStrLn ("removal-edges " <> show (length removalEdges))
putStrLn ("fresh-edges " <> show (length freshEdges))
putStrLn ("only-in-removal " <> show (length onlyRemoval))
putStrLn ("only-in-fresh " <> show (length onlyFresh))
let (invertedRemoval, illegalRemoval) = meshHealth nextTriangulation
(invertedFresh, illegalFresh) = meshHealth fresh
putStrLn ("removal-inverted-faces " <> show invertedRemoval)
putStrLn ("removal-illegal-half-edges " <> show illegalRemoval)
putStrLn ("fresh-inverted-faces " <> show invertedFresh)
putStrLn ("fresh-illegal-half-edges " <> show illegalFresh)
traverse_ (putStrLn . (" removal-only " <>)) (take 8 onlyRemoval)
traverse_ (putStrLn . (" fresh-only " <>)) (take 8 onlyFresh)
step base 0
-- | The geometry around one removal, reconstructed with public API only. For
-- an interior removal the cavity border is the cyclic link of the removed
-- vertex, and nothing but a drain flip can destroy one of those non-star
-- edges: cleanup removes only the star, and swap-remove touches only the
-- removed vertex. The reported hull predicate guards that interpretation.
-- The fan origin 'remeshRing' would pick is also public: 'collectOutgoing'
-- starts at the stored out-edge, so the origin is the destination of that
-- edge.
reportRemovalContext :: Int -> Int -> IO ()
reportRemovalContext pointCount target = do
let allPoints = V.fromList (randomPoints 0x9e3779b97f4a7c15 pointCount)
case allPoints V.!? target of
Nothing -> do
hPutStrLn stderr "removal-context: target index out of range"
exitWith (ExitFailure 2)
Just removedPoint -> do
base <- buildTriangulation <$> require (delaunay unitElementDefaults allPoints)
before <- removeAllByCoordinate base (V.take target allPoints)
after <- removeAllByCoordinate before (V.singleton removedPoint)
fresh <- buildTriangulation <$> require (delaunay unitElementDefaults (V.drop (target + 1) allPoints))
lab <- buildTriangulation <$> require (delaunay unitElementDefaults (V.drop target allPoints))
removedQuery <- require (mkQueryPoint removedPoint)
putStrLn ("lab-clean " <> show (canonicalEdges before == canonicalEdges lab))
case locatePoint before removedQuery of
OnVertex removedVertex -> do
let ring = ringOf before removedVertex
ringPoints = map snd ring
ringMember point = any (== point) ringPoints
rOnHull = any (\(edge, _) -> incidentFace before edge == outerFace) ring
afterPairs = canonicalEdgePairs after
freshPairs = canonicalEdgePairs fresh
borderPairs = cyclicPairs ringPoints
destroyedBorder = [pair | pair <- borderPairs, not (any (== pair) afterPairs)]
removalOnly = [pair | pair <- afterPairs, not (any (== pair) freshPairs)]
freshOnly = [pair | pair <- freshPairs, not (any (== pair) afterPairs)]
putStrLn ("r-on-hull " <> show rOnHull)
putStrLn ("ring-size " <> show (length ringPoints))
traverse_
(\(index, (_, point)) -> putStrLn (" ring " <> show index <> " " <> pointHex point))
(zip [(0 :: Int) ..] ring)
putStrLn ("border-edges-before " <> show (length borderPairs))
putStrLn ("border-edges-destroyed " <> show (length destroyedBorder))
traverse_
(\(p, q) -> putStrLn (" destroyed " <> pointHex p <> " " <> pointHex q))
destroyedBorder
case ringPoints of
[] -> putStrLn "fan-unavailable empty-ring"
fanOriginPoint : ringTail ->
if rOnHull
then putStrLn "fan-unavailable hull-removal"
else do
let fanOrientations =
zipWith3
(\index p q -> (index, p, q, orient2d fanOriginPoint p q))
[(1 :: Int) ..]
ringTail
(drop 1 ringTail)
fanCandidates =
zipWith3
( \(index, priorPoint) current following ->
( index
, current
, orient2d priorPoint following fanOriginPoint == GT
&& orient2d following priorPoint current == GT
, inCircle current fanOriginPoint priorPoint following
)
)
(zip [(2 :: Int) ..] ringTail)
(drop 1 ringTail)
(drop 2 ringTail)
putStrLn ("fan-origin " <> pointHex fanOriginPoint)
putStrLn
( "fan-inverted "
<> show (length [() | (_, _, _, verdict) <- fanOrientations, verdict == LT])
<> " of "
<> show (length fanOrientations)
)
traverse_
( \(index, p, q, verdict) ->
when (verdict == LT) $
putStrLn
( " fan-inverted-at "
<> show index
<> " "
<> pointHex fanOriginPoint
<> " "
<> pointHex p
<> " "
<> pointHex q
)
)
fanOrientations
traverse_
( \(index, current, convex, circle) ->
putStrLn
( " fan-candidate "
<> show index
<> " "
<> pointHex current
<> "->"
<> pointHex fanOriginPoint
<> " convex="
<> show convex
<> " incircle="
<> show circle
)
)
fanCandidates
putStrLn ("only-in-removal " <> show (length removalOnly))
traverse_
( \(p, q) ->
putStrLn
( " removal-only "
<> pointHex p
<> " "
<> pointHex q
<> " ring="
<> show (ringMember p, ringMember q)
)
)
removalOnly
putStrLn ("only-in-fresh " <> show (length freshOnly))
traverse_
( \(p, q) ->
putStrLn
( " fresh-only "
<> pointHex p
<> " "
<> pointHex q
<> " ring="
<> show (ringMember p, ringMember q)
)
)
freshOnly
let inverted = invertedFaces after
putStrLn ("inverted-faces " <> show (length inverted))
traverse_
( \(p, q, r) ->
putStrLn
( " inverted-face "
<> pointHex p
<> " "
<> pointHex q
<> " "
<> pointHex r
<> " ring="
<> show (ringMember p, ringMember q, ringMember r)
)
)
inverted
other -> putStrLn ("locate-unexpected " <> show other)
where
ringOf triangulation vertex =
case vertexOutEdge triangulation vertex of
Nothing -> []
Just start -> walk start start []
where
walk start current acc =
let ringVertex = destination triangulation current
nextEdge = counterClockwise triangulation current
acc' = (current, vertexPoint triangulation ringVertex) : acc
in if nextEdge == start
then reverse acc'
else walk start nextEdge acc'
cyclicPairs (first : second : rest) =
zipWith canonicalPair (first : second : rest) ((second : rest) ++ [first])
cyclicPairs _ = []
canonicalEdgePairs triangulation =
sortOn
(\(p, q) -> (pointKey p, pointKey q))
[ canonicalPair p q
| edge <- undirectedEdges triangulation
, let (from, to) = undirectedEndpoints triangulation edge
p = vertexPoint triangulation from
q = vertexPoint triangulation to
]
canonicalPair p q =
if pointKey p <= pointKey q then (p, q) else (q, p)
pointKey (Point x y) = (x, y)
invertedFaces triangulation =
[ (vertexPoint triangulation a, vertexPoint triangulation b, vertexPoint triangulation c)
| face <- innerFaces triangulation
, Just (a, b, c) <- [innerFaceVertices triangulation face]
, orient2d
(vertexPoint triangulation a)
(vertexPoint triangulation b)
(vertexPoint triangulation c)
== LT
]
-- | Two independent health measures, which together decide what a divergence
-- means. Delaunay's theorem says a /valid/ triangulation with no locally
-- illegal edge is globally Delaunay — so a mesh that differs from the fresh
-- build while carrying no illegal edge cannot be a valid triangulation, and
-- one carrying illegal edges was simply left under-legalized.
meshHealth
:: Triangulation mode vertex directed undirected face
-> (Int, Int)
meshHealth triangulation =
( length [() | face <- innerFaces triangulation, not (isCounterClockwise face)]
, length
[ ()
| face <- innerFaces triangulation
, edge <- faceEdgeList face
, edgeIsIllegal edge
]
)
where
faceEdgeList face = case innerFaceDirectedEdges triangulation face of
Nothing -> []
Just (e0, e1, e2) -> [e0, e1, e2]
isCounterClockwise face = case innerFaceVertices triangulation face of
Nothing -> True
Just (a, b, c) ->
orient2d
(vertexPoint triangulation a)
(vertexPoint triangulation b)
(vertexPoint triangulation c)
== GT
edgeIsIllegal edge =
let twin = reverseEdge edge
in incidentFace triangulation twin /= outerFace
&& inCircle
(vertexPoint triangulation (origin triangulation edge))
(vertexPoint triangulation (destination triangulation edge))
(vertexPoint triangulation (destination triangulation (next triangulation edge)))
(vertexPoint triangulation (destination triangulation (next triangulation twin)))
== GT
-- | The Delaunay triangulation of the survivors, built from scratch. Diffed
-- against a removal result it decides whether that result is the Delaunay
-- triangulation it is required to be.
writeFreshSuffix :: Int -> Int -> FilePath -> IO ()
writeFreshSuffix pointCount dropped path = do
triangulation <-
buildTriangulation
<$> require
( delaunay
unitElementDefaults
(V.drop dropped (V.fromList (randomPoints 0x9e3779b97f4a7c15 pointCount)))
)
writeFile path (unlines (canonicalEdges triangulation))
-- | Symmetric difference of two ascending lists, by merge.
sortedDifference :: Ord a => [a] -> [a] -> ([a], [a])
sortedDifference [] right = ([], right)
sortedDifference left [] = (left, [])
sortedDifference left@(l : ls) right@(r : rs)
| l == r = sortedDifference ls rs
| l < r = let (onlyLeft, onlyRight) = sortedDifference ls right in (l : onlyLeft, onlyRight)
| otherwise = let (onlyLeft, onlyRight) = sortedDifference left rs in (onlyLeft, r : onlyRight)