Started to write code for traversing and computing homology dimension…

… by dimension, instead of streaming all dimensions in the same stream.
appliedtopology · May 13, 2024 · 5d7a5ae · 5d7a5ae
1 parent 8860875
commit 5d7a5ae
Show file tree

Hide file tree

Showing 7 changed files with 271 additions and 86 deletions.
diff --git a/src/main/scala/org/appliedtopology/tda4j/Cube.scala b/src/main/scala/org/appliedtopology/tda4j/Cube.scala
@@ -22,56 +22,58 @@ given Ordering[ElementaryInterval] = Ordering.by(i => i.n)
 given Ordering[ElementaryCube] = Ordering.by(c => c.intervals)
 
 case class ElementaryCube(val intervals: List[ElementaryInterval]) extends Cell[ElementaryCube] {
-  override def boundary[CoefficientT: Fractional]: Chain[ElementaryCube, CoefficientT] = {
-    val chainOps = ChainOps[ElementaryCube, CoefficientT]()
-    import chainOps.{*, given}
-    val fr = summon[Fractional[CoefficientT]]
-    import math.Fractional.Implicits.infixFractionalOps
-
-    // For a cube that decomposes as Q = I*P where I is a single interval, we define
-    // ∂Q = ∂I * P + (-1)^{dim I} I * ∂P
-
-    // Given stuff-already-processed and an I*P decomposition, figure out whether this I changes the accumulated
-    // sign, and produce the ∂I * P + sign pair
-    def process(
-      left: List[ElementaryInterval],
-      current: ElementaryInterval,
-      right: List[ElementaryInterval]
-    ): (CoefficientT, Chain[ElementaryCube, CoefficientT]) = current match {
-      case DegenerateInterval(n) => (fr.one, Chain())
-      case FullInterval(n) =>
-        (
-          fr.negate(fr.one),
-          Chain(
-            ElementaryCube(left ++ (DegenerateInterval(n + 1) :: right)) -> fr.one,
-            ElementaryCube(left ++ (DegenerateInterval(n) :: right)) -> fr.negate(fr.one)
+  override def boundary[CoefficientT: Fractional]: Chain[ElementaryCube, CoefficientT] =
+    if(dim == 0) Chain()
+    else {
+      val chainOps = ChainOps[ElementaryCube, CoefficientT]()
+      import chainOps.{*, given}
+      val fr = summon[Fractional[CoefficientT]]
+      import math.Fractional.Implicits.infixFractionalOps
+
+      // For a cube that decomposes as Q = I*P where I is a single interval, we define
+      // ∂Q = ∂I * P + (-1)^{dim I} I * ∂P
+
+      // Given stuff-already-processed and an I*P decomposition, figure out whether this I changes the accumulated
+      // sign, and produce the ∂I * P + sign pair
+      def process(
+        left: List[ElementaryInterval],
+        current: ElementaryInterval,
+        right: List[ElementaryInterval]
+      ): (CoefficientT, Chain[ElementaryCube, CoefficientT]) = current match {
+        case DegenerateInterval(n) => (fr.one, Chain())
+        case FullInterval(n) =>
+          (
+            fr.negate(fr.one),
+            Chain(
+              ElementaryCube(left ++ (DegenerateInterval(n + 1) :: right)) -> fr.one,
+              ElementaryCube(left ++ (DegenerateInterval(n) :: right)) -> fr.negate(fr.one)
+            )
           )
-        )
-    }
+      }
 
-    @tailrec
-    def boundaryOf(
-      left: List[ElementaryInterval],
-      current: ElementaryInterval,
-      right: List[ElementaryInterval],
-      sign: CoefficientT,
-      acc: Chain[ElementaryCube, CoefficientT]
-    ): Chain[ElementaryCube, CoefficientT] = {
-      val (sgn, newchain) = process(left, current, right)
-      right match {
-        case Nil =>
-          // process current, then return the resulting acc
-          acc + ((sign * sgn) ⊠ newchain)
-        case c :: cs =>
-          // process current, then call with c as new current
-          boundaryOf(left.appended(current), c, cs, sign * sgn, acc + ((sign * sgn) ⊠ newchain))
+      @tailrec
+      def boundaryOf(
+        left: List[ElementaryInterval],
+        current: ElementaryInterval,
+        right: List[ElementaryInterval],
+        sign: CoefficientT,
+        acc: Chain[ElementaryCube, CoefficientT]
+      ): Chain[ElementaryCube, CoefficientT] = {
+        val (sgn, newchain) = process(left, current, right)
+        right match {
+          case Nil =>
+            // process current, then return the resulting acc
+            acc + ((sign * sgn) ⊠ newchain)
+          case c :: cs =>
+            // process current, then call with c as new current
+            boundaryOf(left.appended(current), c, cs, sign * sgn, acc + ((sign * sgn) ⊠ newchain))
+        }
+      }
+      intervals match {
+        case Nil       => Chain()
+        case (i :: is) => boundaryOf(List.empty, i, is, fr.one, Chain())
       }
     }
-    intervals match {
-      case Nil       => Chain()
-      case (i :: is) => boundaryOf(List.empty, i, is, fr.one, Chain())
-    }
-  }
 
   def emb: Int = intervals.size
   def dim: Int = intervals.count(i => i.isInstanceOf[FullInterval])
@@ -83,31 +85,33 @@ case class ElementaryCube(val intervals: List[ElementaryInterval]) extends Cell[
     s"Cubical[${intervals.map(_.toString).mkString("x")}]"
 }
 
-trait CubeStream[FiltrationT : Ordering]
-  extends Filtration[ElementaryCube,FiltrationT]
-  with IterableOnce[ElementaryCube]
-
-
+trait CubeStream[FiltrationT: Ordering] extends CellStream[ElementaryCube, FiltrationT]
 
 object Cubical {
-  def from[T:Filterable:Ordering](grid: Array[Array[T]]): CubeStream[T] = new CubeStream[T] {
+  def from[T: Filterable: Ordering](grid: Array[Array[T]]): CubeStream[T] = new CubeStream[T] {
     val largest = summon[Filterable[T]].largest
-    val smallest =  summon[Filterable[T]].smallest
-
-    val cubes : Map[ElementaryCube,T] = Map.from((
-      for
-        (row,i) <- grid.zipWithIndex
-        (pixel,j) <- row.zipWithIndex
-        firstInterval <- Seq(FullInterval(i), DegenerateInterval(i))
-        secondInterval <- Seq(FullInterval(j), DegenerateInterval(j))
-      yield
-        ElementaryCube(List(firstInterval,secondInterval)) -> pixel
-    ).toList.sorted.reverse) // reversed so that the last element saved for any one cube is the one with lowest T-value
+    val smallest = summon[Filterable[T]].smallest
+
+    override def filtrationOrdering: Ordering[ElementaryCube] =
+      Ordering.by(cubes).orElseBy(_.dim).orElseBy(_.intervals.map{(int) => int.n})
+
+    val cubes: Map[ElementaryCube, T] = Map.from(
+      (
+        for
+          (row, i) <- grid.zipWithIndex
+          (pixel, j) <- row.zipWithIndex
+          firstInterval <- Seq(FullInterval(i), DegenerateInterval(i))
+          secondInterval <- Seq(FullInterval(j), DegenerateInterval(j))
+        yield ElementaryCube(List(firstInterval, secondInterval)) -> pixel
+      ).toList.sorted.reverse
+    ) // reversed so that the last element saved for any one cube is the one with lowest T-value
 
     override def filtrationValue: PartialFunction[ElementaryCube, T] =
       cubes.apply
 
     override def iterator: Iterator[ElementaryCube] =
-      cubes.keys.toList.sorted(using Ordering.by(cubes.apply).orElseBy(_.dim)(using ord=Ordering.Int.reverse)).iterator
+      cubes.keys.toList
+        .sorted(using Ordering.by(cubes.apply).orElseBy(_.dim)(using ord = Ordering.Int.reverse))
+        .iterator
   }
 }
diff --git a/src/main/scala/org/appliedtopology/tda4j/Homology.scala b/src/main/scala/org/appliedtopology/tda4j/Homology.scala
@@ -2,7 +2,7 @@ package org.appliedtopology.tda4j
 
 import org.appliedtopology.tda4j.barcode.PersistenceBar
 
-import collection.mutable
+import collection.{immutable, mutable}
 import scala.annotation.tailrec
 
 //class ReducedSimplicialHomologyContext[VertexT: Ordering, CoefficientT: Fractional, FiltrationT: Ordering]()
@@ -159,5 +159,158 @@ class CellularHomologyContext[CellT <: Cell[CellT]: Ordering, CoefficientT: Frac
       stream.smallest: FiltrationT, // computation done up until filtrationValue
       mutable.ArrayDeque.empty
     ) // torsion part of barcode
+}
+
+class SimplicialHomologyByDimensionContext[VertexT: Ordering, CoefficientT: Fractional]
+    extends ChainOps[Simplex[VertexT], CoefficientT]() {
+  case class HomologyState(
+    cycles: mutable.Map[Simplex[VertexT], Chain[Simplex[VertexT], CoefficientT]],
+    cyclesBornBy: mutable.Map[Simplex[VertexT], Simplex[VertexT]],
+    boundaries: mutable.Map[Simplex[VertexT], Chain[Simplex[VertexT], CoefficientT]],
+    boundariesBornBy: mutable.Map[Simplex[VertexT], Simplex[VertexT]],
+    coboundaries: mutable.Map[Simplex[VertexT], Chain[Simplex[VertexT], CoefficientT]],
+    stream: StratifiedCellStream[Simplex[VertexT], Double],
+    var current: Double,
+    var currentDim: Int,
+    var currentIterator: collection.BufferedIterator[Simplex[VertexT]],
+    barcode: mutable.Map[Int, immutable.Queue[(Double, Double, Chain[Simplex[VertexT], CoefficientT])]]
+  ) {
+    // first off, all vertices are immediately cycles
+    cycles.addAll(stream.iterateDimension(0).map(cell => cell -> Chain(cell)))
+    cyclesBornBy.addAll(cycles.map((cell, chain) => cell -> cell))
+
+    // secondly, we can read off homology completely from a minimal spanning tree
+    val kruskal = new Kruskal[Simplex[VertexT]](
+      cycles.keys.toSeq,
+      { (x: Simplex[VertexT], y: Simplex[VertexT]) => stream.filtrationValue(x.union(y)) }
+    )(using stream.filtrationOrdering)
+
+    kruskal.mstIterator.foreach { (src, tgt) =>
+      val edge: Simplex[VertexT] = src.union(tgt)
+
+      // the edge src -- tgt will connect src to tgt thus removing one of the cycles
+      val dEdge = edge.boundary
+      val dyingVertex = dEdge.leadingCell.get
+      boundaries.addOne(dEdge.leadingCell.get -> dEdge)
+      coboundaries.addOne(dyingVertex, dEdge)
+      barcode(0) =
+        barcode(0).appended((stream.filtrationValue(dyingVertex), stream.filtrationValue(edge), cycles(dyingVertex)))
+      cycles.remove(dyingVertex)
+    }
+
+    kruskal.cyclesIterator.foreach { (src, tgt) =>
+      val edge: Simplex[VertexT] = src.union(tgt)
+
+      // the edge src -- tgt will connect src to tgt thus closing a loop
+      val dEdge = edge.boundary
+      // TODO is it worth it to have a more complex UnionFind that allows us to get the entire path along the MST?
+      val (reduced, reductionLog): (Chain[Simplex[VertexT], CoefficientT], Chain[Simplex[VertexT], CoefficientT]) =
+        reduceBy(dEdge, boundaries)
+      val fr = summon[Fractional[CoefficientT]]
+      val coboundary: Chain[Simplex[VertexT], CoefficientT] =
+        reductionLog.items.foldRight(fr.negate(fr.one) ⊠ Chain(edge)) { (item, acc) =>
+          val (spx, coeff) = item
+          if (coboundaries.contains(spx))
+            acc + coeff ⊠ coboundaries(spx)
+          else
+            acc
+        }
+      cycles(coboundary.leadingCell.get) = coboundary
+      cyclesBornBy(coboundary.leadingCell.get) = edge
+    }
+
+    // setup is done, we should be ready to start dimension 2
+    currentDim = 1
+    current = Double.PositiveInfinity
+
+    @tailrec
+    private def reduceBy(
+                          z: Chain[Simplex[VertexT], CoefficientT],
+                          basis: mutable.Map[Simplex[VertexT], Chain[Simplex[VertexT], CoefficientT]],
+                          reductionLog: Chain[Simplex[VertexT], CoefficientT] = Chain()
+                        )(using
+                          fr: Fractional[CoefficientT]
+                        ): (Chain[Simplex[VertexT], CoefficientT], Chain[Simplex[VertexT], CoefficientT]) =
+      z.leadingCell match {
+        case None => (z, reductionLog)
+        case Some(sigma) =>
+          if (basis.contains(sigma)) {
+            val redCoeff = fr.div(z.leadingCoefficient, basis(sigma).leadingCoefficient)
+            reduceBy(z - redCoeff ⊠ basis(sigma), basis, reductionLog + redCoeff ⊠ Chain(sigma))
+          } else (z, reductionLog)
+      }
+
+    def advanceOne(): Unit =
+      if (currentIterator.hasNext) {
+        val fr = summon[Fractional[CoefficientT]]
+        val sigma = currentIterator.next()
+        val dsigma : Chain[Simplex[VertexT], CoefficientT] = sigma.boundary
+        val (dsigmaReduced, reduction) = reduceBy(dsigma, boundaries)
+        val coboundary = reduction.items.foldRight(fr.negate(fr.one) ⊠ Chain(sigma)) { (next, acc) =>
+          val (spx, coeff) = next
+          if (coboundaries.contains(spx))
+            acc + coeff ⊠ coboundaries(spx)
+          else
+            acc
+        }
+        if (dsigmaReduced.isZero()) {
+          // adding a boundary to a boundary creates a new cycle as sigma + whatever whose boundary eliminated dsigma
+          cycles(coboundary.leadingCell.get) = coboundary
+          cyclesBornBy(coboundary.leadingCell.get) = sigma
+        } else {
+          // we have a new boundary witnessed
+          boundaries(dsigmaReduced.leadingCell.get) = dsigmaReduced
+          boundariesBornBy(dsigmaReduced.leadingCell.get) = sigma
+          coboundaries(dsigmaReduced.leadingCell.get) = coboundary
 
+          val (_, cycleBasis) = reduceBy(dsigmaReduced, cycles)
+          val representativeCycle: Chain[Simplex[VertexT], CoefficientT] = cycleBasis.leadingCell match {
+            case None => Chain()
+            case Some(cell) => cycles(cell)
+          }
+          cycleBasis.leadingCell match {
+            case None => ()
+            case Some(cell) => cycles.remove(cell)
+          }
+
+          val lower: Double = cycleBasis.leadingCell match {
+            case None => Double.NegativeInfinity
+            case Some(spx) =>
+              stream.filtrationValue.orElse{(_) => Double.NegativeInfinity}.compose(cyclesBornBy)(spx)
+          }
+          val upper: Double =
+            stream.filtrationValue.orElse{(_) => Double.PositiveInfinity}(sigma)
+
+          barcode(currentDim) = barcode(currentDim).appended((lower, upper, representativeCycle))
+        }
+        current = stream.filtrationValue.lift(sigma).getOrElse(stream.smallest)
+      } else {
+        currentDim += 1
+        currentIterator = stream
+          .iterateDimension
+          .applyOrElse(currentDim, _ => Iterator.empty)
+          .buffered
+        current = Double.NegativeInfinity
+      }
+
+    def advanceTo(dim: Int, f: Double = Double.PositiveInfinity): Unit =
+      while(currentIterator.hasNext &&
+        currentDim <= dim &&
+        f > current)
+        advanceOne()
+  }
+
+  def persistentHomology(stream: => StratifiedCellStream[Simplex[VertexT], Double]): HomologyState =
+    HomologyState(
+      mutable.Map.empty,
+      mutable.Map.empty,
+      mutable.Map.empty,
+      mutable.Map.empty,
+      mutable.Map.empty,
+      stream,
+      stream.smallest,
+      0,
+      Iterator.empty.buffered,
+      mutable.Map.empty
+    )
 }
diff --git a/src/main/scala/org/appliedtopology/tda4j/SimplexStream.scala b/src/main/scala/org/appliedtopology/tda4j/SimplexStream.scala
@@ -187,3 +187,17 @@ class FilteredSimplexOrdering[VertexT, FiltrationT](
         Ordering.Int.compare(x.size, y.size)
   }
 }
+
+
+
+trait StratifiedCellStream[CellT <: Cell[CellT] :Ordering,FiltrationT:Filterable] extends CellStream[CellT,FiltrationT] {
+  def iterateDimension : PartialFunction[Int, Iterator[CellT]]
+
+  override def iterator: Iterator[CellT] =
+    Iterator.from(0)
+      .filter(iterateDimension.isDefinedAt)
+      .map{ (dim:Int) => iterateDimension.applyOrElse(dim, (d:Int) => Iterator.empty) }
+      .fold(Iterator.empty:Iterator[CellT])((x,y) => x++y)
+}
+
+trait StratifiedSimplexStream[VertexT:Ordering, FiltrationT:Filterable] extends StratifiedCellStream[Simplex[VertexT],FiltrationT] { }
diff --git a/src/main/scala/org/appliedtopology/tda4j/UnionFind.scala b/src/main/scala/org/appliedtopology/tda4j/UnionFind.scala
@@ -28,17 +28,18 @@ class UnionFind[T](vertices: IterableOnce[T]) {
 /** This implementation of Kruskal's algorithm will return two iterators of vertex pairs: the first iterator is a
   * Minimal Spanning Tree in increasing weight order, while the second iterator gives all the non-included
   */
-class Kruskal[T](metricSpace: FiniteMetricSpace[T])(using
+
+class Kruskal[T](elements: Seq[T], distance: (T,T) => Double)(using
   orderingT: Ordering[T]
 ) {
-  val unionFind: UnionFind[T] = UnionFind(metricSpace.elements)
+  val unionFind: UnionFind[T] = UnionFind(elements)
 
   val sortedEdges: List[(Double, unionFind.UFSet, unionFind.UFSet)] =
     (for
       x <- unionFind.sets.keysIterator
       y <- unionFind.sets.keysIterator
       if orderingT.lt(x.label, y.label)
-    yield (metricSpace.distance(x.label, y.label), x, y)).toList.sortWith { (l, r) =>
+    yield (distance(x.label, y.label), x, y)).toList.sortWith { (l, r) =>
       l._1 < r._1
     }
 
@@ -55,6 +56,11 @@ class Kruskal[T](metricSpace: FiniteMetricSpace[T])(using
   val cyclesIterator: Iterator[(T, T)] = lrList._2.iterator
 }
 
+object Kruskal {
+  def apply[T:Ordering](metricSpace: FiniteMetricSpace[T]): Kruskal[T] =
+    new Kruskal(metricSpace.elements.toSeq, metricSpace.distance)
+}
+
 /*
 def KruskalF[T](metricSpace: FiniteMetricSpace[T])(using orderingT: Ordering[T]) = {
   val unionFind: UnionFind[T] = UnionFind(metricSpace.elements)

diff --git a/src/test/scala/org/appliedtopology/tda4j/APISpec.scala b/src/test/scala/org/appliedtopology/tda4j/APISpec.scala
@@ -2,6 +2,8 @@ package org.appliedtopology.tda4j
 
 import org.specs2.mutable
 
+import scala.math.Numeric.DoubleIsFractional
+
 class APISpec extends mutable.Specification {
   """Test case class for developing the non-Scala facing API functionality
     |and the non-expert API functionality""".stripMargin