diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,15 @@ All notable changes to `moonlight-homology` are documented here. +## 0.1.0.3 - 2026-08-29 + +- Add checked finite chain maps, arbitrary forward/backward finite chain + zigzags, and exact rational interval decomposition through one streaming + right-filtration descent. Each induced arrow is consumed once rather than + retained globally. Arrow payloads and interval endpoints are functorial; + invalid shapes, endpoint gluing, chain-map laws, and interval multiplicities + remain typed failures. + ## 0.1.0.2 - 2026-08-22 - Generalize `FilteredFiniteChainComplex` from the fixed binary64 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -29,7 +29,8 @@ - **Exact and spectral sequences.** Filtered spectral families with page-by-page reduction and convergence tracking; exact-sequence helpers; Block–Schur reductions. - **Persistence.** Arbitrary ordered one-parameter birth keys, mod-2 - persistence pairs and closed-sublevel Betti queries; two-parameter vocabulary. + persistence pairs and closed-sublevel Betti queries; checked finite chain + maps and exact rational zigzag intervals; 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 @@ -150,6 +151,9 @@ 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. + `mkFiniteChainMapChecked` admits only boundary-commuting maps, + `mkFiniteChainZigzag` glues arbitrary forward/backward diagrams, and + `rationalZigzagIntervals` returns their exact interval decomposition. `BiPersistencePair` carries the two-parameter case. In `Moonlight.Homology.Persistence`. - **Spectral sequences.** `mkSpectralSource` and `spectralFamilyPages` produce the @@ -212,7 +216,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` | Ordered filtered complexes, mod-2 persistence pairs, and barcode Betti queries. | +| `Moonlight.Homology.Persistence` | Ordered filtered complexes, mod-2 persistence, checked chain maps, and exact rational zigzag intervals. | | `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/bench/Main.hs b/bench/Main.hs --- a/bench/Main.hs +++ b/bench/Main.hs @@ -15,6 +15,9 @@ import Test.Tasty.Bench ( defaultMain, ) +import ZigzagPersistence + ( zigzagPersistenceBenchmarks, + ) main :: IO () main = do @@ -23,5 +26,6 @@ putStrLn (benchmarkNotice includeLarge include100k) defaultMain [ morseSpectralBenchmarks includeLarge, - sparseSpectralBenchmarks includeLarge include100k + sparseSpectralBenchmarks includeLarge include100k, + zigzagPersistenceBenchmarks includeLarge ] diff --git a/bench/topology/ZigzagPersistence.hs b/bench/topology/ZigzagPersistence.hs new file mode 100644 --- /dev/null +++ b/bench/topology/ZigzagPersistence.hs @@ -0,0 +1,139 @@ +module ZigzagPersistence + ( zigzagPersistenceBenchmarks, + ) +where + +import Data.List qualified as List +import Moonlight.Homology.Boundary + ( FiniteChainComplex, + emptyBoundaryIncidence, + emptyBoundaryIncidenceOf, + mkBoundaryEntry, + mkBoundaryIncidence, + ) +import Moonlight.Homology.Boundary.Finite (degreeCardinality, mkFiniteChainComplex) +import Moonlight.Homology.Chain (HomologicalDegree (..)) +import Moonlight.Homology.Persistence + ( FiniteChainMap, + FiniteChainZigzag, + ZigzagArrow (..), + ZigzagInterval (..), + mkFiniteChainMapChecked, + mkFiniteChainZigzag, + rationalZigzagIntervals, + ) +import Test.Tasty.Bench (Benchmark, bench, bgroup, whnf) + +zigzagPersistenceBenchmarks :: Bool -> Benchmark +zigzagPersistenceBenchmarks includeLarge = + case traverse benchmarkCase ([("vertices-9", 9), ("vertices-17", 17)] <> [("vertices-33", 33) | includeLarge]) of + Left fixtureFailure -> bench "invalid-zigzag-fixture" (whnf id fixtureFailure) + Right benchmarks -> bgroup "zigzag-persistence" benchmarks + +benchmarkCase :: (String, Int) -> Either String Benchmark +benchmarkCase (caseName, vertexCount) = do + let seeds = benchmarkIntervalSeeds vertexCount + diagram <- intervalSumDiagram vertexCount seeds + intervals <- firstShow "zigzag preflight failed" (rationalZigzagIntervals diagram) + let recoveredMultiplicity = sum (fmap zigzagIntervalMultiplicity intervals) + if recoveredMultiplicity /= length seeds + then Left ("zigzag preflight recovered " <> show recoveredMultiplicity <> " of " <> show (length seeds) <> " interval summands") + else Right (bench caseName (whnf zigzagChecksum diagram)) + +zigzagChecksum :: FiniteChainZigzag Int -> Int +zigzagChecksum diagram = + either + (const minBound) + ( foldl' + ( \checksum interval -> + checksum * 16777619 + + unHomologicalDegree (zigzagIntervalDegree interval) * 31 + + zigzagIntervalFirst interval * 17 + + zigzagIntervalLast interval * 7 + + zigzagIntervalMultiplicity interval + ) + 2166136261 + ) + (rationalZigzagIntervals diagram) + +type IntervalSeed = (Int, Int, Int) + +benchmarkIntervalSeeds :: Int -> [IntervalSeed] +benchmarkIntervalSeeds vertexCount = + fmap (\seedIndex -> (seedIndex, 0, vertexCount - 1)) [0 .. 3] + <> fmap + ( \seedIndex -> + let firstIndex = seedIndex `mod` vertexCount + maximumSpan = vertexCount - firstIndex + intervalWidth = 1 + (seedIndex * 7 `mod` maximumSpan) + in (seedIndex + 4, firstIndex, firstIndex + intervalWidth - 1) + ) + [0 .. 3 * vertexCount - 1] + +intervalSumDiagram :: Int -> [IntervalSeed] -> Either String (FiniteChainZigzag Int) +intervalSumDiagram vertexCount seeds + | vertexCount <= 0 = Left "zigzag benchmark requires at least one vertex" + | otherwise = do + let activeSeeds stageIndex = + filter (\(_, firstIndex, lastIndex) -> firstIndex <= stageIndex && stageIndex <= lastIndex) seeds + complexAt stageIndex = zeroComplex (length (activeSeeds stageIndex)) + coordinates sourceIndex targetIndex = + [ (sourceCoordinate, targetCoordinate) + | (sourceCoordinate, seed) <- zip [0 :: Int ..] (activeSeeds sourceIndex), + Just targetCoordinate <- [List.elemIndex seed (activeSeeds targetIndex)] + ] + arrowAt arrowIndex = + if even arrowIndex + then + ForwardArrow + <$> coordinateMap + (complexAt arrowIndex) + (complexAt (arrowIndex + 1)) + (coordinates arrowIndex (arrowIndex + 1)) + else + BackwardArrow + <$> coordinateMap + (complexAt (arrowIndex + 1)) + (complexAt arrowIndex) + (coordinates (arrowIndex + 1) arrowIndex) + arrows <- traverse arrowAt [0 .. vertexCount - 2] + firstShow "invalid zigzag fixture" (mkFiniteChainZigzag (complexAt 0) arrows) + +zeroComplex :: Int -> FiniteChainComplex Int +zeroComplex dimensionValue = + mkFiniteChainComplex (HomologicalDegree 0) $ \degreeValue -> + case degreeValue of + HomologicalDegree 0 -> emptyBoundaryIncidenceOf (fromIntegral dimensionValue) 0 + _ -> emptyBoundaryIncidence + +coordinateMap :: + FiniteChainComplex Int -> + FiniteChainComplex Int -> + [(Int, Int)] -> + Either String (FiniteChainMap Int) +coordinateMap sourceComplex targetComplex coordinates = do + degreeZeroMap <- + firstShow + "invalid coordinate map" + ( mkBoundaryIncidence + (fromIntegral (sourceDimension sourceComplex)) + (fromIntegral (sourceDimension targetComplex)) + ( fmap + (\(sourceIndex, targetIndex) -> mkBoundaryEntry (fromIntegral sourceIndex) (fromIntegral targetIndex) (1 :: Int)) + coordinates + ) + ) + firstShow + "invalid chain map" + ( mkFiniteChainMapChecked sourceComplex targetComplex $ \degreeValue -> + case degreeValue of + HomologicalDegree 0 -> degreeZeroMap + _ -> emptyBoundaryIncidence + ) + +sourceDimension :: FiniteChainComplex Int -> Int +sourceDimension finite = degreeCardinality finite (HomologicalDegree 0) + +firstShow :: Show failure => String -> Either failure value -> Either String value +firstShow contextMessage = + either (Left . ((contextMessage <> ": ") <>) . show) Right diff --git a/moonlight-homology.cabal b/moonlight-homology.cabal --- a/moonlight-homology.cabal +++ b/moonlight-homology.cabal @@ -1,10 +1,10 @@ cabal-version: 3.0 name: moonlight-homology -version: 0.1.0.2 +version: 0.1.0.3 homepage: https://github.com/PaleRoses/moonlight bug-reports: https://github.com/PaleRoses/moonlight/issues synopsis: Chain complexes, phase-gated homology interfaces, and spectral scaffolding. -description: Finite chain complexes, validated boundary matrices, field and Smith-normal-form rank backends, phase-gated Betti numbers and spectral sequences, discrete Morse reductions, and persistence helpers. +description: Finite chain complexes, validated boundary matrices, field and Smith-normal-form rank backends, phase-gated Betti numbers and spectral sequences, discrete Morse reductions, ordered persistence, and exact finite zigzag persistence. license: MIT license-file: LICENSE copyright: (c) 2026 Blue Rose @@ -162,6 +162,7 @@ Moonlight.Homology.Pure.Topology.SparseAlgebra Moonlight.Homology.Pure.Topology.Spectral Moonlight.Homology.Pure.Topology.Target + Moonlight.Homology.Pure.Topology.Zigzag Moonlight.Homology.Pure.TopologyObserver Moonlight.Homology.Pure.TopologyView build-depends: @@ -273,6 +274,7 @@ SpectralSpec TestFixtures TopologySpec + ZigzagSpec build-depends: base >= 4.22 && < 5 , containers >= 0.6 && < 0.9 @@ -300,6 +302,7 @@ other-modules: MorseSpectral SparseSpectral + ZigzagPersistence ghc-options: -O2 -rtsopts ghc-prof-options: -fprof-auto-top build-depends: @@ -317,5 +320,5 @@ source-repository this type: git location: https://github.com/PaleRoses/moonlight.git - tag: moonlight-homology-0.1.0.2 + tag: moonlight-homology-0.1.0.3 subdir: moonlight-homology diff --git a/src-matrix/Moonlight/Homology/Pure/Matrix/SparseLinAlg.hs b/src-matrix/Moonlight/Homology/Pure/Matrix/SparseLinAlg.hs --- a/src-matrix/Moonlight/Homology/Pure/Matrix/SparseLinAlg.hs +++ b/src-matrix/Moonlight/Homology/Pure/Matrix/SparseLinAlg.hs @@ -3,7 +3,10 @@ SparseMatrix (..), sparseBoundaryMatrixWith, sparseBoundaryMatrix, + sparseBoundaryColumns, sparseTransposeMatrix, + sparseMatrixVectorProduct, + sparseLinearCombination, SparseRref (..), sparseRref, sparseKernelBasisFromRref, @@ -15,12 +18,15 @@ sparseEchelonBasis, sparseEchelonContains, sparseEchelonRank, + sparseExtendEchelonBasis, sparseIndependentModulo, sparseIndependentModuloWithBasis, sparseSpanRank, SparseCoordinateBasis (..), sparseCoordinateBasis, sparseCoordinatesInBasis, + SparseColumnEchelon (..), + sparseColumnEchelon, compactSparseRow, scaleSparseRow, addScaledSparseRow, @@ -38,10 +44,12 @@ import Data.Kind (Type) import qualified Data.List as List import Data.Ratio (denominator, numerator) +import Data.Vector qualified as Vector import Moonlight.Homology.Boundary.LinAlg ( BoundaryIncidence, boundaryCoefficient, boundaryEntries, + boundaryEntriesBySource, sourceCardinality, sourceIndex, targetCardinality, @@ -93,6 +101,18 @@ sparseBoundaryMatrix = sparseBoundaryMatrixWith fromIntegral +-- | Exact source columns read directly from the incidence owner's cached +-- source cover. Consumers needing an image section avoid constructing a row +-- matrix merely to transpose it again. +sparseBoundaryColumns :: Integral r => BoundaryIncidence r -> Vector.Vector SparseRow +sparseBoundaryColumns = + fmap + ( IntMap.fromDistinctAscList + . fmap (\entry -> (targetIndex entry, fromIntegral (boundaryCoefficient entry))) + ) + . boundaryEntriesBySource +{-# INLINE sparseBoundaryColumns #-} + sparseTransposeMatrix :: SparseMatrix -> SparseMatrix sparseTransposeMatrix matrix = let transposedBuckets = @@ -121,17 +141,45 @@ smColumnCount = newColumnCount } +sparseMatrixVectorProduct :: SparseMatrix -> SparseRow -> SparseRow +sparseMatrixVectorProduct matrix vectorValue = + smRows matrix + & zip [0 :: Int ..] + & fmap (\(rowIndex, rowValue) -> (rowIndex, sparseDotProduct rowValue vectorValue)) + & filter ((/= 0) . snd) + & IntMap.fromDistinctAscList +{-# INLINE sparseMatrixVectorProduct #-} + +-- | The exact linear combination of an indexed vector section. Work is +-- proportional to the selected columns rather than to the ambient matrix. +sparseLinearCombination :: Vector.Vector SparseRow -> SparseRow -> SparseRow +sparseLinearCombination columns = + IntMap.foldlWithKey' + ( \result columnIndex coefficient -> + maybe + result + (addScaledSparseRow coefficient result) + (columns Vector.!? columnIndex) + ) + IntMap.empty +{-# INLINE sparseLinearCombination #-} + +sparseDotProduct :: SparseRow -> SparseRow -> Rational +sparseDotProduct leftRow rightRow = + IntMap.foldlWithKey' + (\total coordinate coefficient -> total + coefficient * IntMap.findWithDefault 0 coordinate rightRow) + 0 + leftRow + type SparseRref :: Type -data SparseRref = SparseRref - { srrefPivots :: ![(Int, SparseRow)], - srrefColumnCount :: !Int +newtype SparseRref = SparseRref + { srrefPivots :: [(Int, SparseRow)] } deriving stock (Eq, Show) type SparseEchelonBasis :: Type -data SparseEchelonBasis = SparseEchelonBasis - { sebColumnCount :: !Int, - sebPivotRows :: !(IntMap SparseRow) +newtype SparseEchelonBasis = SparseEchelonBasis + { sebPivotRows :: IntMap SparseRow } deriving stock (Eq, Show) @@ -143,10 +191,8 @@ deriving stock (Eq, Show) type SparseCoordinateBasis :: Type -data SparseCoordinateBasis = SparseCoordinateBasis - { scbAmbientDimension :: !Int, - scbGeneratorCount :: !Int, - scbPivotRows :: !(IntMap SparseCoordinatePivot) +newtype SparseCoordinateBasis = SparseCoordinateBasis + { scbPivotRows :: IntMap SparseCoordinatePivot } deriving stock (Eq, Show) @@ -157,6 +203,15 @@ } deriving stock (Eq, Show) +-- | A column reduction together with the domain transformations that witness +-- its kernel and image pivots. Pivot coordinates are ascending. +type SparseColumnEchelon :: Type +data SparseColumnEchelon = SparseColumnEchelon + { sparseColumnKernelBasis :: ![SparseRow], + sparseColumnPivotPreimages :: ![(Int, SparseRow)] + } + deriving stock (Eq, Show) + type SparseSupportIndex :: Type data SparseSupportIndex = SparseSupportIndex { ssiColumnRows :: !(IntMap IntSet.IntSet), @@ -200,10 +255,7 @@ -- exactly the order 'canonicalRrefPivots' consumes (latest pivot -- first); no re-reversal is needed on either side. pivots = canonicalRrefPivots (sesSelectedPivots finalState) - in SparseRref - { srrefPivots = pivots, - srrefColumnCount = smColumnCount matrix - } + in SparseRref {srrefPivots = pivots} initialSparseEliminationState :: SparseMatrix -> SparseEliminationState initialSparseEliminationState matrix = @@ -752,53 +804,39 @@ (\columnIndex -> IntMap.findWithDefault IntMap.empty columnIndex selectedColumns) validPivotColumns -sparseIndependentModulo :: Int -> [SparseRow] -> [SparseRow] -> [SparseRow] -sparseIndependentModulo ambientDimension imageBasis kernelBasis = - sparseIndependentModuloWithBasis (sparseEchelonBasis ambientDimension imageBasis) kernelBasis +sparseIndependentModulo :: [SparseRow] -> [SparseRow] -> [SparseRow] +sparseIndependentModulo imageBasis kernelBasis = + sparseIndependentModuloWithBasis (sparseEchelonBasis imageBasis) kernelBasis sparseIndependentModuloWithBasis :: SparseEchelonBasis -> [SparseRow] -> [SparseRow] sparseIndependentModuloWithBasis spanBasis kernelBasis = - let initialSelection = - SparseModuloSelection - { smsSpanBasis = spanBasis, - smsSelectedRows = [] - } - in kernelBasis - & List.foldl' selectIndependentModulo initialSelection - & reverse . smsSelectedRows - -type SparseModuloSelection :: Type -data SparseModuloSelection = SparseModuloSelection - { smsSpanBasis :: !SparseEchelonBasis, - smsSelectedRows :: ![SparseRow] - } - deriving stock (Eq, Show) + sparseExtendEchelonBasis spanBasis kernelBasis + & fst -selectIndependentModulo :: SparseModuloSelection -> SparseRow -> SparseModuloSelection -selectIndependentModulo selection candidateVector = - case adjoinSparseEchelonRow (smsSpanBasis selection) candidateVector of - (Nothing, unchangedBasis) -> - selection {smsSpanBasis = unchangedBasis} - (Just _residualVector, extendedBasis) -> - SparseModuloSelection - { smsSpanBasis = extendedBasis, - smsSelectedRows = candidateVector : smsSelectedRows selection - } +-- | Extend an admitted echelon basis with the independent members of a +-- candidate section. The selected rows retain input order, and the returned +-- basis carries the same extension so callers never repeat the elimination. +sparseExtendEchelonBasis :: SparseEchelonBasis -> [SparseRow] -> ([SparseRow], SparseEchelonBasis) +sparseExtendEchelonBasis spanBasis candidateRows = + let (reversedSelections, extendedBasis) = + List.foldl' selectIndependent ([], spanBasis) candidateRows + in (reverse reversedSelections, extendedBasis) + where + selectIndependent (selectedRows, basis) candidateVector = + case adjoinSparseEchelonRow basis candidateVector of + (Nothing, unchangedBasis) -> (selectedRows, unchangedBasis) + (Just _, extendedBasis) -> (candidateVector : selectedRows, extendedBasis) -sparseEchelonBasis :: Int -> [SparseRow] -> SparseEchelonBasis -sparseEchelonBasis ambientDimension = +sparseEchelonBasis :: [SparseRow] -> SparseEchelonBasis +sparseEchelonBasis = List.foldl' ( \basis rowValue -> snd (adjoinSparseEchelonRow basis rowValue) ) - (emptySparseEchelonBasis ambientDimension) + emptySparseEchelonBasis -emptySparseEchelonBasis :: Int -> SparseEchelonBasis -emptySparseEchelonBasis ambientDimension = - SparseEchelonBasis - { sebColumnCount = ambientDimension, - sebPivotRows = IntMap.empty - } +emptySparseEchelonBasis :: SparseEchelonBasis +emptySparseEchelonBasis = SparseEchelonBasis IntMap.empty sparseEchelonContains :: SparseEchelonBasis -> SparseRow -> Bool sparseEchelonContains basis = @@ -836,24 +874,25 @@ (compactSparseRow rowValue) (sebPivotRows basis) -sparseSpanRank :: Int -> [SparseRow] -> Int -sparseSpanRank ambientDimension vectorList = - sparseEchelonRank (sparseEchelonBasis ambientDimension vectorList) +sparseSpanRank :: [SparseRow] -> Int +sparseSpanRank = sparseEchelonRank . sparseEchelonBasis -sparseCoordinateBasis :: Int -> [SparseRow] -> SparseCoordinateBasis -sparseCoordinateBasis ambientDimension generatorRows = +sparseCoordinateBasis :: [SparseRow] -> SparseCoordinateBasis +sparseCoordinateBasis generatorRows = generatorRows & zip [0 :: Int ..] & List.foldl' - adjoinSparseCoordinateGenerator - SparseCoordinateBasis - { scbAmbientDimension = ambientDimension, - scbGeneratorCount = length generatorRows, - scbPivotRows = IntMap.empty - } + (\basis generator -> snd (extendSparseCoordinateBasis basis generator)) + (SparseCoordinateBasis IntMap.empty) -adjoinSparseCoordinateGenerator :: SparseCoordinateBasis -> (Int, SparseRow) -> SparseCoordinateBasis -adjoinSparseCoordinateGenerator basis (generatorIndex, generatorRow) = +-- | Reduce one generator while carrying its coordinates in the original +-- domain. A dependent generator returns the resulting kernel relation; +-- an independent generator extends the image basis instead. +extendSparseCoordinateBasis :: + SparseCoordinateBasis -> + (Int, SparseRow) -> + (Maybe SparseRow, SparseCoordinateBasis) +extendSparseCoordinateBasis basis (generatorIndex, generatorRow) = let residual = reduceSparseCoordinateGenerator basis @@ -861,18 +900,44 @@ { scrVector = compactSparseRow generatorRow, scrCoordinates = IntMap.singleton generatorIndex 1 } - in case rowLeadingColumn (scrVector residual) of - Nothing -> basis - Just pivotColumn -> + in case IntMap.lookupMin (scrVector residual) of + Nothing -> (Just (scrCoordinates residual), basis) + Just (pivotColumn, pivotCoefficient) -> let pivot = normalizeSparseCoordinatePivot - pivotColumn + pivotCoefficient residual - in basis - { scbPivotRows = - IntMap.insert pivotColumn pivot (scbPivotRows basis) - } + in ( Nothing, + basis + { scbPivotRows = + IntMap.insert pivotColumn pivot (scbPivotRows basis) + } + ) +-- | Reduce a sparse column map once, retaining an exact basis of its kernel +-- and a source preimage for every image pivot. This is the compositional +-- alternative to rebuilding a kernel for every prefix of a target flag. +sparseColumnEchelon :: Vector.Vector SparseRow -> SparseColumnEchelon +sparseColumnEchelon columns = + let initialBasis = SparseCoordinateBasis IntMap.empty + (reversedKernelBasis, imageBasis) = + Vector.ifoldl' + ( \(kernelBasis, basis) columnIndex column -> + case extendSparseCoordinateBasis basis (columnIndex, column) of + (Nothing, extendedBasis) -> (kernelBasis, extendedBasis) + (Just kernelVector, unchangedBasis) -> + (kernelVector : kernelBasis, unchangedBasis) + ) + ([], initialBasis) + columns + in SparseColumnEchelon + { sparseColumnKernelBasis = reverse reversedKernelBasis, + sparseColumnPivotPreimages = + fmap + (\(pivotColumn, pivot) -> (pivotColumn, scpCoordinates pivot)) + (IntMap.toAscList (scbPivotRows imageBasis)) + } + reduceSparseCoordinateGenerator :: SparseCoordinateBasis -> SparseCoordinateResidual -> @@ -891,41 +956,23 @@ eliminateCoordinateGeneratorPivot residual pivotColumn pivot = case IntMap.lookup pivotColumn (scrVector residual) of Nothing -> residual - Just coefficient - | coefficient == 0 -> residual - | otherwise -> - SparseCoordinateResidual - { scrVector = - eliminateColumnFromRow pivotColumn (scpVector pivot) (scrVector residual), - scrCoordinates = - addScaledSparseRow - (negate coefficient) - (scrCoordinates residual) - (scpCoordinates pivot) - } + Just coefficient -> + SparseCoordinateResidual + { scrVector = + addScaledSparseRow (negate coefficient) (scrVector residual) (scpVector pivot), + scrCoordinates = + addScaledSparseRow (negate coefficient) (scrCoordinates residual) (scpCoordinates pivot) + } normalizeSparseCoordinatePivot :: - Int -> + Rational -> SparseCoordinateResidual -> SparseCoordinatePivot -normalizeSparseCoordinatePivot pivotColumn residual = - case IntMap.lookup pivotColumn (scrVector residual) of - Nothing -> - SparseCoordinatePivot - { scpVector = scrVector residual, - scpCoordinates = scrCoordinates residual - } - Just pivotCoefficient - | pivotCoefficient == 0 -> - SparseCoordinatePivot - { scpVector = scrVector residual, - scpCoordinates = scrCoordinates residual - } - | otherwise -> - SparseCoordinatePivot - { scpVector = scaleSparseRow (recip pivotCoefficient) (scrVector residual), - scpCoordinates = scaleSparseRow (recip pivotCoefficient) (scrCoordinates residual) - } +normalizeSparseCoordinatePivot pivotCoefficient residual = + SparseCoordinatePivot + { scpVector = scaleSparseRow (recip pivotCoefficient) (scrVector residual), + scpCoordinates = scaleSparseRow (recip pivotCoefficient) (scrCoordinates residual) + } sparseCoordinatesInBasis :: SparseCoordinateBasis -> SparseRow -> Maybe SparseRow sparseCoordinatesInBasis basis rowValue = @@ -949,18 +996,13 @@ eliminateCoordinateCandidatePivot residual pivotColumn pivot = case IntMap.lookup pivotColumn (scrVector residual) of Nothing -> residual - Just coefficient - | coefficient == 0 -> residual - | otherwise -> - SparseCoordinateResidual - { scrVector = - eliminateColumnFromRow pivotColumn (scpVector pivot) (scrVector residual), - scrCoordinates = - addScaledSparseRow - coefficient - (scrCoordinates residual) - (scpCoordinates pivot) - } + Just coefficient -> + SparseCoordinateResidual + { scrVector = + addScaledSparseRow (negate coefficient) (scrVector residual) (scpVector pivot), + scrCoordinates = + addScaledSparseRow coefficient (scrCoordinates residual) (scpCoordinates pivot) + } scaleSparseRow :: Rational -> SparseRow -> SparseRow scaleSparseRow scalarValue = 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 @@ -124,6 +124,23 @@ finiteAbelianExactOrderElementCount, isPrime, matchesOptional, + ChainMapEndpoint (..), + ZigzagDirection (..), + ZigzagFailure (..), + FiniteChainMap, + finiteChainMapSource, + finiteChainMapTarget, + finiteChainMapAt, + mkFiniteChainMapChecked, + ZigzagArrow (..), + zigzagArrowDirection, + FiniteChainZigzag, + mkFiniteChainZigzag, + finiteChainZigzagComplexes, + finiteChainZigzagArrows, + ZigzagInterval (..), + rationalZigzagIntervals, + zigzagBettiAt, DegreeSelection (..), GradedAggregation (..), GradedQuery (..), 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 @@ -18,6 +18,23 @@ persistentBettiAtMany, persistentBettiAtCriticalValues, mod2PersistenceTopologyWitness, + ChainMapEndpoint (..), + ZigzagDirection (..), + ZigzagFailure (..), + FiniteChainMap, + finiteChainMapSource, + finiteChainMapTarget, + finiteChainMapAt, + mkFiniteChainMapChecked, + ZigzagArrow (..), + zigzagArrowDirection, + FiniteChainZigzag, + mkFiniteChainZigzag, + finiteChainZigzagComplexes, + finiteChainZigzagArrows, + ZigzagInterval (..), + rationalZigzagIntervals, + zigzagBettiAt, ) where @@ -44,6 +61,25 @@ persistentBettiAtMany, persistentBettiAtCriticalValues, mod2PersistenceTopologyWitness, + ) +import Moonlight.Homology.Pure.Topology.Zigzag + ( ChainMapEndpoint (..), + FiniteChainMap, + FiniteChainZigzag, + ZigzagArrow (..), + ZigzagDirection (..), + ZigzagFailure (..), + ZigzagInterval (..), + finiteChainMapAt, + finiteChainMapSource, + finiteChainMapTarget, + finiteChainZigzagArrows, + finiteChainZigzagComplexes, + mkFiniteChainMapChecked, + mkFiniteChainZigzag, + rationalZigzagIntervals, + zigzagArrowDirection, + zigzagBettiAt, ) -- | A cell with two independent filtration parameters. Forward-looking diff --git a/src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Build.hs b/src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Build.hs --- a/src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Build.hs +++ b/src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Build.hs @@ -468,7 +468,7 @@ reducedNumerator <- reduceBasisChecked ambientDimension numeratorBasis reducedDenominator <- reduceBasisChecked ambientDimension denominatorBasis assertDenominatorSubset bidegreeValue ambientDimension reducedNumerator reducedDenominator - let quotientBasis = independentModuloBasis ambientDimension reducedDenominator reducedNumerator + let quotientBasis = independentModuloBasis reducedDenominator reducedNumerator presentation = presentationFromSparseBases (bidegreeTotalDegree bidegreeValue) @@ -509,7 +509,7 @@ [AmbientVector] -> Either HomologyFailure () assertDenominatorSubset bidegreeValue ambientDimension numeratorBasis denominatorBasis = - case firstVectorOutsideSpan ambientDimension numeratorBasis denominatorBasis of + case firstVectorOutsideSpan numeratorBasis denominatorBasis of Nothing -> Right () Just denominatorVector -> Left diff --git a/src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Linear.hs b/src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Linear.hs --- a/src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Linear.hs +++ b/src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Linear.hs @@ -54,6 +54,7 @@ sparseImageBasisOf, sparseIndependentModulo, sparseKernelBasisOf, + sparseMatrixVectorProduct, sparseRowLookup, sparseRowToDense, sparseSpanRank, @@ -284,7 +285,7 @@ reduceBasisChecked :: Int -> [AmbientVector] -> Either HomologyFailure [AmbientVector] reduceBasisChecked ambientDimension basisVectors = do compactBasis <- traverse (validateAmbientVector ambientDimension) basisVectors - pure (independentModuloBasis ambientDimension [] compactBasis) + pure (independentModuloBasis [] compactBasis) intersectionBasisChecked :: Int -> @@ -323,19 +324,19 @@ imageBasisOfMatrix = sparseImageBasisOf -independentModuloBasis :: Int -> [AmbientVector] -> [AmbientVector] -> [AmbientVector] +independentModuloBasis :: [AmbientVector] -> [AmbientVector] -> [AmbientVector] independentModuloBasis = sparseIndependentModulo -firstVectorOutsideSpan :: Int -> [AmbientVector] -> [AmbientVector] -> Maybe AmbientVector -firstVectorOutsideSpan ambientDimension spanBasis candidateVectors = +firstVectorOutsideSpan :: [AmbientVector] -> [AmbientVector] -> Maybe AmbientVector +firstVectorOutsideSpan spanBasis candidateVectors = let echelonBasis = - sparseEchelonBasis ambientDimension spanBasis + sparseEchelonBasis spanBasis in List.find (not . sparseEchelonContains echelonBasis) candidateVectors -spanRankOfBasis :: Int -> [AmbientVector] -> Int +spanRankOfBasis :: [AmbientVector] -> Int spanRankOfBasis = sparseSpanRank @@ -438,18 +439,7 @@ Either HomologyFailure AmbientVector applySparseMatrixChecked matrixValue vectorValue = do compactVector <- validateMatrixVector matrixValue vectorValue - pure - ( smRows matrixValue - & zip [0 :: Int ..] - & List.foldl' - ( \imageRow (rowIndex, rowValue) -> - let coefficientValue = sparseDot rowValue compactVector - in if coefficientValue == 0 - then imageRow - else IntMap.insert rowIndex coefficientValue imageRow - ) - IntMap.empty - ) + pure (sparseMatrixVectorProduct matrixValue compactVector) sparseLinearCombination :: [AmbientVector] -> AmbientVector -> AmbientVector sparseLinearCombination basisVectors coefficients = @@ -525,15 +515,6 @@ negateSparseRow :: AmbientVector -> AmbientVector negateSparseRow = scaleSparseRow (-1) - -sparseDot :: AmbientVector -> AmbientVector -> Rational -sparseDot leftRow rightRow = - IntMap.foldlWithKey' - ( \dotValue columnIndex coefficientValue -> - dotValue + coefficientValue * sparseRowLookup columnIndex rightRow - ) - 0 - leftRow elementAt :: Int -> [a] -> Maybe a elementAt indexValue _ diff --git a/src-topology/Moonlight/Homology/Pure/Topology/Algebra.hs b/src-topology/Moonlight/Homology/Pure/Topology/Algebra.hs --- a/src-topology/Moonlight/Homology/Pure/Topology/Algebra.hs +++ b/src-topology/Moonlight/Homology/Pure/Topology/Algebra.hs @@ -141,7 +141,6 @@ presentationDenominatorBasis = denominatorBasis, presentationCoordinateBasis = sparseCoordinateBasis - ambientDimension (fmap sparseRowFromDense (basisVectors <> denominatorBasis)), presentationRepresentatives = representatives } @@ -173,7 +172,7 @@ sparseMatrixFromRows (matrixColumnCount incomingMatrix) incomingMatrix kernelBasis = sparseKernelBasisOf ambientDimension currentSparse imageBasis = sparseImageBasisOf incomingSparse - quotientBasis = sparseIndependentModulo ambientDimension imageBasis kernelBasis + quotientBasis = sparseIndependentModulo imageBasis kernelBasis in quotientBasis & fmap (sparseVectorToRepresentative degreeValue) @@ -261,7 +260,7 @@ [RepresentativeCycle Rational Int] representativeCyclesOverQPrepared finite preparedBoundaries = dimensionsOf finite - >>= homologyBasisAtPrepared finite preparedBoundaries + >>= homologyBasisAtPrepared preparedBoundaries representativeCocyclesOverQPrepared :: FiniteChainComplex r -> @@ -269,38 +268,34 @@ [RepresentativeCocycle Rational Int] representativeCocyclesOverQPrepared finite preparedBoundaries = dimensionsOf finite - >>= cohomologyBasisAtPrepared finite preparedBoundaries + >>= cohomologyBasisAtPrepared preparedBoundaries homologyBasisAtPrepared :: - FiniteChainComplex r -> IntMap.IntMap RationalBoundaryDecomposition -> HomologicalDegree -> [RepresentativeCycle Rational Int] -homologyBasisAtPrepared finite preparedBoundaries degreeValue@(HomologicalDegree degreeIndex) = - let ambientDimension = cellCountAtDegree finite degreeValue - currentKernel = +homologyBasisAtPrepared preparedBoundaries degreeValue@(HomologicalDegree degreeIndex) = + let currentKernel = rationalBoundaryKernelBasis (rationalBoundaryAt preparedBoundaries degreeValue) incomingImage = rationalBoundaryImageBasis (rationalBoundaryAt preparedBoundaries (HomologicalDegree (degreeIndex + 1))) - in sparseIndependentModulo ambientDimension incomingImage currentKernel + in sparseIndependentModulo incomingImage currentKernel & fmap (sparseVectorToRepresentative degreeValue) cohomologyBasisAtPrepared :: - FiniteChainComplex r -> IntMap.IntMap RationalBoundaryDecomposition -> HomologicalDegree -> [RepresentativeCocycle Rational Int] -cohomologyBasisAtPrepared finite preparedBoundaries degreeValue@(HomologicalDegree degreeIndex) = - let ambientDimension = cellCountAtDegree finite degreeValue - currentKernel = +cohomologyBasisAtPrepared preparedBoundaries degreeValue@(HomologicalDegree degreeIndex) = + let currentKernel = rationalCoboundaryKernelBasis (rationalBoundaryAt preparedBoundaries (HomologicalDegree (degreeIndex + 1))) incomingImage = rationalCoboundaryImageBasis (rationalBoundaryAt preparedBoundaries degreeValue) - in sparseIndependentModulo ambientDimension incomingImage currentKernel + in sparseIndependentModulo incomingImage currentKernel & fmap (sparseVectorToRepresentative degreeValue) rationalBoundaryAt :: diff --git a/src-topology/Moonlight/Homology/Pure/Topology/SparseAlgebra.hs b/src-topology/Moonlight/Homology/Pure/Topology/SparseAlgebra.hs --- a/src-topology/Moonlight/Homology/Pure/Topology/SparseAlgebra.hs +++ b/src-topology/Moonlight/Homology/Pure/Topology/SparseAlgebra.hs @@ -10,6 +10,7 @@ import Data.IntMap.Strict qualified as IntMap import Data.Maybe (mapMaybe) import Data.Set qualified as Set +import Data.Vector qualified as Vector import Moonlight.Homology.Boundary.Finite (FiniteChainComplex, incidenceMatrixAt) import Moonlight.Homology.Pure.Chain ( HomologicalDegree (..), @@ -28,6 +29,7 @@ import Moonlight.Homology.Pure.Matrix.SparseLinAlg ( SparseMatrix (..), SparseRow, + sparseBoundaryColumns, sparseBoundaryMatrix, sparseIndependentModulo, sparseKernelBasisOf, @@ -75,7 +77,9 @@ degreeValue (cellCountAtDegree finite degreeValue) (sparseBoundaryMatrix (incidenceMatrixAt finite degreeValue)) - (sparseBoundaryMatrix (incidenceMatrixAt finite (HomologicalDegree (degreeIndex + 1)))) + ( Vector.toList + (sparseBoundaryColumns (incidenceMatrixAt finite (HomologicalDegree (degreeIndex + 1)))) + ) genericSparseCohomologyBasisAt :: Integral r => @@ -88,7 +92,7 @@ degreeValue (cellCountAtDegree finite degreeValue) (sparseTransposeMatrix (sparseBoundaryMatrix (incidenceMatrixAt finite (HomologicalDegree (degreeIndex + 1))))) - (sparseTransposeMatrix (sparseBoundaryMatrix (incidenceMatrixAt finite degreeValue))) + (smRows (sparseBoundaryMatrix (incidenceMatrixAt finite degreeValue))) graphHomologyZeroRepresentatives :: GraphOneComplex -> [RepresentativeCycle Rational Int] graphHomologyZeroRepresentatives graph = @@ -122,12 +126,11 @@ HomologicalDegree -> Int -> SparseMatrix -> - SparseMatrix -> + [SparseRow] -> [RepresentativeChain Rational Int] -sparseQuotientRepresentatives degreeValue ambientDimension currentMatrix incomingMatrix = +sparseQuotientRepresentatives degreeValue ambientDimension currentMatrix imageGenerators = let kernelBasis = sparseKernelBasisOf ambientDimension currentMatrix - imageGenerators = smRows (sparseTransposeMatrix incomingMatrix) - quotientBasis = sparseIndependentModulo ambientDimension imageGenerators kernelBasis + quotientBasis = sparseIndependentModulo imageGenerators kernelBasis in fmap (sparseVectorToRepresentative degreeValue) quotientBasis sparseVectorToRepresentative :: HomologicalDegree -> SparseRow -> RepresentativeChain Rational Int diff --git a/src-topology/Moonlight/Homology/Pure/Topology/Zigzag.hs b/src-topology/Moonlight/Homology/Pure/Topology/Zigzag.hs new file mode 100644 --- /dev/null +++ b/src-topology/Moonlight/Homology/Pure/Topology/Zigzag.hs @@ -0,0 +1,692 @@ +{-# LANGUAGE DeriveFoldable #-} +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE DeriveTraversable #-} +{-# LANGUAGE LambdaCase #-} + +-- | Exact persistence for finite, non-monotone diagrams of chain complexes. +-- +-- A 'FiniteChainMap' is admitted only after its component maps commute with +-- the two boundary operators. A 'FiniteChainZigzag' then glues checked maps +-- along structurally equal endpoints. Persistence is computed on rational +-- homology by propagating the right filtration of each type-A prefix; quotient +-- layers that fail descent close at that arrow, and the terminal layers are +-- the surviving intervals. +module Moonlight.Homology.Pure.Topology.Zigzag + ( ChainMapEndpoint (..), + ZigzagDirection (..), + ZigzagFailure (..), + FiniteChainMap, + finiteChainMapSource, + finiteChainMapTarget, + finiteChainMapAt, + mkFiniteChainMapChecked, + ZigzagArrow (..), + zigzagArrowDirection, + FiniteChainZigzag, + mkFiniteChainZigzag, + finiteChainZigzagComplexes, + finiteChainZigzagArrows, + ZigzagInterval (..), + rationalZigzagIntervals, + zigzagBettiAt, + ) +where + +import Control.Monad (foldM) +import Data.Bifunctor (first) +import Data.Foldable (traverse_) +import Data.Function ((&)) +import Data.IntMap.Strict qualified as IntMap +import Data.Kind (Type) +import Data.List qualified as List +import Data.List.NonEmpty (NonEmpty (..)) +import Data.Map.Strict (Map) +import Data.Map.Strict qualified as Map +import Data.Vector (Vector) +import Data.Vector qualified as Vector +import Moonlight.Core (Semiring) +import Moonlight.Homology.Boundary.Finite + ( FiniteChainComplex, + degreeCardinality, + incidenceMatrixAt, + maxHomologicalDegree, + validateFiniteChainComplexShape, + ) +import Moonlight.Homology.Boundary.LinAlg + ( BoundaryIncidence, + BoundaryIncidenceShapeError, + composeBoundaryIncidence, + emptyBoundaryIncidenceOf, + sourceCardinality, + targetCardinality, + ) +import Moonlight.Homology.Pure.Chain + ( HomologicalDegree (..), + RepresentativeChain (..), + ) +import Moonlight.Homology.Pure.Failure (HomologyFailure) +import Moonlight.Homology.Pure.Matrix.SparseLinAlg + ( SparseCoordinateBasis, + SparseColumnEchelon (..), + SparseRow, + compactSparseRow, + sparseBoundaryColumns, + sparseCoordinateBasis, + sparseCoordinatesInBasis, + sparseColumnEchelon, + sparseEchelonBasis, + sparseExtendEchelonBasis, + sparseLinearCombination, + ) +import Moonlight.Homology.Pure.Topology.SparseAlgebra + ( sparseHomologyBasisAt, + ) +import Numeric.Natural (Natural) + +type ChainMapEndpoint :: Type +data ChainMapEndpoint + = ChainMapSource + | ChainMapTarget + deriving stock (Eq, Ord, Show) + +type ZigzagDirection :: Type +data ZigzagDirection + = ZigzagForward + | ZigzagBackward + deriving stock (Eq, Ord, Show) + +-- | Every refusal names the exact descent obligation that failed. In +-- particular, absence of a persistence class is a successful empty result, +-- never one of these failures. +type ZigzagFailure :: Type +data ZigzagFailure + = ZigzagComplexInvalid !ChainMapEndpoint !HomologyFailure + | ZigzagMapSourceCardinalityMismatch !HomologicalDegree !Int !Int + | ZigzagMapTargetCardinalityMismatch !HomologicalDegree !Int !Int + | ZigzagMapComponentInvalid !HomologicalDegree !BoundaryIncidenceShapeError + | ZigzagMapCompositionInvalid !HomologicalDegree !BoundaryIncidenceShapeError + | ZigzagChainMapLawViolation !HomologicalDegree + | ZigzagEndpointMismatch !Int !ZigzagDirection + | ZigzagHomologyCoordinatesMissing !Int !HomologicalDegree + | ZigzagNegativeIntervalMultiplicity !HomologicalDegree !Int !Int !Int + deriving stock (Eq, Show) + +-- | A checked degree-preserving chain map. Its source, target, and materialized +-- degree components are retained together so no caller can later pair the map +-- with different complexes. +type FiniteChainMap :: Type -> Type +data FiniteChainMap r = FiniteChainMap + { storedChainMapSource :: !(FiniteChainComplex r), + storedChainMapTarget :: !(FiniteChainComplex r), + storedChainMapComponents :: !(Vector (BoundaryIncidence r)) + } + +finiteChainMapSource :: FiniteChainMap r -> FiniteChainComplex r +finiteChainMapSource = storedChainMapSource + +finiteChainMapTarget :: FiniteChainMap r -> FiniteChainComplex r +finiteChainMapTarget = storedChainMapTarget + +finiteChainMapAt :: FiniteChainMap r -> HomologicalDegree -> BoundaryIncidence r +finiteChainMapAt chainMap (HomologicalDegree degreeIndex) = + if degreeIndex < 0 + then emptyBoundaryIncidenceOf 0 0 + else + case storedChainMapComponents chainMap Vector.!? degreeIndex of + Just component -> component + Nothing -> emptyBoundaryIncidenceOf 0 0 + +-- | Admit a finite chain map after materializing its relevant degrees and +-- proving @d_target . f = f . d_source@ at every positive degree. +mkFiniteChainMapChecked :: + (Eq r, Num r, Semiring r) => + FiniteChainComplex r -> + FiniteChainComplex r -> + (HomologicalDegree -> BoundaryIncidence r) -> + Either ZigzagFailure (FiniteChainMap r) +mkFiniteChainMapChecked sourceComplex targetComplex componentAt = do + first (ZigzagComplexInvalid ChainMapSource) (validateFiniteChainComplexShape sourceComplex) + first (ZigzagComplexInvalid ChainMapTarget) (validateFiniteChainComplexShape targetComplex) + let maximumDegree = max (maximumDegreeOf sourceComplex) (maximumDegreeOf targetComplex) + degreeValues = fmap HomologicalDegree [0 .. maximumDegree] + componentList = fmap componentAt degreeValues + components = Vector.fromList componentList + traverse_ + (uncurry (validateComponentShape sourceComplex targetComplex)) + (zip degreeValues componentList) + traverse_ + (\(degreeValue, precedingComponent, degreeComponent) -> + validateChainMapLaw sourceComplex targetComplex degreeValue precedingComponent degreeComponent + ) + (zip3 (drop 1 degreeValues) componentList (drop 1 componentList)) + pure + FiniteChainMap + { storedChainMapSource = sourceComplex, + storedChainMapTarget = targetComplex, + storedChainMapComponents = components + } + +validateComponentShape :: + FiniteChainComplex r -> + FiniteChainComplex r -> + HomologicalDegree -> + BoundaryIncidence r -> + Either ZigzagFailure () +validateComponentShape sourceComplex targetComplex degreeValue component + | sourceCardinality component /= expectedSource = + Left (ZigzagMapSourceCardinalityMismatch degreeValue expectedSource (sourceCardinality component)) + | targetCardinality component /= expectedTarget = + Left (ZigzagMapTargetCardinalityMismatch degreeValue expectedTarget (targetCardinality component)) + | otherwise = Right () + where + expectedSource = degreeCardinality sourceComplex degreeValue + expectedTarget = degreeCardinality targetComplex degreeValue + +validateChainMapLaw :: + (Eq r, Num r, Semiring r) => + FiniteChainComplex r -> + FiniteChainComplex r -> + HomologicalDegree -> + BoundaryIncidence r -> + BoundaryIncidence r -> + Either ZigzagFailure () +validateChainMapLaw sourceComplex targetComplex degreeValue precedingComponent degreeComponent = do + targetAfterMap <- + first (ZigzagMapCompositionInvalid degreeValue) + ( composeBoundaryIncidence + (finiteBoundaryAt targetComplex degreeValue) + degreeComponent + ) + mapAfterSource <- + first (ZigzagMapCompositionInvalid degreeValue) + ( composeBoundaryIncidence + precedingComponent + (finiteBoundaryAt sourceComplex degreeValue) + ) + if targetAfterMap == mapAfterSource + then Right () + else Left (ZigzagChainMapLawViolation degreeValue) + +-- | Orientation separated from its payload. The same functor carries checked +-- chain maps during authoring and exact linear maps during reduction. +type ZigzagArrow :: Type -> Type +data ZigzagArrow map + = ForwardArrow !map + | BackwardArrow !map + deriving stock (Eq, Ord, Show, Functor, Foldable, Traversable) + +zigzagArrowDirection :: ZigzagArrow r -> ZigzagDirection +zigzagArrowDirection = \case + ForwardArrow _ -> ZigzagForward + BackwardArrow _ -> ZigzagBackward + +type FiniteChainZigzag :: Type -> Type +data FiniteChainZigzag r = FiniteChainZigzag + { storedZigzagFirstComplex :: !(FiniteChainComplex r), + storedZigzagArrows :: !(Vector (ZigzagArrow (FiniteChainMap r))) + } + +-- | Glue a line of already checked maps. A forward arrow is interpreted as +-- @current -> next@; a backward arrow is @current <- next@. +mkFiniteChainZigzag :: + Eq r => + FiniteChainComplex r -> + [ZigzagArrow (FiniteChainMap r)] -> + Either ZigzagFailure (FiniteChainZigzag r) +mkFiniteChainZigzag firstComplex arrows = do + _ <- foldM glueZigzagArrow firstComplex (zip [0 :: Int ..] arrows) + pure + FiniteChainZigzag + { storedZigzagFirstComplex = firstComplex, + storedZigzagArrows = Vector.fromList arrows + } + +glueZigzagArrow :: + Eq r => + FiniteChainComplex r -> + (Int, ZigzagArrow (FiniteChainMap r)) -> + Either ZigzagFailure (FiniteChainComplex r) +glueZigzagArrow currentComplex (arrowIndex, arrow) = + let (requiredCurrent, nextComplex) = finiteArrowEndpoints arrow + in if finiteComplexesAgree currentComplex requiredCurrent + then Right nextComplex + else Left (ZigzagEndpointMismatch arrowIndex (zigzagArrowDirection arrow)) + +finiteArrowEndpoints :: ZigzagArrow (FiniteChainMap r) -> (FiniteChainComplex r, FiniteChainComplex r) +finiteArrowEndpoints = \case + ForwardArrow chainMap -> (finiteChainMapSource chainMap, finiteChainMapTarget chainMap) + BackwardArrow chainMap -> (finiteChainMapTarget chainMap, finiteChainMapSource chainMap) + +finiteChainZigzagComplexes :: FiniteChainZigzag r -> NonEmpty (FiniteChainComplex r) +finiteChainZigzagComplexes zigzag = + storedZigzagFirstComplex zigzag + :| Vector.toList + (fmap (snd . finiteArrowEndpoints) (storedZigzagArrows zigzag)) + +finiteChainZigzagArrows :: FiniteChainZigzag r -> [ZigzagArrow (FiniteChainMap r)] +finiteChainZigzagArrows = Vector.toList . storedZigzagArrows + +zigzagComplexVector :: FiniteChainZigzag r -> Vector (FiniteChainComplex r) +zigzagComplexVector zigzag = + Vector.cons + (storedZigzagFirstComplex zigzag) + (fmap (snd . finiteArrowEndpoints) (storedZigzagArrows zigzag)) + +-- | One indecomposable interval of a zigzag barcode. Both endpoints are +-- inclusive. 'Traversable' transports the same interval from numeric diagram +-- indices to a caller's stage vocabulary without defining a parallel carrier. +type ZigzagInterval :: Type -> Type +data ZigzagInterval endpoint = ZigzagInterval + { zigzagIntervalDegree :: !HomologicalDegree, + zigzagIntervalFirst :: !endpoint, + zigzagIntervalLast :: !endpoint, + zigzagIntervalMultiplicity :: !Int + } + deriving stock (Eq, Ord, Show, Functor, Foldable, Traversable) + +-- | Compute the unique interval decomposition of the rational homology +-- zigzag. Chain homology and every induced adjacent map are prepared once per +-- degree; one right-filtration descent then closes every interval. +rationalZigzagIntervals :: + Integral r => + FiniteChainZigzag r -> + Either ZigzagFailure [ZigzagInterval Int] +rationalZigzagIntervals zigzag = + let complexes = zigzagComplexVector zigzag + arrows = storedZigzagArrows zigzag + maximumDegree = + complexes + & Vector.foldl' + (\currentMaximum complexValue -> max currentMaximum (maximumDegreeOf complexValue)) + 0 + in concat + <$> traverse + ( intervalsAtDegree + (storedZigzagFirstComplex zigzag) + arrows + . HomologicalDegree + ) + [0 .. maximumDegree] + +-- | Betti numbers reconstructed from a barcode at one diagram index. The map +-- is sparse; absent degrees have dimension zero. +zigzagBettiAt :: Int -> [ZigzagInterval Int] -> Map HomologicalDegree Int +zigzagBettiAt diagramIndex = + Map.fromListWith (+) + . fmap + (\interval -> (zigzagIntervalDegree interval, zigzagIntervalMultiplicity interval)) + . filter + ( \interval -> + zigzagIntervalFirst interval <= diagramIndex + && diagramIndex <= zigzagIntervalLast interval + ) + +type HomologyPresentation :: Type +data HomologyPresentation = HomologyPresentation + { homologyBasisVectors :: !(Vector SparseRow), + homologyCoordinateBasis :: !SparseCoordinateBasis + } + +homologyDimension :: HomologyPresentation -> Int +homologyDimension = Vector.length . homologyBasisVectors + +type RationalLinearMap :: Type +data RationalLinearMap = RationalLinearMap + { rationalMapTargetDimension :: !Int, + rationalMapColumns :: !(Vector SparseRow) + } + +-- | One quotient layer of the right filtration on the current endpoint. The +-- vectors form a basis for that layer modulo every preceding layer; the layer +-- order, not numeric birth order, is the zigzag orientation witness. +type RightFiltrationLayer :: Type +data RightFiltrationLayer = RightFiltrationLayer + { rightLayerBirthIndex :: !Int, + rightLayerBasis :: ![SparseRow] + } + +type IntervalMultiplicities :: Type +type IntervalMultiplicities = Map (Int, Int) Int + +intervalsAtDegree :: + Integral r => + FiniteChainComplex r -> + Vector (ZigzagArrow (FiniteChainMap r)) -> + HomologicalDegree -> + Either ZigzagFailure [ZigzagInterval Int] +intervalsAtDegree firstComplex arrows degreeValue = do + let initialPresentation = homologyPresentationAt firstComplex degreeValue + initialFiltration = + [ RightFiltrationLayer + { rightLayerBirthIndex = 0, + rightLayerBasis = standardSparseBasis (homologyDimension initialPresentation) + } + ] + -- Prepare, induce, and descend one arrow at a time: only the adjacent + -- presentations and current filtration remain live. + (_, terminalFiltration, closedIntervals) <- + Vector.ifoldM' + ( \(leftPresentation, currentFiltration, intervals) arrowIndex arrow -> do + let rightPresentation = + homologyPresentationAt + (snd (finiteArrowEndpoints arrow)) + degreeValue + inducedArrow <- + inducedArrowAt + degreeValue + arrowIndex + (leftPresentation, rightPresentation, arrow) + (nextFiltration, nextIntervals) <- + advanceRightFiltration + degreeValue + (currentFiltration, intervals) + (arrowIndex, inducedArrow) + pure (rightPresentation, nextFiltration, nextIntervals) + ) + (initialPresentation, initialFiltration, Map.empty) + arrows + terminalIntervals <- + closeTerminalIntervals + degreeValue + (Vector.length arrows) + terminalFiltration + closedIntervals + pure (intervalsFromMultiplicities degreeValue terminalIntervals) + +homologyPresentationAt :: + Integral r => + FiniteChainComplex r -> + HomologicalDegree -> + HomologyPresentation +homologyPresentationAt finite degreeValue@(HomologicalDegree degreeIndex) + | degreeIndex < 0 || degreeIndex > maximumDegreeOf finite = + HomologyPresentation + { homologyBasisVectors = Vector.empty, + homologyCoordinateBasis = sparseCoordinateBasis [] + } + | otherwise = + let basisRows = + sparseHomologyBasisAt finite degreeValue + & fmap representativeSparseRow + basisVectors = Vector.fromList basisRows + boundaryGenerators = + finiteBoundaryAt finite (HomologicalDegree (degreeIndex + 1)) + & sparseBoundaryColumns + in HomologyPresentation + { homologyBasisVectors = basisVectors, + homologyCoordinateBasis = + sparseCoordinateBasis + (basisRows <> Vector.toList boundaryGenerators) + } + +inducedArrowAt :: + Integral r => + HomologicalDegree -> + Int -> + (HomologyPresentation, HomologyPresentation, ZigzagArrow (FiniteChainMap r)) -> + Either ZigzagFailure (ZigzagArrow RationalLinearMap) +inducedArrowAt degreeValue arrowIndex (leftPresentation, rightPresentation, arrow) = + case arrow of + ForwardArrow chainMap -> + ForwardArrow + <$> inducedHomologyMap arrowIndex degreeValue chainMap leftPresentation rightPresentation + BackwardArrow chainMap -> + BackwardArrow + <$> inducedHomologyMap arrowIndex degreeValue chainMap rightPresentation leftPresentation + +inducedHomologyMap :: + Integral r => + Int -> + HomologicalDegree -> + FiniteChainMap r -> + HomologyPresentation -> + HomologyPresentation -> + Either ZigzagFailure RationalLinearMap +inducedHomologyMap arrowIndex degreeValue chainMap sourcePresentation targetPresentation = do + let rationalColumns = + sparseBoundaryColumns (finiteChainMapAt chainMap degreeValue) + imageCoordinates <- + traverse + ( \sourceCycle -> do + let mappedCycle = sparseLinearCombination rationalColumns sourceCycle + coordinates <- + maybe + (Left (ZigzagHomologyCoordinatesMissing arrowIndex degreeValue)) + Right + ( sparseCoordinatesInBasis + (homologyCoordinateBasis targetPresentation) + mappedCycle + ) + pure + ( IntMap.filterWithKey + (\coordinateIndex _ -> coordinateIndex < homologyDimension targetPresentation) + coordinates + ) + ) + (homologyBasisVectors sourcePresentation) + pure + RationalLinearMap + { rationalMapTargetDimension = homologyDimension targetPresentation, + rationalMapColumns = imageCoordinates + } + +applyRationalLinearMap :: RationalLinearMap -> SparseRow -> SparseRow +applyRationalLinearMap = sparseLinearCombination . rationalMapColumns + +advanceRightFiltration :: + HomologicalDegree -> + ([RightFiltrationLayer], IntervalMultiplicities) -> + (Int, ZigzagArrow RationalLinearMap) -> + Either ZigzagFailure ([RightFiltrationLayer], IntervalMultiplicities) +advanceRightFiltration degreeValue (currentFiltration, intervals) (arrowIndex, arrow) = do + let nextBirthIndex = arrowIndex + 1 + (nextFiltration, survivingLayers) <- + case arrow of + ForwardArrow linearMap -> + Right + ( forwardRightFiltration + nextBirthIndex + linearMap + currentFiltration + ) + BackwardArrow linearMap -> + backwardRightFiltration + arrowIndex + degreeValue + nextBirthIndex + linearMap + currentFiltration + updatedIntervals <- + closeExpiredIntervals + degreeValue + arrowIndex + currentFiltration + survivingLayers + intervals + pure (nextFiltration, updatedIntervals) + +forwardRightFiltration :: + Int -> + RationalLinearMap -> + [RightFiltrationLayer] -> + ([RightFiltrationLayer], [RightFiltrationLayer]) +forwardRightFiltration nextBirthIndex linearMap currentFiltration = + let targetDimension = rationalMapTargetDimension linearMap + (reversedSurvivingLayers, imageBasis) = + List.foldl' + ( \(reversedLayers, accumulatedBasis) layer -> + let mappedVectors = fmap (applyRationalLinearMap linearMap) (rightLayerBasis layer) + (independentImages, extendedBasis) = + sparseExtendEchelonBasis accumulatedBasis mappedVectors + in ( layer {rightLayerBasis = independentImages} : reversedLayers, + extendedBasis + ) + ) + ([], sparseEchelonBasis []) + currentFiltration + survivingLayers = reverse reversedSurvivingLayers + (newLayerBasis, _) = + sparseExtendEchelonBasis imageBasis (standardSparseBasis targetDimension) + nextFiltration = + survivingLayers + <> [ RightFiltrationLayer + { rightLayerBirthIndex = nextBirthIndex, + rightLayerBasis = newLayerBasis + } + ] + in (nextFiltration, survivingLayers) + +backwardRightFiltration :: + Int -> + HomologicalDegree -> + Int -> + RationalLinearMap -> + [RightFiltrationLayer] -> + Either ZigzagFailure ([RightFiltrationLayer], [RightFiltrationLayer]) +backwardRightFiltration arrowIndex degreeValue nextBirthIndex linearMap currentFiltration = do + let targetDimension = rationalMapTargetDimension linearMap + targetFiltrationBasis = currentFiltration >>= rightLayerBasis + targetCoordinates = sparseCoordinateBasis targetFiltrationBasis + imageCoordinateColumns <- + traverse + ( \imageColumn -> + maybe + (Left (ZigzagHomologyCoordinatesMissing arrowIndex degreeValue)) + Right + (sparseCoordinatesInBasis targetCoordinates imageColumn) + ) + (rationalMapColumns linearMap) + let columnEchelon = + sparseColumnEchelon + (fmap (reverseSparseCoordinates targetDimension) imageCoordinateColumns) + ascendingPivotPreimages = + sparseColumnPivotPreimages columnEchelon + & fmap + (\(reversedPivot, preimage) -> (targetDimension - reversedPivot - 1, preimage)) + & reverse + (_, survivingLayers) = + List.mapAccumL + pullbackLayer + (0, ascendingPivotPreimages) + currentFiltration + kernelLayer = + RightFiltrationLayer + { rightLayerBirthIndex = nextBirthIndex, + rightLayerBasis = sparseColumnKernelBasis columnEchelon + } + pure (kernelLayer : survivingLayers, survivingLayers) + where + pullbackLayer (lowerBound, remainingPivots) layer = + let upperBound = lowerBound + length (rightLayerBasis layer) + (layerPivots, laterPivots) = + span ((< upperBound) . fst) remainingPivots + in ( (upperBound, laterPivots), + layer {rightLayerBasis = fmap snd layerPivots} + ) + +reverseSparseCoordinates :: Int -> SparseRow -> SparseRow +reverseSparseCoordinates dimensionValue = + IntMap.fromDistinctAscList + . fmap (\(coordinateIndex, coefficient) -> (dimensionValue - coordinateIndex - 1, coefficient)) + . IntMap.toDescList + +closeExpiredIntervals :: + HomologicalDegree -> + Int -> + [RightFiltrationLayer] -> + [RightFiltrationLayer] -> + IntervalMultiplicities -> + Either ZigzagFailure IntervalMultiplicities +closeExpiredIntervals degreeValue deathIndex currentLayers survivingLayers intervals = + foldM + ( \currentIntervals (currentLayer, survivingLayer) -> + recordIntervalMultiplicity + degreeValue + (rightLayerBirthIndex currentLayer) + deathIndex + (length (rightLayerBasis currentLayer) - length (rightLayerBasis survivingLayer)) + currentIntervals + ) + intervals + (zip currentLayers survivingLayers) + +closeTerminalIntervals :: + HomologicalDegree -> + Int -> + [RightFiltrationLayer] -> + IntervalMultiplicities -> + Either ZigzagFailure IntervalMultiplicities +closeTerminalIntervals degreeValue deathIndex terminalFiltration intervals = + foldM + ( \currentIntervals layer -> + recordIntervalMultiplicity + degreeValue + (rightLayerBirthIndex layer) + deathIndex + (length (rightLayerBasis layer)) + currentIntervals + ) + intervals + terminalFiltration + +recordIntervalMultiplicity :: + HomologicalDegree -> + Int -> + Int -> + Int -> + IntervalMultiplicities -> + Either ZigzagFailure IntervalMultiplicities +recordIntervalMultiplicity degreeValue firstIndex lastIndex multiplicity intervals + | multiplicity < 0 = + Left (ZigzagNegativeIntervalMultiplicity degreeValue firstIndex lastIndex multiplicity) + | multiplicity == 0 = Right intervals + | otherwise = Right (Map.insertWith (+) (firstIndex, lastIndex) multiplicity intervals) + +intervalsFromMultiplicities :: HomologicalDegree -> IntervalMultiplicities -> [ZigzagInterval Int] +intervalsFromMultiplicities degreeValue = + fmap + ( \((firstIndex, lastIndex), multiplicity) -> + ZigzagInterval + { zigzagIntervalDegree = degreeValue, + zigzagIntervalFirst = firstIndex, + zigzagIntervalLast = lastIndex, + zigzagIntervalMultiplicity = multiplicity + } + ) + . Map.toAscList + +standardSparseBasis :: Int -> [SparseRow] +standardSparseBasis dimensionValue = + fmap (\coordinateIndex -> IntMap.singleton coordinateIndex 1) [0 .. dimensionValue - 1] + +representativeSparseRow :: RepresentativeChain Rational Int -> SparseRow +representativeSparseRow representative = + representativeTerms representative + & fmap (\(coefficient, basisIndex) -> (basisIndex, coefficient)) + & IntMap.fromListWith (+) + & compactSparseRow + +finiteComplexesAgree :: Eq r => FiniteChainComplex r -> FiniteChainComplex r -> Bool +finiteComplexesAgree left right = + let maximumDegree = max (maximumDegreeOf left) (maximumDegreeOf right) + in all + ( \degreeIndex -> + finiteBoundaryAt left (HomologicalDegree degreeIndex) + == finiteBoundaryAt right (HomologicalDegree degreeIndex) + ) + [0 .. maximumDegree] + +finiteBoundaryAt :: FiniteChainComplex r -> HomologicalDegree -> BoundaryIncidence r +finiteBoundaryAt finite degreeValue@(HomologicalDegree degreeIndex) = + if degreeIndex >= 0 && degreeIndex <= maximumDegreeOf finite + then incidenceMatrixAt finite degreeValue + else + emptyBoundaryIncidenceOf + (naturalCardinality (degreeCardinality finite degreeValue)) + (naturalCardinality (degreeCardinality finite (HomologicalDegree (degreeIndex - 1)))) + +maximumDegreeOf :: FiniteChainComplex r -> Int +maximumDegreeOf = unHomologicalDegree . maxHomologicalDegree + +naturalCardinality :: Int -> Natural +naturalCardinality = fromIntegral . max 0 diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -28,6 +28,7 @@ triangleCycleComplex, ) import TopologySpec qualified +import ZigzagSpec qualified import Test.Tasty (TestTree, defaultMain, testGroup) import Test.Tasty.HUnit (assertBool, assertEqual, assertFailure, testCase) import qualified Test.Tasty.QuickCheck as QC @@ -57,7 +58,8 @@ determinismHarnessTests, MorseSpec.tests, PresentationSpec.tests, - TopologySpec.tests + TopologySpec.tests, + ZigzagSpec.tests ] emptyIntFiniteComplex :: FiniteChainComplex Int 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.2-inplace"] + [GhcPackageId "moonlight-homology-0.1.0.3-inplace"] expectRight :: Show left => diff --git a/test/sequence/SpectralSpec.hs b/test/sequence/SpectralSpec.hs --- a/test/sequence/SpectralSpec.hs +++ b/test/sequence/SpectralSpec.hs @@ -639,11 +639,7 @@ && all ((== expectedColumnCount) . length) matrixValue matrixRank :: [[Rational]] -> Int -matrixRank matrixValue = - case matrixValue of - [] -> 0 - firstRow : _ -> - sparseSpanRank (length firstRow) (fmap sparseRowFromDense matrixValue) +matrixRank = sparseSpanRank . fmap sparseRowFromDense zeroMatrixValue :: [[Rational]] -> Bool zeroMatrixValue = diff --git a/test/topology/ZigzagSpec.hs b/test/topology/ZigzagSpec.hs new file mode 100644 --- /dev/null +++ b/test/topology/ZigzagSpec.hs @@ -0,0 +1,388 @@ +module ZigzagSpec (tests) where + +import Data.Bifunctor (first) +import Data.Foldable (traverse_) +import Data.List qualified as List +import Data.Map.Strict qualified as Map +import Moonlight.Homology.Boundary + ( BoundaryIncidence, + BoundaryIncidenceShapeError, + FiniteChainComplex, + degreeCardinality, + emptyBoundaryIncidence, + emptyBoundaryIncidenceOf, + mkBoundaryEntry, + mkBoundaryIncidence, + ) +import Moonlight.Homology.Boundary.Finite (mkFiniteChainComplex) +import Moonlight.Homology.Chain (HomologicalDegree (..)) +import Moonlight.Homology.Persistence + ( FiniteChainMap, + ZigzagArrow (..), + ZigzagDirection (..), + ZigzagFailure (..), + ZigzagInterval (..), + mkFiniteChainMapChecked, + mkFiniteChainZigzag, + rationalZigzagIntervals, + zigzagBettiAt, + ) +import Test.Tasty (TestTree, testGroup) +import Test.Tasty.HUnit (Assertion, assertFailure, testCase, (@?=)) + +tests :: TestTree +tests = + testGroup + "exact zigzag persistence" + [ testCase "decomposes a singleton vector space" singletonIntervals, + testCase "preserves one class through a forward identity" forwardIdentityInterval, + testCase "separates classes across a forward zero map" forwardZeroIntervals, + testCase "separates classes across a backward zero map" backwardZeroIntervals, + testCase "decomposes a non-monotone cospan" nonMonotoneCospanIntervals, + testCase "pulls a flag backward through exact cancellation" backwardPullbackCancellation, + testCase "recovers interval sums under every arrow orientation" everyDirectionIntervalSum, + testCase "induces the identity on degree-one homology" degreeOneIdentityInterval, + testCase "treats degrees above a complex maximum as zero" mixedMaximumDegrees, + testCase "reconstructs every vertex Betti number" reconstructVertexBetti, + testCase "rejects a component with the wrong source dimension" rejectMapShape, + testCase "rejects a degree map that does not commute with boundaries" rejectChainMapLaw, + testCase "rejects diagram endpoint mismatches" rejectEndpointMismatch + ] + +singletonIntervals :: Assertion +singletonIntervals = do + diagram <- requireRight "singleton diagram" (mkFiniteChainZigzag (zeroComplex 2) []) + intervals <- requireRight "singleton intervals" (rationalZigzagIntervals diagram) + intervals @?= [ZigzagInterval (HomologicalDegree 0) 0 0 2] + +forwardIdentityInterval :: Assertion +forwardIdentityInterval = do + identityMap <- requireRight "identity map" (coordinateMap (zeroComplex 1) (zeroComplex 1) [(0, 0)]) + diagram <- requireRight "forward diagram" (mkFiniteChainZigzag (zeroComplex 1) [ForwardArrow identityMap]) + intervals <- requireRight "forward intervals" (rationalZigzagIntervals diagram) + intervals @?= [ZigzagInterval (HomologicalDegree 0) 0 1 1] + +forwardZeroIntervals :: Assertion +forwardZeroIntervals = do + zeroMap <- requireRight "forward zero map" (coordinateMap (zeroComplex 1) (zeroComplex 1) []) + diagram <- requireRight "forward zero diagram" (mkFiniteChainZigzag (zeroComplex 1) [ForwardArrow zeroMap]) + intervals <- requireRight "forward zero intervals" (rationalZigzagIntervals diagram) + intervals + @?= [ ZigzagInterval (HomologicalDegree 0) 0 0 1, + ZigzagInterval (HomologicalDegree 0) 1 1 1 + ] + +backwardZeroIntervals :: Assertion +backwardZeroIntervals = do + zeroMap <- requireRight "backward zero map" (coordinateMap (zeroComplex 1) (zeroComplex 1) []) + diagram <- requireRight "backward zero diagram" (mkFiniteChainZigzag (zeroComplex 1) [BackwardArrow zeroMap]) + intervals <- requireRight "backward zero intervals" (rationalZigzagIntervals diagram) + intervals + @?= [ ZigzagInterval (HomologicalDegree 0) 0 0 1, + ZigzagInterval (HomologicalDegree 0) 1 1 1 + ] + +nonMonotoneCospanIntervals :: Assertion +nonMonotoneCospanIntervals = do + leftInclusion <- requireRight "left inclusion" (coordinateMap (zeroComplex 1) (zeroComplex 2) [(0, 0)]) + rightInclusion <- requireRight "right inclusion" (coordinateMap (zeroComplex 1) (zeroComplex 2) [(0, 1)]) + diagram <- + requireRight + "cospan diagram" + ( mkFiniteChainZigzag + (zeroComplex 1) + [ForwardArrow leftInclusion, BackwardArrow rightInclusion] + ) + intervals <- requireRight "cospan intervals" (rationalZigzagIntervals diagram) + intervals + @?= [ ZigzagInterval (HomologicalDegree 0) 0 1 1, + ZigzagInterval (HomologicalDegree 0) 1 2 1 + ] + +backwardPullbackCancellation :: Assertion +backwardPullbackCancellation = do + firstMap <- + requireRight + "filtered inclusion" + (coordinateMap (zeroComplex 1) (zeroComplex 2) [(0, 0)]) + cancellationMap <- + requireRight + "cancelling backward map" + (coordinateMap (zeroComplex 2) (zeroComplex 2) [(0, 0), (0, 1), (1, 1)]) + diagram <- + requireRight + "backward cancellation diagram" + ( mkFiniteChainZigzag + (zeroComplex 1) + [ForwardArrow firstMap, BackwardArrow cancellationMap] + ) + intervals <- requireRight "backward cancellation intervals" (rationalZigzagIntervals diagram) + intervals + @?= [ ZigzagInterval (HomologicalDegree 0) 0 2 1, + ZigzagInterval (HomologicalDegree 0) 1 2 1 + ] + +everyDirectionIntervalSum :: Assertion +everyDirectionIntervalSum = do + let intervalSeeds :: [(Int, Int, Int)] + intervalSeeds = + zipWith + (\seedIndex (firstIndex, lastIndex) -> (seedIndex, firstIndex, lastIndex)) + [0 :: Int ..] + [(0, 3), (0, 1), (1, 2), (2, 3), (1, 1), (1, 1)] + stageSeeds stageIndex = + filter + (\(_, firstIndex, lastIndex) -> firstIndex <= stageIndex && stageIndex <= lastIndex) + intervalSeeds + stageComplex stageIndex = zeroComplex (length (stageSeeds stageIndex)) + coordinates sourceIndexValue targetIndexValue = + [ (sourceCoordinate, targetCoordinate) + | (sourceCoordinate, seed) <- zip [0 :: Int ..] (stageSeeds sourceIndexValue) + , Just targetCoordinate <- [List.elemIndex seed (stageSeeds targetIndexValue)] + ] + arrowAt arrowIndex direction = + let leftIndex = arrowIndex + rightIndex = arrowIndex + 1 + in case direction of + ZigzagForward -> + ForwardArrow + <$> coordinateMap + (stageComplex leftIndex) + (stageComplex rightIndex) + (coordinates leftIndex rightIndex) + ZigzagBackward -> + BackwardArrow + <$> coordinateMap + (stageComplex rightIndex) + (stageComplex leftIndex) + (coordinates rightIndex leftIndex) + expected = + [ ZigzagInterval (HomologicalDegree 0) 0 1 1 + , ZigzagInterval (HomologicalDegree 0) 0 3 1 + , ZigzagInterval (HomologicalDegree 0) 1 1 2 + , ZigzagInterval (HomologicalDegree 0) 1 2 1 + , ZigzagInterval (HomologicalDegree 0) 2 3 1 + ] + assertOrientation directions = do + arrows <- requireRight "oriented interval-sum maps" (traverse (uncurry arrowAt) (zip [0 ..] directions)) + diagram <- requireRight "oriented interval-sum diagram" (mkFiniteChainZigzag (stageComplex 0) arrows) + intervals <- requireRight "oriented interval-sum barcode" (rationalZigzagIntervals diagram) + intervals @?= expected + traverse_ assertOrientation (sequence (replicate 3 [ZigzagForward, ZigzagBackward])) + +degreeOneIdentityInterval :: Assertion +degreeOneIdentityInterval = do + complexValue <- requireRight "circle complex" circleComplex + degreeZeroIdentity <- requireRight "circle vertex identity" (identityIncidence 3) + degreeOneIdentity <- requireRight "circle edge identity" (identityIncidence 3) + identityMap <- + requireRight + "circle chain identity" + ( mkFiniteChainMapChecked complexValue complexValue $ \(HomologicalDegree degreeIndex) -> + case degreeIndex of + 0 -> degreeZeroIdentity + 1 -> degreeOneIdentity + _ -> emptyBoundaryIncidence + ) + diagram <- + requireRight + "circle identity diagram" + (mkFiniteChainZigzag complexValue [ForwardArrow identityMap]) + intervals <- requireRight "circle identity intervals" (rationalZigzagIntervals diagram) + intervals + @?= [ ZigzagInterval (HomologicalDegree 0) 0 1 1 + , ZigzagInterval (HomologicalDegree 1) 0 1 1 + ] + +mixedMaximumDegrees :: Assertion +mixedMaximumDegrees = do + targetComplex <- requireRight "mixed-maximum circle" circleComplex + degreeZeroInclusion <- + requireRight + "mixed-maximum degree-zero inclusion" + ( first (ZigzagMapComponentInvalid (HomologicalDegree 0)) + (mkBoundaryIncidence 1 3 [mkBoundaryEntry 0 0 (1 :: Int)]) + ) + let degreeOneInclusion :: BoundaryIncidence Int + degreeOneInclusion = emptyBoundaryIncidenceOf 0 3 + inclusion <- + requireRight + "mixed-maximum chain inclusion" + ( mkFiniteChainMapChecked poisonedPointComplex targetComplex $ \(HomologicalDegree degreeIndex) -> + case degreeIndex of + 0 -> degreeZeroInclusion + 1 -> degreeOneInclusion + _ -> emptyBoundaryIncidence + ) + diagram <- + requireRight + "mixed-maximum diagram" + (mkFiniteChainZigzag poisonedPointComplex [ForwardArrow inclusion]) + intervals <- requireRight "mixed-maximum intervals" (rationalZigzagIntervals diagram) + intervals + @?= [ ZigzagInterval (HomologicalDegree 0) 0 1 1 + , ZigzagInterval (HomologicalDegree 1) 1 1 1 + ] + +reconstructVertexBetti :: Assertion +reconstructVertexBetti = do + leftInclusion <- requireRight "shared left inclusion" (coordinateMap (zeroComplex 1) (zeroComplex 2) [(0, 0)]) + rightInclusion <- requireRight "shared right inclusion" (coordinateMap (zeroComplex 1) (zeroComplex 2) [(0, 0)]) + diagram <- + requireRight + "shared cospan diagram" + ( mkFiniteChainZigzag + (zeroComplex 1) + [ForwardArrow leftInclusion, BackwardArrow rightInclusion] + ) + intervals <- requireRight "shared cospan intervals" (rationalZigzagIntervals diagram) + intervals + @?= [ ZigzagInterval (HomologicalDegree 0) 0 2 1, + ZigzagInterval (HomologicalDegree 0) 1 1 1 + ] + fmap (`zigzagBettiAt` intervals) [0, 1, 2] + @?= fmap (Map.singleton (HomologicalDegree 0)) [1, 2, 1] + +rejectMapShape :: Assertion +rejectMapShape = + case + mkFiniteChainMapChecked + (zeroComplex 1) + (zeroComplex 1) + (const (emptyBoundaryIncidenceOf 2 1)) + of + Left (ZigzagMapSourceCardinalityMismatch (HomologicalDegree 0) 1 2) -> pure () + Left failure -> assertFailure ("wrong shape refusal: " <> show failure) + Right _ -> assertFailure "malformed chain map was admitted" + +rejectChainMapLaw :: Assertion +rejectChainMapLaw = do + complexValue <- requireRight "interval complex" intervalComplex + degreeZeroIdentity <- requireRight "degree-zero identity" (identityIncidence 2) + case + mkFiniteChainMapChecked + complexValue + complexValue + ( \(HomologicalDegree degreeIndex) -> + case degreeIndex of + 0 -> degreeZeroIdentity + 1 -> emptyBoundaryIncidenceOf 1 1 + _ -> emptyBoundaryIncidence + ) + of + Left (ZigzagChainMapLawViolation (HomologicalDegree 1)) -> pure () + Left failure -> assertFailure ("wrong chain-law refusal: " <> show failure) + Right _ -> assertFailure "noncommuting chain map was admitted" + +rejectEndpointMismatch :: Assertion +rejectEndpointMismatch = do + identityMap <- requireRight "endpoint identity map" (coordinateMap (zeroComplex 1) (zeroComplex 1) [(0, 0)]) + case mkFiniteChainZigzag (zeroComplex 2) [ForwardArrow identityMap] of + Left (ZigzagEndpointMismatch 0 ZigzagForward) -> pure () + Left failure -> assertFailure ("wrong endpoint refusal: " <> show failure) + Right _ -> assertFailure "mismatched diagram endpoint was admitted" + +zeroComplex :: Int -> FiniteChainComplex Int +zeroComplex dimension = + mkFiniteChainComplex + (HomologicalDegree 0) + (const (emptyBoundaryIncidenceOf (fromIntegral dimension) 0)) + +poisonedPointComplex :: FiniteChainComplex Int +poisonedPointComplex = + mkFiniteChainComplex + (HomologicalDegree 0) + ( \(HomologicalDegree degreeIndex) -> + if degreeIndex == 0 + then emptyBoundaryIncidenceOf 1 0 + else emptyBoundaryIncidenceOf 7 6 + ) + +intervalComplex :: Either String (FiniteChainComplex Int) +intervalComplex = do + degreeOneBoundary <- + first show + ( mkBoundaryIncidence + 1 + 2 + [mkBoundaryEntry 0 0 (-1 :: Int), mkBoundaryEntry 0 1 1] + ) + pure + ( mkFiniteChainComplex + (HomologicalDegree 1) + ( \(HomologicalDegree degreeIndex) -> + case degreeIndex of + 0 -> emptyBoundaryIncidenceOf 2 0 + 1 -> degreeOneBoundary + _ -> emptyBoundaryIncidence + ) + ) + +circleComplex :: Either String (FiniteChainComplex Int) +circleComplex = do + degreeOneBoundary <- + first show + ( mkBoundaryIncidence + 3 + 3 + [ mkBoundaryEntry 0 0 (-1 :: Int) + , mkBoundaryEntry 0 1 1 + , mkBoundaryEntry 1 1 (-1) + , mkBoundaryEntry 1 2 1 + , mkBoundaryEntry 2 2 (-1) + , mkBoundaryEntry 2 0 1 + ] + ) + pure + ( mkFiniteChainComplex + (HomologicalDegree 1) + ( \(HomologicalDegree degreeIndex) -> + case degreeIndex of + 0 -> emptyBoundaryIncidenceOf 3 0 + 1 -> degreeOneBoundary + _ -> emptyBoundaryIncidence + ) + ) + +coordinateMap :: + FiniteChainComplex Int -> + FiniteChainComplex Int -> + [(Int, Int)] -> + Either ZigzagFailure (FiniteChainMap Int) +coordinateMap sourceComplex targetComplex coordinates = do + degreeZeroMap <- + first (ZigzagMapComponentInvalid (HomologicalDegree 0)) + ( mkBoundaryIncidence + (fromIntegral (complexZeroDimension sourceComplex)) + (fromIntegral (complexZeroDimension targetComplex)) + ( fmap + ( \(sourceIndexValue, targetIndexValue) -> + mkBoundaryEntry + (fromIntegral sourceIndexValue) + (fromIntegral targetIndexValue) + (1 :: Int) + ) + coordinates + ) + ) + mkFiniteChainMapChecked sourceComplex targetComplex $ \(HomologicalDegree degreeIndex) -> + case degreeIndex of + 0 -> degreeZeroMap + _ -> emptyBoundaryIncidence + +identityIncidence :: Int -> Either BoundaryIncidenceShapeError (BoundaryIncidence Int) +identityIncidence dimension = + mkBoundaryIncidence + (fromIntegral dimension) + (fromIntegral dimension) + ( fmap + (\indexValue -> mkBoundaryEntry (fromIntegral indexValue) (fromIntegral indexValue) (1 :: Int)) + [0 .. dimension - 1] + ) + +complexZeroDimension :: FiniteChainComplex r -> Int +complexZeroDimension complexValue = degreeCardinality complexValue (HomologicalDegree 0) + +requireRight :: Show failure => String -> Either failure value -> IO value +requireRight context = + either (assertFailure . ((context <> ": ") <>) . show) pure