Homology (#49)

* code reorganization
* homology implementation
* use views for ReducedMatrix columns
* fix type instabilities
* sizehint to_skip in matrix assembly
* update docs, fix bugs and rename representatives
* bump version

Project.toml

@@ -1,7 +1,7 @@
 name = "Ripserer"
 uuid = "aa79e827-bd0b-42a8-9f10-2b302677a641"
 authors = ["mtsch <matijacufar@gmail.com>"]
-version = "0.6.4"
+version = "0.7.0"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
@@ -33,7 +33,8 @@ julia = "1"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Random", "SafeTestsets", "Test"]
+test = ["Aqua", "Random", "SafeTestsets", "Suppressor", "Test"]

docs/src/examples/basics.jl

@@ -134,29 +134,24 @@ most_persistent = sort(result_cut[2], by=persistence)[end]
 # that time. This means that we could stop computing when we reach this time and the result
 # should not change. Let's try it out!
 
-@time result_cut_thr = ripserer(cutout_2000, threshold=1.83)
+@time result_cut_thresh_2 = ripserer(cutout_2000, threshold=2)
 nothing; # hide
 
 #
 
-plot(result_cut_thr, title="Persistence Diagram, threshold=1.83")
+plot(result_cut_thresh_2, title="Persistence Diagram, threshold=2")
 
 # Indeed, the result is exactly the same, but it took less than a third of the time to
 # compute.
 
+result_cut_thresh_2 == result_cut
+
 # If we pick a threshold that is too low, we still detect the interval, but its death time
-# will become infinite. Let's demonstrate that with an animation.
-
-anim = @animate for threshold in 1.83:-0.03:0.09
-    result_cut_i = ripserer(cutout_2000, threshold=threshold)
-    pd_plt = plot(result_cut_i,
-                  xticks=0:0.1:1.9,
-                  yticks=0:0.1:1.9,
-                  title="threshold=$threshold")
-    bc_plt = barcode(result_cut_i[2],
-                     xticks=0:0.1:1.9,
-                     ylim=(1, length(result_cut[2])))
-    plot(pd_plt, bc_plt,
-         size=(800, 400))
-end
-mp4(anim, "basics_anim.mp4") # hide
+# becomes infinite.
+
+@time result_cut_thresh_1 = ripserer(cutout_2000, threshold=1)
+nothing; # hide
+
+#
+
+plot(result_cut_thresh_1, title="Persistence Diagram, threshold=1")

docs/src/examples/cocycles.jl

@@ -25,7 +25,7 @@ plot(curve(100), legend=false, title="Curve", aspect_ratio=1, xlab="x", ylab="y"
 
 anim = @animate for i in vcat(3:100, 100:-1:3)
     points = curve(i)
-    res = ripserer(points, representatives=true)
+    res = ripserer(points, reps=true)
 
     plt1 = plot(title="1-dimensional Representative Cocycle",
                 xlab="x",

docs/src/examples/image_sublevel.jl

@@ -30,8 +30,8 @@ heatmap(blackhole_grayscale, aspect_ratio=1, size=(800, 800))
 # !!! note "Note"
 #     Unlike with Rips filtrations, we can use large data sets here because the number of
 #     simplices in a cubical filtration is linear to the size of the input. We still want to
-#     be careful with calculating representatives because calculating representatives for
-#     thousands of persistence intervals can take a while.
+#     be careful with computing representatives because collecting them for thousands of
+#     persistence intervals can take a while.
 
 result_min = ripserer(Cubical(blackhole_grayscale))
 
@@ -52,7 +52,7 @@ plot(result_min, persistence=true)
 
 # Now that we know we have a smaller number of intervals, we can compute the
 # representatives.
-result_min = ripserer(Cubical(blackhole_grayscale), cutoff=0.1, representatives=true)
+result_min = ripserer(Cubical(blackhole_grayscale), cutoff=0.1, reps=true)
 
 # We separate the finite and infinite intervals. The finite one corresponds to the local
 # minimum in the middle of the image and the infinite one corresponds to the global minimum.
@@ -113,7 +113,7 @@ scatter!(result_min[2][2], blackhole_grayscale,
 
 # Let's repeat what we just did, but with the image negated.
 
-result_max = ripserer(Cubical(-blackhole_grayscale), cutoff=0.1, representatives=true)
+result_max = ripserer(Cubical(-blackhole_grayscale), cutoff=0.1, reps=true)
 plot(result_max)
 
 heatmap(blackhole_grayscale, aspect_ratio=1, size=(800, 800))

docs/src/examples/time_series_sublevel.jl

@@ -29,7 +29,7 @@ plot(x, y, xlab="x", ylab="y", legend=false, title="Data")
 # height of an adjacent local maximum. An interval with infinite persistence represents the
 # global minimum.
 
-res = ripserer(Cubical(y), dim_max=0, representatives=true)[1]
+res = ripserer(Cubical(y), dim_max=0, reps=true)[1]
 plot(res)
 
 # We notice there is a lot of noise on the persistence diagram. We can filter it out by

docs/src/index.md

@@ -14,12 +14,19 @@ introduction](https://towardsdatascience.com/persistent-homology-with-examples-1
 
 See the Examples for further info.
 
+Ripserer was created by Matija Čufar. If you used this software in your project, or if
+you have any comments, questions or suggestions, feel free to contact me at
+[matijacufar@gmail.com](mailto:matijacufar@gmail.com).
+
+While this package is fully functional, it is still in development and should not be
+considered stable. Interfaces and internals may still change.
+
 ## Installation
 
-This package is registered. To install it, run the following.
+This package is registered. To install it, simply run the following.
 
 ```julia
-julia> using Pkg
+julia> import Pkg
 julia> Pkg.add("Ripserer")
 ```
 
@@ -50,12 +57,12 @@ Much like Ripser, Ripserer uses the following optimizations to achieve its speed
 * Skip apparent and emergent persistence pairs.
 
 For a detailed description of the algorithm, please see the
-[original article](https://arxiv.org/abs/1908.02518).
+[Ulrich Bauer's article on Ripser](https://arxiv.org/abs/1908.02518).
 
 In general, the performance of Ripserer is very close to
-[Ripser](https://github.com/Ripser/ripser). Depending on the data set, one or the other may
-be faster and the differences are usually small. There are no official benchmarks, because I
-have found benchmarking on my computer or a CI system to be too unreliable.
+[Ripser](https://github.com/Ripser/ripser), within around 30%. Depending on the data set,
+one or the other may be faster. There are no official benchmarks, because I have found
+benchmarking on my computer or a CI system to be too unreliable.
 
 ## Extending
 
@@ -78,3 +85,4 @@ I would like to thank:
 * [@ctralie](https://github.com/ctralie) and [@sauln](https://github.com/sauln) for creating
   [ripser.py](https://github.com/scikit-tda/ripser.py/) which has been a source of
   inspiration.
+* Žiga Virk, for giving ideas and helping with the theoretical side of things.

src/Ripserer.jl

@@ -22,13 +22,13 @@ using ProgressMeter
 using RecipesBase
 using TupleTools
 
-# reexport basics from PersistenceDiagrams
+# reexporting basics these makes Ripserer usable without having to import another package.
 import PersistenceDiagrams: birth, threshold
 export PersistenceDiagram, PersistenceInterval
 export birth, death, persistence, representative, barcode
 
 export Mod
-export AbstractSimplex, diam, coface_type, vertices, coboundary, dim
+export AbstractSimplex, diam, coface_type, face_type, vertices, coboundary, boundary, dim
 export IndexedSimplex, index, Simplex
 export AbstractFiltration, n_vertices, edges, threshold
 export AbstractFlagFiltration, Rips, SparseRips
@@ -46,8 +46,10 @@ include("ripsfiltration.jl")
 include("cubical.jl")
 
 include("chainelement.jl")
-include("reduction_structs.jl")
-include("reduction.jl")
+include("zerodimensional.jl")
+include("reductionmatrix.jl")
+include("ripserer.jl")
+
 include("simplexrecipes.jl")
 
 end

src/abstractsimplex.jl

@@ -12,9 +12,11 @@ accessed by [`dim`](@ref).
 * [`Base.sign(::AbstractSimplex)`](@ref)
 * [`Base.:-(::AbstractSimplex)`](@ref)
 * `Base.isless(::AbstractSimplex, ::AbstractSimplex)`
-* [`coface_type(::AbstractSimplex)`](@ref)
 * [`vertices(::AbstractSimplex)`](@ref)
+* [`coface_type(::AbstractSimplex)`](@ref)
 * [`coboundary(::Any, ::AbstractSimplex)`](@ref)
+* [`face_type(::AbstractSimplex)`](@ref) only required for homology.
+* [`boundary(::Any, ::AbstractSimplex)`](@ref) only required for homology.
 """
 abstract type AbstractSimplex{D, T} end
 
@@ -75,7 +77,7 @@ Base.hash(sx::AbstractSimplex, h::UInt64) = hash(vertices(sx), hash(diam(sx), h)
     coface_type(::AbstractSimplex)
     coface_type(::Type{<:AbstractSimplex})
 
-Get the type of coface a simplex hax. For a `D`-dimensional simplex, this is usually its
+Get the type of simplex's coface. For a `D`-dimensional simplex, this is usually its
 `D+1`-dimensional counterpart. Only the method for the type needs to be implemented.
 
 ```jldoctest
@@ -88,6 +90,23 @@ Simplex{3, Float64, Int}
 """
 coface_type(sx::AbstractSimplex) = coface_type(typeof(sx))
 
+"""
+    face_type(::AbstractSimplex)
+    face_type(::Type{<:AbstractSimplex})
+
+Get the type of a simplex's face. For a `D`-dimensional simplex, this is usually its
+`D-1`-dimensional counterpart. Only the method for the type needs to be implemented.
+
+```jldoctest
+face_type(Simplex{2}((3, 2, 1), 3.2))
+
+# output
+
+Simplex{1, Float64, Int}
+```
+"""
+face_type(sx::AbstractSimplex) = face_type(typeof(sx))
+
 """
     dim(::AbstractSimplex)
     dim(::Type{<:AbstractSimplex})

src/cubical.jl

@@ -46,7 +46,8 @@ index(csx::Cubelet) = csx.index
     return :(vertices(index(csx), Val($(2^D))))
 end
 
-coface_type(::Type{Cubelet{D, T, I}}) where {D, T, I} = Cubelet{D+1, T, I}
+coface_type(::Type{Cubelet{D, T, I}}) where {D, T, I} = Cubelet{D + 1, T, I}
+face_type(::Type{Cubelet{D, T, I}}) where {D, T, I} = Cubelet{D - 1, T, I}
 
 """
     Cubical{T, N} <: AbstractFiltration{T, <:Cubelet{0, T}}
@@ -173,3 +174,39 @@ function Base.iterate(cc::CubeletCoboundary{A, N, C}, (dim, dir)=(1, 1)) where {
     all_vertices = TupleTools.sort(all_vertices, rev=true)
     return coface_type(C)(-dir * index(all_vertices), diameter), (dim, dir)
 end
+
+struct CubeletBoundary{D, C<:Cubelet, N, F<:Cubical{<:Any, N}, K}
+    filtration ::F
+    cubelet    ::C
+    vertices   ::NTuple{K, CartesianIndex{N}}
+end
+
+function CubeletBoundary(
+    filtration::F, cubelet::C
+) where {N, D, F<:Cubical{<:Any, N}, C<:Cubelet{D}}
+
+    K = 2^D
+    vxs = map(v -> CartesianIndices(filtration)[v], vertices(cubelet))
+    return CubeletBoundary{D, C, N, F, K}(filtration, cubelet, vxs)
+end
+
+boundary(filtration::Cubical, cubelet::Cubelet) = CubeletBoundary(filtration, cubelet)
+
+# Idea: split the cube in each dimension, returning two halves depending on dir.
+function Base.iterate(cb::CubeletBoundary{D, C}, (dim, dir)=(1, 1)) where {D, C}
+    if dim > D
+        return nothing
+    else
+        lin_vertices = map(v -> LinearIndices(cb.filtration)[v],
+                           TupleTools.sort(cb.vertices, by=v -> v[dim], rev=true))
+        if dir == 1
+            new_vertices = lin_vertices[1:end÷2]
+            diameter = diam(cb.filtration, new_vertices)
+            return face_type(C)(index(new_vertices), diameter), (dim, -dir)
+        else
+            new_vertices = lin_vertices[end÷2+1:end]
+            diameter = diam(cb.filtration, new_vertices)
+            return face_type(C)(-index(new_vertices), diameter), (dim + 1, 1)
+        end
+    end
+end