moonlight-triangulation 1.3.0.0 → 1.3.0.2
raw patch · 31 files changed
+2689/−365 lines, 31 filesdep ~basedep ~containersdep ~moonlight-triangulationPVP: minor bump suggested
API additions: PVP suggests at least a minor version bump
Dependency ranges changed: base, containers, moonlight-triangulation
API changes (from Hackage documentation)
+ Moonlight.Triangulation: AlphaCircumradiusDivisionFailed :: FaceId -> ExactArithmeticError -> AlphaFiltrationError
+ Moonlight.Triangulation: AlphaEdgeBirthWitnessMissing :: UndirectedEdgeId -> AlphaFiltrationError
+ Moonlight.Triangulation: AlphaFaceDegenerate :: FaceId -> AlphaFiltrationError
+ Moonlight.Triangulation: AlphaFaceEdgeIncidenceInvalid :: FaceId -> UndirectedEdgeId -> [VertexId] -> AlphaFiltrationError
+ Moonlight.Triangulation: AlphaFaceVerticesUnavailable :: FaceId -> AlphaFiltrationError
+ Moonlight.Triangulation: AlphaFiltrationCellSelection :: CellSelectionError -> AlphaFiltrationError
+ Moonlight.Triangulation: AlphaIncidentFaceBirthMissing :: UndirectedEdgeId -> FaceId -> AlphaFiltrationError
+ Moonlight.Triangulation: alphaBirthDenominator :: AlphaBirth -> Integer
+ Moonlight.Triangulation: alphaBirthFromRadiusSquared :: RadiusSquared -> AlphaBirth
+ Moonlight.Triangulation: alphaBirthNumerator :: AlphaBirth -> Integer
+ Moonlight.Triangulation: alphaBirthToDouble :: AlphaBirth -> Double
+ Moonlight.Triangulation: alphaComplexAtBirth :: AlphaBirth -> AlphaFiltration -> Either CellSelectionError ExactCellSet
+ Moonlight.Triangulation: alphaComplexAtRadius :: RadiusSquared -> AlphaFiltration -> Either CellSelectionError ExactCellSet
+ Moonlight.Triangulation: alphaEdgeBirth :: AlphaFiltration -> UndirectedEdgeId -> Maybe AlphaBirth
+ Moonlight.Triangulation: alphaFaceBirth :: AlphaFiltration -> FaceId -> Maybe AlphaBirth
+ Moonlight.Triangulation: alphaFiltration :: DelaunayTriangulation vertex -> Either AlphaFiltrationError AlphaFiltration
+ Moonlight.Triangulation: alphaFiltrationCellSet :: AlphaFiltration -> ExactCellSet
+ Moonlight.Triangulation: alphaFiltrationCriticalBirths :: AlphaFiltration -> [AlphaBirth]
+ Moonlight.Triangulation: alphaVertexBirth :: AlphaFiltration -> VertexId -> Maybe AlphaBirth
+ Moonlight.Triangulation: data AlphaBirth
+ Moonlight.Triangulation: data AlphaFiltration
+ Moonlight.Triangulation: data AlphaFiltrationError
+ Moonlight.Triangulation: radiusSquaredValue :: RadiusSquared -> Double
Files
- CHANGELOG.md +28/−0
- LICENSE +1/−1
- README.md +69/−4
- bench/alpha/Main.hs +128/−0
- bench/delaunay-compare/Moonlight/Triangulation/Bench/DelaunayCompare/Picture.hs +0/−60
- bench/delaunay-compare/README.md +41/−32
- bench/delaunay-compare/results/moonlight-delaunay-compare-small.svg +415/−0
- bench/support/BenchSupport.hs +18/−7
- bindings/README.md +4/−0
- bindings/python/src/moonlight_triangulation/__init__.py +13/−5
- bindings/python/tests/test_binding.py +26/−0
- bindings/typescript/src/index.ts +71/−45
- bindings/typescript/test/binding.test.ts +25/−0
- moonlight-triangulation.cabal +46/−9
- src-build/Moonlight/Triangulation/BulkLoad.hs +83/−61
- src-build/Moonlight/Triangulation/Internal/CircleSweep.hs +185/−23
- src-cell-complex/Moonlight/Triangulation/CellComplex.hs +405/−10
- src-core/Moonlight/Triangulation/Internal/Dyadic.hs +106/−12
- src-core/Moonlight/Triangulation/Internal/ExactRational.hs +47/−3
- src-dcel/Moonlight/Triangulation/Alpha.hs +424/−0
- src-dcel/Moonlight/Triangulation/FloodFillIterator.hs +4/−46
- src-dcel/Moonlight/Triangulation/Internal/CellSet.hs +28/−0
- src-dcel/Moonlight/Triangulation/Internal/Mutable.hs +69/−15
- src-dcel/Moonlight/Triangulation/Internal/Types.hs +31/−0
- src-dcel/Moonlight/Triangulation/Types.hs +4/−0
- src-dcel/Moonlight/Triangulation/Validation.hs +11/−4
- src-public/Moonlight/Triangulation.hs +38/−4
- src-serialize/Moonlight/Triangulation/Serialization.hs +3/−19
- test/cell-complex/Main.hs +321/−4
- test/native/Moonlight/Triangulation/NativeSpec.hs +1/−0
- test/serialization/Moonlight/Triangulation/SerializationSpec.hs +44/−1
CHANGELOG.md view
@@ -6,6 +6,34 @@ The serialization format carries its own version tag, independent of the package version; any change to it is recorded here explicitly. +## 1.3.0.2 - 2026-08-22++* Add an opaque exact alpha filtration over every resident Delaunay vertex,+ edge, and bounded face, exact threshold cell selection, checked integral+ cellular-chain lowering, and direct lowering into Moonlight Homology's+ ordered filtered complex.+* Add exact Gabriel/non-Gabriel birth laws, boundary monotonicity, degenerate+ support, and a cocircular-square persistence witness whose degree-one class+ is born at squared radius 1 and dies at squared radius 2.+* Add a focused benchmark comparing one persistence reduction with recomputing+ homology at every exact critical radius over the same prepared geometry.+* Compute exact edge lengths and face circumradii as aligned dyadic integer+ kernels with one rational normalization at publication, and lower already+ ordered DCEL incidence before Homology checks and seals the chain law.+* Batch every critical-radius Betti profile through Homology's dense rank+ quotient. On the shared 1,000-site Spade diagnostic this cuts Moonlight from+ 941.312 ms and 3.11 GB allocated to a 34.177 ms median and 96.0 MB while+ preserving the exact receipt.+* Remove the duplicate packed radial-geometry arena from bulk Delaunay ingress.+ Coordinates remain in dense planes while a 32-bit vertex-id plane descends+ through prefix radix order and exact equal-prefix normalization.+* Refresh the complete 120-case construction board. Moonlight's mean is below+ plain Spade and Spade hierarchy on all 24 fixtures, with non-overlapping+ reported 2σ on 18 and 22 fixtures respectively.+* Remove promotional advantage callouts from both SVG projections, ship the+ regenerated small SVG with the package, and reference it directly from the+ README instead of displaying a stale remote mirror copy.+ ## 1.3.0.0 * Add the public GHC-9.14 `cell-complex` component. It interprets an admitted
LICENSE view
@@ -1,6 +1,6 @@ MIT License -Copyright (c) 2026 Fable, Blue Rose+Copyright (c) 2026 agent, Blue Rose Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal
README.md view
@@ -37,6 +37,8 @@ | `minkowskiSum` | polygonal convolution | Exact convex or general polygonal Minkowski addition, with work recorded in a receipt. | | `erodeBy` / `openWith` / `closeWith` | regularized morphology | Full-dimensional polygonal erosion and its opening/closing compositions for an origin-anchored convex kernel. | | `fromExactCellSet` | topology interpretation | An opaque `DCELComplex` view of the already-admitted cell inventory; no mesh copy or second validation authority. |+| `alphaFiltration` | exact scale geometry | One handle-indexed birth section for every resident vertex, edge, and bounded face of a Delaunay triangulation. |+| `filteredAlphaComplex` | persistent topology | The same resident cells and exact births lowered into Homology's checked filtered chain complex. | ### Choose the operation @@ -68,7 +70,8 @@ | Measure a closed exact selection or region | `cellValuations` or `regionValuations` | Do not coerce irrational perimeter to an allegedly exact `Double`; inspect its radical expression and certified interval. | | Expand a polygon by another polygon | `minkowskiSum` or `polygonOffset` | Do not add resident mesh vertices pairwise and rebuild; morphology acts on the represented continuum. | | Erode, open, or close a polygonal region | `erodeBy`, `openWith`, or `closeWith` | Do not infer a hidden universe or invent lower-dimensional polygons; the result is regularized to representable 2D area. |-| Keep Delaunay faces below an alpha threshold | `alphaShapeContainsFace threshold mesh` | This is the filled-face filtration, not palette logic and not a claim to expose every simplex of a full alpha-complex carrier. |+| Render only Delaunay faces below one alpha threshold | `alphaShapeContainsFace threshold mesh` | This constant-work-per-face projection deliberately does not build the full filtration. |+| Ask topology questions across many radii | `alphaFiltration mesh`, then `filteredAlphaComplex` and `mod2PersistentPairs` | Do not rebuild disc unions, triangulations, or homology independently at every radius. | | Require construction-independent numbering | `canonicalize` at the observation boundary | Do not canonicalize every intermediate value; it is intentionally global work. | ### Cell-complex and category interpretation@@ -105,6 +108,67 @@ incidence category and normalized nerve, validates their overlap counts, and emits only the checked visualization projection. +### Exact alpha filtration and persistent topology++The full alpha path is a different observation of the same resident Delaunay+DCEL. `alphaFiltration` computes exact squared-radius births once: vertices at+zero, Gabriel edges at one quarter of their squared length, other edges at their+least incident-face witness, and faces at their circumradius squared.+`filteredAlphaComplex` lowers that section to the one Homology persistence+owner. `AlphaBirth` is an opaque reduced rational; equality and ordering never+round through `Double`.++This is the planar alpha-complex shortcut to the Cech filtration: the Delaunay+subcomplex has the same homotopy type as the corresponding union of discs, but+has only the resident planar cells. One reduction therefore yields the barcode+for every critical radius. `persistentBettiAtCriticalValues` sweeps all of+those sublevels through dense complex-relative ranks without rebuilding+geometry, reducing another boundary matrix, or rescanning the barcode per+radius.++The minimal component path names only Delaunay construction, the DCEL alpha+owner, the opt-in lowering component, and Homology. A consumer already using+the broad `moonlight-triangulation` facade may substitute it for `build` and+`dcel` below.++```cabal+build-depends:+ moonlight-triangulation:build >= 1.3 && < 1.4,+ moonlight-triangulation:dcel >= 1.3 && < 1.4,+ moonlight-triangulation:cell-complex >= 1.3 && < 1.4,+ moonlight-homology >= 0.1 && < 0.2+```++```haskell+import Data.Foldable (traverse_)+import qualified Data.Vector as Vector+import Moonlight.Homology.Persistence+ ( filteredCriticalValues+ , mod2PersistentPairs+ , persistentBettiAtCriticalValues+ )+import Moonlight.Triangulation.Alpha (alphaFiltration)+import Moonlight.Triangulation.BulkLoad (delaunayGeometry)+import Moonlight.Triangulation.CellComplex (filteredAlphaComplex)+import Moonlight.Triangulation.Types (Point (..))++main :: IO ()+main = do+ mesh <- requireRight (delaunayGeometry (Vector.fromList [Point (-1) (-1), Point 1 (-1), Point 1 1, Point (-1) 1]))+ filtration <- requireRight (alphaFiltration mesh)+ filtered <- requireRight (filteredAlphaComplex filtration)+ pairs <- requireRight (mod2PersistentPairs filtered)+ profiles <- requireRight (persistentBettiAtCriticalValues filtered pairs)+ traverse_ print (zip (filteredCriticalValues filtered) profiles)++requireRight :: Show failure => Either failure value -> IO value+requireRight = either (fail . show) pure+```++The focused `moonlight-triangulation-alpha-bench` measures exact preparation,+one persistence reduction plus all critical-radius queries, and repeated+threshold homology as separate lanes over the same prebuilt Delaunay mesh.+ ## Foreign bindings The `moonlight-triangulation-c` foreign-library component exposes opaque@@ -185,7 +249,7 @@ all 24 fixture summaries across all five implementations before timing, and derives both pictures from the retained 120-case CSV. -+ The [benchmark README](bench/delaunay-compare/README.md) owns the single retained board, its result, receipt, big-point picture, and reproduction@@ -629,8 +693,9 @@ * `.Cdt` — constraint recovery by conflict strip, with requeue on re-intersection. * `.PointLocation`, `.HintGenerator` — walk location and Delaunay hierarchy hints. * `.IntersectionIterator`, `.FloodFillIterator` — ordered line traversal,- barrier fill, labelled face components, boundary loops, and alpha-face- filtration.+ barrier fill, labelled face components, and boundary loops.+* `.Alpha` — exact vertex/edge/face births, closed alpha subcomplexes, and the+ constant-work bounded-face membership projection. * `.Exact`, `.CellSet`, `.Region`, `.Valuation` — exact rational geometry, closed cell selections, admitted polygonal regions, and intrinsic measures. * `.Overlay`, `.Minkowski` — labelled common refinement, planar Booleans,
+ bench/alpha/Main.hs view
@@ -0,0 +1,128 @@+-- | Measure the exact alpha shortcut against recomputing homology at every+-- critical radius. Delaunay construction is deliberately outside every timed+-- lane: both algorithms consume the same already-built geometry.+module Main (main) where++import BenchSupport+ ( randomPoints+ , requireRight+ , timedProjection+ , timedValue+ )+import Data.Map.Strict qualified as Map+import Data.Vector qualified as Vector+import Moonlight.Homology.Boundary+ ( degreeCardinality+ , maxHomologicalDegree+ , restrictComplex+ )+import Moonlight.Homology.Chain+ ( HomologicalDegree (..)+ , HomologyFailure+ )+import Moonlight.Homology.Persistence+ ( FilteredFiniteChainComplex+ , criticalBettiVectors+ , filteredBaseComplex+ , filteredCellBirths+ , filteredCriticalValues+ , mod2PersistentPairsWithCriticalBettiTable+ )+import Moonlight.Homology.Topology (freeBettiVector)+import Moonlight.Triangulation.Alpha+ ( AlphaBirth+ , alphaFiltration+ )+import Moonlight.Triangulation.BulkLoad (delaunayGeometry)+import Moonlight.Triangulation.CellComplex (filteredAlphaComplex)+import Moonlight.Triangulation.Types (DelaunayTriangulation)++type PreparedAlpha = FilteredFiniteChainComplex AlphaBirth Int++main :: IO ()+main = benchmarkAlphaFiltration 64++benchmarkAlphaFiltration :: Int -> IO ()+benchmarkAlphaFiltration siteCount = do+ triangulation <-+ requireRight+ (delaunayGeometry (Vector.fromList (randomPoints 0x9e3779b97f4a7c15 siteCount)))+ filtered <-+ timedProjection+ (benchmarkLabel siteCount "alpha-filtered-complex")+ preparationReceipt+ (prepareAlpha triangulation)+ let thresholds = filteredCriticalValues filtered+ maximumDegree = maxHomologicalDegreeValue filtered+ persistentProfile <-+ timedValue+ (benchmarkLabel siteCount "persistence-critical-profile")+ ( do+ (_, bettiTable) <-+ requireRight (mod2PersistentPairsWithCriticalBettiTable filtered)+ pure (criticalBettiVectors bettiTable)+ )+ repeatedProfile <-+ timedValue+ (benchmarkLabel siteCount "repeated-threshold-homology")+ (requireRight (traverse (bettiVectorAtThreshold maximumDegree filtered) thresholds))+ if persistentProfile == repeatedProfile+ then+ putStrLn+ ( benchmarkLabel siteCount "agreement"+ <> ": critical-radii="+ <> show (length thresholds)+ <> " profiles="+ <> show (length persistentProfile)+ <> " prepared-cells="+ <> show (preparationReceipt filtered)+ )+ else fail (benchmarkLabel siteCount "profile-mismatch")++benchmarkLabel :: Int -> String -> String+benchmarkLabel siteCount suffix =+ "alpha/n" <> show siteCount <> "/" <> suffix++prepareAlpha+ :: DelaunayTriangulation ()+ -> IO PreparedAlpha+prepareAlpha triangulation = do+ filtration <- requireRight (alphaFiltration triangulation)+ requireRight (filteredAlphaComplex filtration)++preparationReceipt :: PreparedAlpha -> (Int, Int, [Int])+preparationReceipt filtered =+ let finite = filteredBaseComplex filtered+ HomologicalDegree maximumDegree = maxHomologicalDegree finite+ in ( length (filteredCriticalValues filtered)+ , Map.size (filteredCellBirths filtered)+ , fmap+ (degreeCardinality finite . HomologicalDegree)+ [0 .. maximumDegree]+ )++maxHomologicalDegreeValue+ :: FilteredFiniteChainComplex filtration coefficient+ -> Int+maxHomologicalDegreeValue filtered =+ case maxHomologicalDegree (filteredBaseComplex filtered) of+ HomologicalDegree degreeValue -> degreeValue++bettiVectorAtThreshold+ :: Int+ -> FilteredFiniteChainComplex AlphaBirth Int+ -> AlphaBirth+ -> Either HomologyFailure [Int]+bettiVectorAtThreshold maximumDegree filtered threshold = do+ restricted <-+ restrictComplex+ (Map.keysSet (Map.filter (<= threshold) (filteredCellBirths filtered)))+ (filteredBaseComplex filtered)+ pure (padBettiVector maximumDegree (freeBettiVector restricted))++padBettiVector :: Int -> [Int] -> [Int]+padBettiVector maximumDegree bettiVector =+ bettiVector+ <> replicate+ (max 0 (maximumDegree + 1 - length bettiVector))+ 0
bench/delaunay-compare/Moonlight/Triangulation/Bench/DelaunayCompare/Picture.hs view
@@ -293,10 +293,6 @@ , "<text x=\"" <> coordinate (panelLeft - 55) <> "\" y=\"" <> coordinate (plotTop + plotHeight / 2) <> "\" text-anchor=\"middle\" font-family=\"Helvetica, Arial, sans-serif\" font-size=\"12\" fill=\"rgb(20,20,20)\" transform=\"rotate(-90 " <> coordinate (panelLeft - 55) <> " " <> coordinate (plotTop + plotHeight / 2) <> ")\">CPU time (ms)</text>" ] <> concatMap (renderSeries panelLeft panelWidth pointCounts axisMaximum) (picturePanelSeries panel)- <> maybe- []- (renderAdvantageCallout panelLeft panelWidth pointCounts axisMaximum)- (bestReportedSpreadAdvantage (picturePanelSeries panel)) <> ["</g>"] where panelLabel = pointDistributionLabel (picturePanelDistribution panel)@@ -351,62 +347,6 @@ y = observationY observation upperY = pictureY axisMaximum (observationUpperMilliseconds observation) lowerY = pictureY axisMaximum (observationLowerMilliseconds observation)--data ReportedSpreadAdvantage = ReportedSpreadAdvantage- { advantagePointCount :: !Int- , advantageMoonlight :: !BenchmarkObservation- , advantageCompetitor :: !Implementation- , advantageReduction :: !Double- }--bestReportedSpreadAdvantage :: [PictureSeries] -> Maybe ReportedSpreadAdvantage-bestReportedSpreadAdvantage series = foldr selectLargerReduction Nothing reportedSpreadAdvantages- where- moonlightObservations =- concatMap pictureSeriesObservations (filter ((== Moonlight) . pictureSeriesImplementation) series)- spadeSeries =- filter- ((`elem` [Spade, SpadeHierarchy]) . pictureSeriesImplementation)- series- reportedSpreadAdvantages =- [ ReportedSpreadAdvantage- { advantagePointCount = pointCountValue (fixturePointCount moonlightFixture)- , advantageMoonlight = moonlight- , advantageCompetitor = pictureSeriesImplementation competitorSeries- , advantageReduction = 1 - benchmarkMeanPicoseconds moonlight / benchmarkMeanPicoseconds competitor- }- | moonlight <- moonlightObservations- , let moonlightFixture = observationFixture (benchmarkKey moonlight)- , competitorSeries <- spadeSeries- , competitor <- pictureSeriesObservations competitorSeries- , fixturePointCount (observationFixture (benchmarkKey competitor)) == fixturePointCount moonlightFixture- , benchmarkMeanPicoseconds moonlight + benchmarkSpreadPicoseconds moonlight- < benchmarkMeanPicoseconds competitor - benchmarkSpreadPicoseconds competitor- ]- selectLargerReduction candidate = \case- Nothing -> Just candidate- Just incumbent ->- Just- ( if advantageReduction candidate > advantageReduction incumbent- then candidate- else incumbent- )--renderAdvantageCallout :: Double -> Double -> [Int] -> Double -> ReportedSpreadAdvantage -> [String]-renderAdvantageCallout panelLeft panelWidth pointCounts axisMaximum advantage =- [ "<g aria-label=\"Moonlight comparison highlight\">"- , " <line x1=\"" <> coordinate pointX <> "\" y1=\"" <> coordinate pointY <> "\" x2=\"" <> coordinate (labelX + 8) <> "\" y2=\"" <> coordinate (labelY + 18) <> "\" stroke=\"rgb(255,140,0)\" stroke-width=\"1.2\"/>"- , " <rect x=\"" <> coordinate labelX <> "\" y=\"" <> coordinate labelY <> "\" width=\"250\" height=\"38\" rx=\"6\" fill=\"white\" stroke=\"rgb(255,140,0)\" stroke-width=\"1.2\"/>"- , " " <> svgText (labelX + 10) (labelY + 16) "start" 11 "500" (show reductionPercent <> "% less time than " <> implementationLabel (advantageCompetitor advantage))- , " " <> svgText (labelX + 10) (labelY + 31) "start" 10 "400" (pointCountLabel (advantagePointCount advantage) <> " points · reported spreads do not overlap")- , "</g>"- ]- where- pointX = pictureX panelLeft panelWidth pointCounts (advantagePointCount advantage)- pointY = pictureY axisMaximum (observationMeanMilliseconds (advantageMoonlight advantage))- labelX = min (panelLeft + panelWidth - 258) (pointX + 28)- labelY = max (picturePlotTop pictureLayout + 12) (pointY - 58)- reductionPercent = round (100 * advantageReduction advantage) :: Int implementationColor :: Implementation -> String implementationColor = \case
bench/delaunay-compare/README.md view
@@ -39,44 +39,47 @@ - timed cases: `120` - fixture summaries agreeing across all five implementations: `24 / 24`-- Moonlight mean below plain Spade: `12 / 24`-- Moonlight mean below Spade hierarchy: `17 / 24`-- Moonlight below plain Spade with non-overlapping reported 2σ: `0 / 24`-- Moonlight below Spade hierarchy with non-overlapping reported 2σ: `1 / 24`-- lowest fixture mean: `delaunator` `23 / 24`; `cdt` `1 / 24`-- median reported 2σ / mean: Moonlight `7.86%`; plain Spade `6.74%`+- Moonlight mean below plain Spade: `24 / 24`+- Moonlight mean below Spade hierarchy: `24 / 24`+- Moonlight below plain Spade with non-overlapping reported 2σ: `18 / 24`+- Moonlight below Spade hierarchy with non-overlapping reported 2σ: `22 / 24`+- lowest fixture mean: `delaunator` `21 / 24`; `cdt` `3 / 24`+- median reported 2σ / mean: Moonlight `6.88%`; plain Spade `6.67%` Selected means are in milliseconds. Δ is `100 × (Moonlight / competitor − 1)`. | fixture | points | Moonlight | plain Spade | Δ | Spade hierarchy | Δ | |---|---:|---:|---:|---:|---:|---:|-| local insertion | 4,000 | 0.663 | 0.748 | −11.3% | 0.784 | −15.4% |-| local insertion | 6,000 | 1.121 | 1.245 | −10.0% | 1.320 | −15.1% |-| local insertion | 8,000 | 1.552 | 1.736 | −10.6% | 1.824 | −14.9% |-| local insertion | 250,000 | 68.017 | 71.675 | −5.1% | 76.208 | −10.7% |-| uniform | 4,000 | 0.676 | 0.693 | −2.4% | 0.729 | −7.3% |-| uniform | 14,000 | 3.591 | 2.959 | +21.4% | 3.171 | +13.3% |-| uniform | 250,000 | 73.329 | 72.076 | +1.7% | 76.447 | −4.1% |+| local insertion | 4,000 | 0.587 | 0.771 | −23.8% | 0.800 | −26.6% |+| local insertion | 6,000 | 0.971 | 1.300 | −25.3% | 1.365 | −28.9% |+| local insertion | 8,000 | 1.354 | 1.776 | −23.7% | 1.903 | −28.8% |+| local insertion | 250,000 | 60.300 | 72.713 | −17.1% | 77.496 | −22.2% |+| uniform | 4,000 | 0.552 | 0.706 | −21.8% | 0.741 | −25.6% |+| uniform | 14,000 | 2.663 | 3.067 | −13.2% | 3.318 | −19.7% |+| uniform | 250,000 | 60.028 | 72.783 | −17.5% | 77.435 | −22.5% | The measured Moonlight path is the canonical circle-sweep bulk loader. Its-local descent carries an already-read hull angle through candidate search,-rejects incompatible closure sections before paying for exact orientation,-and compiles the exact predicate and repair kernels at their LLVM `-O3`-boundary.+geometry ingress keeps coordinates in immutable planes, moves only a 32-bit+vertex-id plane through three radix passes, and restores the complete+`(distance, x, y, vertex)` order inside equal-prefix fibres before descent.+Fresh dense vertices are committed once, and the sweep consumes reserved cells+without constructing a second packed geometry arena. Receipt identity: -- command completed: `2026-08-17T10:39:30Z`, 297.39 s over 120 cases+- command completed: `2026-08-23T00:00:15Z`, 318.68 s over 120 cases - host: Apple M4 Pro, arm64, macOS 26.5.2, GHC 9.14.1, rustc 1.92.0-- Haskell build: Cabal `-O1` package profile, comparison executable `-O2`, hot- predicate, sweep, and repair modules LLVM `-O3`-- checkout HEAD: `6da6feb1a6`, with seven modified files in the tree, none of- them on the measured path-- CSV: [`moonlight-delaunay-compare-2026-08-17.csv`](results/moonlight-delaunay-compare-2026-08-17.csv), SHA-256- `eba77eb960790f2e1cbe3f82c112e8690721de4e5aa3c69961925670e0822162`+- Haskell build: cabal-install 3.16.1.0, focused `-O1` package profile,+ comparison executable `-O2`, hot predicate, sweep, and repair modules LLVM+ `-O3`+- checkout base HEAD: `41906f44f0`, plus the three optimized runtime-source+ owners named by the manifest+- CSV: [`moonlight-delaunay-compare-2026-08-22.csv`](results/moonlight-delaunay-compare-2026-08-22.csv), SHA-256+ `cce6af4ac7d570970a69124f5648db4e13ca41e91e05383306407a0dca323be7` - benchmark runtime-source hashes:- [`moonlight-delaunay-compare-2026-08-17.source-sha256`](results/moonlight-delaunay-compare-2026-08-17.source-sha256)+ [`moonlight-delaunay-compare-2026-08-22.source-sha256`](results/moonlight-delaunay-compare-2026-08-22.source-sha256), SHA-256+ `069d3bb1639b6d865fe4491b2fbbb9d0cdc8324173e30bebf690ed6592190c20` The Haskell projection in `Moonlight.Triangulation.Bench.DelaunayCompare.Picture` parses the tasty-bench@@ -91,6 +94,7 @@ ```console scripts/safe-cabal.sh run \ moonlight-triangulation:exe:moonlight-triangulation-delaunay-compare \+ --project-file=cabal.project.triangulation-dev \ -- --list-tests ``` @@ -99,6 +103,7 @@ ```console scripts/safe-cabal.sh run \ moonlight-triangulation:exe:moonlight-triangulation-delaunay-compare \+ --project-file=cabal.project.triangulation-dev \ -- --time-mode wall -j1 ``` @@ -116,7 +121,8 @@ ```console caffeinate -i scripts/safe-cabal.sh run \- moonlight-triangulation:exe:moonlight-triangulation-delaunay-compare -- \+ moonlight-triangulation:exe:moonlight-triangulation-delaunay-compare \+ --project-file=cabal.project.triangulation-dev -- \ --stdev 5 --timeout 30s --time-mode cpu -j1 \ --csv foundation/moonlight-triangulation/bench/delaunay-compare/results/moonlight-delaunay-compare-YYYY-MM-DD.csv \ --color never --hide-progress --min-duration-to-report 1h@@ -127,16 +133,19 @@ ```console scripts/safe-cabal.sh run \- moonlight-triangulation:exe:moonlight-triangulation-delaunay-pictures -- \+ moonlight-triangulation:exe:moonlight-triangulation-delaunay-pictures \+ --project-file=cabal.project.triangulation-dev -- \ foundation/moonlight-triangulation/bench/delaunay-compare/results/moonlight-delaunay-compare-YYYY-MM-DD.csv \ foundation/moonlight-triangulation/bench/delaunay-compare/results ``` -The runtime-source manifest beside the CSV is the list of files that entered-timed actions, hashed at run time; picture sources and package-only metadata-are excluded because they do not. The receipt and SVGs stay package-owned-beside the benchmark but are not Cabal package inputs, because taking a-measurement must not rebuild the triangulation library.+The runtime-source manifest beside the CSV is regenerated from the measured+worktree immediately after a completed command. It hashes the closed list of+files that entered timed actions; picture sources and package-only metadata are+excluded because they did not. The receipt, CSV, and big SVG stay+benchmark-owned beside the executable and outside Cabal package inputs. The+small SVG is also the package README projection, so it enters the sdist as+documentation without entering any compiled component. ## Upstream-compatible fixtures
+ bench/delaunay-compare/results/moonlight-delaunay-compare-small.svg view
@@ -0,0 +1,415 @@+<?xml version="1.0" encoding="utf-8"?>+<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 1280.0 720.0" role="img" aria-labelledby="title description">+<title id="title">Moonlight Delaunay construction · small point sets</title>+<desc id="description">CPU construction time for Spade, cdt, delaunator, and Moonlight over the upstream Delaunay comparison fixtures. Lower lines are faster.</desc>+<defs>+ <pattern id="triangulation-mesh" width="54" height="46" patternUnits="userSpaceOnUse">+ <path d="M0 46L27 0L54 46ZM0 46L54 46M27 0L27 46" fill="none" stroke="rgb(255,140,0)" stroke-width="0.7"/>+ </pattern>+ <clipPath id="header-mesh-clip"><rect x="955" y="8" width="290" height="82" rx="20"/></clipPath>+</defs>+<rect width="1280.0" height="720.0" fill="white"/>+<rect x="955" y="8" width="290" height="82" fill="url(#triangulation-mesh)" opacity="0.10" clip-path="url(#header-mesh-clip)"/>+<text x="640.0" y="34.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="24" font-weight="500" fill="rgb(20,20,20)">Moonlight Delaunay construction · small point sets</text>+<text x="640.0" y="61.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="13" font-weight="400" fill="rgb(20,20,20)">Moonlight against Spade's construction referents · CPU time · lower is faster</text>+<g aria-label="spade 2">+ <line x1="105.0" y1="101" x2="139.0" y2="101" stroke="rgb(0,0,139)" stroke-width="2.0"/>+ <circle cx="122.0" cy="101.0" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <text x="148.0" y="105.0" text-anchor="start" font-family="Helvetica, Arial, sans-serif" font-size="12" font-weight="400" fill="rgb(20,20,20)">spade 2</text>+</g>+<g aria-label="spade 2 hierarchy">+ <line x1="315.0" y1="101" x2="349.0" y2="101" stroke="rgb(139,0,139)" stroke-width="2.0" stroke-dasharray="7 4"/>+ <path d="M 332.0 96.7 L 327.9 104.5 L 336.1 104.5 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <text x="358.0" y="105.0" text-anchor="start" font-family="Helvetica, Arial, sans-serif" font-size="12" font-weight="400" fill="rgb(20,20,20)">spade 2 hierarchy</text>+</g>+<g aria-label="cdt">+ <line x1="560.0" y1="101" x2="594.0" y2="101" stroke="rgb(178,34,34)" stroke-width="2.0"/>+ <rect x="573.8" y="97.8" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <text x="603.0" y="105.0" text-anchor="start" font-family="Helvetica, Arial, sans-serif" font-size="12" font-weight="400" fill="rgb(20,20,20)">cdt</text>+</g>+<g aria-label="delaunator">+ <line x1="735.0" y1="101" x2="769.0" y2="101" stroke="rgb(0,139,139)" stroke-width="2.0"/>+ <path d="M 752.0 96.8 L 756.2 101.0 L 752.0 105.2 L 747.8 101.0 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <text x="778.0" y="105.0" text-anchor="start" font-family="Helvetica, Arial, sans-serif" font-size="12" font-weight="400" fill="rgb(20,20,20)">delaunator</text>+</g>+<g aria-label="moonlight-triangulation">+ <line x1="940.0" y1="101" x2="974.0" y2="101" stroke="rgb(255,140,0)" stroke-width="3.2"/>+ <path d="M 957.0 95.6 L 958.7 99.3 L 962.4 101.0 L 958.7 102.7 L 957.0 106.4 L 955.3 102.7 L 951.6 101.0 L 955.3 99.3 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <text x="983.0" y="105.0" text-anchor="start" font-family="Helvetica, Arial, sans-serif" font-size="12" font-weight="400" fill="rgb(20,20,20)">moonlight-triangulation</text>+</g>+<g aria-label="local insertion">+<text x="340.0" y="144.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="15" font-weight="500" fill="rgb(20,20,20)">local insertion</text>+<line x1="70.0" y1="610.0" x2="610.0" y2="610.0" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="60.0" y="614.0" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">0.00</text>+<line x1="70.0" y1="521.6" x2="610.0" y2="521.6" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="60.0" y="525.6" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">1.0</text>+<line x1="70.0" y1="433.2" x2="610.0" y2="433.2" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="60.0" y="437.2" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">2.0</text>+<line x1="70.0" y1="344.8" x2="610.0" y2="344.8" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="60.0" y="348.8" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">3.0</text>+<line x1="70.0" y1="256.4" x2="610.0" y2="256.4" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="60.0" y="260.4" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">4.0</text>+<line x1="70.0" y1="168.0" x2="610.0" y2="168.0" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="60.0" y="172.0" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">5.0</text>+<line x1="70.0" y1="610" x2="70.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="70.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">2k</text>+<line x1="160.0" y1="610" x2="160.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="160.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">4k</text>+<line x1="250.0" y1="610" x2="250.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="250.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">6k</text>+<line x1="340.0" y1="610" x2="340.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="340.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">8k</text>+<line x1="430.0" y1="610" x2="430.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="430.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">10k</text>+<line x1="520.0" y1="610" x2="520.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="520.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">12k</text>+<line x1="610.0" y1="610" x2="610.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="610.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">14k</text>+<rect x="70.0" y="168.0" width="540.0" height="442.0" fill="none" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="340.0" y="664.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="12" font-weight="400" fill="rgb(20,20,20)">input size (points)</text>+<text x="15.0" y="389.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="12" fill="rgb(20,20,20)" transform="rotate(-90 15.0 389.0)">CPU time (ms)</text>+<g aria-label="spade 2">+ <polyline points="70.0,583.1 160.0,541.9 250.0,495.1 340.0,453.0 430.0,407.7 520.0,366.9 610.0,314.3" fill="none" stroke="rgb(0,0,139)" stroke-width="2.0" stroke-linecap="round" stroke-linejoin="round"/>+ <line x1="70.0" y1="581.7" x2="70.0" y2="584.5" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="67.0" y1="581.7" x2="73.0" y2="581.7" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="67.0" y1="584.5" x2="73.0" y2="584.5" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="70.0" cy="583.1" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="160.0" y1="537.6" x2="160.0" y2="546.2" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="157.0" y1="537.6" x2="163.0" y2="537.6" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="157.0" y1="546.2" x2="163.0" y2="546.2" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="160.0" cy="541.9" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="250.0" y1="487.9" x2="250.0" y2="502.2" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="247.0" y1="487.9" x2="253.0" y2="487.9" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="247.0" y1="502.2" x2="253.0" y2="502.2" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="250.0" cy="495.1" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="340.0" y1="445.2" x2="340.0" y2="460.9" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="337.0" y1="445.2" x2="343.0" y2="445.2" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="337.0" y1="460.9" x2="343.0" y2="460.9" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="340.0" cy="453.0" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="430.0" y1="395.0" x2="430.0" y2="420.4" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="427.0" y1="395.0" x2="433.0" y2="395.0" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="427.0" y1="420.4" x2="433.0" y2="420.4" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="430.0" cy="407.7" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="520.0" y1="354.2" x2="520.0" y2="379.6" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="517.0" y1="354.2" x2="523.0" y2="354.2" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="517.0" y1="379.6" x2="523.0" y2="379.6" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="520.0" cy="366.9" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="610.0" y1="291.9" x2="610.0" y2="336.8" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="607.0" y1="291.9" x2="613.0" y2="291.9" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="607.0" y1="336.8" x2="613.0" y2="336.8" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="610.0" cy="314.3" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+</g>+<g aria-label="spade 2 hierarchy">+ <polyline points="70.0,582.7 160.0,539.3 250.0,489.3 340.0,441.8 430.0,391.2 520.0,346.0 610.0,292.9" fill="none" stroke="rgb(139,0,139)" stroke-width="2.0" stroke-linecap="round" stroke-linejoin="round" stroke-dasharray="7 4"/>+ <line x1="70.0" y1="580.0" x2="70.0" y2="585.3" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="67.0" y1="580.0" x2="73.0" y2="580.0" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="67.0" y1="585.3" x2="73.0" y2="585.3" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 70.0 578.4 L 65.9 586.2 L 74.1 586.2 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="160.0" y1="532.3" x2="160.0" y2="546.2" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="157.0" y1="532.3" x2="163.0" y2="532.3" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="157.0" y1="546.2" x2="163.0" y2="546.2" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 160.0 535.0 L 155.9 542.8 L 164.1 542.8 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="250.0" y1="482.4" x2="250.0" y2="496.2" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="247.0" y1="482.4" x2="253.0" y2="482.4" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="247.0" y1="496.2" x2="253.0" y2="496.2" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 250.0 485.0 L 245.9 492.8 L 254.1 492.8 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="340.0" y1="430.2" x2="340.0" y2="453.3" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="337.0" y1="430.2" x2="343.0" y2="430.2" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="337.0" y1="453.3" x2="343.0" y2="453.3" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 340.0 437.5 L 335.9 445.3 L 344.1 445.3 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="430.0" y1="381.5" x2="430.0" y2="400.9" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="427.0" y1="381.5" x2="433.0" y2="381.5" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="427.0" y1="400.9" x2="433.0" y2="400.9" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 430.0 386.9 L 425.9 394.7 L 434.1 394.7 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="520.0" y1="324.9" x2="520.0" y2="367.2" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="517.0" y1="324.9" x2="523.0" y2="324.9" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="517.0" y1="367.2" x2="523.0" y2="367.2" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 520.0 341.7 L 515.9 349.5 L 524.1 349.5 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="610.0" y1="274.3" x2="610.0" y2="311.5" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="607.0" y1="274.3" x2="613.0" y2="274.3" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="607.0" y1="311.5" x2="613.0" y2="311.5" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 610.0 288.6 L 605.9 296.4 L 614.1 296.4 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+</g>+<g aria-label="cdt">+ <polyline points="70.0,593.6 160.0,570.7 250.0,538.0 340.0,506.8 430.0,475.6 520.0,445.6 610.0,411.9" fill="none" stroke="rgb(178,34,34)" stroke-width="2.0" stroke-linecap="round" stroke-linejoin="round"/>+ <line x1="70.0" y1="592.9" x2="70.0" y2="594.2" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="67.0" y1="592.9" x2="73.0" y2="592.9" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="67.0" y1="594.2" x2="73.0" y2="594.2" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="66.8" y="590.4" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="160.0" y1="567.9" x2="160.0" y2="573.5" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="157.0" y1="567.9" x2="163.0" y2="567.9" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="157.0" y1="573.5" x2="163.0" y2="573.5" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="156.8" y="567.5" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="250.0" y1="532.4" x2="250.0" y2="543.6" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="247.0" y1="532.4" x2="253.0" y2="532.4" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="247.0" y1="543.6" x2="253.0" y2="543.6" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="246.8" y="534.8" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="340.0" y1="500.2" x2="340.0" y2="513.3" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="337.0" y1="500.2" x2="343.0" y2="500.2" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="337.0" y1="513.3" x2="343.0" y2="513.3" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="336.8" y="503.6" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="430.0" y1="468.6" x2="430.0" y2="482.6" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="427.0" y1="468.6" x2="433.0" y2="468.6" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="427.0" y1="482.6" x2="433.0" y2="482.6" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="426.8" y="472.4" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="520.0" y1="435.9" x2="520.0" y2="455.3" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="517.0" y1="435.9" x2="523.0" y2="435.9" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="517.0" y1="455.3" x2="523.0" y2="455.3" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="516.8" y="442.4" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="610.0" y1="392.8" x2="610.0" y2="430.9" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="607.0" y1="392.8" x2="613.0" y2="392.8" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="607.0" y1="430.9" x2="613.0" y2="430.9" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="606.8" y="408.7" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+</g>+<g aria-label="delaunator">+ <polyline points="70.0,597.4 160.0,580.1 250.0,554.3 340.0,527.1 430.0,495.7 520.0,461.5 610.0,432.1" fill="none" stroke="rgb(0,139,139)" stroke-width="2.0" stroke-linecap="round" stroke-linejoin="round"/>+ <line x1="70.0" y1="596.7" x2="70.0" y2="598.0" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="67.0" y1="596.7" x2="73.0" y2="596.7" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="67.0" y1="598.0" x2="73.0" y2="598.0" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 70.0 593.2 L 74.2 597.4 L 70.0 601.6 L 65.8 597.4 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="160.0" y1="577.8" x2="160.0" y2="582.4" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="157.0" y1="577.8" x2="163.0" y2="577.8" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="157.0" y1="582.4" x2="163.0" y2="582.4" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 160.0 575.9 L 164.2 580.1 L 160.0 584.3 L 155.8 580.1 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="250.0" y1="550.2" x2="250.0" y2="558.5" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="247.0" y1="550.2" x2="253.0" y2="550.2" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="247.0" y1="558.5" x2="253.0" y2="558.5" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 250.0 550.1 L 254.2 554.3 L 250.0 558.5 L 245.8 554.3 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="340.0" y1="522.0" x2="340.0" y2="532.2" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="337.0" y1="522.0" x2="343.0" y2="522.0" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="337.0" y1="532.2" x2="343.0" y2="532.2" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 340.0 522.9 L 344.2 527.1 L 340.0 531.3 L 335.8 527.1 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="430.0" y1="486.9" x2="430.0" y2="504.5" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="427.0" y1="486.9" x2="433.0" y2="486.9" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="427.0" y1="504.5" x2="433.0" y2="504.5" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 430.0 491.5 L 434.2 495.7 L 430.0 499.9 L 425.8 495.7 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="520.0" y1="453.0" x2="520.0" y2="470.1" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="517.0" y1="453.0" x2="523.0" y2="453.0" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="517.0" y1="470.1" x2="523.0" y2="470.1" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 520.0 457.3 L 524.2 461.5 L 520.0 465.7 L 515.8 461.5 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="610.0" y1="418.1" x2="610.0" y2="446.1" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="607.0" y1="418.1" x2="613.0" y2="418.1" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="607.0" y1="446.1" x2="613.0" y2="446.1" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 610.0 427.9 L 614.2 432.1 L 610.0 436.3 L 605.8 432.1 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+</g>+<g aria-label="moonlight-triangulation">+ <polyline points="70.0,585.5 160.0,558.1 250.0,524.2 340.0,490.3 430.0,456.2 520.0,417.5 610.0,361.2" fill="none" stroke="rgb(255,140,0)" stroke-width="3.2" stroke-linecap="round" stroke-linejoin="round"/>+ <line x1="70.0" y1="583.9" x2="70.0" y2="587.1" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="67.0" y1="583.9" x2="73.0" y2="583.9" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="67.0" y1="587.1" x2="73.0" y2="587.1" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 70.0 580.1 L 71.7 583.8 L 75.4 585.5 L 71.7 587.2 L 70.0 590.9 L 68.3 587.2 L 64.6 585.5 L 68.3 583.8 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="160.0" y1="554.3" x2="160.0" y2="561.9" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="157.0" y1="554.3" x2="163.0" y2="554.3" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="157.0" y1="561.9" x2="163.0" y2="561.9" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 160.0 552.7 L 161.7 556.4 L 165.4 558.1 L 161.7 559.8 L 160.0 563.5 L 158.3 559.8 L 154.6 558.1 L 158.3 556.4 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="250.0" y1="518.4" x2="250.0" y2="529.9" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="247.0" y1="518.4" x2="253.0" y2="518.4" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="247.0" y1="529.9" x2="253.0" y2="529.9" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 250.0 518.8 L 251.7 522.5 L 255.4 524.2 L 251.7 525.9 L 250.0 529.6 L 248.3 525.9 L 244.6 524.2 L 248.3 522.5 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="340.0" y1="480.9" x2="340.0" y2="499.6" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="337.0" y1="480.9" x2="343.0" y2="480.9" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="337.0" y1="499.6" x2="343.0" y2="499.6" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 340.0 484.9 L 341.7 488.6 L 345.4 490.3 L 341.7 492.0 L 340.0 495.7 L 338.3 492.0 L 334.6 490.3 L 338.3 488.6 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="430.0" y1="449.5" x2="430.0" y2="462.9" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="427.0" y1="449.5" x2="433.0" y2="449.5" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="427.0" y1="462.9" x2="433.0" y2="462.9" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 430.0 450.8 L 431.7 454.5 L 435.4 456.2 L 431.7 457.9 L 430.0 461.6 L 428.3 457.9 L 424.6 456.2 L 428.3 454.5 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="520.0" y1="398.7" x2="520.0" y2="436.3" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="517.0" y1="398.7" x2="523.0" y2="398.7" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="517.0" y1="436.3" x2="523.0" y2="436.3" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 520.0 412.1 L 521.7 415.8 L 525.4 417.5 L 521.7 419.2 L 520.0 422.9 L 518.3 419.2 L 514.6 417.5 L 518.3 415.8 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="610.0" y1="346.9" x2="610.0" y2="375.6" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="607.0" y1="346.9" x2="613.0" y2="346.9" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="607.0" y1="375.6" x2="613.0" y2="375.6" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 610.0 355.8 L 611.7 359.5 L 615.4 361.2 L 611.7 362.9 L 610.0 366.6 L 608.3 362.9 L 604.6 361.2 L 608.3 359.5 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+</g>+</g>+<g aria-label="uniform">+<text x="940.0" y="144.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="15" font-weight="500" fill="rgb(20,20,20)">uniform</text>+<line x1="670.0" y1="610.0" x2="1210.0" y2="610.0" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="660.0" y="614.0" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">0.00</text>+<line x1="670.0" y1="521.6" x2="1210.0" y2="521.6" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="660.0" y="525.6" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">1.0</text>+<line x1="670.0" y1="433.2" x2="1210.0" y2="433.2" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="660.0" y="437.2" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">2.0</text>+<line x1="670.0" y1="344.8" x2="1210.0" y2="344.8" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="660.0" y="348.8" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">3.0</text>+<line x1="670.0" y1="256.4" x2="1210.0" y2="256.4" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="660.0" y="260.4" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">4.0</text>+<line x1="670.0" y1="168.0" x2="1210.0" y2="168.0" stroke="rgb(170,170,170)" stroke-width="0.7" stroke-dasharray="2 4"/>+<text x="660.0" y="172.0" text-anchor="end" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">5.0</text>+<line x1="670.0" y1="610" x2="670.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="670.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">2k</text>+<line x1="760.0" y1="610" x2="760.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="760.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">4k</text>+<line x1="850.0" y1="610" x2="850.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="850.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">6k</text>+<line x1="940.0" y1="610" x2="940.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="940.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">8k</text>+<line x1="1030.0" y1="610" x2="1030.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="1030.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">10k</text>+<line x1="1120.0" y1="610" x2="1120.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="1120.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">12k</text>+<line x1="1210.0" y1="610" x2="1210.0" y2="616" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="1210.0" y="635.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">14k</text>+<rect x="670.0" y="168.0" width="540.0" height="442.0" fill="none" stroke="rgb(25,25,25)" stroke-width="1"/>+<text x="940.0" y="664.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="12" font-weight="400" fill="rgb(20,20,20)">input size (points)</text>+<text x="615.0" y="389.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="12" fill="rgb(20,20,20)" transform="rotate(-90 615.0 389.0)">CPU time (ms)</text>+<g aria-label="spade 2">+ <polyline points="670.0,583.7 760.0,547.6 850.0,506.8 940.0,466.9 1030.0,426.3 1120.0,383.7 1210.0,338.9" fill="none" stroke="rgb(0,0,139)" stroke-width="2.0" stroke-linecap="round" stroke-linejoin="round"/>+ <line x1="670.0" y1="582.4" x2="670.0" y2="585.0" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="667.0" y1="582.4" x2="673.0" y2="582.4" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="667.0" y1="585.0" x2="673.0" y2="585.0" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="670.0" cy="583.7" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="760.0" y1="542.4" x2="760.0" y2="552.8" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="757.0" y1="542.4" x2="763.0" y2="542.4" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="757.0" y1="552.8" x2="763.0" y2="552.8" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="760.0" cy="547.6" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="850.0" y1="497.7" x2="850.0" y2="515.8" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="847.0" y1="497.7" x2="853.0" y2="497.7" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="847.0" y1="515.8" x2="853.0" y2="515.8" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="850.0" cy="506.8" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="940.0" y1="463.1" x2="940.0" y2="470.7" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="937.0" y1="463.1" x2="943.0" y2="463.1" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="937.0" y1="470.7" x2="943.0" y2="470.7" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="940.0" cy="466.9" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="1030.0" y1="415.5" x2="1030.0" y2="437.1" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1027.0" y1="415.5" x2="1033.0" y2="415.5" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1027.0" y1="437.1" x2="1033.0" y2="437.1" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="1030.0" cy="426.3" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="1120.0" y1="370.6" x2="1120.0" y2="396.9" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1117.0" y1="370.6" x2="1123.0" y2="370.6" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1117.0" y1="396.9" x2="1123.0" y2="396.9" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="1120.0" cy="383.7" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+ <line x1="1210.0" y1="317.9" x2="1210.0" y2="359.9" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1207.0" y1="317.9" x2="1213.0" y2="317.9" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1207.0" y1="359.9" x2="1213.0" y2="359.9" stroke="rgb(0,0,139)" stroke-width="0.9" opacity="0.48"/>+ <circle cx="1210.0" cy="338.9" r="3.2" fill="white" stroke="rgb(0,0,139)" stroke-width="1.6"/>+</g>+<g aria-label="spade 2 hierarchy">+ <polyline points="670.0,583.3 760.0,544.5 850.0,499.8 940.0,455.4 1030.0,406.4 1120.0,364.4 1210.0,316.7" fill="none" stroke="rgb(139,0,139)" stroke-width="2.0" stroke-linecap="round" stroke-linejoin="round" stroke-dasharray="7 4"/>+ <line x1="670.0" y1="582.1" x2="670.0" y2="584.6" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="667.0" y1="582.1" x2="673.0" y2="582.1" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="667.0" y1="584.6" x2="673.0" y2="584.6" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 670.0 579.0 L 665.9 586.8 L 674.1 586.8 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="760.0" y1="541.1" x2="760.0" y2="547.8" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="757.0" y1="541.1" x2="763.0" y2="541.1" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="757.0" y1="547.8" x2="763.0" y2="547.8" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 760.0 540.2 L 755.9 548.0 L 764.1 548.0 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="850.0" y1="489.8" x2="850.0" y2="509.8" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="847.0" y1="489.8" x2="853.0" y2="489.8" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="847.0" y1="509.8" x2="853.0" y2="509.8" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 850.0 495.5 L 845.9 503.3 L 854.1 503.3 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="940.0" y1="443.6" x2="940.0" y2="467.2" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="937.0" y1="443.6" x2="943.0" y2="443.6" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="937.0" y1="467.2" x2="943.0" y2="467.2" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 940.0 451.1 L 935.9 458.9 L 944.1 458.9 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="1030.0" y1="387.9" x2="1030.0" y2="424.9" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1027.0" y1="387.9" x2="1033.0" y2="387.9" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1027.0" y1="424.9" x2="1033.0" y2="424.9" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 1030.0 402.1 L 1025.9 409.9 L 1034.1 409.9 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="1120.0" y1="344.7" x2="1120.0" y2="384.0" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1117.0" y1="344.7" x2="1123.0" y2="344.7" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1117.0" y1="384.0" x2="1123.0" y2="384.0" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 1120.0 360.1 L 1115.9 367.9 L 1124.1 367.9 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+ <line x1="1210.0" y1="306.8" x2="1210.0" y2="326.6" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1207.0" y1="306.8" x2="1213.0" y2="306.8" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1207.0" y1="326.6" x2="1213.0" y2="326.6" stroke="rgb(139,0,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 1210.0 312.4 L 1205.9 320.2 L 1214.1 320.2 Z" fill="white" stroke="rgb(139,0,139)" stroke-width="1.6"/>+</g>+<g aria-label="cdt">+ <polyline points="670.0,593.7 760.0,570.4 850.0,534.1 940.0,500.7 1030.0,465.7 1120.0,431.7 1210.0,396.1" fill="none" stroke="rgb(178,34,34)" stroke-width="2.0" stroke-linecap="round" stroke-linejoin="round"/>+ <line x1="670.0" y1="592.7" x2="670.0" y2="594.6" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="667.0" y1="592.7" x2="673.0" y2="592.7" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="667.0" y1="594.6" x2="673.0" y2="594.6" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="666.8" y="590.5" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="760.0" y1="569.0" x2="760.0" y2="571.7" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="757.0" y1="569.0" x2="763.0" y2="569.0" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="757.0" y1="571.7" x2="763.0" y2="571.7" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="756.8" y="567.2" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="850.0" y1="529.0" x2="850.0" y2="539.3" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="847.0" y1="529.0" x2="853.0" y2="529.0" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="847.0" y1="539.3" x2="853.0" y2="539.3" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="846.8" y="530.9" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="940.0" y1="491.0" x2="940.0" y2="510.3" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="937.0" y1="491.0" x2="943.0" y2="491.0" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="937.0" y1="510.3" x2="943.0" y2="510.3" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="936.8" y="497.5" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="1030.0" y1="455.0" x2="1030.0" y2="476.4" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1027.0" y1="455.0" x2="1033.0" y2="455.0" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1027.0" y1="476.4" x2="1033.0" y2="476.4" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="1026.8" y="462.5" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="1120.0" y1="420.1" x2="1120.0" y2="443.4" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1117.0" y1="420.1" x2="1123.0" y2="420.1" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1117.0" y1="443.4" x2="1123.0" y2="443.4" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="1116.8" y="428.5" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+ <line x1="1210.0" y1="379.4" x2="1210.0" y2="412.8" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1207.0" y1="379.4" x2="1213.0" y2="379.4" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1207.0" y1="412.8" x2="1213.0" y2="412.8" stroke="rgb(178,34,34)" stroke-width="0.9" opacity="0.48"/>+ <rect x="1206.8" y="392.9" width="6.4" height="6.4" fill="white" stroke="rgb(178,34,34)" stroke-width="1.6"/>+</g>+<g aria-label="delaunator">+ <polyline points="670.0,597.2 760.0,580.7 850.0,558.9 940.0,530.2 1030.0,497.9 1120.0,466.9 1210.0,434.3" fill="none" stroke="rgb(0,139,139)" stroke-width="2.0" stroke-linecap="round" stroke-linejoin="round"/>+ <line x1="670.0" y1="596.0" x2="670.0" y2="598.5" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="667.0" y1="596.0" x2="673.0" y2="596.0" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="667.0" y1="598.5" x2="673.0" y2="598.5" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 670.0 593.0 L 674.2 597.2 L 670.0 601.4 L 665.8 597.2 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="760.0" y1="578.8" x2="760.0" y2="582.6" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="757.0" y1="578.8" x2="763.0" y2="578.8" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="757.0" y1="582.6" x2="763.0" y2="582.6" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 760.0 576.5 L 764.2 580.7 L 760.0 584.9 L 755.8 580.7 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="850.0" y1="554.6" x2="850.0" y2="563.2" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="847.0" y1="554.6" x2="853.0" y2="554.6" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="847.0" y1="563.2" x2="853.0" y2="563.2" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 850.0 554.7 L 854.2 558.9 L 850.0 563.1 L 845.8 558.9 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="940.0" y1="529.0" x2="940.0" y2="531.5" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="937.0" y1="529.0" x2="943.0" y2="529.0" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="937.0" y1="531.5" x2="943.0" y2="531.5" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 940.0 526.0 L 944.2 530.2 L 940.0 534.4 L 935.8 530.2 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="1030.0" y1="487.6" x2="1030.0" y2="508.3" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1027.0" y1="487.6" x2="1033.0" y2="487.6" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1027.0" y1="508.3" x2="1033.0" y2="508.3" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 1030.0 493.7 L 1034.2 497.9 L 1030.0 502.1 L 1025.8 497.9 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="1120.0" y1="461.4" x2="1120.0" y2="472.4" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1117.0" y1="461.4" x2="1123.0" y2="461.4" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1117.0" y1="472.4" x2="1123.0" y2="472.4" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 1120.0 462.7 L 1124.2 466.9 L 1120.0 471.1 L 1115.8 466.9 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+ <line x1="1210.0" y1="419.9" x2="1210.0" y2="448.7" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1207.0" y1="419.9" x2="1213.0" y2="419.9" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1207.0" y1="448.7" x2="1213.0" y2="448.7" stroke="rgb(0,139,139)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 1210.0 430.1 L 1214.2 434.3 L 1210.0 438.5 L 1205.8 434.3 Z" fill="white" stroke="rgb(0,139,139)" stroke-width="1.6"/>+</g>+<g aria-label="moonlight-triangulation">+ <polyline points="670.0,585.5 760.0,561.2 850.0,530.7 940.0,496.8 1030.0,462.3 1120.0,419.3 1210.0,374.6" fill="none" stroke="rgb(255,140,0)" stroke-width="3.2" stroke-linecap="round" stroke-linejoin="round"/>+ <line x1="670.0" y1="584.0" x2="670.0" y2="587.0" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="667.0" y1="584.0" x2="673.0" y2="584.0" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="667.0" y1="587.0" x2="673.0" y2="587.0" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 670.0 580.1 L 671.7 583.8 L 675.4 585.5 L 671.7 587.2 L 670.0 590.9 L 668.3 587.2 L 664.6 585.5 L 668.3 583.8 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="760.0" y1="560.5" x2="760.0" y2="562.0" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="757.0" y1="560.5" x2="763.0" y2="560.5" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="757.0" y1="562.0" x2="763.0" y2="562.0" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 760.0 555.8 L 761.7 559.5 L 765.4 561.2 L 761.7 562.9 L 760.0 566.6 L 758.3 562.9 L 754.6 561.2 L 758.3 559.5 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="850.0" y1="522.8" x2="850.0" y2="538.6" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="847.0" y1="522.8" x2="853.0" y2="522.8" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="847.0" y1="538.6" x2="853.0" y2="538.6" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 850.0 525.3 L 851.7 529.0 L 855.4 530.7 L 851.7 532.4 L 850.0 536.1 L 848.3 532.4 L 844.6 530.7 L 848.3 529.0 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="940.0" y1="485.9" x2="940.0" y2="507.6" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="937.0" y1="485.9" x2="943.0" y2="485.9" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="937.0" y1="507.6" x2="943.0" y2="507.6" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 940.0 491.4 L 941.7 495.1 L 945.4 496.8 L 941.7 498.5 L 940.0 502.2 L 938.3 498.5 L 934.6 496.8 L 938.3 495.1 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="1030.0" y1="451.9" x2="1030.0" y2="472.7" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1027.0" y1="451.9" x2="1033.0" y2="451.9" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1027.0" y1="472.7" x2="1033.0" y2="472.7" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 1030.0 456.9 L 1031.7 460.6 L 1035.4 462.3 L 1031.7 464.0 L 1030.0 467.7 L 1028.3 464.0 L 1024.6 462.3 L 1028.3 460.6 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="1120.0" y1="401.2" x2="1120.0" y2="437.4" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1117.0" y1="401.2" x2="1123.0" y2="401.2" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1117.0" y1="437.4" x2="1123.0" y2="437.4" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 1120.0 413.9 L 1121.7 417.6 L 1125.4 419.3 L 1121.7 421.0 L 1120.0 424.7 L 1118.3 421.0 L 1114.6 419.3 L 1118.3 417.6 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+ <line x1="1210.0" y1="364.2" x2="1210.0" y2="385.0" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1207.0" y1="364.2" x2="1213.0" y2="364.2" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <line x1="1207.0" y1="385.0" x2="1213.0" y2="385.0" stroke="rgb(255,140,0)" stroke-width="0.9" opacity="0.48"/>+ <path d="M 1210.0 369.2 L 1211.7 372.9 L 1215.4 374.6 L 1211.7 376.3 L 1210.0 380.0 L 1208.3 376.3 L 1204.6 374.6 L 1208.3 372.9 Z" fill="rgb(255,140,0)" stroke="white" stroke-width="1"/>+</g>+</g>+<text x="640.0" y="704.0" text-anchor="middle" font-family="Helvetica, Arial, sans-serif" font-size="11" font-weight="400" fill="rgb(20,20,20)">Whiskers show tasty-bench's reported 2× standard-deviation spread.</text>+</svg>
bench/support/BenchSupport.hs view
@@ -6,6 +6,7 @@ -- around a timing distribution cannot express either. module BenchSupport ( timedValue+ , timedProjection , requireRight , randomPoints , latticePoints@@ -24,13 +25,15 @@ , getRTSStats , getRTSStatsEnabled )-import Moonlight.Triangulation- ( FaceId- , Point (Point)+import Moonlight.Triangulation.Dcel+ ( faceVertices+ , vertexPoint+ )+import Moonlight.Triangulation.Handles.HandleDefs (FaceId)+import Moonlight.Triangulation.Types+ ( Point (Point) , Triangulation- , faceVertices , pointX- , vertexPoint ) import System.CPUTime (getCPUTime) @@ -52,12 +55,20 @@ -- as the watermark up to here, and only against a run with the same GC -- schedule; @-G1@ is what makes two schedules comparable. timedValue :: NFData value => String -> IO value -> IO value-timedValue label action = do+timedValue label = timedProjection label id++-- | Time an action while forcing only its authoritative observation. This is+-- for opaque prepared carriers whose public projections have 'NFData'+-- instances even when the carrier itself deliberately does not. The original+-- value is returned for subsequent benchmark lanes.+timedProjection :: NFData observation => String -> (value -> observation) -> IO value -> IO value+timedProjection label observe action = do statsEnabled <- getRTSStatsEnabled allocationStart <- getAllocationCounter wallStart <- getMonotonicTimeNSec cpuStart <- getCPUTime- value <- action >>= evaluate . force+ value <- action+ _ <- evaluate (force (observe value)) cpuEnd <- getCPUTime wallEnd <- getMonotonicTimeNSec allocationEnd <- getAllocationCounter
bindings/README.md view
@@ -44,6 +44,10 @@ remain valid until explicitly freed; every operation publishes a fresh output handle. `ml_runtime_initialize` is idempotent and process-lifetime. +Python and Rust owned handles are thread-affine: create, use, and close each+handle on one thread. This prevents explicit close from racing a foreign call+after the runtime has received the raw handle.+ Batch related mesh edits into one `insert_many` call. Author a whole region with one component list; there is deliberately no mutable builder or per-vertex FFI surface.
bindings/python/src/moonlight_triangulation/__init__.py view
@@ -9,6 +9,7 @@ from enum import IntEnum from fractions import Fraction from pathlib import Path+from threading import get_ident from typing import Final, Self @@ -212,7 +213,7 @@ handle = _produce_handle( self._native, lambda output, obstruction: self._native.library.ml_delaunay_f64(- pointer, len(points), output, obstruction+ pointer, len(coordinates) // 2, output, obstruction ), ) return Mesh(self._native, handle)@@ -228,7 +229,7 @@ self._native, lambda output, obstruction: self._native.library.ml_region_create_f64( coordinate_pointer,- sum(loop_counts),+ len(coordinates) // 2, loop_buffer, len(loop_counts), component_buffer,@@ -244,22 +245,24 @@ handle = _produce_handle( self._native, lambda output, obstruction: self._native.library.ml_structuring_element_create_f64(- pointer, len(points), output, obstruction+ pointer, len(coordinates) // 2, output, obstruction ), ) return StructuringElement(self._native, handle) class _OwnedHandle:- __slots__ = ("_native", "_handle", "_finalizer", "_kind", "__weakref__")+ __slots__ = ("_native", "_handle", "_finalizer", "_kind", "_owner_thread", "__weakref__") def __init__(self, native: _NativeApi, handle: ctypes.c_void_p, free: Callable[[ctypes.c_void_p], None], kind: str) -> None: self._native = native self._handle = handle self._kind = kind+ self._owner_thread = get_ident() self._finalizer = weakref.finalize(self, free, handle) def close(self) -> None:+ self._require_owner_thread() self._finalizer() self._handle = ctypes.c_void_p() @@ -270,10 +273,15 @@ self.close() def _live_handle(self) -> ctypes.c_void_p:+ self._require_owner_thread() if not self._finalizer.alive: raise RuntimeError(f"{self._kind} is closed") return self._handle + def _require_owner_thread(self) -> None:+ if get_ident() != self._owner_thread:+ raise RuntimeError(f"{self._kind} belongs to another thread")+ def _require_same_runtime(self, other: _OwnedHandle) -> None: if self._native is not other._native: raise ValueError("both values must belong to the same Moonlight runtime")@@ -329,7 +337,7 @@ handle = _produce_handle( self._native, lambda output, obstruction: self._native.library.ml_mesh_insert_many_f64(- self._live_handle(), pointer, len(points), output, obstruction+ self._live_handle(), pointer, len(coordinates) // 2, output, obstruction ), ) return Mesh(self._native, handle)
bindings/python/tests/test_binding.py view
@@ -3,6 +3,8 @@ import math import os import unittest+from collections.abc import Sequence+from concurrent.futures import ThreadPoolExecutor from fractions import Fraction from typing import ClassVar @@ -15,6 +17,16 @@ ) +class LengthInflatingPoints(Sequence[tuple[float, float]]):+ def __len__(self) -> int:+ return 1_000++ def __getitem__(self, index: int) -> tuple[float, float]:+ if index == 0:+ return (0.0, 0.0)+ raise IndexError(index)++ class MoonlightBindingTest(unittest.TestCase): engine: ClassVar[Moonlight] @@ -83,6 +95,20 @@ self.assertEqual(raised.exception.code, 1) self.assertEqual(raised.exception.coordinate_error, 1) self.assertEqual(raised.exception.input_index, 1)++ def test_native_count_descends_from_materialized_coordinates(self) -> None:+ mesh = self.engine.delaunay(LengthInflatingPoints())+ self.addCleanup(mesh.close)+ self.assertEqual(mesh.vertex_count, 1)++ def test_owned_handles_are_thread_affine(self) -> None:+ mesh = self.engine.delaunay([(0, 0), (1, 0), (0, 1)])+ self.addCleanup(mesh.close)+ with ThreadPoolExecutor(max_workers=1) as executor:+ foreign_access = executor.submit(lambda: mesh.vertex_count)+ with self.assertRaisesRegex(RuntimeError, "belongs to another thread"):+ foreign_access.result()+ self.assertEqual(mesh.vertex_count, 3) if __name__ == "__main__":
bindings/typescript/src/index.ts view
@@ -3,6 +3,8 @@ const STATUS_OK = 0; const STATUS_BUFFER_TOO_SMALL = 3; const ABI_VERSION = 2;+const bindingCapability: unique symbol = Symbol("moonlight binding capability");+type BindingCapability = typeof bindingCapability; export type Point = readonly [x: number, y: number]; export type Triangle = readonly [first: number, second: number, third: number];@@ -163,6 +165,12 @@ const handleFinalizer = new FinalizationRegistry<FinalizerState>(({ free, handle }) => free(handle)); +function requireBindingCapability(capability: BindingCapability): void {+ if (capability !== bindingCapability) {+ throw new Error("native binding capability is unavailable");+ }+}+ export class MoonlightError extends Error { readonly status: number; readonly code: number;@@ -211,9 +219,10 @@ delaunay(points: readonly Point[]): Mesh { const coordinates = flattenPoints(points); return new Mesh(+ bindingCapability, this.#native, produceHandle((output, obstruction) =>- this.#native.delaunay(coordinates, points.length, output, obstruction),+ this.#native.delaunay(coordinates, coordinates.length / 2, output, obstruction), ), ); }@@ -226,15 +235,16 @@ ); const coordinates = flattenLoops(loops); return new Region(+ bindingCapability, this.#native, produceHandle((output, obstruction) => this.#native.regionCreate( coordinates, coordinates.length / 2, loopPointCounts,- loops.length,+ loopPointCounts.length, componentLoopCounts,- components.length,+ componentLoopCounts.length, output, obstruction, ),@@ -245,9 +255,10 @@ structuringElement(points: readonly Point[]): StructuringElement { const coordinates = flattenPoints(points); return new StructuringElement(+ bindingCapability, this.#native, produceHandle((output, obstruction) =>- this.#native.structuringElementCreate(coordinates, points.length, output, obstruction),+ this.#native.structuringElementCreate(coordinates, coordinates.length / 2, output, obstruction), ), ); }@@ -260,11 +271,13 @@ #handle: NativeHandle | null; protected constructor(+ capability: BindingCapability, native: NativeApi, handle: NativeHandle, free: (handle: NativeHandle) => void, kind: string, ) {+ requireBindingCapability(capability); this.#native = native; this.#free = free; this.#kind = kind;@@ -281,11 +294,13 @@ } } - native(): NativeApi {+ protected native(capability: BindingCapability): NativeApi {+ requireBindingCapability(capability); return this.#native; } - nativeHandle(): NativeHandle {+ nativeHandle(capability: BindingCapability): NativeHandle {+ requireBindingCapability(capability); if (this.#handle === null) { throw new Error(`${this.#kind} is closed`); }@@ -300,16 +315,16 @@ } export class Mesh extends OwnedHandle {- constructor(native: NativeApi, handle: NativeHandle) {- super(native, handle, native.meshFree, "mesh");+ constructor(capability: BindingCapability, native: NativeApi, handle: NativeHandle) {+ super(capability, native, handle, native.meshFree, "mesh"); } vertexCount(): number {- return this.count(this.native().vertexCount);+ return this.count(this.native(bindingCapability).vertexCount); } triangleCount(): number {- return this.count(this.native().triangleCount);+ return this.count(this.native(bindingCapability).triangleCount); } vertices(): readonly Point[] {@@ -317,7 +332,7 @@ const coordinates = new Float64Array(count * 2); const written: NativeInteger[] = [0]; const obstruction: NativeObstruction = {};- const status = this.native().copyVertices(this.nativeHandle(), coordinates, count, written, obstruction);+ const status = this.native(bindingCapability).copyVertices(this.nativeHandle(bindingCapability), coordinates, count, written, obstruction); checkStatus(status, obstruction); return Array.from({ length: toSafeNumber(written[0]) }, (_unused, index): Point => [ coordinates[index * 2] ?? 0,@@ -330,7 +345,7 @@ const triangles = new Uint32Array(count * 3); const written: NativeInteger[] = [0]; const obstruction: NativeObstruction = {};- const status = this.native().copyTriangles(this.nativeHandle(), triangles, count, written, obstruction);+ const status = this.native(bindingCapability).copyTriangles(this.nativeHandle(bindingCapability), triangles, count, written, obstruction); checkStatus(status, obstruction); return Array.from({ length: toSafeNumber(written[0]) }, (_unused, index): Triangle => [ triangles[index * 3] ?? 0,@@ -342,35 +357,37 @@ insertMany(points: readonly Point[]): Mesh { const coordinates = flattenPoints(points); return new Mesh(- this.native(),+ bindingCapability,+ this.native(bindingCapability), produceHandle((output, obstruction) =>- this.native().insertMany(this.nativeHandle(), coordinates, points.length, output, obstruction),+ this.native(bindingCapability).insertMany(this.nativeHandle(bindingCapability), coordinates, coordinates.length / 2, output, obstruction), ), ); } siteUnion(other: Mesh): Mesh {- return this.binary(other, this.native().meshSiteUnion);+ return this.binary(other, this.native(bindingCapability).meshSiteUnion); } siteIntersection(other: Mesh): Mesh {- return this.binary(other, this.native().meshSiteIntersection);+ return this.binary(other, this.native(bindingCapability).meshSiteIntersection); } siteDifference(other: Mesh): Mesh {- return this.binary(other, this.native().meshSiteDifference);+ return this.binary(other, this.native(bindingCapability).meshSiteDifference); } siteSymmetricDifference(other: Mesh): Mesh {- return this.binary(other, this.native().meshSiteSymmetricDifference);+ return this.binary(other, this.native(bindingCapability).meshSiteSymmetricDifference); } private binary(other: Mesh, operation: BinaryOperation): Mesh { this.requireSameRuntime(other); return new Mesh(- this.native(),+ bindingCapability,+ this.native(bindingCapability), produceHandle((output, obstruction) =>- operation(this.nativeHandle(), other.nativeHandle(), output, obstruction),+ operation(this.nativeHandle(bindingCapability), other.nativeHandle(bindingCapability), output, obstruction), ), ); }@@ -378,15 +395,15 @@ private count(operation: CountOperation): number { const output: NativeInteger[] = [0]; const obstruction: NativeObstruction = {};- const status = operation(this.nativeHandle(), output, obstruction);+ const status = operation(this.nativeHandle(bindingCapability), output, obstruction); checkStatus(status, obstruction); return toSafeNumber(output[0]); } } export class Region extends OwnedHandle {- constructor(native: NativeApi, handle: NativeHandle) {- super(native, handle, native.regionFree, "region");+ constructor(capability: BindingCapability, native: NativeApi, handle: NativeHandle) {+ super(capability, native, handle, native.regionFree, "region"); } components(): readonly PolygonComponent[] {@@ -395,8 +412,8 @@ const loopPointOffsets = new BigUint64Array(loopCount + 1); const componentLoopOffsets = new BigUint64Array(componentCount + 1); const obstruction: NativeObstruction = {};- const status = this.native().regionCopy(- this.nativeHandle(),+ const status = this.native(bindingCapability).regionCopy(+ this.nativeHandle(bindingCapability), coordinates, pointCount, loopPointOffsets,@@ -431,54 +448,59 @@ locate([x, y]: Point): RegionLocation { const location = [0]; const obstruction: NativeObstruction = {};- const status = this.native().regionLocate(this.nativeHandle(), x, y, location, obstruction);+ const status = this.native(bindingCapability).regionLocate(this.nativeHandle(bindingCapability), x, y, location, obstruction); checkStatus(status, obstruction); return locationCode(location[0]); } union(other: Region): Region {- return this.binary(other, this.native().regionUnion);+ return this.binary(other, this.native(bindingCapability).regionUnion); } intersection(other: Region): Region {- return this.binary(other, this.native().regionIntersection);+ return this.binary(other, this.native(bindingCapability).regionIntersection); } difference(other: Region): Region {- return this.binary(other, this.native().regionDifference);+ return this.binary(other, this.native(bindingCapability).regionDifference); } symmetricDifference(other: Region): Region {- return this.binary(other, this.native().regionSymmetricDifference);+ return this.binary(other, this.native(bindingCapability).regionSymmetricDifference); } minkowskiSum(other: Region): readonly [Region, MinkowskiReceipt] { this.requireSameRuntime(other);- return this.morph(this.native().regionMinkowskiSum, this.nativeHandle(), other.nativeHandle());+ return this.morph(+ this.native(bindingCapability).regionMinkowskiSum,+ this.nativeHandle(bindingCapability),+ other.nativeHandle(bindingCapability),+ ); } offset(element: StructuringElement): readonly [Region, MinkowskiReceipt] {- return this.withElement(element, this.native().regionOffset);+ return this.withElement(element, this.native(bindingCapability).regionOffset); } inset(element: StructuringElement): readonly [Region, MinkowskiReceipt] {- return this.withElement(element, this.native().regionInset);+ return this.withElement(element, this.native(bindingCapability).regionInset); } open(element: StructuringElement): readonly [Region, MinkowskiReceipt] {- return this.withElement(element, this.native().regionOpen);+ return this.withElement(element, this.native(bindingCapability).regionOpen); } closeWith(element: StructuringElement): readonly [Region, MinkowskiReceipt] {- return this.withElement(element, this.native().regionClose);+ return this.withElement(element, this.native(bindingCapability).regionClose); } private binary(other: Region, operation: BinaryOperation): Region { this.requireSameRuntime(other); return new Region(- this.native(),+ bindingCapability,+ this.native(bindingCapability), produceHandle((output, obstruction) =>- operation(this.nativeHandle(), other.nativeHandle(), output, obstruction),+ operation(this.nativeHandle(bindingCapability), other.nativeHandle(bindingCapability), output, obstruction), ), ); }@@ -488,7 +510,11 @@ operation: MorphologyOperation, ): readonly [Region, MinkowskiReceipt] { this.requireSameRuntime(element);- return this.morph(operation, element.nativeHandle(), this.nativeHandle());+ return this.morph(+ operation,+ element.nativeHandle(bindingCapability),+ this.nativeHandle(bindingCapability),+ ); } private morph(@@ -499,7 +525,7 @@ const [handle, receipt] = produceMorphology((output, nativeReceipt, obstruction) => operation(first, second, output, nativeReceipt, obstruction), );- return [new Region(this.native(), handle), receipt];+ return [new Region(bindingCapability, this.native(bindingCapability), handle), receipt]; } private counts(): readonly [number, number, number] {@@ -507,8 +533,8 @@ const loopCount: NativeInteger[] = [0]; const pointCount: NativeInteger[] = [0]; const obstruction: NativeObstruction = {};- const status = this.native().regionCounts(- this.nativeHandle(),+ const status = this.native(bindingCapability).regionCounts(+ this.nativeHandle(bindingCapability), componentCount, loopCount, pointCount,@@ -525,8 +551,8 @@ const perimeterLower = [0]; const perimeterUpper = [0]; const obstruction: NativeObstruction = {};- const status = this.native().regionMeasure(- this.nativeHandle(),+ const status = this.native(bindingCapability).regionMeasure(+ this.nativeHandle(bindingCapability), euler, areaRatio, capacity,@@ -553,8 +579,8 @@ } export class StructuringElement extends OwnedHandle {- constructor(native: NativeApi, handle: NativeHandle) {- super(native, handle, native.structuringElementFree, "structuring element");+ constructor(capability: BindingCapability, native: NativeApi, handle: NativeHandle) {+ super(capability, native, handle, native.structuringElementFree, "structuring element"); } }
bindings/typescript/test/binding.test.ts view
@@ -76,3 +76,28 @@ failure.inputIndex === 1n, ); });++test("native count descends from materialized coordinates", () => {+ const engine = new Moonlight();+ const points = new Proxy<readonly (readonly [number, number])[]>([[0, 0]], {+ get(target, property, receiver) {+ return property === "length" ? 1_000 : Reflect.get(target, property, receiver);+ },+ });+ const mesh = engine.delaunay(points);+ assert.equal(mesh.vertexCount(), 1);+ mesh.close();+});++test("raw native handles require the module capability", () => {+ const engine = new Moonlight();+ const mesh = engine.delaunay([[0, 0]]);+ const capabilityGuarded = mesh as unknown as {+ nativeHandle(capability: symbol): object;+ };+ assert.throws(+ () => capabilityGuarded.nativeHandle(Symbol("forged capability")),+ /native binding capability is unavailable/,+ );+ mesh.close();+});
moonlight-triangulation.cabal view
@@ -1,6 +1,6 @@ cabal-version: 3.4 name: moonlight-triangulation-version: 1.3.0.0+version: 1.3.0.2 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@@ -20,13 +20,15 @@ hints; a bounded-concurrency interpreter for the join tournament; and a versioned binary serialization surface. On GHC 9.14, a public cell-complex component interprets- admitted exact cell selections for Homology and Category.+ admitted exact cell selections for Homology and Category+ and lowers exact Delaunay alpha filtrations into persistent+ homology. Failure is values: every refusal names its witness. license: MIT license-file: LICENSE-author: Fable, Blue Rose+author: Blue Rose maintainer: rosaliafialkova@gmail.com-copyright: (c) 2026 Fable, Blue Rose+copyright: (c) 2026 Blue Rose category: Geometry, Math homepage: https://github.com/PaleRoses/moonlight bug-reports: https://github.com/PaleRoses/moonlight/issues@@ -38,11 +40,13 @@ extra-doc-files: README.md CHANGELOG.md+ bench/delaunay-compare/results/moonlight-delaunay-compare-small.svg extra-source-files: weeder.toml bench/delaunay-compare/README.md- -- Receipts and their SVG projections remain beside the benchmark, but are- -- not package inputs: taking a new measurement must not rebuild the library.+ -- The retained CSV, source receipt, and big-point SVG remain beside the+ -- benchmark but are not package inputs. The small README projection is the+ -- sole published chart and is declared as package documentation above. bench/delaunay-compare/rust/Cargo.lock bench/delaunay-compare/rust/Cargo.toml bench/delaunay-compare/rust/src/lib.rs@@ -69,7 +73,7 @@ source-repository this type: git location: https://github.com/PaleRoses/moonlight.git- tag: moonlight-triangulation-1.3.0.0+ tag: moonlight-triangulation-1.3.0.2 subdir: moonlight-triangulation flag warnings-as-errors@@ -133,6 +137,7 @@ visibility: public hs-source-dirs: src-dcel exposed-modules:+ Moonlight.Triangulation.Alpha Moonlight.Triangulation.Types Moonlight.Triangulation.Math Moonlight.Triangulation.Exact@@ -190,8 +195,9 @@ -- An ExactCellSet already carries every handle and closure witness needed to -- interpret it as a generic finite cell complex. This component is the--- canonical owner of that interpretation and deliberately remains outside the--- portable geometry tower because its Homology dependency requires GHC 9.14.+-- canonical owner of that interpretation and of the alpha-birth lowering into+-- Homology. It deliberately remains outside the portable geometry tower+-- because its Homology dependency requires GHC 9.14. library cell-complex import: shared-properties visibility: public@@ -206,8 +212,10 @@ base >= 4.22 && < 5 , containers >= 0.6 && < 0.9 , moonlight-algebra:abstract >= 0.1 && < 0.2+ , moonlight-homology >= 0.1 && < 0.2 , moonlight-homology:cell-complex >= 0.1 && < 0.2 , moonlight-triangulation:dcel+ , vector >= 0.13 && < 0.14 library build import: shared-properties@@ -490,7 +498,10 @@ buildable: False build-depends: base >= 4.22 && < 5+ , containers >= 0.6 && < 0.9+ , moonlight-homology >= 0.1 && < 0.2 , moonlight-homology:cell-complex >= 0.1 && < 0.2+ , moonlight-triangulation:core , moonlight-triangulation:build , moonlight-triangulation:cell-complex , moonlight-triangulation:dcel@@ -620,6 +631,7 @@ common triangulation-benchmark-support-slice other-modules: BenchSupport+ build-depends: moonlight-triangulation:dcel common triangulation-build-benchmark-slice other-modules: Moonlight.Triangulation.BuildBench@@ -713,6 +725,31 @@ , moonlight-category:simplicial >= 1.1.0.0 && < 1.2 , moonlight-homology:cell-category >= 0.1 && < 0.2 , moonlight-homology:cell-complex >= 0.1 && < 0.2+ , moonlight-triangulation:build+ , moonlight-triangulation:cell-complex+ , moonlight-triangulation:dcel+ , vector >= 0.13 && < 0.14++-- One exact Delaunay alpha filtration versus independently restricting and+-- reducing the same chain complex at every critical radius. Construction of+-- the shared Delaunay geometry is outside all timed lanes.+benchmark moonlight-triangulation-alpha-bench+ import:+ shared-properties,+ triangulation-benchmark-support-slice+ type: exitcode-stdio-1.0+ main-is: Main.hs+ hs-source-dirs:+ bench/alpha+ bench/support+ ghc-options: -threaded -rtsopts "-with-rtsopts=-T"+ if impl(ghc < 9.14)+ buildable: False+ build-depends:+ base >= 4.22 && < 5+ , containers >= 0.6 && < 0.9+ , deepseq >= 1.5 && < 1.6+ , moonlight-homology >= 0.1 && < 0.2 , moonlight-triangulation:build , moonlight-triangulation:cell-complex , moonlight-triangulation:dcel
src-build/Moonlight/Triangulation/BulkLoad.hs view
@@ -48,7 +48,11 @@ , newOperationState , setCounter )-import Moonlight.Triangulation.Internal.CircleSweep (circleSweepInsert, radiallyOrderArena)+import Moonlight.Triangulation.Internal.CircleSweep+ ( circleSweepInsert+ , radiallyOrderArena+ , radiallyOrderGeometry+ ) import Moonlight.Triangulation.Internal.PointIndex ( MutablePointIndex , emptyPointIndex@@ -157,99 +161,117 @@ let !inputScale = recip (fromIntegral inputCount) !inputCenterX = inputSumX * inputScale !inputCenterY = inputSumY * inputScale- assignGeometryRadialDistances inputArena inputCenterX inputCenterY- orderedInputArena <- radiallyOrderArena inputArena- defaultedMutable <-+ (distanceArena, pointXMutable, pointYMutable, inputOrder) =+ MUV.unzip4 inputArena+ assignGeometryRadialDistances+ distanceArena+ pointXMutable+ pointYMutable+ inputCenterX+ inputCenterY+ distances <- U.unsafeFreeze distanceArena+ pointXs <- U.unsafeFreeze pointXMutable+ pointYs <- U.unsafeFreeze pointYMutable+ orderedInputArena <-+ radiallyOrderGeometry distances pointXs pointYs inputOrder+ defaultedOutcome <- newMutableDcelWithVertexDefault () unitElementDefaults (planarDcelCapacity inputCount)- let !mutable = defaultedVertexDcel defaultedMutable- operation <- newOperationState (halfEdgeCapacity mutable)- unique <- appendSortedGeometry defaultedMutable inputArena- orderedArena <-- if unique == inputCount- then pure orderedInputArena- else do- let !scale = recip (fromIntegral unique)- !uniqueArena = MUV.unsafeSlice 0 unique inputArena- (sumX, sumY) <- sumGeometryArena uniqueArena- recenterGeometryArena- (sumX * scale)- (sumY * scale)- uniqueArena- radiallyOrderArena uniqueArena- inserted <- circleSweepInsert mutable operation orderedArena- case inserted of+ case defaultedOutcome of Left failure -> pure (Left failure)- Right _ -> freezeTriangulation mutable+ Right defaultedMutable -> do+ let !mutable = defaultedVertexDcel defaultedMutable+ !dense = defaultedDenseVertexDcel defaultedMutable+ operation <- newOperationState (halfEdgeCapacity mutable)+ unique <- appendSortedGeometry defaultedMutable pointXs pointYs inputOrder+ orderedArena <-+ if unique == inputCount+ then pure orderedInputArena+ else do+ let !scale = recip (fromIntegral unique)+ !uniqueOrder = MUV.unsafeSlice 0 unique inputOrder+ (sumX, sumY) <- sumGeometryVertices dense uniqueOrder+ uniqueArena <-+ fillRadialArena+ mutable+ (sumX * scale)+ (sumY * scale)+ (MUV.unsafeRead uniqueOrder)+ unique+ radiallyOrderArena uniqueArena+ inserted <- circleSweepInsert mutable operation orderedArena+ case inserted of+ Left failure -> pure (Left failure)+ Right _ -> freezeTriangulation mutable where !inputCount = V.length coordinates assignGeometryRadialDistances :: forall s- . MUV.MVector s (Double, Double, Double, Word32)+ . MUV.MVector s Double+ -> MUV.MVector s Double+ -> MUV.MVector s Double -> Double -> Double -> ST s ()- assignGeometryRadialDistances arena centerX centerY =+ assignGeometryRadialDistances distances pointXs pointYs centerX centerY = MUV.imapM_- (\index (_, x, y, input) -> do+ (\index _ -> do+ x <- MUV.unsafeRead pointXs index+ y <- MUV.unsafeRead pointYs index let !deltaX = centerX - x !deltaY = centerY - y- MUV.unsafeWrite arena index (deltaX * deltaX + deltaY * deltaY, x, y, input)+ MUV.unsafeWrite distances index (deltaX * deltaX + deltaY * deltaY) )- arena+ distances appendSortedGeometry :: forall s . DefaultedVertexDcel s () () () ()- -> MUV.MVector s (Double, Double, Double, Word32)+ -> U.Vector Double+ -> U.Vector Double+ -> MUV.MVector s Word32 -> ST s Int- appendSortedGeometry defaultedMutable arena+ appendSortedGeometry defaultedMutable pointXs pointYs order | inputCount == 0 = pure 0 | otherwise = do- (distance, x, y, _) <- MUV.unsafeRead arena 0- first <- appendDefaultVertexCoordinates defaultedMutable x y- MUV.unsafeWrite arena 0 (distance, x, y, fromIntegral first)+ firstInput <- MUV.unsafeRead order 0+ let !x = U.unsafeIndex pointXs (fromIntegral firstInput)+ !y = U.unsafeIndex pointYs (fromIntegral firstInput)+ denseWriteFreshDefaultPoint defaultedMutable 0 x y+ MUV.unsafeWrite order 0 0 (_, _, unique) <- MUV.foldM'- (\(!previousX, !previousY, !uniqueCount) (nextDistance, nextX, nextY, _) ->- if nextX == previousX && nextY == previousY- then pure (previousX, previousY, uniqueCount)- else do- vertex <- appendDefaultVertexCoordinates defaultedMutable nextX nextY- MUV.unsafeWrite arena uniqueCount (nextDistance, nextX, nextY, fromIntegral vertex)- pure (nextX, nextY, uniqueCount + 1)+ (\(!previousX, !previousY, !uniqueCount) nextInput ->+ let !nextX = U.unsafeIndex pointXs (fromIntegral nextInput)+ !nextY = U.unsafeIndex pointYs (fromIntegral nextInput)+ in if nextX == previousX && nextY == previousY+ then pure (previousX, previousY, uniqueCount)+ else do+ denseWriteFreshDefaultPoint defaultedMutable uniqueCount nextX nextY+ MUV.unsafeWrite order uniqueCount (fromIntegral uniqueCount)+ pure (nextX, nextY, uniqueCount + 1) ) (x, y, 1)- (MUV.unsafeSlice 1 (inputCount - 1) arena)- pure unique+ (MUV.unsafeSlice 1 (inputCount - 1) order)+ unique <$ denseCommitFreshDefaultPoints defaultedMutable unique - sumGeometryArena+ sumGeometryVertices :: forall s- . MUV.MVector s (Double, Double, Double, Word32)+ . DenseMutableDcel s () () () ()+ -> MUV.MVector s Word32 -> ST s (Double, Double)- sumGeometryArena arena =+ sumGeometryVertices dense vertices = MUV.foldM'- (\(!sumX, !sumY) (_, x, y, _) -> pure (sumX + x, sumY + y))- (0, 0)- arena-- recenterGeometryArena- :: forall s- . Double- -> Double- -> MUV.MVector s (Double, Double, Double, Word32)- -> ST s ()- recenterGeometryArena centerX centerY arena =- MUV.imapM_- (\index (_, x, y, vertex) -> do- let !deltaX = centerX - x- !deltaY = centerY - y- MUV.unsafeWrite arena index (deltaX * deltaX + deltaY * deltaY, x, y, vertex)+ (\(!sumX, !sumY) raw -> do+ x <- denseReadPointX dense (fromIntegral raw)+ y <- denseReadPointY dense (fromIntegral raw)+ pure (sumX + x, sumY + y) )- arena+ (0, 0)+ vertices -- | Canonical construction from separate geometry and annotation sources. -- The coordinate vector remains the only geometry in ingress; payloads never
src-build/Moonlight/Triangulation/Internal/CircleSweep.hs view
@@ -7,15 +7,19 @@ module Moonlight.Triangulation.Internal.CircleSweep ( RadiallyOrderedArena , radiallyOrderArena+ , radiallyOrderGeometry , circleSweepInsert ) where import Control.Monad (forM_, when) import Control.Monad.ST (ST)-import Data.Bits (xor)+import Data.Bits ((.&.), shiftL, shiftR, xor) import qualified Data.Vector.Algorithms.Intro as Intro+import qualified Data.Vector.Algorithms.Radix as Radix+import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as MUV-import Data.Word (Word32)+import Data.Word (Word32, Word64)+import GHC.Float (castDoubleToWord64) import Moonlight.Triangulation.Handles.HandleDefs (DirectedEdgeId (..)) import Moonlight.Triangulation.Insertion (insertExistingVertex) import Moonlight.Triangulation.Internal.DcelOperations.Hull@@ -110,26 +114,155 @@ !(LegalizationArena s) {-# UNPACK #-} !SweepCellCursor --- | A packed sweep arena after its radial keys have descended to one total--- order. The constructor is private: the circle sweep consumes the proof and--- therefore never pays to establish the same ordering twice.+-- | Vertex ids after their radial keys have descended to one total order. The+-- constructor is private: the circle sweep consumes the proof and therefore+-- never pays to establish the same ordering twice. Coordinates remain owned+-- by the DCEL rather than being duplicated in the sweep input. newtype RadiallyOrderedArena s = RadiallyOrderedArena- (MUV.MVector s (Double, Double, Double, Word32))+ (MUV.MVector s Word32) radiallyOrderArena :: MUV.MVector s (Double, Double, Double, Word32) -> ST s (RadiallyOrderedArena s) radiallyOrderArena arena = do Intro.sort arena- pure (RadiallyOrderedArena arena)+ let (_, _, _, vertices) = MUV.unzip4 arena+ pure (RadiallyOrderedArena vertices) {-# INLINE radiallyOrderArena #-} +-- | Order the geometry-only ingress by moving only its existing vertex-id+-- plane. Stable radix descent first orders the high 36 bits of each binary64+-- distance. Every equal-prefix fibre is then normalized by the complete+-- historical @(distance, x, y, vertex)@ order, recovering the discarded low+-- 28 bits before the fibres glue into the authoritative radial order. The key+-- coordinates remain in their immutable planes, so the sole linear scratch+-- plane is four bytes per site rather than a second copy of every radial+-- record.+--+-- Coordinate admission bounds every squared distance to a finite,+-- non-negative binary64 value, whose unsigned word order is numeric order.+radiallyOrderGeometry+ :: U.Vector Double+ -> U.Vector Double+ -> U.Vector Double+ -> MUV.MVector s Word32+ -> ST s (RadiallyOrderedArena s)+radiallyOrderGeometry distances pointXs pointYs vertices = do+ Radix.sortBy+ radialPassCount+ radialBucketCount+ (radialDistanceBucket distances)+ vertices+ when (MUV.length vertices > 1) $ do+ first <- MUV.unsafeRead vertices 0+ let !firstPrefix = geometryDistancePrefix distances first+ (!finalFibreStart, _) <-+ MUV.ifoldM'+ (\(!fibreStart, !distancePrefix) relativeIndex vertex ->+ let !nextPrefix = geometryDistancePrefix distances vertex+ !index = relativeIndex + 1+ in if nextPrefix == distancePrefix+ then pure (fibreStart, distancePrefix)+ else do+ sortGeometryPrefixFibre+ distances+ pointXs+ pointYs+ vertices+ fibreStart+ index+ pure (index, nextPrefix)+ )+ (0, firstPrefix)+ (MUV.unsafeSlice 1 (MUV.length vertices - 1) vertices)+ sortGeometryPrefixFibre+ distances+ pointXs+ pointYs+ vertices+ finalFibreStart+ (MUV.length vertices)+ pure (RadiallyOrderedArena vertices)+{-# INLINE radiallyOrderGeometry #-}++radialPassCount :: Int+radialPassCount = 3++radialBucketCount :: Int+radialBucketCount = 1 `shiftL` radialPassWidth++radialPassWidth :: Int+radialPassWidth = 12++radialPrefixDiscardedBits :: Int+radialPrefixDiscardedBits = 64 - radialPassCount * radialPassWidth++radialBucketMask :: Word64+radialBucketMask = fromIntegral (radialBucketCount - 1)++radialDistanceBucket :: U.Vector Double -> Int -> Word32 -> Int+radialDistanceBucket distances pass vertex =+ fromIntegral+ ( (castDoubleToWord64 (geometryDistance distances vertex)+ `shiftR` (radialPrefixDiscardedBits + radialPassWidth * pass)+ )+ .&. radialBucketMask+ )+{-# INLINE radialDistanceBucket #-}++geometryDistancePrefix :: U.Vector Double -> Word32 -> Word64+geometryDistancePrefix distances vertex =+ castDoubleToWord64 (geometryDistance distances vertex)+ `shiftR` radialPrefixDiscardedBits+{-# INLINE geometryDistancePrefix #-}++geometryDistance :: U.Vector Double -> Word32 -> Double+geometryDistance distances vertex = U.unsafeIndex distances (fromIntegral vertex)+{-# INLINE geometryDistance #-}++compareGeometryVertex+ :: U.Vector Double+ -> U.Vector Double+ -> U.Vector Double+ -> Word32+ -> Word32+ -> Ordering+compareGeometryVertex distances pointXs pointYs left right =+ case compare (geometryDistance distances left) (geometryDistance distances right) of+ EQ ->+ case compare (coordinateAt pointXs left) (coordinateAt pointXs right) of+ EQ ->+ case compare (coordinateAt pointYs left) (coordinateAt pointYs right) of+ EQ -> compare left right+ yOrder -> yOrder+ xOrder -> xOrder+ distanceOrder -> distanceOrder+{-# INLINE compareGeometryVertex #-}++coordinateAt :: U.Vector Double -> Word32 -> Double+coordinateAt coordinates vertex = U.unsafeIndex coordinates (fromIntegral vertex)+{-# INLINE coordinateAt #-}++sortGeometryPrefixFibre+ :: U.Vector Double+ -> U.Vector Double+ -> U.Vector Double+ -> MUV.MVector s Word32+ -> Int+ -> Int+ -> ST s ()+sortGeometryPrefixFibre distances pointXs pointYs vertices start end =+ when (end - start > 1) $+ Intro.sortBy+ (compareGeometryVertex distances pointXs pointYs)+ (MUV.unsafeSlice start (end - start) vertices)+{-# INLINE sortGeometryPrefixFibre #-}+ noOuterEdge :: Word32 noOuterEdge = maxBound --- | Circle sweep over one mutable DCEL, consuming one packed radial arena of--- @(squaredDistance, x, y, vertex)@ records. The arena is sorted in place and--- then read directly — no decorated freeze, no undecoration pass. Insertions+-- | 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 -- full convexity, so construction does not repeatedly pay for global -- convexity that no intermediate observer can see.@@ -145,9 +278,10 @@ {-# INLINE circleSweepInsert #-} circleSweepInsertDense- :: DenseMutableDcel s vertex directed undirected face+ :: forall s vertex directed undirected face+ . DenseMutableDcel s vertex directed undirected face -> OperationState s- -> MUV.MVector s (Double, Double, Double, Word32)+ -> MUV.MVector s Word32 -> ST s (Either BuildError Int) circleSweepInsertDense dense operation arena | MUV.length arena == 0 = pure (Right 0)@@ -170,10 +304,21 @@ case builtHull of Left failure -> pure (Left failure) Right (!hull, !initialHullIndex) -> do- skipped <- MUV.new (ordered - seedCount) initialLegalizationArena <- legalizationArena operation let !initialCursor = initialSweepCellCursor reserved- inserted <- insertRemaining reserved initialCursor initialLegalizationArena hull initialHullIndex ordered seedCount skipped 0 0 0 0+ inserted <-+ insertRemaining+ reserved+ initialCursor+ initialLegalizationArena+ hull+ initialHullIndex+ ordered+ seedCount+ (0 :: Int)+ (0 :: Int)+ (0 :: Int)+ (0 :: Int) case inserted of Left failure -> pure (Left failure) Right (!skippedCount, !fastCount, !flips, !maxDepth, !sweepCursor, !sweepArena) -> do@@ -209,7 +354,7 @@ maxCounter operation CounterLegalizationMaxStack maxDepth addCounter operation CounterEdgeFlips terminalFlips maxCounter operation CounterLegalizationMaxStack terminalMaxDepth- insertedSkipped <- insertSkipped skipped skippedCount 0+ insertedSkipped <- insertSkipped skippedCount 0 pure (seedCount <$ insertedSkipped) where !mutable = denseMutableOwner dense@@ -217,7 +362,7 @@ insertSeed !index | index >= MUV.length arena = pure (Right index) | otherwise = do- (_, _, _, raw) <- MUV.unsafeRead arena index+ raw <- MUV.unsafeRead arena index result <- insertExistingVertex @'ProbeOff mutable operation (fromIntegral raw) case result of Left failure -> pure (Left failure)@@ -227,11 +372,26 @@ then pure (Right (index + 1)) else insertSeed (index + 1) - insertRemaining !reserved !cursor !candidateArena !hull !hullIndex !ordered !index !skipped !skippedCount !fastCount !flips !maxDepth+ insertRemaining+ :: ReservedSweepCells s vertex directed undirected face+ -> SweepCellCursor+ -> LegalizationArena s+ -> Hull s+ -> HullIndex s+ -> Int+ -> Int+ -> Int+ -> Int+ -> Int+ -> Int+ -> ST s (Either BuildError (Int, Int, Int, Int, SweepCellCursor, LegalizationArena s))+ insertRemaining !reserved !cursor !candidateArena !hull !hullIndex !ordered !index !skippedCount !fastCount !flips !maxDepth | index >= ordered = pure (Right (skippedCount, fastCount, flips, maxDepth, cursor, candidateArena)) | otherwise = do- (_, queryXWide, queryYWide, raw) <- MUV.unsafeRead arena index+ raw <- MUV.unsafeRead arena index let !vertex = fromIntegral raw+ queryXWide <- denseReadPointX dense vertex+ queryYWide <- denseReadPointY dense vertex let !queryAngle = pseudoAngle (hullCenterX hull) (hullCenterY hull) queryXWide queryYWide edge <- hullCandidate dense hull hullIndex queryAngle queryXWide queryYWide@@ -269,21 +429,23 @@ nextHullIndex ordered (index + 1)- skipped skippedCount (fastCount + 1) (flips + newFlips) (max maxDepth newMaxDepth) else do- MUV.unsafeWrite skipped skippedCount raw- insertRemaining reserved cursor candidateArena hull hullIndex ordered (index + 1) skipped (skippedCount + 1) fastCount flips maxDepth+ -- The radial order is consumed monotonically, so its processed+ -- prefix is dead. Compact deferred vertices into that prefix+ -- instead of allocating a second linear arena.+ MUV.unsafeWrite arena skippedCount raw+ insertRemaining reserved cursor candidateArena hull hullIndex ordered (index + 1) (skippedCount + 1) fastCount flips maxDepth - insertSkipped skipped !count = go+ insertSkipped !count = go where go !index | index >= count = pure (Right ()) | otherwise = do- vertex <- fromIntegral <$> MUV.unsafeRead skipped index+ vertex <- fromIntegral <$> MUV.unsafeRead arena index inserted <- insertExistingVertex @'ProbeOff mutable operation vertex case inserted of Left failure -> pure (Left failure)
src-cell-complex/Moonlight/Triangulation/CellComplex.hs view
@@ -1,27 +1,68 @@ {-# LANGUAGE EmptyDataDeriving #-} -- | The admitted Moonlight triangulation cell section as a generic--- 'CellComplex2D'. The 'ExactCellSet' remains the semantic owner: this module--- only supplies the incidence interpretation required by downstream topology.+-- 'CellComplex2D', an integral cellular chain complex, and an exact filtered+-- alpha complex. The 'ExactCellSet' remains the semantic owner: this module+-- supplies only the incidence and Homology interpretations required by+-- downstream topology. module Moonlight.Triangulation.CellComplex ( DCELComplex,- DCELError,+ DCELError (..), fromExactCellSet,+ finiteChainComplex,+ filteredAlphaComplex, ) where +import Data.Bifunctor (first) import Data.IntMap.Strict qualified as IntMap import Data.IntSet qualified as IntSet+import Data.List qualified as List+import Data.Map.Strict (Map)+import Data.Map.Strict qualified as Map+import Data.Vector qualified as Vector import Moonlight.Algebra.Pure.Orientation (Orientation (..))+import Moonlight.Homology.Boundary+ ( BoundaryIncidence+ , BoundaryIncidenceShapeError+ , BoundaryEntry+ , FiniteChainComplex+ , emptyBoundaryIncidence+ , emptyBoundaryIncidenceOf+ , mkBoundaryEntryFromInts+ , mkBoundaryIncidenceFromOrderedColumns+ , mkBoundaryIncidenceFromOrderedEntries+ , mkFiniteChainComplexChecked+ , targetIndex+ )+import Moonlight.Homology.Chain+ ( HomologicalDegree (..)+ , HomologyFailure+ )+import Moonlight.Homology.Persistence+ ( FilteredFiniteChainComplex+ , mkFilteredFiniteChainComplex+ ) import Moonlight.Homology.Pure.Topology.CellComplex- ( CellComplex2D (..),- CellTypes (..),- OrientedEdge (..),- ValidateComplex2D (..),+ ( CellComplex2D (..)+ , CellRef (..)+ , CellTypes (..)+ , OrientedEdge (..)+ , ValidateComplex2D (..) )+import Moonlight.Homology.Topology (BasisCellRef (..))+import Moonlight.Triangulation.Alpha+ ( AlphaBirth+ , AlphaFiltration+ , alphaEdgeBirth+ , alphaFaceBirth+ , alphaFiltrationCellSet+ , alphaVertexBirth+ ) import Moonlight.Triangulation.Dcel qualified as Dcel import Moonlight.Triangulation.Handles.HandleDefs- ( FaceId (..),+ ( DirectedEdgeId,+ FaceId (..), UndirectedEdgeId (..), VertexId (..), asUndirected,@@ -36,13 +77,367 @@ newtype DCELComplex = DCELComplex ExactCellSet -- | 'ExactCellSet' construction discharges every closure and handle--- obligation before this view exists, so there are no residual validation--- failures for the adapter to manufacture.+-- obligation before this view exists. Incidence materialization, Homology's+-- independent chain-law seal, and the exact-birth join retain typed+-- obstructions at their respective boundaries.+type DCELCellRef = CellRef VertexId UndirectedEdgeId FaceId+ data DCELError+ = DCELBoundaryCellMissing !DCELCellRef !DCELCellRef+ | DCELBoundaryIncidenceInvalid !BoundaryIncidenceShapeError+ | DCELChainComplexInvalid !HomologyFailure+ | DCELAlphaBirthMissing !DCELCellRef+ | DCELFilteredComplexInvalid !HomologyFailure deriving stock (Eq, Show) fromExactCellSet :: ExactCellSet -> DCELComplex fromExactCellSet = DCELComplex++-- | Canonical integral cellular chains in ascending resident-handle order.+-- Degree one uses target minus source; degree two uses the DCEL's oriented+-- face boundary. The Homology boundary seals the result only after checking+-- shape and @d . d = 0@; no unchecked chain constructor crosses the public+-- package boundary.+finiteChainComplex :: DCELComplex -> Either DCELError (FiniteChainComplex Int)+finiteChainComplex complexValue =+ let basis = dcelBasis complexValue+ in finiteChainComplexWithBasis complexValue basis++-- | Lower the exact alpha section into Homology without converting its birth+-- order through binary64. Persistence remains wholly owned by Homology.+filteredAlphaComplex+ :: AlphaFiltration+ -> Either DCELError (FilteredFiniteChainComplex AlphaBirth Int)+filteredAlphaComplex filtration =+ case alphaFiltrationCellSet filtration of+ ExactCellSet triangulation _ _ _ -> do+ finite <- residentAlphaFiniteChainComplex triangulation+ vertexBirthAssignments <-+ traverse+ ( residentAlphaBirthAssignment+ CellVertexRef+ residentVertexBasisRef+ (alphaVertexBirth filtration)+ )+ (vertexHandlesOf triangulation)+ edgeBirthAssignments <-+ traverse+ ( residentAlphaBirthAssignment+ CellEdgeRef+ residentEdgeBasisRef+ (alphaEdgeBirth filtration)+ )+ (undirectedEdgesOf triangulation)+ faceBirthAssignments <-+ traverse+ ( residentAlphaBirthAssignment+ CellFaceRef+ residentFaceBasisRef+ (alphaFaceBirth filtration)+ )+ (innerFacesOf triangulation)+ first DCELFilteredComplexInvalid+ ( mkFilteredFiniteChainComplex+ finite+ (vertexBirthAssignments <> edgeBirthAssignments <> faceBirthAssignments)+ )++-- | The opaque 'AlphaFiltration' constructor admits the entire resident DCEL,+-- whose handle ranges are dense. This local section therefore lowers those+-- handles directly to basis indices while retaining Homology's independent+-- shape and nilpotence seal. Sparse 'ExactCellSet' values continue through the+-- generic map-indexed 'finiteChainComplex' path.+residentAlphaFiniteChainComplex+ :: Triangulation mode vertex directed undirected face+ -> Either DCELError (FiniteChainComplex Int)+residentAlphaFiniteChainComplex triangulation = do+ let degreeOneColumns =+ Vector.generate+ (Dcel.numUndirectedEdges triangulation)+ (residentEdgeBoundaryEntries triangulation . UndirectedEdgeId . fromIntegral)+ degreeTwoColumns =+ Vector.imap+ residentFaceBoundaryEntries+ (Dcel.innerFaceDirectedEdgeTriples triangulation)+ degreeOneBoundary <-+ first DCELBoundaryIncidenceInvalid+ ( mkBoundaryIncidenceFromOrderedColumns+ (fromIntegral (Dcel.numUndirectedEdges triangulation))+ (fromIntegral (Dcel.numVertices triangulation))+ degreeOneColumns+ )+ degreeTwoBoundary <-+ first DCELBoundaryIncidenceInvalid+ ( mkBoundaryIncidenceFromOrderedColumns+ (fromIntegral (Dcel.numInnerFaces triangulation))+ (fromIntegral (Dcel.numUndirectedEdges triangulation))+ degreeTwoColumns+ )+ let degreeZeroBoundary =+ emptyBoundaryIncidenceOf+ (fromIntegral (Dcel.numVertices triangulation))+ 0+ boundaryAt (HomologicalDegree degreeValue) =+ case degreeValue of+ 0 -> degreeZeroBoundary+ 1 -> degreeOneBoundary+ 2 -> degreeTwoBoundary+ _ -> emptyBoundaryIncidence+ first DCELChainComplexInvalid+ (mkFiniteChainComplexChecked (HomologicalDegree 2) boundaryAt)++residentEdgeBoundaryEntries+ :: Triangulation mode vertex directed undirected face+ -> UndirectedEdgeId+ -> [BoundaryEntry Int]+residentEdgeBoundaryEntries triangulation edgeValue =+ let (sourceVertex, targetVertex) = Dcel.undirectedEndpoints triangulation edgeValue+ sourceEntry =+ mkBoundaryEntryFromInts+ (basisIndexOfEdge edgeValue)+ (basisIndexOfVertex sourceVertex)+ (-1)+ targetEntry =+ mkBoundaryEntryFromInts+ (basisIndexOfEdge edgeValue)+ (basisIndexOfVertex targetVertex)+ 1+ in if targetIndex sourceEntry <= targetIndex targetEntry+ then [sourceEntry, targetEntry]+ else [targetEntry, sourceEntry]++residentFaceBoundaryEntries+ :: Int+ -> (DirectedEdgeId, DirectedEdgeId, DirectedEdgeId)+ -> [BoundaryEntry Int]+residentFaceBoundaryEntries faceIndexValue (firstEdge, secondEdge, thirdEdge) =+ sortThreeBoundaryEntries+ (residentFaceBoundaryEntry faceIndexValue firstEdge)+ (residentFaceBoundaryEntry faceIndexValue secondEdge)+ (residentFaceBoundaryEntry faceIndexValue thirdEdge)++sortThreeBoundaryEntries+ :: BoundaryEntry Int+ -> BoundaryEntry Int+ -> BoundaryEntry Int+ -> [BoundaryEntry Int]+sortThreeBoundaryEntries firstEntry secondEntry thirdEntry =+ let (firstLow, firstHigh) = orderedBoundaryPair firstEntry secondEntry+ (secondLow, finalHigh) = orderedBoundaryPair firstHigh thirdEntry+ (finalLow, finalMiddle) = orderedBoundaryPair firstLow secondLow+ in [finalLow, finalMiddle, finalHigh]++orderedBoundaryPair+ :: BoundaryEntry Int+ -> BoundaryEntry Int+ -> (BoundaryEntry Int, BoundaryEntry Int)+orderedBoundaryPair firstEntry secondEntry =+ if targetIndex firstEntry <= targetIndex secondEntry+ then (firstEntry, secondEntry)+ else (secondEntry, firstEntry)++residentFaceBoundaryEntry :: Int -> DirectedEdgeId -> BoundaryEntry Int+residentFaceBoundaryEntry faceIndexValue directedEdge =+ mkBoundaryEntryFromInts+ faceIndexValue+ (basisIndexOfEdge (asUndirected directedEdge))+ (if isNormalized directedEdge then 1 else -1)++residentAlphaBirthAssignment+ :: (cell -> DCELCellRef)+ -> (cell -> BasisCellRef)+ -> (cell -> Maybe AlphaBirth)+ -> cell+ -> Either DCELError (BasisCellRef, AlphaBirth)+residentAlphaBirthAssignment cellReference basisReference birthAt cell =+ maybe+ (Left (DCELAlphaBirthMissing (cellReference cell)))+ (Right . (,) (basisReference cell))+ (birthAt cell)++residentVertexBasisRef :: VertexId -> BasisCellRef+residentVertexBasisRef vertexValue =+ BasisCellRef (HomologicalDegree 0) (basisIndexOfVertex vertexValue)++residentEdgeBasisRef :: UndirectedEdgeId -> BasisCellRef+residentEdgeBasisRef edgeValue =+ BasisCellRef (HomologicalDegree 1) (basisIndexOfEdge edgeValue)++residentFaceBasisRef :: FaceId -> BasisCellRef+residentFaceBasisRef faceValue =+ BasisCellRef (HomologicalDegree 2) (basisIndexOfFace faceValue)++basisIndexOfVertex :: VertexId -> Int+basisIndexOfVertex (VertexId rawVertex) = fromIntegral rawVertex++basisIndexOfEdge :: UndirectedEdgeId -> Int+basisIndexOfEdge (UndirectedEdgeId rawEdge) = fromIntegral rawEdge++basisIndexOfFace :: FaceId -> Int+basisIndexOfFace (FaceId rawFace) = fromIntegral rawFace - 1++undirectedEdgesOf+ :: Triangulation mode vertex directed undirected face+ -> [UndirectedEdgeId]+undirectedEdgesOf triangulation =+ fmap (UndirectedEdgeId . fromIntegral) [0 .. Dcel.numUndirectedEdges triangulation - 1]++vertexHandlesOf+ :: Triangulation mode vertex directed undirected face+ -> [VertexId]+vertexHandlesOf triangulation =+ fmap (VertexId . fromIntegral) [0 .. Dcel.numVertices triangulation - 1]++innerFacesOf+ :: Triangulation mode vertex directed undirected face+ -> [FaceId]+innerFacesOf triangulation =+ fmap (FaceId . fromIntegral) [1 .. Dcel.numFaces triangulation - 1]++data DCELBasis = DCELBasis+ { dcelVertexBasis :: !(Map VertexId BasisCellRef)+ , dcelEdgeBasis :: !(Map UndirectedEdgeId BasisCellRef)+ , dcelFaceBasis :: !(Map FaceId BasisCellRef)+ }++dcelBasis :: DCELComplex -> DCELBasis+dcelBasis complexValue =+ DCELBasis+ { dcelVertexBasis = basisMap 0 (vertices complexValue)+ , dcelEdgeBasis = basisMap 1 (edges complexValue)+ , dcelFaceBasis = basisMap 2 (faces complexValue)+ }++basisMap :: Ord cell => Int -> [cell] -> Map cell BasisCellRef+basisMap degreeValue cells =+ Map.fromAscList+ ( zipWith+ (\indexValue cell -> (cell, BasisCellRef (HomologicalDegree degreeValue) indexValue))+ [0 ..]+ cells+ )++finiteChainComplexWithBasis+ :: DCELComplex+ -> DCELBasis+ -> Either DCELError (FiniteChainComplex Int)+finiteChainComplexWithBasis complexValue basis = do+ degreeOneEntries <-+ concat+ <$> traverse+ (edgeBoundaryEntries complexValue basis)+ (Map.toAscList (dcelEdgeBasis basis))+ degreeTwoEntries <-+ concat+ <$> traverse+ (faceBoundaryEntries complexValue basis)+ (Map.toAscList (dcelFaceBasis basis))+ degreeOneBoundary <-+ first DCELBoundaryIncidenceInvalid+ ( mkBoundaryIncidenceFromOrderedEntries+ (fromIntegral (Map.size (dcelEdgeBasis basis)))+ (fromIntegral (Map.size (dcelVertexBasis basis)))+ degreeOneEntries+ )+ degreeTwoBoundary <-+ first DCELBoundaryIncidenceInvalid+ ( mkBoundaryIncidenceFromOrderedEntries+ (fromIntegral (Map.size (dcelFaceBasis basis)))+ (fromIntegral (Map.size (dcelEdgeBasis basis)))+ degreeTwoEntries+ )+ let degreeZeroBoundary =+ emptyBoundaryIncidenceOf+ (fromIntegral (Map.size (dcelVertexBasis basis)))+ 0+ boundaryAt :: HomologicalDegree -> BoundaryIncidence Int+ boundaryAt (HomologicalDegree degreeValue) =+ case degreeValue of+ 0 -> degreeZeroBoundary+ 1 -> degreeOneBoundary+ 2 -> degreeTwoBoundary+ _ -> emptyBoundaryIncidence+ first DCELChainComplexInvalid+ (mkFiniteChainComplexChecked (HomologicalDegree 2) boundaryAt)++edgeBoundaryEntries+ :: DCELComplex+ -> DCELBasis+ -> (UndirectedEdgeId, BasisCellRef)+ -> Either DCELError [BoundaryEntry Int]+edgeBoundaryEntries complexValue basis (edgeValue, edgeBasisRef) = do+ let (sourceVertex, targetVertex) = edgeBoundary complexValue edgeValue+ sourceCell = CellEdgeRef edgeValue+ sourceBasisRef <-+ requireBoundaryCell+ sourceCell+ (CellVertexRef sourceVertex)+ sourceVertex+ (dcelVertexBasis basis)+ targetBasisRef <-+ requireBoundaryCell+ sourceCell+ (CellVertexRef targetVertex)+ targetVertex+ (dcelVertexBasis basis)+ pure+ ( List.sortOn+ targetIndex+ [ mkBoundaryEntryFromInts+ (cellIndex edgeBasisRef)+ (cellIndex sourceBasisRef)+ (-1)+ , mkBoundaryEntryFromInts+ (cellIndex edgeBasisRef)+ (cellIndex targetBasisRef)+ 1+ ]+ )++faceBoundaryEntries+ :: DCELComplex+ -> DCELBasis+ -> (FaceId, BasisCellRef)+ -> Either DCELError [BoundaryEntry Int]+faceBoundaryEntries complexValue basis (faceValue, faceBasisRef) =+ List.sortOn targetIndex+ <$> traverse+ boundaryEntry+ (faceBoundary complexValue faceValue)+ where+ boundaryEntry orientedBoundary = do+ let edgeValue = orientedEdge orientedBoundary+ edgeBasisRef <-+ requireBoundaryCell+ (CellFaceRef faceValue)+ (CellEdgeRef edgeValue)+ edgeValue+ (dcelEdgeBasis basis)+ pure+ ( mkBoundaryEntryFromInts+ (cellIndex faceBasisRef)+ (cellIndex edgeBasisRef)+ (orientationCoefficient (edgeOrientation orientedBoundary))+ )++orientationCoefficient :: Orientation -> Int+orientationCoefficient orientation =+ case orientation of+ Positive -> 1+ Negative -> -1++requireBoundaryCell+ :: Ord cell+ => DCELCellRef+ -> DCELCellRef+ -> cell+ -> Map cell BasisCellRef+ -> Either DCELError BasisCellRef+requireBoundaryCell sourceCell targetCell cell basis =+ maybe+ (Left (DCELBoundaryCellMissing sourceCell targetCell))+ Right+ (Map.lookup cell basis) instance CellTypes DCELComplex where type Vertex DCELComplex = VertexId
src-core/Moonlight/Triangulation/Internal/Dyadic.hs view
@@ -7,6 +7,8 @@ ( exactOrientDet , exactOrientSignDouble , exactInCircleDet+ , exactQuarterSquaredDistanceRational+ , exactCircumradiusSquaredRational , exactCircumradiusSquaredWithin , exactBarycentricDeterminants , exactDiametralDot@@ -23,6 +25,14 @@ ) import qualified Data.List as List import Data.Word (Word64)+import GHC.Exts (Int (I#), (+#))+import GHC.Integer.Logarithms (integerLog2#)+import Moonlight.Triangulation.Internal.ExactRational+ ( ExactArithmeticError+ , ExactRational+ , exactRationalFromDyadic+ , exactRationalFromDyadicRatio+ ) #if WORD_SIZE_IN_BITS == 64 import GHC.Exts ( Double (D#)@@ -45,7 +55,6 @@ , word2Int# , word64ToWord# , (*#)- , (+#) , (-#) , (<=#) , (==#)@@ -74,6 +83,55 @@ | otherwise = mantissa `shiftL` (sourcePower - power) {-# INLINE alignDecoded #-} +decodedExponentFloor :: Decoded -> Int+decodedExponentFloor (mantissa, power)+ | mantissa == 0 = 0+ | otherwise = min 0 power+{-# INLINE decodedExponentFloor #-}++commonExponent4 :: Decoded -> Decoded -> Decoded -> Decoded -> Int+commonExponent4 firstValue secondValue thirdValue fourthValue =+ min+ (min (decodedExponentFloor firstValue) (decodedExponentFloor secondValue))+ (min (decodedExponentFloor thirdValue) (decodedExponentFloor fourthValue))+{-# INLINE commonExponent4 #-}++commonExponent6+ :: Decoded+ -> Decoded+ -> Decoded+ -> Decoded+ -> Decoded+ -> Decoded+ -> Int+commonExponent6 firstValue secondValue thirdValue fourthValue fifthValue sixthValue =+ min+ ( commonExponent4+ firstValue+ secondValue+ thirdValue+ fourthValue+ )+ (min (decodedExponentFloor fifthValue) (decodedExponentFloor sixthValue))+{-# INLINE commonExponent6 #-}++aligned4+ :: Double -> Double -> Double -> Double+ -> (Int, Integer, Integer, Integer, Integer)+aligned4 a b c d =+ let !da = decodeFloat a+ !db = decodeFloat b+ !dc = decodeFloat c+ !dd = decodeFloat d+ !power = commonExponent4 da db dc dd+ in ( power+ , alignDecoded power da+ , alignDecoded power db+ , alignDecoded power dc+ , alignDecoded power dd+ )+{-# INLINE aligned4 #-}+ aligned6 :: Double -> Double -> Double -> Double -> Double -> Double -> (Int, Integer, Integer, Integer, Integer, Integer, Integer)@@ -84,7 +142,7 @@ !dd = decodeFloat d !de = decodeFloat e !df = decodeFloat f- !power = commonExponent [da, db, dc, dd, de, df]+ !power = commonExponent6 da db dc dd de df in ( power , alignDecoded power da , alignDecoded power db@@ -127,6 +185,48 @@ !bcy = iby - icy in acx * bcy - acy * bcx +-- | One quarter of the exact squared distance between two finite binary64+-- points. Alpha's Gabriel-edge birth needs this value directly, so the dyadic+-- exponent absorbs the quarter before the single rational normalization.+exactQuarterSquaredDistanceRational+ :: Double -> Double -> Double -> Double -> ExactRational+exactQuarterSquaredDistanceRational ax ay bx by =+ let (!coordinatePower, !iax, !iay, !ibx, !iby) = aligned4 ax ay bx by+ !deltaX = ibx - iax+ !deltaY = iby - iay+ !squaredDistance = deltaX * deltaX + deltaY * deltaY+ in exactRationalFromDyadic squaredDistance (2 * coordinatePower - 2)+{-# INLINE exactQuarterSquaredDistanceRational #-}++-- | Exact squared circumradius of three finite binary64 points. The local+-- dyadic section computes the integer numerator and denominator without+-- constructing intermediate ratios; descent normalizes the authoritative+-- result exactly once. A degenerate triangle is the typed zero-divisor+-- obstruction from 'ExactRational'.+exactCircumradiusSquaredRational+ :: Double -> Double -> Double -> Double -> Double -> Double+ -> Either ExactArithmeticError ExactRational+exactCircumradiusSquaredRational ax ay bx by cx cy =+ let (!coordinatePower, !iax, !iay, !ibx, !iby, !icx, !icy) =+ aligned6 ax ay bx by cx cy+ !abx = ibx - iax+ !aby = iby - iay+ !acx = icx - iax+ !acy = icy - iay+ !bcx = icx - ibx+ !bcy = icy - iby+ !abSquared = abx * abx + aby * aby+ !acSquared = acx * acx + acy * acy+ !bcSquared = bcx * bcx + bcy * bcy+ !determinant = abx * acy - aby * acx+ !radiusNumerator = abSquared * acSquared * bcSquared+ !radiusDenominator = 4 * determinant * determinant+ in exactRationalFromDyadicRatio+ radiusNumerator+ radiusDenominator+ (2 * coordinatePower)+{-# INLINE exactCircumradiusSquaredRational #-}+ -- | Exact closed comparison of squared circumradius with a finite, -- non-negative binary64 threshold. No constructed circumcenter participates. exactCircumradiusSquaredWithin@@ -346,16 +446,10 @@ in if sameSign then magnitude else negate magnitude integerBitLength :: Integer -> Int-integerBitLength = go 0- where- go !bits value- | value <= 0xffffffff = bits + wordBitLength value- | otherwise = go (bits + 32) (value `shiftR` 32)-- wordBitLength = count 0- count :: Int -> Integer -> Int- count !bits 0 = bits- count !bits value = count (bits + 1) (value `shiftR` 1)+integerBitLength value+ | value <= 0 = 0+ | otherwise = I# (integerLog2# value +# 1#)+{-# INLINE integerBitLength #-} -- --------------------------------------------------------------------------- -- Fixed-precision exact orient sign for Double.
src-core/Moonlight/Triangulation/Internal/ExactRational.hs view
@@ -2,6 +2,7 @@ {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE DerivingStrategies #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE MagicHash #-} -- | Normalized exact rational arithmetic without geometric dependencies. module Moonlight.Triangulation.Internal.ExactRational@@ -11,6 +12,7 @@ , exactRationalFromDouble , exactRationalFromFiniteDouble , exactRationalFromDyadic+ , exactRationalFromDyadicRatio , exactRationalNumerator , exactRationalDenominator , exactRationalIsZero@@ -19,10 +21,13 @@ ) where import Control.DeepSeq (NFData)-import Data.Bits (shiftL)+import Data.Bits ((.&.), shiftL, shiftR) import Data.Ratio (Ratio, (%)) import qualified Data.Ratio as Ratio import GHC.Generics (Generic)+import GHC.Exts (Int (I#))+import GHC.Integer.Logarithms (integerLog2#)+import GHC.Real (Ratio ((:%))) -- | A checked wrapper around a reduced ratio with a strictly positive -- denominator. 'Ratio' owns normalization, including the unique zero@@ -65,9 +70,48 @@ -- | Construct @numerator * 2^power@ without a partial denominator path. exactRationalFromDyadic :: Integer -> Int -> ExactRational exactRationalFromDyadic numerator power- | power >= 0 = fromInteger (numerator `shiftL` power)- | otherwise = ExactRational (numerator % (1 `shiftL` negate power))+ | numerator == 0 = ExactRational (0 :% 1)+ | power >= 0 = ExactRational ((numerator `shiftL` power) :% 1)+ | otherwise =+ let denominatorPower = negate power+ removableFactor =+ min denominatorPower (integerTrailingZeroBits numerator)+ reducedNumerator = numerator `shiftR` removableFactor+ reducedDenominator =+ 1 `shiftL` (denominatorPower - removableFactor)+ in ExactRational (reducedNumerator :% reducedDenominator) {-# INLINE exactRationalFromDyadic #-}++-- | The denominator of a dyadic rational is a power of two, so its complete+-- normalization requires only the numerator's least set bit. Calling generic+-- rational GCD here merely rediscovers that closed arithmetic fact.+integerTrailingZeroBits :: Integer -> Int+integerTrailingZeroBits value =+ let magnitude = abs value+ in I# (integerLog2# (magnitude .&. negate magnitude))+{-# INLINE integerTrailingZeroBits #-}++-- | Construct @(numerator / denominator) * 2^power@ with one final rational+-- normalization. This is the publication boundary for exact dyadic kernels:+-- their integer arithmetic must not pay a greatest-common-divisor reduction+-- after every intermediate operation.+exactRationalFromDyadicRatio+ :: Integer+ -> Integer+ -> Int+ -> Either ExactArithmeticError ExactRational+exactRationalFromDyadicRatio _ 0 _ = Left ExactZeroDivisor+exactRationalFromDyadicRatio 0 _ _ = Right (ExactRational (0 :% 1))+exactRationalFromDyadicRatio numerator denominator power =+ let numeratorFactor = integerTrailingZeroBits numerator+ denominatorFactor = integerTrailingZeroBits denominator+ oddNumerator = numerator `shiftR` numeratorFactor+ oddDenominator = denominator `shiftR` denominatorFactor+ residualPower = power + numeratorFactor - denominatorFactor+ in if residualPower >= 0+ then exactRational (oddNumerator `shiftL` residualPower) oddDenominator+ else exactRational oddNumerator (oddDenominator `shiftL` negate residualPower)+{-# INLINE exactRationalFromDyadicRatio #-} -- | Read the reduced numerator. exactRationalNumerator :: ExactRational -> Integer
+ src-dcel/Moonlight/Triangulation/Alpha.hs view
@@ -0,0 +1,424 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveAnyClass #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DerivingStrategies #-}++-- | Exact alpha births over one resident Delaunay triangulation. Geometry+-- remains owned by the DCEL and 'ExactCellSet'; this module supplies the+-- handle-indexed filtration section used by topology consumers.+module Moonlight.Triangulation.Alpha+ ( AlphaBirth+ , alphaBirthNumerator+ , alphaBirthDenominator+ , alphaBirthToDouble+ , alphaBirthFromRadiusSquared+ , AlphaFiltration+ , AlphaFiltrationError (..)+ , alphaFiltration+ , alphaFiltrationCellSet+ , alphaVertexBirth+ , alphaEdgeBirth+ , alphaFaceBirth+ , alphaFiltrationCriticalBirths+ , alphaComplexAtBirth+ , alphaComplexAtRadius+ , alphaShapeContainsFace+ ) where++import Control.DeepSeq (NFData)+import Data.Bifunctor (first)+import Data.Bits (shiftL)+import Data.List.NonEmpty qualified as NonEmpty+import Data.Set qualified as Set+import Data.Vector (Vector)+import Data.Vector qualified as Vector+import Data.Word (Word64)+import GHC.Generics (Generic)+import Moonlight.Triangulation.CellSet+ ( CellSelectionError+ , exactCellSet+ )+import Moonlight.Triangulation.Dcel qualified as Dcel+import Moonlight.Triangulation.Handles.HandleDefs+ ( FaceId (..)+ , UndirectedEdgeId (..)+ , VertexId (..)+ , directedPair+ )+import Moonlight.Triangulation.Internal.CellSet+ ( ExactCellSet (..)+ , residentExactCellSet+ )+import Moonlight.Triangulation.Internal.Dyadic+ ( exactCircumradiusSquaredRational+ , exactQuarterSquaredDistanceRational+ , integerBitLength+ , integerRatioToDouble+ )+import Moonlight.Triangulation.Internal.ExactRational+ ( ExactArithmeticError (..)+ , ExactRational+ , exactRationalDenominator+ , exactRationalFromFiniteDouble+ , exactRationalNumerator+ )+import Moonlight.Triangulation.Internal.Representation+ ( DelaunayTriangulation+ , Triangulation+ )+import Moonlight.Triangulation.Math (inDiametralCircle)+import Moonlight.Triangulation.Scalar (circumradiusSquaredWithinCoordinates)+import Moonlight.Triangulation.Types+ ( ConstraintMode (Unconstrained)+ , Point (..)+ , RadiusSquared+ , radiusSquaredValue+ )++-- | One exact squared-radius birth. The constructor is withheld so negative+-- exact rationals cannot masquerade as geometric radii.+data AlphaBirth = AlphaBirth !AlphaOrderKey !ExactRational+ deriving stock (Generic)+ deriving anyclass (NFData)++data AlphaOrderKey+ = NonPositiveAlphaOrder+ | PositiveAlphaOrder !Int !Word64+ deriving stock (Eq, Ord, Generic)+ deriving anyclass (NFData)++instance Eq AlphaBirth where+ AlphaBirth _ leftValue == AlphaBirth _ rightValue = leftValue == rightValue++instance Ord AlphaBirth where+ compare (AlphaBirth leftKey leftValue) (AlphaBirth rightKey rightValue) =+ case compare leftKey rightKey of+ EQ -> compare leftValue rightValue+ distinctKeyOrder -> distinctKeyOrder++instance Show AlphaBirth where+ showsPrec precedence (AlphaBirth _ value) =+ showParen+ (precedence > 10)+ (showString "AlphaBirth " . showsPrec 11 value)++-- | Numerator of the reduced exact squared-radius birth.+alphaBirthNumerator :: AlphaBirth -> Integer+alphaBirthNumerator (AlphaBirth _ value) = exactRationalNumerator value++-- | Positive denominator of the reduced exact squared-radius birth.+alphaBirthDenominator :: AlphaBirth -> Integer+alphaBirthDenominator (AlphaBirth _ value) = exactRationalDenominator value++-- | Binary64 display projection. Equality and filtration order must use the+-- 'AlphaBirth' itself, not this rounded view.+alphaBirthToDouble :: AlphaBirth -> Double+alphaBirthToDouble birth =+ integerRatioToDouble+ (alphaBirthNumerator birth)+ (alphaBirthDenominator birth)++-- | Embed an admitted binary64 threshold into the exact alpha order without+-- loss. Every finite binary64 value is a dyadic rational.+alphaBirthFromRadiusSquared :: RadiusSquared -> AlphaBirth+alphaBirthFromRadiusSquared =+ alphaBirthFromExact . exactRationalFromFiniteDouble . radiusSquaredValue++-- | The exact handle-indexed filtration of one Delaunay DCEL. The cell set is+-- the only geometry inventory; the maps are derived birth sections over those+-- resident handles.+data AlphaFiltration = AlphaFiltration+ { alphaFiltrationCells :: !ExactCellSet+ , alphaVertexBirths :: !(Vector AlphaBirth)+ , alphaEdgeBirths :: !(Vector AlphaBirth)+ , alphaFaceBirths :: !(Vector AlphaBirth)+ }++data FaceAlphaSection = FaceAlphaSection+ { faceAlphaBirth :: !AlphaBirth+ , faceAlphaVertices :: !(VertexId, VertexId, VertexId)+ }++-- | Typed obstruction to deriving the exact birth section from a resident+-- Delaunay triangulation. Every incidence failure retains the affected handle.+data AlphaFiltrationError+ = AlphaFiltrationCellSelection !CellSelectionError+ | AlphaFaceVerticesUnavailable !FaceId+ | AlphaFaceDegenerate !FaceId+ | AlphaCircumradiusDivisionFailed !FaceId !ExactArithmeticError+ | AlphaFaceEdgeIncidenceInvalid !FaceId !UndirectedEdgeId ![VertexId]+ | AlphaIncidentFaceBirthMissing !UndirectedEdgeId !FaceId+ | AlphaEdgeBirthWitnessMissing !UndirectedEdgeId+ deriving stock (Eq, Show, Generic)+ deriving anyclass (NFData)++-- | Construct every exact alpha birth once. Faces descend first because a+-- non-Gabriel edge is born at the least incident-face circumradius.+alphaFiltration+ :: DelaunayTriangulation vertex+ -> Either AlphaFiltrationError AlphaFiltration+alphaFiltration triangulation = do+ cellSet <-+ first AlphaFiltrationCellSelection+ (residentExactCellSet triangulation)+ let vertexBirths = Vector.replicate (Dcel.numVertices triangulation) zeroAlphaBirth+ faceSections <-+ Vector.imapM+ (\faceIndex -> exactFaceSection triangulation (FaceId (fromIntegral (faceIndex + 1))))+ (Dcel.innerFaceVertexTriples triangulation)+ edgeBirths <-+ Vector.generateM+ (Dcel.numUndirectedEdges triangulation)+ (exactEdgeBirth faceSections triangulation . UndirectedEdgeId . fromIntegral)+ let faceBirths = fmap faceAlphaBirth faceSections+ pure+ AlphaFiltration+ { alphaFiltrationCells = cellSet+ , alphaVertexBirths = vertexBirths+ , alphaEdgeBirths = edgeBirths+ , alphaFaceBirths = faceBirths+ }++-- | The authoritative resident cells underlying this filtration. This is a+-- read-only projection rather than an exported record field, so callers cannot+-- splice births from one triangulation onto the cells of another.+alphaFiltrationCellSet :: AlphaFiltration -> ExactCellSet+alphaFiltrationCellSet = alphaFiltrationCells++-- | 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++-- | 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++-- | Exact birth of a resident bounded face, when the handle belongs to this+-- filtration.+alphaFaceBirth :: AlphaFiltration -> FaceId -> Maybe AlphaBirth+alphaFaceBirth filtration =+ faceBirthAt (alphaFaceBirths filtration)++-- | Distinct critical births in exact ascending order.+alphaFiltrationCriticalBirths :: AlphaFiltration -> [AlphaBirth]+alphaFiltrationCriticalBirths filtration =+ Set.toAscList+ ( Set.fromList+ ( Vector.toList (alphaVertexBirths filtration)+ <> Vector.toList (alphaEdgeBirths filtration)+ <> Vector.toList (alphaFaceBirths filtration)+ )+ )++-- | Select the closed alpha subcomplex at one exact birth. The resident DCEL+-- is retained; only its downward-closed cell section changes.+alphaComplexAtBirth+ :: AlphaBirth+ -> AlphaFiltration+ -> Either CellSelectionError ExactCellSet+alphaComplexAtBirth threshold filtration =+ case alphaFiltrationCells filtration of+ ExactCellSet triangulation _ _ _ ->+ exactCellSet+ triangulation+ (selectedHandles (VertexId . fromIntegral) (alphaVertexBirths filtration))+ (selectedHandles (UndirectedEdgeId . fromIntegral) (alphaEdgeBirths filtration))+ (selectedHandles (FaceId . fromIntegral . (+ 1)) (alphaFaceBirths filtration))+ where+ selectedHandles :: (Int -> handle) -> Vector AlphaBirth -> [handle]+ selectedHandles handleAt =+ Vector.ifoldr+ (\indexValue birthValue selected -> if birthValue <= threshold then handleAt indexValue : selected else selected)+ []++alphaComplexAtRadius+ :: RadiusSquared+ -> AlphaFiltration+ -> Either CellSelectionError ExactCellSet+alphaComplexAtRadius radius =+ alphaComplexAtBirth (alphaBirthFromRadiusSquared radius)++-- | Membership of a bounded face in the closed alpha shape. Exact dyadic+-- comparison makes equality independent of circumcenter rounding. This+-- specialized observation does not construct the full filtration.+alphaShapeContainsFace+ :: RadiusSquared+ -> Triangulation 'Unconstrained vertex directed undirected face+ -> FaceId+ -> Bool+alphaShapeContainsFace radius triangulation =+ maybe False withinRadius . Dcel.innerFaceVertices triangulation+ where+ threshold = radiusSquaredValue radius+ withinRadius (firstVertex, secondVertex, thirdVertex) =+ let Point ax ay = Dcel.vertexPoint triangulation firstVertex+ Point bx by = Dcel.vertexPoint triangulation secondVertex+ Point cx cy = Dcel.vertexPoint triangulation thirdVertex+ in circumradiusSquaredWithinCoordinates threshold ax ay bx by cx cy++zeroAlphaBirth :: AlphaBirth+zeroAlphaBirth = alphaBirthFromExact 0++exactFaceSection+ :: DelaunayTriangulation vertex+ -> FaceId+ -> (VertexId, VertexId, VertexId)+ -> Either AlphaFiltrationError FaceAlphaSection+exactFaceSection triangulation face verticesValue@(firstVertex, secondVertex, thirdVertex) =+ let Point ax ay = Dcel.vertexPoint triangulation firstVertex+ Point bx by = Dcel.vertexPoint triangulation secondVertex+ Point cx cy = Dcel.vertexPoint triangulation thirdVertex+ in case exactCircumradiusSquaredRational ax ay bx by cx cy of+ Left ExactZeroDivisor -> Left (AlphaFaceDegenerate face)+ Left arithmeticError -> Left (AlphaCircumradiusDivisionFailed face arithmeticError)+ Right exactBirth ->+ Right+ FaceAlphaSection+ { faceAlphaBirth = alphaBirthFromExact exactBirth+ , faceAlphaVertices = verticesValue+ }++exactEdgeBirth+ :: Vector FaceAlphaSection+ -> DelaunayTriangulation vertex+ -> UndirectedEdgeId+ -> Either AlphaFiltrationError AlphaBirth+exactEdgeBirth faceSections triangulation edge = do+ let (fromVertex, toVertex) = Dcel.undirectedEndpoints triangulation edge+ fromPoint = Dcel.vertexPoint triangulation fromVertex+ toPoint = Dcel.vertexPoint triangulation toVertex+ incidentFaces = innerIncidentFaces triangulation edge+ oppositeVertices <-+ traverse+ (oppositeVertexAcross faceSections triangulation edge fromVertex toVertex)+ incidentFaces+ let oppositePoints = fmap (Dcel.vertexPoint triangulation) oppositeVertices+ if all (diametralWitnessOutside fromPoint toPoint) oppositePoints+ then pure (alphaBirthFromExact (quarterSquaredPointDistance fromPoint toPoint))+ else leastIncidentFaceBirth edge faceSections incidentFaces++innerIncidentFaces+ :: DelaunayTriangulation vertex+ -> UndirectedEdgeId+ -> [FaceId]+innerIncidentFaces triangulation edge =+ let (forward, backward) = directedPair edge+ forwardFace = Dcel.incidentFace triangulation forward+ backwardFace = Dcel.incidentFace triangulation backward+ in case (forwardFace == Dcel.outerFace, backwardFace == Dcel.outerFace) of+ (True, True) -> []+ (False, True) -> [forwardFace]+ (True, False) -> [backwardFace]+ (False, False)+ | forwardFace == backwardFace -> [forwardFace]+ | forwardFace < backwardFace -> [forwardFace, backwardFace]+ | otherwise -> [backwardFace, forwardFace]++oppositeVertexAcross+ :: Vector FaceAlphaSection+ -> DelaunayTriangulation vertex+ -> UndirectedEdgeId+ -> VertexId+ -> VertexId+ -> FaceId+ -> Either AlphaFiltrationError VertexId+oppositeVertexAcross faceSections triangulation edge fromVertex toVertex face =+ case faceSectionAt faceSections face of+ Nothing ->+ Left+ (AlphaFaceEdgeIncidenceInvalid face edge (Dcel.faceVertices triangulation face))+ Just faceSection ->+ let verticesValue@(firstVertex, secondVertex, thirdVertex) = faceAlphaVertices faceSection+ outsideEdge vertex = vertex /= fromVertex && vertex /= toVertex+ in case+ ( outsideEdge firstVertex+ , outsideEdge secondVertex+ , outsideEdge thirdVertex+ )+ of+ (True, False, False) -> Right firstVertex+ (False, True, False) -> Right secondVertex+ (False, False, True) -> Right thirdVertex+ _ ->+ Left+ ( AlphaFaceEdgeIncidenceInvalid+ face+ edge+ (filter outsideEdge (triangleVertices verticesValue))+ )++diametralWitnessOutside :: Point -> Point -> Point -> Bool+diametralWitnessOutside firstPoint secondPoint =+ not . inDiametralCircle firstPoint secondPoint++quarterSquaredPointDistance :: Point -> Point -> ExactRational+quarterSquaredPointDistance (Point ax ay) (Point bx by) =+ exactQuarterSquaredDistanceRational ax ay bx by++leastIncidentFaceBirth+ :: UndirectedEdgeId+ -> Vector FaceAlphaSection+ -> [FaceId]+ -> Either AlphaFiltrationError AlphaBirth+leastIncidentFaceBirth edge faceSections incidentFaces = do+ births <-+ traverse+ (\face -> maybe (Left (AlphaIncidentFaceBirthMissing edge face)) (Right . faceAlphaBirth) (faceSectionAt faceSections face))+ incidentFaces+ case NonEmpty.nonEmpty births of+ Nothing -> Left (AlphaEdgeBirthWitnessMissing edge)+ Just nonEmptyBirths -> Right (minimum nonEmptyBirths)++faceSectionAt :: Vector FaceAlphaSection -> FaceId -> Maybe FaceAlphaSection+faceSectionAt faceSections (FaceId rawFace)+ | rawFace == 0 = Nothing+ | otherwise = faceSections Vector.!? (fromIntegral rawFace - 1)++faceBirthAt :: Vector AlphaBirth -> FaceId -> Maybe AlphaBirth+faceBirthAt faceBirths (FaceId rawFace)+ | rawFace == 0 = Nothing+ | otherwise = faceBirths Vector.!? (fromIntegral rawFace - 1)++triangleVertices :: (VertexId, VertexId, VertexId) -> [VertexId]+triangleVertices (firstVertex, secondVertex, thirdVertex) =+ [firstVertex, secondVertex, thirdVertex]++alphaBirthFromExact :: ExactRational -> AlphaBirth+alphaBirthFromExact exactValue =+ AlphaBirth (alphaOrderKey exactValue) exactValue++-- | A monotone exact fixed-point projection of the positive rational. Distinct+-- keys prove order; collisions descend to the authoritative rational. One+-- division at birth replaces repeated cross-products during every downstream+-- ordered-map comparison.+alphaOrderKey :: ExactRational -> AlphaOrderKey+alphaOrderKey exactValue =+ let numerator = exactRationalNumerator exactValue+ denominator = exactRationalDenominator exactValue+ in if numerator <= 0+ then NonPositiveAlphaOrder+ else+ let candidateExponent = integerBitLength numerator - integerBitLength denominator+ binaryExponent =+ if candidateExponent >= 0+ then+ if numerator < denominator `shiftL` candidateExponent+ then candidateExponent - 1+ else candidateExponent+ else+ if numerator `shiftL` negate candidateExponent < denominator+ then candidateExponent - 1+ else candidateExponent+ scale = alphaOrderFractionBits - binaryExponent+ scaledMantissa =+ if scale >= 0+ then (numerator `shiftL` scale) `quot` denominator+ else numerator `quot` (denominator `shiftL` negate scale)+ in PositiveAlphaOrder binaryExponent (fromIntegral scaledMantissa)++alphaOrderFractionBits :: Int+alphaOrderFractionBits = 16
src-dcel/Moonlight/Triangulation/FloodFillIterator.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE DataKinds #-} {-# LANGUAGE DeriveAnyClass #-} {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE DerivingStrategies #-}@@ -9,7 +8,6 @@ ( DistanceMetric (..) , CircleMetric , CircleMetricError (..)- , RadiusSquaredError (..) , RectangleMetric , RectangleMetricError (..) , circleMetric@@ -38,9 +36,6 @@ , labelledRegionBoundaries , componentBoundaryLoops , componentBoundary- , RadiusSquared- , mkRadiusSquared- , alphaShapeContainsFace ) where import Control.DeepSeq (NFData)@@ -65,7 +60,6 @@ ) import Moonlight.Triangulation.Math import Moonlight.Triangulation.PointLocation-import Moonlight.Triangulation.Scalar (circumradiusSquaredWithinCoordinates) import Moonlight.Triangulation.Types -- | One non-empty connected set of equally labelled bounded face indices.@@ -122,26 +116,6 @@ !DirectedEdgeId !DirectedEdgeId --- | An admitted finite, non-negative squared radius.-newtype RadiusSquared = RadiusSquared Double- deriving stock (Eq, Ord, Show, Generic)- deriving anyclass (NFData)---- | Typed refusal shared by every squared-radius constructor.-data RadiusSquaredError- = NonFiniteRadiusSquared !NonFiniteValue- | NegativeRadiusSquared !Double- deriving stock (Eq, Ord, Show)---- | Admit a finite, non-negative squared radius.-mkRadiusSquared :: Double -> Either RadiusSquaredError RadiusSquared-mkRadiusSquared value =- case classifyNonFinite value of- Just nonFinite -> Left (NonFiniteRadiusSquared nonFinite)- Nothing- | value < 0 -> Left (NegativeRadiusSquared value)- | otherwise -> Right (RadiusSquared value)- -- | A query shape that can admit points, test edges, and supply a location -- seed. class DistanceMetric metric where@@ -190,10 +164,10 @@ Right (RectangleMetric queryLower queryUpper queryCenter) instance DistanceMetric CircleMetric where- metricContainsPoint (CircleMetric center (RadiusSquared radiusSquared)) point =- squaredDistanceWide (queryPointValue center) point <= radiusSquared- metricIntersectsEdge (CircleMetric center (RadiusSquared radiusSquared)) from to =- segmentDistanceSquaredWide from to (queryPointValue center) <= radiusSquared+ metricContainsPoint (CircleMetric center radius) point =+ squaredDistanceWide (queryPointValue center) point <= radiusSquaredValue radius+ metricIntersectsEdge (CircleMetric center radius) from to =+ segmentDistanceSquaredWide from to (queryPointValue center) <= radiusSquaredValue radius metricStartPoint (CircleMetric center _) = center instance DistanceMetric RectangleMetric where@@ -612,22 +586,6 @@ winding previousVertex current nextVertex = orient2d (point previousVertex) (point current) (point nextVertex) key vertex = (point vertex, vertex)---- | Membership of a bounded face in the closed alpha shape. Exact dyadic--- comparison makes equality independent of circumcenter rounding.-alphaShapeContainsFace- :: RadiusSquared- -> Triangulation 'Unconstrained vertex directed undirected face- -> FaceId- -> Bool-alphaShapeContainsFace (RadiusSquared threshold) triangulation =- maybe False withinRadius . innerFaceVertices triangulation- where- withinRadius (firstVertex, secondVertex, thirdVertex) =- let Point ax ay = vertexPoint triangulation firstVertex- Point bx by = vertexPoint triangulation secondVertex- Point cx cy = vertexPoint triangulation thirdVertex- in circumradiusSquaredWithinCoordinates threshold ax ay bx by cx cy -- | Inner faces separated from the outer face by an even minimum number of -- barriers. A 0–1 BFS floods freely within one depth before crossing a barrier,
src-dcel/Moonlight/Triangulation/Internal/CellSet.hs view
@@ -10,6 +10,7 @@ ( ExactCellSet (..) , CellSelectionError (..) , exactCellSet+ , residentExactCellSet , closeFaceCellSet , closeExactCellSetWith , exactCellSetVertexCount@@ -91,6 +92,33 @@ (ordinaryExactPoint triangulation) triangulation (ClosedCellIds verticesSet edgesSet facesSet)++-- | Seal the complete resident DCEL without rechecking ranges or closure+-- already established by its dense handle families. Exact coordinates are+-- still admitted at this representation boundary; only impossible handle+-- refusals and repeated incidence walks disappear.+residentExactCellSet+ :: Triangulation mode vertex directed undirected face+ -> Either CellSelectionError ExactCellSet+residentExactCellSet triangulation = do+ let vertexIndices = [0 .. numVertices triangulation - 1]+ edgeIndices = [0 .. numUndirectedEdges triangulation - 1]+ faceIndices = [1 .. numFaces triangulation - 1]+ exactPoints <-+ IntMap.fromDistinctAscList+ <$> traverse+ ( \vertexIndex -> do+ point <- ordinaryExactPoint triangulation (VertexId (fromIntegral vertexIndex))+ pure (vertexIndex, point)+ )+ vertexIndices+ pure+ ( ExactCellSet+ triangulation+ exactPoints+ (IntSet.fromDistinctAscList edgeIndices)+ (IntSet.fromDistinctAscList faceIndices)+ ) -- | Close a bounded face selection over all resident boundary cells. Face -- handles are author input and are checked; constructed closure is not then
src-dcel/Moonlight/Triangulation/Internal/Mutable.hs view
@@ -18,6 +18,7 @@ , DefaultedVertexDcel , newMutableDcelWithVertexDefault , defaultedVertexDcel+ , defaultedDenseVertexDcel , thawTriangulation , thawTriangulationDense , freezeTriangulation@@ -84,6 +85,8 @@ , readConstraint , denseReadPointX , denseReadPointY+ , denseWriteFreshDefaultPoint+ , denseCommitFreshDefaultPoints , denseReadOrigin , denseWriteOrigin , denseReadNext@@ -362,26 +365,32 @@ newMutableDcel :: ElementDefaults directed undirected face -> DcelCapacity -> ST s (MutableDcel s vertex directed undirected face) newMutableDcel defaults capacity = newMutableDcelFrom DenseTransaction Nothing defaults capacity Nothing --- | A fresh mutable DCEL whose vertex payload plane is uniformly filled.--- The witness is what permits geometry-only ingress to extend that plane--- without materializing one boxed unit value per site.-newtype DefaultedVertexDcel s vertex directed undirected face = DefaultedVertexDcel- { defaultedVertexDcel :: MutableDcel s vertex directed undirected face+-- | A fresh dense mutable DCEL whose vertex payload plane is uniformly filled.+-- The witness keeps the flat section tied to its canonical owner and permits+-- geometry-only ingress to extend the payload plane without materializing one+-- boxed value per site.+data DefaultedVertexDcel s vertex directed undirected face = DefaultedVertexDcel+ { defaultedVertexDcel :: !(MutableDcel s vertex directed undirected face)+ , defaultedDenseVertexDcel :: !(DenseMutableDcel s vertex directed undirected face) } newMutableDcelWithVertexDefault :: vertex -> ElementDefaults directed undirected face -> DcelCapacity- -> ST s (DefaultedVertexDcel s vertex directed undirected face)-newMutableDcelWithVertexDefault vertexDefault defaults capacity =- DefaultedVertexDcel- <$> newMutableDcelFrom- DenseTransaction- (Just vertexDefault)- defaults- capacity- Nothing+ -> ST s (Either BuildError (DefaultedVertexDcel s vertex directed undirected face))+newMutableDcelWithVertexDefault vertexDefault defaults capacity = do+ mutable <-+ newMutableDcelFrom+ DenseTransaction+ (Just vertexDefault)+ defaults+ capacity+ Nothing+ pure $+ case denseMutableDcel mutable of+ Nothing -> Left CircleSweepRequiresDenseStorage+ Just dense -> Right (DefaultedVertexDcel mutable dense) -- | Open a local-edit transaction: copy-on-write pages, publication -- proportional to dirtied pages. The section for singleton persistent verbs.@@ -788,7 +797,7 @@ -> Double -> Double -> ST s Int-appendDefaultVertexCoordinates (DefaultedVertexDcel mutable) x y = do+appendDefaultVertexCoordinates DefaultedVertexDcel{defaultedVertexDcel = mutable} x y = do slot <- nextVertexSlot mutable initializeVertexCoordinatesAtSlot mutable slot x y {-# INLINE appendDefaultVertexCoordinates #-}@@ -1326,6 +1335,51 @@ denseReadPointY DenseMutableDcel{dmdPointY} = readFlatMutable dmdPointY {-# INLINE denseReadPointX #-} {-# INLINE denseReadPointY #-}++-- | Materialize one coordinate slot in a fresh dense point section. The+-- caller owns the capacity witness and supplies the final cardinality once;+-- the point index remains the fresh owner's missing derived view until+-- publication rebuilds it from these authoritative coordinate planes.+denseWriteFreshDefaultPoint+ :: DefaultedVertexDcel s vertex directed undirected face+ -> Int+ -> Double+ -> Double+ -> ST s ()+denseWriteFreshDefaultPoint+ DefaultedVertexDcel+ { defaultedDenseVertexDcel =+ dense@DenseMutableDcel+ { dmdOwner = MutableDcel{mdInitialPointCount}+ , dmdPointX+ , dmdPointY+ , dmdNewConnected+ }+ }+ vertex+ x+ y = do+ writeFlatMutable dmdPointX vertex x+ writeFlatMutable dmdPointY vertex y+ denseWriteVertexOut dense vertex (-1)+ writeFlatMutable dmdNewConnected (vertex - mdInitialPointCount) 0+{-# INLINE denseWriteFreshDefaultPoint #-}++-- | Glue a completed fresh point section into the mutable owner's logical+-- extent. Per-site writes remain local; the aggregate cardinality descends+-- once after the section is complete.+denseCommitFreshDefaultPoints+ :: DefaultedVertexDcel s vertex directed undirected face+ -> Int+ -> ST s ()+denseCommitFreshDefaultPoints+ DefaultedVertexDcel+ { defaultedDenseVertexDcel = DenseMutableDcel+ { dmdOwner = MutableDcel{mdPointCount}+ }+ }+ = writeSTRef mdPointCount+{-# INLINE denseCommitFreshDefaultPoints #-} denseReadOrigin :: DenseMutableDcel s vertex directed undirected face -> Int -> ST s Int denseReadOrigin DenseMutableDcel{dmdHalfTopology} index = fromIntegral <$> readFlatMutable dmdHalfTopology (4 * index)
src-dcel/Moonlight/Triangulation/Internal/Types.hs view
@@ -31,6 +31,10 @@ , CoordinateError (..) , NonFiniteValue (..) , classifyNonFinite+ , RadiusSquared+ , RadiusSquaredError (..)+ , mkRadiusSquared+ , radiusSquaredValue , BuildError (..) , Location (..) , LocationHint (..)@@ -229,6 +233,33 @@ | isInfinite value && value < 0 = Just ValueNegativeInfinity | isInfinite value = Just ValuePositiveInfinity | otherwise = Nothing++-- | An admitted finite, non-negative squared radius. This belongs to the+-- shared geometric vocabulary because both circle queries and alpha+-- filtration consume it.+newtype RadiusSquared = RadiusSquared Double+ deriving stock (Eq, Ord, Show, Generic)+ deriving anyclass (NFData)++-- | Typed refusal shared by every squared-radius consumer.+data RadiusSquaredError+ = NonFiniteRadiusSquared !NonFiniteValue+ | NegativeRadiusSquared !Double+ deriving stock (Eq, Ord, Show, Generic)+ deriving anyclass (NFData)++-- | Admit a finite, non-negative squared radius.+mkRadiusSquared :: Double -> Either RadiusSquaredError RadiusSquared+mkRadiusSquared value =+ case classifyNonFinite value of+ Just nonFinite -> Left (NonFiniteRadiusSquared nonFinite)+ Nothing+ | value < 0 -> Left (NegativeRadiusSquared value)+ | otherwise -> Right (RadiusSquared value)++-- | The admitted binary64 threshold.+radiusSquaredValue :: RadiusSquared -> Double+radiusSquaredValue (RadiusSquared value) = value -- | Total construction and rewrite obstruction surface. data BuildError
src-dcel/Moonlight/Triangulation/Types.hs view
@@ -22,6 +22,10 @@ , CoordinateError (..) , NonFiniteValue (..) , classifyNonFinite+ , RadiusSquared+ , RadiusSquaredError (..)+ , mkRadiusSquared+ , radiusSquaredValue , BuildError (..) , Location (..) , LocationHint (..)
src-dcel/Moonlight/Triangulation/Validation.hs view
@@ -76,7 +76,10 @@ , flag /= 0 && flag /= 1 ] ++ [CachedConstraintCountMismatch (triConstraintCount triangulation) actualConstraintCount | triConstraintCount triangulation /= actualConstraintCount]- ++ [CachedConstraintIndexMismatch | triConstraintEdges triangulation /= indexedConstraintEdges]+ ++ [ CachedConstraintIndexMismatch+ | constraintLength == edgeCount+ , triConstraintEdges triangulation /= indexedConstraintEdges+ ] ++ [MissingOuterFace | facesCount == 0] ++ [FacePayloadLengthMismatch facePayloadCount facesCount | facePayloadCount /= facesCount] @@ -134,7 +137,7 @@ ] edgeViolations- | not halfArraysEqual || odd halfCount = []+ | not halfArraysEqual || odd halfCount || not (null rangeViolations) = [] | otherwise = concatMap validateEdge (directedEdges triangulation) validateEdge edge@(DirectedEdgeId raw) =@@ -163,7 +166,9 @@ validEdge (DirectedEdgeId value) = fromIntegral value < halfCount - faceViolations = concatMap validateFace (allFaces triangulation)+ faceViolations+ | not (null rangeViolations) = []+ | otherwise = concatMap validateFace (allFaces triangulation) validateFace face@(FaceId _) = case adjacentEdge triangulation face of Nothing@@ -180,7 +185,9 @@ , length faceVertexIds /= 3 || length (nub faceVertexIds) /= 3 ] - vertexViolations = concatMap validateVertex (vertices triangulation)+ vertexViolations+ | vertexOutCount /= verticesCount || not (null rangeViolations) = []+ | otherwise = concatMap validateVertex (vertices triangulation) validateVertex vertex = case vertexOutEdge triangulation vertex of Nothing | verticesCount <= 1 -> []
src-public/Moonlight/Triangulation.hs view
@@ -202,6 +202,22 @@ , RadiusSquared , RadiusSquaredError (..) , mkRadiusSquared+ , radiusSquaredValue+ , AlphaBirth+ , alphaBirthNumerator+ , alphaBirthDenominator+ , alphaBirthToDouble+ , alphaBirthFromRadiusSquared+ , AlphaFiltration+ , AlphaFiltrationError (..)+ , alphaFiltration+ , alphaFiltrationCellSet+ , alphaVertexBirth+ , alphaEdgeBirth+ , alphaFaceBirth+ , alphaFiltrationCriticalBirths+ , alphaComplexAtBirth+ , alphaComplexAtRadius , alphaShapeContainsFace -- * Exact planar regions — authoritative rational geometry, labelled@@ -358,15 +374,30 @@ , vertexPoint , vertexPoints )+import Moonlight.Triangulation.Alpha+ ( AlphaBirth+ , AlphaFiltration+ , AlphaFiltrationError (..)+ , alphaBirthDenominator+ , alphaBirthFromRadiusSquared+ , alphaBirthNumerator+ , alphaBirthToDouble+ , alphaComplexAtBirth+ , alphaComplexAtRadius+ , alphaEdgeBirth+ , alphaFaceBirth+ , alphaFiltration+ , alphaFiltrationCellSet+ , alphaFiltrationCriticalBirths+ , alphaShapeContainsFace+ , alphaVertexBirth+ ) import Moonlight.Triangulation.FloodFillIterator ( BoundaryLoop , BoundaryObstruction (..) , BoundaryOrientation (..) , FaceComponent- , RadiusSquared , RegionBoundary- , RadiusSquaredError (..)- , alphaShapeContainsFace , boundaryLoopOrientation , boundaryLoopVertices , componentBoundary@@ -374,7 +405,6 @@ , faceComponentFaces , faceComponents , facesAtEvenBarrierDepth- , mkRadiusSquared , regionBoundaryHoleLoops , regionBoundaryOuterLoop )@@ -507,10 +537,14 @@ , Point (..) , PointValidationError (..) , QueryPoint+ , RadiusSquared+ , RadiusSquaredError (..) , RefinementParameters (..) , SiteRelation (..) , defaultRefinementParameters+ , mkRadiusSquared , queryPointValue+ , radiusSquaredValue , unitElementDefaults ) import Moonlight.Triangulation.Interpolation
src-serialize/Moonlight/Triangulation/Serialization.hs view
@@ -60,29 +60,13 @@ import Moonlight.Triangulation.Handles.HandleDefs import Moonlight.Triangulation.Internal.Representation import Moonlight.Triangulation.Internal.Types-import Moonlight.Triangulation.Math (mkQueryPoint)+import Moonlight.Triangulation.Math (canonicalCoordinate, mkQueryPoint) import Moonlight.Triangulation.Validation (validateTriangulation) instance Binary (Point) where put (Point x y) = putDoublebe x >> putDoublebe y get = Point <$> getDoublebe <*> getDoublebe -instance Binary VertexId where- put (VertexId value) = putWord32be value- get = VertexId <$> getWord32be--instance Binary FaceId where- put (FaceId value) = putWord32be value- get = FaceId <$> getWord32be--instance Binary DirectedEdgeId where- put (DirectedEdgeId value) = putWord32be value- get = DirectedEdgeId <$> getWord32be--instance Binary UndirectedEdgeId where- put (UndirectedEdgeId value) = putWord32be value- get = UndirectedEdgeId <$> getWord32be- -- | The finite resource envelope admitted by the canonical decoder. The byte -- budget bounds the complete input before parsing; the element budget bounds -- the total number of library-owned serialized section elements before any@@ -293,8 +277,8 @@ faceCount = fromIntegral encodedFaceCount cachedConstraintCount = fromIntegral encodedConstraintCount defaults <- ElementDefaults <$> lift get <*> lift get <*> lift get- pointXs <- U.replicateM vertexCount (lift getDoublebe)- pointYs <- U.replicateM vertexCount (lift getDoublebe)+ pointXs <- U.replicateM vertexCount (canonicalCoordinate <$> lift getDoublebe)+ pointYs <- U.replicateM vertexCount (canonicalCoordinate <$> lift getDoublebe) vertexDefault <- lift get vertexDataVector <- V.replicateM vertexCount (lift get) vertexOut <- U.replicateM vertexCount (lift getWord32be)
test/cell-complex/Main.hs view
@@ -1,8 +1,36 @@ module Main (main) where +import Data.Bits (shiftL)+import Data.Function ((&)) import Data.Foldable (traverse_)+import Data.List qualified as List+import Data.Map.Strict qualified as Map+import Data.Maybe (mapMaybe)+import Data.Ratio ((%)) import Data.Vector qualified as Vector-import Moonlight.Triangulation.CellComplex (DCELComplex, fromExactCellSet)+import Moonlight.Homology.Boundary+ ( BoundaryEntry+ , boundaryCoefficient+ , boundaryEntries+ , degreeCardinality+ , incidenceMatrixAt+ , sourceIndex+ , targetIndex+ )+import Moonlight.Homology.Chain+ ( HomologicalDegree (..)+ , PersistencePair (..)+ )+import Moonlight.Homology.Persistence+ ( FilteredFiniteChainComplex+ , filteredBaseComplex+ , filteredCellBirths+ , mod2PersistentPairs+ )+import Moonlight.Homology.Topology+ ( BasisCellRef (..)+ , freeBettiVector+ ) import Moonlight.Homology.Pure.Topology.CellComplex ( CellComplex2D (..), CellTypes (..),@@ -11,6 +39,17 @@ eulerCharacteristic, isBoundaryEdge, )+import Moonlight.Triangulation.Alpha+ ( AlphaBirth+ , AlphaFiltration+ , alphaBirthDenominator+ , alphaBirthNumerator+ , alphaEdgeBirth+ , alphaFaceBirth+ , alphaFiltration+ , alphaFiltrationCellSet+ , alphaVertexBirth+ ) import Moonlight.Triangulation.BulkLoad (delaunayGeometry) import Moonlight.Triangulation.CellSet ( ExactCellSet,@@ -19,8 +58,34 @@ exactCellSetFaceCount, exactCellSetVertexCount, )-import Moonlight.Triangulation.Handles.Iterators.FixedIterators (innerFaces)-import Moonlight.Triangulation.Types (Point (..))+import Moonlight.Triangulation.CellComplex+ ( DCELComplex+ , filteredAlphaComplex+ , finiteChainComplex+ , fromExactCellSet+ )+import Moonlight.Triangulation.Dcel+ ( undirectedEndpoints+ , vertexPoint+ )+import Moonlight.Triangulation.Handles.HandleDefs+ ( UndirectedEdgeId+ )+import Moonlight.Triangulation.Handles.Iterators.FixedIterators+ ( innerFaces+ , undirectedEdges+ )+import Moonlight.Triangulation.Internal.Dyadic (exactDiametralDot)+import Moonlight.Triangulation.Internal.ExactRational+ ( exactRational+ , exactRationalFromDyadic+ , exactRationalFromDyadicRatio+ )+import Moonlight.Triangulation.Math (inDiametralCircle)+import Moonlight.Triangulation.Types+ ( DelaunayTriangulation+ , Point (..)+ ) import Test.Tasty (TestTree, defaultMain, testGroup) import Test.Tasty.HUnit ((@?=), Assertion, assertBool, assertFailure, testCase) @@ -33,7 +98,15 @@ "ExactCellSet bridge" [ testCase "preserves the admitted cell inventory" preserveCellInventory, testCase "preserves downward-closed incidence" preserveClosedIncidence,- testCase "marks the triangular exterior as absent" preserveExteriorAdjacency+ testCase "marks the triangular exterior as absent" preserveExteriorAdjacency,+ testCase "lowers oriented DCEL incidence to a checked integral chain complex" lowerIntegralChainComplex,+ testCase "assigns exact Gabriel and non-Gabriel alpha births" exactAlphaBirths,+ testCase "preserves exact dyadic rational normalization" exactDyadicRationalNormalization,+ testCase "orders cached alpha keys by their exact rational births" alphaBirthOrderMatchesExactRationals,+ testCase "agrees with the exact diametral predicate" diametralPredicateMatchesExactSign,+ testCase "keeps every boundary cell no later than its coface" alphaBirthsAreFaceMonotone,+ testCase "finds the square alpha hole from radius one to radius two" squareAlphaPersistence,+ testCase "keeps point and collinear alpha filtrations total" degenerateSupportFiltrations ] preserveCellInventory :: Assertion@@ -60,6 +133,250 @@ assertBool "every edge of a single selected triangle has one exterior incident face" (all (isBoundaryEdge complexValue) (edges complexValue))++lowerIntegralChainComplex :: Assertion+lowerIntegralChainComplex =+ withTriangleComplex $ \_ complexValue -> do+ finite <- requireRight "triangle finite chain lowering" (finiteChainComplex complexValue)+ freeBettiVector finite @?= [1, 0, 0]++exactAlphaBirths :: Assertion+exactAlphaBirths = do+ triangulation <- requireDelaunay "obtuse triangle" obtuseTrianglePoints+ filtration <- requireRight "obtuse triangle alpha filtration" (alphaFiltration triangulation)+ let vertexBirthRatios =+ vertices (fromExactCellSet (alphaFiltrationCellSet filtration))+ & fmap (fmap alphaBirthRatio . alphaVertexBirth filtration)+ faceBirthRatios =+ faces (fromExactCellSet (alphaFiltrationCellSet filtration))+ & fmap (fmap alphaBirthRatio . alphaFaceBirth filtration)+ assertBool "every vertex is born at zero" (all (== Just (0, 1)) vertexBirthRatios)+ faceBirthRatios @?= [Just (25, 4)]+ longEdge <- requireSome "obtuse triangle long edge" (edgeBetween triangulation (Point (-2) 0) (Point 2 0))+ fmap alphaBirthRatio (alphaEdgeBirth filtration longEdge) @?= Just (25, 4)+ let shortEdgeBirths =+ edges (fromExactCellSet (alphaFiltrationCellSet filtration))+ & filter (/= longEdge)+ & fmap (fmap alphaBirthRatio . alphaEdgeBirth filtration)+ shortEdgeBirths @?= [Just (5, 4), Just (5, 4)]++exactDyadicRationalNormalization :: Assertion+exactDyadicRationalNormalization = do+ traverse_ assertDyadic+ [ (0, -200)+ , (12, -5)+ , (-12, -5)+ , (3, 4)+ , (2 ^ (80 :: Int) + 8, -70)+ ]+ traverse_ assertDyadicRatio+ [ (0, 7, -200)+ , (96, 40, -11)+ , (-96, 40, 9)+ , (45, -28, -3)+ , (2 ^ (120 :: Int) + 24, 2 ^ (75 :: Int) + 12, -51)+ ]+ where+ assertDyadic (numerator, power) =+ Right (exactRationalFromDyadic numerator power)+ @?= ( if power >= 0+ then exactRational (numerator `shiftL` power) 1+ else exactRational numerator (1 `shiftL` negate power)+ )+ assertDyadicRatio (numerator, denominator, power) =+ exactRationalFromDyadicRatio numerator denominator power+ @?= ( if power >= 0+ then exactRational (numerator `shiftL` power) denominator+ else exactRational numerator (denominator `shiftL` negate power)+ )++alphaBirthOrderMatchesExactRationals :: Assertion+alphaBirthOrderMatchesExactRationals = do+ filtrations <-+ traverse+ (\(label, points) -> requireDelaunay label points >>= requireRight (label <> " alpha filtration") . alphaFiltration)+ [ ("exact-order obtuse triangle", obtuseTrianglePoints)+ , ("exact-order square", squarePoints)+ ]+ let births = filtrations >>= allAlphaBirths+ traverse_+ ( \(leftBirth, rightBirth) ->+ compare leftBirth rightBirth+ @?= compare (alphaBirthRational leftBirth) (alphaBirthRational rightBirth)+ )+ [(leftBirth, rightBirth) | leftBirth <- births, rightBirth <- births]++diametralPredicateMatchesExactSign :: Assertion+diametralPredicateMatchesExactSign =+ traverse_+ ( \(firstPoint, secondPoint, witnessPoint) ->+ inDiametralCircle firstPoint secondPoint witnessPoint+ @?= exactDiametralMembership firstPoint secondPoint witnessPoint+ )+ [ (Point (-1) 0, Point 1 0, Point 0 0)+ , (Point (-1) 0, Point 1 0, Point 0 2)+ , (Point (-1) 0, Point 1 0, Point 0 1)+ , (Point (-1) 0, Point 1 0, Point 0 (1 + encodeFloat 1 (-52)))+ , (Point (-1e40) 0, Point 1e40 0, Point 0 1e40)+ ]++allAlphaBirths :: AlphaFiltration -> [AlphaBirth]+allAlphaBirths filtration =+ let complexValue = fromExactCellSet (alphaFiltrationCellSet filtration)+ in mapMaybe (alphaVertexBirth filtration) (vertices complexValue)+ <> mapMaybe (alphaEdgeBirth filtration) (edges complexValue)+ <> mapMaybe (alphaFaceBirth filtration) (faces complexValue)++alphaBirthRational :: AlphaBirth -> Rational+alphaBirthRational birth =+ alphaBirthNumerator birth % alphaBirthDenominator birth++exactDiametralMembership :: Point -> Point -> Point -> Bool+exactDiametralMembership (Point ax ay) (Point bx by) (Point px py) =+ exactDiametralDot ax ay bx by px py <= 0++alphaBirthsAreFaceMonotone :: Assertion+alphaBirthsAreFaceMonotone = do+ triangulation <- requireDelaunay "alpha monotonicity square" squarePoints+ filtration <- requireRight "alpha monotonicity filtration" (alphaFiltration triangulation)+ filtered <- requireRight "alpha monotonicity filtered complex" (filteredAlphaComplex filtration)+ traverse_+ (assertBoundaryDegreeMonotone filtered)+ [HomologicalDegree 1, HomologicalDegree 2]++assertBoundaryDegreeMonotone+ :: FilteredFiniteChainComplex AlphaBirth Int+ -> HomologicalDegree+ -> Assertion+assertBoundaryDegreeMonotone filtered degreeValue =+ traverse_+ (assertBoundaryEntryMonotone filtered degreeValue)+ ( boundaryEntries+ (incidenceMatrixAt (filteredBaseComplex filtered) degreeValue)+ )++assertBoundaryEntryMonotone+ :: FilteredFiniteChainComplex AlphaBirth Int+ -> HomologicalDegree+ -> BoundaryEntry Int+ -> Assertion+assertBoundaryEntryMonotone filtered degreeValue@(HomologicalDegree degreeIndex) entryValue = do+ let sourceCell =+ BasisCellRef+ { cellDegree = degreeValue+ , cellIndex = sourceIndex entryValue+ }+ targetCell =+ BasisCellRef+ { cellDegree = HomologicalDegree (degreeIndex - 1)+ , cellIndex = targetIndex entryValue+ }+ births = filteredCellBirths filtered+ sourceBirth <- requireSome "source alpha birth" (Map.lookup sourceCell births)+ targetBirth <- requireSome "boundary alpha birth" (Map.lookup targetCell births)+ assertBool+ ("boundary birth exceeds coface birth for " <> show (targetCell, sourceCell, boundaryCoefficient entryValue))+ (targetBirth <= sourceBirth)++squareAlphaPersistence :: Assertion+squareAlphaPersistence = do+ triangulation <- requireDelaunay "cocircular square" squarePoints+ filtration <- requireRight "cocircular square alpha filtration" (alphaFiltration triangulation)+ filtered <- requireRight "cocircular square filtered chain complex" (filteredAlphaComplex filtration)+ pairs <- requireRight "cocircular square persistence" (mod2PersistentPairs filtered)+ let positiveDegreeOnePairs =+ pairs+ & filter ((== HomologicalDegree 1) . persistenceDegree)+ & mapMaybe positiveBirthAndDeath+ positiveDegreeOnePairs @?= [((1, 1), (2, 1))]++degenerateSupportFiltrations :: Assertion+degenerateSupportFiltrations =+ traverse_+ assertSupport+ [ ("two sites", Vector.fromList [Point 0 0, Point 2 0], (2, 1, 0))+ , ("three collinear sites", Vector.fromList [Point 0 0, Point 1 0, Point 3 0], (3, 2, 0))+ ]+ where+ assertSupport (label, points, expectedCounts) = do+ triangulation <- requireDelaunay label points+ filtration <- requireRight (label <> " alpha filtration") (alphaFiltration triangulation)+ filtered <- requireRight (label <> " filtered chain complex") (filteredAlphaComplex filtration)+ let finite = filteredBaseComplex filtered+ counts =+ ( degreeCardinality finite (HomologicalDegree 0)+ , degreeCardinality finite (HomologicalDegree 1)+ , degreeCardinality finite (HomologicalDegree 2)+ )+ counts @?= expectedCounts++positiveBirthAndDeath+ :: PersistencePair AlphaBirth+ -> Maybe ((Integer, Integer), (Integer, Integer))+positiveBirthAndDeath pairValue =+ case persistenceDeath pairValue of+ Just deathBirth+ | persistenceBirth pairValue < deathBirth ->+ Just+ ( alphaBirthRatio (persistenceBirth pairValue)+ , alphaBirthRatio deathBirth+ )+ _ -> Nothing++alphaBirthRatio :: AlphaBirth -> (Integer, Integer)+alphaBirthRatio birth =+ (alphaBirthNumerator birth, alphaBirthDenominator birth)++edgeBetween+ :: DelaunayTriangulation ()+ -> Point+ -> Point+ -> Maybe UndirectedEdgeId+edgeBetween triangulation firstPoint secondPoint =+ List.find hasEndpoints (undirectedEdges triangulation)+ where+ hasEndpoints edgeValue =+ let (fromVertex, toVertex) = undirectedEndpoints triangulation edgeValue+ fromPoint = vertexPoint triangulation fromVertex+ toPoint = vertexPoint triangulation toVertex+ in (fromPoint == firstPoint && toPoint == secondPoint)+ || (fromPoint == secondPoint && toPoint == firstPoint)++requireDelaunay+ :: String+ -> Vector.Vector Point+ -> IO (DelaunayTriangulation ())+requireDelaunay label points =+ requireRight (label <> " Delaunay construction") (delaunayGeometry points)++requireRight :: Show failure => String -> Either failure value -> IO value+requireRight label result =+ case result of+ Left failureValue -> assertFailure (label <> " failed: " <> show failureValue)+ Right value -> pure value++requireSome :: String -> Maybe value -> IO value+requireSome label maybeValue =+ case maybeValue of+ Nothing -> assertFailure (label <> " was absent")+ Just value -> pure value++obtuseTrianglePoints :: Vector.Vector Point+obtuseTrianglePoints =+ Vector.fromList+ [ Point (-2) 0+ , Point 2 0+ , Point 0 1+ ]++squarePoints :: Vector.Vector Point+squarePoints =+ Vector.fromList+ [ Point (-1) (-1)+ , Point 1 (-1)+ , Point 1 1+ , Point (-1) 1+ ] assertSelectedEdgeEndpoints :: DCELComplex -> Edge DCELComplex -> Assertion assertSelectedEdgeEndpoints complexValue edgeValue = do
test/native/Moonlight/Triangulation/NativeSpec.hs view
@@ -31,6 +31,7 @@ import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as MUV import Data.Word (Word32, Word64)+import Moonlight.Triangulation.Alpha import Moonlight.Triangulation.BulkLoad import Moonlight.Triangulation.Cdt import Moonlight.Triangulation.Dcel
test/serialization/Moonlight/Triangulation/SerializationSpec.hs view
@@ -10,13 +10,15 @@ import Control.DeepSeq (NFData) import Control.Monad (unless) import Data.Binary (Binary)-import Data.Binary.Put (putWord16be, putWord64be, runPut)+import Data.Binary.Put (putWord16be, putWord32be, putWord64be, putWord8, runPut) import qualified Data.ByteString.Lazy as BL import Data.Foldable (traverse_)+import Data.List.NonEmpty (NonEmpty (..)) import qualified Data.Vector as V import Data.Word (Word16, Word64) import GHC.Generics (Generic) import Moonlight.Triangulation+import Moonlight.Triangulation.Handles.HandleDefs (DirectedEdgeId (..), FaceId (..)) import Moonlight.Triangulation.Serialization import Moonlight.Triangulation.Types (KnownConstraintMode) import Support (assertEqual, assertValid, requireRight)@@ -28,7 +30,9 @@ testConstrainedRoundTrip testIndependentPayloadGeometryRoundTrip testPointPayloadRoundTrip+ testCanonicalizesSerializedSignedZero testRejectsHostileStructuralPrefixes+ testRejectsOutOfRangeOuterFaceReference testRejectsCorruption putStrLn "all serialization tests passed" @@ -131,6 +135,32 @@ (Point 13 17) (vertexData original vertex) assertSerializationRoundTrip "point payload serialization" original +-- The wire format admits every binary64 bit pattern, but resident point identity+-- admits only canonical zero. Mutate the first singleton coordinate to negative+-- zero and require decode to restore the canonical representation.+testCanonicalizesSerializedSignedZero :: IO ()+testCanonicalizesSerializedSignedZero = do+ original <- requireRight "signed-zero source" (delaunayGeometry (V.singleton (Point 0 0)))+ let bytes = encodeTriangulation original+ structuralPrefixSize = 8 + 2 + 1 + 1 + 4 * 8+ signedZeroBytes =+ BL.concat+ [ BL.take structuralPrefixSize bytes+ , BL.singleton 0x80+ , BL.drop (structuralPrefixSize + 1) bytes+ ]+ decoded <-+ requireRight+ "signed-zero decode"+ ( decodeTriangulation testDecodingBudget trustedBinaryPayloadDecoders signedZeroBytes+ :: Either SerializationError (Triangulation 'Unconstrained () () () ())+ )+ vertex <- case vertices decoded of+ [onlyVertex] -> pure onlyVertex+ unexpected -> fail ("signed-zero fixture produced " <> show (length unexpected) <> " vertices")+ assertEqual "signed-zero coordinate canonicalization" (Point 0 0) (vertexPoint decoded vertex)+ assertEqual "signed-zero canonical re-encoding" bytes (encodeTriangulation decoded)+ assertSerializationRoundTrip :: ( KnownConstraintMode mode , Binary vertex@@ -189,6 +219,19 @@ ) , ("version 5", testDecodingBudget, structuralPrefix 5 0 0 1 0, UnsupportedFormatVersion 5) ]++-- Range validation must obstruct before relational validation dereferences the+-- hostile handle. Fully comparing the error forces its complete violation+-- inventory and therefore pins the absence of a lazy unsafe-index crash.+testRejectsOutOfRangeOuterFaceReference :: IO ()+testRejectsOutOfRangeOuterFaceReference =+ assertDecodeFailure+ "out-of-range outer-face reference"+ testDecodingBudget+ (structuralPrefix 6 0 0 1 0 <> runPut (putWord32be 0 >> putWord8 0))+ ( DecodedInvariantViolations+ (FaceAdjacentOutOfRange (FaceId 0) (DirectedEdgeId 0) 0 :| [])+ ) assertDecodeFailure :: String