diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,25 @@ All notable changes to `moonlight-homology` are documented here. +## 0.1.0.2 - 2026-08-22 + +- Generalize `FilteredFiniteChainComplex` from the fixed binary64 + `FiltrationValue` to any ordered birth-key type, preserving exact geometric + filtration order through the canonical persistence reducer. Its constructor + is now sealed; callers retain read-only projections and construct only + through the coverage- and monotonicity-checking boundary. +- Add `persistentBettiAt`, the batched `persistentBettiAtMany`, and + `persistentBettiAtCriticalValues`. Filtered complexes cache exact critical + equivalence classes and dense complex-relative ranks, so all critical Betti + profiles are one ordered event sweep rather than one barcode scan per birth. +- Lower mod-two persistence through dense integer indices, `IntMap`/`IntSet` + columns, degree-wise clearing, and a sealed union-find specialization for + graph boundaries, with the general sparse reducer retained for incompatible + boundary fibers. +- Keep the checked finite-chain constructor as the sole public entrance while + deciding sparse boundary-composition zero directly, without materializing a + throwaway product matrix on the valid path. + ## 0.1.0.1 - 2026-08-21 - Rebuilt the multi-library Haddock archive with external Algebra, Core, and diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -28,8 +28,8 @@ backend is a compile error. - **Exact and spectral sequences.** Filtered spectral families with page-by-page reduction and convergence tracking; exact-sequence helpers; Block–Schur reductions. -- **Persistence.** One- and two-parameter filtered complexes and mod-2 persistence - pairs. +- **Persistence.** Arbitrary ordered one-parameter birth keys, mod-2 + persistence pairs and closed-sublevel Betti queries; two-parameter vocabulary. - **Discrete Morse theory.** Acyclic matchings that reduce a complex to its critical cells while preserving homology. - **Topological carriers.** Cell complexes, graph 1-skeletons, Reeb/macro-scaffold @@ -49,37 +49,30 @@ each oriented edge to `head − tail`; every degree above 1 is empty. ```haskell -{-# LANGUAGE DataKinds #-} - -import Moonlight.Homology +import Moonlight.Homology (FiniteChainComplex) +import Moonlight.Homology.Presentation -circle :: Either BoundaryIncidenceShapeError (FiniteChainComplex Rational) -circle = do - d1 <- - mkBoundaryIncidence 3 3 - [ mkBoundaryEntry 0 0 (-1), mkBoundaryEntry 0 1 1, - mkBoundaryEntry 1 1 (-1), mkBoundaryEntry 1 2 1, - mkBoundaryEntry 2 2 (-1), mkBoundaryEntry 2 0 1 - ] - pure $ - mkFiniteChainComplex (HomologicalDegree 1) $ \degree -> - case degree of - HomologicalDegree 1 -> d1 - HomologicalDegree 0 -> emptyBoundaryIncidenceOf 3 0 - _ -> emptyBoundaryIncidence +circle :: Either ChainBuildError (FiniteChainComplex Rational) +circle = + compileChain + ChainSpec + { chainCellCounts = [3, 3], + chainBoundaries = + [ [ (0, 0, -1), (0, 1, 1), + (1, 1, -1), (1, 2, 1), + (2, 2, -1), (2, 0, 1) + ] + ] + } ``` -`mkBoundaryIncidence sourceDim targetDim entries` builds the validated matrix of ∂ₙ, -the boundary map from the `sourceDim` cells of degree *n* to the `targetDim` cells of -degree *n − 1*. Each `mkBoundaryEntry source target coefficient` is one nonzero -incidence: the degree-*n* cell `source` contains the degree-*(n − 1)* cell `target` in -its boundary with that `coefficient`. Edge 0, entries `(0, 0, -1)` and +`chainCellCounts` lists cell counts from degree zero upward; +`chainBoundaries` supplies one sparse `(source, target, coefficient)` section for +each positive degree. Edge 0, entries `(0, 0, -1)` and `(0, 1, 1)`, encodes ∂(edge₀) = vertex₁ − vertex₀, running from vertex 0 (tail, `-1`) to vertex 1 (head, `+1`); edges 1 and 2 close the loop v₀ → v₁ → v₂ → v₀. Construction -is total: mismatched shapes are rejected as `BoundaryIncidenceShapeError`. `emptyBoundaryIncidenceOf sourceDim targetDim` is the zero map of that shape -(here ∂₀, three vertices to nothing), and `emptyBoundaryIncidence` the empty map used -above the top degree. `mkFiniteChainComplex topDegree atDegree` then assembles the -complex from its boundary at each degree: a top degree and one `∂` per degree. +is total: `compileChain` reports malformed incidence, shape mismatch, or failure of +∂ ∘ ∂ = 0 through `ChainBuildError`; the unchecked matrix constructor remains private. ### Betti numbers over a field @@ -109,20 +102,16 @@ with the 2-cell attached by a degree-2 map, giving H₁(RP²) = ℤ/2. ```haskell -{-# LANGUAGE DataKinds #-} - import Moonlight.Homology +import Moonlight.Homology.Presentation -realProjectivePlane :: Either BoundaryIncidenceShapeError (FiniteChainComplex Integer) -realProjectivePlane = do - d2 <- mkBoundaryIncidence 1 1 [mkBoundaryEntry 0 0 2] - pure $ - mkFiniteChainComplex (HomologicalDegree 2) $ \degree -> - case degree of - HomologicalDegree 2 -> d2 - HomologicalDegree 1 -> emptyBoundaryIncidenceOf 1 1 - HomologicalDegree 0 -> emptyBoundaryIncidenceOf 1 0 - _ -> emptyBoundaryIncidence +realProjectivePlane :: Either ChainBuildError (FiniteChainComplex Integer) +realProjectivePlane = + compileChain + ChainSpec + { chainCellCounts = [1, 1, 1], + chainBoundaries = [[], [(0, 0, 2)]] + } integralHomology :: FiniteChainComplex Integer -> Either HomologyFailure [HomologyGroup Integer] @@ -154,8 +143,15 @@ `Moonlight.Homology` umbrella, or from the narrower module noted below. - **Persistence.** `mkFilteredFiniteChainComplex` builds a filtered complex; - `mod2PersistentPairs` reads its birth/death pairs, and `BiPersistencePair` carries - the two-parameter case. In `Moonlight.Homology.Persistence`. + its birth key may be any ordered type, while `FiltrationValue` remains the + binary64 convenience specialization. `mod2PersistentPairs` reads exact + birth/death pairs, `persistentBettiAt` answers one closed-sublevel query, and + `persistentBettiAtMany` sweeps an arbitrary threshold family without + rescanning the barcode. For every admitted critical value, + `persistentBettiAtCriticalValues` uses the filtered complex's dense derived + ranks while retaining the exact births as the public authority. + `BiPersistencePair` carries the two-parameter case. In + `Moonlight.Homology.Persistence`. - **Spectral sequences.** `mkSpectralSource` and `spectralFamilyPages` produce the page-by-page family; `spectralFamilyLimitPage`, `spectralFamilyStableFrom`, and `convergenceDepth` track convergence. In `Moonlight.Homology.Sequence`. @@ -216,7 +212,7 @@ | `Moonlight.Homology.Backend` | The `HomologyBackend` dispatcher: Smith / rational / GF(2). | | `Moonlight.Homology.Sequence` | Exact and spectral sequences, Block–Schur reductions, graph spectral helpers. | | `Moonlight.Homology.Topology` | Cell complexes, graph skeletons, macro-scaffolds, discrete Morse, persistence values, observers, and constraints. | -| `Moonlight.Homology.Persistence` | Filtered complexes and mod-2 persistence pairs. | +| `Moonlight.Homology.Persistence` | Ordered filtered complexes, mod-2 persistence pairs, and barcode Betti queries. | | `Moonlight.Homology.Pure.Topology.CellComplex` | Generic two-dimensional cell incidence; requires `moonlight-homology:cell-complex`. | | `Moonlight.Homology.Pure.Topology.CellCategory` | Finite incidence category for a `CellComplex2D`; requires `moonlight-homology:cell-category`. | | `Moonlight.Homology.Effect.Laws` | Boundary-nilpotence and reduction law harnesses. | diff --git a/moonlight-homology.cabal b/moonlight-homology.cabal --- a/moonlight-homology.cabal +++ b/moonlight-homology.cabal @@ -1,6 +1,6 @@ cabal-version: 3.0 name: moonlight-homology -version: 0.1.0.1 +version: 0.1.0.2 homepage: https://github.com/PaleRoses/moonlight bug-reports: https://github.com/PaleRoses/moonlight/issues synopsis: Chain complexes, phase-gated homology interfaces, and spectral scaffolding. @@ -90,6 +90,7 @@ , moonlight-linalg:moonlight-linalg-dense >= 0.1 && < 0.2 , moonlight-pale:diagnostic >= 0.1 && < 0.2 , moonlight-homology:moonlight-homology-chain + , vector >= 0.13 && < 0.14 -- The incidence interface is independently useful to mesh and visualization -- consumers. Keep it outside the matrix/diagnostic topology closure. @@ -289,6 +290,7 @@ , tasty-hunit >= 0.10 && < 0.11 , tasty-quickcheck >= 0.10 && < 0.12 , text >= 2.0 && < 2.2 + , vector >= 0.13 && < 0.14 benchmark moonlight-homology-bench import: shared-properties @@ -315,5 +317,5 @@ source-repository this type: git location: https://github.com/PaleRoses/moonlight.git - tag: moonlight-homology-0.1.0.1 + tag: moonlight-homology-0.1.0.2 subdir: moonlight-homology diff --git a/src-matrix/Moonlight/Homology/Boundary/Finite.hs b/src-matrix/Moonlight/Homology/Boundary/Finite.hs --- a/src-matrix/Moonlight/Homology/Boundary/Finite.hs +++ b/src-matrix/Moonlight/Homology/Boundary/Finite.hs @@ -21,17 +21,17 @@ import qualified Data.Map.Strict as Map import Data.Set (Set) import qualified Data.Set as Set +import Data.Vector qualified as Vector import Moonlight.Core (Semiring) import Moonlight.Homology.Boundary.LinAlg ( BoundaryIncidence, boundaryCoefficient, - boundaryEntries, - composeBoundaryIncidence, + boundaryEntriesBySource, + boundaryIncidenceCompositionIsZero, emptyBoundaryIncidence, emptyBoundaryIncidenceOf, materializeIncidenceBoundary, sourceCardinality, - sourceIndex, targetCardinality, targetIndex, ) @@ -69,13 +69,13 @@ Int -> Either HomologyFailure () adjacentNilpotenceAt finite degreeIndex = - composeBoundaryIncidence + boundaryIncidenceCompositionIsZero (incidenceMatrixAt finite (HomologicalDegree degreeIndex)) (incidenceMatrixAt finite (HomologicalDegree (degreeIndex + 1))) & either (Left . InvalidBoundaryIncidence . show) - ( \composed -> - if null (boundaryEntries composed) + ( \compositionIsZero -> + if compositionIsZero then Right () else Left (ChainComplexNilpotenceViolation degreeIndex) ) @@ -184,15 +184,9 @@ ) | otherwise = let incidence = incidenceMatrixAt finiteComplex (HomologicalDegree dimensionValue) - entriesBySource = - foldl' - ( \accumulator entryValue -> - Map.insertWith (<>) (sourceIndex entryValue) [entryValue] accumulator - ) - Map.empty - (boundaryEntries incidence) + entriesBySource = boundaryEntriesBySource incidence restrictedEntriesOf sourceCell = - Map.findWithDefault [] (cellIndex sourceCell) entriesBySource + maybe [] id (entriesBySource Vector.!? cellIndex sourceCell) & mapMaybe ( \entryValue -> let targetCell = diff --git a/src-matrix/Moonlight/Homology/Boundary/LinAlg.hs b/src-matrix/Moonlight/Homology/Boundary/LinAlg.hs --- a/src-matrix/Moonlight/Homology/Boundary/LinAlg.hs +++ b/src-matrix/Moonlight/Homology/Boundary/LinAlg.hs @@ -12,8 +12,10 @@ sourceCardinality, targetCardinality, boundaryEntries, + boundaryEntriesBySource, mkBoundaryIncidence, mkBoundaryIncidenceFromOrderedEntries, + mkBoundaryIncidenceFromOrderedColumns, overlapBoundaryIncidence, emptyBoundaryIncidence, emptyBoundaryIncidenceOf, @@ -23,6 +25,7 @@ boundaryIncidenceApply, transposeBoundaryIncidence, composeBoundaryIncidence, + boundaryIncidenceCompositionIsZero, boundaryIncidenceDiagonal, addBoundaryIncidence, mapBoundaryCoefficients, @@ -43,10 +46,15 @@ ) where +import Control.Monad (foldM) +import Control.Monad.ST (ST, runST) import Data.Function ((&)) import Data.Kind (Type) +import Data.List qualified as List import qualified Data.Map.Strict as Map import Data.Maybe (listToMaybe, mapMaybe) +import Data.Vector qualified as Vector +import Data.Vector.Mutable qualified as MutableVector import Moonlight.Core (AdditiveMonoid (..), Semiring) import Moonlight.Homology.Pure.Failure (HomologyFailure (..), HomologyLaw (..)) import Numeric.Natural (Natural) @@ -93,10 +101,36 @@ data BoundaryIncidence r = BoundaryIncidence { sourceCardinality :: Int, targetCardinality :: Int, - boundaryEntries :: [BoundaryEntry r] + boundaryEntries :: [BoundaryEntry r], + boundaryEntriesBySourceCache :: Vector.Vector [BoundaryEntry r] } - deriving stock (Eq, Show) +instance Eq r => Eq (BoundaryIncidence r) where + left == right = + sourceCardinality left == sourceCardinality right + && targetCardinality left == targetCardinality right + && boundaryEntries left == boundaryEntries right + +instance Show r => Show (BoundaryIncidence r) where + showsPrec precedence incidence = + showParen + (precedence > 10) + ( showString "BoundaryIncidence {sourceCardinality = " + . shows (sourceCardinality incidence) + . showString ", targetCardinality = " + . shows (targetCardinality incidence) + . showString ", boundaryEntries = " + . shows (boundaryEntries incidence) + . showString "}" + ) + +-- | Canonical entries sectioned by source column. The projection is cached +-- lazily by the abstract incidence owner, so chain-law checking, restriction, +-- graph extraction, and persistence share one grouping without making it +-- semantic state. +boundaryEntriesBySource :: BoundaryIncidence r -> Vector.Vector [BoundaryEntry r] +boundaryEntriesBySource = boundaryEntriesBySourceCache + mkBoundaryIncidence :: (Eq r, Semiring r) => Natural -> Natural -> [BoundaryEntry r] -> Either BoundaryIncidenceShapeError (BoundaryIncidence r) mkBoundaryIncidence sourceCardinalityValue targetCardinalityValue entries = let sourceDimension = fromIntegral sourceCardinalityValue @@ -160,6 +194,10 @@ = OrderedEntriesOutOfOrder | OrderedCanonicalization !(Maybe (BoundaryEntry r)) ![BoundaryEntry r] +data OrderedColumnState + = OrderedColumnInvalid + | OrderedColumnPrefix !(Maybe Int) + canonicalizeOrderedEntries :: (Eq r, Semiring r) => [BoundaryEntry r] -> Maybe [BoundaryEntry r] canonicalizeOrderedEntries entries = finalizeOrderedCanonicalization @@ -218,7 +256,55 @@ (sourceIndex left, targetIndex left) (sourceIndex right, targetIndex right) +-- | Admit source-indexed canonical columns without flattening and regrouping +-- them. A noncanonical cover descends to the general ordered-entry constructor; +-- the authoritative 'BoundaryIncidence' and every typed shape obstruction stay +-- identical. +mkBoundaryIncidenceFromOrderedColumns :: + (Eq r, Semiring r) => + Natural -> + Natural -> + Vector.Vector [BoundaryEntry r] -> + Either BoundaryIncidenceShapeError (BoundaryIncidence r) +mkBoundaryIncidenceFromOrderedColumns sourceCardinalityValue targetCardinalityValue columns = + let sourceDimension = fromIntegral sourceCardinalityValue + targetDimension = fromIntegral targetCardinalityValue + flattenedEntries = concat (Vector.toList columns) + in if Vector.length columns == sourceDimension + && Vector.and (Vector.imap (columnIsCanonical targetDimension) columns) + then + Right + ( uncheckedBoundaryIncidenceFromColumns + sourceDimension + targetDimension + flattenedEntries + columns + ) + else + mkBoundaryIncidenceFromOrderedEntries + sourceCardinalityValue + targetCardinalityValue + flattenedEntries +columnIsCanonical :: (Eq r, Semiring r) => Int -> Int -> [BoundaryEntry r] -> Bool +columnIsCanonical targetDimension sourceIndexValue entries = + case List.foldl' advanceOrderedColumn (OrderedColumnPrefix Nothing) entries of + OrderedColumnInvalid -> False + OrderedColumnPrefix _ -> True + where + advanceOrderedColumn OrderedColumnInvalid _ = OrderedColumnInvalid + advanceOrderedColumn (OrderedColumnPrefix previousTarget) entry = + let currentTarget = targetIndex entry + targetFollows = maybe True (< currentTarget) previousTarget + in if sourceIndex entry == sourceIndexValue + && currentTarget >= 0 + && currentTarget < targetDimension + && boundaryCoefficient entry /= zero + && targetFollows + then OrderedColumnPrefix (Just currentTarget) + else OrderedColumnInvalid + + emptyBoundaryIncidence :: BoundaryIncidence r emptyBoundaryIncidence = uncheckedBoundaryIncidence 0 0 [] @@ -361,6 +447,100 @@ (fromIntegral (targetCardinality left)) composedEntries +-- | Decide whether a sparse composite is zero without materializing the +-- composite matrix. The local columns are the cover: each upper source is +-- reduced independently, and the result descends exactly when every overlap +-- coefficient cancels. +boundaryIncidenceCompositionIsZero :: + (Eq r, Num r, Semiring r) => + BoundaryIncidence r -> + BoundaryIncidence r -> + Either BoundaryIncidenceShapeError Bool +boundaryIncidenceCompositionIsZero left right + | targetCardinality right /= sourceCardinality left = + Left + ( BoundaryIncidenceShapeMismatch + (sourceCardinality left) + (targetCardinality left) + (sourceCardinality right) + (targetCardinality right) + ) + | Vector.all null (boundaryEntriesBySource left) + || Vector.all null (boundaryEntriesBySource right) = Right True + | otherwise = + let leftTermsBySource = boundaryEntriesBySource left + rightTermsBySource = boundaryEntriesBySource right + in Right + ( runST $ do + coefficientsByTarget <- + MutableVector.replicate (targetCardinality left) zero + Vector.foldM + (checkCompositeColumn coefficientsByTarget leftTermsBySource) + True + rightTermsBySource + ) + +checkCompositeColumn :: + (Eq r, Num r, Semiring r) => + MutableVector.MVector s r -> + Vector.Vector [BoundaryEntry r] -> + Bool -> + [BoundaryEntry r] -> + ST s Bool +checkCompositeColumn coefficientsByTarget leftTermsBySource precedingColumnsAreZero rightTerms = + if precedingColumnsAreZero + then do + touchedTargets <- + foldM + (accumulateRightTerm coefficientsByTarget leftTermsBySource) + [] + rightTerms + foldM + (inspectAndClearTarget coefficientsByTarget) + True + touchedTargets + else pure False + +accumulateRightTerm :: + Num r => + MutableVector.MVector s r -> + Vector.Vector [BoundaryEntry r] -> + [Int] -> + BoundaryEntry r -> + ST s [Int] +accumulateRightTerm coefficientsByTarget leftTermsBySource touchedTargets rightEntry = + foldM + (accumulateLeftTerm coefficientsByTarget (boundaryCoefficient rightEntry)) + touchedTargets + (Vector.unsafeIndex leftTermsBySource (targetIndex rightEntry)) + +accumulateLeftTerm :: + Num r => + MutableVector.MVector s r -> + r -> + [Int] -> + BoundaryEntry r -> + ST s [Int] +accumulateLeftTerm coefficientsByTarget rightCoefficient touchedTargets leftEntry = do + let targetIndexValue = targetIndex leftEntry + accumulatedCoefficient <- MutableVector.unsafeRead coefficientsByTarget targetIndexValue + MutableVector.unsafeWrite + coefficientsByTarget + targetIndexValue + (accumulatedCoefficient + boundaryCoefficient leftEntry * rightCoefficient) + pure (targetIndexValue : touchedTargets) + +inspectAndClearTarget :: + (Eq r, Semiring r) => + MutableVector.MVector s r -> + Bool -> + Int -> + ST s Bool +inspectAndClearTarget coefficientsByTarget precedingTargetsAreZero targetIndexValue = do + coefficientValue <- MutableVector.unsafeRead coefficientsByTarget targetIndexValue + MutableVector.unsafeWrite coefficientsByTarget targetIndexValue zero + pure (precedingTargetsAreZero && coefficientValue == zero) + boundaryIncidenceDiagonal :: Num r => BoundaryIncidence r -> Map.Map Int r boundaryIncidenceDiagonal incidence = boundaryEntries incidence @@ -551,10 +731,30 @@ uncheckedBoundaryIncidence :: Int -> Int -> [BoundaryEntry r] -> BoundaryIncidence r uncheckedBoundaryIncidence sourceDimension targetDimension entries = + uncheckedBoundaryIncidenceFromColumns + sourceDimension + targetDimension + entries + ( fmap reverse + ( Vector.accum + (flip (:)) + (Vector.replicate sourceDimension []) + (fmap (\entry -> (sourceIndex entry, entry)) entries) + ) + ) + +uncheckedBoundaryIncidenceFromColumns :: + Int -> + Int -> + [BoundaryEntry r] -> + Vector.Vector [BoundaryEntry r] -> + BoundaryIncidence r +uncheckedBoundaryIncidenceFromColumns sourceDimension targetDimension entries columns = BoundaryIncidence { sourceCardinality = sourceDimension, targetCardinality = targetDimension, - boundaryEntries = entries + boundaryEntries = entries, + boundaryEntriesBySourceCache = columns } firstOutOfBoundsEntry :: Int -> Int -> [BoundaryEntry r] -> Maybe (BoundaryEntry r) diff --git a/src-matrix/Moonlight/Homology/Pure/Rank/Field.hs b/src-matrix/Moonlight/Homology/Pure/Rank/Field.hs --- a/src-matrix/Moonlight/Homology/Pure/Rank/Field.hs +++ b/src-matrix/Moonlight/Homology/Pure/Rank/Field.hs @@ -14,7 +14,6 @@ import Data.IntMap.Strict (IntMap) import Data.IntMap.Strict qualified as IntMap import Data.Kind (Type) -import Data.Maybe (listToMaybe) import Moonlight.Core (Semiring, mkCapability) import Moonlight.Homology.Boundary.Finite ( FiniteChainComplex, @@ -26,7 +25,7 @@ import Moonlight.Homology.Boundary.LinAlg ( BoundaryIncidence, BoundaryIncidenceShapeError, - boundaryEntries, + boundaryIncidenceCompositionIsZero, composeBoundaryIncidence, ) import Moonlight.Homology.Pure.Matrix.Reducer @@ -217,19 +216,24 @@ HomologicalDegree -> Either (FieldHomologyFailure coeff) () validateNilpotenceAt finite degreeValue@(HomologicalDegree degreeInt) = - case - composeBoundaryIncidence - (incidenceMatrixAt finite (HomologicalDegree (degreeInt - 1))) - (incidenceMatrixAt finite degreeValue) - of - Left shapeError -> - Left (FieldHomologyBoundaryShapeFailed degreeValue shapeError) - Right composite -> - case listToMaybe (boundaryEntries composite) of - Nothing -> + let lowerBoundary = + incidenceMatrixAt finite (HomologicalDegree (degreeInt - 1)) + upperBoundary = incidenceMatrixAt finite degreeValue + in case + boundaryIncidenceCompositionIsZero + lowerBoundary + upperBoundary + of + Left shapeError -> + Left (FieldHomologyBoundaryShapeFailed degreeValue shapeError) + Right True -> Right () - Just _ -> - Left (FieldHomologyNonNilpotent degreeValue composite) + Right False -> + case composeBoundaryIncidence lowerBoundary upperBoundary of + Left shapeError -> + Left (FieldHomologyBoundaryShapeFailed degreeValue shapeError) + Right composite -> + Left (FieldHomologyNonNilpotent degreeValue composite) {-# INLINEABLE validateNilpotenceAt #-} rationalRank :: diff --git a/src-public/Moonlight/Homology.hs b/src-public/Moonlight/Homology.hs --- a/src-public/Moonlight/Homology.hs +++ b/src-public/Moonlight/Homology.hs @@ -36,6 +36,7 @@ boundaryEntries, mkBoundaryIncidence, mkBoundaryIncidenceFromOrderedEntries, + mkBoundaryIncidenceFromOrderedColumns, overlapBoundaryIncidence, emptyBoundaryIncidence, emptyBoundaryIncidenceOf, @@ -448,9 +449,21 @@ BiFilteredCell (..), BiPersistencePair (..), FiltrationValue (..), - FilteredFiniteChainComplex (..), + FilteredFiniteChainComplex, + filteredBaseComplex, + filteredCellBirths, + filteredCriticalValues, mkFilteredFiniteChainComplex, mod2PersistentPairs, + CriticalBettiTable, + criticalBettiDegreeCount, + criticalBettiRankCount, + criticalBettiTableValues, + criticalBettiVectors, + mod2PersistentPairsWithCriticalBettiTable, + persistentBettiAt, + persistentBettiAtMany, + persistentBettiAtCriticalValues, mod2PersistenceTopologyWitness, Bound (..), TargetBetti (..), diff --git a/src-public/Moonlight/Homology/Boundary.hs b/src-public/Moonlight/Homology/Boundary.hs --- a/src-public/Moonlight/Homology/Boundary.hs +++ b/src-public/Moonlight/Homology/Boundary.hs @@ -32,6 +32,7 @@ boundaryEntries, mkBoundaryIncidence, mkBoundaryIncidenceFromOrderedEntries, + mkBoundaryIncidenceFromOrderedColumns, overlapBoundaryIncidence, emptyBoundaryIncidence, emptyBoundaryIncidenceOf, diff --git a/src-public/Moonlight/Homology/Persistence.hs b/src-public/Moonlight/Homology/Persistence.hs --- a/src-public/Moonlight/Homology/Persistence.hs +++ b/src-public/Moonlight/Homology/Persistence.hs @@ -2,9 +2,21 @@ ( BiFilteredCell (..), BiPersistencePair (..), FiltrationValue (..), - FilteredFiniteChainComplex (..), + FilteredFiniteChainComplex, + filteredBaseComplex, + filteredCellBirths, + filteredCriticalValues, mkFilteredFiniteChainComplex, mod2PersistentPairs, + CriticalBettiTable, + criticalBettiDegreeCount, + criticalBettiRankCount, + criticalBettiTableValues, + criticalBettiVectors, + mod2PersistentPairsWithCriticalBettiTable, + persistentBettiAt, + persistentBettiAtMany, + persistentBettiAtCriticalValues, mod2PersistenceTopologyWitness, ) where @@ -13,10 +25,24 @@ import Moonlight.Homology.Pure.Chain (HomologicalDegree) import Moonlight.Homology.Topology (BasisCellRef) import Moonlight.Homology.Pure.Filtration (FiltrationValue (..)) -import Moonlight.Homology.Pure.Topology.Core (FilteredFiniteChainComplex (..)) +import Moonlight.Homology.Pure.Topology.Core + ( FilteredFiniteChainComplex, + filteredBaseComplex, + filteredCellBirths, + filteredCriticalValues, + ) import Moonlight.Homology.Pure.Topology.Persistence ( mkFilteredFiniteChainComplex, mod2PersistentPairs, + CriticalBettiTable, + criticalBettiDegreeCount, + criticalBettiRankCount, + criticalBettiTableValues, + criticalBettiVectors, + mod2PersistentPairsWithCriticalBettiTable, + persistentBettiAt, + persistentBettiAtMany, + persistentBettiAtCriticalValues, mod2PersistenceTopologyWitness, ) diff --git a/src-topology/Moonlight/Homology/Pure/Topology.hs b/src-topology/Moonlight/Homology/Pure/Topology.hs --- a/src-topology/Moonlight/Homology/Pure/Topology.hs +++ b/src-topology/Moonlight/Homology/Pure/Topology.hs @@ -126,6 +126,15 @@ representativeToVector, exactTopologyWitness, mod2PersistentPairs, + CriticalBettiTable, + criticalBettiDegreeCount, + criticalBettiRankCount, + criticalBettiTableValues, + criticalBettiVectors, + mod2PersistentPairsWithCriticalBettiTable, + persistentBettiAt, + persistentBettiAtMany, + persistentBettiAtCriticalValues, mod2PersistenceTopologyWitness, TopologyTarget (..), TargetViolation (..), @@ -218,7 +227,7 @@ type TopologyWitnessSeed :: Type -> Type data TopologyWitnessSeed r - = GraphTopologySeed Graph1Skeleton (Maybe (FilteredFiniteChainComplex r)) (Maybe ScalarPotentialField) Int + = GraphTopologySeed Graph1Skeleton (Maybe (FilteredFiniteChainComplex FiltrationValue r)) (Maybe ScalarPotentialField) Int | FiniteTopologySeed (FiniteChainComplex r) (TopologyObservationConfig r) observeTopologyWitnessSeed :: @@ -234,7 +243,7 @@ observeGraphTopologyWitness :: Integral r => - Maybe (FilteredFiniteChainComplex r) -> + Maybe (FilteredFiniteChainComplex FiltrationValue r) -> Maybe ScalarPotentialField -> Int -> Graph1Skeleton -> diff --git a/src-topology/Moonlight/Homology/Pure/Topology/BlockSchur.hs b/src-topology/Moonlight/Homology/Pure/Topology/BlockSchur.hs --- a/src-topology/Moonlight/Homology/Pure/Topology/BlockSchur.hs +++ b/src-topology/Moonlight/Homology/Pure/Topology/BlockSchur.hs @@ -17,13 +17,13 @@ import Data.Bifunctor (first) import Data.Foldable (traverse_) -import Data.IntMap.Strict qualified as IntMap import Data.Kind (Type) import Data.List (transpose) import Data.Map.Strict (Map) import Data.Map.Strict qualified as Map import Data.Maybe (mapMaybe) import Data.Set qualified as Set +import Data.Vector qualified as Vector import Moonlight.Core (Semiring) import Moonlight.Homology.Boundary.Finite ( FiniteChainComplex, @@ -38,6 +38,7 @@ BoundaryIncidenceShapeError (..), boundaryCoefficient, boundaryEntries, + boundaryEntriesBySource, composeBoundaryIncidence, emptyBoundaryIncidenceOf, mkBoundaryEntry, @@ -690,7 +691,7 @@ } boundaryOf :: - Map HomologicalDegree (IntMap.IntMap [BoundaryEntry coefficient]) -> + Map HomologicalDegree (Vector.Vector [BoundaryEntry coefficient]) -> BasisCellRef -> [(coefficient, BasisCellRef)] boundaryOf groupedBoundaryEntries basisRef = @@ -699,26 +700,17 @@ | degreeInt <= 0 -> [] | otherwise -> let targetDegree = HomologicalDegree (degreeInt - 1) - degreeEntries = Map.findWithDefault IntMap.empty (HomologicalDegree degreeInt) groupedBoundaryEntries + degreeEntries = Map.findWithDefault Vector.empty (HomologicalDegree degreeInt) groupedBoundaryEntries in [ (boundaryCoefficient entry, basisRefAt targetDegree (targetIndex entry)) - | entry <- IntMap.findWithDefault [] (cellIndex basisRef) degreeEntries + | entry <- maybe [] id (degreeEntries Vector.!? cellIndex basisRef) ] boundaryEntriesByDegree :: FiniteChainComplex coefficient -> - Map HomologicalDegree (IntMap.IntMap [BoundaryEntry coefficient]) + Map HomologicalDegree (Vector.Vector [BoundaryEntry coefficient]) boundaryEntriesByDegree complex = Map.fromList ( fmap - (\degreeValue -> (degreeValue, entriesBySource (incidenceMatrixAt complex degreeValue))) + (\degreeValue -> (degreeValue, boundaryEntriesBySource (incidenceMatrixAt complex degreeValue))) (dimensionsOf complex) ) - -entriesBySource :: BoundaryIncidence coefficient -> IntMap.IntMap [BoundaryEntry coefficient] -entriesBySource incidence = - foldr - ( \entryValue -> - IntMap.insertWith (<>) (sourceIndex entryValue) [entryValue] - ) - IntMap.empty - (boundaryEntries incidence) diff --git a/src-topology/Moonlight/Homology/Pure/Topology/Core.hs b/src-topology/Moonlight/Homology/Pure/Topology/Core.hs --- a/src-topology/Moonlight/Homology/Pure/Topology/Core.hs +++ b/src-topology/Moonlight/Homology/Pure/Topology/Core.hs @@ -24,8 +24,6 @@ enumerateFromZero, toRationalFromIntegral, mapMaybeWithLookup, - symmetricDifference, - lowIndex, ) where @@ -46,10 +44,25 @@ import Moonlight.Homology.Pure.Matrix.Shape (cellCountAtDegree, dimensionsOf) -type FilteredFiniteChainComplex :: Type -> Type -data FilteredFiniteChainComplex r = FilteredFiniteChainComplex +-- | A finite chain complex indexed by an arbitrary ordered filtration key. +-- The key is deliberately independent of the chain coefficient: exact +-- geometric filtrations must not be rounded through 'FiltrationValue' merely +-- to reach the persistence reducer. +type FilteredFiniteChainComplex :: Type -> Type -> Type +data FilteredFiniteChainComplex filtration r = FilteredFiniteChainComplex { filteredBaseComplex :: FiniteChainComplex r, - filteredCellBirths :: Map.Map BasisCellRef FiltrationValue + filteredCellBirths :: Map.Map BasisCellRef filtration, + -- | Exact critical-value equivalence classes in ascending order. + filteredCriticalValues :: ![filtration], + -- | Exact critical value to its dense, complex-relative rank. + filteredBirthValueRanks :: !(Map.Map filtration Int), + -- | Dense ranks local to this filtered complex. They are a derived + -- reduction index, never a semantic replacement for the exact birth. + filteredCellBirthRanks :: Map.Map BasisCellRef Int, + -- | Cells in the canonical @(birth rank, degree, basis index)@ reduction + -- order. Construction derives this once from the exact birth quotient; + -- reducers must not rediscover it through map lookups and a second sort. + filteredRankedCells :: ![(BasisCellRef, filtration, Int)] } type GraphEdge :: Type @@ -188,13 +201,3 @@ mapMaybeWithLookup :: Ord key => Map.Map key value -> [key] -> [value] mapMaybeWithLookup valueMap = mapMaybe (flip Map.lookup valueMap) - -symmetricDifference :: Ord a => Set.Set a -> Set.Set a -> Set.Set a -symmetricDifference left right = - (left `Set.difference` right) `Set.union` (right `Set.difference` left) - -lowIndex :: Set.Set Int -> Maybe Int -lowIndex rowSet = - case Set.maxView rowSet of - Nothing -> Nothing - Just (maximumIndexValue, _) -> Just maximumIndexValue diff --git a/src-topology/Moonlight/Homology/Pure/Topology/Graph/Skeleton.hs b/src-topology/Moonlight/Homology/Pure/Topology/Graph/Skeleton.hs --- a/src-topology/Moonlight/Homology/Pure/Topology/Graph/Skeleton.hs +++ b/src-topology/Moonlight/Homology/Pure/Topology/Graph/Skeleton.hs @@ -17,10 +17,10 @@ import Data.Function ((&)) import Data.Graph qualified as Graph -import Data.IntMap.Strict qualified as IntMap import Data.Map.Strict qualified as Map import Data.Set qualified as Set import Data.Tree qualified as Tree +import Data.Vector qualified as Vector import Moonlight.Core (Semiring) import Moonlight.Homology.Boundary.Finite ( FiniteChainComplex, @@ -32,6 +32,7 @@ BoundaryIncidence, boundaryCoefficient, boundaryEntries, + boundaryEntriesBySource, emptyBoundaryIncidence, emptyBoundaryIncidenceOf, materializeIncidenceBoundary, @@ -130,20 +131,11 @@ orientedUnitEdgeSupports :: Integral r => BoundaryIncidence r -> Either GraphSkeletonExtractionFailure [(Int, Int)] orientedUnitEdgeSupports edgeBoundary = enumerateFromZero (sourceCardinality edgeBoundary) - & traverse (unitEdgeSupport (entriesBySource edgeBoundary)) - -entriesBySource :: BoundaryIncidence r -> IntMap.IntMap [BoundaryEntry r] -entriesBySource incidence = - boundaryEntries incidence - & foldr - ( \entryValue -> - IntMap.insertWith (<>) (sourceIndex entryValue) [entryValue] - ) - IntMap.empty + & traverse (unitEdgeSupport (boundaryEntriesBySource edgeBoundary)) -unitEdgeSupport :: Integral r => IntMap.IntMap [BoundaryEntry r] -> Int -> Either GraphSkeletonExtractionFailure (Int, Int) +unitEdgeSupport :: Integral r => Vector.Vector [BoundaryEntry r] -> Int -> Either GraphSkeletonExtractionFailure (Int, Int) unitEdgeSupport groupedEntries edgeIndex = - case traverse signedUnitTarget (filter ((/= 0) . boundaryCoefficient) (IntMap.findWithDefault [] edgeIndex groupedEntries)) of + case traverse signedUnitTarget (filter ((/= 0) . boundaryCoefficient) (maybe [] id (groupedEntries Vector.!? edgeIndex))) of Just [(leftTarget, leftSign), (rightTarget, rightSign)] | leftTarget /= rightTarget && leftSign + rightSign == 0 -> Right (leftTarget, rightTarget) diff --git a/src-topology/Moonlight/Homology/Pure/Topology/Observation.hs b/src-topology/Moonlight/Homology/Pure/Topology/Observation.hs --- a/src-topology/Moonlight/Homology/Pure/Topology/Observation.hs +++ b/src-topology/Moonlight/Homology/Pure/Topology/Observation.hs @@ -5,12 +5,13 @@ where import Data.Kind (Type) +import Moonlight.Homology.Pure.Filtration (FiltrationValue) import Moonlight.Homology.Pure.Topology.Core (FilteredFiniteChainComplex) import Moonlight.Homology.Pure.Topology.MacroScaffold (ScalarPotentialField) type TopologyObservationConfig :: Type -> Type data TopologyObservationConfig r = TopologyObservationConfig - { observationFiltration :: Maybe (FilteredFiniteChainComplex r), + { observationFiltration :: Maybe (FilteredFiniteChainComplex FiltrationValue r), observationPotential :: Maybe ScalarPotentialField, observationLowModeCount :: Int } diff --git a/src-topology/Moonlight/Homology/Pure/Topology/Persistence.hs b/src-topology/Moonlight/Homology/Pure/Topology/Persistence.hs --- a/src-topology/Moonlight/Homology/Pure/Topology/Persistence.hs +++ b/src-topology/Moonlight/Homology/Pure/Topology/Persistence.hs @@ -1,373 +1,1128 @@ module Moonlight.Homology.Pure.Topology.Persistence ( mkFilteredFiniteChainComplex, mod2PersistentPairs, - mod2PersistenceTopologyWitness, - mod2PersistentBoundaryColumn, - persistentPairs, - persistenceTopologyWitness, - persistentBoundaryColumn, - persistenceEssentialBirths, - reducePersistentColumn, - reduceBoundaryColumn, - materializeFinitePersistencePair, - materializeEssentialPersistencePair, - orderedFilteredCells, - validateBirthCoverage, - validateBirthExactness, - validateFiltrationMonotonicity, - ) -where - -import Data.Function ((&)) -import Data.IntMap.Strict qualified as IntMap -import Data.Kind (Type) -import qualified Data.List as List -import qualified Data.Map.Strict as Map -import Data.Maybe (mapMaybe) -import qualified Data.Set as Set -import Moonlight.Homology.Boundary.Finite - ( FiniteChainComplex, - incidenceMatrixAt, - ) -import Moonlight.Homology.Boundary.LinAlg - ( BoundaryEntry, - BoundaryIncidence, - boundaryCoefficient, - boundaryEntries, - sourceIndex, - targetIndex, - ) -import Moonlight.Homology.Pure.Chain - ( HomologicalDegree (..), - PersistencePair (..), - TopologyWitness (..), - decrementDegree, - emptyTopologyWitness, - ) -import Moonlight.Homology.Pure.Failure (HomologyFailure (..)) -import Moonlight.Homology.Pure.Topology.Core - -type OrderedFilteredCell :: Type -data OrderedFilteredCell = OrderedFilteredCell - { orderedCellIdentity :: BasisCellRef, - orderedCellBirth :: FiltrationValue - } - deriving stock (Eq, Show) - -type PersistenceState :: Type -data PersistenceState = PersistenceState - { persistenceLowColumns :: !(Map.Map Int (Set.Set Int)), - persistencePairsByIndex :: ![(Int, Int)], - persistenceCreators :: !(Set.Set Int) - } - -emptyPersistenceState :: PersistenceState -emptyPersistenceState = - PersistenceState - { persistenceLowColumns = Map.empty, - persistencePairsByIndex = [], - persistenceCreators = Set.empty - } - -mkFilteredFiniteChainComplex :: - Integral r => - FiniteChainComplex r -> - [(BasisCellRef, FiltrationValue)] -> - Either HomologyFailure (FilteredFiniteChainComplex r) -mkFilteredFiniteChainComplex finite births = do - let birthMap = Map.fromList births - validateBirthUniqueness births birthMap - validateBirthCoverage finite birthMap - validateBirthExactness finite birthMap - validateFiltrationMonotonicity finite birthMap - pure - FilteredFiniteChainComplex - { filteredBaseComplex = finite, - filteredCellBirths = birthMap - } - -mod2PersistentPairs :: - Integral r => - FilteredFiniteChainComplex r -> - Either HomologyFailure [PersistencePair FiltrationValue] -mod2PersistentPairs filtered = do - let orderedCells = orderedFilteredCells filtered - orderedCellByIndex = - orderedCells - & zip [0 :: Int ..] - & Map.fromList - globalIndexByCell = - orderedCells - & zip [0 :: Int ..] - & fmap (\(globalIndexValue, orderedCell) -> (orderedCellIdentity orderedCell, globalIndexValue)) - & Map.fromList - groupedBoundaryEntries = boundaryEntriesByDegree (filteredBaseComplex filtered) - boundaryColumns = fmap (mod2PersistentBoundaryColumnFromIndex groupedBoundaryEntries globalIndexByCell) orderedCells - stateAfterReduction = foldl' reducePersistentColumn emptyPersistenceState (zip [0 :: Int ..] boundaryColumns) - finitePairs = - persistencePairsByIndex stateAfterReduction - & reverse - & mapMaybe (uncurry (materializeFinitePersistencePair orderedCellByIndex)) - essentialPairs = - persistenceEssentialBirths stateAfterReduction - & Set.toAscList - & mapMaybe (materializeEssentialPersistencePair orderedCellByIndex) - pure (finitePairs <> essentialPairs) - -persistentPairs :: - Integral r => - FilteredFiniteChainComplex r -> - Either HomologyFailure [PersistencePair FiltrationValue] -{-# DEPRECATED persistentPairs "Use mod2PersistentPairs — this computes mod-2 persistence despite its Integral constraint" #-} -persistentPairs = mod2PersistentPairs - -mod2PersistenceTopologyWitness :: - Integral r => - FilteredFiniteChainComplex r -> - Either HomologyFailure (TopologyWitness scaffold spectral FiltrationValue coefficient basis) -mod2PersistenceTopologyWitness filtered = do - pairs <- mod2PersistentPairs filtered - pure - emptyTopologyWitness - { topologyPersistencePairs = pairs - } - -persistenceTopologyWitness :: - Integral r => - FilteredFiniteChainComplex r -> - Either HomologyFailure (TopologyWitness scaffold spectral FiltrationValue coefficient basis) -{-# DEPRECATED persistenceTopologyWitness "Use mod2PersistenceTopologyWitness — this computes mod-2 persistence despite its Integral constraint" #-} -persistenceTopologyWitness = mod2PersistenceTopologyWitness - -persistenceEssentialBirths :: PersistenceState -> Set.Set Int -persistenceEssentialBirths stateValue = - let pairedBirths = persistencePairsByIndex stateValue & fmap fst & Set.fromList - in persistenceCreators stateValue `Set.difference` pairedBirths - -reducePersistentColumn :: PersistenceState -> (Int, Set.Set Int) -> PersistenceState -reducePersistentColumn stateValue (columnIndexValue, initialColumn) = - let reducedColumn = reduceBoundaryColumn (persistenceLowColumns stateValue) initialColumn - in case lowIndex reducedColumn of - Nothing -> - stateValue - { persistenceCreators = Set.insert columnIndexValue (persistenceCreators stateValue) - } - Just lowValue -> - stateValue - { persistenceLowColumns = Map.insert lowValue reducedColumn (persistenceLowColumns stateValue), - persistencePairsByIndex = (lowValue, columnIndexValue) : persistencePairsByIndex stateValue - } - -reduceBoundaryColumn :: Map.Map Int (Set.Set Int) -> Set.Set Int -> Set.Set Int -reduceBoundaryColumn lowColumns columnValue = - case lowIndex columnValue >>= (`Map.lookup` lowColumns) of - Nothing -> columnValue - Just pivotColumn -> reduceBoundaryColumn lowColumns (symmetricDifference columnValue pivotColumn) - -materializeFinitePersistencePair :: - Map.Map Int OrderedFilteredCell -> - Int -> - Int -> - Maybe (PersistencePair FiltrationValue) -materializeFinitePersistencePair orderedCellByIndex birthIndexValue deathIndexValue = - case (Map.lookup birthIndexValue orderedCellByIndex, Map.lookup deathIndexValue orderedCellByIndex) of - (Just birthCell, Just deathCell) -> - Just - PersistencePair - { persistenceDegree = cellDegree (orderedCellIdentity birthCell), - persistenceBirth = orderedCellBirth birthCell, - persistenceDeath = Just (orderedCellBirth deathCell) - } - _ -> Nothing - -materializeEssentialPersistencePair :: - Map.Map Int OrderedFilteredCell -> - Int -> - Maybe (PersistencePair FiltrationValue) -materializeEssentialPersistencePair orderedCellByIndex birthIndexValue = - Map.lookup birthIndexValue orderedCellByIndex - & fmap - ( \birthCell -> - PersistencePair - { persistenceDegree = cellDegree (orderedCellIdentity birthCell), - persistenceBirth = orderedCellBirth birthCell, - persistenceDeath = Nothing - } - ) - -orderedFilteredCells :: FilteredFiniteChainComplex r -> [OrderedFilteredCell] -orderedFilteredCells filtered = - allBasisCellRefs (filteredBaseComplex filtered) - & mapMaybe - ( \cellRefValue -> - fmap (OrderedFilteredCell cellRefValue) - (Map.lookup cellRefValue (filteredCellBirths filtered)) - ) - & List.sortOn - ( \orderedCell -> - let degreeValue = cellDegree (orderedCellIdentity orderedCell) - in ( orderedCellBirth orderedCell, - unHomologicalDegree degreeValue, - cellIndex (orderedCellIdentity orderedCell) - ) - ) - -mod2PersistentBoundaryColumn :: - Integral r => - FilteredFiniteChainComplex r -> - Map.Map BasisCellRef Int -> - OrderedFilteredCell -> - Set.Set Int -mod2PersistentBoundaryColumn filtered globalIndexByCell orderedCell = - let cellRefValue = orderedCellIdentity orderedCell - degreeValue = cellDegree cellRefValue - incidence = incidenceMatrixAt (filteredBaseComplex filtered) degreeValue - in mod2PersistentBoundaryColumnFromEntries (entriesBySource incidence) globalIndexByCell orderedCell - -mod2PersistentBoundaryColumnFromIndex :: - Integral r => - Map.Map HomologicalDegree (IntMap.IntMap [BoundaryEntry r]) -> - Map.Map BasisCellRef Int -> - OrderedFilteredCell -> - Set.Set Int -mod2PersistentBoundaryColumnFromIndex groupedBoundaryEntries globalIndexByCell orderedCell = - let cellRefValue = orderedCellIdentity orderedCell - degreeEntries = Map.findWithDefault IntMap.empty (cellDegree cellRefValue) groupedBoundaryEntries - in mod2PersistentBoundaryColumnFromEntries degreeEntries globalIndexByCell orderedCell - -mod2PersistentBoundaryColumnFromEntries :: - Integral r => - IntMap.IntMap [BoundaryEntry r] -> - Map.Map BasisCellRef Int -> - OrderedFilteredCell -> - Set.Set Int -mod2PersistentBoundaryColumnFromEntries groupedEntries globalIndexByCell orderedCell = - let cellRefValue = orderedCellIdentity orderedCell - degreeValue = cellDegree cellRefValue - in IntMap.findWithDefault [] (cellIndex cellRefValue) groupedEntries - & filter (\entry -> odd (abs (boundaryCoefficient entry))) - & fmap - ( \entry -> - BasisCellRef - { cellDegree = decrementDegree degreeValue, - cellIndex = targetIndex entry - } - ) - & mapMaybeWithLookup globalIndexByCell - & Set.fromList - -boundaryEntriesByDegree :: - FiniteChainComplex r -> - Map.Map HomologicalDegree (IntMap.IntMap [BoundaryEntry r]) -boundaryEntriesByDegree finite = - dimensionsOf finite - & fmap (\degreeValue -> (degreeValue, entriesBySource (incidenceMatrixAt finite degreeValue))) - & Map.fromList - -entriesBySource :: BoundaryIncidence r -> IntMap.IntMap [BoundaryEntry r] -entriesBySource incidence = - boundaryEntries incidence - & foldr - ( \entryValue -> - IntMap.insertWith (<>) (sourceIndex entryValue) [entryValue] - ) - IntMap.empty - -persistentBoundaryColumn :: - Integral r => - FilteredFiniteChainComplex r -> - Map.Map BasisCellRef Int -> - OrderedFilteredCell -> - Set.Set Int -{-# DEPRECATED persistentBoundaryColumn "Use mod2PersistentBoundaryColumn — this computes mod-2 boundary despite its Integral constraint" #-} -persistentBoundaryColumn = mod2PersistentBoundaryColumn - -validateBirthUniqueness :: - [(BasisCellRef, FiltrationValue)] -> - Map.Map BasisCellRef FiltrationValue -> - Either HomologyFailure () -validateBirthUniqueness births birthMap = - if length births == Map.size birthMap - then Right () - else Left (InvalidTopologyInput "duplicate birth assignments for the same cell") - -validateBirthCoverage :: - FiniteChainComplex r -> - Map.Map BasisCellRef FiltrationValue -> - Either HomologyFailure () -validateBirthCoverage finite birthMap = - allBasisCellRefs finite - & List.find (\cellRefValue -> Map.notMember cellRefValue birthMap) - & maybe (Right ()) missingCellFailure - where - missingCellFailure :: BasisCellRef -> Either HomologyFailure () - missingCellFailure cellRefValue = - Left - ( InvalidTopologyInput - ( "missing filtration value for cell " - <> show cellRefValue - ) - ) - -validateBirthExactness :: - FiniteChainComplex r -> - Map.Map BasisCellRef FiltrationValue -> - Either HomologyFailure () -validateBirthExactness finite birthMap = - let basisSet = Set.fromList (allBasisCellRefs finite) - extraKeys = Map.keysSet birthMap `Set.difference` basisSet - in if Set.null extraKeys - then Right () - else - Left - ( InvalidTopologyInput - ( "birth map contains cells absent from the chain complex: " - <> show (Set.toList extraKeys) - ) - ) - -validateFiltrationMonotonicity :: - Integral r => - FiniteChainComplex r -> - Map.Map BasisCellRef FiltrationValue -> - Either HomologyFailure () -validateFiltrationMonotonicity finite birthMap = - filtrationViolations - & List.find (const True) - & maybe (Right ()) (Left . InvalidTopologyInput) - where - filtrationViolations = - dimensionsOf finite - >>= ( \degreeValue@(HomologicalDegree degreeIndex) -> - if degreeIndex <= 0 - then [] - else - let incidence = incidenceMatrixAt finite degreeValue - in boundaryEntries incidence - & filter (\entry -> boundaryCoefficient entry /= 0) - & mapMaybe - ( \entry -> - let sourceCell = - BasisCellRef - { cellDegree = degreeValue, - cellIndex = sourceIndex entry - } - targetCell = - BasisCellRef - { cellDegree = decrementDegree degreeValue, - cellIndex = targetIndex entry - } - in case (Map.lookup sourceCell birthMap, Map.lookup targetCell birthMap) of - (Just sourceBirth, Just targetBirth) -> - if targetBirth <= sourceBirth - then Nothing - else - Just - ( "filtration violates face monotonicity for " - <> show sourceCell - <> " -> " - <> show targetCell - ) - _ -> Nothing - ) - ) + CriticalBettiTable, + criticalBettiDegreeCount, + criticalBettiRankCount, + criticalBettiTableValues, + criticalBettiVectors, + mod2PersistentPairsWithCriticalBettiTable, + persistentBettiAt, + persistentBettiAtMany, + persistentBettiAtCriticalValues, + mod2PersistenceTopologyWitness, + mod2PersistentBoundaryColumn, + persistentPairs, + persistenceTopologyWitness, + persistentBoundaryColumn, + persistenceEssentialBirths, + materializeFinitePersistencePair, + materializeEssentialPersistencePair, + orderedFilteredCells, + validateBirthCoverage, + validateBirthExactness, + validateFiltrationMonotonicity, + ) +where + +import Control.Monad (foldM) +import Control.Monad.ST (ST, runST) +import Data.Foldable (traverse_) +import Data.Function ((&)) +import Data.IntMap.Strict qualified as IntMap +import Data.IntSet qualified as IntSet +import Data.Kind (Type) +import qualified Data.List as List +import qualified Data.Map.Strict as Map +import Data.Maybe (mapMaybe) +import qualified Data.Set as Set +import Data.Vector qualified as Vector +import Data.Vector.Mutable qualified as BoxedMutableVector +import Data.Vector.Unboxed qualified as UnboxedVector +import Data.Vector.Unboxed.Mutable qualified as MutableVector +import Moonlight.Homology.Boundary.Finite + ( FiniteChainComplex, + incidenceMatrixAt, + ) +import Moonlight.Homology.Boundary.LinAlg + ( BoundaryEntry, + BoundaryIncidence, + boundaryCoefficient, + boundaryEntriesBySource, + sourceIndex, + targetIndex, + ) +import Moonlight.Homology.Pure.Chain + ( HomologicalDegree (..), + PersistencePair (..), + TopologyWitness (..), + decrementDegree, + emptyTopologyWitness, + ) +import Moonlight.Homology.Pure.Failure (HomologyFailure (..)) +import Moonlight.Homology.Pure.Topology.Core + +type OrderedFilteredCell :: Type -> Type +data OrderedFilteredCell filtration = OrderedFilteredCell + { orderedCellIdentity :: BasisCellRef, + orderedCellBirth :: filtration, + orderedCellBirthRank :: Int + } + deriving stock (Eq, Show) + +-- | Dense Betti profiles over the finite complex's critical-rank and degree +-- cover. Values are stored rank-major in one unboxed vector; the constructor +-- is withheld so the two dimensions cannot disagree with its storage. +data CriticalBettiTable = CriticalBettiTable + { criticalBettiDegreeCount :: !Int, + criticalBettiRankCount :: !Int, + criticalBettiTableValues :: !(UnboxedVector.Vector Int) + } + deriving stock (Eq, Show) + +-- | Materialize the familiar rank-indexed Betti vectors from the dense table. +criticalBettiVectors :: CriticalBettiTable -> [[Int]] +criticalBettiVectors table = + fmap profileAt [0 .. criticalBettiRankCount table - 1] + where + degreeCount = criticalBettiDegreeCount table + values = criticalBettiTableValues table + profileAt rankValue = + values + & UnboxedVector.drop (rankValue * degreeCount) + & UnboxedVector.take degreeCount + & UnboxedVector.toList + +type PersistenceState :: Type +data PersistenceState = PersistenceState + { persistencePairsByIndex :: ![(Int, Int)], + persistencePairedBirths :: !IntSet.IntSet, + persistenceCreators :: !IntSet.IntSet + } + +type BettiSweepPoint :: Type +data BettiSweepPoint = BettiSweepPoint + { bettiSweepDelta :: !(Map.Map HomologicalDegree Int), + bettiSweepRequested :: !Bool + } + +emptyPersistenceState :: PersistenceState +emptyPersistenceState = + PersistenceState + { persistencePairsByIndex = [], + persistencePairedBirths = IntSet.empty, + persistenceCreators = IntSet.empty + } + +mkFilteredFiniteChainComplex :: + (Integral r, Ord filtration) => + FiniteChainComplex r -> + [(BasisCellRef, filtration)] -> + Either HomologyFailure (FilteredFiniteChainComplex filtration r) +mkFilteredFiniteChainComplex finite births = do + let birthMap = Map.fromList births + validateBirthUniqueness births birthMap + denseBirthsByDegree <- denseBirthSection finite birthMap + let expectedCellCount = Vector.sum (fmap Vector.length denseBirthsByDegree) + if Map.size birthMap == expectedCellCount + then Right () + else validateBirthExactness finite birthMap + let birthSections = + Map.foldlWithKey' + (\sections cellRefValue birthValue -> Map.insertWith (<>) birthValue [cellRefValue] sections) + Map.empty + birthMap + & Map.toAscList + rankedBirthSections = zip [0 :: Int ..] birthSections + criticalValues = fmap (fst . snd) rankedBirthSections + birthValueRanks = + rankedBirthSections + & fmap (\(rankValue, (birthValue, _)) -> (birthValue, rankValue)) + & Map.fromDistinctAscList + cellBirthRanks = + rankedBirthSections + >>= ( \(rankValue, (_, cellsAtBirth)) -> + fmap (\cellRefValue -> (cellRefValue, rankValue)) cellsAtBirth + ) + & Map.fromList + rankedCells = + rankedBirthSections + >>= ( \(rankValue, (birthValue, cellsAtBirth)) -> + reverse cellsAtBirth + & fmap (\cellRefValue -> (cellRefValue, birthValue, rankValue)) + ) + validateDenseFiltrationMonotonicity finite denseBirthsByDegree + pure + FilteredFiniteChainComplex + { filteredBaseComplex = finite, + filteredCellBirths = birthMap, + filteredCriticalValues = criticalValues, + filteredBirthValueRanks = birthValueRanks, + filteredCellBirthRanks = cellBirthRanks, + filteredRankedCells = rankedCells + } + +mod2PersistentPairs :: + Integral r => + FilteredFiniteChainComplex filtration r -> + Either HomologyFailure [PersistencePair filtration] +mod2PersistentPairs filtered = + let (orderedCellByIndex, stateAfterReduction) = reduceMod2Persistence filtered + in Right (materializePersistencePairs orderedCellByIndex stateAfterReduction) + +-- | Reduce once, then derive both the exact barcode and all critical-rank +-- Betti profiles from the same indexed witness. The fused observation avoids +-- throwing away proven cell ranks only to recover them through exact-value +-- map lookups in 'persistentBettiAtCriticalValues'. +mod2PersistentPairsWithCriticalBettiTable :: + Integral r => + FilteredFiniteChainComplex filtration r -> + Either + HomologyFailure + ( [PersistencePair filtration] + , CriticalBettiTable + ) +mod2PersistentPairsWithCriticalBettiTable filtered = + let (orderedCellByIndex, stateAfterReduction) = reduceMod2Persistence filtered + pairs = materializePersistencePairs orderedCellByIndex stateAfterReduction + criticalBettiTable = + criticalBettiTableFromRankedEvents + (length (dimensionsOf (filteredBaseComplex filtered))) + (length (filteredCriticalValues filtered)) + (persistenceStateRankedEvents orderedCellByIndex stateAfterReduction) + in Right (pairs, criticalBettiTable) + +reduceMod2Persistence :: + Integral r => + FilteredFiniteChainComplex filtration r -> + (Vector.Vector (OrderedFilteredCell filtration), PersistenceState) +reduceMod2Persistence filtered = + let orderedCells = orderedFilteredCells filtered + orderedCellByIndex = + orderedCells + & Vector.fromList + globalIndicesByDegree = + denseGlobalIndicesByDegree + (filteredBaseComplex filtered) + (zip [0 :: Int ..] orderedCells) + groupedBoundaryEntries = boundaryEntriesByDegree (filteredBaseComplex filtered) + degreeCount = length (dimensionsOf (filteredBaseComplex filtered)) + indexedCellsByDegree = + zip [0 :: Int ..] orderedCells + & mapMaybe (indexedCellDegreeUpdate degreeCount) + & Vector.accum + (flip (:)) + (Vector.replicate degreeCount []) + & fmap reverse + higherDegreeColumns = + indexedCellsByDegree + & Vector.drop 2 + & fmap + ( fmap + ( \(columnIndexValue, orderedCell) -> + ( columnIndexValue + , mod2PersistentBoundaryColumnFromIndex + groupedBoundaryEntries + globalIndicesByDegree + orderedCell + ) + ) + ) + & Vector.toList + & reverse + degreeOneColumns = + maybe [] id (indexedCellsByDegree Vector.!? 1) + & fmap + ( \(columnIndexValue, orderedCell) -> + ( columnIndexValue + , mod2PersistentBoundaryTargetsFromIndex + groupedBoundaryEntries + globalIndicesByDegree + orderedCell + ) + ) + degreeZeroIndices = + maybe [] id (indexedCellsByDegree Vector.!? 0) + & fmap fst + stateAfterReduction = + reducePersistentColumnsByDegree + (length orderedCells) + higherDegreeColumns + degreeOneColumns + degreeZeroIndices + in (orderedCellByIndex, stateAfterReduction) + +indexedCellDegreeUpdate :: + Int -> + (Int, OrderedFilteredCell filtration) -> + Maybe (Int, (Int, OrderedFilteredCell filtration)) +indexedCellDegreeUpdate degreeCount indexedCell@(_, orderedCell) = + let degreeIndex = + unHomologicalDegree (cellDegree (orderedCellIdentity orderedCell)) + in if degreeIndex < 0 || degreeIndex >= degreeCount + then Nothing + else Just (degreeIndex, indexedCell) + +materializePersistencePairs :: + Vector.Vector (OrderedFilteredCell filtration) -> + PersistenceState -> + [PersistencePair filtration] +materializePersistencePairs orderedCellByIndex stateAfterReduction = + let finitePairs = + persistencePairsByIndex stateAfterReduction + & List.sortOn snd + & mapMaybe (uncurry (materializeFinitePersistencePair orderedCellByIndex)) + essentialPairs = + persistenceEssentialBirths stateAfterReduction + & IntSet.toAscList + & mapMaybe (materializeEssentialPersistencePair orderedCellByIndex) + in finitePairs <> essentialPairs + +persistenceStateRankedEvents :: + Vector.Vector (OrderedFilteredCell filtration) -> + PersistenceState -> + [(Int, (Int, Int))] +persistenceStateRankedEvents orderedCellByIndex stateAfterReduction = + let finiteEvents = + persistencePairsByIndex stateAfterReduction + >>= ( \(birthIndexValue, deathIndexValue) -> + case + ( orderedCellByIndex Vector.!? birthIndexValue + , orderedCellByIndex Vector.!? deathIndexValue + ) + of + (Just birthCell, Just deathCell) -> + rankedIntervalEvents + (cellDegree (orderedCellIdentity birthCell)) + (orderedCellBirthRank birthCell) + (Just (orderedCellBirthRank deathCell)) + _ -> [] + ) + essentialEvents = + persistenceEssentialBirths stateAfterReduction + & IntSet.toAscList + >>= ( \birthIndexValue -> + maybe + [] + ( \birthCell -> + rankedIntervalEvents + (cellDegree (orderedCellIdentity birthCell)) + (orderedCellBirthRank birthCell) + Nothing + ) + (orderedCellByIndex Vector.!? birthIndexValue) + ) + in finiteEvents <> essentialEvents + +-- | Read the Betti numbers of one closed filtration sublevel from an already +-- reduced barcode. The returned map is sparse: absent degrees have Betti +-- number zero. Finite persistence intervals are half-open, so a class whose +-- death equals the requested threshold is no longer present after every cell +-- at that threshold has entered the complex. +persistentBettiAt :: + Ord filtration => + filtration -> + [PersistencePair filtration] -> + Map.Map HomologicalDegree Int +persistentBettiAt threshold = + Map.fromListWith (+) + . fmap (\pairValue -> (persistenceDegree pairValue, 1)) + . filter (persistencePairAliveAt threshold) + +-- | Read several closed filtration sublevels in one ordered sweep over the +-- barcode. Results retain the supplied threshold order and multiplicity. +-- Unlike mapping 'persistentBettiAt', this does not rescan every persistence +-- interval for every requested threshold. +-- +-- Finite intervals use the same half-open semantics as 'persistentBettiAt'. A +-- birth and death at the same threshold cancel before that threshold is +-- observed. Malformed finite intervals whose death does not follow their birth +-- contribute nowhere, matching the pointwise query. +persistentBettiAtMany :: + Ord filtration => + [filtration] -> + [PersistencePair filtration] -> + [Map.Map HomologicalDegree Int] +persistentBettiAtMany thresholds pairs = + let requestedPoints = + thresholds + & fmap (\threshold -> (threshold, requestedBettiSweepPoint)) + & Map.fromListWith mergeBettiSweepPoints + eventPoints = + foldl' + insertPersistencePairSweepPoints + Map.empty + pairs + sweepPoints = Map.unionWith mergeBettiSweepPoints eventPoints requestedPoints + observedByThreshold = + sweepPoints + & Map.mapAccumWithKey applyBettiSweepPoint Map.empty + & snd + & Map.mapMaybe id + in fmap + (\threshold -> Map.findWithDefault Map.empty threshold observedByThreshold) + thresholds + +-- | Read every exact critical sublevel of one admitted filtered complex. The +-- complex's cached exact-value quotient transports pair endpoints to dense +-- ranks, performs the sweep over 'Int', then returns profiles aligned with +-- 'filteredCriticalValues'. Exact births remain authoritative; ranks never +-- escape this derived observation. +persistentBettiAtCriticalValues :: + Ord filtration => + FilteredFiniteChainComplex filtration r -> + [PersistencePair filtration] -> + Either HomologyFailure [Map.Map HomologicalDegree Int] +persistentBettiAtCriticalValues filtered pairs = do + rankedPairs <- traverse (rankPersistencePair filtered) pairs + pure + ( criticalProfilesFromRankedEvents + (length (filteredCriticalValues filtered)) + (rankedPairs >>= rankedPersistencePairEvents) + ) + +criticalProfilesFromRankedEvents :: + Int -> + [(Int, (Int, Int))] -> + [Map.Map HomologicalDegree Int] +criticalProfilesFromRankedEvents criticalRankCount rankedEvents = + let eventDeltas = + Vector.accum + insertRankedBettiEvent + (Vector.replicate criticalRankCount IntMap.empty) + rankedEvents + (_, criticalProfiles) = + List.mapAccumL + applyRankedBettiEvents + IntMap.empty + (Vector.toList eventDeltas) + in criticalProfiles + +criticalBettiTableFromRankedEvents :: + Int -> + Int -> + [(Int, (Int, Int))] -> + CriticalBettiTable +criticalBettiTableFromRankedEvents degreeCount criticalRankCount rankedEvents = + let eventDeltas = + Vector.accum + insertRankedBettiEvent + (Vector.replicate criticalRankCount IntMap.empty) + rankedEvents + degreeIndices = [0 .. degreeCount - 1] + values = + UnboxedVector.create $ do + mutableValues <- + MutableVector.new (max 0 (degreeCount * criticalRankCount)) + _ <- + Vector.ifoldM' + ( \currentBetti rankValue rankEvents -> do + let nextBetti = nextRankedBetti currentBetti rankEvents + traverse_ + ( \degreeValue -> + MutableVector.write + mutableValues + (rankValue * degreeCount + degreeValue) + (IntMap.findWithDefault 0 degreeValue nextBetti) + ) + degreeIndices + pure nextBetti + ) + IntMap.empty + eventDeltas + pure mutableValues + in CriticalBettiTable + { criticalBettiDegreeCount = degreeCount, + criticalBettiRankCount = criticalRankCount, + criticalBettiTableValues = values + } + +rankPersistencePair :: + Ord filtration => + FilteredFiniteChainComplex filtration r -> + PersistencePair filtration -> + Either HomologyFailure (PersistencePair Int) +rankPersistencePair filtered pairValue = do + birthRank <- + maybe + (Left (InvalidTopologyInput "persistence-pair birth is absent from the filtered complex")) + Right + (Map.lookup (persistenceBirth pairValue) (filteredBirthValueRanks filtered)) + deathRank <- + traverse + ( \deathValue -> + maybe + (Left (InvalidTopologyInput "persistence-pair death is absent from the filtered complex")) + Right + (Map.lookup deathValue (filteredBirthValueRanks filtered)) + ) + (persistenceDeath pairValue) + pure + PersistencePair + { persistenceDegree = persistenceDegree pairValue, + persistenceBirth = birthRank, + persistenceDeath = deathRank + } + +rankedPersistencePairEvents :: + PersistencePair Int -> + [(Int, (Int, Int))] +rankedPersistencePairEvents pairValue = + rankedIntervalEvents + (persistenceDegree pairValue) + (persistenceBirth pairValue) + (persistenceDeath pairValue) + +rankedIntervalEvents :: + HomologicalDegree -> + Int -> + Maybe Int -> + [(Int, (Int, Int))] +rankedIntervalEvents degreeValue birthRank deathRank = + case deathRank of + Nothing -> [rankedBettiEventAt 1 birthRank degreeValue] + Just finiteDeathRank -> + if birthRank < finiteDeathRank + then + [ rankedBettiEventAt 1 birthRank degreeValue + , rankedBettiEventAt (-1) finiteDeathRank degreeValue + ] + else [] + +rankedBettiEventAt :: + Int -> + Int -> + HomologicalDegree -> + (Int, (Int, Int)) +rankedBettiEventAt deltaValue criticalRank degreeValue = + (criticalRank, (unHomologicalDegree degreeValue, deltaValue)) + +insertRankedBettiEvent :: + IntMap.IntMap Int -> + (Int, Int) -> + IntMap.IntMap Int +insertRankedBettiEvent eventDeltas (degreeValue, deltaValue) = + IntMap.insertWith + (+) + degreeValue + deltaValue + eventDeltas + +applyRankedBettiEvents :: + IntMap.IntMap Int -> + IntMap.IntMap Int -> + (IntMap.IntMap Int, Map.Map HomologicalDegree Int) +applyRankedBettiEvents currentBetti eventDeltas = + let nextBetti = nextRankedBetti currentBetti eventDeltas + publicBetti = + nextBetti + & IntMap.foldrWithKey + ( \degreeValue countValue remainingDegrees -> + (HomologicalDegree degreeValue, countValue) : remainingDegrees + ) + [] + & Map.fromDistinctAscList + in (nextBetti, publicBetti) + +nextRankedBetti :: + IntMap.IntMap Int -> + IntMap.IntMap Int -> + IntMap.IntMap Int +nextRankedBetti currentBetti eventDeltas = + IntMap.unionWith (+) currentBetti eventDeltas + & IntMap.filter (/= 0) + +requestedBettiSweepPoint :: BettiSweepPoint +requestedBettiSweepPoint = + BettiSweepPoint + { bettiSweepDelta = Map.empty, + bettiSweepRequested = True + } + +insertPersistencePairSweepPoints :: + Ord filtration => + Map.Map filtration BettiSweepPoint -> + PersistencePair filtration -> + Map.Map filtration BettiSweepPoint +insertPersistencePairSweepPoints eventPoints pairValue = + case persistenceDeath pairValue of + Nothing -> insertPersistencePairDelta 1 (persistenceBirth pairValue) pairValue eventPoints + Just deathValue -> + if persistenceBirth pairValue < deathValue + then + eventPoints + & insertPersistencePairDelta 1 (persistenceBirth pairValue) pairValue + & insertPersistencePairDelta (-1) deathValue pairValue + else eventPoints + +insertPersistencePairDelta :: + Ord filtration => + Int -> + filtration -> + PersistencePair filtration -> + Map.Map filtration BettiSweepPoint -> + Map.Map filtration BettiSweepPoint +insertPersistencePairDelta deltaValue threshold pairValue = + Map.insertWith + mergeBettiSweepPoints + threshold + (persistencePairDelta deltaValue pairValue) + +persistencePairDelta :: Int -> PersistencePair filtration -> BettiSweepPoint +persistencePairDelta deltaValue pairValue = + BettiSweepPoint + { bettiSweepDelta = Map.singleton (persistenceDegree pairValue) deltaValue, + bettiSweepRequested = False + } + +mergeBettiSweepPoints :: BettiSweepPoint -> BettiSweepPoint -> BettiSweepPoint +mergeBettiSweepPoints leftPoint rightPoint = + BettiSweepPoint + { bettiSweepDelta = Map.unionWith (+) (bettiSweepDelta leftPoint) (bettiSweepDelta rightPoint), + bettiSweepRequested = bettiSweepRequested leftPoint || bettiSweepRequested rightPoint + } + +applyBettiSweepPoint :: + Map.Map HomologicalDegree Int -> + filtration -> + BettiSweepPoint -> + (Map.Map HomologicalDegree Int, Maybe (Map.Map HomologicalDegree Int)) +applyBettiSweepPoint currentBetti _ pointValue = + let nextBetti = + Map.unionWith (+) currentBetti (bettiSweepDelta pointValue) + & Map.filter (/= 0) + observation = + if bettiSweepRequested pointValue + then Just nextBetti + else Nothing + in (nextBetti, observation) + +persistencePairAliveAt :: + Ord filtration => + filtration -> + PersistencePair filtration -> + Bool +persistencePairAliveAt threshold pairValue = + persistenceBirth pairValue <= threshold + && maybe True (threshold <) (persistenceDeath pairValue) + +persistentPairs :: + Integral r => + FilteredFiniteChainComplex filtration r -> + Either HomologyFailure [PersistencePair filtration] +{-# DEPRECATED persistentPairs "Use mod2PersistentPairs — this computes mod-2 persistence despite its Integral constraint" #-} +persistentPairs = mod2PersistentPairs + +mod2PersistenceTopologyWitness :: + Integral r => + FilteredFiniteChainComplex filtration r -> + Either HomologyFailure (TopologyWitness scaffold spectral filtration coefficient basis) +mod2PersistenceTopologyWitness filtered = do + pairs <- mod2PersistentPairs filtered + pure + emptyTopologyWitness + { topologyPersistencePairs = pairs + } + +persistenceTopologyWitness :: + Integral r => + FilteredFiniteChainComplex filtration r -> + Either HomologyFailure (TopologyWitness scaffold spectral filtration coefficient basis) +{-# DEPRECATED persistenceTopologyWitness "Use mod2PersistenceTopologyWitness — this computes mod-2 persistence despite its Integral constraint" #-} +persistenceTopologyWitness = mod2PersistenceTopologyWitness + +persistenceEssentialBirths :: PersistenceState -> IntSet.IntSet +persistenceEssentialBirths stateValue = + persistenceCreators stateValue + `IntSet.difference` persistencePairedBirths stateValue + +-- | Reduce one homological degree after all higher-degree pairings are known. +-- A cell already paired as the low simplex of a higher-dimensional column is +-- a proven creator, so the standard persistence clearing theorem lets us +-- replace its column by zero rather than rediscovering that fact by reduction. +reducePersistentDegree :: + Int -> + PersistenceState -> + [(Int, IntSet.IntSet)] -> + PersistenceState +reducePersistentDegree totalCellCount stateValue indexedColumns = + let clearedColumnIndices = persistencePairedBirths stateValue + in runST $ do + pivotColumns <- BoxedMutableVector.replicate totalCellCount Nothing + foldM + (reduceOrClearPersistentColumn pivotColumns clearedColumnIndices) + stateValue + indexedColumns + +reduceOrClearPersistentColumn :: + BoxedMutableVector.MVector s (Maybe IntSet.IntSet) -> + IntSet.IntSet -> + PersistenceState -> + (Int, IntSet.IntSet) -> + ST s PersistenceState +reduceOrClearPersistentColumn pivotColumns clearedColumnIndices stateValue indexedColumn@(columnIndexValue, _) = + if IntSet.member columnIndexValue clearedColumnIndices + then + pure + stateValue + { persistenceCreators = IntSet.insert columnIndexValue (persistenceCreators stateValue) + } + else reducePersistentColumn pivotColumns stateValue indexedColumn + +-- | Reduce one sparse column against the dense pivot section for its +-- homological degree. Global row identifiers are assigned from @[0 .. n)@ +-- before boundary materialization, so the fresh arena is a total local view; +-- it cannot escape this degree's sealed 'ST' descent. +reducePersistentColumn :: + BoxedMutableVector.MVector s (Maybe IntSet.IntSet) -> + PersistenceState -> + (Int, IntSet.IntSet) -> + ST s PersistenceState +reducePersistentColumn pivotColumns stateValue (columnIndexValue, initialColumn) = do + reducedColumn <- reduceBoundaryColumn pivotColumns initialColumn + case IntSet.lookupMax reducedColumn of + Nothing -> + pure + stateValue + { persistenceCreators = IntSet.insert columnIndexValue (persistenceCreators stateValue) + } + Just lowValue -> do + BoxedMutableVector.write pivotColumns lowValue (Just reducedColumn) + pure + stateValue + { persistencePairsByIndex = (lowValue, columnIndexValue) : persistencePairsByIndex stateValue, + persistencePairedBirths = IntSet.insert lowValue (persistencePairedBirths stateValue) + } + +-- | Descend through the degree cover. Higher-dimensional columns use the +-- general sparse reducer; the graph boundary uses its equivalent elder-rule +-- union-find section when every mod-two one-cell boundary has cardinality zero +-- or two. Any non-graph fiber is a typed local incompatibility with that +-- specialization and glues back to the general reducer instead. +reducePersistentColumnsByDegree :: + Int -> + [[(Int, IntSet.IntSet)]] -> + [(Int, [Int])] -> + [Int] -> + PersistenceState +reducePersistentColumnsByDegree totalCellCount higherDegreeColumns degreeOneColumns degreeZeroIndices = + let stateAfterHigherDegrees = + foldl' + (reducePersistentDegree totalCellCount) + emptyPersistenceState + higherDegreeColumns + clearedDegreeOneIndices = + persistencePairedBirths stateAfterHigherDegrees + graphReduction = + reduceGraphBoundaryDegree + totalCellCount + clearedDegreeOneIndices + degreeZeroIndices + degreeOneColumns + in case graphReduction of + Just (graphPairs, graphCreators) -> + stateAfterHigherDegrees + { persistencePairsByIndex = graphPairs <> persistencePairsByIndex stateAfterHigherDegrees, + persistencePairedBirths = + IntSet.union + (IntSet.fromList (fmap fst graphPairs)) + (persistencePairedBirths stateAfterHigherDegrees), + persistenceCreators = IntSet.union graphCreators (persistenceCreators stateAfterHigherDegrees) + } + Nothing -> + let sparseDegreeOneColumns = + fmap + (\(columnIndexValue, targets) -> (columnIndexValue, IntSet.fromList targets)) + degreeOneColumns + degreeZeroColumns = + fmap (\columnIndexValue -> (columnIndexValue, IntSet.empty)) degreeZeroIndices + in foldl' + (reducePersistentDegree totalCellCount) + stateAfterHigherDegrees + [sparseDegreeOneColumns, degreeZeroColumns] + +reduceGraphBoundaryDegree :: + Int -> + IntSet.IntSet -> + [Int] -> + [(Int, [Int])] -> + Maybe ([(Int, Int)], IntSet.IntSet) +reduceGraphBoundaryDegree totalCellCount clearedEdgeIndices vertexIndices edgeColumns = + runST $ do + parents <- MutableVector.generate totalCellCount id + foldM + (reduceGraphBoundaryColumn parents clearedEdgeIndices) + (Just ([], IntSet.fromList vertexIndices)) + edgeColumns + +reduceGraphBoundaryColumn :: + MutableVector.MVector s Int -> + IntSet.IntSet -> + Maybe ([(Int, Int)], IntSet.IntSet) -> + (Int, [Int]) -> + ST s (Maybe ([(Int, Int)], IntSet.IntSet)) +reduceGraphBoundaryColumn _ _ Nothing _ = pure Nothing +reduceGraphBoundaryColumn parents clearedEdgeIndices (Just (pairsByIndex, creators)) (edgeIndex, boundaryTargets) = + if IntSet.member edgeIndex clearedEdgeIndices + then pure (Just (pairsByIndex, IntSet.insert edgeIndex creators)) + else + case boundaryTargets of + [] -> pure (Just (pairsByIndex, IntSet.insert edgeIndex creators)) + [firstVertex, secondVertex] -> do + firstRoot <- findGraphComponentRoot parents firstVertex + secondRoot <- findGraphComponentRoot parents secondVertex + case (firstRoot, secondRoot) of + (Just firstRootIndex, Just secondRootIndex) -> + if firstRootIndex == secondRootIndex + then pure (Just (pairsByIndex, IntSet.insert edgeIndex creators)) + else do + let olderRoot = min firstRootIndex secondRootIndex + youngerRoot = max firstRootIndex secondRootIndex + MutableVector.write parents youngerRoot olderRoot + pure (Just ((youngerRoot, edgeIndex) : pairsByIndex, creators)) + _ -> pure Nothing + _ -> pure Nothing + +findGraphComponentRoot :: + MutableVector.MVector s Int -> + Int -> + ST s (Maybe Int) +findGraphComponentRoot parents indexValue = + if indexValue < 0 || indexValue >= MutableVector.length parents + then pure Nothing + else do + parentIndex <- MutableVector.read parents indexValue + if parentIndex == indexValue + then pure (Just indexValue) + else do + rootIndex <- findGraphComponentRoot parents parentIndex + traverse + (\rootValue -> MutableVector.write parents indexValue rootValue >> pure rootValue) + rootIndex + +reduceBoundaryColumn :: + BoxedMutableVector.MVector s (Maybe IntSet.IntSet) -> + IntSet.IntSet -> + ST s IntSet.IntSet +reduceBoundaryColumn pivotColumns columnValue = + case IntSet.lookupMax columnValue of + Nothing -> pure columnValue + Just lowValue -> do + pivotColumn <- BoxedMutableVector.read pivotColumns lowValue + case pivotColumn of + Nothing -> pure columnValue + Just existingPivot -> + reduceBoundaryColumn + pivotColumns + (IntSet.symmetricDifference columnValue existingPivot) + +materializeFinitePersistencePair :: + Vector.Vector (OrderedFilteredCell filtration) -> + Int -> + Int -> + Maybe (PersistencePair filtration) +materializeFinitePersistencePair orderedCellByIndex birthIndexValue deathIndexValue = + case (orderedCellByIndex Vector.!? birthIndexValue, orderedCellByIndex Vector.!? deathIndexValue) of + (Just birthCell, Just deathCell) -> + Just + PersistencePair + { persistenceDegree = cellDegree (orderedCellIdentity birthCell), + persistenceBirth = orderedCellBirth birthCell, + persistenceDeath = Just (orderedCellBirth deathCell) + } + _ -> Nothing + +materializeEssentialPersistencePair :: + Vector.Vector (OrderedFilteredCell filtration) -> + Int -> + Maybe (PersistencePair filtration) +materializeEssentialPersistencePair orderedCellByIndex birthIndexValue = + orderedCellByIndex Vector.!? birthIndexValue + & fmap + ( \birthCell -> + PersistencePair + { persistenceDegree = cellDegree (orderedCellIdentity birthCell), + persistenceBirth = orderedCellBirth birthCell, + persistenceDeath = Nothing + } + ) + +orderedFilteredCells :: + FilteredFiniteChainComplex filtration r -> + [OrderedFilteredCell filtration] +orderedFilteredCells filtered = + filteredRankedCells filtered + & fmap + ( \(cellRefValue, birthValue, birthRank) -> + OrderedFilteredCell + { orderedCellIdentity = cellRefValue, + orderedCellBirth = birthValue, + orderedCellBirthRank = birthRank + } + ) + +mod2PersistentBoundaryColumn :: + Integral r => + FilteredFiniteChainComplex filtration r -> + Map.Map BasisCellRef Int -> + OrderedFilteredCell filtration -> + IntSet.IntSet +mod2PersistentBoundaryColumn filtered globalIndexByCell orderedCell = + let cellRefValue = orderedCellIdentity orderedCell + degreeValue = cellDegree cellRefValue + incidence = incidenceMatrixAt (filteredBaseComplex filtered) degreeValue + in mod2PersistentBoundaryColumnFromEntries (boundaryEntriesBySource incidence) globalIndexByCell orderedCell + +mod2PersistentBoundaryColumnFromIndex :: + Integral r => + Vector.Vector (Vector.Vector [BoundaryEntry r]) -> + Vector.Vector (Vector.Vector (Maybe Int)) -> + OrderedFilteredCell filtration -> + IntSet.IntSet +mod2PersistentBoundaryColumnFromIndex groupedBoundaryEntries globalIndicesByDegree = + IntSet.fromList + . mod2PersistentBoundaryTargetsFromIndex groupedBoundaryEntries globalIndicesByDegree + +mod2PersistentBoundaryTargetsFromIndex :: + Integral r => + Vector.Vector (Vector.Vector [BoundaryEntry r]) -> + Vector.Vector (Vector.Vector (Maybe Int)) -> + OrderedFilteredCell filtration -> + [Int] +mod2PersistentBoundaryTargetsFromIndex groupedBoundaryEntries globalIndicesByDegree orderedCell = + let cellRefValue = orderedCellIdentity orderedCell + degreeValue = unHomologicalDegree (cellDegree cellRefValue) + degreeEntries = + maybe Vector.empty id (groupedBoundaryEntries Vector.!? degreeValue) + targetIndices = + maybe + Vector.empty + id + (globalIndicesByDegree Vector.!? (degreeValue - 1)) + in maybe [] id (degreeEntries Vector.!? cellIndex cellRefValue) + & filter (\entry -> odd (abs (boundaryCoefficient entry))) + & mapMaybe (\entry -> targetIndices Vector.!? (targetIndex entry) >>= id) + +mod2PersistentBoundaryColumnFromEntries :: + Integral r => + Vector.Vector [BoundaryEntry r] -> + Map.Map BasisCellRef Int -> + OrderedFilteredCell filtration -> + IntSet.IntSet +mod2PersistentBoundaryColumnFromEntries groupedEntries globalIndexByCell orderedCell = + let cellRefValue = orderedCellIdentity orderedCell + degreeValue = cellDegree cellRefValue + in maybe [] id (groupedEntries Vector.!? cellIndex cellRefValue) + & filter (\entry -> odd (abs (boundaryCoefficient entry))) + & fmap + ( \entry -> + BasisCellRef + { cellDegree = decrementDegree degreeValue, + cellIndex = targetIndex entry + } + ) + & mapMaybeWithLookup globalIndexByCell + & IntSet.fromList + +boundaryEntriesByDegree :: + FiniteChainComplex r -> + Vector.Vector (Vector.Vector [BoundaryEntry r]) +boundaryEntriesByDegree finite = + dimensionsOf finite + & fmap + ( \degreeValue -> + boundaryEntriesBySource (incidenceMatrixAt finite degreeValue) + ) + & Vector.fromList + +denseGlobalIndicesByDegree :: + FiniteChainComplex r -> + [(Int, OrderedFilteredCell filtration)] -> + Vector.Vector (Vector.Vector (Maybe Int)) +denseGlobalIndicesByDegree finite indexedCells = + let degreeCount = length (dimensionsOf finite) + updatesByDegree = + Vector.accum + (flip (:)) + (Vector.replicate degreeCount []) + (mapMaybe (globalIndexUpdate degreeCount) indexedCells) + in Vector.imap + ( \degreeIndex updates -> + Vector.accum + (\_ globalIndexValue -> Just globalIndexValue) + ( Vector.replicate + (cellCountAtDegree finite (HomologicalDegree degreeIndex)) + Nothing + ) + updates + ) + updatesByDegree + where + globalIndexUpdate :: + Int -> + (Int, OrderedFilteredCell filtration) -> + Maybe (Int, (Int, Int)) + globalIndexUpdate degreeCountValue (globalIndexValue, orderedCell) = + let cellRefValue = orderedCellIdentity orderedCell + degreeIndex = unHomologicalDegree (cellDegree cellRefValue) + in if degreeIndex < 0 || degreeIndex >= degreeCountValue + then Nothing + else Just (degreeIndex, (cellIndex cellRefValue, globalIndexValue)) + +persistentBoundaryColumn :: + Integral r => + FilteredFiniteChainComplex filtration r -> + Map.Map BasisCellRef Int -> + OrderedFilteredCell filtration -> + IntSet.IntSet +{-# DEPRECATED persistentBoundaryColumn "Use mod2PersistentBoundaryColumn — this computes mod-2 boundary despite its Integral constraint" #-} +persistentBoundaryColumn = mod2PersistentBoundaryColumn + +validateBirthUniqueness :: + [(BasisCellRef, filtration)] -> + Map.Map BasisCellRef filtration -> + Either HomologyFailure () +validateBirthUniqueness births birthMap = + if length births == Map.size birthMap + then Right () + else Left (InvalidTopologyInput "duplicate birth assignments for the same cell") + +validateBirthCoverage :: + FiniteChainComplex r -> + Map.Map BasisCellRef filtration -> + Either HomologyFailure () +validateBirthCoverage finite birthMap = + () <$ denseBirthSection finite birthMap + +validateBirthExactness :: + FiniteChainComplex r -> + Map.Map BasisCellRef filtration -> + Either HomologyFailure () +validateBirthExactness finite birthMap = + let basisSet = Set.fromList (allBasisCellRefs finite) + extraKeys = Map.keysSet birthMap `Set.difference` basisSet + in if Set.null extraKeys + then Right () + else + Left + ( InvalidTopologyInput + ( "birth map contains cells absent from the chain complex: " + <> show (Set.toList extraKeys) + ) + ) + +validateFiltrationMonotonicity :: + (Integral r, Ord filtration) => + FiniteChainComplex r -> + Map.Map BasisCellRef filtration -> + Either HomologyFailure () +validateFiltrationMonotonicity finite birthMap = do + denseBirthsByDegree <- denseBirthSection finite birthMap + validateDenseFiltrationMonotonicity finite denseBirthsByDegree + +denseBirthSection :: + FiniteChainComplex r -> + Map.Map BasisCellRef filtration -> + Either HomologyFailure (Vector.Vector (Vector.Vector filtration)) +denseBirthSection finite birthMap = + dimensionsOf finite + & traverse + ( \degreeValue -> + Vector.generateM + (cellCountAtDegree finite degreeValue) + ( \indexValue -> + let cellRefValue = BasisCellRef degreeValue indexValue + in maybe + ( Left + ( InvalidTopologyInput + ("missing filtration value for cell " <> show cellRefValue) + ) + ) + Right + (Map.lookup cellRefValue birthMap) + ) + ) + & fmap Vector.fromList + +validateDenseFiltrationMonotonicity :: + (Integral r, Ord filtration) => + FiniteChainComplex r -> + Vector.Vector (Vector.Vector filtration) -> + Either HomologyFailure () +validateDenseFiltrationMonotonicity finite denseBirthsByDegree = + dimensionsOf finite + & mapMaybe + ( \degreeValue -> + firstFiltrationViolation + denseBirthsByDegree + degreeValue + (incidenceMatrixAt finite degreeValue) + ) + & List.find (const True) + & maybe (Right ()) (Left . InvalidTopologyInput) + +firstFiltrationViolation :: + (Integral r, Ord filtration) => + Vector.Vector (Vector.Vector filtration) -> + HomologicalDegree -> + BoundaryIncidence r -> + Maybe String +firstFiltrationViolation denseBirthsByDegree degreeValue@(HomologicalDegree degreeIndex) incidence + | degreeIndex <= 0 = Nothing + | otherwise = + Vector.foldr + (firstColumnFiltrationViolation denseBirthsByDegree degreeValue) + Nothing + (boundaryEntriesBySource incidence) + +firstColumnFiltrationViolation :: + (Integral r, Ord filtration) => + Vector.Vector (Vector.Vector filtration) -> + HomologicalDegree -> + [BoundaryEntry r] -> + Maybe String -> + Maybe String +firstColumnFiltrationViolation denseBirthsByDegree degreeValue entries remainingViolation = + case List.find (const True) (mapMaybe (filtrationViolationAtEntry denseBirthsByDegree degreeValue) entries) of + Just violation -> Just violation + Nothing -> remainingViolation + +filtrationViolationAtEntry :: + (Integral r, Ord filtration) => + Vector.Vector (Vector.Vector filtration) -> + HomologicalDegree -> + BoundaryEntry r -> + Maybe String +filtrationViolationAtEntry denseBirthsByDegree degreeValue entry + | boundaryCoefficient entry == 0 = Nothing + | otherwise = + let sourceCell = + BasisCellRef + { cellDegree = degreeValue, + cellIndex = sourceIndex entry + } + targetCell = + BasisCellRef + { cellDegree = decrementDegree degreeValue, + cellIndex = targetIndex entry + } + in case (denseBirthAt denseBirthsByDegree sourceCell, denseBirthAt denseBirthsByDegree targetCell) of + (Just sourceBirth, Just targetBirth) + | targetBirth > sourceBirth -> + Just + ( "filtration violates face monotonicity for " + <> show sourceCell + <> " -> " + <> show targetCell + ) + _ -> Nothing + +denseBirthAt :: + Vector.Vector (Vector.Vector filtration) -> + BasisCellRef -> + Maybe filtration +denseBirthAt denseBirthsByDegree cellRefValue = + case cellDegree cellRefValue of + HomologicalDegree degreeIndex -> + denseBirthsByDegree Vector.!? degreeIndex + >>= (\birthsAtDegree -> birthsAtDegree Vector.!? cellIndex cellRefValue) diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -6,9 +6,13 @@ import Data.Kind (Type) import Data.List.NonEmpty (toList) import qualified Data.Map.Strict as Map +import Data.Vector qualified as Vector import Moonlight.Core (mkCapability) import Moonlight.Homology import Moonlight.Homology.Boundary.Finite (mkFiniteChainComplex) +import Moonlight.Homology.Boundary.LinAlg + ( boundaryEntriesBySource, + ) import Moonlight.Homology.Effect.Laws import Moonlight.Homology.Effect.Determinism import BlockSchurSpec qualified @@ -131,6 +135,41 @@ "unordered entries still use canonical semantics" (mkBoundaryIncidence 2 2 entries) (mkBoundaryIncidenceFromOrderedEntries 2 2 entries), + testCase "glues ordered columns to the canonical source projection" $ do + let columns = + Vector.fromList + [ [mkBoundaryEntry 0 0 (2 :: Int), mkBoundaryEntry 0 1 3], + [mkBoundaryEntry 1 1 5] + ] + flattenedEntries = concat (Vector.toList columns) + assertEqual + "ordered columns retain the same canonical sections" + (fmap boundaryEntriesBySource (mkBoundaryIncidence 2 2 flattenedEntries)) + (fmap boundaryEntriesBySource (mkBoundaryIncidenceFromOrderedColumns 2 2 columns)), + testCase "canonicalizes a noncanonical column cover through the general owner" $ do + let columns = + Vector.fromList + [ [mkBoundaryEntry 0 1 (3 :: Int), mkBoundaryEntry 0 0 2], + [mkBoundaryEntry 1 1 5] + ] + flattenedEntries = concat (Vector.toList columns) + assertEqual + "noncanonical columns descend to the canonical source sections" + (fmap boundaryEntriesBySource (mkBoundaryIncidence 2 2 flattenedEntries)) + (fmap boundaryEntriesBySource (mkBoundaryIncidenceFromOrderedColumns 2 2 columns)), + testCase "preserves typed bounds obstructions for malformed columns" $ + assertEqual + "out-of-bounds column entry" + (Left (BoundaryIncidenceEntryOutOfBounds 0 2 2 2)) + ( mkBoundaryIncidenceFromOrderedColumns + 2 + 2 + ( Vector.fromList + [ [mkBoundaryEntry 0 2 (1 :: Int)], + [] + ] + ) + ), testCase "rejects out-of-bounds entries before fast construction" $ assertEqual "out-of-bounds entry" diff --git a/test/facade/CompileFailSpec.hs b/test/facade/CompileFailSpec.hs --- a/test/facade/CompileFailSpec.hs +++ b/test/facade/CompileFailSpec.hs @@ -155,7 +155,7 @@ -- fixture verdict. fixturePackageIds :: [GhcPackageSpec] fixturePackageIds = - [GhcPackageId "moonlight-homology-0.1.0.1-inplace"] + [GhcPackageId "moonlight-homology-0.1.0.2-inplace"] expectRight :: Show left => diff --git a/test/topology/TopologySpec.hs b/test/topology/TopologySpec.hs --- a/test/topology/TopologySpec.hs +++ b/test/topology/TopologySpec.hs @@ -14,6 +14,7 @@ import qualified Moonlight.Homology as H import qualified Moonlight.Homology.Boundary.Finite as H (mkFiniteChainComplex) import Moonlight.Homology.Pure.Topology.Algebra (mkQuotientPresentation) +import Moonlight.Homology.Pure.Topology.Core qualified as HomologyCore import Moonlight.Homology.Pure.Topology.Harmonic (harmonicBasisAt) import TestFixtures ( intervalComplex, @@ -42,6 +43,10 @@ testCase "exact witness recovers torsion and exact classes for Moore complex" testExactWitnessMoore, testCase "projective plane cellular attaching map exposes Z/2 torsion" testProjectivePlaneTorsionAnchor, testCase "persistent witness tracks the essential loop birth" testPersistentTriangleLoop, + testCase "persistence retains a non-binary64 exact filtration order" testExactOrderedPersistence, + testCase "critical Betti tables cancel same-rank interval pairs" testCriticalBettiSameRankCancellation, + testCase "persistence falls back for non-graph-shaped degree-one columns" testPersistentGraphFallback, + testCase "graph specialization preserves higher-degree clearing" testPersistentTetrahedronClearing, testCase "graph witness extracts scaffold and low modes on interval" testGraphWitnessInterval, testCase "graph witness remains convergent on branched five-vertex skeletons" testGraphWitnessBranchedFiveVertexSkeleton, testCase "graph witness seed enriches graph scaffolds with exact loop data" testGraphWitnessSeedTriangle, @@ -263,6 +268,156 @@ length essentialLoops @?= 1 fmap H.persistenceBirth essentialLoops @?= [H.FiltrationValue 2.0] +newtype ExactTestBirth = ExactTestBirth Rational + deriving stock (Eq, Ord, Show) + +testExactOrderedPersistence :: Assertion +testExactOrderedPersistence = do + filteredComplex <- + expectRight + (H.mkFilteredFiniteChainComplex triangleCycleComplex exactTriangleFiltration) + persistencePairs <- expectRight (H.mod2PersistentPairs filteredComplex) + assertCriticalBettiTableAgrees filteredComplex persistencePairs + let criticalValues = fmap ExactTestBirth [0, 1, 2] + H.filteredCriticalValues filteredComplex @?= criticalValues + HomologyCore.filteredBirthValueRanks filteredComplex + @?= Map.fromList (zip criticalValues [0, 1, 2]) + fmap + (\(cellRefValue, _) -> Map.lookup cellRefValue (HomologyCore.filteredCellBirthRanks filteredComplex)) + exactTriangleFiltration + @?= fmap Just [0, 0, 0, 1, 1, 2] + let essentialLoops = + persistencePairs + & filter (\pairValue -> H.persistenceDegree pairValue == H.HomologicalDegree 1) + & filter (isNothing . H.persistenceDeath) + fmap H.persistenceBirth essentialLoops @?= [ExactTestBirth 2] + H.persistentBettiAt (ExactTestBirth 1) persistencePairs + @?= Map.fromList [(H.HomologicalDegree 0, 1)] + H.persistentBettiAt (ExactTestBirth 2) persistencePairs + @?= Map.fromList + [ (H.HomologicalDegree 0, 1), + (H.HomologicalDegree 1, 1) + ] + let requestedThresholds = fmap ExactTestBirth [2, 0, 1, 2] + H.persistentBettiAtMany requestedThresholds persistencePairs + @?= fmap (`H.persistentBettiAt` persistencePairs) requestedThresholds + criticalProfiles <- + expectRight (H.persistentBettiAtCriticalValues filteredComplex persistencePairs) + criticalProfiles + @?= fmap (`H.persistentBettiAt` persistencePairs) criticalValues + let malformedPairs = + [ H.PersistencePair + { H.persistenceDegree = H.HomologicalDegree 0, + H.persistenceBirth = ExactTestBirth 1, + H.persistenceDeath = Just (ExactTestBirth 1) + }, + H.PersistencePair + { H.persistenceDegree = H.HomologicalDegree 1, + H.persistenceBirth = ExactTestBirth 2, + H.persistenceDeath = Just (ExactTestBirth 1) + } + ] + H.persistentBettiAtMany requestedThresholds malformedPairs + @?= fmap (`H.persistentBettiAt` malformedPairs) requestedThresholds + +testCriticalBettiSameRankCancellation :: Assertion +testCriticalBettiSameRankCancellation = do + filteredComplex <- + expectRight + ( H.mkFilteredFiniteChainComplex + intervalComplex + (filtrationByDegree intervalComplex (const (ExactTestBirth 0))) + ) + (persistencePairs, criticalBettiTable) <- + expectRight (H.mod2PersistentPairsWithCriticalBettiTable filteredComplex) + persistencePairs + @?= [ H.PersistencePair + { H.persistenceDegree = H.HomologicalDegree 0, + H.persistenceBirth = ExactTestBirth 0, + H.persistenceDeath = Just (ExactTestBirth 0) + }, + H.PersistencePair + { H.persistenceDegree = H.HomologicalDegree 0, + H.persistenceBirth = ExactTestBirth 0, + H.persistenceDeath = Nothing + } + ] + H.criticalBettiVectors criticalBettiTable @?= [[1, 0]] + +testPersistentGraphFallback :: Assertion +testPersistentGraphFallback = do + singleEndpointGraph <- + expectRight (boundaryOneComplex 2 1 [boundaryEntry 0 0 1]) + filteredComplex <- + expectRight + ( H.mkFilteredFiniteChainComplex + singleEndpointGraph + ( filtrationByDegree singleEndpointGraph $ \degreeValue -> + ExactTestBirth + (if degreeValue == H.HomologicalDegree 0 then 0 else 1) + ) + ) + (persistencePairs, criticalBettiTable) <- + expectRight (H.mod2PersistentPairsWithCriticalBettiTable filteredComplex) + persistencePairs + @?= [ H.PersistencePair + { H.persistenceDegree = H.HomologicalDegree 0, + H.persistenceBirth = ExactTestBirth 0, + H.persistenceDeath = Just (ExactTestBirth 1) + }, + H.PersistencePair + { H.persistenceDegree = H.HomologicalDegree 0, + H.persistenceBirth = ExactTestBirth 0, + H.persistenceDeath = Nothing + } + ] + H.criticalBettiVectors criticalBettiTable @?= [[2, 0], [1, 0]] + +testPersistentTetrahedronClearing :: Assertion +testPersistentTetrahedronClearing = do + filteredComplex <- + expectRight + ( H.mkFilteredFiniteChainComplex + tetrahedronBoundaryComplex + ( filtrationByDegree + tetrahedronBoundaryComplex + (ExactTestBirth . fromIntegral . H.unHomologicalDegree) + ) + ) + persistencePairs <- expectRight (H.mod2PersistentPairs filteredComplex) + assertCriticalBettiTableAgrees filteredComplex persistencePairs + let pairsAt degreeValue = + filter ((== degreeValue) . H.persistenceDegree) persistencePairs + degreeOnePairs = pairsAt (H.HomologicalDegree 1) + degreeTwoPairs = pairsAt (H.HomologicalDegree 2) + length degreeOnePairs @?= 3 + fmap H.persistenceBirth degreeOnePairs @?= replicate 3 (ExactTestBirth 1) + fmap H.persistenceDeath degreeOnePairs @?= replicate 3 (Just (ExactTestBirth 2)) + fmap H.persistenceBirth degreeTwoPairs @?= [ExactTestBirth 2] + fmap H.persistenceDeath degreeTwoPairs @?= [Nothing] + +assertCriticalBettiTableAgrees :: + (Integral r, Ord filtration, Show filtration) => + H.FilteredFiniteChainComplex filtration r -> + [H.PersistencePair filtration] -> + Assertion +assertCriticalBettiTableAgrees filteredComplex persistencePairs = do + (fusedPairs, criticalBettiTable) <- + expectRight (H.mod2PersistentPairsWithCriticalBettiTable filteredComplex) + criticalProfiles <- + expectRight (H.persistentBettiAtCriticalValues filteredComplex persistencePairs) + let H.HomologicalDegree maximumDegree = + H.maxHomologicalDegree (H.filteredBaseComplex filteredComplex) + degreeValues = fmap H.HomologicalDegree [0 .. maximumDegree] + denseProfiles = + fmap + (\profile -> fmap (\degreeValue -> Map.findWithDefault 0 degreeValue profile) degreeValues) + criticalProfiles + fusedPairs @?= persistencePairs + H.criticalBettiDegreeCount criticalBettiTable @?= maximumDegree + 1 + H.criticalBettiRankCount criticalBettiTable @?= length criticalProfiles + H.criticalBettiVectors criticalBettiTable @?= denseProfiles + testGraphWitnessInterval :: Assertion testGraphWitnessInterval = withIntervalObservation $ \intervalObservationValue -> do @@ -856,6 +1011,31 @@ (H.BasisCellRef {H.cellDegree = H.HomologicalDegree 1, H.cellIndex = 1}, H.FiltrationValue 1.0), (H.BasisCellRef {H.cellDegree = H.HomologicalDegree 1, H.cellIndex = 2}, H.FiltrationValue 2.0) ] + +exactTriangleFiltration :: [(H.BasisCellRef, ExactTestBirth)] +exactTriangleFiltration = + [ (H.BasisCellRef {H.cellDegree = H.HomologicalDegree 0, H.cellIndex = 0}, ExactTestBirth 0), + (H.BasisCellRef {H.cellDegree = H.HomologicalDegree 0, H.cellIndex = 1}, ExactTestBirth 0), + (H.BasisCellRef {H.cellDegree = H.HomologicalDegree 0, H.cellIndex = 2}, ExactTestBirth 0), + (H.BasisCellRef {H.cellDegree = H.HomologicalDegree 1, H.cellIndex = 0}, ExactTestBirth 1), + (H.BasisCellRef {H.cellDegree = H.HomologicalDegree 1, H.cellIndex = 1}, ExactTestBirth 1), + (H.BasisCellRef {H.cellDegree = H.HomologicalDegree 1, H.cellIndex = 2}, ExactTestBirth 2) + ] + +filtrationByDegree :: + H.FiniteChainComplex coefficient -> + (H.HomologicalDegree -> filtration) -> + [(H.BasisCellRef, filtration)] +filtrationByDegree finite birthAtDegree = + let H.HomologicalDegree maximumDegree = H.maxHomologicalDegree finite + in concatMap + ( \degreeIndex -> + let degreeValue = H.HomologicalDegree degreeIndex + in fmap + (\cellRefValue -> (cellRefValue, birthAtDegree degreeValue)) + (H.finiteChainBasisRefsAtDegree finite degreeValue) + ) + [0 .. maximumDegree] withIntervalObservation :: (H.TopologyObservationConfig Integer -> Assertion) -> Assertion withIntervalObservation assertion = do