moonlight-homology 0.1.0.2 → 0.1.0.3
raw patch · 17 files changed
+1488/−177 lines, 17 filesPVP: major bump suggested
API removals or changes: PVP suggests a major version bump
API changes (from Hackage documentation)
+ Moonlight.Homology: BackwardArrow :: map -> ZigzagArrow map
+ Moonlight.Homology: ChainMapSource :: ChainMapEndpoint
+ Moonlight.Homology: ChainMapTarget :: ChainMapEndpoint
+ Moonlight.Homology: ForwardArrow :: map -> ZigzagArrow map
+ Moonlight.Homology: ZigzagBackward :: ZigzagDirection
+ Moonlight.Homology: ZigzagChainMapLawViolation :: HomologicalDegree -> ZigzagFailure
+ Moonlight.Homology: ZigzagComplexInvalid :: ChainMapEndpoint -> HomologyFailure -> ZigzagFailure
+ Moonlight.Homology: ZigzagEndpointMismatch :: Int -> ZigzagDirection -> ZigzagFailure
+ Moonlight.Homology: ZigzagForward :: ZigzagDirection
+ Moonlight.Homology: ZigzagHomologyCoordinatesMissing :: Int -> HomologicalDegree -> ZigzagFailure
+ Moonlight.Homology: ZigzagInterval :: HomologicalDegree -> endpoint -> endpoint -> Int -> ZigzagInterval endpoint
+ Moonlight.Homology: ZigzagMapComponentInvalid :: HomologicalDegree -> BoundaryIncidenceShapeError -> ZigzagFailure
+ Moonlight.Homology: ZigzagMapCompositionInvalid :: HomologicalDegree -> BoundaryIncidenceShapeError -> ZigzagFailure
+ Moonlight.Homology: ZigzagMapSourceCardinalityMismatch :: HomologicalDegree -> Int -> Int -> ZigzagFailure
+ Moonlight.Homology: ZigzagMapTargetCardinalityMismatch :: HomologicalDegree -> Int -> Int -> ZigzagFailure
+ Moonlight.Homology: ZigzagNegativeIntervalMultiplicity :: HomologicalDegree -> Int -> Int -> Int -> ZigzagFailure
+ Moonlight.Homology: [zigzagIntervalDegree] :: ZigzagInterval endpoint -> HomologicalDegree
+ Moonlight.Homology: [zigzagIntervalFirst] :: ZigzagInterval endpoint -> endpoint
+ Moonlight.Homology: [zigzagIntervalLast] :: ZigzagInterval endpoint -> endpoint
+ Moonlight.Homology: [zigzagIntervalMultiplicity] :: ZigzagInterval endpoint -> Int
+ Moonlight.Homology: data ChainMapEndpoint
+ Moonlight.Homology: data FiniteChainMap r
+ Moonlight.Homology: data FiniteChainZigzag r
+ Moonlight.Homology: data ZigzagArrow map
+ Moonlight.Homology: data ZigzagDirection
+ Moonlight.Homology: data ZigzagFailure
+ Moonlight.Homology: data ZigzagInterval endpoint
+ Moonlight.Homology: finiteChainMapAt :: FiniteChainMap r -> HomologicalDegree -> BoundaryIncidence r
+ Moonlight.Homology: finiteChainMapSource :: FiniteChainMap r -> FiniteChainComplex r
+ Moonlight.Homology: finiteChainMapTarget :: FiniteChainMap r -> FiniteChainComplex r
+ Moonlight.Homology: finiteChainZigzagArrows :: FiniteChainZigzag r -> [ZigzagArrow (FiniteChainMap r)]
+ Moonlight.Homology: finiteChainZigzagComplexes :: FiniteChainZigzag r -> NonEmpty (FiniteChainComplex r)
+ Moonlight.Homology: mkFiniteChainMapChecked :: (Eq r, Num r, Semiring r) => FiniteChainComplex r -> FiniteChainComplex r -> (HomologicalDegree -> BoundaryIncidence r) -> Either ZigzagFailure (FiniteChainMap r)
+ Moonlight.Homology: mkFiniteChainZigzag :: Eq r => FiniteChainComplex r -> [ZigzagArrow (FiniteChainMap r)] -> Either ZigzagFailure (FiniteChainZigzag r)
+ Moonlight.Homology: rationalZigzagIntervals :: Integral r => FiniteChainZigzag r -> Either ZigzagFailure [ZigzagInterval Int]
+ Moonlight.Homology: zigzagArrowDirection :: ZigzagArrow r -> ZigzagDirection
+ Moonlight.Homology: zigzagBettiAt :: Int -> [ZigzagInterval Int] -> Map HomologicalDegree Int
+ Moonlight.Homology.Persistence: BackwardArrow :: map -> ZigzagArrow map
+ Moonlight.Homology.Persistence: ChainMapSource :: ChainMapEndpoint
+ Moonlight.Homology.Persistence: ChainMapTarget :: ChainMapEndpoint
+ Moonlight.Homology.Persistence: ForwardArrow :: map -> ZigzagArrow map
+ Moonlight.Homology.Persistence: ZigzagBackward :: ZigzagDirection
+ Moonlight.Homology.Persistence: ZigzagChainMapLawViolation :: HomologicalDegree -> ZigzagFailure
+ Moonlight.Homology.Persistence: ZigzagComplexInvalid :: ChainMapEndpoint -> HomologyFailure -> ZigzagFailure
+ Moonlight.Homology.Persistence: ZigzagEndpointMismatch :: Int -> ZigzagDirection -> ZigzagFailure
+ Moonlight.Homology.Persistence: ZigzagForward :: ZigzagDirection
+ Moonlight.Homology.Persistence: ZigzagHomologyCoordinatesMissing :: Int -> HomologicalDegree -> ZigzagFailure
+ Moonlight.Homology.Persistence: ZigzagInterval :: HomologicalDegree -> endpoint -> endpoint -> Int -> ZigzagInterval endpoint
+ Moonlight.Homology.Persistence: ZigzagMapComponentInvalid :: HomologicalDegree -> BoundaryIncidenceShapeError -> ZigzagFailure
+ Moonlight.Homology.Persistence: ZigzagMapCompositionInvalid :: HomologicalDegree -> BoundaryIncidenceShapeError -> ZigzagFailure
+ Moonlight.Homology.Persistence: ZigzagMapSourceCardinalityMismatch :: HomologicalDegree -> Int -> Int -> ZigzagFailure
+ Moonlight.Homology.Persistence: ZigzagMapTargetCardinalityMismatch :: HomologicalDegree -> Int -> Int -> ZigzagFailure
+ Moonlight.Homology.Persistence: ZigzagNegativeIntervalMultiplicity :: HomologicalDegree -> Int -> Int -> Int -> ZigzagFailure
+ Moonlight.Homology.Persistence: [zigzagIntervalDegree] :: ZigzagInterval endpoint -> HomologicalDegree
+ Moonlight.Homology.Persistence: [zigzagIntervalFirst] :: ZigzagInterval endpoint -> endpoint
+ Moonlight.Homology.Persistence: [zigzagIntervalLast] :: ZigzagInterval endpoint -> endpoint
+ Moonlight.Homology.Persistence: [zigzagIntervalMultiplicity] :: ZigzagInterval endpoint -> Int
+ Moonlight.Homology.Persistence: data ChainMapEndpoint
+ Moonlight.Homology.Persistence: data FiniteChainMap r
+ Moonlight.Homology.Persistence: data FiniteChainZigzag r
+ Moonlight.Homology.Persistence: data ZigzagArrow map
+ Moonlight.Homology.Persistence: data ZigzagDirection
+ Moonlight.Homology.Persistence: data ZigzagFailure
+ Moonlight.Homology.Persistence: data ZigzagInterval endpoint
+ Moonlight.Homology.Persistence: finiteChainMapAt :: FiniteChainMap r -> HomologicalDegree -> BoundaryIncidence r
+ Moonlight.Homology.Persistence: finiteChainMapSource :: FiniteChainMap r -> FiniteChainComplex r
+ Moonlight.Homology.Persistence: finiteChainMapTarget :: FiniteChainMap r -> FiniteChainComplex r
+ Moonlight.Homology.Persistence: finiteChainZigzagArrows :: FiniteChainZigzag r -> [ZigzagArrow (FiniteChainMap r)]
+ Moonlight.Homology.Persistence: finiteChainZigzagComplexes :: FiniteChainZigzag r -> NonEmpty (FiniteChainComplex r)
+ Moonlight.Homology.Persistence: mkFiniteChainMapChecked :: (Eq r, Num r, Semiring r) => FiniteChainComplex r -> FiniteChainComplex r -> (HomologicalDegree -> BoundaryIncidence r) -> Either ZigzagFailure (FiniteChainMap r)
+ Moonlight.Homology.Persistence: mkFiniteChainZigzag :: Eq r => FiniteChainComplex r -> [ZigzagArrow (FiniteChainMap r)] -> Either ZigzagFailure (FiniteChainZigzag r)
+ Moonlight.Homology.Persistence: rationalZigzagIntervals :: Integral r => FiniteChainZigzag r -> Either ZigzagFailure [ZigzagInterval Int]
+ Moonlight.Homology.Persistence: zigzagArrowDirection :: ZigzagArrow r -> ZigzagDirection
+ Moonlight.Homology.Persistence: zigzagBettiAt :: Int -> [ZigzagInterval Int] -> Map HomologicalDegree Int
- Moonlight.Homology: sparseQuotientRepresentatives :: HomologicalDegree -> Int -> SparseMatrix -> SparseMatrix -> [RepresentativeChain Rational Int]
+ Moonlight.Homology: sparseQuotientRepresentatives :: HomologicalDegree -> Int -> SparseMatrix -> [SparseRow] -> [RepresentativeChain Rational Int]
- Moonlight.Homology.Topology: sparseQuotientRepresentatives :: HomologicalDegree -> Int -> SparseMatrix -> SparseMatrix -> [RepresentativeChain Rational Int]
+ Moonlight.Homology.Topology: sparseQuotientRepresentatives :: HomologicalDegree -> Int -> SparseMatrix -> [SparseRow] -> [RepresentativeChain Rational Int]
Files
- CHANGELOG.md +9/−0
- README.md +6/−2
- bench/Main.hs +5/−1
- bench/topology/ZigzagPersistence.hs +139/−0
- moonlight-homology.cabal +6/−3
- src-matrix/Moonlight/Homology/Pure/Matrix/SparseLinAlg.hs +157/−115
- src-public/Moonlight/Homology.hs +17/−0
- src-public/Moonlight/Homology/Persistence.hs +36/−0
- src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Build.hs +2/−2
- src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Linear.hs +8/−27
- src-topology/Moonlight/Homology/Pure/Topology/Algebra.hs +9/−14
- src-topology/Moonlight/Homology/Pure/Topology/SparseAlgebra.hs +9/−6
- src-topology/Moonlight/Homology/Pure/Topology/Zigzag.hs +692/−0
- test/Main.hs +3/−1
- test/facade/CompileFailSpec.hs +1/−1
- test/sequence/SpectralSpec.hs +1/−5
- test/topology/ZigzagSpec.hs +388/−0
CHANGELOG.md view
@@ -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
README.md view
@@ -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. |
bench/Main.hs view
@@ -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 ]
+ bench/topology/ZigzagPersistence.hs view
@@ -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
moonlight-homology.cabal view
@@ -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
src-matrix/Moonlight/Homology/Pure/Matrix/SparseLinAlg.hs view
@@ -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 =
src-public/Moonlight/Homology.hs view
@@ -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 (..),
src-public/Moonlight/Homology/Persistence.hs view
@@ -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
src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Build.hs view
@@ -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
src-sequence/Moonlight/Homology/Pure/Sequence/Spectral/Linear.hs view
@@ -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 _
src-topology/Moonlight/Homology/Pure/Topology/Algebra.hs view
@@ -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 ::
src-topology/Moonlight/Homology/Pure/Topology/SparseAlgebra.hs view
@@ -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
+ src-topology/Moonlight/Homology/Pure/Topology/Zigzag.hs view
@@ -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
test/Main.hs view
@@ -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
test/facade/CompileFailSpec.hs view
@@ -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 =>
test/sequence/SpectralSpec.hs view
@@ -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 =
+ test/topology/ZigzagSpec.hs view
@@ -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