Taking Cubical Homology Seriously (#72)

* add option to turn emergent pairs off, add columns_to_reduce

* improve usability of postprocess_interval

* fix emergent checking

* remove emergent switch

* code organization

* major refactoring: rename diam to birth, make index mandatory in AbstractSimplex

* new cubical homology implementation

* update docs

* go back to using cohomology terminology in the main algorithm

* enable emergnet pairs for homology

* fix plotting recipe

* remove dead code and improve tests

* bump version

Project.toml

@@ -1,7 +1,7 @@
 name = "Ripserer"
 uuid = "aa79e827-bd0b-42a8-9f10-2b302677a641"
 authors = ["mtsch <matijacufar@gmail.com>"]
-version = "0.13.2"
+version = "0.14.0"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"

benchmark/bench_image.jl

@@ -4,14 +4,14 @@ using FileIO
 using Ripserer
 hubble1024 = getfield.(load(joinpath(@__DIR__, "hubble1024px.jpg")), :r)
 
-# TODO: add bonsai when supported by new version of Cubical
-# bonsai256 = Array{UInt8}(undef, (256, 256, 256))
-# read!(joinpath(@__DIR__, "bonsai.raw"), bonsai256)
-# bonsai64 = bonsai256[1:4:end, 1:4:end, 1:4:end]
+bonsai256 = Array{UInt8}(undef, (256, 256, 256))
+read!(joinpath(@__DIR__, "bonsai.raw"), bonsai256)
+bonsai64 = bonsai256[1:4:end, 1:4:end, 1:4:end]
 
 suite = BenchmarkGroup()
 
-suite["hubble"] = @benchmarkable ripserer(Cubical{Int128}($hubble1024))
+suite["hubble"] = @benchmarkable ripserer(Cubical($hubble1024))
+suite["bonsai"] = @benchmarkable ripserer(Cubical($bonsai64))
 
 end
 BenchImage.suite

docs/src/api/filtrations.md

@@ -36,6 +36,14 @@ Ripserer.unsafe_simplex
 Ripserer.unsafe_cofacet
 ```
 
+```@docs
+Ripserer.columns_to_reduce
+```
+
+```@docs
+Ripserer.emergent_pairs
+```
+
 ```@docs
 Ripserer.postprocess_interval
 ```

docs/src/api/ripserer.md

@@ -29,7 +29,7 @@ Simplex
 ```
 
 ```@docs
-Cubelet
+Cube
 ```
 
 ```@docs

docs/src/api/simplices.md

@@ -9,7 +9,11 @@ dim(::Ripserer.AbstractSimplex)
 ```
 
 ```@docs
-diam(::Ripserer.AbstractSimplex)
+birth(::Ripserer.AbstractSimplex)
+```
+
+```@docs
+index(::Ripserer.AbstractSimplex)
 ```
 
 ```@docs
@@ -31,11 +35,3 @@ Ripserer.coboundary(::Any, ::Ripserer.AbstractSimplex)
 ```@docs
 Ripserer.boundary(::Any, ::Ripserer.AbstractSimplex)
 ```
-
-```@docs
-Ripserer.IndexedSimplex
-```
-
-```@docs
-Ripserer.index(::Ripserer.AbstractSimplex)
-```

docs/src/examples/sublevelset.jl

@@ -68,11 +68,14 @@ end
 
 # To get the locations of the minima, extract the critical simplices from intervals. Since
 # simplices act like collections of vertex indices, and zero-dimensional representatives
-# have a single point each, we can use `only` to extract them.
+# have a single point each, we can use `only` to extract them. Cubical simplices store
+# vertices as `CartesianIndex`, so we need to index into the x-values with them in order to
+# plot them.
 
-x_min = [only(birth_simplex(int)) for int in result]
+min_indices = [only(birth_simplex(int)) for int in result]
+min_x = eachindex(curve)[min_indices]
 
-scatter!(plt, x_min, curve[x_min], color=1:6, markershape=:star)
+scatter!(plt, min_x, curve[min_x], color=1:6, markershape=:star)
 plot(plt, plot(result, markercolor=1:6, markeralpha=1))
 
 
@@ -104,6 +107,7 @@ plot(filtered_result, persistence=true)
 
 # We could achieve the same result by simply filtering the resulting diagram.
 
+@assert filter(int -> persistence(int) > 1, noisy_result) == filtered_result # hide
 filter(int -> persistence(int) > 1, noisy_result) == filtered_result
 
 # Keep in mind, however, that for very large diagrams, collecting representatives on all of
@@ -122,9 +126,10 @@ for interval in filter(isfinite, filtered_result)
     plot!(plt, interval, noisy_curve, seriestype=:path)
 end
 
-x_min = [first(birth_simplex(int)) for int in filtered_result]
+min_indices = [first(birth_simplex(int)) for int in filtered_result]
+min_x = eachindex(curve)[min_indices]
 
-scatter!(plt, x_min, noisy_curve[x_min], color=eachindex(filtered_result), markershape=:star)
+scatter!(plt, min_x, noisy_curve[min_x], color=eachindex(filtered_result), markershape=:star)
 plot(plt, plot(filtered_result, markercolor=eachindex(filtered_result), markeralpha=1))
 
 # The result we got is very similar even though the data set was very noisy.
@@ -136,29 +141,22 @@ plot(plt, plot(filtered_result, markercolor=eachindex(filtered_result), markeral
 # Instead of looking for local minima, let's look for local maxima. To do that, we have to
 # invert the image.
 
-result = ripserer(Cubical(1 .- blackhole))
-
-# !!! tip "Overflow"
-#     Ripserer stores vertices in large integers. In some cases, `Int64` is
-#     not large enough to represent them. This is particularly likely to happen with
-#     `Cubical`. When it does, `ripserer` will throw an `OverflowError`. If that happens,
-#     try constructing the `Cubical` filtration with `Cubical{Int128}(data)`.
-
+result = ripserer(Cubical(-blackhole))
 plot(blackhole_plot, plot(result))
 
 # We notice there are quite a lot of intervals along the diagonal. These correspond to local
 # geometry of the image, so we are not interested in them right now. To filter them out, we
 # set a cutoff.
 
-result = ripserer(Cubical(1 .- blackhole), cutoff=0.1)
+result = ripserer(Cubical(-blackhole), cutoff=0.1)
 plot(blackhole_plot, plot(result))
 
 # Like earlier, we can show the local extrema in the image. We will show a different way to
 # plot them. Let's use the `threshold` argument with `plot`, which only keeps parts of the
 # representative with a diameter equal to or lower than `threshold`. If we needed a strict
 # `<`, we could use `threshold_strict`.
 
-result = ripserer(Cubical(1 .- blackhole), cutoff=0.1, reps=true)
+result = ripserer(Cubical(-blackhole), cutoff=0.1, reps=true)
 plt = plot(blackhole_image, title="Black Hole")
 for interval in result[1]
     plot!(plt, interval, blackhole, threshold=birth(interval))
@@ -188,7 +186,7 @@ plot!(plt, only(result[2]), blackhole, label="H₁ cocycle", color=1)
 #     Ripserer currently can't compute infinite intervals in dimensions
 #     higher than zero with persistent homology.
 
-result = ripserer(Cubical(1 .- blackhole), cutoff=0.1, reps=true, cohomology=false)
+result = ripserer(Cubical(-blackhole), cutoff=0.1, reps=true, cohomology=false)
 plot!(plt, only(result[2]), blackhole, label="H₁ cycle", color=3)
 
 # We have successfully found the hole in a black hole.

src/Ripserer.jl

@@ -10,6 +10,7 @@ module Ripserer
 using Compat
 
 using Base: @propagate_inbounds, @pure
+using Base.Cartesian
 using LinearAlgebra
 using SparseArrays
 
@@ -29,24 +30,24 @@ export
     birth, death, persistence, representative, birth_simplex, death_simplex, barcode
 
 export
-    Mod, Simplex, Cubelet, vertices, diam, dim, threshold, simplex, coefficient,
+    Mod, Simplex, Cube, index, vertices, dim, threshold, simplex, coefficient,
     Rips, SparseRips, Cubical, ripserer
 
 
 include("primefield.jl")
 
-include("abstractsimplex.jl")
 include("simplex.jl")
 
 include("abstractfiltration.jl")
 include("ripsfiltration.jl")
-include("cubical.jl")
 
 include("chainelement.jl")
 include("zerodimensional.jl")
 include("reductionmatrix.jl")
 include("main.jl")
 
+include("cubical.jl")
+
 include("simplexrecipes.jl")
 
 end

src/abstractfiltration.jl

@@ -15,6 +15,8 @@ A filtration is used to find the edges in filtration and to create simplices. An
 * [`unsafe_cofacet`](@ref)`(::AbstractFiltration, simplex, vertices, vertex[, sign, edges])`
 * [`birth(::AbstractFiltration, v)`](@ref)
 * [`threshold(::AbstractFiltration)`](@ref)
+* [`columns_to_reduce(::AbstractFiltration)`](@ref)
+* [`emergent_pairs(::AbstractFiltration)`](@ref)
 * [`postprocess_interval(::AbstractFiltration, ::Any)`](@ref)
 """
 abstract type AbstractFiltration{I, T} end
@@ -125,10 +127,35 @@ can have. Defaults to `Inf`.
 """
 threshold(::AbstractFiltration) = Inf
 
+"""
+    columns_to_reduce(::AbstractFilration, prev_column_itr)
+
+List all columns to reduce in next dimension, possibly computing it from previous
+columns. Default implementation uses [`coboundary`](@ref) with the all cofacets parameter
+set to `Val(false)`.
+"""
+function columns_to_reduce(flt::AbstractFiltration, prev_column_itr)
+    return Iterators.flatten(
+        imap(σ -> coboundary(flt, σ, Val(false)), prev_column_itr)
+    )
+end
+
+"""
+    emergent_pairs(::AbstractFiltration)
+
+Perform the emergent pairs optimization. Default to returning `true`. Should be set to
+`false` for a filtration type that is unable to produce (co)boundary simplices in the
+correct order.
+"""
+emergent_pairs(::AbstractFiltration) = true
+
 """
     postprocess_interval(::AbstractFiltration, interval)
 
-This function is called on each resulting persistence interval. The default implementation
-does nothing.
+This function is called on each resulting persistence interval. If returns `nothing`, the
+interval is skipped. The default implementation returns the unchanged interval.
 """
 postprocess_interval(::AbstractFiltration, interval) = interval
+
+_postprocess_interval(flt, interval) = postprocess_interval(flt, interval)
+_postprocess_interval(flt, ::Nothing) = nothing