moonlight-triangulation 1.4.0.5 → 1.5.0.0

raw patch · 84 files changed

+7861/−2818 lines, 84 filesdep ~basedep ~containersdep ~moonlight-triangulationPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: base, containers, moonlight-triangulation

API changes (from Hackage documentation)

Files

• CHANGELOG.md +28/−0
• README.md +16/−2
• bench/alpha/Main.hs +78/−3
• bench/power/Moonlight/Triangulation/PowerBench.hs +248/−2
• bench/spade-compare/hs/Main.hs +3/−7
• bench/spade-compare/hs/Moonlight/Triangulation/Bench/SpadeCompare/Gates.hs +2/−4
• bench/spade-compare/hs/Moonlight/Triangulation/Bench/SpadeCompare/Hierarchy.hs +6/−1
• bench/spade-compare/hs/Moonlight/Triangulation/Bench/SpadeCompare/Interpolation.hs +4/−3
• bench/spade-compare/hs/Moonlight/Triangulation/Bench/SpadeCompare/Intersection.hs +4/−6
• bench/spade-compare/hs/Moonlight/Triangulation/Bench/SpadeCompare/Lane.hs +1/−1
• bench/spade-compare/hs/Moonlight/Triangulation/Bench/SpadeCompare/Removal.hs +8/−8
• docs/README.md +2/−0
• docs/regular-site-algebra.md +86/−0
• docs/weighted-alpha.md +32/−0
• moonlight-triangulation.cabal +48/−11
• src-build/Moonlight/Triangulation/Insertion.hs +2/−2
• src-build/Moonlight/Triangulation/Internal/BulkLoad.hs +16/−11
• src-build/Moonlight/Triangulation/Internal/Cdt/Admission.hs +1/−1
• src-build/Moonlight/Triangulation/Internal/Cdt/Combinators.hs +0/−13
• src-build/Moonlight/Triangulation/Internal/Cdt/Corridor.hs +9/−13
• src-build/Moonlight/Triangulation/Internal/Cdt/Corridor/Trace.hs +11/−11
• src-build/Moonlight/Triangulation/Internal/Cdt/Recovery.hs +4/−5
• src-build/Moonlight/Triangulation/Internal/Cdt/Region.hs +2/−3
• src-build/Moonlight/Triangulation/Internal/Cdt/Segment.hs +9/−9
• src-build/Moonlight/Triangulation/Internal/Cdt/Split.hs +6/−8
• src-build/Moonlight/Triangulation/Internal/Cdt/Union.hs +4/−2
• src-build/Moonlight/Triangulation/Internal/CircleSweep.hs +31/−9
• src-build/Moonlight/Triangulation/Internal/Excision.hs +1/−1
• src-build/Moonlight/Triangulation/Internal/Join/Seam.hs +4/−1
• src-build/Moonlight/Triangulation/Internal/Location.hs +6/−12
• src-build/Moonlight/Triangulation/Internal/Refinement.hs +13/−11
• src-build/Moonlight/Triangulation/Internal/Session.hs +56/−17
• src-build/Moonlight/Triangulation/Refinement.hs +14/−19
• src-build/Moonlight/Triangulation/Session.hs +1/−0
• src-cell-complex/Moonlight/Triangulation/LabelledComplex.hs +131/−0
• src-core/Moonlight/Triangulation/Internal/ExactRational.hs +8/−0
• src-core/Moonlight/Triangulation/Scalar.hs +0/−64
• src-dcel/Moonlight/Triangulation/Alpha.hs +14/−11
• src-dcel/Moonlight/Triangulation/Dcel.hs +34/−34
• src-dcel/Moonlight/Triangulation/Exact.hs +8/−1
• src-dcel/Moonlight/Triangulation/Handles/Dynamic.hs +11/−11
• src-dcel/Moonlight/Triangulation/Internal/CellSet.hs +20/−17
• src-dcel/Moonlight/Triangulation/Internal/DcelOperations/Chain.hs +156/−1
• src-dcel/Moonlight/Triangulation/Internal/DcelOperations/FlipRule.hs +4/−1
• src-dcel/Moonlight/Triangulation/Internal/DcelOperations/Hull.hs +1/−1
• src-dcel/Moonlight/Triangulation/Internal/DcelOperations/Normalize.hs +0/−1
• src-dcel/Moonlight/Triangulation/Internal/FaceProbe.hs +1/−3
• src-dcel/Moonlight/Triangulation/Internal/FloodFillIterator.hs +36/−32
• src-dcel/Moonlight/Triangulation/Internal/HandleDefs.hs +19/−0
• src-dcel/Moonlight/Triangulation/Internal/Mutable.hs +12/−12
• src-dcel/Moonlight/Triangulation/Internal/OperationState.hs +2/−2
• src-dcel/Moonlight/Triangulation/Internal/Predicates.hs +143/−0
• src-dcel/Moonlight/Triangulation/Internal/Representation.hs +4/−6
• src-dcel/Moonlight/Triangulation/Internal/Validation.hs +33/−35
• src-dcel/Moonlight/Triangulation/IntersectionIterator.hs +92/−35
• src-dcel/Moonlight/Triangulation/Math.hs +93/−64
• src-dcel/Moonlight/Triangulation/PointLocation.hs +3/−3
• src-dcel/Moonlight/Triangulation/Simplex.hs +226/−0
• src-dcel/Moonlight/Triangulation/Valuation.hs +138/−5
• src-dual/Moonlight/Triangulation/Interpolation.hs +15/−18
• src-dual/Moonlight/Triangulation/Voronoi.hs +3/−20
• src-dual/Moonlight/Triangulation/Voronoi/Handles.hs +4/−19
• src-overlay/Moonlight/Triangulation/Embedding.hs +3/−3
• src-public/Moonlight/Triangulation.hs +74/−0
• src-public/Moonlight/Triangulation/HintGenerator.hs +8/−7
• src-public/Moonlight/Triangulation/Internal/PowerDiagram/Edit.hs +1376/−0
• src-public/Moonlight/Triangulation/Internal/PowerDiagram/Generator.hs +264/−0
• src-public/Moonlight/Triangulation/Internal/PowerDiagram/Hull.hs +925/−0
• src-public/Moonlight/Triangulation/Internal/PowerDiagram/Locality.hs +235/−0
• src-public/Moonlight/Triangulation/Internal/PowerDiagram/Model.hs +558/−0
• src-public/Moonlight/Triangulation/Internal/PowerDiagram/Projection.hs +601/−0
• src-public/Moonlight/Triangulation/Internal/PowerDiagram/Section.hs +489/−0
• src-public/Moonlight/Triangulation/Internal/RegularTriangulation.hs +0/−1111
• src-public/Moonlight/Triangulation/Overlay.hs +1/−0
• src-public/Moonlight/Triangulation/PowerDiagram.hs +22/−913
• src-public/Moonlight/Triangulation/RegularAlpha.hs +393/−0
• src-zigzag/Moonlight/Triangulation/Zigzag.hs +75/−129
• test/algebra/Moonlight/Triangulation/PowerDiagramSpec.hs +558/−10
• test/algebra/Moonlight/Triangulation/ValuationSpec.hs +65/−1
• test/cell-complex/Main.hs +56/−1
• test/native/Moonlight/Triangulation/ExactEmbeddingSpec.hs +1/−1
• test/native/Moonlight/Triangulation/FilteredPredicateOptimizationSpec.hs +15/−12
• test/native/Moonlight/Triangulation/NativeSpec.hs +162/−34
• test/public-components/Main.hs +13/−4

@@ -10,6 +10,34 @@ The serialization format carries its own version tag, independent of the package version; any change to it is recorded here explicitly. +## Unreleased++## 1.5.0.0 - 2026-08-30++* Make the complete stable-labelled `PowerSite` section authoritative in+  `RegularTriangulation`. Add pure insertion, removal, and batch reweighting+  with exhaustive typed disposition transitions, plus bounded power-cell+  construction that reuses an already normalized regular value.+* Collapse edit publication to changed-site support plus exact disposition+  transitions. A normalized face/edge/star section supports exact slope walks,+  insertion conflict cavities, interior removal cavities, and atomic reweight+  composition. Local publication certifies incidence, cavity boundary, and+  lifted convexity; unsupported degeneracies reconstruct once through the exact+  batch owner. On the retained GHC 9.14.1 `-O1` 512-site fixture, full-result+  insertion falls from 250.7 ms / 304.7 MB to 0.53 ms / 0.61 MB, removal from+  246.4 ms / 303.5 MB to 0.72 ms / 0.87 MB, and a general one-site reweight+  from 99.4 ms / 119.2 MB to 22.8 ms / 34.2 MB. A common shift of all 512+  weights takes 0.22 ms / 0.27 MB without changing incidence. Remove the+  transient peak-hull counter+  from the persistent topology receipt: unlike its remaining fields, it was a+  construction-history measurement rather than a property of the value.+* Make batch-reweight validation proportional to the patch rather than the+  resident site set, and descend changed sites without an association-list+  copy.+* Add signed exact weighted-alpha filtration over resident regular topology, a+  generic labelled planar-complex lowering into Homology, and focused weighted+  preparation and persistence benchmark lanes.+ ## 1.4.0.5 - 2026-08-29  * Add the GHC-9.14 `zigzag` sublibrary. Stable labelled activation points become

@@ -62,10 +62,14 @@ | `minkowskiSum` / `polygonOffset` | Regions must add / expand | Two regions / element and region | Region plus receipt / `MinkowskiError` | | `erodeBy` / `openWith` / `closeWith` | Regularized morphology is required | Structuring element, region | Region plus receipt / `MinkowskiError` | | `regularTriangulation` | Weighted sites need unbounded topology and exact dual geometry | Nonempty `PowerSite` family | Opaque regular topology, dispositions, segments/rays/lines, and receipt / `PowerDiagramError` |+| `insertRegularSite` / `removeRegularSite` | A persistent regular value gains or loses one stable-labelled site | `PowerSite` / label and `RegularTriangulation` | New regular value, changed-site support, and exact disposition transitions / `RegularEditError` |+| `reweightRegularSites` | Stable sites receive one optimizer or simulation weight update | `Map label PowerWeight`, `RegularTriangulation` | New regular value, changed-site support, and exact disposition transitions / `RegularEditError` | | `boundedPowerDiagram` | Labelled weighted sites need exact cells inside one finite convex domain | `ConvexPolygon`, nonempty `PowerSite` family | Total labelled dispositions plus rational-width receipt / `PowerDiagramError` |+| `boundedPowerDiagramFromRegular` | An admitted regular value already owns the weighted-site topology | `ConvexPolygon`, `RegularTriangulation` | Exact bounded cells without rebuilding the lifted hull / `PowerDiagramError` | | `upperEnvelope` | Exact affine forms need their labelled argmax decomposition in a polygonal window | `PolygonComponent`, `Map label AffineForm` | `PlanarLayer (Maybe label)` / `UpperEnvelopeError` | | `alphaShapeContainsFace` | One face at one radius is enough | `RadiusSquared`, mesh, `FaceId` | `Bool` | | `alphaFiltration` | Every critical radius matters | Delaunay mesh | `AlphaFiltration` / `AlphaFiltrationError` |+| `regularAlphaFiltration` | Weighted regular topology needs signed exact births | `RegularTriangulation` | `RegularAlphaFiltration` / `RegularAlphaError` | | `fromExactCellSet` | Exact cells need a generic complex view | `ExactCellSet` | `DCELComplex` | | `filteredAlphaComplex` | Alpha births must enter persistence | `AlphaFiltration` | Filtered chain complex / `DCELError` | | `activationAlphaZigzag` | Non-nested labelled depths need checked correspondence | `NonEmpty (ActivationSlice depth label)` | Admitted complexes, adjacent-union witnesses, and checked zigzag / `ActivationZigzagError` |@@ -75,23 +79,33 @@ [Interpret cell sets as incidence categories.](docs/category-observatory/README.md#cell-complex-and-category-interpretation) · [Compute persistent alpha topology.](docs/persistence-rose/README.md#exact-alpha-filtration-and-persistent-topology) · [Track topology across non-nested activation depths.](docs/activation-zigzag.md)+· [Compute weighted alpha persistence.](docs/weighted-alpha.md)  ## Exact power cells and affine envelopes  `PowerWeight` admits a finite signed additive offset through its own `PowerWeightError`; `PowerSite` admits its binary64 position once. `regularTriangulation` constructs the exact lifted upper hull once. It retains-every label's visibility disposition and publishes regular faces plus weighted-dual segments, rays, full lines, or collapsed degenerate edges.+every admitted labelled site as its semantic value; visibility, regular faces,+and weighted-dual segments, rays, full lines, or collapsed degenerate edges are+sealed derived projections. Pure insertion, removal, and batch reweighting+return exhaustive typed disposition transitions. Exact conflict-cavity and+face-star descent update only affected topology; an internal typed refusal+falls back once to the canonical batch constructor. `boundedPowerDiagram` then clips only regular-neighbour inequalities against the `ConvexPolygon`; globally hidden sites need no clip, while lower-dimensional sites use the complete HPI oracle. It retains one full-dimensional, lower-dimensional, empty, coincident-equivalent, or coincident-dominated result per label. `powerDiagramPlanarLayer` is only the derived full-dimensional view.+When a regular value is already available, `boundedPowerDiagramFromRegular`+reuses its normalized lifted hull.  `upperEnvelope` maps exact `c0 + cx*x + cy*y` forms into the same owner. Its planar result contains only two-dimensional winning regions, with `Nothing` outside the window; nonconvex windows and holes use the canonical overlay.++[Edit and reweight the complete regular-site section.](docs/regular-site-algebra.md)+· [Lower weighted topology into exact persistence.](docs/weighted-alpha.md)  ## Foreign bindings 

@@ -10,6 +10,8 @@   , timedValue   ) import Data.Map.Strict qualified as Map+import Data.List.NonEmpty qualified as NonEmpty+import Data.Foldable (traverse_) import Data.Vector qualified as Vector import Moonlight.Homology.Boundary   ( degreeCardinality@@ -35,12 +37,29 @@   ) import Moonlight.Triangulation.BulkLoad (delaunayGeometry) import Moonlight.Triangulation.CellComplex (filteredAlphaComplex)-import Moonlight.Triangulation.Types (DelaunayTriangulation)+import Moonlight.Triangulation.LabelledComplex (filteredPlanarComplex)+import Moonlight.Triangulation.PowerDiagram+  ( PowerSite+  , RegularTriangulation+  , powerSite+  , powerWeight+  , regularTriangulation+  , regularTriangulationReceipt+  )+import Moonlight.Triangulation.RegularAlpha+  ( PowerAlphaBirth+  , regularAlphaBirths+  , regularAlphaComplex+  , regularAlphaFiltration+  )+import Moonlight.Triangulation.Types (DelaunayTriangulation, Point)  type PreparedAlpha = FilteredFiniteChainComplex AlphaBirth Int  main :: IO ()-main = benchmarkAlphaFiltration 64+main = do+  benchmarkAlphaFiltration 64+  traverse_ benchmarkWeightedAlphaFiltration [64, 512, 2048]  benchmarkAlphaFiltration :: Int -> IO () benchmarkAlphaFiltration siteCount = do@@ -90,7 +109,9 @@   filtration <- requireRight (alphaFiltration triangulation)   requireRight (filteredAlphaComplex filtration) -preparationReceipt :: PreparedAlpha -> (Int, Int, [Int])+preparationReceipt+  :: FilteredFiniteChainComplex filtration coefficient+  -> (Int, Int, [Int]) preparationReceipt filtered =   let finite = filteredBaseComplex filtered       HomologicalDegree maximumDegree = maxHomologicalDegree finite@@ -126,3 +147,57 @@     <> replicate       (max 0 (maximumDegree + 1 - length bettiVector))       0++type PreparedWeightedAlpha = FilteredFiniteChainComplex PowerAlphaBirth Int++benchmarkWeightedAlphaFiltration :: Int -> IO ()+benchmarkWeightedAlphaFiltration siteCount = do+  sites <-+    traverse prepareWeightedSite+      (zip [0 ..] (randomPoints 0x517cc1b727220a95 siteCount))+  submitted <-+    maybe (fail "weighted alpha benchmark requires at least one site") pure+      (NonEmpty.nonEmpty sites)+  regular <-+    timedProjection+      (benchmarkLabel siteCount "regular-topology")+      regularTriangulationReceipt+      (fst <$> requireRight (regularTriangulation submitted))+  filtered <-+    timedProjection+      (benchmarkLabel siteCount "weighted-filtered-complex")+      preparationReceipt+      (prepareWeightedAlpha regular)+  profiles <-+    timedValue+      (benchmarkLabel siteCount "weighted-persistence-critical-profile")+      ( do+          (_, bettiTable) <-+            requireRight (mod2PersistentPairsWithCriticalBettiTable filtered)+          pure (criticalBettiVectors bettiTable)+      )+  putStrLn+    ( benchmarkLabel siteCount "weighted-profile"+        <> ": critical-births="+        <> show (length (filteredCriticalValues filtered))+        <> " profiles="+        <> show (length profiles)+        <> " prepared-cells="+        <> show (preparationReceipt filtered)+    )++prepareWeightedSite :: (Int, Point) -> IO (PowerSite Int)+prepareWeightedSite (label, point) = do+  weight <- requireRight (powerWeight (fromIntegral (label `mod` 11 - 5) / 256))+  requireRight (powerSite label point weight)++prepareWeightedAlpha+  :: RegularTriangulation Int+  -> IO PreparedWeightedAlpha+prepareWeightedAlpha regular = do+  filtration <- requireRight (regularAlphaFiltration regular)+  requireRight+    ( filteredPlanarComplex+        (regularAlphaComplex filtration)+        (regularAlphaBirths filtration)+    )

@@ -14,8 +14,12 @@ import Control.DeepSeq (force) import Control.Exception (evaluate) import Data.Foldable (traverse_)+import qualified Data.List as List import Data.List.NonEmpty (NonEmpty (..)) import qualified Data.List.NonEmpty as NonEmpty+import qualified Data.Map.Strict as Map+import qualified Data.Set as Set+import qualified Data.Vector as Vector import GHC.Conc.Sync (getAllocationCounter, setAllocationCounter) import Moonlight.Triangulation   ( BoundedPowerDiagram@@ -24,9 +28,13 @@   , PowerDiagramError   , PowerDiagramReceipt   , PowerSite+  , PowerWeight+  , RegularEditError+  , RegularEditResult   , RegularTriangulation   , RegularTriangulationReceipt   , boundedPowerDiagram+  , boundedPowerDiagramFromRegular   , convexPolygon   , exactPoint   , powerDiagramActiveBoundaries@@ -52,8 +60,22 @@   , powerDiagramRegularFaces   , powerDiagramSubmittedSiteConstraints   , powerSite+  , powerSiteLabel+  , powerSitePosition+  , powerSiteWeight   , powerWeight+  , powerWeightExact+  , powerWeightFromExact+  , insertRegularSite+  , removeRegularSite+  , reweightRegularSites   , regularTriangulation+  , regularEditChangedSites+  , regularEditTriangulation+  , regularEditTransitions+  , regularSiteDisposition+  , regularSites+  , RegularSiteDisposition (RegularSiteHidden)   ) import System.Mem (performGC) import Test.Tasty.Bench (Benchmark, bench, bgroup, defaultMain, nf)@@ -68,8 +90,22 @@   { powerBenchLabel :: !String   , powerBenchDomain :: !ConvexPolygon   , powerBenchSites :: !(NonEmpty (PowerSite Int))+  , powerBenchRegular :: !(RegularTriangulation Int)+  , powerBenchInsertionBase :: !(RegularTriangulation Int)+  , powerBenchInsertedSite :: !(PowerSite Int)+  , powerBenchReweight :: !PowerWeight+  , powerBenchCommonShift :: !(Map.Map Int PowerWeight)+  , powerBenchCoincidentLabel :: !Int+  , powerBenchCoincidentRegular :: !(RegularTriangulation Int)+  , powerBenchCoincidentReweight :: !PowerWeight+  , powerBenchHiddenEdit :: !(Maybe HiddenEditFixture)   } +data HiddenEditFixture = HiddenEditFixture+  { hiddenEditLabel :: !Int+  , hiddenEditLowerWeight :: !PowerWeight+  }+ type PowerResult =   Either (PowerDiagramError Int) (BoundedPowerDiagram Int, PowerDiagramReceipt) @@ -81,8 +117,14 @@     (PowerDiagramError Int)     (RegularTriangulation Int, RegularTriangulationReceipt) +type RegularEditBenchmarkResult =+  Either (RegularEditError Int) (RegularEditResult Int)++type RegularEditObservation =+  Either (RegularEditError Int) (Int, Int, Maybe (RegularSiteDisposition Int))+ siteCounts :: [Int]-siteCounts = [16, 64, 128, 169, 256, 512]+siteCounts = [16, 64, 128, 169, 256, 512, 2048]  families :: [PowerBenchFamily] families = [WellConditioned, NearParallel, CollinearSlopes]@@ -107,12 +149,57 @@     traverse (prepareSite family) (zip [0 ..] (fixturePoints family siteCount))   case NonEmpty.nonEmpty sites of     Nothing -> fail "power benchmark fixture requires at least one site"-    Just nonEmptySites ->+    Just nonEmptySites -> do+      regular <- requireRight (regularTriangulation nonEmptySites)+      insertionBaseSites <-+        case NonEmpty.nonEmpty (NonEmpty.init nonEmptySites) of+          Nothing -> fail "power edit benchmark requires at least two sites"+          Just baseSites -> pure baseSites+      insertionBase <- requireRight (regularTriangulation insertionBaseSites)+      replacementWeight <- requireRight (powerWeight 17)+      subordinateWeight <- requireRight (powerWeight (-1048576))+      subordinateReweight <- requireRight (powerWeight (-2097152))+      let representativeSite = NonEmpty.head nonEmptySites+          commonShift =+            Map.fromList+              [ ( powerSiteLabel site+                , powerWeightFromExact (powerWeightExact (powerSiteWeight site) + 1)+                )+              | site <- NonEmpty.toList nonEmptySites+              ]+      subordinateSite <-+        requireRight+          (powerSite siteCount (powerSitePosition representativeSite) subordinateWeight)+      subordinateInsertion <-+        requireRight (insertRegularSite subordinateSite (fst regular))+      let preparedRegular = fst regular+          hiddenEdit =+            fmap+              (\site ->+                 HiddenEditFixture+                   { hiddenEditLabel = powerSiteLabel site+                   , hiddenEditLowerWeight = subordinateWeight+                   })+              ( List.find+                  (\site ->+                     regularSiteDisposition (powerSiteLabel site) preparedRegular+                       == Just RegularSiteHidden)+                  (regularSites preparedRegular)+              )       pure         PowerBenchCase           { powerBenchLabel = familyLabel family <> "/n=" <> show siteCount           , powerBenchDomain = domain           , powerBenchSites = nonEmptySites+          , powerBenchRegular = preparedRegular+          , powerBenchInsertionBase = fst insertionBase+          , powerBenchInsertedSite = NonEmpty.last nonEmptySites+          , powerBenchReweight = replacementWeight+          , powerBenchCommonShift = commonShift+          , powerBenchCoincidentLabel = powerSiteLabel subordinateSite+          , powerBenchCoincidentRegular = regularEditTriangulation subordinateInsertion+          , powerBenchCoincidentReweight = subordinateReweight+          , powerBenchHiddenEdit = hiddenEdit           }  fixturePoints :: PowerBenchFamily -> Int -> [Point]@@ -200,8 +287,53 @@   , bench       (powerBenchLabel benchmarkCase <> "/construct-diagram-only")       (nf constructDiagramOnly benchmarkCase)+  , bench+      (powerBenchLabel benchmarkCase <> "/clip-prepared-regular")+      (nf constructFromRegular benchmarkCase)+  , bench+      (powerBenchLabel benchmarkCase <> "/insert-local-observation")+      (nf observeInsertion benchmarkCase)+  , bench+      (powerBenchLabel benchmarkCase <> "/insert-and-normalize")+      (nf insertAndNormalize benchmarkCase)+  , bench+      (powerBenchLabel benchmarkCase <> "/remove-local-observation")+      (nf observeRemoval benchmarkCase)+  , bench+      (powerBenchLabel benchmarkCase <> "/remove-and-normalize")+      (nf removeAndNormalize benchmarkCase)+  , bench+      (powerBenchLabel benchmarkCase <> "/reweight-local-observation")+      (nf observeReweight benchmarkCase)+  , bench+      (powerBenchLabel benchmarkCase <> "/reweight-and-normalize")+      (nf reweightAndNormalize benchmarkCase)+  , bench+      (powerBenchLabel benchmarkCase <> "/common-weight-shift")+      (nf shiftEveryWeight benchmarkCase)   ]+    <> topologyPreservingBenchmarks benchmarkCase +topologyPreservingBenchmarks :: PowerBenchCase -> [Benchmark]+topologyPreservingBenchmarks benchmarkCase =+  [ bench+      (powerBenchLabel benchmarkCase <> "/remove-coincident-subordinate")+      (nf removeCoincidentSubordinate benchmarkCase)+  , bench+      (powerBenchLabel benchmarkCase <> "/reweight-coincident-subordinate")+      (nf reweightCoincidentSubordinate benchmarkCase)+  ]+    <> maybe [] hiddenBenchmarks (powerBenchHiddenEdit benchmarkCase)+ where+  hiddenBenchmarks hidden =+    [ bench+        (powerBenchLabel benchmarkCase <> "/remove-hidden")+        (nf (removeHidden benchmarkCase) hidden)+    , bench+        (powerBenchLabel benchmarkCase <> "/reweight-hidden-downward")+        (nf (reweightHiddenDownward benchmarkCase) hidden)+    ]+ construct :: PowerBenchCase -> PowerResult construct benchmarkCase =   boundedPowerDiagram@@ -213,6 +345,120 @@  constructRegularTopology :: PowerBenchCase -> RegularResult constructRegularTopology = regularTriangulation . powerBenchSites++constructFromRegular :: PowerBenchCase -> PowerResult+constructFromRegular benchmarkCase =+  boundedPowerDiagramFromRegular+    (powerBenchDomain benchmarkCase)+    (powerBenchRegular benchmarkCase)++insertAndNormalize+  :: PowerBenchCase+  -> RegularEditBenchmarkResult+insertAndNormalize benchmarkCase =+  insertRegularSite+    (powerBenchInsertedSite benchmarkCase)+    (powerBenchInsertionBase benchmarkCase)++observeInsertion+  :: PowerBenchCase+  -> RegularEditObservation+observeInsertion benchmarkCase =+  observeRegularEdit+    (powerSiteLabel (powerBenchInsertedSite benchmarkCase))+    (insertAndNormalize benchmarkCase)++removeAndNormalize+  :: PowerBenchCase+  -> RegularEditBenchmarkResult+removeAndNormalize benchmarkCase =+  removeRegularSite+    (powerSiteLabel (powerBenchInsertedSite benchmarkCase))+    (powerBenchRegular benchmarkCase)++observeRemoval+  :: PowerBenchCase+  -> RegularEditObservation+observeRemoval benchmarkCase =+  observeRegularEdit+    (powerSiteLabel (powerBenchInsertedSite benchmarkCase))+    (removeAndNormalize benchmarkCase)++reweightAndNormalize+  :: PowerBenchCase+  -> RegularEditBenchmarkResult+reweightAndNormalize benchmarkCase =+  reweightRegularSites+    ( Map.singleton+        (powerSiteLabel (powerBenchInsertedSite benchmarkCase))+        (powerBenchReweight benchmarkCase)+    )+    (powerBenchRegular benchmarkCase)++observeReweight+  :: PowerBenchCase+  -> RegularEditObservation+observeReweight benchmarkCase =+  observeRegularEdit+    (powerSiteLabel (powerBenchInsertedSite benchmarkCase))+    (reweightAndNormalize benchmarkCase)++observeRegularEdit+  :: Int+  -> RegularEditBenchmarkResult+  -> RegularEditObservation+observeRegularEdit focus =+  fmap+    (\editResult ->+       ( Set.size (regularEditChangedSites editResult)+       , Vector.length (regularEditTransitions editResult)+       , regularSiteDisposition+           focus+           (regularEditTriangulation editResult)+       ))++shiftEveryWeight+  :: PowerBenchCase+  -> RegularEditBenchmarkResult+shiftEveryWeight benchmarkCase =+  reweightRegularSites+    (powerBenchCommonShift benchmarkCase)+    (powerBenchRegular benchmarkCase)++removeCoincidentSubordinate+  :: PowerBenchCase+  -> RegularEditBenchmarkResult+removeCoincidentSubordinate benchmarkCase =+  removeRegularSite+    (powerBenchCoincidentLabel benchmarkCase)+    (powerBenchCoincidentRegular benchmarkCase)++reweightCoincidentSubordinate+  :: PowerBenchCase+  -> RegularEditBenchmarkResult+reweightCoincidentSubordinate benchmarkCase =+  reweightRegularSites+    ( Map.singleton+        (powerBenchCoincidentLabel benchmarkCase)+        (powerBenchCoincidentReweight benchmarkCase)+    )+    (powerBenchCoincidentRegular benchmarkCase)++removeHidden+  :: PowerBenchCase+  -> HiddenEditFixture+  -> RegularEditBenchmarkResult+removeHidden benchmarkCase fixture =+  removeRegularSite (hiddenEditLabel fixture) (powerBenchRegular benchmarkCase)++reweightHiddenDownward+  :: PowerBenchCase+  -> HiddenEditFixture+  -> RegularEditBenchmarkResult+reweightHiddenDownward benchmarkCase fixture =+  reweightRegularSites+    (Map.singleton (hiddenEditLabel fixture) (hiddenEditLowerWeight fixture))+    (powerBenchRegular benchmarkCase)  data SharedPowerAllocation = SharedPowerAllocation   { sharedFullAllocation :: !Integer

@@ -104,13 +104,9 @@   SweepAngleCollapseLane -> do     points <- evaluate (force (V.fromList (nearCollinearPoints 0x9e3779b97f4a7c15 first)))     measure (buildTriangulation <$> require (delaunay unitElementDefaults points))-  -- Exactly collinear input, which never builds a face, so both loaders fall-  -- out of their sweep and degrade to plain incremental insertion for the-  -- whole load — and neither has a sub-linear answer for locating against a-  -- face-less mesh. This side scans every vertex and then every half-edge per-  -- insertion; spade collects and sorts all vertices per insertion. Both are-  -- quadratic in the load, so the ratio reports which quadratic costs more-  -- rather than whether one exists, and spade's carries the extra log factor.+  -- Exactly collinear input never builds a face. Moonlight recognizes that+  -- section after radial admission and glues its complete DCEL chain once;+  -- Spade still locates against the growing face-less mesh per insertion.   DegenerateLineLane -> do     points <- evaluate (force (V.fromList (exactlyCollinearPoints first)))     measure (buildTriangulation <$> require (delaunay unitElementDefaults points))

@@ -224,10 +224,8 @@   -- Spade timing rows. Their canonical egress is retained for the hard gate.   writePersistentGate directory 1000   writePersistentRemovalGate directory 10000 2500-  -- The remaining cliff lanes. Each times an entry point of this side that is-  -- expected to lose badly, which is worth nothing unless the two are first-  -- shown to compute the same thing — a quadratic path and a linear one that-  -- disagree are not a comparison.+  -- The remaining diagnostic gates. Timing is meaningful only after both+  -- sides have been shown to compute the same canonical result.   writeConstraintIncrementalGate directory 8000 800   writeConstraintSplitGate directory 1000   writeSweepAngleGate directory 2000

@@ -8,6 +8,7 @@ import Moonlight.Triangulation import Moonlight.Triangulation.BulkLoad import Moonlight.Triangulation.HintGenerator+import Moonlight.Triangulation.Math (validatePoint) import Moonlight.Triangulation.Session import Moonlight.Triangulation.Bench.SpadeCompare.Support import Moonlight.Triangulation.Bench.SpadeCompare.Input (randomPoints)@@ -27,7 +28,11 @@   -> Point   -> Session s Point () () () HierarchyHint hierarchyStep hierarchy point = do-  (vertex, disposition) <- insertVertexAt point point+  query <- either refuse pure (validatePoint Nothing point)+  (vertex, disposition) <-+    case hierarchyHint hierarchy query of+      Just (VertexHint seed) -> insertVertexAtNearVertex seed point point+      _ -> insertVertexAt point point   either refuse pure (updateHierarchyAfterInsertion hierarchy point vertex disposition)  -- | Arrival-order insertion with the hierarchy maintained alongside it.

@@ -16,7 +16,8 @@ import Moonlight.Triangulation import Moonlight.Triangulation.Dcel import Moonlight.Triangulation.Interpolation-import Moonlight.Triangulation.Math+import Moonlight.Triangulation.Internal.Predicates qualified as Predicates+import Moonlight.Triangulation.Math (isFinite) import Moonlight.Triangulation.PointLocation import Moonlight.Triangulation.Types import Moonlight.Triangulation.Bench.SpadeCompare.Support@@ -93,7 +94,7 @@               let !pa = vertexPoint triangulation a                   !pb = vertexPoint triangulation b                   !pc = vertexPoint triangulation c-                  !verdict = inCircle pa pb pc query+                  !verdict = Predicates.inCircle pa pb pc query                   !margins' = if insideMargin pa pb pc query then margins + 1 else margins                in if verdict == GT                     then go (adjacent triangulation face ++ rest) (face : seen) (tests + 1) margins'@@ -266,7 +267,7 @@       case innerFaceVertices triangulation face of         Nothing -> False         Just (a, b, c) ->-          inCircle+          Predicates.inCircle             (vertexPoint triangulation a)             (vertexPoint triangulation b)             (vertexPoint triangulation c)

@@ -15,16 +15,14 @@ -- | Query endpoints pulled inward. The generator fills [-1, 1]^2, so a point -- with both coordinates inside [-0.5, 0.5] is strictly inside the hull of a -- few thousand of them and locating it lands on a face. A chord that starts--- OUTSIDE the hull takes a different route on each side — spade walks the hull--- from the located hull edge while this side falls back to a scan over every--- element — so those chords are the "intersection-outside" cliff lane and are--- kept out of this one deliberately.+-- OUTSIDE the hull makes endpoint choice and local hull entry observable, so+-- those chords remain the separate "intersection-outside" cliff lane. interiorQueryPoints :: Word64 -> Int -> [Point] interiorQueryPoints seed count =   [Point (0.5 * x) (0.5 * y) | Point x y <- randomPoints seed count] --- | Query origins pushed well clear of the hull, so locating one is certain to--- report the outside of the convex hull rather than a face.+-- | Query origins pushed away from the dense centre so the corpus repeatedly+-- exercises outside-hull entry rather than only face-to-face corridors. exteriorQueryPoints :: Word64 -> Int -> [Point] exteriorQueryPoints seed count =   [Point (4 * x) (4 * y) | Point x y <- randomPoints seed count]

@@ -270,6 +270,7 @@   , parity VoronoiSweepLane 1000 0 "Voronoi dual sweep, 1k sites"   , parity DcelWalkLane 2000 0 "DCEL traversal, 2k points"   , parity IntersectionLane 10000 500 "Line intersection, 500 chords / 10k"+  , parity DegenerateLineLane 2000 0 "Exactly collinear bulk chain, 2k"   , cliff ConstraintIncrementalLane 8000 800 "Constraint one at a time, 800 / 8k"   , cliff ConstraintSplitLane 1000 0 "Constraint split, 1k crossings"   -- This deliberately stays at 1k while its parity twin also reaches 10k.@@ -277,7 +278,6 @@   , cliff HierarchyDuplicateLane 10000 500 "Hierarchy duplicate insert, 500 / 10k"   , cliff HierarchyRemovalLane 10000 250 "Hierarchy removal, 250 / 10k"   , cliff SweepAngleCollapseLane 2000 0 "Sweep hull index at one angle, 2k"-  , cliff DegenerateLineLane 2000 0 "Exactly collinear load, face-less locate, 2k"   , cliff IntersectionOutsideLane 2000 100 "Line intersection from outside the hull, 100 / 2k"   ]  where

