packages feed

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 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