diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -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
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -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
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -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,
diff --git a/bench/alpha/Main.hs b/bench/alpha/Main.hs
new file mode 100644
--- /dev/null
+++ b/bench/alpha/Main.hs
@@ -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
diff --git a/bench/delaunay-compare/Moonlight/Triangulation/Bench/DelaunayCompare/Picture.hs b/bench/delaunay-compare/Moonlight/Triangulation/Bench/DelaunayCompare/Picture.hs
--- a/bench/delaunay-compare/Moonlight/Triangulation/Bench/DelaunayCompare/Picture.hs
+++ b/bench/delaunay-compare/Moonlight/Triangulation/Bench/DelaunayCompare/Picture.hs
@@ -293,10 +293,6 @@
, " 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)"
]
<> concatMap (renderSeries panelLeft panelWidth pointCounts axisMaximum) (picturePanelSeries panel)
- <> maybe
- []
- (renderAdvantageCallout panelLeft panelWidth pointCounts axisMaximum)
- (bestReportedSpreadAdvantage (picturePanelSeries panel))
<> [""]
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 =
- [ ""
- , " coordinate pointX <> "\" y1=\"" <> coordinate pointY <> "\" x2=\"" <> coordinate (labelX + 8) <> "\" y2=\"" <> coordinate (labelY + 18) <> "\" stroke=\"rgb(255,140,0)\" stroke-width=\"1.2\"/>"
- , " 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")
- , ""
- ]
- 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
diff --git a/bench/delaunay-compare/README.md b/bench/delaunay-compare/README.md
--- a/bench/delaunay-compare/README.md
+++ b/bench/delaunay-compare/README.md
@@ -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
diff --git a/bench/delaunay-compare/results/moonlight-delaunay-compare-small.svg b/bench/delaunay-compare/results/moonlight-delaunay-compare-small.svg
new file mode 100644
--- /dev/null
+++ b/bench/delaunay-compare/results/moonlight-delaunay-compare-small.svg
@@ -0,0 +1,415 @@
+
+
diff --git a/bench/support/BenchSupport.hs b/bench/support/BenchSupport.hs
--- a/bench/support/BenchSupport.hs
+++ b/bench/support/BenchSupport.hs
@@ -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
diff --git a/bindings/README.md b/bindings/README.md
--- a/bindings/README.md
+++ b/bindings/README.md
@@ -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.
diff --git a/bindings/python/src/moonlight_triangulation/__init__.py b/bindings/python/src/moonlight_triangulation/__init__.py
--- a/bindings/python/src/moonlight_triangulation/__init__.py
+++ b/bindings/python/src/moonlight_triangulation/__init__.py
@@ -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)
diff --git a/bindings/python/tests/test_binding.py b/bindings/python/tests/test_binding.py
--- a/bindings/python/tests/test_binding.py
+++ b/bindings/python/tests/test_binding.py
@@ -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__":
diff --git a/bindings/typescript/src/index.ts b/bindings/typescript/src/index.ts
--- a/bindings/typescript/src/index.ts
+++ b/bindings/typescript/src/index.ts
@@ -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(({ 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");
}
}
diff --git a/bindings/typescript/test/binding.test.ts b/bindings/typescript/test/binding.test.ts
--- a/bindings/typescript/test/binding.test.ts
+++ b/bindings/typescript/test/binding.test.ts
@@ -76,3 +76,28 @@
failure.inputIndex === 1n,
);
});
+
+test("native count descends from materialized coordinates", () => {
+ const engine = new Moonlight();
+ const points = new Proxy([[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();
+});
diff --git a/moonlight-triangulation.cabal b/moonlight-triangulation.cabal
--- a/moonlight-triangulation.cabal
+++ b/moonlight-triangulation.cabal
@@ -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
diff --git a/src-build/Moonlight/Triangulation/BulkLoad.hs b/src-build/Moonlight/Triangulation/BulkLoad.hs
--- a/src-build/Moonlight/Triangulation/BulkLoad.hs
+++ b/src-build/Moonlight/Triangulation/BulkLoad.hs
@@ -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
diff --git a/src-build/Moonlight/Triangulation/Internal/CircleSweep.hs b/src-build/Moonlight/Triangulation/Internal/CircleSweep.hs
--- a/src-build/Moonlight/Triangulation/Internal/CircleSweep.hs
+++ b/src-build/Moonlight/Triangulation/Internal/CircleSweep.hs
@@ -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)
diff --git a/src-cell-complex/Moonlight/Triangulation/CellComplex.hs b/src-cell-complex/Moonlight/Triangulation/CellComplex.hs
--- a/src-cell-complex/Moonlight/Triangulation/CellComplex.hs
+++ b/src-cell-complex/Moonlight/Triangulation/CellComplex.hs
@@ -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
diff --git a/src-core/Moonlight/Triangulation/Internal/Dyadic.hs b/src-core/Moonlight/Triangulation/Internal/Dyadic.hs
--- a/src-core/Moonlight/Triangulation/Internal/Dyadic.hs
+++ b/src-core/Moonlight/Triangulation/Internal/Dyadic.hs
@@ -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.
diff --git a/src-core/Moonlight/Triangulation/Internal/ExactRational.hs b/src-core/Moonlight/Triangulation/Internal/ExactRational.hs
--- a/src-core/Moonlight/Triangulation/Internal/ExactRational.hs
+++ b/src-core/Moonlight/Triangulation/Internal/ExactRational.hs
@@ -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
diff --git a/src-dcel/Moonlight/Triangulation/Alpha.hs b/src-dcel/Moonlight/Triangulation/Alpha.hs
new file mode 100644
--- /dev/null
+++ b/src-dcel/Moonlight/Triangulation/Alpha.hs
@@ -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
diff --git a/src-dcel/Moonlight/Triangulation/FloodFillIterator.hs b/src-dcel/Moonlight/Triangulation/FloodFillIterator.hs
--- a/src-dcel/Moonlight/Triangulation/FloodFillIterator.hs
+++ b/src-dcel/Moonlight/Triangulation/FloodFillIterator.hs
@@ -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,
diff --git a/src-dcel/Moonlight/Triangulation/Internal/CellSet.hs b/src-dcel/Moonlight/Triangulation/Internal/CellSet.hs
--- a/src-dcel/Moonlight/Triangulation/Internal/CellSet.hs
+++ b/src-dcel/Moonlight/Triangulation/Internal/CellSet.hs
@@ -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
diff --git a/src-dcel/Moonlight/Triangulation/Internal/Mutable.hs b/src-dcel/Moonlight/Triangulation/Internal/Mutable.hs
--- a/src-dcel/Moonlight/Triangulation/Internal/Mutable.hs
+++ b/src-dcel/Moonlight/Triangulation/Internal/Mutable.hs
@@ -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)
diff --git a/src-dcel/Moonlight/Triangulation/Internal/Types.hs b/src-dcel/Moonlight/Triangulation/Internal/Types.hs
--- a/src-dcel/Moonlight/Triangulation/Internal/Types.hs
+++ b/src-dcel/Moonlight/Triangulation/Internal/Types.hs
@@ -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
diff --git a/src-dcel/Moonlight/Triangulation/Types.hs b/src-dcel/Moonlight/Triangulation/Types.hs
--- a/src-dcel/Moonlight/Triangulation/Types.hs
+++ b/src-dcel/Moonlight/Triangulation/Types.hs
@@ -22,6 +22,10 @@
, CoordinateError (..)
, NonFiniteValue (..)
, classifyNonFinite
+ , RadiusSquared
+ , RadiusSquaredError (..)
+ , mkRadiusSquared
+ , radiusSquaredValue
, BuildError (..)
, Location (..)
, LocationHint (..)
diff --git a/src-dcel/Moonlight/Triangulation/Validation.hs b/src-dcel/Moonlight/Triangulation/Validation.hs
--- a/src-dcel/Moonlight/Triangulation/Validation.hs
+++ b/src-dcel/Moonlight/Triangulation/Validation.hs
@@ -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 -> []
diff --git a/src-public/Moonlight/Triangulation.hs b/src-public/Moonlight/Triangulation.hs
--- a/src-public/Moonlight/Triangulation.hs
+++ b/src-public/Moonlight/Triangulation.hs
@@ -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
diff --git a/src-serialize/Moonlight/Triangulation/Serialization.hs b/src-serialize/Moonlight/Triangulation/Serialization.hs
--- a/src-serialize/Moonlight/Triangulation/Serialization.hs
+++ b/src-serialize/Moonlight/Triangulation/Serialization.hs
@@ -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)
diff --git a/test/cell-complex/Main.hs b/test/cell-complex/Main.hs
--- a/test/cell-complex/Main.hs
+++ b/test/cell-complex/Main.hs
@@ -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
diff --git a/test/native/Moonlight/Triangulation/NativeSpec.hs b/test/native/Moonlight/Triangulation/NativeSpec.hs
--- a/test/native/Moonlight/Triangulation/NativeSpec.hs
+++ b/test/native/Moonlight/Triangulation/NativeSpec.hs
@@ -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
diff --git a/test/serialization/Moonlight/Triangulation/SerializationSpec.hs b/test/serialization/Moonlight/Triangulation/SerializationSpec.hs
--- a/test/serialization/Moonlight/Triangulation/SerializationSpec.hs
+++ b/test/serialization/Moonlight/Triangulation/SerializationSpec.hs
@@ -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