@@ -13,7 +13,7 @@ import Moonlight.Triangulation import Moonlight.Triangulation.Dcel import Moonlight.Triangulation.Internal.HandleDefs-import Moonlight.Triangulation.Math+import Moonlight.Triangulation.Internal.Predicates qualified as Predicates import Moonlight.Triangulation.PointLocation import Moonlight.Triangulation.Removal import Moonlight.Triangulation.Session@@ -221,7 +221,7 @@                 else do                   let fanOrientations =                         zipWith3-                          (\index p q -> (index, p, q, orient2d fanOriginPoint p q))+                          (\index p q -> (index, p, q, Predicates.orient2d fanOriginPoint p q))                           [(1 :: Int) ..]                           ringTail                           (drop 1 ringTail)@@ -230,9 +230,9 @@                           ( \(index, priorPoint) current following ->                               ( index                               , current-                              , orient2d priorPoint following fanOriginPoint == GT-                                  && orient2d following priorPoint current == GT-                              , inCircle current fanOriginPoint priorPoint following+                              , Predicates.orient2d priorPoint following fanOriginPoint == GT+                                  && Predicates.orient2d following priorPoint current == GT+                              , Predicates.inCircle current fanOriginPoint priorPoint following                               )                           )                           (zip [(2 :: Int) ..] ringTail)@@ -356,7 +356,7 @@     [ (vertexPoint triangulation a, vertexPoint triangulation b, vertexPoint triangulation c)     | face <- innerFaces triangulation     , Just (a, b, c) <- [innerFaceVertices triangulation face]-    , orient2d+    , Predicates.orient2d         (vertexPoint triangulation a)         (vertexPoint triangulation b)         (vertexPoint triangulation c)@@ -388,7 +388,7 @@   isCounterClockwise face = case innerFaceVertices triangulation face of     Nothing -> True     Just (a, b, c) ->-      orient2d+      Predicates.orient2d         (vertexPoint triangulation a)         (vertexPoint triangulation b)         (vertexPoint triangulation c)@@ -397,7 +397,7 @@   edgeIsIllegal edge =     let twin = reverseEdge edge      in incidentFace triangulation twin /= outerFace-          && inCircle+          && Predicates.inCircle             (vertexPoint triangulation (origin triangulation edge))             (vertexPoint triangulation (destination triangulation edge))             (vertexPoint triangulation (destination triangulation (next triangulation edge)))

@@ -9,6 +9,8 @@ | Which library token should I depend on? | [Consumer libraries](../README.md#consumer-libraries) | | What algebraic laws do meshes satisfy? | [Package algebraic contract](../README.md#algebraic-contract) | | How are exact power cells and affine envelopes represented? | [Package affine-cell guide](../README.md#exact-power-cells-and-affine-envelopes) |+| How do weighted sites survive insertion, removal, and reweighting? | [Persistent regular-site algebra](./regular-site-algebra.md) |+| How do weighted regular cells become an exact alpha filtration? | [Exact weighted alpha filtration](./weighted-alpha.md) | | What are the set-operation types? | [Package type signatures](../README.md#set-operation-types) | | Which execution schedule publishes each result? | [Mesh publication guide](./mesh-publication.md) | | How do I build and validate components locally? | [Development guide](./development.md) |

@@ -0,0 +1,86 @@+# Persistent regular-site algebra++`Moonlight.Triangulation.PowerDiagram` owns the complete weighted-site value.+Use this surface when site identity must survive visibility changes or repeated+weight updates.++## Construct and edit++A `PowerSite label` combines a stable identity, an admitted planar position,+and an exact signed `PowerWeight`. Construction retains every site, including+hidden and coincident sites:++```haskell+regularTriangulation sites++editRegular newSite oldLabel weightPatch regular = do+  inserted <- insertRegularSite newSite regular+  removed <-+    removeRegularSite oldLabel (regularEditTriangulation inserted)+  reweightRegularSites+    weightPatch+    (regularEditTriangulation removed)+```++Each edit returns `RegularEditResult`:++- `regularEditTriangulation` is the new authoritative value.+- `regularEditChangedSites` is the exact support changed by the edit; it is+  empty for idempotent insertion, removal, or reweighting.+- `regularEditTransitions` carries every resulting disposition transition.++Read structural counts from the returned value with+`regularTriangulationReceipt`; the edit result does not duplicate them.++`RegularSiteTransition` distinguishes appearance, disappearance, and a change+between any two `RegularSiteDisposition` values. Thus removing or down-weighting+a visible site can report a previously hidden site becoming visible without a+caller-maintained hidden-site table.++## Identity and obstructions++The label is site identity. Re-inserting the same labelled position and weight+is idempotent. Reusing a label for different geometry returns+`RegularEditSiteConflict`; changing only weights belongs in+`reweightRegularSites`. A batch weight patch is applied once and returns+`RegularEditUnknownSites` before construction if any label is absent.++Distinct labels may occupy the same position. Their exact affine forms remain+separate submissions, and the result explicitly classifies them as visible,+coincident-equivalent, or coincident-dominated rather than silently deleting an+identity.++`emptyRegularTriangulation` is the insertion identity and the lawful result of+removing the final site. `regularSite`, `regularSites`, and `regularSiteCount`+read the canonical site section; faces, edges, neighbours, and visibility are+derived from it.++## Reuse the normalized topology++```haskell+(diagram, receipt) <- boundedPowerDiagramFromRegular domain regular+```++This clips regular-neighbour inequalities without reconstructing the lifted+hull. `boundedPowerDiagram domain sites` remains the one-shot entrance when the+editable regular value is not otherwise needed.++Edits descend through a persistent exact index. Coincident subordinates and+hidden sites update directly. A distinct insertion walks the slope+triangulation, discovers its lifted conflict cavity, glues the exact horizon,+and reattaches only hidden generators supported by removed faces. Interior+removal retriangulates the removed face star from its boundary and attached+hidden candidates. A single general reweight composes removal with insertion+atomically; a uniform shift of every weight translates the resident affine+section without changing topology. Multi-site patches use the direct+topology-preserving cases when exact, and otherwise normalize once through the+batch owner rather than accumulating local edits.++Every local publication checks manifold incidence, exact cavity boundary, and+strict lifted convexity. Degenerate, boundary-changing, or unsupported cases+decline to the exact batch constructor used by `regularTriangulation`; no+uncertified local result is published. On the retained 512-site fixture this+collapses insertion and removal from whole-hull reconstruction to+sub-millisecond face-star work.++[Compute signed weighted alpha persistence from the same value.](./weighted-alpha.md)

@@ -0,0 +1,32 @@+# Exact weighted alpha filtration++`Moonlight.Triangulation.RegularAlpha` derives a signed exact alpha filtration+from the canonical `RegularTriangulation`. It reuses resident regular faces,+edges, weighted-dual points, rays, and lines; no second weighted geometry is+constructed.++```haskell+filtration <- regularAlphaFiltration regular+filtered <-+  filteredPlanarComplex+    (regularAlphaComplex filtration)+    (regularAlphaBirths filtration)+pairs <- mod2PersistentPairs filtered+```++`PowerAlphaBirth` is an opaque exact rational. Negative births are lawful: a+positive site weight can make a simplex present before power radius zero.+Hidden and coincident-subordinate sites remain in the editable regular value but+contribute no simplex. `filteredPlanarComplex` is the only lowering boundary;+Homology owns chain and filtration validation.++Use ordinary `alphaFiltration` for an unweighted Delaunay DCEL. Use+`regularAlphaFiltration` when weights, exact dual rays, or repeated regular-site+edits are semantic inputs.++```console+cabal bench moonlight-triangulation-alpha-bench+```++[Return to the operation index.](../README.md#operations)+· [Edit the regular site section.](./regular-site-algebra.md)

@@ -1,6 +1,6 @@ cabal-version:       3.4 name:                moonlight-triangulation-version:             1.4.0.5+version:             1.5.0.0 synopsis:            Delaunay meshes and exact planar-region algebra. description:         Delaunay and constrained Delaunay triangulation as a lawful                      finite-set algebra, together with exact rational planar@@ -45,6 +45,8 @@   docs/README.md   docs/development.md   docs/activation-zigzag.md+  docs/regular-site-algebra.md+  docs/weighted-alpha.md   docs/examples/README.md   docs/mesh-publication.md   docs/category-observatory/README.md@@ -108,7 +110,7 @@ source-repository this   type:     git   location: https://github.com/PaleRoses/moonlight.git-  tag:      moonlight-triangulation-1.4.0.5+  tag:      moonlight-triangulation-1.5.0.0   subdir:   moonlight-triangulation  flag warnings-as-errors@@ -139,11 +141,13 @@     ghc-options: -optlo-enable-new-pm=0  -- ── the tower ────────────────────────────────────────────────────────────────--- The implementation tower is core <- dcel-internal <- build-internal. Public--- dcel and build components are restrictions of those owners: they reexport--- lawful modules without copying implementations or publishing trusted--- constructors. Dual branches directly from the DCEL owner; the apex owns only--- the genuine overlaps. The graph is acyclic and uses no @.hs-boot@ files.+-- The implementation tower separates the public exact scalar core from the+-- package-private paged-storage kernel. Both feed dcel-internal, followed by+-- build-internal. Public dcel and build components are restrictions of those+-- owners: they reexport lawful modules without copying implementations or+-- publishing trusted constructors. Dual branches directly from the DCEL owner;+-- the apex owns only the genuine overlaps. The graph is acyclic and uses no+-- @.hs-boot@ files.  library core   import: shared-properties@@ -154,6 +158,19 @@     Moonlight.Triangulation.LineSideInfo     Moonlight.Triangulation.Internal.Dyadic     Moonlight.Triangulation.Internal.ExactRational+  build-depends:+    base >= 4.19 && < 5+    , deepseq >= 1.5 && < 1.6+  ghc-options: -fexpose-all-unfoldings++-- Trusted paged storage is shared by the immutable DCEL and construction+-- owners. It is private because raw page indices and unsafe vector operations+-- are valid only after those owners have admitted array bounds and page order.+library storage-internal+  import: shared-properties+  visibility: private+  hs-source-dirs: src-core+  exposed-modules:     Moonlight.Triangulation.Internal.PageDirectory     Moonlight.Triangulation.Internal.Paged     Moonlight.Triangulation.Internal.BoxedPaged@@ -176,6 +193,7 @@     Moonlight.Triangulation.Types     Moonlight.Triangulation.Math     Moonlight.Triangulation.Exact+    Moonlight.Triangulation.Simplex     Moonlight.Triangulation.CellSet     Moonlight.Triangulation.Region     Moonlight.Triangulation.Valuation@@ -204,6 +222,7 @@     Moonlight.Triangulation.Internal.ExactSegmentEvents     Moonlight.Triangulation.Internal.Region.Publication     Moonlight.Triangulation.Internal.Region.Types+    Moonlight.Triangulation.Internal.Predicates     Moonlight.Triangulation.Internal.SegmentRelation     Moonlight.Triangulation.Internal.Types     Moonlight.Triangulation.Internal.Representation@@ -232,6 +251,7 @@     , vector >= 0.13 && < 0.14     , vector-algorithms >= 0.9 && < 0.10     , moonlight-triangulation:core+    , moonlight-triangulation:storage-internal   ghc-options: -fexpose-all-unfoldings  -- Public immutable-DCEL and exact-geometry surface. Every module is the module@@ -245,6 +265,7 @@     , Moonlight.Triangulation.Types     , Moonlight.Triangulation.Math     , Moonlight.Triangulation.Exact+    , Moonlight.Triangulation.Simplex     , Moonlight.Triangulation.CellSet     , Moonlight.Triangulation.Region     , Moonlight.Triangulation.Valuation@@ -279,6 +300,7 @@     TypeFamilies   exposed-modules:     Moonlight.Triangulation.CellComplex+    Moonlight.Triangulation.LabelledComplex   if impl(ghc < 9.14)     buildable: False   build-depends:@@ -359,6 +381,7 @@     , vector-algorithms >= 0.9 && < 0.10     , moonlight-triangulation:core     , moonlight-triangulation:dcel-internal+    , moonlight-triangulation:storage-internal   ghc-options: -fexpose-all-unfoldings  -- Public construction surface over the private kernel. Trusted admitted-input@@ -411,6 +434,7 @@     , vector >= 0.13 && < 0.14     , moonlight-triangulation:core     , moonlight-triangulation:dcel-internal+    , moonlight-triangulation:storage-internal   -- The @Binary@ instances for the identifier and point types are orphans by   -- construction: the types belong to @dcel@ and the class to @binary@, and the   -- whole purpose of this component is that neither one has to know about the@@ -456,6 +480,7 @@     , moonlight-triangulation:core     , moonlight-triangulation:dcel-internal     , moonlight-triangulation:build-internal+    , moonlight-triangulation:storage-internal  library   import: shared-properties@@ -467,6 +492,7 @@     Moonlight.Triangulation.Minkowski     Moonlight.Triangulation.Overlay     Moonlight.Triangulation.PowerDiagram+    Moonlight.Triangulation.RegularAlpha   -- The default library is the broad surface. Scoped/owning handles are the safe   -- path; Dcel is the explicitly unchecked fixed-index kernel. These are the   -- exact focused-component modules, not facade copies.@@ -506,7 +532,13 @@   other-modules:     Moonlight.Triangulation.Internal.Minkowski.Convex     Moonlight.Triangulation.Internal.Minkowski.Types-    Moonlight.Triangulation.Internal.RegularTriangulation+    Moonlight.Triangulation.Internal.PowerDiagram.Edit+    Moonlight.Triangulation.Internal.PowerDiagram.Generator+    Moonlight.Triangulation.Internal.PowerDiagram.Hull+    Moonlight.Triangulation.Internal.PowerDiagram.Locality+    Moonlight.Triangulation.Internal.PowerDiagram.Model+    Moonlight.Triangulation.Internal.PowerDiagram.Projection+    Moonlight.Triangulation.Internal.PowerDiagram.Section   build-depends:     base >= 4.19 && < 5     , containers >= 0.8 && < 0.9@@ -517,6 +549,7 @@     , moonlight-triangulation:build-internal     , moonlight-triangulation:dual     , moonlight-triangulation:overlay-internal+    , moonlight-triangulation:storage-internal   ghc-options: -fexpose-all-unfoldings  library ffi@@ -553,7 +586,7 @@     moonlight_triangulation.h   build-depends:     base >= 4.19 && < 5-    , moonlight-triangulation:ffi >= 1.4 && < 1.5+    , moonlight-triangulation:ffi >= 1.5 && < 1.6   ghc-options: -threaded   if os(windows)     options: standalone@@ -571,7 +604,7 @@     base >= 4.19 && < 5     , directory >= 1.3 && < 1.4     , filepath >= 1.4 && < 1.6-    , moonlight-triangulation:ffi >= 1.4 && < 1.5+    , moonlight-triangulation:ffi >= 1.5 && < 1.6  -- ── test slices ────────────────────────────────────────────────────────────── @@ -616,6 +649,7 @@     , moonlight-triangulation:build-internal     , moonlight-triangulation:dual     , moonlight-triangulation:overlay-internal+    , moonlight-triangulation:storage-internal  common triangulation-serialization-test-slice   other-modules: Moonlight.Triangulation.SerializationSpec@@ -963,6 +997,7 @@     , moonlight-triangulation:core     , moonlight-triangulation:dcel-internal     , moonlight-triangulation:build-internal+    , moonlight-triangulation:storage-internal  common triangulation-region-benchmark-slice   other-modules: Moonlight.Triangulation.RegionBench@@ -1180,6 +1215,7 @@     , containers >= 0.6 && < 0.9     , deepseq >= 1.5 && < 1.6     , moonlight-homology >= 0.1.0.2 && < 0.2+    , moonlight-triangulation     , moonlight-triangulation:build-internal     , moonlight-triangulation:cell-complex     , moonlight-triangulation:dcel-internal@@ -1200,7 +1236,8 @@   other-modules:     Moonlight.Triangulation.PowerBench   build-depends:-    tasty-bench >= 0.3 && < 0.6+    containers >= 0.6 && < 0.9+    , tasty-bench >= 0.3 && < 0.6  -- Exact affine-envelope construction, balanced n-ary refinement, and direct -- cell analytics against the public compositions callers would otherwise

@@ -1,5 +1,4 @@ {-# LANGUAGE AllowAmbiguousTypes #-}-{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-}@@ -35,7 +34,8 @@   , addCounter   ) import Moonlight.Triangulation.Internal.Probe (KnownProbe, Probe (..))-import Moonlight.Triangulation.Math+import Moonlight.Triangulation.Internal.Predicates+import Moonlight.Triangulation.Math (validatePoint) import Moonlight.Triangulation.Types  -- | Admit a raw site, then delegate to the shared admitted-site interpreter.

@@ -35,7 +35,12 @@ import qualified Data.Vector.Unboxed.Mutable as MUV import Data.Word (Word32) import Moonlight.Triangulation.Dcel (numInnerFaces, numVertices)-import Moonlight.Triangulation.Internal.HandleDefs (DirectedEdgeId (..), FaceId (..), VertexId (..))+import Moonlight.Triangulation.Internal.HandleDefs+  ( VertexId (..)+  , directedEdgeIdIndex+  , faceIdIndex+  , vertexIdIndex+  ) import Moonlight.Triangulation.Internal.BoxedPaged (boxedFromVector, boxedUpdate, emptyBoxedPaged) import Moonlight.Triangulation.Internal.Capacity (ensureCapacity) import Moonlight.Triangulation.Insertion (insertExistingVertexAtLocation, insertVertexAtPoint)@@ -66,7 +71,7 @@ import Moonlight.Triangulation.PointLocation (locatePointWithHint) import Moonlight.Triangulation.Internal.Probe (Probe (..)) import Moonlight.Triangulation.Internal.Representation (Triangulation (..))-import Moonlight.Triangulation.Internal.PackedIndex (noIndex)+import Moonlight.Triangulation.Internal.PackedIndex (noIndex, unpackIndex) import Moonlight.Triangulation.Internal.Paged (TransactionShape (DenseTransaction, LocalTransaction), emptyPaged, fromVector) import Moonlight.Triangulation.Internal.Transaction (runTransaction) import Moonlight.Triangulation.Types@@ -279,8 +284,8 @@   sumGeometryVertices dense vertices =     MUV.foldM'       (\(!sumX, !sumY) raw -> do-         x <- denseReadPointX dense (fromIntegral raw)-         y <- denseReadPointY dense (fromIntegral raw)+         x <- denseReadPointX dense (unpackIndex raw)+         y <- denseReadPointY dense (unpackIndex raw)          pure (sumX + x, sumY + y)       )       (0, 0)@@ -701,13 +706,13 @@                        let thawedLocation =                              case located of                                EmptyTriangulation -> MutableEmpty-                               OnEdge (DirectedEdgeId raw) -> MutableOnEdge (fromIntegral raw)-                               InFace (FaceId raw) -> MutableInFace (fromIntegral raw)+                               OnEdge edge -> MutableOnEdge (directedEdgeIdIndex edge)+                               InFace face -> MutableInFace (faceIdIndex face)                                -- The frozen locator emits no edge only for a                                -- singleton mesh. Its mutable interpreter ignores                                -- this sentinel while constructing the second vertex.                                OutsideConvexHull Nothing -> MutableOutsideHull 0-                               OutsideConvexHull (Just (DirectedEdgeId raw)) -> MutableOutsideHull (fromIntegral raw)+                               OutsideConvexHull (Just edge) -> MutableOutsideHull (directedEdgeIdIndex edge)                        ((vertex, Inserted) <$) <$> insertExistingVertexAtLocation @'ProbeOff mutable operation vertex thawedLocation              case inserted of                Left failure -> pure (Left failure)@@ -746,12 +751,12 @@   -> LocationStats   -> vertex   -> InsertionResult mode vertex directed undirected face-replaceResidentPayload triangulation resident@(VertexId raw) walked vertexData =+replaceResidentPayload triangulation resident walked vertexData =   InsertionResult     { insertionTriangulation =         triangulation           { triVertexData =-              boxedUpdate (fromIntegral raw) vertexData (triVertexData triangulation)+              boxedUpdate (vertexIdIndex resident) vertexData (triVertexData triangulation)           }     , insertionVertex = resident     , insertionDisposition = AlreadyPresent@@ -884,8 +889,8 @@   arena <- MUV.new count   forM_ [0 .. count - 1] $ \index -> do     raw <- lookupId index-    x <- readPointX mutable (fromIntegral raw)-    y <- readPointY mutable (fromIntegral raw)+    x <- readPointX mutable (unpackIndex raw)+    y <- readPointY mutable (unpackIndex raw)     let !wideX = x         !wideY = y         !deltaX = centerX - wideX

@@ -12,9 +12,9 @@ import qualified Data.IntSet as IntSet import qualified Moonlight.Triangulation.Dcel as Dcel import Moonlight.Triangulation.Internal.HandleDefs+import Moonlight.Triangulation.Internal.Predicates (segmentsProperlyCross) import Moonlight.Triangulation.Internal.Representation import Moonlight.Triangulation.Internal.Types-import Moonlight.Triangulation.Math  -- | Immutable admission descends through the sparse constrained-edge section, -- not through every ordinary edge in the requested corridor. Constraint

@@ -1,4 +1,3 @@-{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DataKinds #-}  -- | Bounded monadic folding, refusable transaction sequencing, and the handle@@ -7,14 +6,8 @@   ( foldWhileM   , bindMutable   , asConstraintStep-  , vertexInt-  , directedInt   ) where -import Moonlight.Triangulation.Internal.HandleDefs-  ( DirectedEdgeId (..)-  , VertexId (..)-  ) import Moonlight.Triangulation.Internal.Cdt.Types (CdtError (..)) import Moonlight.Triangulation.Internal.Types (BuildError) @@ -58,9 +51,3 @@ asConstraintStep :: Functor f => f (Either BuildError a) -> f (Either (CdtError) a) asConstraintStep = fmap (either (Left . CdtBuildError) Right) {-# INLINE asConstraintStep #-}--vertexInt :: VertexId -> Int-vertexInt (VertexId value) = fromIntegral value--directedInt :: DirectedEdgeId -> Int-directedInt (DirectedEdgeId value) = fromIntegral value

@@ -22,11 +22,7 @@ import Data.Bits (xor) import Moonlight.Triangulation.Internal.HandleDefs import Moonlight.Triangulation.IntersectionIterator (Intersection (..))-import Moonlight.Triangulation.Internal.Cdt.Combinators-  ( directedInt-  , foldWhileM-  , vertexInt-  )+import Moonlight.Triangulation.Internal.Cdt.Combinators (foldWhileM) import Moonlight.Triangulation.Internal.Cdt.Corridor.Trace (nextMutableIntersection) import Moonlight.Triangulation.Internal.Cdt.Types import Moonlight.Triangulation.Internal.Growable@@ -78,8 +74,8 @@   -> ST s (Either (CdtError) MutableConstraintScan) scanMutableConstraint workspace@ConstraintWorkspace{constraintProgramWords, constraintWalkBudget} mutable from to = do   clearGrowable constraintProgramWords-  lineFrom <- pointAt mutable (vertexInt from)-  lineTo <- pointAt mutable (vertexInt to)+  lineFrom <- pointAt mutable (vertexIdIndex from)+  lineTo <- pointAt mutable (vertexIdIndex to)   initialCursor <- beginMutableRecover workspace 0 from 0 False   walked <-     foldWhileM@@ -112,11 +108,11 @@ advanceMutablePlan workspace mutable lineFrom lineTo target (MutablePlanActive event cursor) _ =   case event of     EdgeIntersection directed -> do-      constrained <- readConstraint mutable (directedInt directed)+      constrained <- readConstraint mutable (directedEdgeIdIndex directed)       if constrained         then pure (MutablePlanBlocked directed)         else do-          writeConstraintProgram workspace (mutableCursorWriteAt cursor) (directedInt directed)+          writeConstraintProgram workspace (mutableCursorWriteAt cursor) (directedEdgeIdIndex directed)           continueMutablePlan             workspace             mutable@@ -159,10 +155,10 @@               False           continueMutablePlan workspace mutable lineFrom lineTo target event nextCursor     EdgeOverlap rawDirected -> do-      let rawEdge = directedInt rawDirected+      let rawEdge = directedEdgeIdIndex rawDirected       rawOrigin <- readOrigin mutable rawEdge       rawDestination <- readOrigin mutable (rawEdge `xor` 1)-      let current = vertexInt (mutableCursorAt cursor)+      let current = vertexIdIndex (mutableCursorAt cursor)           oriented             | rawOrigin == current = Just rawEdge             | rawDestination == current = Just (rawEdge `xor` 1)@@ -227,7 +223,7 @@   -> ST s MutableProgramCursor beginMutableRecover workspace header from pieceCount afterOverlap = do   writeConstraintProgram workspace header recoverProgramTag-  writeConstraintProgram workspace (header + 1) (vertexInt from)+  writeConstraintProgram workspace (header + 1) (vertexIdIndex from)   writeConstraintProgram workspace (header + 2) 0   writeConstraintProgram workspace (header + 3) 0   pure@@ -246,7 +242,7 @@   -> MutableProgramCursor   -> ST s MutableProgramCursor finishMutableRecover workspace to cursor = do-  writeConstraintProgram workspace (mutableCursorHeader cursor + 2) (vertexInt to)+  writeConstraintProgram workspace (mutableCursorHeader cursor + 2) (vertexIdIndex to)   writeConstraintProgram workspace (mutableCursorHeader cursor + 3) (mutableCursorConflictCount cursor)   pure     cursor

@@ -19,15 +19,15 @@ import Data.Bits (xor) import Moonlight.Triangulation.Internal.HandleDefs import Moonlight.Triangulation.IntersectionIterator-import Moonlight.Triangulation.Internal.Cdt.Combinators-  ( directedInt-  , foldWhileM-  , vertexInt-  )+import Moonlight.Triangulation.Internal.Cdt.Combinators (foldWhileM) import Moonlight.Triangulation.Internal.Mutable+import Moonlight.Triangulation.Internal.Predicates+  ( onClosedSegment+  , orient2d+  , orient2dCoordinates+  ) import Moonlight.Triangulation.Internal.Types-import Moonlight.Triangulation.Math-import Moonlight.Triangulation.Scalar (orient2dCoordinates)+import Moonlight.Triangulation.Math (projectionFactor)  data MutableVertexOut   = MutableVertexOutHull@@ -53,7 +53,7 @@ nextMutableIntersection mutable lineFrom lineTo current =   case current of     EdgeIntersection directed -> do-      edgeOut <- traceMutableDirectionOutOfEdge mutable (directedInt directed) lineFrom lineTo+      edgeOut <- traceMutableDirectionOutOfEdge mutable (directedEdgeIdIndex directed) lineFrom lineTo       pure $ case edgeOut of         MutableEdgeOutHull -> Nothing         MutableEdgeOutVertex vertex ->@@ -62,11 +62,11 @@           Just (EdgeIntersection (DirectedEdgeId (fromIntegral edge)))         MutableEdgeOutNone -> Nothing     VertexIntersection vertex -> do-      currentPoint <- pointAt mutable (vertexInt vertex)+      currentPoint <- pointAt mutable (vertexIdIndex vertex)       if currentPoint == lineTo         then pure Nothing         else do-          vertexOut <- traceMutableDirectionOutOfVertex mutable (vertexInt vertex) lineTo+          vertexOut <- traceMutableDirectionOutOfVertex mutable (vertexIdIndex vertex) lineTo           case vertexOut of             MutableVertexOutHull -> pure Nothing             MutableVertexOutOverlap edge ->@@ -82,7 +82,7 @@     EdgeOverlap directed       | lineFrom == lineTo -> pure Nothing       | otherwise -> do-          destination <- readOrigin mutable (directedInt directed `xor` 1)+          destination <- readOrigin mutable (directedEdgeIdIndex directed `xor` 1)           destinationPoint <- pointAt mutable destination           pure             ( if onClosedSegment lineFrom lineTo destinationPoint

@@ -21,7 +21,6 @@ import Moonlight.Triangulation.Internal.HandleDefs import Moonlight.Triangulation.Internal.Cdt.Combinators   ( foldWhileM-  , vertexInt   ) import Moonlight.Triangulation.Internal.Cdt.Corridor   ( constraintWorkspaceFor@@ -45,7 +44,7 @@   , writeScratch   ) import Moonlight.Triangulation.Internal.Types-import Moonlight.Triangulation.Scalar (orient2dCoordinates)+import Moonlight.Triangulation.Internal.Predicates (orient2dCoordinates)  data ConflictRecovery = ConflictRecovery   { recoveredConstraintEdge :: {-# UNPACK #-} !Int@@ -238,8 +237,8 @@   -> Int   -> ST s (Either (CdtError) ConflictRecovery) resolveConflictStrip workspace mutable operation from to conflictStart stripLength = do-  fromPoint <- pointAt mutable (vertexInt from)-  toPoint <- pointAt mutable (vertexInt to)+  fromPoint <- pointAt mutable (vertexIdIndex from)+  toPoint <- pointAt mutable (vertexIdIndex to)   traverse_     (\index -> do       edge <- readConstraintProgram workspace (conflictStart + index)@@ -290,7 +289,7 @@               )           )     | pendingCount <= 0 = do-        direct <- findMutableEdge mutable (vertexInt from) (vertexInt to)+        direct <- findMutableEdge mutable (vertexIdIndex from) (vertexIdIndex to)         case direct of           Just edge -> pure (Right (edge, flipped))           Nothing ->

@@ -34,8 +34,7 @@ boundedRegionFaces triangulation =   [ face   | face <- innerFaces triangulation-  , let FaceId raw = face-  , not (IntSet.member (fromIntegral raw) outside)+  , not (IntSet.member (faceIdIndex face) outside)   ]  where-  outside = IntSet.fromList [fromIntegral raw | FaceId raw <- outerRegionFaces triangulation]+  outside = IntSet.fromList (fmap faceIdIndex (outerRegionFaces triangulation))

@@ -106,8 +106,8 @@     -> ST s (Either BuildError VertexId)   placeOrReuseConstraintEndpoint mutable operation existing point payload =     case existing of-      Just vertex@(VertexId rawVertex) -> do-        writeVertexData mutable (fromIntegral rawVertex) payload+      Just vertex -> do+        writeVertexData mutable (vertexIdIndex vertex) payload         pure (Right vertex)       Nothing -> placeConstraintEndpoint mutable operation Nothing point payload @@ -122,13 +122,13 @@   -> VertexId   -> vertex   -> Either (CdtError) (ConstraintResult vertex directed undirected face)-applyConstraintToExistingEndpoints triangulation from@(VertexId rawFrom) fromPayload to@(VertexId rawTo) toPayload =+applyConstraintToExistingEndpoints triangulation from fromPayload to toPayload =   do     (request, frozen) <-       runUnmeasuredTransaction CdtBuildError LocalTransaction triangulation 0 $ \mutable operation -> do         programWords <- newGrowableWord32 256-        writeVertexData mutable (fromIntegral rawFrom) fromPayload-        writeVertexData mutable (fromIntegral rawTo) toPayload+        writeVertexData mutable (vertexIdIndex from) fromPayload+        writeVertexData mutable (vertexIdIndex to) toPayload         recoverConstraintRequest programWords mutable operation from to     pure (publishConstraintResult frozen request) @@ -207,8 +207,8 @@   :: Triangulation 'Constrained vertex directed undirected face   -> UndirectedEdgeId   -> Either (CdtError) (Triangulation 'Constrained vertex directed undirected face)-removeConstraintEdge triangulation edge@(UndirectedEdgeId raw)-  | fromIntegral raw >= edgeCount =+removeConstraintEdge triangulation edge+  | undirectedEdgeIdIndex edge >= edgeCount =       Left (ConstraintEdgeIndexOutOfRange edge edgeCount)   -- An edge carrying no constraint has nothing to retire, and answering that   -- without thawing is the difference between O(1) and a republished mesh.@@ -228,9 +228,9 @@   -> OperationState s   -> UndirectedEdgeId   -> ST s (Either (CdtError) ())-retireConstraintEdge mutable operation edge@(UndirectedEdgeId raw) = do+retireConstraintEdge mutable operation edge = do   halfEdges <- directedEdgeCount mutable-  let !directed = fromIntegral raw * 2+  let !directed = directedEdgeIdIndex (normalizedDirected edge)   if directed >= halfEdges     then pure (Left (ConstraintEdgeIndexOutOfRange edge (halfEdges `quot` 2)))     else do

@@ -18,8 +18,6 @@ import Moonlight.Triangulation.Internal.Cdt.Combinators   ( asConstraintStep   , bindMutable-  , directedInt-  , vertexInt   ) import Moonlight.Triangulation.Internal.Cdt.Corridor   ( constraintWorkspaceFor@@ -336,10 +334,10 @@       else divideWithin cursor crossing    divideWithin cursor crossing = do-    let !crossingEdge = directedInt crossing+    let !crossingEdge = directedEdgeIdIndex crossing         !oldConstraint = asUndirected crossing-    segmentFrom <- pointAt mutable (vertexInt (splitCursorAt cursor))-    segmentTo <- pointAt mutable (vertexInt target)+    segmentFrom <- pointAt mutable (vertexIdIndex (splitCursorAt cursor))+    segmentTo <- pointAt mutable (vertexIdIndex target)     oldFrom <- readOrigin mutable crossingEdge     oldTo <- readOrigin mutable (crossingEdge `xor` 1)     edgeFrom <- pointAt mutable oldFrom@@ -405,9 +403,9 @@   -> (VertexId, VertexId)   -> VertexId   -> ST s Bool-mutableEdgeWasSplit mutable (UndirectedEdgeId raw) (from, to) splitVertex = do+mutableEdgeWasSplit mutable edge (from, to) splitVertex = do   halfEdges <- directedEdgeCount mutable-  let !directed = fromIntegral raw * 2+  let !directed = directedEdgeIdIndex (normalizedDirected edge)   stillConstrained <-     if directed >= halfEdges       then pure False@@ -426,7 +424,7 @@   -> VertexId   -> ST s Bool mutableConstraintBetween mutable from to = do-  found <- findMutableEdge mutable (vertexInt from) (vertexInt to)+  found <- findMutableEdge mutable (vertexIdIndex from) (vertexIdIndex to)   case found of     Nothing -> pure False     Just edge -> readConstraint mutable edge

@@ -1,4 +1,3 @@-{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE DerivingStrategies #-} {-# LANGUAGE ScopedTypeVariables #-}@@ -63,12 +62,15 @@   ) import Moonlight.Triangulation.Internal.BoxedPaged (boxedFill) import Moonlight.Triangulation.Internal.Paged (TransactionShape (DenseTransaction, LocalTransaction))+import Moonlight.Triangulation.Internal.Predicates+  ( onClosedSegment+  , segmentsProperlyCross+  ) import Moonlight.Triangulation.Internal.Representation import Moonlight.Triangulation.Internal.Mutable (MutableDcel) import Moonlight.Triangulation.Internal.OperationState (OperationState) import Moonlight.Triangulation.Internal.Transaction (runTransaction) import Moonlight.Triangulation.Internal.Types-import Moonlight.Triangulation.Math  canonicalSegment :: Point -> Point -> CanonicalSegment canonicalSegment from to

@@ -36,6 +36,10 @@ import Moonlight.Triangulation.Internal.DcelOperations.CandidateArena   ( seedGenericPairInArena   )+import Moonlight.Triangulation.Internal.DcelOperations.Chain+  ( connectAscendingLineArena+  , prepareAscendingLineArena+  ) import Moonlight.Triangulation.Internal.DcelOperations.Normalize   ( LegalizationDrain (..)   , drainDenseUnconstrainedGenericLegalization@@ -43,6 +47,7 @@ import Moonlight.Triangulation.Internal.Mutable   ( DenseMutableDcel   , MutableDcel+  , connectedCount   , denseMutableDcel   , denseMutableOwner   , denseFaceEdges@@ -68,8 +73,9 @@   , storeLegalizationArena   , writeScratch   )+import Moonlight.Triangulation.Internal.PackedIndex (noIndex, unpackIndex) import Moonlight.Triangulation.Internal.Probe (Probe (..))-import Moonlight.Triangulation.Scalar (orient2dCoordinates)+import Moonlight.Triangulation.Internal.Predicates (orient2dCoordinates) import Moonlight.Triangulation.Types (BuildError (..))  -- | The hull is an angular index over the DCEL outer-face cycle. The cycle@@ -258,9 +264,6 @@       (MUV.unsafeSlice start (end - start) vertices) {-# INLINE sortGeometryPrefixFibre #-} -noOuterEdge :: Word32-noOuterEdge = maxBound- -- | Circle sweep over one mutable DCEL, consuming one radial vertex order. -- Coordinates are read from the authoritative dense DCEL. Insertions -- initially close only acute hull turns. One terminal Graham pass restores@@ -286,6 +289,25 @@ circleSweepInsertDense dense operation arena   | MUV.length arena == 0 = pure (Right 0)   | otherwise = do+      residentCount <- connectedCount mutable+      ascendingLine <-+        if residentCount == 0+          then prepareAscendingLineArena dense arena+          else pure Nothing+      case ascendingLine of+        Just line -> connectAscendingLineArena dense operation line+        Nothing -> circleSweepInsertArea dense operation arena+ where+  !mutable = denseMutableOwner dense+{-# INLINE circleSweepInsertDense #-}++circleSweepInsertArea+  :: forall s vertex directed undirected face+   . DenseMutableDcel s vertex directed undirected face+  -> OperationState s+  -> MUV.MVector s Word32+  -> ST s (Either BuildError Int)+circleSweepInsertArea dense operation arena = do       let !ordered = MUV.length arena       seed <- insertSeed 0       case seed of@@ -363,7 +385,7 @@     | index >= MUV.length arena = pure (Right index)     | otherwise = do         raw <- MUV.unsafeRead arena index-        result <- insertExistingVertex @'ProbeOff mutable operation (fromIntegral raw)+        result <- insertExistingVertex @'ProbeOff mutable operation (unpackIndex raw)         case result of           Left failure -> pure (Left failure)           Right () -> do@@ -389,7 +411,7 @@     | index >= ordered = pure (Right (skippedCount, fastCount, flips, maxDepth, cursor, candidateArena))     | otherwise = do         raw <- MUV.unsafeRead arena index-        let !vertex = fromIntegral raw+        let !vertex = unpackIndex raw         queryXWide <- denseReadPointX dense vertex         queryYWide <- denseReadPointY dense vertex         let !queryAngle =@@ -554,7 +576,7 @@ rebuildBuckets :: DenseMutableDcel s vertex directed undirected face -> Hull s -> Int -> Int -> ST s (HullIndex s) rebuildBuckets dense hull active requested = do   let !count = max 1 (min (hullBucketCapacity hull) requested)-  buckets <- MUV.replicate count noOuterEdge+  buckets <- MUV.replicate count noIndex   start <- denseReadFaceEdge dense 0   let go !remaining !edge         | remaining <= 0 = pure ()@@ -673,9 +695,9 @@   let !count = MUV.length buckets       !bucket = bucketFor count queryAngle   raw <- MUV.unsafeRead buckets bucket-  if raw == noOuterEdge+  if raw == noIndex     then denseReadFaceEdge dense 0-    else adjustFromBoundary bucket (fromIntegral raw)+    else adjustFromBoundary bucket (unpackIndex raw)  where   adjustFromBoundary !bucket !initial = do     initialAngle <- readAngle hull initial

@@ -26,7 +26,7 @@   , writeScratch   ) import Moonlight.Triangulation.Internal.Types-import Moonlight.Triangulation.Math (orient2d)+import Moonlight.Triangulation.Internal.Predicates (orient2d)  -- | The proved outgoing section of one ordinary removal. Its edges live in the -- operation scratch arena in counter-clockwise order; the record carries only

@@ -96,7 +96,10 @@   , RefinementParameters (..)   , unitElementDefaults   )-import Moonlight.Triangulation.Scalar (inCircleCoordinates, orient2dCoordinates)+import Moonlight.Triangulation.Internal.Predicates+  ( inCircleCoordinates+  , orient2dCoordinates+  ) import Moonlight.Triangulation.Math (distance, midpoint, validatePoint) import Moonlight.Triangulation.Internal.Transaction (runTransactionWithPublication) 

@@ -1,5 +1,4 @@ {-# LANGUAGE BangPatterns #-}-{-# LANGUAGE NamedFieldPuns #-}  module Moonlight.Triangulation.Internal.Location   ( MutableLocation (..)@@ -26,7 +25,11 @@   ( FaceId (..)   , UndirectedEdgeId (..)   )+import Moonlight.Triangulation.Internal.Predicates import Moonlight.Triangulation.Math+  ( canonicalPoint+  , squaredDistanceWide+  ) import Moonlight.Triangulation.Types  data MutableLocation@@ -89,7 +92,7 @@  -- | Greedy first-improvement descent through vertex neighbours: hop to the -- first neighbour strictly closer to the query until none improves. Distance--- is a heuristic only — exact containment belongs to 'walk' — so the widened+-- is a heuristic only — exact containment belongs to 'walkWithBarrier' — so the widened -- comparison can lengthen the path but never move the located answer. -- -- Strictly decreasing distance bounds the hops and nothing else: no vertex is@@ -99,7 +102,7 @@ -- step per outgoing half-edge, so a budget past 'directedEdgeCount' cannot be -- exhausted while the links are well formed. Exhausting it means they are not, -- and yields the vertex in hand rather than an error — a truncated descent can--- only lengthen 'walk', never move what it finds.+-- only lengthen 'walkWithBarrier', never move what it finds. descendToNearest   :: MutableDcel s vertex directed undirected face   -> OperationState s@@ -229,15 +232,6 @@                 _ -> go (edge + 1) (Just (terminal, candidate))    complementOne = -2--walk-  :: MutableDcel s vertex directed undirected face-  -> OperationState s-  -> Point-  -> Int-  -> ST s (Either BuildError MutableLocation)-walk mutable operation query initialFace = do-  walkWithBarrier Nothing mutable operation query initialFace  walkWithBarrier   :: Maybe IntSet.IntSet

@@ -85,16 +85,18 @@   , readScratch   , writeScratch   )-import Moonlight.Triangulation.Internal.PackedIndex (noIndex, packIndex)+import Moonlight.Triangulation.Internal.PackedIndex (noIndex, packIndex, unpackIndex)+import Moonlight.Triangulation.Internal.Predicates+  ( inCircle+  , inDiametralCircle+  , orient2d+  ) import Moonlight.Triangulation.Internal.Probe (Probe (..)) import Moonlight.Triangulation.Math   ( canonicalPoint   , circumcenter-  , inCircle-  , inDiametralCircle   , isFinite   , midpoint-  , orient2d   , squaredDistance   , squaredDistanceWide   , triangleArea@@ -950,7 +952,7 @@     count <- growableLength faces     traverse_       (\index -> do-         face <- fromIntegral <$> readGrowable faces index+         face <- unpackIndex <$> readGrowable faces index          retirePermittedFace workspace face       )       [0 .. count - 1]@@ -960,7 +962,7 @@     count <- growableLength faces     traverse       (\index -> do-         face <- fromIntegral <$> readGrowable faces index+         face <- unpackIndex <$> readGrowable faces index          signature <- mutableFaceSignature face          pure (face, signature)       )@@ -1239,7 +1241,7 @@   chainBoundary Workspace{wsCavityEpoch, wsCavityFaceMarks, wsCavityVertexMarks, wsCavityBoundary, wsCavityChain} count = do     epoch <- readSTRef wsCavityEpoch     walkBudget <- directedEdgeCount mutable-    start <- fromIntegral <$> readGrowable wsCavityBoundary 0+    start <- unpackIndex <$> readGrowable wsCavityBoundary 0     clearGrowable wsCavityChain     pushGrowable wsCavityChain (packIndex start)     walk epoch (walkBudget + 2) start start (count - 1)@@ -1300,7 +1302,7 @@           spoke !index !pair             | index >= boundaryCount = pure ()             | otherwise = do-                boundaryEdge <- fromIntegral <$> readGrowable wsCavityChain index+                boundaryEdge <- unpackIndex <$> readGrowable wsCavityChain index                 face <- fanFace faceBase index                 nextPair <-                   if index + 1 >= boundaryCount@@ -1329,7 +1331,7 @@       -- full shouldFlip pass would ask again what it already answered.       cocircular <- growableLength wsCavityCocircular       when (cocircular > 0) $ do-        ties <- traverse (\index -> fromIntegral <$> readGrowable wsCavityCocircular index) [0 .. cocircular - 1]+        ties <- traverse (\index -> unpackIndex <$> readGrowable wsCavityCocircular index) [0 .. cocircular - 1]         legalizeEdges mutable operation ties       admitStar workspace vertex       let offer :: Int -> ST s ()@@ -1345,12 +1347,12 @@      fanFace :: Int -> Int -> ST s Int     fanFace faceBase index-      | index < cavityFaceCount = fromIntegral <$> readGrowable wsCavityFaces index+      | index < cavityFaceCount = unpackIndex <$> readGrowable wsCavityFaces index       | otherwise = pure (faceBase + (index - cavityFaceCount))      fanPair :: Int -> Int -> ST s Int     fanPair edgeBase index-      | index < internalCount = fromIntegral <$> readGrowable wsCavityInternal index+      | index < internalCount = unpackIndex <$> readGrowable wsCavityInternal index       | otherwise = pure ((edgeBase + 2 * (index - internalCount)) `quot` 2)    -- A candidate outside the hull (only reachable when the hull is not

@@ -12,6 +12,7 @@   , insertVertex   , insertVertexAt   , insertVertexAtNear+  , insertVertexAtNearVertex   , insertAdmittedVertexAtCombining   , removeAt   , removeAtNear@@ -123,10 +124,7 @@   :: Point   -> vertex   -> Session s vertex directed undirected face (VertexId, InsertionDisposition)-insertVertexAt point payload = Session $ \mesh operation ->-  fmap-    (fmap (\(vertex, disposition) -> (VertexId (fromIntegral vertex), disposition)))-    (insertPointCombining (\_ replacement -> replacement) Nothing mesh operation point payload)+insertVertexAt = insertVertexAtFrom (const (pure Nothing))  -- | 'insertVertexAt' with the walk seeded at a face the caller vouches for -- -- typically the face a locate on the just-published value settled on, which@@ -138,17 +136,39 @@   -> Point   -> vertex   -> Session s vertex directed undirected face (VertexId, InsertionDisposition)-insertVertexAtNear (FaceId rawSeed) point payload = Session $ \mesh operation ->+insertVertexAtNear seed =+  insertVertexAtFrom (const (pure (Just (faceIdIndex seed))))++-- | 'insertVertexAt' with the exact walk restricted first to an inner face in+-- the supplied vertex's star. This is the native boundary for hierarchy and+-- nearest-neighbour hints, which answer vertices rather than faces. The seed+-- remains non-authoritative: an invalid handle or a vertex with no incident+-- inner face degrades to the ordinary unhinted descent.+insertVertexAtNearVertex+  :: VertexId+  -> Point+  -> vertex+  -> Session s vertex directed undirected face (VertexId, InsertionDisposition)+insertVertexAtNearVertex seed = insertVertexAtFrom (`interiorFaceNearVertex` seed)++insertVertexAtFrom+  :: (MutableDcel s vertex directed undirected face -> ST s (Maybe Int))+  -> Point+  -> vertex+  -> Session s vertex directed undirected face (VertexId, InsertionDisposition)+insertVertexAtFrom seedFor point payload = Session $ \mesh operation -> do+  seed <- seedFor mesh   fmap     (fmap (\(vertex, disposition) -> (VertexId (fromIntegral vertex), disposition)))     (insertPointCombining        (\_ replacement -> replacement)-       (Just (fromIntegral rawSeed))+       seed        mesh        operation        point        payload     )+{-# INLINE insertVertexAtFrom #-}  -- | Insert an admitted site while combining an occupied annotation. This is -- package-internal descent from a boundary that already owns coordinate@@ -205,9 +225,9 @@   :: VertexId   -> Point   -> Session s vertex directed undirected face (Maybe (RemovalOutcome vertex))-removeAtNear (VertexId rawGuess) point = Session $ \mesh operation -> do+removeAtNear guessVertex point = Session $ \mesh operation -> do   indexed <- identityIndexActive mesh-  let !guess = fromIntegral rawGuess+  let !guess = vertexIdIndex guessVertex       !query = canonicalPoint point   if indexed     then removeIndexed mesh operation query@@ -220,16 +240,35 @@           if stored == query             then fmap (fmap Just) (removeMutableOutcome mesh operation guess)             else do-              outgoing <- readVertexOut mesh guess-              let interiorFace edge fallback = do-                    face <- readFace mesh edge-                    if face > 0 then pure (Just face) else fallback-              hint <--                if outgoing < 0-                  then pure Nothing-                  else interiorFace outgoing (interiorFace (outgoing `xor` 1) (pure Nothing))+              hint <- interiorFaceAtVertex mesh guess               walkAndRemove mesh operation hint query +interiorFaceNearVertex+  :: MutableDcel s vertex directed undirected face+  -> VertexId+  -> ST s (Maybe Int)+interiorFaceNearVertex mesh vertex@(VertexId rawVertex) = do+  vertices <- pointCount mesh+  if toInteger rawVertex >= toInteger vertices+    then pure Nothing+    else interiorFaceAtVertex mesh (vertexIdIndex vertex)++interiorFaceAtVertex+  :: MutableDcel s vertex directed undirected face+  -> Int+  -> ST s (Maybe Int)+interiorFaceAtVertex mesh vertex = do+  outgoing <- readVertexOut mesh vertex+  if outgoing < 0+    then pure Nothing+    else do+      primary <- innerFaceOf outgoing+      maybe (innerFaceOf (outgoing `xor` 1)) (pure . Just) primary+ where+  innerFaceOf edge = do+    face <- readFace mesh edge+    pure (if face > 0 then Just face else Nothing)+ walkAndRemove   :: MutableDcel s vertex directed undirected face   -> OperationState s@@ -314,7 +353,7 @@     then pure (Left (RemovalVertexOutOfRange requested vertices))     else do       activatePointIndex mesh-      removeMutableOutcome mesh operation (fromIntegral raw)+      removeMutableOutcome mesh operation (vertexIdIndex requested)  removeMutableOutcome   :: MutableDcel s vertex directed undirected face

@@ -37,7 +37,9 @@   , UndirectedEdgeId (..)   , asUndirected   , directedPair+  , faceIdIndex   , reverseEdge+  , undirectedEdgeIdIndex   ) import Moonlight.Triangulation.Internal.Mutable import Moonlight.Triangulation.Internal.Paged@@ -120,15 +122,12 @@       [ pair       | rawFace <- IntSet.toAscList permitted       , edge <- faceDirectedEdges triangulation (FaceId (fromIntegral rawFace))-      , let pair = fromIntegral (unUndirected (asUndirected edge))+      , let pair = undirectedEdgeIdIndex (asUndirected edge)       , let adjacent = incidentFace triangulation (reverseEdge edge)       , adjacent /= FaceId 0-      , IntSet.notMember (fromIntegral (unFace adjacent)) permitted+      , IntSet.notMember (faceIdIndex adjacent) permitted       ] -  unUndirected (UndirectedEdgeId raw) = raw-  unFace (FaceId raw) = raw- -- | Refine exactly one checked local section. Interface edges are installed as -- transaction-local legalization barriers and removed before publication; -- the receipt and collar proof inspect only the admitted local closure.@@ -223,8 +222,8 @@         if domain == Nothing && refineExcludeOuterFaces parameters           then             IntSet.fromList-              [ fromIntegral raw-              | FaceId raw <-+              [ faceIdIndex face+              | face <-                   facesAtEvenBarrierDepth                     triangulation                     (isConstraintEdge triangulation)@@ -373,12 +372,10 @@   -> IntSet.IntSet localDomainEdges triangulation faces =   IntSet.fromList-    [ fromIntegral (unUndirected (asUndirected edge))+    [ undirectedEdgeIdIndex (asUndirected edge)     | rawFace <- IntSet.toAscList faces     , edge <- faceDirectedEdges triangulation (FaceId (fromIntegral rawFace))     ]- where-  unUndirected (UndirectedEdgeId raw) = raw  localEdgeSignature   :: Triangulation mode vertex directed undirected face@@ -436,8 +433,8 @@       , let edge = UndirectedEdgeId (fromIntegral rawPair)       , let (forward, backward) = directedPair edge       , rawFace <--          [ fromIntegral (unFace (incidentFace before forward))-          , fromIntegral (unFace (incidentFace before backward))+          [ faceIdIndex (incidentFace before forward)+          , faceIdIndex (incidentFace before backward)           ]       , rawFace > 0       , IntSet.notMember rawFace permitted@@ -464,7 +461,7 @@         | protectedFace `elem` afterFaces ->             case filter (/= protectedFace) afterFaces of               [oppositeFace]-                | IntSet.member (fromIntegral (unFace oppositeFace)) finalPermitted ->+                | IntSet.member (faceIdIndex oppositeFace) finalPermitted ->                     Right (edge, protectedFace, oppositeFace)                 | otherwise ->                     Left@@ -476,11 +473,9 @@       _ -> Left RefinementDomainTopologyChanged    protected face =-    let raw = fromIntegral (unFace face)+    let raw = faceIdIndex face      in raw > 0 && IntSet.notMember raw permitted -  unFace (FaceId raw) = raw- -- | Refuse outer or absent faces rather than silently treating an invalid -- topology witness as an empty repair. Face zero is the outer face and has no -- active refinement equation.@@ -494,10 +489,10 @@   totalFaces = numFaces triangulation    validateFace :: FaceId -> Either BuildError Int-  validateFace face@(FaceId raw)-    | raw == 0 || toInteger raw >= toInteger totalFaces =+  validateFace face+    | face == FaceId 0 || toInteger (unFaceId face) >= toInteger totalFaces =         Left (RefinementSeedFaceNotActive face totalFaces)-    | otherwise = Right (fromIntegral raw)+    | otherwise = Right (faceIdIndex face)  -- | Refuse parameters no mesh can satisfy. validateRefinementParameters :: RefinementParameters -> Either BuildError ()

@@ -6,6 +6,7 @@   , insertVertex   , insertVertexAt   , insertVertexAtNear+  , insertVertexAtNearVertex   , removeAt   , removeAtNear   , removeManyAt

@@ -0,0 +1,131 @@+-- | Homology interpretations of canonical labelled planar complexes.+--+-- 'PlanarComplex' owns simplex identity and downward closure. This module+-- materializes that admitted section once; Homology remains the sole owner of+-- chain-law and filtration-law validation.+module Moonlight.Triangulation.LabelledComplex+  ( PlanarComplexInterpretationError (..)+  , finitePlanarComplex+  , filteredPlanarComplex+  )+where++import Data.Bifunctor (first)+import Data.Map.Strict (Map)+import Data.Map.Strict qualified as Map+import Data.Vector qualified as Vector+import Moonlight.Homology.Boundary+  ( BoundaryIncidence+  , FiniteChainComplex+  , emptyBoundaryIncidence+  , emptyBoundaryIncidenceOf+  , materializeIncidenceBoundary+  , mkFiniteChainComplexChecked+  )+import Moonlight.Homology.Chain+  ( HomologicalDegree (..)+  , HomologyFailure+  )+import Moonlight.Homology.Persistence+  ( FilteredFiniteChainComplex+  , mkFilteredFiniteChainComplex+  )+import Moonlight.Homology.Topology (BasisCellRef (..))+import Moonlight.Triangulation.Simplex+  ( PlanarComplex+  , PlanarSimplex+  , SimplexDimension (..)+  , planarComplexBasisAt+  , planarSimplexBoundary+  )++-- | Typed obstructions at the labelled-complex/Homology boundary.+data PlanarComplexInterpretationError label+  = PlanarComplexChainInvalid !HomologyFailure+  | PlanarComplexBirthUnexpected !(PlanarSimplex label)+  | PlanarComplexFiltrationInvalid !HomologyFailure+  deriving stock (Eq, Show)++-- | Canonical integral chains in ascending labelled-simplex order.+finitePlanarComplex+  :: Ord label+  => PlanarComplex label+  -> Either (PlanarComplexInterpretationError label) (FiniteChainComplex Int)+finitePlanarComplex complexValue = do+  degreeOneBoundary <-+    first PlanarComplexChainInvalid+      ( materializeIncidenceBoundary+          planarSimplexBoundary+          (basisList SimplexDimension1 complexValue)+          (basisList SimplexDimension0 complexValue)+      )+  degreeTwoBoundary <-+    first PlanarComplexChainInvalid+      ( materializeIncidenceBoundary+          planarSimplexBoundary+          (basisList SimplexDimension2 complexValue)+          (basisList SimplexDimension1 complexValue)+      )+  let degreeZeroBoundary :: BoundaryIncidence Int+      degreeZeroBoundary =+        emptyBoundaryIncidenceOf+          (fromIntegral (Vector.length (planarComplexBasisAt SimplexDimension0 complexValue)))+          0+      boundaryAt (HomologicalDegree degreeIndex) =+        case degreeIndex of+          0 -> degreeZeroBoundary+          1 -> degreeOneBoundary+          2 -> degreeTwoBoundary+          _ -> emptyBoundaryIncidence+  first PlanarComplexChainInvalid+    (mkFiniteChainComplexChecked (HomologicalDegree 2) boundaryAt)++-- | Add exact births to the admitted complex. Homology validates uniqueness,+-- total coverage, and every boundary inequality; this interpreter deliberately+-- does not duplicate those laws.+filteredPlanarComplex+  :: (Ord label, Ord filtration)+  => PlanarComplex label+  -> Map (PlanarSimplex label) filtration+  -> Either+      (PlanarComplexInterpretationError label)+      (FilteredFiniteChainComplex filtration Int)+filteredPlanarComplex complexValue births = do+  finite <- finitePlanarComplex complexValue+  assignments <- traverse birthAssignment (Map.toAscList births)+  first PlanarComplexFiltrationInvalid+    (mkFilteredFiniteChainComplex finite assignments)+ where+  basisReferences = planarBasisReferences complexValue+  birthAssignment (simplex, birthValue) =+    maybe+      (Left (PlanarComplexBirthUnexpected simplex))+      (Right . (,birthValue))+      (Map.lookup simplex basisReferences)++planarBasisReferences+  :: Ord label+  => PlanarComplex label+  -> Map (PlanarSimplex label) BasisCellRef+planarBasisReferences complexValue =+  Map.unions+    [ referencesAt 0 SimplexDimension0+    , referencesAt 1 SimplexDimension1+    , referencesAt 2 SimplexDimension2+    ]+ where+  referencesAt degreeValue dimension =+    Map.fromDistinctAscList+      ( zipWith+          (\basisIndex simplex ->+              (simplex, BasisCellRef (HomologicalDegree degreeValue) basisIndex)+          )+          [0 ..]+          (basisList dimension complexValue)+      )++basisList+  :: SimplexDimension+  -> PlanarComplex label+  -> [PlanarSimplex label]+basisList dimension = Vector.toList . planarComplexBasisAt dimension

@@ -13,6 +13,7 @@   , exactRationalFromFiniteDouble   , exactRationalFromDyadic   , exactRationalFromDyadicRatio+  , exactRationalFromNormalizedRatio   , exactRationalNumerator   , exactRationalDenominator   , exactRationalBitWidth@@ -114,6 +115,13 @@         then exactRational (oddNumerator `shiftL` residualPower) oddDenominator         else exactRational oddNumerator (oddDenominator `shiftL` negate residualPower) {-# INLINE exactRationalFromDyadicRatio #-}++-- | Internal bridge from the normalized carrier owned by @Data.Ratio@. This+-- exists for statically nonzero rational constants in exact kernels; public+-- callers continue through 'exactRational'.+exactRationalFromNormalizedRatio :: Ratio Integer -> ExactRational+exactRationalFromNormalizedRatio = ExactRational+{-# INLINE exactRationalFromNormalizedRatio #-}  -- | Read the reduced numerator. exactRationalNumerator :: ExactRational -> Integer

@@ -9,8 +9,6 @@   , scalarUnitRoundoff   , scalarCcwErrorBound   , scalarInCircleErrorBound-  , orient2dCoordinates-  , inCircleCoordinates   , circumradiusSquaredWithinCoordinates   , BinaryFormat   , formatRadix@@ -24,8 +22,6 @@  import Moonlight.Triangulation.Internal.Dyadic   ( exactCircumradiusSquaredWithin-  , exactInCircleDet-  , exactOrientSignDouble   )  -- | The coordinate component of canonical point identity. IEEE signed zeros@@ -81,21 +77,6 @@ scalarInCircleErrorBound :: Double scalarInCircleErrorBound = 1.1102230246251577e-15 --- | Exact orientation ordering of three binary64 coordinate pairs.-orient2dCoordinates-  :: Double -> Double -> Double -> Double -> Double -> Double-  -> Ordering-orient2dCoordinates = filteredOrient2dDouble-{-# INLINE orient2dCoordinates #-}---- | Exact in-circle ordering of four binary64 coordinate pairs.-inCircleCoordinates-  :: Double -> Double -> Double -> Double-  -> Double -> Double -> Double -> Double-  -> Ordering-inCircleCoordinates = filteredInCircle scalarInCircleErrorBound-{-# INLINE inCircleCoordinates #-}- -- | Closed exact circumradius membership at a finite, non-negative binary64 -- threshold. Invalid thresholds and collinear triples are outside. circumradiusSquaredWithinCoordinates@@ -133,51 +114,6 @@ isFinite :: Double -> Bool isFinite value = value - value == 0 {-# INLINE isFinite #-}---- The binary64 kernel pairs the approximation test with the--- fixed-precision exact sign, which answers the dyadic determinant's sign in--- machine words rather than allocated Integers whenever the exponent spread--- allows, and defers to the dyadic determinant when it does not.-filteredOrient2dDouble-  :: Double -> Double -> Double -> Double -> Double -> Double -> Ordering-filteredOrient2dDouble ax ay bx by cx cy-  | abs determinant > errorBound * determinantSum = compare determinant 0-  | otherwise = exactOrientSignDouble ax ay bx by cx cy- where-  errorBound = 3.3306690738754716e-16-  !left = (ax - cx) * (by - cy)-  !right = (ay - cy) * (bx - cx)-  !determinant = left - right-  !determinantSum = abs left + abs right-{-# INLINE filteredOrient2dDouble #-}--filteredInCircle-  :: Double-  -> Double -> Double -> Double -> Double-  -> Double -> Double -> Double -> Double-  -> Ordering-filteredInCircle errorBound ax ay bx by cx cy dx dy-  | abs determinant > errorBound * permanent = compare determinant 0-  | otherwise = compare (exactInCircleDet ax ay bx by cx cy dx dy) 0- where-  !adx = ax - dx-  !ady = ay - dy-  !bdx = bx - dx-  !bdy = by - dy-  !cdx = cx - dx-  !cdy = cy - dy-  !abdet = adx * bdy - bdx * ady-  !bcdet = bdx * cdy - cdx * bdy-  !cadet = cdx * ady - adx * cdy-  !alift = adx * adx + ady * ady-  !blift = bdx * bdx + bdy * bdy-  !clift = cdx * cdx + cdy * cdy-  !determinant = alift * bcdet + blift * cadet + clift * abdet-  !permanent =-    (abs (bdx * cdy) + abs (cdx * bdy)) * alift-      + (abs (cdx * ady) + abs (adx * cdy)) * blift-      + (abs (adx * bdy) + abs (bdx * ady)) * clift-{-# INLINE filteredInCircle #-}  -- | The smallest coordinate magnitude the exact predicates accept. minimumAllowedCoordinate :: Double

@@ -44,6 +44,9 @@   , UndirectedEdgeId (..)   , VertexId (..)   , directedPair+  , faceIdIndex+  , undirectedEdgeIdIndex+  , vertexIdIndex   ) import Moonlight.Triangulation.Internal.CellSet   ( ExactCellSet (..)@@ -66,7 +69,7 @@   ( DelaunayTriangulation   , Triangulation   )-import Moonlight.Triangulation.Math (inDiametralCircle)+import Moonlight.Triangulation.Internal.Predicates (inDiametralCircle) import Moonlight.Triangulation.Scalar (circumradiusSquaredWithinCoordinates) import Moonlight.Triangulation.Types   ( ConstraintMode (Unconstrained)@@ -188,14 +191,14 @@ -- | Exact birth of a resident vertex, when the handle belongs to this -- filtration. alphaVertexBirth :: AlphaFiltration -> VertexId -> Maybe AlphaBirth-alphaVertexBirth filtration (VertexId rawVertex) =-  alphaVertexBirths filtration Vector.!? fromIntegral rawVertex+alphaVertexBirth filtration vertex =+  alphaVertexBirths filtration Vector.!? vertexIdIndex vertex  -- | Exact birth of a resident undirected edge, when the handle belongs to this -- filtration. alphaEdgeBirth :: AlphaFiltration -> UndirectedEdgeId -> Maybe AlphaBirth-alphaEdgeBirth filtration (UndirectedEdgeId rawEdge) =-  alphaEdgeBirths filtration Vector.!? fromIntegral rawEdge+alphaEdgeBirth filtration edge =+  alphaEdgeBirths filtration Vector.!? undirectedEdgeIdIndex edge  -- | Exact birth of a resident bounded face, when the handle belongs to this -- filtration.@@ -374,14 +377,14 @@     Just nonEmptyBirths -> Right (minimum nonEmptyBirths)  faceSectionAt :: Vector FaceAlphaSection -> FaceId -> Maybe FaceAlphaSection-faceSectionAt faceSections (FaceId rawFace)-  | rawFace == 0 = Nothing-  | otherwise = faceSections Vector.!? (fromIntegral rawFace - 1)+faceSectionAt faceSections face+  | face == Dcel.outerFace = Nothing+  | otherwise = faceSections Vector.!? (faceIdIndex face - 1)  faceBirthAt :: Vector AlphaBirth -> FaceId -> Maybe AlphaBirth-faceBirthAt faceBirths (FaceId rawFace)-  | rawFace == 0 = Nothing-  | otherwise = faceBirths Vector.!? (fromIntegral rawFace - 1)+faceBirthAt faceBirths face+  | face == Dcel.outerFace = Nothing+  | otherwise = faceBirths Vector.!? (faceIdIndex face - 1)  triangleVertices :: (VertexId, VertexId, VertexId) -> [VertexId] triangleVertices (firstVertex, secondVertex, thirdVertex) =

@@ -97,8 +97,8 @@  -- | Authoritative geometric position of an admitted vertex handle. vertexPoint :: Triangulation mode vertex directed undirected face -> VertexId -> Point-vertexPoint triangulation (VertexId vertex) =-  let !index = fromIntegral vertex+vertexPoint triangulation vertex =+  let !index = vertexIdIndex vertex    in Point (pagedUnsafeIndex (triPointX triangulation) index) (pagedUnsafeIndex (triPointY triangulation) index) {-# INLINE vertexPoint #-} @@ -110,26 +110,26 @@  -- | Annotation carried by an admitted vertex handle. vertexData :: Triangulation mode vertex directed undirected face -> VertexId -> vertex-vertexData triangulation (VertexId vertex) =-  boxedUnsafeIndex (triVertexData triangulation) (fromIntegral vertex)+vertexData triangulation vertex =+  boxedUnsafeIndex (triVertexData triangulation) (vertexIdIndex vertex) {-# INLINE vertexData #-}  -- | Annotation carried by an admitted directed edge. directedEdgeData :: Triangulation mode vertex directed undirected face -> DirectedEdgeId -> directed-directedEdgeData triangulation (DirectedEdgeId edge) =-  boxedUnsafeIndex (triDirectedData triangulation) (fromIntegral edge)+directedEdgeData triangulation edge =+  boxedUnsafeIndex (triDirectedData triangulation) (directedEdgeIdIndex edge) {-# INLINE directedEdgeData #-}  -- | Annotation carried by an admitted undirected edge. undirectedEdgeData :: Triangulation mode vertex directed undirected face -> UndirectedEdgeId -> undirected-undirectedEdgeData triangulation (UndirectedEdgeId edge) =-  boxedUnsafeIndex (triUndirectedData triangulation) (fromIntegral edge)+undirectedEdgeData triangulation edge =+  boxedUnsafeIndex (triUndirectedData triangulation) (undirectedEdgeIdIndex edge) {-# INLINE undirectedEdgeData #-}  -- | Annotation carried by an admitted face. faceData :: Triangulation mode vertex directed undirected face -> FaceId -> face-faceData triangulation (FaceId face) =-  boxedUnsafeIndex (triFaceData triangulation) (fromIntegral face)+faceData triangulation face =+  boxedUnsafeIndex (triFaceData triangulation) (faceIdIndex face) {-# INLINE faceData #-}  -- | Replace one vertex annotation without changing geometry or topology.@@ -138,8 +138,8 @@   -> VertexId   -> vertex   -> Triangulation mode vertex directed undirected face-setVertexData triangulation (VertexId vertex) payload =-  triangulation{triVertexData = boxedUpdate (fromIntegral vertex) payload (triVertexData triangulation)}+setVertexData triangulation vertex payload =+  triangulation{triVertexData = boxedUpdate (vertexIdIndex vertex) payload (triVertexData triangulation)}  -- | Replace one directed-edge annotation without changing geometry or topology. setDirectedEdgeData@@ -147,8 +147,8 @@   -> DirectedEdgeId   -> directed   -> Triangulation mode vertex directed undirected face-setDirectedEdgeData triangulation (DirectedEdgeId edge) payload =-  triangulation{triDirectedData = boxedUpdate (fromIntegral edge) payload (triDirectedData triangulation)}+setDirectedEdgeData triangulation edge payload =+  triangulation{triDirectedData = boxedUpdate (directedEdgeIdIndex edge) payload (triDirectedData triangulation)}  -- | Replace one undirected-edge annotation without changing geometry or topology. setUndirectedEdgeData@@ -156,8 +156,8 @@   -> UndirectedEdgeId   -> undirected   -> Triangulation mode vertex directed undirected face-setUndirectedEdgeData triangulation (UndirectedEdgeId edge) payload =-  triangulation{triUndirectedData = boxedUpdate (fromIntegral edge) payload (triUndirectedData triangulation)}+setUndirectedEdgeData triangulation edge payload =+  triangulation{triUndirectedData = boxedUpdate (undirectedEdgeIdIndex edge) payload (triUndirectedData triangulation)}  -- | Replace one face annotation without changing geometry or topology. setFaceData@@ -165,25 +165,25 @@   -> FaceId   -> face   -> Triangulation mode vertex directed undirected face-setFaceData triangulation (FaceId face) payload =-  triangulation{triFaceData = boxedUpdate (fromIntegral face) payload (triFaceData triangulation)}+setFaceData triangulation face payload =+  triangulation{triFaceData = boxedUpdate (faceIdIndex face) payload (triFaceData triangulation)}  -- | One outgoing directed edge, if the vertex is connected. vertexOutEdge :: Triangulation mode vertex directed undirected face -> VertexId -> Maybe DirectedEdgeId-vertexOutEdge triangulation (VertexId vertex) =-  DirectedEdgeId . fromIntegral <$> unpackOptionalIndex (pagedUnsafeIndex (triVertexOut triangulation) (fromIntegral vertex))+vertexOutEdge triangulation vertex =+  DirectedEdgeId . fromIntegral <$> unpackOptionalIndex (pagedUnsafeIndex (triVertexOut triangulation) (vertexIdIndex vertex)) {-# INLINE vertexOutEdge #-}  -- | One boundary edge of a face, if the face has a boundary. adjacentEdge :: Triangulation mode vertex directed undirected face -> FaceId -> Maybe DirectedEdgeId-adjacentEdge triangulation (FaceId face) =-  DirectedEdgeId . fromIntegral <$> unpackOptionalIndex (pagedUnsafeIndex (triFaceEdge triangulation) (fromIntegral face))+adjacentEdge triangulation face =+  DirectedEdgeId . fromIntegral <$> unpackOptionalIndex (pagedUnsafeIndex (triFaceEdge triangulation) (faceIdIndex face)) {-# INLINE adjacentEdge #-}  -- | Origin vertex of an admitted directed edge. origin :: Triangulation mode vertex directed undirected face -> DirectedEdgeId -> VertexId-origin triangulation (DirectedEdgeId edge) =-  VertexId (pagedUnsafeIndex (triHalfTopology triangulation) (4 * fromIntegral edge))+origin triangulation edge =+  VertexId (pagedUnsafeIndex (triHalfTopology triangulation) (4 * directedEdgeIdIndex edge)) {-# INLINE origin #-}  -- | Destination vertex of an admitted directed edge.@@ -193,20 +193,20 @@  -- | Next directed edge around the incident face. next :: Triangulation mode vertex directed undirected face -> DirectedEdgeId -> DirectedEdgeId-next triangulation (DirectedEdgeId edge) =-  DirectedEdgeId (pagedUnsafeIndex (triHalfTopology triangulation) (4 * fromIntegral edge + 1))+next triangulation edge =+  DirectedEdgeId (pagedUnsafeIndex (triHalfTopology triangulation) (4 * directedEdgeIdIndex edge + 1)) {-# INLINE next #-}  -- | Previous directed edge around the incident face. previous :: Triangulation mode vertex directed undirected face -> DirectedEdgeId -> DirectedEdgeId-previous triangulation (DirectedEdgeId edge) =-  DirectedEdgeId (pagedUnsafeIndex (triHalfTopology triangulation) (4 * fromIntegral edge + 2))+previous triangulation edge =+  DirectedEdgeId (pagedUnsafeIndex (triHalfTopology triangulation) (4 * directedEdgeIdIndex edge + 2)) {-# INLINE previous #-}  -- | Face on the left of an admitted directed edge. incidentFace :: Triangulation mode vertex directed undirected face -> DirectedEdgeId -> FaceId-incidentFace triangulation (DirectedEdgeId edge) =-  FaceId (pagedUnsafeIndex (triHalfTopology triangulation) (4 * fromIntegral edge + 3))+incidentFace triangulation edge =+  FaceId (pagedUnsafeIndex (triHalfTopology triangulation) (4 * directedEdgeIdIndex edge + 3)) {-# INLINE incidentFace #-}  -- | Whether the directed edge is incident to the unbounded face.@@ -223,8 +223,8 @@  -- | Whether the edge belongs to the constrained-edge section. isConstraintEdge :: Triangulation mode vertex directed undirected face -> UndirectedEdgeId -> Bool-isConstraintEdge triangulation (UndirectedEdgeId edge) =-  pagedUnsafeIndex (triConstraint triangulation) (fromIntegral edge) /= (0 :: Word8)+isConstraintEdge triangulation edge =+  pagedUnsafeIndex (triConstraint triangulation) (undirectedEdgeIdIndex edge) /= (0 :: Word8) {-# INLINE isConstraintEdge #-}  -- | Number of constrained undirected edges.@@ -254,8 +254,8 @@  -- | The three directed edges of a bounded triangular face. innerFaceDirectedEdges :: Triangulation mode vertex directed undirected face -> FaceId -> Maybe (DirectedEdgeId, DirectedEdgeId, DirectedEdgeId)-innerFaceDirectedEdges triangulation face@(FaceId rawFace)-  | face == outerFace || fromIntegral rawFace >= numFaces triangulation = Nothing+innerFaceDirectedEdges triangulation face+  | face == outerFace || faceIdIndex face >= numFaces triangulation = Nothing   | otherwise = do       e0 <- adjacentEdge triangulation face       let !e1 = next triangulation e0

@@ -4,7 +4,14 @@  -- | Exact rational planar geometry over admitted binary64 points. module Moonlight.Triangulation.Exact-  ( ExactPoint+  ( ExactRational+  , ExactArithmeticError (..)+  , exactRational+  , exactRationalNumerator+  , exactRationalDenominator+  , exactRationalBitWidth+  , exactDivide+  , ExactPoint   , exactPoint   , exactPointCoordinates   , exactPointCross

@@ -2,7 +2,6 @@ {-# LANGUAGE EmptyDataDecls #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE KindSignatures #-}-{-# LANGUAGE NamedFieldPuns #-} {-# LANGUAGE RoleAnnotations #-}  -- | Owning handles that prevent identifiers from crossing mesh boundaries.@@ -66,6 +65,7 @@  import Moonlight.Triangulation.Dcel qualified as Dcel import Moonlight.Triangulation.Internal.HandleDefs+import Moonlight.Triangulation.Internal.Predicates qualified as Admitted import Moonlight.Triangulation.LineSideInfo (LineSideInfo) import Moonlight.Triangulation.Math qualified as Math import Moonlight.Triangulation.Types@@ -129,8 +129,8 @@   :: Triangulation mode vertex directed undirected face   -> VertexId   -> Maybe (VertexHandle mode vertex directed undirected face)-vertexHandle triangulation vertex@(VertexId raw)-  | fromIntegral raw < Dcel.numVertices triangulation = Just (VertexHandle triangulation vertex)+vertexHandle triangulation vertex+  | vertexIdIndex vertex < Dcel.numVertices triangulation = Just (VertexHandle triangulation vertex)   | otherwise = Nothing  -- | Admit a directed-edge identifier into a triangulation.@@ -138,8 +138,8 @@   :: Triangulation mode vertex directed undirected face   -> DirectedEdgeId   -> Maybe (DirectedEdgeHandle mode vertex directed undirected face)-directedEdgeHandle triangulation edge@(DirectedEdgeId raw)-  | fromIntegral raw < Dcel.numDirectedEdges triangulation = Just (DirectedEdgeHandle triangulation edge)+directedEdgeHandle triangulation edge+  | directedEdgeIdIndex edge < Dcel.numDirectedEdges triangulation = Just (DirectedEdgeHandle triangulation edge)   | otherwise = Nothing  -- | Admit an undirected-edge identifier into a triangulation.@@ -147,8 +147,8 @@   :: Triangulation mode vertex directed undirected face   -> UndirectedEdgeId   -> Maybe (UndirectedEdgeHandle mode vertex directed undirected face)-undirectedEdgeHandle triangulation edge@(UndirectedEdgeId raw)-  | fromIntegral raw < Dcel.numUndirectedEdges triangulation = Just (UndirectedEdgeHandle triangulation edge)+undirectedEdgeHandle triangulation edge+  | undirectedEdgeIdIndex edge < Dcel.numUndirectedEdges triangulation = Just (UndirectedEdgeHandle triangulation edge)   | otherwise = Nothing  -- | Admit a face identifier that may denote the outer face.@@ -156,8 +156,8 @@   :: Triangulation mode vertex directed undirected face   -> FaceId   -> Maybe (FaceHandle PossiblyOuterTag mode vertex directed undirected face)-faceHandle triangulation face@(FaceId raw)-  | fromIntegral raw < Dcel.numFaces triangulation = Just (FaceHandle triangulation (FixedFaceHandle face))+faceHandle triangulation face+  | faceIdIndex face < Dcel.numFaces triangulation = Just (FaceHandle triangulation (FixedFaceHandle face))   | otherwise = Nothing  -- | Admit a face identifier while proving that it is bounded.@@ -293,11 +293,11 @@ -- | Exact side of the owning edge's oriented line. directedEdgeSideQuery   :: DirectedEdgeHandle mode vertex directed undirected face-  -> Point+  -> QueryPoint   -> LineSideInfo directedEdgeSideQuery edge query =   let (from, to) = directedEdgePositions edge-   in Math.sideQuery from to query+   in Admitted.sideQuery from to (queryPointValue query)  -- | Undirected-edge annotation through an owning handle. undirectedEdgeDataH :: UndirectedEdgeHandle mode vertex directed undirected face -> undirected

@@ -42,6 +42,9 @@   , UndirectedEdgeId (..)   , VertexId (..)   , asUndirected+  , faceIdIndex+  , undirectedEdgeIdIndex+  , vertexIdIndex   ) import Moonlight.Triangulation.Internal.Representation (Triangulation) import Moonlight.Triangulation.Internal.Types (PointValidationError)@@ -107,9 +110,9 @@   exactPoints <-     IntMap.fromDistinctAscList       <$> traverse-        ( \vertexIndex -> do-            point <- ordinaryExactPoint triangulation (VertexId (fromIntegral vertexIndex))-            pure (vertexIndex, point)+        ( \residentIndex -> do+            point <- ordinaryExactPoint triangulation (VertexId (fromIntegral residentIndex))+            pure (residentIndex, point)         )         vertexIndices   pure@@ -226,15 +229,15 @@   traverse_ validateEdge selectedEdges   traverse_ validateFace selectedFaces  where-  validateVertex vertex@(VertexId raw)-    | fromIntegral raw < numVertices triangulation = Right ()+  validateVertex vertex+    | vertexIdIndex vertex < numVertices triangulation = Right ()     | otherwise = Left (CellVertexOutOfRange vertex (numVertices triangulation))-  validateEdge edge@(UndirectedEdgeId raw)-    | fromIntegral raw < numUndirectedEdges triangulation = Right ()+  validateEdge edge+    | undirectedEdgeIdIndex edge < numUndirectedEdges triangulation = Right ()     | otherwise = Left (CellEdgeOutOfRange edge (numUndirectedEdges triangulation))-  validateFace face@(FaceId raw)-    | raw == 0 = Left CellOuterFaceSelected-    | fromIntegral raw < numFaces triangulation = Right ()+  validateFace face+    | faceIdIndex face == 0 = Left CellOuterFaceSelected+    | faceIdIndex face < numFaces triangulation = Right ()     | otherwise = Left (CellFaceOutOfRange face (numFaces triangulation))  validateEdgeClosure@@ -246,8 +249,8 @@   traverse_ requireVertex [from, to]  where   (from, to) = undirectedEndpoints triangulation edge-  requireVertex vertex@(VertexId raw)-    | IntSet.member (fromIntegral raw) selectedVertices = Right ()+  requireVertex vertex+    | IntSet.member (vertexIdIndex vertex) selectedVertices = Right ()     | otherwise = Left (CellEdgeBoundaryMissing edge vertex)  validateFaceClosure@@ -258,18 +261,18 @@ validateFaceClosure triangulation selectedEdges face =   traverse_ (requireEdge . asUndirected) (faceDirectedEdges triangulation face)  where-  requireEdge edge@(UndirectedEdgeId raw)-    | IntSet.member (fromIntegral raw) selectedEdges = Right ()+  requireEdge edge+    | IntSet.member (undirectedEdgeIdIndex edge) selectedEdges = Right ()     | otherwise = Left (CellFaceEdgeMissing face edge)  vertexSet :: [VertexId] -> IntSet.IntSet-vertexSet = IntSet.fromList . map (\(VertexId raw) -> fromIntegral raw)+vertexSet = IntSet.fromList . map vertexIdIndex  edgeSet :: [UndirectedEdgeId] -> IntSet.IntSet-edgeSet = IntSet.fromList . map (\(UndirectedEdgeId raw) -> fromIntegral raw)+edgeSet = IntSet.fromList . map undirectedEdgeIdIndex  faceSet :: [FaceId] -> IntSet.IntSet-faceSet = IntSet.fromList . map (\(FaceId raw) -> fromIntegral raw)+faceSet = IntSet.fromList . map faceIdIndex  exactCellSetVertexCount :: ExactCellSet -> Int exactCellSetVertexCount (ExactCellSet _ selected _ _) = IntMap.size selected

@@ -11,6 +11,8 @@   , setupSecondVertex   , splitLineEdge   , extendLine+  , prepareAscendingLineArena+  , connectAscendingLineArena   , collectLineChain   , lineToArea   ) where@@ -18,14 +20,28 @@ import Control.Monad (forM_, when) import Control.Monad.ST (ST) import Data.STRef (writeSTRef)+import qualified Data.Vector.Unboxed.Mutable as MUV+import Data.Word (Word32) import Moonlight.Triangulation.Internal.HandleDefs (VertexId (..)) import Moonlight.Triangulation.Internal.DcelOperations.Legalize (legalizeScratch) import Moonlight.Triangulation.Internal.DcelOperations.Twin (reverseIndex)+import Moonlight.Triangulation.Internal.PackedIndex (unpackIndex) import Moonlight.Triangulation.Internal.Mutable   ( MutableDcel (..)+  , DenseMutableDcel   , addEdge   , addEdgeBlock   , addFaceBlock+  , denseAddEdgeBlock+  , denseCommitFreshConnections+  , denseMarkFreshConnected+  , denseMutableOwner+  , denseReadPointX+  , denseReadPointY+  , denseWriteFace+  , denseWriteFaceEdge+  , denseLinkEdges+  , denseWriteOrigin   , directedEdgeCount   , edgeOriginPoint   , ensureCellCapacity@@ -56,7 +72,7 @@   ) import Moonlight.Triangulation.Internal.Probe (KnownProbe) import Moonlight.Triangulation.Internal.Types (BuildError (..))-import Moonlight.Triangulation.Math (orient2d)+import Moonlight.Triangulation.Internal.Predicates (orient2d, orient2dCoordinates)  setupFirstVertex :: MutableDcel s vertex directed undirected face -> Int -> ST s () setupFirstVertex mutable vertex = markConnected mutable vertex (-1)@@ -146,6 +162,145 @@           markConnected mutable vertex newTwin           addCounter operation CounterLineExtensions 1           pure (Right ())++-- | A consumed radial arena whose sites have been proved collinear and placed+-- in ascending coordinate order. The constructor stays here with the chain+-- topology it licenses; the sweep can only obtain the witness through+-- 'prepareAscendingLineArena'.+newtype AscendingLineArena s = AscendingLineArena (MUV.MVector s Word32)++-- | Recognize a one-dimensional fresh section and derive its chain order from+-- the radial order already paid for by bulk admission. The nearest radial site+-- is a lawful pivot: distances increase monotonically on each side of it, so+-- reversing the lower fibre and retaining the upper fibre yields total line+-- order without another sort.+--+-- The third point is inspected before the remaining fold. General-position+-- input therefore pays one exact orientation, not a speculative linear scan.+prepareAscendingLineArena+  :: forall s vertex directed undirected face+   . DenseMutableDcel s vertex directed undirected face+  -> MUV.MVector s Word32+  -> ST s (Maybe (AscendingLineArena s))+prepareAscendingLineArena dense radial+  | MUV.length radial < 3 = pure Nothing+  | otherwise = do+      first <- MUV.unsafeRead radial 0+      second <- MUV.unsafeRead radial 1+      third <- MUV.unsafeRead radial 2+      firstX <- coordinateX first+      firstY <- coordinateY first+      secondX <- coordinateX second+      secondY <- coordinateY second+      thirdX <- coordinateX third+      thirdY <- coordinateY third+      let !useX = firstX /= secondX+          pointCoordinate :: Double -> Double -> Double+          pointCoordinate x y = if useX then x else y+          coordinateAt :: Word32 -> ST s Double+          coordinateAt = if useX then coordinateX else coordinateY+          !pivotCoordinate = pointCoordinate firstX firstY+          !initialLowerCount =+            fromEnum (pointCoordinate secondX secondY < pivotCoordinate)+              + fromEnum (pointCoordinate thirdX thirdY < pivotCoordinate)+      if orient2dCoordinates firstX firstY secondX secondY thirdX thirdY /= EQ+        then pure Nothing+        else do+          (!collinear, !lowerCount) <-+            MUV.foldM'+              (\(!admitted, !count) raw ->+                 if admitted+                   then do+                     x <- coordinateX raw+                     y <- coordinateY raw+                     pure+                       ( orient2dCoordinates firstX firstY secondX secondY x y == EQ+                       , count + fromEnum (pointCoordinate x y < pivotCoordinate)+                       )+                   else pure (False, count)+              )+              (True, initialLowerCount)+              (MUV.unsafeSlice 3 (MUV.length radial - 3) radial)+          if collinear+            then Just <$> orderAlongLine coordinateAt pivotCoordinate lowerCount+            else pure Nothing+ where+  coordinateX :: Word32 -> ST s Double+  coordinateX = denseReadPointX dense . fromIntegral+  coordinateY :: Word32 -> ST s Double+  coordinateY = denseReadPointY dense . fromIntegral++  orderAlongLine+    :: (Word32 -> ST s Double)+    -> Double+    -> Int+    -> ST s (AscendingLineArena s)+  orderAlongLine coordinateAt pivotCoordinate lowerCount = do+    ascending <- MUV.new (MUV.length radial)+    _ <-+      MUV.foldM'+        (\(!lowerSeen, !upperSeen) raw -> do+           coordinate <- coordinateAt raw+           if coordinate < pivotCoordinate+             then do+               MUV.unsafeWrite ascending (lowerCount - lowerSeen - 1) raw+               pure (lowerSeen + 1, upperSeen)+             else do+               MUV.unsafeWrite ascending (lowerCount + upperSeen) raw+               pure (lowerSeen, upperSeen + 1)+        )+        (0, 0)+        radial+    pure (AscendingLineArena ascending)++-- | Glue a proved ascending line as one outer-face cycle. Every site and cell+-- is fresh, so the dense section initializes each topology slot exactly once+-- and commits connectivity once after the local pieces agree on their shared+-- endpoints.+connectAscendingLineArena+  :: DenseMutableDcel s vertex directed undirected face+  -> OperationState s+  -> AscendingLineArena s+  -> ST s (Either BuildError Int)+connectAscendingLineArena dense operation (AscendingLineArena vertices) = do+  let !vertexCount = MUV.length vertices+      !segmentCount = vertexCount - 1+      mutable = denseMutableOwner dense+  capacity <- ensureCellCapacity mutable segmentCount 0+  case capacity of+    Left obstruction -> pure (Left obstruction)+    Right () -> do+      edgeBase <- denseAddEdgeBlock dense segmentCount+      MUV.imapM_+        (\index fromRaw -> do+           toRaw <- MUV.unsafeRead vertices (index + 1)+           let !forward = edgeBase + 2 * index+               !backward = forward + 1+               !nextForward = if index == segmentCount - 1 then backward else forward + 2+               !nextBackward = if index == 0 then edgeBase else backward - 2+           denseWriteOrigin dense forward (unpackIndex fromRaw)+           denseWriteOrigin dense backward (unpackIndex toRaw)+           denseLinkEdges dense forward nextForward+           denseLinkEdges dense backward nextBackward+           denseWriteFace dense forward 0+           denseWriteFace dense backward 0+        )+        (MUV.unsafeSlice 0 segmentCount vertices)+      MUV.imapM_+        (\index raw ->+           denseMarkFreshConnected+             dense+             (unpackIndex raw)+             ( if index < segmentCount+                 then edgeBase + 2 * index+                 else edgeBase + 2 * segmentCount - 1+             )+        )+        vertices+      denseCommitFreshConnections dense vertexCount+      denseWriteFaceEdge dense 0 edgeBase+      addCounter operation CounterLineExtensions (max 0 (vertexCount - 2))+      pure (Right vertexCount)  collectLineChain :: MutableDcel s vertex directed undirected face -> OperationState s -> ST s (Either BuildError Int) collectLineChain mutable operation = do

@@ -25,7 +25,10 @@   , readPointY   , readPrevious   )-import Moonlight.Triangulation.Scalar (inCircleCoordinates, orient2dCoordinates)+import Moonlight.Triangulation.Internal.Predicates+  ( inCircleCoordinates+  , orient2dCoordinates+  )  -- | Whether the flip rule may fire on a quadrilateral @a b c d@ under a law, -- where @ab@ is the diagonal and @c@, @d@ the apexes opposite it. The

@@ -92,7 +92,7 @@   ) import Moonlight.Triangulation.Internal.Probe (KnownProbe) import Moonlight.Triangulation.Types (BuildError (..), Point (..))-import Moonlight.Triangulation.Scalar (orient2dCoordinates)+import Moonlight.Triangulation.Internal.Predicates (orient2dCoordinates)  -- | Proof that the remaining circle-sweep program fits the particular mutable -- arena carried here. The constructor is private: only 'reserveSweepCells' can

@@ -1,7 +1,6 @@ {-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE DataKinds #-}-{-# LANGUAGE NamedFieldPuns #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} {-# OPTIONS_GHC -O3 -fllvm -optlo-O3 -optlc-O3 #-}

@@ -1,11 +1,9 @@-{-# LANGUAGE BangPatterns #-}- module Moonlight.Triangulation.Internal.FaceProbe   ( BoundaryProbe (..)   , probeBoundary   ) where -import Moonlight.Triangulation.Math (onClosedSegment, orient2d)+import Moonlight.Triangulation.Internal.Predicates (onClosedSegment, orient2d) import Moonlight.Triangulation.Types (Point)  -- | Classification of one oriented boundary of a triangular face. A

@@ -61,7 +61,12 @@   ( simplifyBoundaryCycle   , traceOrientedBoundaryCircuits   )+import Moonlight.Triangulation.Internal.Predicates import Moonlight.Triangulation.Math+  ( mkQueryPoint+  , segmentDistanceSquaredWide+  , squaredDistanceWide+  ) import Moonlight.Triangulation.PointLocation import Moonlight.Triangulation.Types @@ -233,8 +238,8 @@   edges = edgesInShape triangulation metric   set = List.foldl' addEndpoints IntSet.empty edges   addEndpoints acc edge =-    let (VertexId from, VertexId to) = undirectedEndpoints triangulation edge-     in IntSet.insert (fromIntegral from) (IntSet.insert (fromIntegral to) acc)+    let (from, to) = undirectedEndpoints triangulation edge+     in IntSet.insert (vertexIdIndex from) (IntSet.insert (vertexIdIndex to) acc)   candidateVertices     | numVertices triangulation == 1 = [VertexId 0]     | otherwise = map (VertexId . fromIntegral) (IntSet.toAscList set)@@ -267,7 +272,7 @@   let (faces, accepted, rejected) = go initialStack initialVisited IntSet.empty IntSet.empty []    in (reverse faces, accepted, rejected)  where-  valid face@(FaceId value) = face /= outerFace && fromIntegral value < numFaces triangulation+  valid face = face /= outerFace && faceIdIndex face < numFaces triangulation   (initialStack, initialVisited) = List.foldl' enqueueStart ([], IntSet.empty) starts    enqueueStart state face@@ -285,8 +290,8 @@      in go stack' visited' accepted' rejected' (face : result)    expand (stack, visited, accepted, rejected) edge =-    let undirected@(UndirectedEdgeId raw) = asUndirected edge-        edgeIndex = fromIntegral raw+    let undirected = asUndirected edge+        edgeIndex = undirectedEdgeIdIndex undirected         adjacent = incidentFace triangulation (reverseEdge edge)         edgeAdmission           | IntSet.member edgeIndex accepted = (True, accepted, rejected)@@ -300,11 +305,11 @@             else (stack, visited)      in (stack', visited', accepted', rejected') -  enqueue face@(FaceId value) (stack, visited)+  enqueue face (stack, visited)     | IntSet.member index visited = (stack, visited)     | otherwise = (face : stack, IntSet.insert index visited)    where-    index = fromIntegral value+    index = faceIdIndex face  -- | Connected components of equally labelled bounded faces. Labels are -- evaluated once; the component carrier is the same 'IntSet' used by descent.@@ -394,9 +399,8 @@         )       else (queued, unvisited)    where-    adjacent@(FaceId adjacentRaw) =-      incidentFace triangulation (reverseEdge edge)-    adjacentIndex = fromIntegral adjacentRaw+    adjacent = incidentFace triangulation (reverseEdge edge)+    adjacentIndex = faceIdIndex adjacent {-# INLINE faceComponentsFromLabels #-}  -- | Descend every equally labelled bounded-face component through the one@@ -477,8 +481,7 @@             )         Nothing -> Right (FaceComponent indices)  where-  indices = IntSet.fromList (fmap faceIndex (Set.toAscList selectedFaces))-  faceIndex (FaceId raw) = fromIntegral raw+  indices = IntSet.fromList (fmap faceIdIndex (Set.toAscList selectedFaces))  componentBoundaryGraph   :: Triangulation mode vertex directed undirected face@@ -518,11 +521,11 @@               firstEdge               secondEdge           SimpleBoundaryGraph edges outgoing ->-            case IntMap.lookup vertexIndex outgoing of+            case IntMap.lookup residentIndex outgoing of               Nothing ->                 SimpleBoundaryGraph                   (insertEdge edges)-                  (IntMap.insert vertexIndex edge outgoing)+                  (IntMap.insert residentIndex edge outgoing)               Just previousEdge ->                 PinchedBoundaryGraph                   (insertEdge edges)@@ -530,12 +533,10 @@                   previousEdge                   edge    where-    FaceId adjacent = incidentFace triangulation (reverseEdge edge)-    adjacentIndex = fromIntegral adjacent-    vertex@(VertexId rawVertex) = origin triangulation edge-    vertexIndex = fromIntegral rawVertex-    DirectedEdgeId rawEdge = edge-    insertEdge = IntSet.insert (fromIntegral rawEdge)+    adjacentIndex = faceIdIndex (incidentFace triangulation (reverseEdge edge))+    vertex = origin triangulation edge+    residentIndex = vertexIdIndex vertex+    insertEdge = IntSet.insert (directedEdgeIdIndex edge)  regionBoundaryFromLoops   :: [BoundaryLoop]@@ -600,15 +601,13 @@         descend (loop : loops) remaining    traceCycle start current unvisited reversedVertices =-    let DirectedEdgeId rawCurrent = current-        remaining = IntSet.delete (fromIntegral rawCurrent) unvisited+    let remaining = IntSet.delete (directedEdgeIdIndex current) unvisited         accumulated = origin triangulation current : reversedVertices-        VertexId rawTarget = destination triangulation current-     in case IntMap.lookup (fromIntegral rawTarget) outgoingByVertex of+        target = destination triangulation current+     in case IntMap.lookup (vertexIdIndex target) outgoingByVertex of           Just successor             | successor == start -> Right (reverse accumulated, remaining)-            | let DirectedEdgeId rawSuccessor = successor-            , IntSet.member (fromIntegral rawSuccessor) remaining ->+            | IntSet.member (directedEdgeIdIndex successor) remaining ->                 traceCycle start successor remaining accumulated             | otherwise -> Left (BoundaryCycleDidNotClose start successor)           Nothing -> Left (BoundaryCycleDidNotClose start current)@@ -666,9 +665,7 @@ barrierDepthLayers triangulation isBarrier =   map (filter (/= outerFace)) (layers IntSet.empty [outerFace])  where-  known (FaceId value) = fromIntegral value < numFaces triangulation-  key :: FaceId -> Int-  key (FaceId value) = fromIntegral value+  known face = faceIdIndex face < numFaces triangulation    layers visited frontier = case flood visited [] frontier of     ([], _) -> []@@ -676,10 +673,10 @@    flood visited acc [] = (reverse acc, visited)   flood visited acc (face : rest)-    | not (known face) || IntSet.member (key face) visited = flood visited acc rest+    | not (known face) || IntSet.member (faceIdIndex face) visited = flood visited acc rest     | otherwise =         flood-          (IntSet.insert (key face) visited)+          (IntSet.insert (faceIdIndex face) visited)           (face : acc)           (neighbours (not . isBarrier) face <> rest) @@ -710,7 +707,14 @@     OnEdge edge -> filter (/= outerFace) [incidentFace triangulation edge, incidentFace triangulation (reverseEdge edge)]     OnVertex vertex ->       intSetToFaces-        (List.foldl' (\set edge -> let FaceId value = incidentFace triangulation edge in if value == 0 then set else IntSet.insert (fromIntegral value) set) IntSet.empty (vertexOutgoingEdges triangulation vertex))+        ( List.foldl'+            (\set edge ->+              let face = incidentFace triangulation edge+               in if face == outerFace then set else IntSet.insert (faceIdIndex face) set+            )+            IntSet.empty+            (vertexOutgoingEdges triangulation vertex)+        )     OutsideConvexHull _ ->       [ incidentFace triangulation (reverseEdge edge)       | edge <- faceDirectedEdges triangulation outerFace

@@ -4,9 +4,13 @@ -- its index complement, so orientation is a bit rather than a lookup. module Moonlight.Triangulation.Internal.HandleDefs   ( VertexId (..)+  , vertexIdIndex   , FaceId (..)+  , faceIdIndex   , DirectedEdgeId (..)+  , directedEdgeIdIndex   , UndirectedEdgeId (..)+  , undirectedEdgeIdIndex   , reverseEdge   , asUndirected   , normalizedDirected@@ -18,6 +22,7 @@ import Control.DeepSeq (NFData) import Data.Bits (shiftL, shiftR, xor, (.&.)) import Data.Word (Word32)+import Moonlight.Triangulation.Internal.PackedIndex (unpackIndex)  -- | Index of a vertex in the immutable DCEL. newtype VertexId = VertexId { unVertexId :: Word32 }@@ -38,6 +43,20 @@ newtype UndirectedEdgeId = UndirectedEdgeId { unUndirectedEdgeId :: Word32 }   deriving stock (Show)   deriving newtype (Eq, Ord, NFData)++-- | Project compact identifiers into resident vector indices.+vertexIdIndex :: VertexId -> Int+vertexIdIndex = unpackIndex . unVertexId+{-# INLINE vertexIdIndex #-}+faceIdIndex :: FaceId -> Int+faceIdIndex = unpackIndex . unFaceId+{-# INLINE faceIdIndex #-}+directedEdgeIdIndex :: DirectedEdgeId -> Int+directedEdgeIdIndex = unpackIndex . unDirectedEdgeId+{-# INLINE directedEdgeIdIndex #-}+undirectedEdgeIdIndex :: UndirectedEdgeId -> Int+undirectedEdgeIdIndex = unpackIndex . unUndirectedEdgeId+{-# INLINE undirectedEdgeIdIndex #-}  -- | Select the opposite orientation of the same undirected edge. reverseEdge :: DirectedEdgeId -> DirectedEdgeId

@@ -128,7 +128,7 @@   , VertexId (..)   ) import Moonlight.Triangulation.Internal.BoxedPaged-import Moonlight.Triangulation.Internal.PackedIndex (noIndex, packIndex)+import Moonlight.Triangulation.Internal.PackedIndex (noIndex, packIndex, unpackIndex) import Moonlight.Triangulation.Internal.Paged import Moonlight.Triangulation.Internal.PointIndex   ( MutablePointIndex@@ -1427,13 +1427,13 @@ {-# INLINE denseCommitFreshDefaultPoints #-}  denseReadOrigin :: DenseMutableDcel s vertex directed undirected face -> Int -> ST s Int-denseReadOrigin DenseMutableDcel{dmdHalfTopology} index = fromIntegral <$> readFlatMutable dmdHalfTopology (4 * index)+denseReadOrigin DenseMutableDcel{dmdHalfTopology} index = unpackIndex <$> readFlatMutable dmdHalfTopology (4 * index) denseReadNext :: DenseMutableDcel s vertex directed undirected face -> Int -> ST s Int-denseReadNext DenseMutableDcel{dmdHalfTopology} index = fromIntegral <$> readFlatMutable dmdHalfTopology (4 * index + 1)+denseReadNext DenseMutableDcel{dmdHalfTopology} index = unpackIndex <$> readFlatMutable dmdHalfTopology (4 * index + 1) denseReadPrevious :: DenseMutableDcel s vertex directed undirected face -> Int -> ST s Int-denseReadPrevious DenseMutableDcel{dmdHalfTopology} index = fromIntegral <$> readFlatMutable dmdHalfTopology (4 * index + 2)+denseReadPrevious DenseMutableDcel{dmdHalfTopology} index = unpackIndex <$> readFlatMutable dmdHalfTopology (4 * index + 2) denseReadFace :: DenseMutableDcel s vertex directed undirected face -> Int -> ST s Int-denseReadFace DenseMutableDcel{dmdHalfTopology} index = fromIntegral <$> readFlatMutable dmdHalfTopology (4 * index + 3)+denseReadFace DenseMutableDcel{dmdHalfTopology} index = unpackIndex <$> readFlatMutable dmdHalfTopology (4 * index + 3) {-# INLINE denseReadOrigin #-} {-# INLINE denseReadNext #-} {-# INLINE denseReadPrevious #-}@@ -1482,7 +1482,7 @@ denseReadFaceEdge :: DenseMutableDcel s vertex directed undirected face -> Int -> ST s Int denseReadFaceEdge DenseMutableDcel{dmdFaceEdge} index = do   value <- readFlatMutable dmdFaceEdge index-  pure (if value == noIndex then -1 else fromIntegral value)+  pure (if value == noIndex then -1 else unpackIndex value) {-# INLINE denseReadFaceEdge #-}  denseWriteFaceEdge :: DenseMutableDcel s vertex directed undirected face -> Int -> Int -> ST s ()@@ -1496,13 +1496,13 @@ {-# INLINE denseReadConstraint #-}  readOrigin :: MutableDcel s vertex directed undirected face -> Int -> ST s Int-readOrigin MutableDcel{mdHalfTopology} index = fromIntegral <$> readPaged mdHalfTopology (4 * index)+readOrigin MutableDcel{mdHalfTopology} index = unpackIndex <$> readPaged mdHalfTopology (4 * index) readNext :: MutableDcel s vertex directed undirected face -> Int -> ST s Int-readNext MutableDcel{mdHalfTopology} index = fromIntegral <$> readPaged mdHalfTopology (4 * index + 1)+readNext MutableDcel{mdHalfTopology} index = unpackIndex <$> readPaged mdHalfTopology (4 * index + 1) readPrevious :: MutableDcel s vertex directed undirected face -> Int -> ST s Int-readPrevious MutableDcel{mdHalfTopology} index = fromIntegral <$> readPaged mdHalfTopology (4 * index + 2)+readPrevious MutableDcel{mdHalfTopology} index = unpackIndex <$> readPaged mdHalfTopology (4 * index + 2) readFace :: MutableDcel s vertex directed undirected face -> Int -> ST s Int-readFace MutableDcel{mdHalfTopology} index = fromIntegral <$> readPaged mdHalfTopology (4 * index + 3)+readFace MutableDcel{mdHalfTopology} index = unpackIndex <$> readPaged mdHalfTopology (4 * index + 3) {-# INLINE readOrigin #-} {-# INLINE readNext #-} {-# INLINE readPrevious #-}@@ -1524,7 +1524,7 @@ readVertexOut :: MutableDcel s vertex directed undirected face -> Int -> ST s Int readVertexOut MutableDcel{mdVertexOut} index = do   value <- readPaged mdVertexOut index-  pure (if value == noIndex then -1 else fromIntegral value)+  pure (if value == noIndex then -1 else unpackIndex value) {-# INLINE readVertexOut #-}  writeVertexOut :: MutableDcel s vertex directed undirected face -> Int -> Int -> ST s ()@@ -1535,7 +1535,7 @@ readFaceEdge :: MutableDcel s vertex directed undirected face -> Int -> ST s Int readFaceEdge MutableDcel{mdFaceEdge} index = do   value <- readPaged mdFaceEdge index-  pure (if value == noIndex then -1 else fromIntegral value)+  pure (if value == noIndex then -1 else unpackIndex value) {-# INLINE readFaceEdge #-}  writeFaceEdge :: MutableDcel s vertex directed undirected face -> Int -> Int -> ST s ()

@@ -37,7 +37,7 @@   , readGrowable   , writeGrowable   )-import Moonlight.Triangulation.Internal.PackedIndex (packIndex)+import Moonlight.Triangulation.Internal.PackedIndex (packIndex, unpackIndex) import Moonlight.Triangulation.Internal.Types (BuildStats (..))  -- | The typed candidate section owned by one operation. The newtype prevents@@ -146,7 +146,7 @@ {-# INLINE writeScratch #-}  readScratch :: OperationState s -> Int -> ST s Int-readScratch OperationState{osScratchArena} index = fromIntegral <$> readGrowable osScratchArena index+readScratch OperationState{osScratchArena} index = unpackIndex <$> readGrowable osScratchArena index {-# INLINE readScratch #-}  addCounter :: OperationState s -> Counter -> Int -> ST s ()

@@ -0,0 +1,143 @@+{-# LANGUAGE BangPatterns #-}+{-# OPTIONS_GHC -O3 -fllvm -optlo-O3 -optlc-O3 #-}++-- | Exact predicates over points already admitted by a construction or query+-- boundary.  This module is package-private: public callers must carry+-- 'QueryPoint' evidence through "Moonlight.Triangulation.Math" instead of+-- presenting arbitrary IEEE-754 payloads to an exact classifier.+module Moonlight.Triangulation.Internal.Predicates+  ( orient2dCoordinates+  , inCircleCoordinates+  , orient2d+  , sideQuery+  , inCircle+  , onClosedSegment+  , segmentRelation+  , segmentsProperlyCross+  , segmentsIntersect+  , inDiametralCircle+  ) where++import Moonlight.Triangulation.Internal.Dyadic+  ( exactDiametralDot+  , exactInCircleDet+  , exactOrientSignDouble+  )+import Moonlight.Triangulation.Internal.SegmentRelation+  ( SegmentRelation (..)+  , segmentRelationWith+  )+import Moonlight.Triangulation.LineSideInfo (LineSideInfo, fromOrdering)+import Moonlight.Triangulation.Scalar+  ( isFinite+  , scalarCcwErrorBound+  , scalarInCircleErrorBound+  )+import Moonlight.Triangulation.Internal.Types (Point (..))++-- The coordinate kernels are package-private because their raw 'Double'+-- arguments carry no finite-range evidence. Every public exact classifier+-- descends through 'QueryPoint'; internal hot paths call these only after their+-- construction boundary has admitted the stored coordinates.+orient2dCoordinates+  :: Double -> Double -> Double -> Double -> Double -> Double+  -> Ordering+orient2dCoordinates ax ay bx by cx cy+  | abs determinant > scalarCcwErrorBound * determinantSum = compare determinant 0+  | otherwise = exactOrientSignDouble ax ay bx by cx cy+ where+  !left = (ax - cx) * (by - cy)+  !right = (ay - cy) * (bx - cx)+  !determinant = left - right+  !determinantSum = abs left + abs right+{-# INLINE orient2dCoordinates #-}++inCircleCoordinates+  :: Double -> Double -> Double -> Double+  -> Double -> Double -> Double -> Double+  -> Ordering+inCircleCoordinates ax ay bx by cx cy dx dy+  | abs determinant > scalarInCircleErrorBound * permanent = compare determinant 0+  | otherwise = compare (exactInCircleDet ax ay bx by cx cy dx dy) 0+ where+  !adx = ax - dx+  !ady = ay - dy+  !bdx = bx - dx+  !bdy = by - dy+  !cdx = cx - dx+  !cdy = cy - dy+  !abdet = adx * bdy - bdx * ady+  !bcdet = bdx * cdy - cdx * bdy+  !cadet = cdx * ady - adx * cdy+  !alift = adx * adx + ady * ady+  !blift = bdx * bdx + bdy * bdy+  !clift = cdx * cdx + cdy * cdy+  !determinant = alift * bcdet + blift * cadet + clift * abdet+  !permanent =+    (abs (bdx * cdy) + abs (cdx * bdy)) * alift+      + (abs (cdx * ady) + abs (adx * cdy)) * blift+      + (abs (adx * bdy) + abs (bdx * ady)) * clift+{-# INLINE inCircleCoordinates #-}++-- | Exact orientation of points whose coordinates have already been admitted.+orient2d :: Point -> Point -> Point -> Ordering+orient2d (Point ax ay) (Point bx by) (Point cx cy) =+  orient2dCoordinates ax ay bx by cx cy+{-# INLINE orient2d #-}++-- | Exact side of an admitted oriented line.+sideQuery :: Point -> Point -> Point -> LineSideInfo+sideQuery from to query = fromOrdering (orient2d from to query)+{-# INLINE sideQuery #-}++-- | Ordering of the admitted oriented incircle determinant.+inCircle :: Point -> Point -> Point -> Point -> Ordering+inCircle+  (Point ax ay)+  (Point bx by)+  (Point cx cy)+  (Point dx dy) =+    inCircleCoordinates ax ay bx by cx cy dx dy+{-# INLINE inCircle #-}++-- | Whether an admitted point lies on an admitted closed segment.+onClosedSegment :: Point -> Point -> Point -> Bool+onClosedSegment from@(Point fromX fromY) to@(Point toX toY) query@(Point queryX queryY) =+  orient2d from to query == EQ+    && queryX >= min fromX toX+    && queryX <= max fromX toX+    && queryY >= min fromY toY+    && queryY <= max fromY toY+{-# INLINE onClosedSegment #-}++-- | Exact relation between admitted closed segments.+segmentRelation :: Point -> Point -> Point -> Point -> SegmentRelation+segmentRelation = segmentRelationWith (==) compare orient2d onClosedSegment++-- | Whether two admitted closed segments share any point.+segmentsIntersect :: Point -> Point -> Point -> Point -> Bool+segmentsIntersect firstFrom firstTo secondFrom secondTo =+  segmentRelation firstFrom firstTo secondFrom secondTo /= SegmentsDisjoint++-- | Whether two admitted segments cross away from their endpoints.+segmentsProperlyCross :: Point -> Point -> Point -> Point -> Bool+segmentsProperlyCross firstFrom firstTo secondFrom secondTo =+  opposite (orient2d firstFrom firstTo secondFrom) (orient2d firstFrom firstTo secondTo)+    && opposite (orient2d secondFrom secondTo firstFrom) (orient2d secondFrom secondTo firstTo)+ where+  opposite LT GT = True+  opposite GT LT = True+  opposite _ _ = False++-- | Whether an admitted point lies in the closed diametral disk of an admitted+-- segment.+inDiametralCircle :: Point -> Point -> Point -> Bool+inDiametralCircle (Point ax ay) (Point bx by) (Point px py)+  | isFinite dot && abs dot > scalarCcwErrorBound * dotSum = dot < 0+  | otherwise = exactDiametralDot ax ay bx by px py <= 0+ where+  !left = (ax - px) * (bx - px)+  !right = (ay - py) * (by - py)+  !dot = left + right+  !dotSum = abs left + abs right+{-# INLINE inDiametralCircle #-}

@@ -54,6 +54,7 @@   , UndirectedEdgeId (..)   , VertexId (..)   )+import Moonlight.Triangulation.Internal.PackedIndex (unpackIndex, unpackOptionalIndex) import Moonlight.Triangulation.Internal.BoxedPaged (BoxedPaged, boxedDefaulted, boxedFromVector, boxedToVector) import Moonlight.Triangulation.Internal.Paged   ( Paged@@ -161,7 +162,7 @@       closesAtStart = do         lastEdge <- terminal         pure-          ( fromIntegral+          ( unpackIndex               (pagedUnsafeIndex (triHalfTopology triangulation) (4 * lastEdge + 1))               == start           )@@ -178,10 +179,7 @@     :: Triangulation mode vertex directed undirected face     -> Maybe Int   frontierStart mesh =-    let raw = pagedUnsafeIndex (triFaceEdge mesh) 0-     in if raw == (maxBound :: Word32)-          then Nothing-          else Just (fromIntegral raw)+    unpackOptionalIndex (pagedUnsafeIndex (triFaceEdge mesh) 0)    frontierStep     :: Int@@ -233,7 +231,7 @@     -> (Double, Double)   frontierPoint mesh rawEdge =     let rawVertex = pagedUnsafeIndex (triHalfTopology mesh) (4 * rawEdge)-        vertex = fromIntegral rawVertex+        vertex = unpackIndex rawVertex      in ( pagedUnsafeIndex (triPointX mesh) vertex         , pagedUnsafeIndex (triPointY mesh) vertex         )

@@ -1,5 +1,3 @@-{-# LANGUAGE BangPatterns #-}- -- | Discharge: the invariants the constructors guarantee, checkable on a value -- built by any route. module Moonlight.Triangulation.Internal.Validation@@ -22,8 +20,12 @@ import Moonlight.Triangulation.Dcel import Moonlight.Triangulation.Internal.HandleDefs import Moonlight.Triangulation.Handles.Iterators.FixedIterators (allFaces, directedEdges, undirectedEdges, vertices)-import Moonlight.Triangulation.Internal.PackedIndex (noIndex)+import Moonlight.Triangulation.Internal.PackedIndex (noIndex, unpackIndex)+import Moonlight.Triangulation.Internal.Predicates import Moonlight.Triangulation.Math+  ( squaredDistance+  , triangleArea+  ) import Moonlight.Triangulation.Internal.Representation import Moonlight.Triangulation.Internal.Types @@ -109,42 +111,42 @@     [ EdgeOriginOutOfRange (DirectedEdgeId (fromIntegral index)) (VertexId value) verticesCount     | index <- [0 .. halfCount - 1]     , let value = pagedUnsafeIndex (triHalfTopology triangulation) (4 * index)-    , fromIntegral value >= verticesCount+    , unpackIndex value >= verticesCount     ]       ++ [ EdgeNextOutOfRange (DirectedEdgeId (fromIntegral index)) (DirectedEdgeId value) halfCount          | index <- [0 .. halfCount - 1]          , let value = pagedUnsafeIndex (triHalfTopology triangulation) (4 * index + 1)-         , fromIntegral value >= halfCount+         , unpackIndex value >= halfCount          ]       ++ [ EdgePreviousOutOfRange (DirectedEdgeId (fromIntegral index)) (DirectedEdgeId value) halfCount          | index <- [0 .. halfCount - 1]          , let value = pagedUnsafeIndex (triHalfTopology triangulation) (4 * index + 2)-         , fromIntegral value >= halfCount+         , unpackIndex value >= halfCount          ]       ++ [ EdgeFaceOutOfRange (DirectedEdgeId (fromIntegral index)) (FaceId value) facesCount          | index <- [0 .. halfCount - 1]          , let value = pagedUnsafeIndex (triHalfTopology triangulation) (4 * index + 3)-         , fromIntegral value >= facesCount+         , unpackIndex value >= facesCount          ]       ++ [ VertexOutgoingOutOfRange (VertexId (fromIntegral index)) (DirectedEdgeId value) halfCount          | index <- [0 .. pagedLength (triVertexOut triangulation) - 1]          , let value = pagedUnsafeIndex (triVertexOut triangulation) index          , value /= noIndex-         , fromIntegral value >= halfCount+         , unpackIndex value >= halfCount          ]       ++ [ FaceAdjacentOutOfRange (FaceId (fromIntegral index)) (DirectedEdgeId value) halfCount          | index <- [0 .. facesCount - 1]          , let value = pagedUnsafeIndex (triFaceEdge triangulation) index          , value /= noIndex-         , fromIntegral value >= halfCount+         , unpackIndex value >= halfCount          ]    edgeViolations     | not halfArraysEqual || odd halfCount || not (null rangeViolations) = []     | otherwise = concatMap validateEdge (directedEdges triangulation) -  validateEdge edge@(DirectedEdgeId raw) =-    let index = fromIntegral raw+  validateEdge edge =+    let index = directedEdgeIdIndex edge         nextEdge = next triangulation edge         previousEdge = previous triangulation edge         twinEdge = reverseEdge edge@@ -167,7 +169,7 @@             else []      in if index < halfCount then local ++ innerCycle else [] -  validEdge (DirectedEdgeId value) = fromIntegral value < halfCount+  validEdge edge = directedEdgeIdIndex edge < halfCount    faceViolations     | not (null rangeViolations) = []@@ -257,7 +259,7 @@       , rawPair < numUndirectedEdges triangulation       , let edge = UndirectedEdgeId (fromIntegral rawPair)       , let (forward, backward) = directedPair edge-      , rawFace <- fmap faceRaw [incidentFace triangulation forward, incidentFace triangulation backward]+      , rawFace <- fmap faceIdIndex [incidentFace triangulation forward, incidentFace triangulation backward]       , rawFace > 0       ]   selectedFaces = IntSet.union admittedFaces collarFaces@@ -267,7 +269,7 @@       [ rawEdge       | rawFace <- IntSet.toAscList selectedFaces       , edge <- faceEdgesBounded (FaceId (fromIntegral rawFace))-      , rawEdge <- [edgeRaw edge, edgeRaw (reverseEdge edge)]+      , rawEdge <- [directedEdgeIdIndex edge, directedEdgeIdIndex (reverseEdge edge)]       ]   selectedInterfaceEdges =     IntSet.fromList@@ -284,10 +286,10 @@       | rawEdge <- IntSet.toAscList selectedEdges       , rawEdge >= 0       , rawEdge < halfCount-      , edgeRaw (reverseEdge (DirectedEdgeId (fromIntegral rawEdge))) < halfCount+      , directedEdgeIdIndex (reverseEdge (DirectedEdgeId (fromIntegral rawEdge))) < halfCount       , rawVertex <--          [ vertexRaw (origin triangulation (DirectedEdgeId (fromIntegral rawEdge)))-          , vertexRaw (destination triangulation (DirectedEdgeId (fromIntegral rawEdge)))+          [ vertexIdIndex (origin triangulation (DirectedEdgeId (fromIntegral rawEdge)))+          , vertexIdIndex (destination triangulation (DirectedEdgeId (fromIntegral rawEdge)))           ]       , rawVertex < verticesCount       ]@@ -299,7 +301,7 @@     , let edge = DirectedEdgeId (fromIntegral rawEdge)     , rawEdge < halfCount     , let vertex = origin triangulation edge-    , vertexRaw vertex >= verticesCount+    , vertexIdIndex vertex >= verticesCount     ]       ++ [ EdgeNextOutOfRange edge nextEdge halfCount          | rawEdge <- IntSet.toAscList selectedEdges@@ -307,7 +309,7 @@          , rawEdge < halfCount          , let edge = DirectedEdgeId (fromIntegral rawEdge)          , let nextEdge = next triangulation edge-         , edgeRaw nextEdge >= halfCount+         , directedEdgeIdIndex nextEdge >= halfCount          ]       ++ [ EdgePreviousOutOfRange edge previousEdge halfCount          | rawEdge <- IntSet.toAscList selectedEdges@@ -315,7 +317,7 @@          , rawEdge < halfCount          , let edge = DirectedEdgeId (fromIntegral rawEdge)          , let previousEdge = previous triangulation edge-         , edgeRaw previousEdge >= halfCount+         , directedEdgeIdIndex previousEdge >= halfCount          ]       ++ [ EdgeFaceOutOfRange edge face facesCount          | rawEdge <- IntSet.toAscList selectedEdges@@ -323,7 +325,7 @@          , rawEdge < halfCount          , let edge = DirectedEdgeId (fromIntegral rawEdge)          , let face = incidentFace triangulation edge-         , faceRaw face >= facesCount+         , faceIdIndex face >= facesCount          ]       ++ [ FaceAdjacentOutOfRange face edge halfCount          | rawFace <- IntSet.toAscList selectedFaces@@ -331,7 +333,7 @@          , rawFace < facesCount          , let face = FaceId (fromIntegral rawFace)          , Just edge <- [adjacentEdge triangulation face]-         , edgeRaw edge >= halfCount+         , directedEdgeIdIndex edge >= halfCount          ]    edgeViolations =@@ -344,9 +346,9 @@             nextEdge = next triangulation edge             previousEdge = previous triangulation edge             twinEdge = reverseEdge edge-            nextValid = edgeRaw nextEdge < halfCount-            previousValid = edgeRaw previousEdge < halfCount-            nextNextValid = nextValid && edgeRaw (next triangulation nextEdge) < halfCount+            nextValid = directedEdgeIdIndex nextEdge < halfCount+            previousValid = directedEdgeIdIndex previousEdge < halfCount+            nextNextValid = nextValid && directedEdgeIdIndex (next triangulation nextEdge) < halfCount             innerCycle =               if incidentFace triangulation edge /= outerFace && nextValid && previousValid && nextNextValid                 then [InnerFaceNotTriangularAtEdge edge | next triangulation (next triangulation nextEdge) /= edge]@@ -371,7 +373,7 @@        in case adjacentEdge triangulation face of           Nothing -> [FaceMissingAdjacentEdge face]           Just edge-            | edgeRaw edge >= halfCount -> []+            | directedEdgeIdIndex edge >= halfCount -> []             | otherwise ->                 let representedFace = incidentFace triangulation edge                     (faceEdges, faceVertices') = triangleEdgesAndVertices face@@ -387,7 +389,7 @@    triangleEdgesAndVertices face =     let edges = faceEdgesBounded face-        vertices' = fmap (vertexRaw . origin triangulation) edges+        vertices' = fmap (vertexIdIndex . origin triangulation) edges      in (edges, vertices')    triangleClosed edges =@@ -404,7 +406,7 @@             | verticesCount > 1 -> [ConnectedVertexMissingOutgoing vertex]             | otherwise -> []           Just edge-            | edgeRaw edge >= halfCount ->+            | directedEdgeIdIndex edge >= halfCount ->                 [VertexOutgoingOutOfRange vertex edge halfCount]             | otherwise ->                 [ VertexOutgoingOriginMismatch vertex edge actualOrigin@@ -451,20 +453,16 @@     adjacentEdges = case adjacentEdge triangulation face of       Nothing -> []       Just start-        | edgeRaw start >= halfCount -> []+        | directedEdgeIdIndex start >= halfCount -> []         | otherwise ->             let second = next triangulation start-             in if edgeRaw second >= halfCount+             in if directedEdgeIdIndex second >= halfCount                   then [start]                   else                     let third = next triangulation second-                     in if edgeRaw third >= halfCount+                     in if directedEdgeIdIndex third >= halfCount                           then [start, second]                           else [start, second, third]--  edgeRaw (DirectedEdgeId raw) = fromIntegral raw-  vertexRaw (VertexId raw) = fromIntegral raw-  faceRaw (FaceId raw) = fromIntegral raw  -- | Every edge whose circumcircle is not empty. validateDelaunay :: Triangulation mode vertex directed undirected face -> [InvariantViolation]

@@ -17,8 +17,12 @@ import Data.Void (Void, absurd) import Moonlight.Triangulation.Dcel import Moonlight.Triangulation.Internal.HandleDefs+import Moonlight.Triangulation.Internal.Predicates import Moonlight.Triangulation.Handles.Iterators.FixedIterators (undirectedEdges, vertices) import Moonlight.Triangulation.Math+  ( projectionFactor+  , squaredDistanceWide+  ) import Moonlight.Triangulation.PointLocation import Moonlight.Triangulation.Types @@ -38,16 +42,20 @@ lineIntersections triangulation queryFrom queryTo =   let !from = queryPointValue queryFrom       !to = queryPointValue queryTo-   in case firstIntersection triangulation queryFrom queryTo of-        Nothing -> []-        Just first -> walkIntersections triangulation from to first+   in if endpointNearerToWalkSeed triangulation to from+        then case firstIntersection triangulation queryTo queryFrom of+          Nothing -> []+          Just first -> walkIntersections AgainstQuery triangulation to from first+        else case firstIntersection triangulation queryFrom queryTo of+          Nothing -> []+          Just first -> walkIntersections AlongQuery triangulation from to first  -- | 'lineIntersections' between two existing vertices. lineIntersectionsBetweenVertices :: Triangulation mode vertex directed undirected face -> VertexId -> VertexId -> [Intersection] lineIntersectionsBetweenVertices triangulation fromVertex toVertex =   let from = vertexPoint triangulation fromVertex       to = vertexPoint triangulation toVertex-   in walkIntersections triangulation from to (VertexIntersection fromVertex)+   in walkIntersections AlongQuery triangulation from to (VertexIntersection fromVertex)  -- | The directed edges a crossing list cuts. conflictingEdges :: [Intersection] -> [DirectedEdgeId]@@ -127,16 +135,51 @@ -- where the answer is 'Void', the left branch is erased outright. {-# INLINE foldCorridor #-} -walkIntersections :: Triangulation mode vertex directed undirected face -> Point -> Point -> Intersection -> [Intersection]-walkIntersections triangulation from to first =+data QueryDirection = AlongQuery | AgainstQuery++walkIntersections+  :: QueryDirection+  -> Triangulation mode vertex directed undirected face+  -> Point+  -> Point+  -> Intersection+  -> [Intersection]+walkIntersections direction triangulation from to first =   case foldCorridor triangulation from to first collect [] of-    Nothing -> exactIntersectionScan triangulation from to+    Nothing -> orientExact (exactIntersectionScan triangulation from to)     Just (Left impossible) -> absurd impossible-    Just (Right events) -> reverse events+    Just (Right events) -> orientCollected events  where   collect :: [Intersection] -> Intersection -> Either Void [Intersection]   collect accumulated event = Right (event : accumulated)+  orientExact = case direction of+    AlongQuery -> id+    AgainstQuery -> map reverseIntersection . reverse+  orientCollected = case direction of+    AlongQuery -> reverse+    AgainstQuery -> map reverseIntersection +reverseIntersection :: Intersection -> Intersection+reverseIntersection event = case event of+  VertexIntersection vertex -> VertexIntersection vertex+  EdgeIntersection edge -> EdgeIntersection (reverseEdge edge)+  EdgeOverlap edge -> EdgeOverlap (reverseEdge edge)++-- The unhinted point locator starts at face one. Choosing the segment endpoint+-- decisively nearer one of that face's vertices shortens the same exact walk+-- without making ordinary interior chords pay to reverse similar paths.+endpointNearerToWalkSeed+  :: Triangulation mode vertex directed undirected face+  -> Point+  -> Point+  -> Bool+endpointNearerToWalkSeed triangulation candidate incumbent =+  case innerFaceVertices triangulation (FaceId 1) of+    Just (anchor, _, _) ->+      let !seed = vertexPoint triangulation anchor+       in 4 * squaredDistanceWide seed candidate < squaredDistanceWide seed incumbent+    Nothing -> False+ nextIntersection :: Triangulation mode vertex directed undirected face -> Point -> Point -> Intersection -> Maybe Intersection nextIntersection triangulation lineFrom lineTo current = case current of   EdgeIntersection edge -> case traceDirectionOutOfEdge triangulation edge lineFrom lineTo of@@ -203,11 +246,10 @@   -- event. A ring that outran its budget, an absent locator edge, and   -- endpoints that leave the parameter comparison without a total order all   -- keep the exact scan.-  firstFromOutside (Just edge)-    | finiteEndpoint lineFrom && finiteEndpoint lineTo =-        case ringEntryEvent triangulation lineFrom lineTo edge of-          Just entry -> eventValue <$> entry-          Nothing -> firstFromScan+  firstFromOutside (Just edge) =+    case hullEntryIntersection triangulation lineFrom lineTo edge of+      Just entry -> entry+      Nothing -> firstFromScan   firstFromOutside _ = firstFromScan    firstFromScan = case exactIntersectionScan triangulation lineFrom lineTo of@@ -370,35 +412,50 @@   high = min 1 (max projectedA projectedB)   oriented = if orient2d a b lineTo == LT then reverseEdge directed else directed --- | The earliest event carried by the outer-face ring the given edge sits on,--- walked from that edge. 'Nothing' reports a ring that outran its budget.-ringEntryEvent+-- | Restrict an outside segment to the hull edge supplied by point location,+-- then follow only the side on which its entry can lie. The outer 'Maybe'+-- witnesses a hull walk that closed within its topology budget; the inner one+-- is the geometric answer. An exhausted or impossible local section descends+-- to the exact scan owned by 'firstIntersection'.+hullEntryIntersection   :: Triangulation mode vertex directed undirected face   -> Point   -> Point   -> DirectedEdgeId-  -> Maybe (Maybe (Event))-ringEntryEvent triangulation lineFrom lineTo entry =-  go (numDirectedEdges triangulation + 1) entry False Nothing+  -> Maybe (Maybe Intersection)+hullEntryIntersection triangulation lineFrom lineTo entry+  | orient2d entryFrom entryTo lineFrom == EQ =+      Just+        ( Just+            ( VertexIntersection+                ( if squaredDistanceWide entryTo lineFrom < squaredDistanceWide entryFrom lineFrom+                    then destination triangulation entry+                    else origin triangulation entry+                )+            )+        )+  | otherwise = go (numDirectedEdges triangulation + 1) entry  where-  go !remaining !edge !departed !earliest+  entryFrom = vertexPoint triangulation (origin triangulation entry)+  entryTo = vertexPoint triangulation (destination triangulation entry)++  go !remaining !edge     | remaining <= 0 = Nothing-    | departed && edge == entry = Just earliest     | otherwise =-        let !stepped =-              keepEarliest (edgeEvent triangulation lineFrom lineTo (asUndirected edge)) $-                keepEarliest (vertexEvent triangulation lineFrom lineTo (origin triangulation edge)) earliest-         in go (remaining - 1) (next triangulation edge) True stepped--keepEarliest :: Maybe (Event) -> Maybe (Event) -> Maybe (Event)-keepEarliest Nothing held = held-keepEarliest candidate Nothing = candidate-keepEarliest candidate@(Just proposed) held@(Just incumbent)-  | compareEvent proposed incumbent == LT = candidate-  | otherwise = held--finiteEndpoint :: Point -> Bool-finiteEndpoint point = isFinite (pointX point) && isFinite (pointY point)+        let !edgeFrom = vertexPoint triangulation (origin triangulation edge)+            !edgeTo = vertexPoint triangulation (destination triangulation edge)+            !fromSide = orient2d lineFrom lineTo edgeFrom+            !toSide = orient2d lineFrom lineTo edgeTo+         in if orient2d edgeFrom edgeTo lineTo == GT+              then Just Nothing+              else case (fromSide == GT, toSide /= LT) of+                (True, True) -> go (remaining - 1) (previous triangulation edge)+                (False, False) -> go (remaining - 1) (next triangulation edge)+                (False, True)+                  | toSide == EQ -> Just (Just (VertexIntersection (destination triangulation edge)))+                  | fromSide == EQ -> Just (Just (VertexIntersection (origin triangulation edge)))+                  | otherwise -> Just (Just (EdgeIntersection (reverseEdge edge)))+                (True, False) -> Nothing  segmentIntersectionParameter :: Point -> Point -> Point -> Point -> Double segmentIntersectionParameter (Point ax ay) (Point bx by) (Point cx cy) (Point dx dy)

@@ -18,6 +18,8 @@   , segmentDistanceSquared   , segmentDistanceSquaredWide   , distance+  , negateVector+  , normalizeVector   , midpoint   , centroid   , triangleArea@@ -39,22 +41,19 @@  import Moonlight.Triangulation.Internal.Dyadic   ( exactBarycentricDeterminants-  , exactDiametralDot   , integerRatioToDouble   )+import Moonlight.Triangulation.Internal.Predicates qualified as Admitted import Moonlight.Triangulation.Internal.SegmentRelation   ( SegmentRelation (..)   , allSegmentRelations-  , segmentRelationWith   )-import Moonlight.Triangulation.LineSideInfo (LineSideInfo, fromOrdering)+import Moonlight.Triangulation.LineSideInfo (LineSideInfo) import Moonlight.Triangulation.Scalar   ( canonicalScalarZero-  , inCircleCoordinates   , isFinite   , maximumAllowedCoordinate   , minimumAllowedCoordinate-  , orient2dCoordinates   , scalarCcwErrorBound   ) import Moonlight.Triangulation.Internal.Types@@ -65,42 +64,50 @@   , QueryPoint (..)   ) --- | Exact relation between two closed segments.+-- | Exact relation between two admitted closed segments. segmentRelation-  :: Point-  -> Point-  -> Point-  -> Point+  :: QueryPoint+  -> QueryPoint+  -> QueryPoint+  -> QueryPoint   -> SegmentRelation-segmentRelation a b c d =-  segmentRelationWith (==) compare orient2d onClosedSegment a b c d+segmentRelation firstFrom firstTo secondFrom secondTo =+  Admitted.segmentRelation+    (queryPointValue firstFrom)+    (queryPointValue firstTo)+    (queryPointValue secondFrom)+    (queryPointValue secondTo) --- | Whether two closed segments share any point.+-- | Whether two admitted closed segments share any point. segmentsIntersect-  :: Point-  -> Point-  -> Point-  -> Point+  :: QueryPoint+  -> QueryPoint+  -> QueryPoint+  -> QueryPoint   -> Bool-segmentsIntersect a b c d = segmentRelation a b c d /= SegmentsDisjoint+segmentsIntersect firstFrom firstTo secondFrom secondTo =+  Admitted.segmentsIntersect+    (queryPointValue firstFrom)+    (queryPointValue firstTo)+    (queryPointValue secondFrom)+    (queryPointValue secondTo)  -- | The proper-crossing section of 'segmentRelation'. Consumers which reject -- only that constructor need not compute the collinear and endpoint-touch -- distinctions required by the complete ADT after either side already proves -- separation. segmentsProperlyCross-  :: Point-  -> Point-  -> Point-  -> Point+  :: QueryPoint+  -> QueryPoint+  -> QueryPoint+  -> QueryPoint   -> Bool-segmentsProperlyCross a b c d =-  opposite (orient2d a b c) (orient2d a b d)-    && opposite (orient2d c d a) (orient2d c d b)- where-  opposite LT GT = True-  opposite GT LT = True-  opposite _ _ = False+segmentsProperlyCross firstFrom firstTo secondFrom secondTo =+  Admitted.segmentsProperlyCross+    (queryPointValue firstFrom)+    (queryPointValue firstTo)+    (queryPointValue secondFrom)+    (queryPointValue secondTo)  -- | Classify a coordinate outside the exact-predicate input domain. validateCoordinate :: Double -> Maybe CoordinateError@@ -156,15 +163,22 @@   (ax - cx) * (by - cy) - (ay - cy) * (bx - cx) {-# INLINE orientDetApprox #-} --- | Exact orientation ordering of three points.-orient2d :: Point -> Point -> Point -> Ordering-orient2d (Point ax ay) (Point bx by) (Point cx cy) =-  orient2dCoordinates ax ay bx by cx cy+-- | Exact orientation ordering of three admitted points.+orient2d :: QueryPoint -> QueryPoint -> QueryPoint -> Ordering+orient2d first second third =+  Admitted.orient2d+    (queryPointValue first)+    (queryPointValue second)+    (queryPointValue third) {-# INLINE orient2d #-} --- | Exact side of an oriented line.-sideQuery :: Point -> Point -> Point -> LineSideInfo-sideQuery a b point = fromOrdering (orient2d a b point)+-- | Exact side of an admitted oriented line.+sideQuery :: QueryPoint -> QueryPoint -> QueryPoint -> LineSideInfo+sideQuery from to query =+  Admitted.sideQuery+    (queryPointValue from)+    (queryPointValue to)+    (queryPointValue query) {-# INLINE sideQuery #-}  -- | Fast approximate oriented in-circle determinant.@@ -191,26 +205,26 @@   !clift = cdx * cdx + cdy * cdy {-# INLINE inCircleDetApprox #-} --- | Ordering of the oriented incircle determinant. For a counter-clockwise--- triangle, 'GT' means the fourth point lies strictly inside its circumcircle.-inCircle-  :: Point -> Point -> Point -> Point -> Ordering+-- | Ordering of the admitted oriented incircle determinant. For a+-- counter-clockwise triangle, 'GT' means the fourth point lies strictly inside+-- its circumcircle. inCircle-  (Point ax ay)-  (Point bx by)-  (Point cx cy)-  (Point dx dy) =-    inCircleCoordinates ax ay bx by cx cy dx dy+  :: QueryPoint -> QueryPoint -> QueryPoint -> QueryPoint -> Ordering+inCircle first second third query =+  Admitted.inCircle+    (queryPointValue first)+    (queryPointValue second)+    (queryPointValue third)+    (queryPointValue query) {-# INLINE inCircle #-} --- | Whether a point lies on a closed segment.-onClosedSegment :: Point -> Point -> Point -> Bool-onClosedSegment a@(Point ax ay) b@(Point bx by) query@(Point qx qy) =-  orient2d a b query == EQ-    && qx >= min ax bx-    && qx <= max ax bx-    && qy >= min ay by-    && qy <= max ay by+-- | Whether an admitted point lies on an admitted closed segment.+onClosedSegment :: QueryPoint -> QueryPoint -> QueryPoint -> Bool+onClosedSegment from to query =+  Admitted.onClosedSegment+    (queryPointValue from)+    (queryPointValue to)+    (queryPointValue query) {-# INLINE onClosedSegment #-}  -- | Squared Euclidean distance.@@ -253,6 +267,24 @@ distance left right = sqrt (squaredDistance left right) {-# INLINE distance #-} +-- | Negate a point interpreted as a vector.+negateVector :: Point -> Point+negateVector (Point x y) = Point (-x) (-y)+{-# INLINE negateVector #-}++-- | Normalize a point interpreted as a vector without overflowing its length.+normalizeVector :: Point -> Point+normalizeVector (Point x y)+  | scale == 0 = Point 0 0+  | otherwise =+      let !scaledX = x / scale+          !scaledY = y / scale+          !length' = sqrt (scaledX * scaledX + scaledY * scaledY)+       in Point (scaledX / length') (scaledY / length')+ where+  !scale = max (abs x) (abs y)+{-# INLINE normalizeVector #-}+ -- | Midpoint of two points. midpoint :: Point -> Point -> Point midpoint (Point ax ay) (Point bx by) = Point (0.5 * ax + 0.5 * bx) (0.5 * ay + 0.5 * by)@@ -385,8 +417,8 @@         !right = (py - sy) * (rx - sx)      in (left - right, abs left + abs right) --- | Whether a point lies in the closed diametral disk of a segment, decided by--- the sign of @(a-p)·(b-p)@.+-- | Whether an admitted point lies in the closed diametral disk of an admitted+-- segment, decided by the sign of @(a-p)·(b-p)@. -- -- The approximation is two coordinate differences and one product per term -- combined additively, which is the arithmetic shape 'orient2d' is analysed@@ -394,15 +426,12 @@ -- addition does not depend on the sign of that combination. The orientation -- coefficient therefore transfers unchanged, and the exact dot product remains -- the oracle for the uncertain band.-inDiametralCircle :: Point -> Point -> Point -> Bool-inDiametralCircle (Point ax ay) (Point bx by) (Point px py)-  | isFinite dot && abs dot > scalarCcwErrorBound * dotSum = dot < 0-  | otherwise = exactDiametralDot ax ay bx by px py <= 0- where-  !left = (ax - px) * (bx - px)-  !right = (ay - py) * (by - py)-  !dot = left + right-  !dotSum = abs left + abs right+inDiametralCircle :: QueryPoint -> QueryPoint -> QueryPoint -> Bool+inDiametralCircle from to query =+  Admitted.inDiametralCircle+    (queryPointValue from)+    (queryPointValue to)+    (queryPointValue query) {-# INLINE inDiametralCircle #-}  -- | Projection parameter of a point onto an oriented segment line.

@@ -12,7 +12,7 @@ import Moonlight.Triangulation.Internal.HandleDefs import Moonlight.Triangulation.Handles.Iterators.FixedIterators import Moonlight.Triangulation.Internal.FaceProbe-import Moonlight.Triangulation.Math+import Moonlight.Triangulation.Internal.Predicates import Moonlight.Triangulation.Types  -- | Locate an admitted point without an initial topology hint.@@ -118,8 +118,8 @@ chooseStart :: Triangulation mode vertex directed undirected face -> Maybe LocationHint -> FaceId chooseStart triangulation hint =   case hint of-    Just (FaceHint face@(FaceId index))-      | index > 0 && fromIntegral index < numFaces triangulation -> face+    Just (FaceHint face)+      | face /= outerFace && faceIdIndex face < numFaces triangulation -> face     Just (VertexHint vertex@(VertexId index))       | toInteger index < toInteger (numVertices triangulation) ->           case [face | edge <- vertexOutgoingEdges triangulation vertex, let face = incidentFace triangulation edge, face /= outerFace] of

@@ -0,0 +1,226 @@+{-# LANGUAGE DeriveAnyClass #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DerivingStrategies #-}+{-# LANGUAGE TypeApplications #-}++-- | Canonical labelled simplices and downward-closed planar complexes.+-- Labels, rather than resident DCEL handles, own identity across independently+-- constructed meshes.+module Moonlight.Triangulation.Simplex+  ( SimplexDimension (..)+  , PlanarSimplex+  , PlanarSimplexError (..)+  , planarVertex+  , planarEdge+  , planarFace+  , planarSimplexDimension+  , planarSimplexVertices+  , planarSimplexBoundary+  , PlanarComplex+  , PlanarComplexError (..)+  , planarComplex+  , planarComplexCells+  , planarComplexBasisAt+  , unionPlanarComplex+  , intersectPlanarComplex+  )+where++import Control.DeepSeq (NFData)+import Data.List qualified as List+import Data.List.NonEmpty (NonEmpty (..))+import Data.Set (Set)+import Data.Set qualified as Set+import Data.Vector (Vector)+import Data.Vector qualified as Vector+import GHC.Generics (Generic)++-- | The closed dimension vocabulary for planar simplicial complexes.+data SimplexDimension+  = SimplexDimension0+  | SimplexDimension1+  | SimplexDimension2+  deriving stock (Eq, Ord, Show, Enum, Bounded, Generic)+  deriving anyclass (NFData)++-- | One canonical simplex. Edge and face labels are strictly ascending.+data PlanarSimplex label+  = PlanarVertex !label+  | PlanarEdge !label !label+  | PlanarFace !label !label !label+  deriving stock (Eq, Ord, Show, Generic)+  deriving anyclass (NFData)++-- | A repeated vertex would lower the requested simplex's dimension.+data PlanarSimplexError label+  = PlanarSimplexRepeatedVertex !label+  deriving stock (Eq, Ord, Show, Generic)+  deriving anyclass (NFData)++planarVertex :: label -> PlanarSimplex label+planarVertex = PlanarVertex++planarEdge+  :: Ord label+  => label+  -> label+  -> Either (PlanarSimplexError label) (PlanarSimplex label)+planarEdge firstLabel secondLabel =+  case compare firstLabel secondLabel of+    LT -> Right (PlanarEdge firstLabel secondLabel)+    EQ -> Left (PlanarSimplexRepeatedVertex firstLabel)+    GT -> Right (PlanarEdge secondLabel firstLabel)++planarFace+  :: Ord label+  => label+  -> label+  -> label+  -> Either (PlanarSimplexError label) (PlanarSimplex label)+planarFace firstLabel secondLabel thirdLabel =+  let (firstLow, firstHigh) = ascendingPair firstLabel secondLabel+      (secondLow, finalHigh) = ascendingPair firstHigh thirdLabel+      (finalLow, finalMiddle) = ascendingPair firstLow secondLow+   in if finalLow == finalMiddle+        then Left (PlanarSimplexRepeatedVertex finalLow)+        else+          if finalMiddle == finalHigh+            then Left (PlanarSimplexRepeatedVertex finalMiddle)+            else Right (PlanarFace finalLow finalMiddle finalHigh)++planarSimplexDimension :: PlanarSimplex label -> SimplexDimension+planarSimplexDimension simplex =+  case simplex of+    PlanarVertex _ -> SimplexDimension0+    PlanarEdge _ _ -> SimplexDimension1+    PlanarFace _ _ _ -> SimplexDimension2++planarSimplexVertices :: PlanarSimplex label -> NonEmpty label+planarSimplexVertices simplex =+  case simplex of+    PlanarVertex label -> label :| []+    PlanarEdge firstLabel secondLabel -> firstLabel :| [secondLabel]+    PlanarFace firstLabel secondLabel thirdLabel ->+      firstLabel :| [secondLabel, thirdLabel]++-- | The canonical oriented simplicial boundary.+planarSimplexBoundary+  :: Num coefficient+  => PlanarSimplex label+  -> [(coefficient, PlanarSimplex label)]+planarSimplexBoundary simplex =+  case simplex of+    PlanarVertex _ -> []+    PlanarEdge firstLabel secondLabel ->+      [ (-1, PlanarVertex firstLabel)+      , (1, PlanarVertex secondLabel)+      ]+    PlanarFace firstLabel secondLabel thirdLabel ->+      [ (1, PlanarEdge secondLabel thirdLabel)+      , (-1, PlanarEdge firstLabel thirdLabel)+      , (1, PlanarEdge firstLabel secondLabel)+      ]++-- | A requested cell whose boundary is absent from the same complex.+data PlanarComplexError label+  = PlanarComplexBoundaryMissing+      !(PlanarSimplex label)+      !(PlanarSimplex label)+  deriving stock (Eq, Ord, Show, Generic)+  deriving anyclass (NFData)++-- | An admitted downward-closed complex and its canonical dense bases.+data PlanarComplex label = PlanarComplex+  { storedPlanarComplexCells :: !(Set (PlanarSimplex label))+  , storedPlanarVertexBasis :: !(Vector (PlanarSimplex label))+  , storedPlanarEdgeBasis :: !(Vector (PlanarSimplex label))+  , storedPlanarFaceBasis :: !(Vector (PlanarSimplex label))+  }+  deriving stock (Show, Generic)+  deriving anyclass (NFData)++instance Eq label => Eq (PlanarComplex label) where+  left == right =+    storedPlanarComplexCells left == storedPlanarComplexCells right++planarComplex+  :: Ord label+  => Set (PlanarSimplex label)+  -> Either (PlanarComplexError label) (PlanarComplex label)+planarComplex cells =+  case List.find (\(_, boundary) -> Set.notMember boundary cells) boundaryPairs of+    Just (simplex, boundary) ->+      Left (PlanarComplexBoundaryMissing simplex boundary)+    Nothing -> Right (admittedPlanarComplex cells)+ where+  boundaryPairs =+    [ (simplex, boundary)+    | simplex <- Set.toAscList cells+    , (_, boundary) <- planarSimplexBoundary @Int simplex+    ]++planarComplexCells :: PlanarComplex label -> Set (PlanarSimplex label)+planarComplexCells = storedPlanarComplexCells++planarComplexBasisAt+  :: SimplexDimension+  -> PlanarComplex label+  -> Vector (PlanarSimplex label)+planarComplexBasisAt dimension complexValue =+  case dimension of+    SimplexDimension0 -> storedPlanarVertexBasis complexValue+    SimplexDimension1 -> storedPlanarEdgeBasis complexValue+    SimplexDimension2 -> storedPlanarFaceBasis complexValue++-- | Union preserves downward closure, so no second validation is needed.+unionPlanarComplex+  :: Ord label+  => PlanarComplex label+  -> PlanarComplex label+  -> PlanarComplex label+unionPlanarComplex left right =+  admittedPlanarComplex+    (Set.union (planarComplexCells left) (planarComplexCells right))++-- | Intersection also preserves downward closure.+intersectPlanarComplex+  :: Ord label+  => PlanarComplex label+  -> PlanarComplex label+  -> PlanarComplex label+intersectPlanarComplex left right =+  admittedPlanarComplex+    (Set.intersection (planarComplexCells left) (planarComplexCells right))++admittedPlanarComplex :: Set (PlanarSimplex label) -> PlanarComplex label+admittedPlanarComplex cells =+  let (vertices, edges, faces) =+        foldr partitionSimplex ([], [], []) (Set.toAscList cells)+   in PlanarComplex+        { storedPlanarComplexCells = cells+        , storedPlanarVertexBasis = Vector.fromList vertices+        , storedPlanarEdgeBasis = Vector.fromList edges+        , storedPlanarFaceBasis = Vector.fromList faces+        }++partitionSimplex+  :: PlanarSimplex label+  -> ( [PlanarSimplex label]+     , [PlanarSimplex label]+     , [PlanarSimplex label]+     )+  -> ( [PlanarSimplex label]+     , [PlanarSimplex label]+     , [PlanarSimplex label]+     )+partitionSimplex simplex (vertices, edges, faces) =+  case simplex of+    PlanarVertex _ -> (simplex : vertices, edges, faces)+    PlanarEdge _ _ -> (vertices, simplex : edges, faces)+    PlanarFace _ _ _ -> (vertices, edges, simplex : faces)++ascendingPair :: Ord value => value -> value -> (value, value)+ascendingPair firstValue secondValue =+  if firstValue <= secondValue+    then (firstValue, secondValue)+    else (secondValue, firstValue)

@@ -11,6 +11,15 @@   , eulerCharacteristicValue   , ExactArea   , exactAreaValue+  , ExactPlanarMoments+  , exactPlanarMeasure+  , exactPlanarFirstX+  , exactPlanarFirstY+  , exactPlanarSecondXX+  , exactPlanarSecondXY+  , exactPlanarSecondYY+  , scaleExactPlanarMoments+  , polygonComponentMoments   , polygonComponentArea   , ExactLengthTerm   , lengthCoefficient@@ -68,6 +77,8 @@   , UndirectedEdgeId (..)   , VertexId (..)   , directedPair+  , faceIdIndex+  , vertexIdIndex   ) import Moonlight.Triangulation.Internal.CellSet   ( ExactCellSet (..)@@ -84,6 +95,7 @@   , exactRationalDenominator   , exactRationalFromDyadic   , exactRationalFromFiniteDouble+  , exactRationalFromNormalizedRatio   , exactRationalIsZero   , exactRationalNumerator   , exactSignum@@ -128,12 +140,126 @@ exactAreaValue :: ExactArea -> ExactRational exactAreaValue (ExactArea value) = value +-- | Exact moments of a bounded planar measure through total degree two.+data ExactPlanarMoments = ExactPlanarMoments+  { exactPlanarMeasure :: !ExactRational+    -- ^ Total signed measure.+  , exactPlanarFirstX :: !ExactRational+    -- ^ Raw first moment integral of @x@.+  , exactPlanarFirstY :: !ExactRational+    -- ^ Raw first moment integral of @y@.+  , exactPlanarSecondXX :: !ExactRational+    -- ^ Raw second moment integral of @x^2@.+  , exactPlanarSecondXY :: !ExactRational+    -- ^ Raw mixed moment integral of @x*y@.+  , exactPlanarSecondYY :: !ExactRational+    -- ^ Raw second moment integral of @y^2@.+  }+  deriving stock (Eq, Ord, Show, Generic)+  deriving anyclass (NFData)++instance Semigroup ExactPlanarMoments where+  left <> right =+    ExactPlanarMoments+      { exactPlanarMeasure = exactPlanarMeasure left + exactPlanarMeasure right+      , exactPlanarFirstX = exactPlanarFirstX left + exactPlanarFirstX right+      , exactPlanarFirstY = exactPlanarFirstY left + exactPlanarFirstY right+      , exactPlanarSecondXX = exactPlanarSecondXX left + exactPlanarSecondXX right+      , exactPlanarSecondXY = exactPlanarSecondXY left + exactPlanarSecondXY right+      , exactPlanarSecondYY = exactPlanarSecondYY left + exactPlanarSecondYY right+      }++instance Monoid ExactPlanarMoments where+  mempty = ExactPlanarMoments 0 0 0 0 0 0++-- | Scale every moment by one exact coefficient.+scaleExactPlanarMoments+  :: ExactRational+  -> ExactPlanarMoments+  -> ExactPlanarMoments+scaleExactPlanarMoments scalar moments =+  ExactPlanarMoments+    { exactPlanarMeasure = scalar * exactPlanarMeasure moments+    , exactPlanarFirstX = scalar * exactPlanarFirstX moments+    , exactPlanarFirstY = scalar * exactPlanarFirstY moments+    , exactPlanarSecondXX = scalar * exactPlanarSecondXX moments+    , exactPlanarSecondXY = scalar * exactPlanarSecondXY moments+    , exactPlanarSecondYY = scalar * exactPlanarSecondYY moments+    }++-- | Exact area and raw moments of one admitted component. The six unscaled+-- boundary sums are accumulated together and normalized once.+polygonComponentMoments :: PolygonComponent -> ExactPlanarMoments+polygonComponentMoments component =+  let raw =+        foldMap+          loopMomentSums+          (polygonOuterLoop component : polygonHoleLoops component)+   in ExactPlanarMoments+        { exactPlanarMeasure = oneHalf * rawDoubleArea raw+        , exactPlanarFirstX = oneSixth * rawFirstX raw+        , exactPlanarFirstY = oneSixth * rawFirstY raw+        , exactPlanarSecondXX = oneTwelfth * rawSecondXX raw+        , exactPlanarSecondXY = oneTwentyFourth * rawSecondXY raw+        , exactPlanarSecondYY = oneTwelfth * rawSecondYY raw+        }+ -- | Exact unsigned area of one already-admitted polygon component.  Winding -- and hole containment were discharged by 'polygonComponent', so this -- observation performs no second geometric validation. polygonComponentArea :: PolygonComponent -> ExactArea polygonComponentArea = ExactArea . (oneHalf *) . componentDoubleArea +data RawPlanarMomentSums = RawPlanarMomentSums+  { rawDoubleArea :: !ExactRational+  , rawFirstX :: !ExactRational+  , rawFirstY :: !ExactRational+  , rawSecondXX :: !ExactRational+  , rawSecondXY :: !ExactRational+  , rawSecondYY :: !ExactRational+  }++instance Semigroup RawPlanarMomentSums where+  left <> right =+    RawPlanarMomentSums+      { rawDoubleArea = rawDoubleArea left + rawDoubleArea right+      , rawFirstX = rawFirstX left + rawFirstX right+      , rawFirstY = rawFirstY left + rawFirstY right+      , rawSecondXX = rawSecondXX left + rawSecondXX right+      , rawSecondXY = rawSecondXY left + rawSecondXY right+      , rawSecondYY = rawSecondYY left + rawSecondYY right+      }++instance Monoid RawPlanarMomentSums where+  mempty = RawPlanarMomentSums 0 0 0 0 0 0++loopMomentSums :: ExactLoop -> RawPlanarMomentSums+loopMomentSums (ExactLoop points) =+  List.foldl' accumulateMomentEdge mempty (cyclePairs points)++accumulateMomentEdge+  :: RawPlanarMomentSums+  -> (ExactPoint, ExactPoint)+  -> RawPlanarMomentSums+accumulateMomentEdge accumulated (from, to) =+  let (fromX, fromY) = exactPointCoordinates from+      (toX, toY) = exactPointCoordinates to+      cross = exactPointCross from to+   in RawPlanarMomentSums+        { rawDoubleArea = rawDoubleArea accumulated + cross+        , rawFirstX = rawFirstX accumulated + (fromX + toX) * cross+        , rawFirstY = rawFirstY accumulated + (fromY + toY) * cross+        , rawSecondXX =+            rawSecondXX accumulated+              + (fromX * fromX + fromX * toX + toX * toX) * cross+        , rawSecondXY =+            rawSecondXY accumulated+              + (2 * fromX * fromY + fromX * toY + toX * fromY + 2 * toX * toY) * cross+        , rawSecondYY =+            rawSecondYY accumulated+              + (fromY * fromY + fromY * toY + toY * toY) * cross+        }+ data ExactLengthTerm = ExactLengthTerm   { lengthCoefficient :: !ExactRational   , squaredLength :: !ExactRational@@ -280,9 +406,7 @@ cellEdgeLengthContribution triangulation points selectedFaces edge = do   let (fromVertex, toVertex) = undirectedEndpoints triangulation edge       (forward, backward) = directedPair edge-      selected face =-        let FaceId raw = face-         in IntSet.member (fromIntegral raw) selectedFaces+      selected face = IntSet.member (faceIdIndex face) selectedFaces       coefficient = case (selected (incidentFace triangulation forward), selected (incidentFace triangulation backward)) of         (False, False) -> 1         (True, True) -> 0@@ -295,11 +419,11 @@   :: IntMap.IntMap ExactPoint   -> VertexId   -> Either ValuationError ExactPoint-cellPoint points vertex@(VertexId raw) =+cellPoint points vertex =   maybe     (Left (ValuationCoordinateMissing vertex))     Right-    (IntMap.lookup (fromIntegral raw) points)+    (IntMap.lookup (vertexIdIndex vertex) points)  triangleDoubleArea :: ExactPoint -> ExactPoint -> ExactPoint -> ExactRational triangleDoubleArea firstPoint secondPoint thirdPoint =@@ -660,3 +784,12 @@  oneHalf :: ExactRational oneHalf = exactRationalFromDyadic 1 (-1)++oneSixth :: ExactRational+oneSixth = exactRationalFromNormalizedRatio (1 Ratio.% 6)++oneTwelfth :: ExactRational+oneTwelfth = exactRationalFromNormalizedRatio (1 Ratio.% 12)++oneTwentyFourth :: ExactRational+oneTwentyFourth = exactRationalFromNormalizedRatio (1 Ratio.% 24)

@@ -3,7 +3,6 @@ {-# LANGUAGE DeriveAnyClass #-} {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE DerivingStrategies #-}-{-# LANGUAGE NamedFieldPuns #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies #-} @@ -54,14 +53,16 @@   ( DirectedEdgeId (..)   , asUndirected   , FaceId (..)+  , faceIdIndex   , VertexId (..)   , reverseEdge+  , vertexIdIndex   ) import Moonlight.Triangulation.Internal.InterpolationWorkspace+import Moonlight.Triangulation.Internal.Predicates (inCircle) import Moonlight.Triangulation.Math   ( barycentricCoordinates   , circumcenter-  , inCircle   , isFinite   , projectionFactor   , squaredDistanceWide@@ -135,8 +136,8 @@   !startDistance = squaredDistanceWide query (vertexPoint triangulation start)   !bound = numDirectedEdges triangulation + numVertices triangulation + 1 -  validateHint (Just vertex@(VertexId index))-    | fromIntegral index < numVertices triangulation = vertex+  validateHint (Just vertex)+    | vertexIdIndex vertex < numVertices triangulation = vertex   validateHint _ = VertexId 0    -- A vertex's Voronoi cell is the intersection of the half-planes its@@ -335,8 +336,8 @@   -> Point   -> FaceId   -> ST s (Int, Int)-discoverCavity workspace generation query (FaceId rawStart) = do-  MUV.unsafeWrite (nnFaceSeenMarks workspace) (fromIntegral rawStart) generation+discoverCavity workspace generation query startFace@(FaceId rawStart) = do+  MUV.unsafeWrite (nnFaceSeenMarks workspace) (faceIdIndex startFace) generation   MUV.unsafeWrite (nnFaceQueue workspace) 0 rawStart   go 1 0 0  where@@ -352,7 +353,7 @@         if not inside           then go slot cavitySize tests'           else do-            MUV.unsafeWrite (nnFaceMarks workspace) (fromIntegral rawFace) generation+            MUV.unsafeWrite (nnFaceMarks workspace) (faceIdIndex face) generation             MUV.unsafeWrite (nnCavityFaces workspace) cavitySize rawFace             nextQueue <- pushNeighbors triangulation workspace generation slot face             go nextQueue (cavitySize + 1) tests'@@ -395,7 +396,7 @@  where   pushOne !size edge =     let !adjacent@(FaceId raw) = incidentFace triangulation (reverseEdge edge)-        !index = fromIntegral raw+        !index = faceIdIndex adjacent      in if adjacent == outerFace           then pure size           else do@@ -427,10 +428,10 @@             appendIfBoundary size2 e2         goFaces (index + 1) nextSize   appendIfBoundary !size edge = do-    let adjacent@(FaceId rawAdjacent) = incidentFace triangulation (reverseEdge edge)+    let adjacent = incidentFace triangulation (reverseEdge edge)     outside <- if adjacent == outerFace       then pure True-      else (/= generation) <$> MUV.unsafeRead (nnFaceMarks workspace) (fromIntegral rawAdjacent)+      else (/= generation) <$> MUV.unsafeRead (nnFaceMarks workspace) (faceIdIndex adjacent)     if outside       then case edge of         DirectedEdgeId raw -> MUV.unsafeWrite (nnBoundaryEdges workspace) size raw >> pure (size + 1)@@ -452,7 +453,7 @@     | index >= count = pure ()     | otherwise = do         rawEdge <- MUV.unsafeRead (nnBoundaryEdges workspace) index-        let !from = vertexIndex (origin triangulation (DirectedEdgeId rawEdge))+        let !from = vertexIdIndex (origin triangulation (DirectedEdgeId rawEdge))         MUV.unsafeWrite (nnOriginMarks workspace) from generation         MUV.unsafeWrite (nnOriginEdge workspace) from rawEdge         install generation (index + 1)@@ -461,7 +462,7 @@         if currentRaw == firstRaw then pure count else pure 0     | otherwise = do         MUV.unsafeWrite (nnOrderedEdges workspace) index currentRaw-        let !to = vertexIndex (destination triangulation (DirectedEdgeId currentRaw))+        let !to = vertexIdIndex (destination triangulation (DirectedEdgeId currentRaw))         marked <- MUV.unsafeRead (nnOriginMarks workspace) to         if marked /= generation           then pure 0@@ -591,7 +592,7 @@   -> Point   -> FaceId   -> ST s (Maybe (Point))-cachedFaceCircumcenter workspace generation query face@(FaceId rawFace) = do+cachedFaceCircumcenter workspace generation query face = do   stamp <- MUV.unsafeRead (nnCircumcenterMarks workspace) slot   if stamp == generation     then do@@ -606,7 +607,7 @@         MUV.unsafeWrite (nnCircumcenterMarks workspace) slot generation         pure (Just center)  where-  !slot = fromIntegral rawFace+  !slot = faceIdIndex face  faceCircumcenterRelative   :: Triangulation mode vertex directed undirected face@@ -884,10 +885,6 @@         , alpha + c1Weight * radiusSquared         , beta + weight * radiusSquared         )--vertexIndex :: VertexId -> Int-vertexIndex (VertexId raw) = fromIntegral raw-{-# INLINE vertexIndex #-}  subtractPoint :: Point -> Point -> Point subtractPoint (Point ax ay) (Point bx by) = Point (ax - bx) (ay - by)

@@ -260,34 +260,17 @@       VoronoiSegment <$> faceCircumcenter triangulation fromFace <*> faceCircumcenter triangulation toFace   | innerFrom = do       start <- faceCircumcenter triangulation fromFace-      pure (VoronoiRay start (normalize (voronoiDirectionVector triangulation edge)))+      pure (VoronoiRay start (normalizeVector (voronoiDirectionVector triangulation edge)))   | innerTo = do       end <- faceCircumcenter triangulation toFace-      pure (VoronoiRay end (normalize (negatePoint (voronoiDirectionVector triangulation edge))))+      pure (VoronoiRay end (normalizeVector (negateVector (voronoiDirectionVector triangulation edge))))   | otherwise =       let !center = midpoint             (vertexPoint triangulation (origin triangulation primal))             (vertexPoint triangulation (destination triangulation primal))-       in Just (VoronoiLine center (normalize (voronoiDirectionVector triangulation edge)))+       in Just (VoronoiLine center (normalizeVector (voronoiDirectionVector triangulation edge)))  where   !fromFace = incidentFace triangulation primal   !toFace = incidentFace triangulation (reverseEdge primal)   !innerFrom = fromFace /= outerFace   !innerTo = toFace /= outerFace--normalize :: Point -> Point-normalize (Point x y)-  | scale == 0 = Point 0 0-  | otherwise =-      let !scaledX = x / scale-          !scaledY = y / scale-          !length' = sqrt (scaledX * scaledX + scaledY * scaledY)-       in Point (scaledX / length') (scaledY / length')- where-  !scale = max (abs x) (abs y)--negatePoint :: Point -> Point-negatePoint (Point x y) = Point (-x) (-y)--{-# INLINE normalize #-}-{-# INLINE negatePoint #-}

@@ -1,4 +1,3 @@-{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DataKinds #-}  -- | Zero-copy Voronoi views that retain the owning triangulation.@@ -37,7 +36,7 @@   ) where  import Moonlight.Triangulation.Handles.Dynamic-import Moonlight.Triangulation.Math (midpoint)+import Moonlight.Triangulation.Math (midpoint, negateVector, normalizeVector) import Moonlight.Triangulation.Types import Moonlight.Triangulation.Voronoi @@ -269,27 +268,13 @@       VoronoiSegment <$> innerFaceCircumcenter fromFace <*> innerFaceCircumcenter toFace     (InnerVoronoiVertexHandle fromFace, OuterVoronoiVertexHandle _) -> do       start <- innerFaceCircumcenter fromFace-      pure (VoronoiRay start (normalizePoint (voronoiEdgeDirectionH edge)))+      pure (VoronoiRay start (normalizeVector (voronoiEdgeDirectionH edge)))     (OuterVoronoiVertexHandle _, InnerVoronoiVertexHandle toFace) -> do       end <- innerFaceCircumcenter toFace-      pure (VoronoiRay end (normalizePoint (negatePoint (voronoiEdgeDirectionH edge))))+      pure (VoronoiRay end (normalizeVector (negateVector (voronoiEdgeDirectionH edge))))     (OuterVoronoiVertexHandle _, OuterVoronoiVertexHandle _) ->       let (from, to) = directedEdgePositions primal-       in Just (VoronoiLine (midpoint from to) (normalizePoint (voronoiEdgeDirectionH edge)))--normalizePoint :: Point -> Point-normalizePoint (Point x y)-  | scale == 0 = Point 0 0-  | otherwise =-      let !scaledX = x / scale-          !scaledY = y / scale-          !length' = sqrt (scaledX * scaledX + scaledY * scaledY)-       in Point (scaledX / length') (scaledY / length')- where-  !scale = max (abs x) (abs y)--negatePoint :: Point -> Point-negatePoint (Point x y) = Point (-x) (-y)+       in Just (VoronoiLine (midpoint from to) (normalizeVector (voronoiEdgeDirectionH edge)))  -- The handle verbs cross the same component boundary as the fixed-index ones -- and carry their unfoldings for the same reason. 'voronoiVertexPositionH'

@@ -47,7 +47,7 @@   , unorderedPairs   ) import Moonlight.Triangulation.Internal.ExactRational (ExactRational)-import qualified Moonlight.Triangulation.Math as Math+import qualified Moonlight.Triangulation.Internal.Predicates as Admitted   ( orient2d   , segmentRelation   )@@ -459,7 +459,7 @@   exactRelation =     relationFor exactSegmentRelation projectedExactPoint firstSegment secondSegment   roundedRelation =-    relationFor Math.segmentRelation roundedProjectedPoint firstSegment secondSegment+    relationFor Admitted.segmentRelation roundedProjectedPoint firstSegment secondSegment  relationFor   :: (point -> point -> point -> point -> SegmentRelation)@@ -507,7 +507,7 @@   -> Ordering compareRoundedAround center left right =   compareRoundedAroundWith-    ( Math.orient2d+    ( Admitted.orient2d         (roundedVertexPoint center)         (roundedVertexPoint left)         (roundedVertexPoint right)

@@ -119,6 +119,7 @@   , insertVertex   , insertVertexAt   , insertVertexAtNear+  , insertVertexAtNearVertex   , removeAt   , removeAtNear   , removeManyAt@@ -493,7 +494,16 @@   , EulerCharacteristic   , eulerCharacteristicValue   , ExactArea+  , ExactPlanarMoments   , exactAreaValue+  , exactPlanarMeasure+  , exactPlanarFirstX+  , exactPlanarFirstY+  , exactPlanarSecondXX+  , exactPlanarSecondXY+  , exactPlanarSecondYY+  , scaleExactPlanarMoments+  , polygonComponentMoments   , polygonComponentArea   , ExactLengthTerm   , lengthCoefficient@@ -533,11 +543,13 @@   , PowerWeight   , PowerWeightError (..)   , powerWeight+  , powerWeightFromExact   , powerWeightExact   , PowerSite   , powerSite   , powerSiteLabel   , powerSitePosition+  , powerSiteExactPosition   , powerSiteWeight   , PowerCellDisposition (..)   , RegularSiteDisposition (..)@@ -549,14 +561,26 @@   , regularEdgeLabels   , regularEdgeDual   , RegularTriangulation+  , emptyRegularTriangulation   , regularTriangulation+  , regularSite+  , regularSites+  , regularSiteCount   , regularSiteDisposition   , regularFaces   , regularEdges   , regularNeighbours+  , regularTriangulationReceipt   , RegularTriangulationReceipt (..)+  , RegularSiteTransition (..)+  , RegularEditError (..)+  , RegularEditResult (..)+  , insertRegularSite+  , removeRegularSite+  , reweightRegularSites   , BoundedPowerDiagram   , boundedPowerDiagram+  , boundedPowerDiagramFromRegular   , powerCellDisposition   , powerCellDispositions   , powerDiagramPlanarLayer@@ -569,6 +593,18 @@   , AffineForm (..)   , UpperEnvelopeError (..)   , upperEnvelope+  , PowerAlphaBirth+  , powerAlphaBirthExact+  , powerAlphaBirthNumerator+  , powerAlphaBirthDenominator+  , RegularAlphaFiltration+  , RegularAlphaError (..)+  , regularAlphaFiltration+  , regularAlphaComplex+  , regularAlphaBirths+  , regularAlphaSimplexBirth+  , regularAlphaCriticalBirths+  , regularAlphaComplexAtBirth      -- * Discharge — the invariants the constructors guarantee, checkable on a     -- value built by any route; every violation is a value carrying its witness@@ -1023,13 +1059,19 @@   , PowerWeight   , PowerWeightError (..)   , RegularEdge+  , RegularEditError (..)+  , RegularEditResult (..)   , RegularFace   , RegularSiteDisposition (..)+  , RegularSiteTransition (..)   , RegularTriangulation   , RegularTriangulationReceipt (..)   , RegularTopologyError (..)   , UpperEnvelopeError (..)   , boundedPowerDiagram+  , boundedPowerDiagramFromRegular+  , emptyRegularTriangulation+  , insertRegularSite   , powerCellDisposition   , powerCellDispositions   , powerDiagramPlanarLayer@@ -1039,8 +1081,10 @@   , powerSite   , powerSiteLabel   , powerSitePosition+  , powerSiteExactPosition   , powerSiteWeight   , powerWeight+  , powerWeightFromExact   , powerWeightExact   , regularEdgeDual   , regularEdgeLabels@@ -1049,10 +1093,30 @@   , regularFaceLabels   , regularFaces   , regularNeighbours+  , regularSite+  , regularSiteCount   , regularSiteDisposition+  , regularSites   , regularTriangulation+  , regularTriangulationReceipt+  , removeRegularSite+  , reweightRegularSites   , upperEnvelope   )+import Moonlight.Triangulation.RegularAlpha+  ( PowerAlphaBirth+  , RegularAlphaError (..)+  , RegularAlphaFiltration+  , powerAlphaBirthDenominator+  , powerAlphaBirthExact+  , powerAlphaBirthNumerator+  , regularAlphaBirths+  , regularAlphaComplex+  , regularAlphaComplexAtBirth+  , regularAlphaCriticalBirths+  , regularAlphaFiltration+  , regularAlphaSimplexBirth+  ) import Moonlight.Triangulation.Payload   ( PayloadTraversal   , directedPayloads@@ -1122,6 +1186,7 @@   , insertVertex   , insertVertexAt   , insertVertexAtNear+  , insertVertexAtNearVertex   , refuse   , removeAt   , removeAtNear@@ -1141,6 +1206,7 @@   ( CertifiedInterval (..)   , EulerCharacteristic   , ExactArea+  , ExactPlanarMoments   , ExactLengthExpression   , ExactLengthMeasurement   , ExactLengthTerm@@ -1150,7 +1216,14 @@   , cellValuations   , eulerCharacteristicValue   , exactAreaValue+  , exactPlanarFirstX+  , exactPlanarFirstY+  , exactPlanarMeasure+  , exactPlanarSecondXX+  , exactPlanarSecondXY+  , exactPlanarSecondYY   , polygonComponentArea+  , polygonComponentMoments   , exactLengthBounds   , exactLengthExpression   , exactLengthTerms@@ -1158,6 +1231,7 @@   , planarValuationsPerimeter   , regionPerimeter   , regionValuations+  , scaleExactPlanarMoments   , squaredLength   , valuationArea   , valuationEuler

@@ -32,7 +32,11 @@ import Data.Word (Word32) import Moonlight.Triangulation.BulkLoad (delaunay, insert) import Moonlight.Triangulation.Dcel (numUndirectedEdges, numVertices, undirectedEndpoints, vertexPoint)-import Moonlight.Triangulation.Internal.HandleDefs (UndirectedEdgeId (..), VertexId (..))+import Moonlight.Triangulation.Internal.HandleDefs+  ( UndirectedEdgeId (..)+  , VertexId (..)+  , vertexIdIndex+  ) import Moonlight.Triangulation.Interpolation (nearestNeighbor) import Moonlight.Triangulation.Math (canonicalPoint, validatePoint) import Moonlight.Triangulation.Removal (RemovalOutcome, removalTriangulation, removeVertex)@@ -239,7 +243,7 @@   case disposition of     AlreadyPresent -> Right hierarchy     Inserted-      | vertexIndex vertex /= baseIndex ->+      | vertexIdIndex vertex /= baseIndex ->           Left             ( HierarchyInsertionHandleMismatch                 (VertexId (fromIntegral baseIndex))@@ -277,7 +281,7 @@           )     | otherwise = do         inserted <- insert level point-        if vertexIndex (insertionVertex inserted) /= numVertices level+        if vertexIdIndex (insertionVertex inserted) /= numVertices level           then             Left               ( HierarchyInsertionHandleMismatch@@ -348,7 +352,7 @@   -- The freed slot is the removed vertex's own index, and it is where the   -- former last vertex now stands. A removal that took the last vertex frees   -- no slot and reports none, and its removed index is that last index.-  !freedSlot = maybe surviving (vertexIndex . fst) swap+  !freedSlot = maybe surviving (vertexIdIndex . fst) swap   !vacatedIndex = surviving    -- The nesting law is stated here as an equation for the same reason the@@ -473,9 +477,6 @@ samplePopulation count divisor   | count <= 0 = 0   | otherwise = 1 + (count - 1) `quot` divisor--vertexIndex :: VertexId -> Int-vertexIndex (VertexId value) = fromIntegral value  -- The hierarchy is consulted once per query from another package, so its -- polymorphic entries expose their unfoldings for the same reason the search

@@ -0,0 +1,1376 @@+-- | Pure topology-changing regular edits by exact conflict-cavity descent and gluing.+module Moonlight.Triangulation.Internal.PowerDiagram.Edit where++import Control.Monad (when)+import Data.Bifunctor (first)+import qualified Data.Foldable as Foldable+import qualified Data.List as List+import qualified Data.List.NonEmpty as NonEmpty+import qualified Data.Map.Strict as Map+import Data.Map.Strict (Map)+import Data.Maybe (isJust)+import qualified Data.Set as Set+import Data.Set (Set)+import qualified Data.Vector as Vector+import Moonlight.Triangulation.Exact+  ( exactOrient2d+  , exactPoint+  )+import Moonlight.Triangulation.Internal.BoundaryCycle (orderedPair)+import Moonlight.Triangulation.Internal.ExactRational (ExactRational)+import Moonlight.Triangulation.Internal.PowerDiagram.Generator+import Moonlight.Triangulation.Internal.PowerDiagram.Hull (regularGeneratorTopology)+import Moonlight.Triangulation.Internal.PowerDiagram.Locality+import Moonlight.Triangulation.Internal.PowerDiagram.Model+import Moonlight.Triangulation.Internal.PowerDiagram.Section++-- | Insert one stable-labelled site. Repeating the identical site is+-- idempotent; reusing its label for different geometry is a typed conflict.+insertRegularSite+  :: Ord label+  => PowerSite label+  -> RegularTriangulation label+  -> Either+      (RegularEditError label)+      (RegularEditResult label)+insertRegularSite site triangulation =+  let label = powerSiteLabel site+   in case Map.lookup label (storedRegularSites triangulation) of+        Nothing ->+          withRegularEditFallback+            triangulation+            (Set.singleton label)+            ( normalizeRegularEdit+                triangulation+                (Map.insert label site (storedRegularSites triangulation))+                (Set.singleton label)+            )+            (insertRegularSiteLocally site triangulation)+        Just resident+          | resident == site ->+              Right (unchangedRegularEdit triangulation)+          | otherwise -> Left (RegularEditSiteConflict label resident site)++-- | Remove one labelled site. Absence is an idempotent successful result.+removeRegularSite+  :: Ord label+  => label+  -> RegularTriangulation label+  -> Either+      (RegularEditError label)+      (RegularEditResult label)+removeRegularSite label triangulation =+  case Map.lookup label (storedRegularSites triangulation) of+    Nothing ->+      Right (unchangedRegularEdit triangulation)+    Just _ ->+      case removeTopologyPreservingSite label triangulation of+        Just (edited, transitions) ->+          Right+            ( RegularEditResult+                edited+                (Set.singleton label)+                transitions+            )+        Nothing ->+          withRegularEditFallback+            triangulation+            (Set.singleton label)+            ( normalizeRegularEdit+                triangulation+                (Map.delete label (storedRegularSites triangulation))+                (Set.singleton label)+            )+            (removeRegularSiteLocally label triangulation)++-- | Replace many weights in one normalization. Every label must already+-- belong to the triangulation; an empty or unchanged patch is idempotent.+reweightRegularSites+  :: Ord label+  => Map label PowerWeight+  -> RegularTriangulation label+  -> Either+      (RegularEditError label)+      (RegularEditResult label)+reweightRegularSites requested triangulation =+  case NonEmpty.nonEmpty unknownLabels of+    Just unknown -> Left (RegularEditUnknownSites unknown)+    Nothing+      | Map.null changedSites ->+          Right (unchangedRegularEdit triangulation)+      | otherwise ->+          case reweightTopologyPreservingSites changedLabels changedSites reweightedSites triangulation of+            Just (edited, transitions) ->+              Right+                ( RegularEditResult+                    edited+                    changedLabels+                    transitions+                )+            Nothing ->+              case commonRegularWeightShift requested sites of+                Just shift ->+                  Right+                    ( translateRegularWeights+                        shift+                        reweightedSites+                        changedLabels+                        triangulation+                    )+                Nothing ->+                  case Map.elems changedSites of+                    [revised] ->+                      withRegularEditFallback+                        triangulation+                        changedLabels+                        normalizeGlobally+                        (reweightRegularSiteLocally revised triangulation)+                    _ -> normalizeGlobally+ where+  sites = storedRegularSites triangulation+  unknownLabels = Map.keys (Map.difference requested sites)+  changedSites =+    Map.mapMaybeWithKey+      (\label weight ->+          case Map.lookup label sites of+            Just site+              | powerSiteWeight site /= weight ->+                  Just site{powerSiteWeight = weight}+            _ -> Nothing+      )+      requested+  reweightedSites = Map.union changedSites sites+  changedLabels = Map.keysSet changedSites+  normalizeGlobally =+    normalizeRegularEdit triangulation reweightedSites changedLabels++withRegularEditFallback+  :: Ord label+  => RegularTriangulation label+  -> Set label+  -> Either (RegularEditError label) (RegularEditResult label)+  -> Either obstruction (RegularTriangulation label, Set label)+  -> Either (RegularEditError label) (RegularEditResult label)+withRegularEditFallback before changed fallback =+  either (const fallback) (Right . publishLocalRegularEdit before changed)+{-# INLINE withRegularEditFallback #-}++publishLocalRegularEdit+  :: Ord label+  => RegularTriangulation label+  -> Set label+  -> (RegularTriangulation label, Set label)+  -> RegularEditResult label+publishLocalRegularEdit before changed (edited, transitionSupport) =+  RegularEditResult+    edited+    changed+    ( dispositionTransitionsOn+        transitionSupport+        (regularDispositionSection before)+        (regularDispositionSection edited)+    )+{-# INLINE publishLocalRegularEdit #-}++insertRegularSiteLocally+  :: Ord label+  => PowerSite label+  -> RegularTriangulation label+  -> Either+      (LocalRegularEditObstruction label)+      (RegularTriangulation label, Set label)+insertRegularSiteLocally site triangulation = do+  section <-+    maybe (Left LocalRegularSectionUnavailable) Right+      (storedRegularSection triangulation)+  let (generator, _) = prepareExactPowerGenerator site+      slope = exactGeneratorSlope generator+  case Map.lookup slope (sectionSlopeRepresentatives section) of+    Just representative ->+      insertCoincidentRegularSite site generator representative section triangulation+    Nothing ->+      insertDistinctRegularSite site generator section triangulation++insertCoincidentRegularSite+  :: Ord label+  => PowerSite label+  -> ExactPowerGenerator label+  -> label+  -> RegularSection label+  -> RegularTriangulation label+  -> Either+      (LocalRegularEditObstruction label)+      (RegularTriangulation label, Set label)+insertCoincidentRegularSite site generator representativeLabel section triangulation = do+  representative <-+    requireGeneratorFromMap (sectionGenerators section) representativeLabel+  let label = powerSiteLabel site+  case+    compareCoincidentPriority+      (exactPowerGeneratorConstant generator, label)+      (exactPowerGeneratorConstant representative, representativeLabel) of+    LT ->+      let coincident = classifyCoincidentGenerator representative generator+          disposition = coincidentRegularDisposition coincident+          editedSection =+            section+              { sectionCoincidentDispositions =+                  Map.insert label coincident (sectionCoincidentDispositions section)+              , sectionRegularDispositions =+                  Map.insert label disposition (sectionRegularDispositions section)+              , sectionRegularStars =+                  Map.insert label emptyRegularSiteStar (sectionRegularStars section)+              }+       in Right+            ( replaceRegularSiteSection+                (Map.insert label site (storedRegularSites triangulation))+                editedSection+                triangulation+            , Set.singleton label+            )+    _ -> Left (LocalRegularSlopeCoincidence representativeLabel)++insertDistinctRegularSite+  :: Ord label+  => PowerSite label+  -> ExactPowerGenerator label+  -> RegularSection label+  -> RegularTriangulation label+  -> Either+      (LocalRegularEditObstruction label)+      (RegularTriangulation label, Set label)+insertDistinctRegularSite site generator section triangulation = do+  locality <-+    maybe (Left LocalRegularLocalityUnavailable) Right+      (sectionRegularLocality section)+  location <-+    maybe (Left LocalRegularSlopeLocationFailed) Right+      (locateRegularSlope (localityFaceSeed locality) section generator)+  case location of+    RegularSlopeOutside edgeKey -> Left (LocalRegularSlopeOutside edgeKey)+    RegularSlopeInside faceKey -> do+      volume <- regularFaceVolume section faceKey generator+      case compare volume 0 of+        LT -> insertHiddenRegularSite site generator faceKey locality section triangulation+        EQ -> Left (LocalRegularDegenerateFace faceKey)+        GT -> insertVisibleRegularSite site generator faceKey locality section triangulation++insertHiddenRegularSite+  :: Ord label+  => PowerSite label+  -> ExactPowerGenerator label+  -> RegularFaceKey label+  -> RegularLocality label+  -> RegularSection label+  -> RegularTriangulation label+  -> Either+      (LocalRegularEditObstruction label)+      (RegularTriangulation label, Set label)+insertHiddenRegularSite site generator supportFace locality section triangulation =+  let label = powerSiteLabel site+      receipt = sectionRegularReceipt section+      editedLocality = attachRegularHidden label supportFace locality+      editedSection =+        section+          { sectionGenerators = Map.insert label generator (sectionGenerators section)+          , sectionSlopeRepresentatives =+              Map.insert (exactGeneratorSlope generator) label (sectionSlopeRepresentatives section)+          , sectionRegularDispositions =+              Map.insert label RegularSiteHidden (sectionRegularDispositions section)+          , sectionRegularStars =+              Map.insert label emptyRegularSiteStar (sectionRegularStars section)+          , sectionRegularReceipt =+              receipt+                { generatorRegularInputSites = generatorRegularInputSites receipt + 1+                , generatorRegularHiddenSites = generatorRegularHiddenSites receipt + 1+                }+          , sectionRegularLocality = Just editedLocality+          }+   in Right+        ( replaceRegularSiteSection+            (Map.insert label site (storedRegularSites triangulation))+            editedSection+            triangulation+        , Set.singleton label+        )++insertVisibleRegularSite+  :: Ord label+  => PowerSite label+  -> ExactPowerGenerator label+  -> RegularFaceKey label+  -> RegularLocality label+  -> RegularSection label+  -> RegularTriangulation label+  -> Either+      (LocalRegularEditObstruction label)+      (RegularTriangulation label, Set label)+insertVisibleRegularSite site generator seedFace locality section triangulation = do+  conflictFaces <- regularConflictCavity section generator seedFace+  horizon <- regularConflictHorizon section conflictFaces+  when (Map.null horizon) (Left LocalRegularEmptyHorizon)+  newFaces <- traverse (regularHorizonFace section generator) (Map.elems horizon)+  let label = powerSiteLabel site+      generators = Map.insert label generator (sectionGenerators section)+  facePatch <-+    applyRegularFacePatch+      generators+      conflictFaces+      newFaces+      section+  certifyRegularFacePatch generators facePatch+  let changedTopologyLabels = Set.insert label (patchedRegularLabels facePatch)+      provisionalDispositions =+        Set.foldl'+          (publishPatchedRegularDisposition (patchedRegularStars facePatch))+          (Map.insert label RegularSiteVisible (sectionRegularDispositions section))+          changedTopologyLabels+      newlyHidden =+        Set.filter+          (becameRegularHidden (sectionRegularDispositions section) provisionalDispositions)+          changedTopologyLabels+      hiddenOnConflict =+        Set.unions+          [ Map.findWithDefault Set.empty faceKey (localitySupportHidden locality)+          | faceKey <- Set.toAscList conflictFaces+          ]+      hiddenToAttach = Set.union newlyHidden hiddenOnConflict+      receipt = sectionRegularReceipt section+      provisionalSection =+        section+          { sectionGenerators = generators+          , sectionSlopeRepresentatives =+              Map.insert (exactGeneratorSlope generator) label (sectionSlopeRepresentatives section)+          , sectionRegularDispositions = provisionalDispositions+          , sectionRegularFaces = patchedRegularFaces facePatch+          , sectionRegularEdges = patchedRegularEdges facePatch+          , sectionRegularStars = patchedRegularStars facePatch+          , sectionRegularReceipt =+              receipt+                { generatorRegularInputSites = generatorRegularInputSites receipt + 1+                , generatorRegularVisibleSites =+                    generatorRegularVisibleSites receipt + 1 - Set.size newlyHidden+                , generatorRegularHiddenSites =+                    generatorRegularHiddenSites receipt + Set.size newlyHidden+                , generatorRegularFaceCount = Map.size (patchedRegularFaces facePatch)+                , generatorRegularEdgeCount = Map.size (patchedRegularEdges facePatch)+                }+          , sectionRegularLocality = Nothing+          }+  editedLocality <-+    patchRegularLocality locality conflictFaces hiddenToAttach provisionalSection+  let editedSection = provisionalSection{sectionRegularLocality = Just editedLocality}+  pure+    ( replaceRegularSiteSection+        (Map.insert label site (storedRegularSites triangulation))+        editedSection+        triangulation+    , Set.insert label newlyHidden+    )+++regularFaceVolume+  :: Ord label+  => RegularSection label+  -> RegularFaceKey label+  -> ExactPowerGenerator label+  -> Either (LocalRegularEditObstruction label) ExactRational+regularFaceVolume section faceKey candidate = do+  face <- requireFaceFromMap (sectionRegularFaces section) faceKey+  (firstGenerator, secondGenerator, thirdGenerator) <-+    regularFaceGenerators (sectionGenerators section) face+  pure+    ( exactGeneratorLiftedVolume+        firstGenerator+        secondGenerator+        thirdGenerator+        candidate+    )++regularConflictCavity+  :: Ord label+  => RegularSection label+  -> ExactPowerGenerator label+  -> RegularFaceKey label+  -> Either+      (LocalRegularEditObstruction label)+      (Set (RegularFaceKey label))+regularConflictCavity section candidate seed =+  let initial =+        RegularConflictDescent+          { conflictPendingFaces = Set.singleton seed+          , conflictVisitedFaces = Set.empty+          , conflictVisibleFaces = Set.empty+          , conflictObstruction = Nothing+          }+      descended =+        Map.foldl'+          (\current _ -> advanceRegularConflict section candidate current)+          initial+          (sectionRegularFaces section)+   in case conflictObstruction descended of+        Just obstruction -> Left obstruction+        Nothing+          | not (Set.null (conflictPendingFaces descended)) ->+              Left LocalRegularSlopeLocationFailed+          | Set.null (conflictVisibleFaces descended) ->+              Left (LocalRegularDegenerateFace seed)+          | otherwise -> Right (conflictVisibleFaces descended)++advanceRegularConflict+  :: Ord label+  => RegularSection label+  -> ExactPowerGenerator label+  -> RegularConflictDescent label+  -> RegularConflictDescent label+advanceRegularConflict _ _ descent+  | Set.null (conflictPendingFaces descent)+      || isJust (conflictObstruction descent) = descent+advanceRegularConflict section candidate descent =+  case Set.minView (conflictPendingFaces descent) of+    Nothing -> descent+    Just (faceKey, remaining) ->+      case regularFaceVolume section faceKey candidate of+        Left obstruction ->+          descent+            { conflictPendingFaces = remaining+            , conflictObstruction = Just obstruction+            }+        Right volume ->+          let visited = Set.insert faceKey (conflictVisitedFaces descent)+           in case compare volume 0 of+                LT ->+                  descent+                    { conflictPendingFaces = remaining+                    , conflictVisitedFaces = visited+                    }+                EQ ->+                  descent+                    { conflictPendingFaces = remaining+                    , conflictVisitedFaces = visited+                    , conflictObstruction = Just (LocalRegularDegenerateFace faceKey)+                    }+                GT ->+                  case regularFaceNeighbours section faceKey of+                    Left obstruction ->+                      descent+                        { conflictPendingFaces = remaining+                        , conflictVisitedFaces = visited+                        , conflictObstruction = Just obstruction+                        }+                    Right neighbours ->+                      descent+                        { conflictPendingFaces =+                            Set.union remaining (Set.difference neighbours visited)+                        , conflictVisitedFaces = visited+                        , conflictVisibleFaces =+                            Set.insert faceKey (conflictVisibleFaces descent)+                        }+++regularConflictHorizon+  :: Ord label+  => RegularSection label+  -> Set (RegularFaceKey label)+  -> Either+      (LocalRegularEditObstruction label)+      (Map (RegularEdgeKey label) (label, label))+regularConflictHorizon section conflictFaces = do+  boundarySections <-+    traverse+      (regularConflictFaceHorizon section conflictFaces)+      (Set.toAscList conflictFaces)+  pure (Map.unions boundarySections)++regularConflictFaceHorizon+  :: Ord label+  => RegularSection label+  -> Set (RegularFaceKey label)+  -> RegularFaceKey label+  -> Either+      (LocalRegularEditObstruction label)+      (Map (RegularEdgeKey label) (label, label))+regularConflictFaceHorizon section conflictFaces faceKey = do+  face <- requireFaceFromMap (sectionRegularFaces section) faceKey+  associations <-+    traverse+      (classifyConflictEdge section conflictFaces)+      (regularFaceDirectedEdgeKeys face)+  pure (Map.fromList (concat associations))++classifyConflictEdge+  :: Ord label+  => RegularSection label+  -> Set (RegularFaceKey label)+  -> (label, label)+  -> Either+      (LocalRegularEditObstruction label)+      [(RegularEdgeKey label, (label, label))]+classifyConflictEdge section conflictFaces directed@(fromLabel, toLabel) = do+  let edgeKey = orderedPair fromLabel toLabel+  incident <- requireRegularEdgeIncidence section edgeKey+  if Set.size incident > 2+    then Left (LocalRegularNonManifoldEdge edgeKey (Set.size incident))+    else+      if Set.size (Set.intersection incident conflictFaces) == 1+        then Right [(edgeKey, directed)]+        else Right []++regularHorizonFace+  :: Ord label+  => RegularSection label+  -> ExactPowerGenerator label+  -> (label, label)+  -> Either (LocalRegularEditObstruction label) (RegularFace label)+regularHorizonFace section candidate (fromLabel, toLabel) = do+  fromGenerator <- requireGeneratorFromMap (sectionGenerators section) fromLabel+  toGenerator <- requireGeneratorFromMap (sectionGenerators section) toLabel+  let candidateLabel = exactPowerGeneratorLabel candidate+      faceKey = sortedRegularFaceLabels fromLabel toLabel candidateLabel+      fromPoint = uncurry exactPoint (exactGeneratorSlope fromGenerator)+      toPoint = uncurry exactPoint (exactGeneratorSlope toGenerator)+      candidatePoint = uncurry exactPoint (exactGeneratorSlope candidate)+  case exactOrient2d fromPoint toPoint candidatePoint of+    GT -> do+      dual <-+        first LocalRegularTopologyObstruction+          (exactGeneratorFaceDual fromGenerator toGenerator candidate)+      pure (canonicalRegularFace fromLabel toLabel candidateLabel dual)+    _ -> Left (LocalRegularDegenerateFace faceKey)++applyRegularFacePatch+  :: Ord label+  => Map label (ExactPowerGenerator label)+  -> Set (RegularFaceKey label)+  -> [RegularFace label]+  -> RegularSection label+  -> Either (LocalRegularEditObstruction label) (RegularFacePatch label)+applyRegularFacePatch generators removedFaceKeys newFaces section = do+  removedFaces <-+    traverse+      (requireFaceFromMap (sectionRegularFaces section))+      (Set.toAscList removedFaceKeys)+  let newFaceSection =+        Map.fromList [(regularFaceKey face, face) | face <- newFaces]+      faces =+        Map.union newFaceSection+          (Map.withoutKeys (sectionRegularFaces section) removedFaceKeys)+      touchedEdges =+        Set.unions+          [ Set.unions (fmap (Set.fromList . regularFaceCanonicalEdgeKeys) removedFaces)+          , Set.unions (fmap (Set.fromList . regularFaceCanonicalEdgeKeys) newFaces)+          ]+      newIncidence =+        List.foldl'+          (Map.unionWith Set.union)+          Map.empty+          (fmap regularFaceEdgeIncidence newFaces)+      affectedLabels =+        Set.unions+          [ Set.unions (fmap regularFaceLabelSet removedFaces)+          , Set.unions (fmap regularFaceLabelSet newFaces)+          ]+  edges <-+    Foldable.foldlM+      (patchRegularEdge generators faces removedFaceKeys newIncidence)+      (sectionRegularEdges section)+      touchedEdges+  let stars =+        patchRegularStars+          affectedLabels+          removedFaceKeys+          newFaceSection+          touchedEdges+          edges+          section+  pure+    RegularFacePatch+      { patchedRegularFaces = faces+      , patchedRegularEdges = edges+      , patchedRegularStars = stars+      , patchedRegularLabels = affectedLabels+      , patchedRegularTouchedEdges = touchedEdges+      }++patchRegularEdge+  :: Ord label+  => Map label (ExactPowerGenerator label)+  -> Map (RegularFaceKey label) (RegularFace label)+  -> Set (RegularFaceKey label)+  -> Map (RegularEdgeKey label) (Set (RegularFaceKey label))+  -> Map (RegularEdgeKey label) (RegularEdgeSection label)+  -> RegularEdgeKey label+  -> Either+      (LocalRegularEditObstruction label)+      (Map (RegularEdgeKey label) (RegularEdgeSection label))+patchRegularEdge generators faces removedFaceKeys newIncidence edges edgeKey =+  let oldIncidence =+        maybe Set.empty sectionRegularEdgeFaces (Map.lookup edgeKey edges)+      incidence =+        Set.union+          (Set.difference oldIncidence removedFaceKeys)+          (Map.findWithDefault Set.empty edgeKey newIncidence)+   in if Set.null incidence+        then Right (Map.delete edgeKey edges)+        else do+          edgeSection <- publishRegularEdgeSection generators faces edgeKey incidence+          pure (Map.insert edgeKey edgeSection edges)++publishRegularEdgeSection+  :: Ord label+  => Map label (ExactPowerGenerator label)+  -> Map (RegularFaceKey label) (RegularFace label)+  -> RegularEdgeKey label+  -> Set (RegularFaceKey label)+  -> Either (LocalRegularEditObstruction label) (RegularEdgeSection label)+publishRegularEdgeSection generators faces edgeKey@(firstLabel, secondLabel) incidence =+  case Set.toAscList incidence of+    [faceKey] -> do+      face <- requireFaceFromMap faces faceKey+      (fromLabel, toLabel, oppositeLabel) <-+        requireFaceEdgeWitness edgeKey face+      fromGenerator <- requireGeneratorFromMap generators fromLabel+      toGenerator <- requireGeneratorFromMap generators toLabel+      oppositeGenerator <- requireGeneratorFromMap generators oppositeLabel+      ray <-+        first LocalRegularTopologyObstruction+          ( exactGeneratorBoundaryDualRay+              (regularFaceDualPoint face)+              fromGenerator+              toGenerator+              oppositeGenerator+          )+      pure+        ( RegularEdgeSection+            (RegularEdge firstLabel secondLabel (UnboundedPowerDual ray))+            incidence+        )+    [firstFaceKey, secondFaceKey] -> do+      firstFace <- requireFaceFromMap faces firstFaceKey+      secondFace <- requireFaceFromMap faces secondFaceKey+      let firstDual = regularFaceDualPoint firstFace+          secondDual = regularFaceDualPoint secondFace+      dual <-+        fmap publishGeneratorDualGeometry+          ( first LocalRegularTopologyObstruction+              ( exactGeneratorDualBetween+                  firstLabel+                  secondLabel+                  firstDual+                  secondDual+              )+          )+      pure+        ( RegularEdgeSection+            (RegularEdge firstLabel secondLabel dual)+            incidence+        )+    incidentFaces ->+      Left (LocalRegularNonManifoldEdge edgeKey (length incidentFaces))++patchRegularStars+  :: Ord label+  => Set label+  -> Set (RegularFaceKey label)+  -> Map (RegularFaceKey label) (RegularFace label)+  -> Set (RegularEdgeKey label)+  -> Map (RegularEdgeKey label) (RegularEdgeSection label)+  -> RegularSection label+  -> Map label (RegularSiteStar label)+patchRegularStars affectedLabels removedFaceKeys newFaceSection touchedEdges edges section =+  Set.foldl' patchStar (sectionRegularStars section) affectedLabels+ where+  affectedSupport = Map.fromSet (const ()) affectedLabels+  touchedStars =+    regularSiteStars+      affectedSupport+      Map.empty+      ( Map.restrictKeys+          (Map.union edges (sectionRegularEdges section))+          touchedEdges+      )+  addedStars =+    regularSiteStars+      affectedSupport+      newFaceSection+      (Map.restrictKeys edges touchedEdges)+  patchStar stars label =+    let oldStar = Map.findWithDefault emptyRegularSiteStar label stars+        touchedStar = Map.findWithDefault emptyRegularSiteStar label touchedStars+        addedStar = Map.findWithDefault emptyRegularSiteStar label addedStars+        incidentFaces =+          Set.union+            (Set.difference (sectionIncidentFaces oldStar) removedFaceKeys)+            (sectionIncidentFaces addedStar)+        retainedNeighbours =+          Set.difference+            (sectionSiteNeighbours oldStar)+            (sectionSiteNeighbours touchedStar)+     in Map.insert+          label+          ( RegularSiteStar+              incidentFaces+              (Set.union retainedNeighbours (sectionSiteNeighbours addedStar))+          )+          stars++certifyRegularFacePatch+  :: Ord label+  => Map label (ExactPowerGenerator label)+  -> RegularFacePatch label+  -> Either (LocalRegularEditObstruction label) ()+certifyRegularFacePatch generators patch =+  Foldable.traverse_+    (\edgeKey ->+       case Map.lookup edgeKey (patchedRegularEdges patch) of+         Nothing -> Right ()+         Just edgeSection ->+           certifyRegularEdgeConvexity+             generators+             (patchedRegularFaces patch)+             (edgeKey, edgeSection))+    (patchedRegularTouchedEdges patch)++certifyRegularEdgeConvexity+  :: Ord label+  => Map label (ExactPowerGenerator label)+  -> Map (RegularFaceKey label) (RegularFace label)+  -> (RegularEdgeKey label, RegularEdgeSection label)+  -> Either (LocalRegularEditObstruction label) ()+certifyRegularEdgeConvexity generators faces (edgeKey, edgeSection) =+  case Set.toAscList (sectionRegularEdgeFaces edgeSection) of+    [_] -> Right ()+    [firstFaceKey, secondFaceKey] -> do+      firstFace <- requireFaceFromMap faces firstFaceKey+      secondFace <- requireFaceFromMap faces secondFaceKey+      (_, _, secondOpposite) <- requireFaceEdgeWitness edgeKey secondFace+      secondOppositeGenerator <- requireGeneratorFromMap generators secondOpposite+      (firstGenerator, secondGenerator, thirdGenerator) <-+        regularFaceGenerators generators firstFace+      if exactGeneratorLiftedVolume+          firstGenerator+          secondGenerator+          thirdGenerator+          secondOppositeGenerator+          < 0+        then Right ()+        else Left (LocalRegularNonConvexEdge edgeKey)+    incident -> Left (LocalRegularNonManifoldEdge edgeKey (length incident))++publishPatchedRegularDisposition+  :: Ord label+  => Map label (RegularSiteStar label)+  -> Map label (RegularSiteDisposition label)+  -> label+  -> Map label (RegularSiteDisposition label)+publishPatchedRegularDisposition stars dispositions label =+  case Map.lookup label dispositions of+    Just RegularSiteVisible+      | maybe True (Set.null . sectionIncidentFaces) (Map.lookup label stars) ->+          Map.insert label RegularSiteHidden dispositions+    _ -> dispositions++becameRegularHidden+  :: Ord label+  => Map label (RegularSiteDisposition label)+  -> Map label (RegularSiteDisposition label)+  -> label+  -> Bool+becameRegularHidden before after label =+  Map.lookup label before == Just RegularSiteVisible+    && Map.lookup label after == Just RegularSiteHidden++requireRegularEdgeIncidence+  :: Ord label+  => RegularSection label+  -> RegularEdgeKey label+  -> Either+      (LocalRegularEditObstruction label)+      (Set (RegularFaceKey label))+requireRegularEdgeIncidence section edgeKey =+  maybe+    (Left (LocalRegularNonManifoldEdge edgeKey 0))+    (Right . sectionRegularEdgeFaces)+    (Map.lookup edgeKey (sectionRegularEdges section))++regularFaceGenerators+  :: Ord label+  => Map label (ExactPowerGenerator label)+  -> RegularFace label+  -> Either+      (LocalRegularEditObstruction label)+      ( ExactPowerGenerator label+      , ExactPowerGenerator label+      , ExactPowerGenerator label+      )+regularFaceGenerators generators face = do+  let (firstLabel, secondLabel, thirdLabel) = regularFaceLabels face+  (,,)+    <$> requireGeneratorFromMap generators firstLabel+    <*> requireGeneratorFromMap generators secondLabel+    <*> requireGeneratorFromMap generators thirdLabel++removeRegularSiteLocally+  :: Ord label+  => label+  -> RegularTriangulation label+  -> Either+      (LocalRegularEditObstruction label)+      (RegularTriangulation label, Set label)+removeRegularSiteLocally label triangulation = do+  section <-+    maybe (Left LocalRegularSectionUnavailable) Right+      (storedRegularSection triangulation)+  locality <-+    maybe (Left LocalRegularLocalityUnavailable) Right+      (sectionRegularLocality section)+  disposition <-+    maybe (Left (LocalRegularGeneratorMissing label)) Right+      (Map.lookup label (sectionRegularDispositions section))+  when+    (disposition /= RegularSiteVisible || hasCoincidentAliases label section)+    (Left (LocalRegularBoundaryRemoval label))+  star <-+    maybe (Left (LocalRegularGeneratorMissing label)) Right+      (Map.lookup label (sectionRegularStars section))+  let removedFaces = sectionIncidentFaces star+  when (Set.null removedFaces) (Left (LocalRegularBoundaryRemoval label))+  Foldable.traverse_+    (requireInteriorRemovalEdge section label)+    (sectionSiteNeighbours star)+  horizon <- regularConflictHorizon section removedFaces+  when (Map.null horizon) (Left LocalRegularEmptyHorizon)+  let hiddenCandidates =+        Set.unions+          [ Map.findWithDefault Set.empty faceKey (localitySupportHidden locality)+          | faceKey <- Set.toAscList removedFaces+          ]+      horizonLabels =+        Set.fromList+          [ endpoint+          | (firstLabel, secondLabel) <- Map.keys horizon+          , endpoint <- [firstLabel, secondLabel]+          ]+      candidateLabels = Set.union horizonLabels hiddenCandidates+      generators = Map.delete label (sectionGenerators section)+  candidateGenerators <-+    traverse+      (requireGeneratorFromMap generators)+      (Set.toAscList candidateLabels)+  nonEmptyCandidates <-+    maybe (Left LocalRegularCandidateSectionEmpty) Right+      (NonEmpty.nonEmpty candidateGenerators)+  candidateTopology <-+    first LocalRegularTopologyObstruction+      (regularGeneratorTopology (DistinctSlopeGenerators nonEmptyCandidates))+  when+    ( any+        ((== RegularGeneratorLowerDimensional) . snd)+        (NonEmpty.toList (generatorRegularDispositions candidateTopology))+    )+    (Left LocalRegularCavityBoundaryMismatch)+  removedFaceValues <-+    traverse+      (requireFaceFromMap (sectionRegularFaces section))+      (Set.toAscList removedFaces)+  let proposedFaces = fmap publishGeneratorFace (generatorRegularFaces candidateTopology)+      newFaces =+        filter+          (regularFaceCentroidInside generators removedFaceValues)+          proposedFaces+      proposedBoundary = regularFaceSectionBoundary newFaces+  when+    (proposedBoundary /= Map.keysSet horizon)+    (Left LocalRegularCavityBoundaryMismatch)+  facePatch <- applyRegularFacePatch generators removedFaces newFaces section+  certifyRegularFacePatch generators facePatch+  let candidateDispositions =+        Set.foldl'+          (publishRemovalCandidateDisposition (patchedRegularStars facePatch))+          (Map.delete label (sectionRegularDispositions section))+          candidateLabels+      exposed =+        Set.filter+          (becameRegularVisible (sectionRegularDispositions section) candidateDispositions)+          candidateLabels+      stillHidden =+        Set.filter+          (\candidate -> Map.lookup candidate candidateDispositions == Just RegularSiteHidden)+          hiddenCandidates+      receipt = sectionRegularReceipt section+      provisionalSection =+        section+          { sectionGenerators = generators+          , sectionSlopeRepresentatives =+              Map.filter (/= label) (sectionSlopeRepresentatives section)+          , sectionRegularDispositions = candidateDispositions+          , sectionRegularFaces = patchedRegularFaces facePatch+          , sectionRegularEdges = patchedRegularEdges facePatch+          , sectionRegularStars = Map.delete label (patchedRegularStars facePatch)+          , sectionRegularReceipt =+              receipt+                { generatorRegularInputSites = generatorRegularInputSites receipt - 1+                , generatorRegularVisibleSites =+                    generatorRegularVisibleSites receipt - 1 + Set.size exposed+                , generatorRegularHiddenSites =+                    generatorRegularHiddenSites receipt - Set.size exposed+                , generatorRegularFaceCount = Map.size (patchedRegularFaces facePatch)+                , generatorRegularEdgeCount = Map.size (patchedRegularEdges facePatch)+                }+          , sectionRegularLocality = Nothing+          }+  editedLocality <-+    patchRegularLocality locality removedFaces stillHidden provisionalSection+  let editedSection = provisionalSection{sectionRegularLocality = Just editedLocality}+  pure+    ( replaceRegularSiteSection+        (Map.delete label (storedRegularSites triangulation))+        editedSection+        triangulation+    , Set.insert label exposed+    )++requireInteriorRemovalEdge+  :: Ord label+  => RegularSection label+  -> label+  -> label+  -> Either (LocalRegularEditObstruction label) ()+requireInteriorRemovalEdge section label neighbour = do+  let edgeKey = orderedPair label neighbour+  incidence <- requireRegularEdgeIncidence section edgeKey+  if Set.size incidence == 2+    then Right ()+    else Left (LocalRegularBoundaryRemoval label)++regularFaceCentroidInside+  :: Ord label+  => Map label (ExactPowerGenerator label)+  -> [RegularFace label]+  -> RegularFace label+  -> Bool+regularFaceCentroidInside generators domainFaces candidate =+  case regularFaceSlopeSum generators candidate of+    Left _ -> False+    Right slopeSum -> any (regularSlopeSumInsideFace generators slopeSum) domainFaces++regularFaceSlopeSum+  :: Ord label+  => Map label (ExactPowerGenerator label)+  -> RegularFace label+  -> Either+      (LocalRegularEditObstruction label)+      (ExactRational, ExactRational)+regularFaceSlopeSum generators face = do+  (firstGenerator, secondGenerator, thirdGenerator) <-+    regularFaceGenerators generators face+  let (firstX, firstY) = exactGeneratorSlope firstGenerator+      (secondX, secondY) = exactGeneratorSlope secondGenerator+      (thirdX, thirdY) = exactGeneratorSlope thirdGenerator+  pure (firstX + secondX + thirdX, firstY + secondY + thirdY)++regularSlopeSumInsideFace+  :: Ord label+  => Map label (ExactPowerGenerator label)+  -> (ExactRational, ExactRational)+  -> RegularFace label+  -> Bool+regularSlopeSumInsideFace generators slopeSum =+  not+    . any (generatorEdgeExcludesScaled 3 slopeSum generators)+    . regularFaceDirectedEdgeKeys++regularFaceSectionBoundary+  :: Ord label+  => [RegularFace label]+  -> Set (RegularEdgeKey label)+regularFaceSectionBoundary =+  Map.keysSet+    . Map.filter ((== 1) . Set.size)+    . List.foldl'+        (Map.unionWith Set.union)+        Map.empty+    . fmap regularFaceEdgeIncidence++publishRemovalCandidateDisposition+  :: Ord label+  => Map label (RegularSiteStar label)+  -> Map label (RegularSiteDisposition label)+  -> label+  -> Map label (RegularSiteDisposition label)+publishRemovalCandidateDisposition stars dispositions label =+  Map.insert+    label+    ( if maybe True (Set.null . sectionIncidentFaces) (Map.lookup label stars)+        then RegularSiteHidden+        else RegularSiteVisible+    )+    dispositions++becameRegularVisible+  :: Ord label+  => Map label (RegularSiteDisposition label)+  -> Map label (RegularSiteDisposition label)+  -> label+  -> Bool+becameRegularVisible before after label =+  Map.lookup label before == Just RegularSiteHidden+    && Map.lookup label after == Just RegularSiteVisible++reweightRegularSiteLocally+  :: Ord label+  => PowerSite label+  -> RegularTriangulation label+  -> Either+      (LocalRegularEditObstruction label)+      (RegularTriangulation label, Set label)+reweightRegularSiteLocally revised triangulation = do+  let label = powerSiteLabel revised+  (removed, removalSupport) <- removeAnyRegularSiteLocally label triangulation+  (inserted, insertionSupport) <- insertRegularSiteLocally revised removed+  pure (inserted, Set.union removalSupport insertionSupport)++commonRegularWeightShift+  :: Ord label+  => Map label PowerWeight+  -> Map label (PowerSite label)+  -> Maybe ExactRational+commonRegularWeightShift requested sites+  | Map.size requested /= Map.size sites = Nothing+  | otherwise = do+      (firstLabel, firstWeight) <- Map.lookupMin requested+      firstSite <- Map.lookup firstLabel sites+      let shift =+            powerWeightExact firstWeight+              - powerWeightExact (powerSiteWeight firstSite)+      if Map.isSubmapOfBy (powerSiteHasWeightShift shift) requested sites+        then Just shift+        else Nothing++powerSiteHasWeightShift+  :: ExactRational+  -> PowerWeight+  -> PowerSite label+  -> Bool+powerSiteHasWeightShift shift weight site =+  powerWeightExact weight+    - powerWeightExact (powerSiteWeight site)+    == shift++translateRegularWeights+  :: ExactRational+  -> Map label (PowerSite label)+  -> Set label+  -> RegularTriangulation label+  -> RegularEditResult label+translateRegularWeights shift sites changedLabels triangulation =+  let shiftedSection =+        fmap+          (\section ->+             section+               { sectionGenerators =+                   fmap+                     (translateExactPowerGenerator shift)+                     (sectionGenerators section)+               })+          (storedRegularSection triangulation)+   in RegularEditResult+        { regularEditTriangulation =+            triangulation+              { storedRegularSites = sites+              , storedRegularSection = shiftedSection+              }+        , regularEditChangedSites = changedLabels+        , regularEditTransitions = Vector.empty+        }++translateExactPowerGenerator+  :: ExactRational+  -> ExactPowerGenerator label+  -> ExactPowerGenerator label+translateExactPowerGenerator shift generator =+  generator+    { exactPowerGeneratorConstant =+        exactPowerGeneratorConstant generator + shift+    }++removeAnyRegularSiteLocally+  :: Ord label+  => label+  -> RegularTriangulation label+  -> Either+      (LocalRegularEditObstruction label)+      (RegularTriangulation label, Set label)+removeAnyRegularSiteLocally label triangulation =+  case removeTopologyPreservingSite label triangulation of+    Just (edited, _) -> Right (edited, Set.singleton label)+    Nothing -> removeRegularSiteLocally label triangulation++-- Exact topology-preserving edits descend through the existing regular+-- section. A failed classification is not an error; it selects the canonical+-- global normalization below.++classifyCoincidentSubordinate+  :: Ord label+  => PowerSite label+  -> PowerSite label+  -> Maybe (CoincidentGeneratorDisposition label)+classifyCoincidentSubordinate representative candidate =+  if powerSiteQueryPoint candidate /= powerSiteQueryPoint representative+    then Nothing+    else+      case+        compareCoincidentPriority+          (powerSiteWeight candidate, powerSiteLabel candidate)+          (powerSiteWeight representative, powerSiteLabel representative) of+        LT ->+          Just+            ( classifyCoincidentValue+                (powerSiteLabel representative)+                (powerSiteWeight representative)+                (powerSiteWeight candidate)+            )+        _ -> Nothing++removeTopologyPreservingSite+  :: Ord label+  => label+  -> RegularTriangulation label+  -> Maybe+      ( RegularTriangulation label+      , Vector.Vector (RegularSiteTransition label)+      )+removeTopologyPreservingSite label triangulation = do+  section <- storedRegularSection triangulation+  disposition <- Map.lookup label (sectionRegularDispositions section)+  editedSection <-+    case disposition of+      RegularSiteCoincidentEquivalentTo _ -> removeCoincident section+      RegularSiteCoincidentDominatedBy _ -> removeCoincident section+      RegularSiteHidden+        | not (hasCoincidentAliases label section) -> do+            generator <- Map.lookup label (sectionGenerators section)+            let receipt = sectionRegularReceipt section+            pure+              section+                { sectionGenerators = Map.delete label (sectionGenerators section)+                , sectionSlopeRepresentatives =+                    Map.delete (exactGeneratorSlope generator) (sectionSlopeRepresentatives section)+                , sectionRegularDispositions =+                    Map.delete label (sectionRegularDispositions section)+                , sectionRegularStars = Map.delete label (sectionRegularStars section)+                , sectionRegularReceipt =+                    receipt+                      { generatorRegularInputSites = generatorRegularInputSites receipt - 1+                      , generatorRegularHiddenSites = generatorRegularHiddenSites receipt - 1+                      }+                , sectionRegularLocality = fmap (detachRegularHidden label) (sectionRegularLocality section)+                }+      _ -> Nothing+  let edited =+        replaceRegularSiteSection+          (Map.delete label (storedRegularSites triangulation))+          editedSection+          triangulation+  pure+    ( edited+    , Vector.singleton (RegularSiteDisappeared label disposition)+    )+ where+  removeCoincident section = do+    _ <- Map.lookup label (sectionCoincidentDispositions section)+    pure+      section+        { sectionCoincidentDispositions =+            Map.delete label (sectionCoincidentDispositions section)+        , sectionRegularDispositions =+            Map.delete label (sectionRegularDispositions section)+        , sectionRegularStars = Map.delete label (sectionRegularStars section)+        }++data RegularReweightDescent label = RegularReweightDescent+  { descendedCoincidentDispositions :: !(Map label (CoincidentGeneratorDisposition label))+  , descendedGenerators :: !(Map label (ExactPowerGenerator label))+  , descendedRegularDispositions :: !(Map label (RegularSiteDisposition label))+  , descendedTransitions :: ![RegularSiteTransition label]+  }++reweightTopologyPreservingSites+  :: Ord label+  => Set label+  -> Map label (PowerSite label)+  -> Map label (PowerSite label)+  -> RegularTriangulation label+  -> Maybe+      ( RegularTriangulation label+      , Vector.Vector (RegularSiteTransition label)+      )+reweightTopologyPreservingSites changedLabels changedSites reweightedSites triangulation = do+  section <- storedRegularSection triangulation+  let coincidentRepresentatives =+        Map.foldl'+          (flip (Set.insert . coincidentRepresentativeLabel))+          Set.empty+          (sectionCoincidentDispositions section)+  descent <-+    Foldable.foldlM+      ( descendTopologyPreservingReweight+          changedLabels+          reweightedSites+          coincidentRepresentatives+          triangulation+      )+      RegularReweightDescent+        { descendedCoincidentDispositions = sectionCoincidentDispositions section+        , descendedGenerators = sectionGenerators section+        , descendedRegularDispositions = sectionRegularDispositions section+        , descendedTransitions = []+        }+      changedSites+  let editedSection =+        section+          { sectionCoincidentDispositions = descendedCoincidentDispositions descent+          , sectionGenerators = descendedGenerators descent+          , sectionRegularDispositions = descendedRegularDispositions descent+          }+      edited =+        replaceRegularSiteSection+          reweightedSites+          editedSection+          triangulation+  pure (edited, Vector.fromList (reverse (descendedTransitions descent)))++descendTopologyPreservingReweight+  :: Ord label+  => Set label+  -> Map label (PowerSite label)+  -> Set label+  -> RegularTriangulation label+  -> RegularReweightDescent label+  -> PowerSite label+  -> Maybe (RegularReweightDescent label)+descendTopologyPreservingReweight changedLabels reweightedSites coincidentRepresentatives triangulation descent revised = do+  resident <- Map.lookup label (storedRegularSites triangulation)+  disposition <-+    storedRegularSection triangulation+      >>= Map.lookup label . sectionRegularDispositions+  case regularCoincidentRepresentative disposition of+    Just representative+      | Set.notMember representative changedLabels -> do+          representativeSite <- Map.lookup representative reweightedSites+          coincident <- classifyCoincidentSubordinate representativeSite revised+          let after = coincidentRegularDisposition coincident+          pure+            descent+              { descendedCoincidentDispositions =+                  Map.insert label coincident (descendedCoincidentDispositions descent)+              , descendedRegularDispositions =+                  Map.insert label after (descendedRegularDispositions descent)+              , descendedTransitions =+                  if disposition == after+                    then descendedTransitions descent+                    else RegularSiteTransitioned label disposition after : descendedTransitions descent+              }+    _+      | disposition == RegularSiteHidden+          && powerSiteWeight revised < powerSiteWeight resident+          && Set.notMember label coincidentRepresentatives ->+          let (generator, _) = prepareExactPowerGenerator revised+           in Just+                descent+                  { descendedGenerators =+                      Map.insert label generator (descendedGenerators descent)+                  }+      | otherwise -> Nothing+ where+  label = powerSiteLabel revised++hasCoincidentAliases+  :: Eq label+  => label+  -> RegularSection label+  -> Bool+hasCoincidentAliases label =+  any ((== label) . coincidentRepresentativeLabel)+    . Map.elems+    . sectionCoincidentDispositions++regularCoincidentRepresentative :: RegularSiteDisposition label -> Maybe label+regularCoincidentRepresentative disposition =+  case disposition of+    RegularSiteCoincidentEquivalentTo label -> Just label+    RegularSiteCoincidentDominatedBy label -> Just label+    _ -> Nothing++coincidentRepresentativeLabel :: CoincidentGeneratorDisposition label -> label+coincidentRepresentativeLabel disposition =+  case disposition of+    CoincidentGeneratorEquivalentTo label -> label+    CoincidentGeneratorDominatedBy label -> label++normalizeRegularEdit+  :: Ord label+  => RegularTriangulation label+  -> Map label (PowerSite label)+  -> Set label+  -> Either+      (RegularEditError label)+      (RegularEditResult label)+normalizeRegularEdit previous sites changedSites = do+  triangulation <-+    first RegularEditTopologyFailed+      (constructRegularTriangulation sites)+  pure+    ( RegularEditResult+        triangulation+        changedSites+        ( dispositionTransitions+            (regularDispositionSection previous)+            (regularDispositionSection triangulation)+        )+    )++unchangedRegularEdit+  :: RegularTriangulation label+  -> RegularEditResult label+unchangedRegularEdit triangulation =+  RegularEditResult triangulation Set.empty Vector.empty+++dispositionTransitions+  :: Ord label+  => Map label (RegularSiteDisposition label)+  -> Map label (RegularSiteDisposition label)+  -> Vector.Vector (RegularSiteTransition label)+dispositionTransitions previous current =+  Vector.fromList+    ( Map.elems+        ( Map.mergeWithKey+            changedDispositionTransition+            (Map.mapWithKey RegularSiteDisappeared)+            (Map.mapWithKey RegularSiteAppeared)+            previous+            current+        )+    )++dispositionTransitionsOn+  :: Ord label+  => Set label+  -> Map label (RegularSiteDisposition label)+  -> Map label (RegularSiteDisposition label)+  -> Vector.Vector (RegularSiteTransition label)+dispositionTransitionsOn support previous current =+  dispositionTransitions+    (Map.restrictKeys previous support)+    (Map.restrictKeys current support)++changedDispositionTransition+  :: Eq label+  => label+  -> RegularSiteDisposition label+  -> RegularSiteDisposition label+  -> Maybe (RegularSiteTransition label)+changedDispositionTransition label before after+  | before == after = Nothing+  | otherwise = Just (RegularSiteTransitioned label before after)

@@ -0,0 +1,264 @@+{-# LANGUAGE DeriveAnyClass #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DerivingStrategies #-}++-- | Exact affine generators and their dual geometry. This module owns the+-- algebra shared by batch hull construction, local edits, and bounded+-- projection; it does not choose a construction schedule.+module Moonlight.Triangulation.Internal.PowerDiagram.Generator+  ( ExactPowerGenerator (..)+  , DistinctSlopeGenerators (..)+  , RegularGeneratorDisposition (..)+  , GeneratorRegularFace (..)+  , GeneratorDualGeometry (..)+  , GeneratorRegularEdge (..)+  , GeneratorRegularTopology (..)+  , GeneratorRegularReceipt (..)+  , RegularTopologyError (..)+  , exactGeneratorAxis+  , exactGeneratorLiftedVolume+  , exactGeneratorFaceDual+  , exactGeneratorBoundaryDualRay+  , exactGeneratorDualBetween+  ) where++import Control.DeepSeq (NFData)+import Data.Bifunctor (first)+import Data.List.NonEmpty (NonEmpty)+import GHC.Generics (Generic)+import Moonlight.Triangulation.Exact+  ( ExactAffineLine+  , ExactGeometryError+  , ExactHalfPlaneError+  , ExactIntersectionError+  , ExactPoint+  , ExactRay+  , ExactSegment+  , ExactVector (..)+  , exactAffineLine+  , exactAffineLineCoefficients+  , exactAffineLineIntersection+  , exactRay+  , exactSegment+  )+import Moonlight.Triangulation.Internal.ExactRational (ExactRational)+import Moonlight.Triangulation.Minkowski (MinkowskiError)+data ExactPowerGenerator label = ExactPowerGenerator+  { exactPowerGeneratorLabel :: !label+  , exactPowerGeneratorXCoefficient :: !ExactRational+  , exactPowerGeneratorYCoefficient :: !ExactRational+  , exactPowerGeneratorConstant :: !ExactRational+  }+  deriving stock (Eq, Ord, Show, Generic)+  deriving anyclass (NFData)++-- | Package-internal witness that coincident affine slopes were resolved once.+newtype DistinctSlopeGenerators label =+  DistinctSlopeGenerators (NonEmpty (ExactPowerGenerator label))++data RegularGeneratorDisposition+  = RegularGeneratorVisible+  | RegularGeneratorLowerDimensional+  | RegularGeneratorHidden+  deriving stock (Eq, Ord, Show, Generic)+  deriving anyclass (NFData)++data GeneratorRegularFace label = GeneratorRegularFace+  { generatorRegularFaceFirst :: !label+  , generatorRegularFaceSecond :: !label+  , generatorRegularFaceThird :: !label+  , generatorRegularFaceDualPoint :: !ExactPoint+  }+  deriving stock (Eq, Ord, Show, Generic)+  deriving anyclass (NFData)++data GeneratorDualGeometry+  = GeneratorDualSegment !ExactSegment+  | GeneratorDualRay !ExactRay+  | GeneratorDualLine !ExactAffineLine+  | GeneratorDualCollapsed !ExactPoint+  deriving stock (Eq, Ord, Show, Generic)+  deriving anyclass (NFData)++data GeneratorRegularEdge label = GeneratorRegularEdge+  { generatorRegularEdgeFirst :: !(ExactPowerGenerator label)+  , generatorRegularEdgeSecond :: !(ExactPowerGenerator label)+  , generatorRegularEdgeDual :: !GeneratorDualGeometry+  }+  deriving stock (Eq, Ord, Show, Generic)+  deriving anyclass (NFData)++data GeneratorRegularTopology label = GeneratorRegularTopology+  { generatorRegularDispositions :: !(NonEmpty (ExactPowerGenerator label, RegularGeneratorDisposition))+  , generatorRegularFaces :: ![GeneratorRegularFace label]+  , generatorRegularEdges :: ![GeneratorRegularEdge label]+  , generatorRegularReceipt :: !GeneratorRegularReceipt+  }+  deriving stock (Eq, Show, Generic)+  deriving anyclass (NFData)++data GeneratorRegularReceipt = GeneratorRegularReceipt+  { generatorRegularInputSites :: !Int+  , generatorRegularVisibleSites :: !Int+  , generatorRegularLowerDimensionalSites :: !Int+  , generatorRegularHiddenSites :: !Int+  , generatorRegularFaceCount :: !Int+  , generatorRegularEdgeCount :: !Int+  }+  deriving stock (Eq, Ord, Show, Generic)+  deriving anyclass (NFData)++-- | Exact slope-hull, incidence, and dual-construction obstructions.+data RegularTopologyError label+  = RegularTopologySlopeHullInvalid !MinkowskiError+  | RegularTopologySlopeVertexMissing !ExactPoint+  | RegularTopologyAxisInvalid !label !label !ExactHalfPlaneError+  | RegularTopologyDualVertexInvalid !label !label !label !ExactIntersectionError+  | RegularTopologyDualSegmentInvalid !label !label !ExactGeometryError+  | RegularTopologyDualRayInvalid !label !label !ExactGeometryError+  | RegularTopologyFlatBoundary !label !label !label+  | RegularTopologyNonManifoldEdge !label !label !Int+  deriving stock (Eq, Show, Generic)+  deriving anyclass (NFData)++exactGeneratorAxis+  :: ExactPowerGenerator label+  -> ExactPowerGenerator label+  -> Either ExactHalfPlaneError ExactAffineLine+exactGeneratorAxis firstGenerator secondGenerator =+  exactAffineLine+    ( exactPowerGeneratorXCoefficient firstGenerator+        - exactPowerGeneratorXCoefficient secondGenerator+    )+    ( exactPowerGeneratorYCoefficient firstGenerator+        - exactPowerGeneratorYCoefficient secondGenerator+    )+    ( exactPowerGeneratorConstant firstGenerator+        - exactPowerGeneratorConstant secondGenerator+    )++-- | Signed exact lifted volume of four affine generators. For a+-- counter-clockwise slope face, a positive result lies above its supporting+-- plane and therefore belongs to the insertion conflict cavity.+exactGeneratorLiftedVolume+  :: ExactPowerGenerator label+  -> ExactPowerGenerator label+  -> ExactPowerGenerator label+  -> ExactPowerGenerator label+  -> ExactRational+exactGeneratorLiftedVolume firstGenerator secondGenerator thirdGenerator fourthGenerator =+  let (firstX, firstY, firstZ) = generatorCoordinates firstGenerator+      (secondX, secondY, secondZ) = generatorCoordinates secondGenerator+      (thirdX, thirdY, thirdZ) = generatorCoordinates thirdGenerator+      (fourthX, fourthY, fourthZ) = generatorCoordinates fourthGenerator+      secondX' = secondX - firstX+      secondY' = secondY - firstY+      secondZ' = secondZ - firstZ+      thirdX' = thirdX - firstX+      thirdY' = thirdY - firstY+      thirdZ' = thirdZ - firstZ+      fourthX' = fourthX - firstX+      fourthY' = fourthY - firstY+      fourthZ' = fourthZ - firstZ+   in secondX' * (thirdY' * fourthZ' - thirdZ' * fourthY')+        - secondY' * (thirdX' * fourthZ' - thirdZ' * fourthX')+        + secondZ' * (thirdX' * fourthY' - thirdY' * fourthX')++-- | Exact dual vertex shared by three affine generators.+exactGeneratorFaceDual+  :: ExactPowerGenerator label+  -> ExactPowerGenerator label+  -> ExactPowerGenerator label+  -> Either (RegularTopologyError label) ExactPoint+exactGeneratorFaceDual firstGenerator secondGenerator thirdGenerator = do+  firstAxis <-+    first+      (RegularTopologyAxisInvalid (exactPowerGeneratorLabel firstGenerator) (exactPowerGeneratorLabel secondGenerator))+      (exactGeneratorAxis firstGenerator secondGenerator)+  secondAxis <-+    first+      (RegularTopologyAxisInvalid (exactPowerGeneratorLabel firstGenerator) (exactPowerGeneratorLabel thirdGenerator))+      (exactGeneratorAxis firstGenerator thirdGenerator)+  first+    ( RegularTopologyDualVertexInvalid+        (exactPowerGeneratorLabel firstGenerator)+        (exactPowerGeneratorLabel secondGenerator)+        (exactPowerGeneratorLabel thirdGenerator)+    )+    (exactAffineLineIntersection firstAxis secondAxis)++-- | Exact outward dual ray of an oriented boundary edge and its incident+-- counter-clockwise face.+exactGeneratorBoundaryDualRay+  :: ExactPoint+  -> ExactPowerGenerator label+  -> ExactPowerGenerator label+  -> ExactPowerGenerator label+  -> Either (RegularTopologyError label) ExactRay+exactGeneratorBoundaryDualRay originPoint firstGenerator secondGenerator thirdGenerator = do+  axis <-+    first+      ( RegularTopologyAxisInvalid+          (exactPowerGeneratorLabel firstGenerator)+          (exactPowerGeneratorLabel secondGenerator)+      )+      (exactGeneratorAxis firstGenerator secondGenerator)+  let (axisX, axisY, _) = exactAffineLineCoefficients axis+      candidate = ExactVector axisY (negate axisX)+      advantageX =+        exactPowerGeneratorXCoefficient firstGenerator+          - exactPowerGeneratorXCoefficient thirdGenerator+      advantageY =+        exactPowerGeneratorYCoefficient firstGenerator+          - exactPowerGeneratorYCoefficient thirdGenerator+      ExactVector candidateX candidateY = candidate+      derivative = advantageX * candidateX + advantageY * candidateY+      direction =+        if derivative > 0+          then candidate+          else ExactVector (negate candidateX) (negate candidateY)+  if derivative == 0+    then+      Left+        ( RegularTopologyFlatBoundary+            (exactPowerGeneratorLabel firstGenerator)+            (exactPowerGeneratorLabel secondGenerator)+            (exactPowerGeneratorLabel thirdGenerator)+        )+    else+      first+        ( RegularTopologyDualRayInvalid+            (exactPowerGeneratorLabel firstGenerator)+            (exactPowerGeneratorLabel secondGenerator)+        )+        (exactRay originPoint direction)++-- | Canonical undirected dual geometry between two incident face vertices.+exactGeneratorDualBetween+  :: label+  -> label+  -> ExactPoint+  -> ExactPoint+  -> Either (RegularTopologyError label) GeneratorDualGeometry+exactGeneratorDualBetween firstLabel secondLabel firstPoint secondPoint+  | firstPoint == secondPoint = Right (GeneratorDualCollapsed firstPoint)+  | otherwise =+      fmap GeneratorDualSegment+        ( first+            (RegularTopologyDualSegmentInvalid firstLabel secondLabel)+            ( uncurry exactSegment+                ( if firstPoint <= secondPoint+                    then (firstPoint, secondPoint)+                    else (secondPoint, firstPoint)+                )+            )+        )++generatorCoordinates+  :: ExactPowerGenerator label+  -> (ExactRational, ExactRational, ExactRational)+generatorCoordinates generator =+  ( exactPowerGeneratorXCoefficient generator+  , exactPowerGeneratorYCoefficient generator+  , exactPowerGeneratorConstant generator+  )

@@ -0,0 +1,925 @@+-- | Exact upper-hull descent for affine generators. The public power owner+-- supplies admitted, distinct slopes; this module returns only topology and+-- exact dual geometry.+module Moonlight.Triangulation.Internal.PowerDiagram.Hull+  ( regularGeneratorTopology+  ) where++import Data.Bifunctor (first)+import qualified Data.Foldable as Foldable+import qualified Data.IntMap.Strict as IntMap+import qualified Data.IntSet as IntSet+import qualified Data.List as List+import Data.List.NonEmpty (NonEmpty (..))+import qualified Data.List.NonEmpty as NonEmpty+import Data.Map.Strict (Map)+import qualified Data.Map.Strict as Map+import Data.Maybe (mapMaybe)+import Data.Set (Set)+import qualified Data.Set as Set+import Moonlight.Triangulation.Exact+  ( ExactAffineLine+  , ExactPoint+  , ExactRay+  , exactPoint+  , exactPointCoordinates+  )+import Moonlight.Triangulation.Internal.BoundaryCycle (consecutivePairs)+import Moonlight.Triangulation.Internal.ExactRational (ExactRational)+import Moonlight.Triangulation.Internal.Minkowski.Convex (convexHullPolygon)+import Moonlight.Triangulation.Internal.PowerDiagram.Generator+import Moonlight.Triangulation.Minkowski (convexPolygonPoints)++-- | Construct exact regular topology from distinct affine slopes. Full-rank+-- input uses one conflict graph; collinear slopes use the one-dimensional+-- upper chain; coplanar lifted input reuses the exact planar hull owner.+regularGeneratorTopology+  :: Ord label+  => DistinctSlopeGenerators label+  -> Either (RegularTopologyError label) (GeneratorRegularTopology label)+regularGeneratorTopology (DistinctSlopeGenerators generators) =+  let points = admittedLiftedPoints generators+   in case points of+        singleton :| [] ->+          pure (topologyFromParts points (Set.singleton singleton) Set.empty [] [])+        firstPoint :| secondPoint : remaining ->+          case List.find (not . slopeCollinear firstPoint secondPoint) remaining of+            Nothing -> collinearSlopeTopology points firstPoint secondPoint+            Just thirdPoint ->+              case List.find (not . liftedCoplanar firstPoint secondPoint thirdPoint) remaining of+                Nothing -> coplanarLiftedTopology points+                Just fourthPoint -> fullRankTopology points firstPoint secondPoint thirdPoint fourthPoint++data LiftedPoint label = LiftedPoint+  { liftedPointId :: !Int+  , liftedPointGenerator :: !(ExactPowerGenerator label)+  }++instance Eq (LiftedPoint label) where+  left == right = liftedPointId left == liftedPointId right++instance Ord (LiftedPoint label) where+  compare left right = compare (liftedPointId left) (liftedPointId right)++admittedLiftedPoints+  :: Ord label+  => NonEmpty (ExactPowerGenerator label)+  -> NonEmpty (LiftedPoint label)+admittedLiftedPoints generators =+  let ordered = NonEmpty.sortWith exactPowerGeneratorLabel generators+   in NonEmpty.zipWith+        (\identifier generator -> LiftedPoint identifier generator)+        (0 :| [1 ..])+        ordered++liftedCoordinates+  :: LiftedPoint label+  -> (ExactRational, ExactRational, ExactRational)+liftedCoordinates point =+  let generator = liftedPointGenerator point+   in ( exactPowerGeneratorXCoefficient generator+      , exactPowerGeneratorYCoefficient generator+      , exactPowerGeneratorConstant generator+      )++slopePoint :: LiftedPoint label -> ExactPoint+slopePoint point =+  let (coordinateX, coordinateY, _) = liftedCoordinates point+   in exactPoint coordinateX coordinateY++slopeCollinear+  :: LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> Bool+slopeCollinear firstPoint secondPoint thirdPoint =+  slopeOrientation firstPoint secondPoint thirdPoint == 0++slopeOrientation+  :: LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> ExactRational+slopeOrientation firstPoint secondPoint thirdPoint =+  let (firstX, firstY, _) = liftedCoordinates firstPoint+      (secondX, secondY, _) = liftedCoordinates secondPoint+      (thirdX, thirdY, _) = liftedCoordinates thirdPoint+   in (secondX - firstX) * (thirdY - firstY)+        - (secondY - firstY) * (thirdX - firstX)++liftedCoplanar+  :: LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> Bool+liftedCoplanar firstPoint secondPoint thirdPoint fourthPoint =+  liftedVolume firstPoint secondPoint thirdPoint fourthPoint == 0++liftedVolume+  :: LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> ExactRational+liftedVolume firstPoint secondPoint thirdPoint fourthPoint =+  exactGeneratorLiftedVolume+    (liftedPointGenerator firstPoint)+    (liftedPointGenerator secondPoint)+    (liftedPointGenerator thirdPoint)+    (liftedPointGenerator fourthPoint)++data HullFace label = HullFace+  { hullFaceId :: !Int+  , hullFaceFirst :: !(LiftedPoint label)+  , hullFaceSecond :: !(LiftedPoint label)+  , hullFaceThird :: !(LiftedPoint label)+  , hullFaceOutside :: !(Set (LiftedPoint label))+  , hullFaceCoplanar :: !(Set (LiftedPoint label))+  }++type HullEdge = (Int, Int)++data HullState label = HullState+  { hullFaces :: !(IntMap.IntMap (HullFace label))+  , hullEdgeFaces :: !(Map HullEdge IntSet.IntSet)+  , hullConflictFaces :: !IntSet.IntSet+  , hullNextFaceId :: !Int+  }++fullRankTopology+  :: NonEmpty (LiftedPoint label)+  -> LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> Either (RegularTopologyError label) (GeneratorRegularTopology label)+fullRankTopology points firstPoint secondPoint thirdPoint fourthPoint = do+  let interior = liftedInteriorSum firstPoint secondPoint thirdPoint fourthPoint+      simplex = Set.fromList [firstPoint, secondPoint, thirdPoint, fourthPoint]+      initialFaces =+        zipWith+          (orientedHullFace interior)+          [0 ..]+          [ (firstPoint, secondPoint, thirdPoint)+          , (firstPoint, fourthPoint, secondPoint)+          , (firstPoint, thirdPoint, fourthPoint)+          , (secondPoint, fourthPoint, thirdPoint)+          ]+      initialState =+        List.foldl'+          insertHullFace+          (HullState IntMap.empty Map.empty IntSet.empty 4)+          initialFaces+      seededState =+        Foldable.foldl'+          (flip assignPointToHull)+          initialState+          (filter (`Set.notMember` simplex) (NonEmpty.toList points))+      finalState = descendHull interior seededState+      upperFaces = filter ((> 0) . hullFaceNormalZ) (IntMap.elems (hullFaces finalState))+  (canonicalFaces, visible, lowerDimensional) <- canonicalUpperFaces upperFaces+  (faces, edges) <- topologyFromHullFaces canonicalFaces+  pure+    (topologyFromParts points visible lowerDimensional faces edges)++liftedInteriorSum+  :: LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> (ExactRational, ExactRational, ExactRational)+liftedInteriorSum firstPoint secondPoint thirdPoint fourthPoint =+  let (firstX, firstY, firstZ) = liftedCoordinates firstPoint+      (secondX, secondY, secondZ) = liftedCoordinates secondPoint+      (thirdX, thirdY, thirdZ) = liftedCoordinates thirdPoint+      (fourthX, fourthY, fourthZ) = liftedCoordinates fourthPoint+   in ( firstX + secondX + thirdX + fourthX+      , firstY + secondY + thirdY + fourthY+      , firstZ + secondZ + thirdZ + fourthZ+      )++orientedHullFace+  :: (ExactRational, ExactRational, ExactRational)+  -> Int+  -> (LiftedPoint label, LiftedPoint label, LiftedPoint label)+  -> HullFace label+orientedHullFace interior identifier (firstPoint, secondPoint, thirdPoint) =+  let provisional = HullFace identifier firstPoint secondPoint thirdPoint Set.empty Set.empty+   in if hullFaceVolumeAtInterior interior provisional < 0+        then provisional+        else HullFace identifier firstPoint thirdPoint secondPoint Set.empty Set.empty++hullFaceVolumeAtInterior+  :: (ExactRational, ExactRational, ExactRational)+  -> HullFace label+  -> ExactRational+hullFaceVolumeAtInterior (sumX, sumY, sumZ) face =+  let (firstX, firstY, firstZ) = liftedCoordinates (hullFaceFirst face)+      (secondX, secondY, secondZ) = liftedCoordinates (hullFaceSecond face)+      (thirdX, thirdY, thirdZ) = liftedCoordinates (hullFaceThird face)+      secondX' = secondX - firstX+      secondY' = secondY - firstY+      secondZ' = secondZ - firstZ+      thirdX' = thirdX - firstX+      thirdY' = thirdY - firstY+      thirdZ' = thirdZ - firstZ+      interiorX' = sumX - 4 * firstX+      interiorY' = sumY - 4 * firstY+      interiorZ' = sumZ - 4 * firstZ+   in secondX' * (thirdY' * interiorZ' - thirdZ' * interiorY')+        - secondY' * (thirdX' * interiorZ' - thirdZ' * interiorX')+        + secondZ' * (thirdX' * interiorY' - thirdY' * interiorX')++liftedVolumeAtFace :: HullFace label -> LiftedPoint label -> ExactRational+liftedVolumeAtFace face =+  liftedVolume+    (hullFaceFirst face)+    (hullFaceSecond face)+    (hullFaceThird face)++hullFaceNormalZ :: HullFace label -> ExactRational+hullFaceNormalZ face =+  slopeOrientation+    (hullFaceFirst face)+    (hullFaceSecond face)+    (hullFaceThird face)++hullFacePoints :: HullFace label -> [LiftedPoint label]+hullFacePoints face =+  [hullFaceFirst face, hullFaceSecond face, hullFaceThird face]++hullFaceDirectedEdges+  :: HullFace label+  -> [(LiftedPoint label, LiftedPoint label)]+hullFaceDirectedEdges face =+  [ (hullFaceFirst face, hullFaceSecond face)+  , (hullFaceSecond face, hullFaceThird face)+  , (hullFaceThird face, hullFaceFirst face)+  ]++hullEdgeKey :: LiftedPoint label -> LiftedPoint label -> HullEdge+hullEdgeKey firstPoint secondPoint =+  let firstId = liftedPointId firstPoint+      secondId = liftedPointId secondPoint+   in if firstId <= secondId then (firstId, secondId) else (secondId, firstId)++insertHullFace :: HullState label -> HullFace label -> HullState label+insertHullFace state face =+  state+    { hullFaces = IntMap.insert (hullFaceId face) face (hullFaces state)+    , hullEdgeFaces =+        List.foldl'+          (\incidence (fromPoint, toPoint) ->+             Map.insertWith+               IntSet.union+               (hullEdgeKey fromPoint toPoint)+               (IntSet.singleton (hullFaceId face))+               incidence)+          (hullEdgeFaces state)+          (hullFaceDirectedEdges face)+    }++assignPointToHull :: LiftedPoint label -> HullState label -> HullState label+assignPointToHull point state =+  attachSelectedPoint point state+    (selectHullFace point (IntMap.elems (hullFaces state)))++data HullFaceSelection label+  = HullFaceInterior+  | HullFaceCoplanar !(HullFace label)+  | HullFaceVisible !(HullFace label) !ExactRational++selectHullFace+  :: LiftedPoint label+  -> [HullFace label]+  -> HullFaceSelection label+selectHullFace point = List.foldl' choose HullFaceInterior+ where+  choose selected face =+    let volume = liftedVolumeAtFace face point+     in case compare volume 0 of+          LT -> selected+          EQ ->+            case selected of+              HullFaceInterior -> HullFaceCoplanar face+              _ -> selected+          GT ->+            case selected of+              HullFaceVisible current currentVolume+                | volume < currentVolume+                    || (volume == currentVolume && hullFaceId current < hullFaceId face) ->+                    selected+              _ -> HullFaceVisible face volume++attachSelectedPoint+  :: LiftedPoint label+  -> HullState label+  -> HullFaceSelection label+  -> HullState label+attachSelectedPoint point state selection =+  case selection of+    HullFaceInterior -> state+    HullFaceCoplanar face ->+      state+        { hullFaces =+            IntMap.adjust+              (\selected -> selected{hullFaceCoplanar = Set.insert point (hullFaceCoplanar selected)})+              (hullFaceId face)+              (hullFaces state)+        }+    HullFaceVisible face _ ->+      state+        { hullFaces =+            IntMap.adjust+              (\selected -> selected{hullFaceOutside = Set.insert point (hullFaceOutside selected)})+              (hullFaceId face)+              (hullFaces state)+        , hullConflictFaces = IntSet.insert (hullFaceId face) (hullConflictFaces state)+        }++descendHull+  :: (ExactRational, ExactRational, ExactRational)+  -> HullState label+  -> HullState label+descendHull interior state =+  case nextHullExpansion state of+    Nothing -> state+    Just (seedFace, apex) ->+      let visible = visibleHullPatch state apex (IntSet.singleton (hullFaceId seedFace)) [hullFaceId seedFace]+          visibleFaces =+            mapMaybe (`IntMap.lookup` hullFaces state) (IntSet.toAscList visible)+          horizon = hullHorizon state visible visibleFaces+          candidates =+            Set.difference+              ( Set.delete apex+                  ( Set.unions+                      [ Set.unions+                          [ hullFaceOutside face+                          , hullFaceCoplanar face+                          , Set.fromList (hullFacePoints face)+                          ]+                      | face <- visibleFaces+                      ]+                  )+              )+              ( Set.fromList+                  [ point+                  | (firstPoint, secondPoint) <- Map.elems horizon+                  , point <- [firstPoint, secondPoint]+                  ]+              )+          withoutVisible = removeHullFaces visibleFaces state+          (withHorizon, newFaceIds) =+            Map.foldl'+              (insertHorizonFace interior apex)+              (withoutVisible, [])+              horizon+          redistributed =+            Foldable.foldl'+              (assignPointToFaces newFaceIds)+              withHorizon+              candidates+       in descendHull interior redistributed++nextHullExpansion+  :: HullState label+  -> Maybe (HullFace label, LiftedPoint label)+nextHullExpansion state = do+  identifier <- IntSet.lookupMin (hullConflictFaces state)+  face <- IntMap.lookup identifier (hullFaces state)+  apex <- farthestOutsidePoint face+  pure (face, apex)++farthestOutsidePoint :: HullFace label -> Maybe (LiftedPoint label)+farthestOutsidePoint face =+  case Set.minView (hullFaceOutside face) of+    Nothing -> Nothing+    Just (initial, remaining) ->+      Just+        ( fst+            ( Foldable.foldl'+                (\(selected, selectedVolume) candidate ->+                   let candidateVolume = liftedVolumeAtFace face candidate+                    in if candidateVolume > selectedVolume+                         then (candidate, candidateVolume)+                         else (selected, selectedVolume))+                (initial, liftedVolumeAtFace face initial)+                remaining+            )+        )++visibleHullPatch+  :: HullState label+  -> LiftedPoint label+  -> IntSet.IntSet+  -> [Int]+  -> IntSet.IntSet+visibleHullPatch _ _ visited [] = visited+visibleHullPatch state apex visited (identifier : remaining) =+  case IntMap.lookup identifier (hullFaces state) of+    Nothing -> visibleHullPatch state apex visited remaining+    Just face ->+      let neighbours =+            IntSet.unions+              [ Map.findWithDefault IntSet.empty (hullEdgeKey fromPoint toPoint) (hullEdgeFaces state)+              | (fromPoint, toPoint) <- hullFaceDirectedEdges face+              ]+          newlyVisible =+            IntSet.filter+              (\candidateId ->+                 IntSet.notMember candidateId visited+                   && maybe+                     False+                     ((> 0) . (`liftedVolumeAtFace` apex))+                     (IntMap.lookup candidateId (hullFaces state)))+              neighbours+       in visibleHullPatch+            state+            apex+            (IntSet.union visited newlyVisible)+            (IntSet.toList newlyVisible <> remaining)++hullHorizon+  :: HullState label+  -> IntSet.IntSet+  -> [HullFace label]+  -> Map HullEdge (LiftedPoint label, LiftedPoint label)+hullHorizon state visible visibleFaces =+  Map.fromList+    [ (hullEdgeKey fromPoint toPoint, (fromPoint, toPoint))+    | face <- visibleFaces+    , (fromPoint, toPoint) <- hullFaceDirectedEdges face+    , not+        ( IntSet.null+            ( IntSet.difference+                (Map.findWithDefault IntSet.empty (hullEdgeKey fromPoint toPoint) (hullEdgeFaces state))+                visible+            )+        )+    ]++removeHullFaces :: [HullFace label] -> HullState label -> HullState label+removeHullFaces removed state = List.foldl' removeHullFace state removed++removeHullFace :: HullState label -> HullFace label -> HullState label+removeHullFace state face =+  state+    { hullFaces = IntMap.delete identifier (hullFaces state)+    , hullEdgeFaces =+        List.foldl'+          (removeHullFaceFromEdge identifier)+          (hullEdgeFaces state)+          (hullFaceDirectedEdges face)+    , hullConflictFaces = IntSet.delete identifier (hullConflictFaces state)+    }+ where+  identifier = hullFaceId face++removeHullFaceFromEdge+  :: Int+  -> Map HullEdge IntSet.IntSet+  -> (LiftedPoint label, LiftedPoint label)+  -> Map HullEdge IntSet.IntSet+removeHullFaceFromEdge identifier incidence (fromPoint, toPoint) =+  Map.update+    (\incident ->+       let retained = IntSet.delete identifier incident+        in if IntSet.null retained then Nothing else Just retained)+    (hullEdgeKey fromPoint toPoint)+    incidence++insertHorizonFace+  :: (ExactRational, ExactRational, ExactRational)+  -> LiftedPoint label+  -> (HullState label, [Int])+  -> (LiftedPoint label, LiftedPoint label)+  -> (HullState label, [Int])+insertHorizonFace interior apex (state, identifiers) (firstPoint, secondPoint) =+  let identifier = hullNextFaceId state+      face = orientedHullFace interior identifier (firstPoint, secondPoint, apex)+      nextState =+        (insertHullFace state face)+          { hullNextFaceId = identifier + 1 }+   in (nextState, identifier : identifiers)++assignPointToFaces+  :: [Int]+  -> HullState label+  -> LiftedPoint label+  -> HullState label+assignPointToFaces identifiers state point =+  attachSelectedPoint point state+    (selectHullFace point (mapMaybe (`IntMap.lookup` hullFaces state) identifiers))++-- | Glue the triangular implementation facets on each exact supporting plane,+-- then recover only that plane's extreme slope vertices. This removes+-- simplex/horizon diagonals from the visibility decision while retaining a+-- deterministic triangulation for incidence and collapsed duals.+canonicalUpperFaces+  :: [HullFace label]+  -> Either+      (RegularTopologyError label)+      ( [(HullFace label, ExactPoint)]+      , Set (LiftedPoint label)+      , Set (LiftedPoint label)+      )+canonicalUpperFaces faces = do+  associations <- traverse upperFaceAssociation faces+  sections <-+    traverse canonicalUpperFaceGroup+      ( Map.elems+          ( List.foldl'+              (\groups (planeKey, dualFace) ->+                 Map.insertWith (<>) planeKey (dualFace :| []) groups)+              Map.empty+              associations+          )+      )+  let (sectionFaces, visible, lowerDimensional) = mconcat sections+      canonicalFaces =+        zipWith+          (\identifier (face, dualPoint) ->+             (face{hullFaceId = identifier}, dualPoint))+          [0 ..]+          sectionFaces+  pure (canonicalFaces, visible, lowerDimensional)++canonicalUpperFaceGroup+  :: NonEmpty (HullFace label, ExactPoint)+  -> Either+      (RegularTopologyError label)+      ( [(HullFace label, ExactPoint)]+      , Set (LiftedPoint label)+      , Set (LiftedPoint label)+      )+canonicalUpperFaceGroup ((face, dualPoint) :| [])+  | Set.null (hullFaceCoplanar face) =+      Right ([(face, dualPoint)], Set.fromList (hullFacePoints face), Set.empty)+canonicalUpperFaceGroup ((firstFace, dualPoint) :| remainingDualFaces) = do+  let groupFaces = firstFace : fmap fst remainingDualFaces+      groupPoints =+        Set.unions+          [ Set.union+              (Set.fromList (hullFacePoints face))+              (hullFaceCoplanar face)+          | face <- groupFaces+          ]+      initialPoint = hullFaceFirst firstFace+      points = initialPoint :| Set.toAscList (Set.delete initialPoint groupPoints)+      pointsBySlope =+        Map.fromList [(slopePoint point, point) | point <- NonEmpty.toList points]+  hull <-+    first RegularTopologySlopeHullInvalid+      (convexHullPolygon (fmap slopePoint points))+  visiblePoints <-+    traverse+      (\point ->+         maybe+           (Left (RegularTopologySlopeVertexMissing point))+           Right+           (Map.lookup point pointsBySlope))+      (NonEmpty.toList (convexPolygonPoints hull))+  let visible = Set.fromList visiblePoints+  pure+    ( fmap (\fanFace -> (fanFace, dualPoint)) (coplanarFanFaces visiblePoints)+    , visible+    , Set.difference groupPoints visible+    )++hullFacePlaneKey+  :: HullFace label+  -> ExactPoint+  -> (ExactRational, ExactRational, ExactRational)+hullFacePlaneKey face dualPoint =+  let (dualX, dualY) = exactPointCoordinates dualPoint+      generator = liftedPointGenerator (hullFaceFirst face)+      support =+        exactPowerGeneratorXCoefficient generator * dualX+          + exactPowerGeneratorYCoefficient generator * dualY+          + exactPowerGeneratorConstant generator+   in (dualX, dualY, support)++upperFaceAssociation+  :: HullFace label+  -> Either+      (RegularTopologyError label)+      ( (ExactRational, ExactRational, ExactRational)+      , (HullFace label, ExactPoint)+      )+upperFaceAssociation face = do+  (_, dualPoint) <- hullFaceDual face+  pure (hullFacePlaneKey face dualPoint, (face, dualPoint))++coplanarLiftedTopology+  :: NonEmpty (LiftedPoint label)+  -> Either (RegularTopologyError label) (GeneratorRegularTopology label)+coplanarLiftedTopology points = do+  hull <-+    first RegularTopologySlopeHullInvalid+      (convexHullPolygon (fmap slopePoint points))+  let pointsBySlope = Map.fromList [(slopePoint point, point) | point <- NonEmpty.toList points]+  hullPoints <-+    traverse+      (\point -> maybe (Left (RegularTopologySlopeVertexMissing point)) Right (Map.lookup point pointsBySlope))+      (NonEmpty.toList (convexPolygonPoints hull))+  let visible = Set.fromList hullPoints+      lowerDimensional = Set.difference (Set.fromList (NonEmpty.toList points)) visible+      hullFaces' = coplanarFanFaces hullPoints+  dualFaces <- coplanarFanDualFaces hullFaces'+  (faces, edges) <- topologyFromHullFaces dualFaces+  pure+    (topologyFromParts points visible lowerDimensional faces edges)++coplanarFanDualFaces+  :: [HullFace label]+  -> Either+      (RegularTopologyError label)+      [(HullFace label, ExactPoint)]+coplanarFanDualFaces [] = Right []+coplanarFanDualFaces faces@(firstFace : _) = do+  (_, dualPoint) <- hullFaceDual firstFace+  pure (fmap (\face -> (face, dualPoint)) faces)++coplanarFanFaces :: [LiftedPoint label] -> [HullFace label]+coplanarFanFaces (firstPoint : secondPoint : thirdPoint : remaining) =+  zipWith+    (\identifier (leftPoint, rightPoint) ->+       HullFace identifier firstPoint leftPoint rightPoint Set.empty Set.empty)+    [0 ..]+    (consecutivePairs (secondPoint : thirdPoint : remaining))+coplanarFanFaces _ = []++collinearSlopeTopology+  :: NonEmpty (LiftedPoint label)+  -> LiftedPoint label+  -> LiftedPoint label+  -> Either (RegularTopologyError label) (GeneratorRegularTopology label)+collinearSlopeTopology points firstPoint secondPoint = do+  let direction = slopeDifference firstPoint secondPoint+      ordered = List.sortOn (slopeParameter direction) (NonEmpty.toList points)+      visibleList = upperSlopeChain direction ordered+      visible = Set.fromList visibleList+      lowerDimensional =+        Set.difference+          (upperChainPointSet direction visibleList ordered)+          visible+  edges <- traverse fullLineEdge (consecutivePairs visibleList)+  pure (topologyFromParts points visible lowerDimensional [] edges)++slopeDifference+  :: LiftedPoint label+  -> LiftedPoint label+  -> (ExactRational, ExactRational)+slopeDifference firstPoint secondPoint =+  let (firstX, firstY, _) = liftedCoordinates firstPoint+      (secondX, secondY, _) = liftedCoordinates secondPoint+   in (secondX - firstX, secondY - firstY)++slopeParameter+  :: (ExactRational, ExactRational)+  -> LiftedPoint label+  -> ExactRational+slopeParameter (directionX, directionY) point =+  let (coordinateX, coordinateY, _) = liftedCoordinates point+   in directionX * coordinateX + directionY * coordinateY++upperSlopeChain+  :: (ExactRational, ExactRational)+  -> [LiftedPoint label]+  -> [LiftedPoint label]+upperSlopeChain direction = reverse . List.foldl' insertUpper []+ where+  insertUpper+    :: [LiftedPoint label]+    -> LiftedPoint label+    -> [LiftedPoint label]+  insertUpper (secondPoint : firstPoint : remaining) candidate+    | slopeChainTurn direction firstPoint secondPoint candidate >= 0 =+        insertUpper (firstPoint : remaining) candidate+  insertUpper retained candidate = candidate : retained++slopeChainTurn+  :: (ExactRational, ExactRational)+  -> LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> ExactRational+slopeChainTurn direction firstPoint secondPoint thirdPoint =+  let firstT = slopeParameter direction firstPoint+      secondT = slopeParameter direction secondPoint+      thirdT = slopeParameter direction thirdPoint+      (_, _, firstConstant) = liftedCoordinates firstPoint+      (_, _, secondConstant) = liftedCoordinates secondPoint+      (_, _, thirdConstant) = liftedCoordinates thirdPoint+   in (secondT - firstT) * (thirdConstant - firstConstant)+        - (secondConstant - firstConstant) * (thirdT - firstT)++upperChainPointSet+  :: (ExactRational, ExactRational)+  -> [LiftedPoint label]+  -> [LiftedPoint label]+  -> Set (LiftedPoint label)+upperChainPointSet direction chain =+  snd+    . List.foldl'+        classifyPoint+        (consecutivePairs chain, Set.empty)+ where+  classifyPoint+    :: ([(LiftedPoint label, LiftedPoint label)], Set (LiftedPoint label))+    -> LiftedPoint label+    -> ([(LiftedPoint label, LiftedPoint label)], Set (LiftedPoint label))+  classifyPoint (remainingEdges, members) point =+    let pointParameter = slopeParameter direction point+        candidateEdges =+          List.dropWhile+            ((< pointParameter) . slopeParameter direction . snd)+            remainingEdges+        retainedMembers =+          case candidateEdges of+            (firstPoint, secondPoint) : _+              | slopeParameter direction firstPoint <= pointParameter+                  && slopeChainTurn direction firstPoint secondPoint point == 0 ->+                  Set.insert point members+            _ -> members+     in (candidateEdges, retainedMembers)++fullLineEdge+  :: (LiftedPoint label, LiftedPoint label)+  -> Either (RegularTopologyError label) (GeneratorRegularEdge label)+fullLineEdge (firstPoint, secondPoint) = do+  axis <- generatorAxis firstPoint secondPoint+  pure (regularEdge firstPoint secondPoint (GeneratorDualLine axis))++topologyFromHullFaces+  :: [(HullFace label, ExactPoint)]+  -> Either+      (RegularTopologyError label)+      ([GeneratorRegularFace label], [GeneratorRegularEdge label])+topologyFromHullFaces dualFaces = do+  let dualByFace = IntMap.fromList [(hullFaceId face, dual) | (face, dual) <- dualFaces]+      edgeIncidence =+        List.foldl'+          (\incidence (face, _) ->+             List.foldl'+               (\current (firstPoint, secondPoint, thirdPoint) ->+                  Map.insertWith+                    (<>)+                    (hullEdgeKey firstPoint secondPoint)+                    ((face, firstPoint, secondPoint, thirdPoint) :| [])+                    current)+               incidence+               (faceEdgesWithOpposite face))+          Map.empty+          dualFaces+  edges <- traverse (dualEdge dualByFace) (Map.elems edgeIncidence)+  pure+    ( [ GeneratorRegularFace+          (pointLabel (hullFaceFirst face))+          (pointLabel (hullFaceSecond face))+          (pointLabel (hullFaceThird face))+          dual+      | (face, dual) <- dualFaces+      ]+    , edges+    )++hullFaceDual+  :: HullFace label+  -> Either (RegularTopologyError label) (HullFace label, ExactPoint)+hullFaceDual face = do+  let firstPoint = hullFaceFirst face+      secondPoint = hullFaceSecond face+      thirdPoint = hullFaceThird face+  dual <-+    exactGeneratorFaceDual+      (liftedPointGenerator firstPoint)+      (liftedPointGenerator secondPoint)+      (liftedPointGenerator thirdPoint)+  pure (face, dual)++faceEdgesWithOpposite+  :: HullFace label+  -> [(LiftedPoint label, LiftedPoint label, LiftedPoint label)]+faceEdgesWithOpposite face =+  [ (hullFaceFirst face, hullFaceSecond face, hullFaceThird face)+  , (hullFaceSecond face, hullFaceThird face, hullFaceFirst face)+  , (hullFaceThird face, hullFaceFirst face, hullFaceSecond face)+  ]++dualEdge+  :: IntMap.IntMap ExactPoint+  -> NonEmpty (HullFace label, LiftedPoint label, LiftedPoint label, LiftedPoint label)+  -> Either (RegularTopologyError label) (GeneratorRegularEdge label)+dualEdge dualByFace incidences =+  case incidences of+    (face, firstPoint, secondPoint, thirdPoint) :| [] -> do+      dual <-+        maybe+          (Left (RegularTopologyFlatBoundary (pointLabel firstPoint) (pointLabel secondPoint) (pointLabel thirdPoint)))+          Right+          (IntMap.lookup (hullFaceId face) dualByFace)+      ray <- boundaryDualRay dual firstPoint secondPoint thirdPoint+      pure (regularEdge firstPoint secondPoint (GeneratorDualRay ray))+    (firstFace, firstPoint, secondPoint, _) :| [(secondFace, _, _, _)] -> do+        firstDual <- requiredDual firstPoint secondPoint firstFace+        secondDual <- requiredDual firstPoint secondPoint secondFace+        dual <-+          exactGeneratorDualBetween+            (pointLabel firstPoint)+            (pointLabel secondPoint)+            firstDual+            secondDual+        pure (regularEdge firstPoint secondPoint dual)+    (_, firstPoint, secondPoint, _) :| _ ->+      Left+        ( RegularTopologyNonManifoldEdge+            (pointLabel firstPoint)+            (pointLabel secondPoint)+            (NonEmpty.length incidences)+        )+ where+  requiredDual+    :: LiftedPoint label+    -> LiftedPoint label+    -> HullFace label+    -> Either (RegularTopologyError label) ExactPoint+  requiredDual firstPoint secondPoint face =+    maybe+      (Left (RegularTopologyFlatBoundary (pointLabel firstPoint) (pointLabel secondPoint) (pointLabel (hullFaceThird face))))+      Right+      (IntMap.lookup (hullFaceId face) dualByFace)++boundaryDualRay+  :: ExactPoint+  -> LiftedPoint label+  -> LiftedPoint label+  -> LiftedPoint label+  -> Either (RegularTopologyError label) ExactRay+boundaryDualRay originPoint firstPoint secondPoint thirdPoint = do+  exactGeneratorBoundaryDualRay+    originPoint+    (liftedPointGenerator firstPoint)+    (liftedPointGenerator secondPoint)+    (liftedPointGenerator thirdPoint)++generatorAxis+  :: LiftedPoint label+  -> LiftedPoint label+  -> Either (RegularTopologyError label) ExactAffineLine+generatorAxis firstPoint secondPoint =+  first+    (RegularTopologyAxisInvalid (pointLabel firstPoint) (pointLabel secondPoint))+    (exactGeneratorAxis (liftedPointGenerator firstPoint) (liftedPointGenerator secondPoint))++regularEdge+  :: LiftedPoint label+  -> LiftedPoint label+  -> GeneratorDualGeometry+  -> GeneratorRegularEdge label+regularEdge firstPoint secondPoint dual =+  if liftedPointId firstPoint <= liftedPointId secondPoint+    then GeneratorRegularEdge (liftedPointGenerator firstPoint) (liftedPointGenerator secondPoint) dual+    else GeneratorRegularEdge (liftedPointGenerator secondPoint) (liftedPointGenerator firstPoint) dual++pointLabel :: LiftedPoint label -> label+pointLabel = exactPowerGeneratorLabel . liftedPointGenerator++topologyFromParts+  :: NonEmpty (LiftedPoint label)+  -> Set (LiftedPoint label)+  -> Set (LiftedPoint label)+  -> [GeneratorRegularFace label]+  -> [GeneratorRegularEdge label]+  -> GeneratorRegularTopology label+topologyFromParts points visible lowerDimensional faces edges =+  let dispositions =+        fmap+          (\point ->+             ( liftedPointGenerator point+             , if Set.member point visible+                 then RegularGeneratorVisible+                 else+                   if Set.member point lowerDimensional+                     then RegularGeneratorLowerDimensional+                     else RegularGeneratorHidden+             ))+          points+      visibleCount = Set.size visible+      lowerCount = Set.size lowerDimensional+      inputCount = NonEmpty.length points+   in GeneratorRegularTopology+        { generatorRegularDispositions = dispositions+        , generatorRegularFaces = faces+        , generatorRegularEdges = edges+        , generatorRegularReceipt =+            GeneratorRegularReceipt+              { generatorRegularInputSites = inputCount+              , generatorRegularVisibleSites = visibleCount+              , generatorRegularLowerDimensionalSites = lowerCount+              , generatorRegularHiddenSites = inputCount - visibleCount - lowerCount+              , generatorRegularFaceCount = length faces+              , generatorRegularEdgeCount = length edges+              }+        }