New Way Of Handling Representatives (#55)

* New Way Of Handling Representatives

* Fix critical bug in cubical

Project.toml

@@ -1,7 +1,7 @@
 name = "Ripserer"
 uuid = "aa79e827-bd0b-42a8-9f10-2b302677a641"
 authors = ["mtsch <matijacufar@gmail.com>"]
-version = "0.8.2"
+version = "0.9.0"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
@@ -24,7 +24,7 @@ DataStructures = "0.17"
 Distances = "0.8, 0.9"
 Hungarian = "0.6"
 IterTools = "1"
-PersistenceDiagrams = "0.4"
+PersistenceDiagrams = "0.5"
 ProgressMeter = "1"
 RecipesBase = "1"
 StaticArrays = "0.12"

docs/src/examples/image_sublevel.jl

@@ -33,7 +33,7 @@ heatmap(blackhole_grayscale, aspect_ratio=1, size=(800, 800))
 #     be careful with computing representatives because collecting them for thousands of
 #     persistence intervals can take a while.
 
-result_min = ripserer(Cubical(blackhole_grayscale))
+result_min = ripserer(Cubical{128}(blackhole_grayscale))
 
 # We plot the diagram with the `persistence` argument. This plots persistence vs birth
 # instead of death vs birth.
@@ -44,15 +44,15 @@ plot(result_min, persistence=true)
 # remove the noise by supplying the `cutoff` argument to `ripserer`. This removes all
 # intervals with persistence strictly smaller than cutoff.
 
-result_min = ripserer(Cubical(blackhole_grayscale), cutoff=0.1)
+result_min = ripserer(Cubical{128}(blackhole_grayscale), cutoff=0.1)
 
 #
 
 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, reps=true)
+result_min = ripserer(Cubical{128}(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, reps=true)
+result_max = ripserer(Cubical{128}(-blackhole_grayscale), cutoff=0.1, reps=true)
 plot(result_max)
 
 heatmap(blackhole_grayscale, aspect_ratio=1, size=(800, 800))

src/Ripserer.jl

@@ -25,17 +25,19 @@ using TupleTools
 
 # 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, face_type, vertices, coboundary, boundary, dim
-export IndexedSimplex, index, Simplex
-export AbstractFiltration, n_vertices, edges, threshold
-export AbstractFlagFiltration, Rips, SparseRips
-export Cubelet, Cubical
-export simplex, coefficient
-export ripserer
+export
+    PersistenceDiagram, PersistenceInterval, RepresentativeInterval,
+    birth, death, persistence, representative, birth_simplex, death_simplex, barcode
+
+export
+    Mod,
+    AbstractSimplex, diam, coface_type, face_type, vertices, coboundary, boundary, dim,
+    IndexedSimplex, index, Simplex,
+    AbstractFiltration, n_vertices, edges, threshold,
+    AbstractFlagFiltration, Rips, SparseRips,
+    Cubelet, Cubical,
+    simplex, coefficient,
+    ripserer
 
 include("primefield.jl")
 

src/chainelement.jl

@@ -59,6 +59,10 @@ function Base.show(io::IO, ce::AbstractChainElement)
     print(io, simplex(ce), " => ", coefficient(ce))
 end
 
+function Base.convert(::Type{C}, simplex::S) where {S, C<:AbstractChainElement{S}}
+    C(simplex)
+end
+
 """
     chain_element_type(simplex, coefficient)
 

src/cubical.jl

@@ -154,9 +154,7 @@ function Base.iterate(
 ) where {A, N, I, C}
     # If not all indices in a given dimension are equal, we can't create a coface by
     # expanding in that direction.
-    diameter = missing
-    new_vertices = cc.vertices
-    while ismissing(diameter)
+    while true
         while dim ≤ N && !all_equal_in_dim(dim, cc.vertices)
             dim += one(I)
         end
@@ -170,14 +168,14 @@ function Base.iterate(
 
         dim += I(dir == -1)
         dir *= -one(I)
-    end
 
-    if ismissing(diameter)
-        return nothing
-    else
-        all_vertices = sort(map(v -> I(LinearIndices(cc.filtration)[v]),
-                                vcat(cc.vertices, new_vertices)), rev=true)
-        return coface_type(C)(-dir * index(all_vertices), diameter), (dim, dir)
+        if !ismissing(diameter)
+            all_vertices = sort(map(v -> I(LinearIndices(cc.filtration)[v]),
+                                    vcat(cc.vertices, new_vertices)), rev=true)
+            if A || all_vertices[end÷2+1:end] == LinearIndices(cc.filtration)[cc.vertices]
+                return coface_type(C)(-dir * index(all_vertices), diameter), (dim, dir)
+            end
+        end
     end
 end
 

src/reductionmatrix.jl

@@ -309,42 +309,50 @@ function reduce_column!(matrix::ReductionMatrix, column_simplex)
     return pivot
 end
 
-function birth_death(matrix::ReductionMatrix{true}, simplex, pivot)
-    return diam(simplex), isnothing(pivot) ? Inf : diam(pivot)
+function birth_death(matrix::ReductionMatrix{true}, column, pivot)
+    return diam(column), isnothing(pivot) ? Inf : diam(pivot)
 end
-function birth_death(matrix::ReductionMatrix{false}, simplex, pivot)
-    return isnothing(pivot) ? Inf : diam(pivot), diam(simplex)
+function birth_death(matrix::ReductionMatrix{false}, column, pivot)
+    return isnothing(pivot) ? Inf : diam(pivot), diam(column)
 end
 
 # Interval with no representative, (co)homology.
 function interval(
-    ::Type{PersistenceInterval{Nothing}}, matrix, simplex, pivot, cutoff
+    ::Type{PersistenceInterval}, matrix::ReductionMatrix, column, pivot, cutoff
 )
-    birth, death = birth_death(matrix, simplex, pivot)
+    birth, death = birth_death(matrix, column, pivot)
     return death - birth > cutoff ? PersistenceInterval(birth, death) : nothing
 end
 # With representative, cohomology.
 function interval(
-    ::Type{PersistenceInterval{V}}, matrix::ReductionMatrix{true}, simplex, pivot, cutoff
-) where V<:AbstractVector
-    birth, death = birth_death(matrix, simplex, pivot)
+    ::Type{R}, matrix::ReductionMatrix{true}, column, pivot, cutoff
+) where R<:RepresentativeInterval
+    birth, death = birth_death(matrix, column, pivot)
 
     if !isfinite(death)
-        return PersistenceInterval(birth, death, V())
+        return R(
+            PersistenceInterval(birth, death), column, nothing, eltype(matrix.reduced)[]
+        )
     elseif death - birth > cutoff
-        return PersistenceInterval(birth, death, collect(matrix.reduced[pivot]))
+        return R(
+            PersistenceInterval(birth, death),
+            column, simplex(pivot), collect(matrix.reduced[pivot])
+        )
     else
         return nothing
     end
 end
 # With representative, homology.
 function interval(
-    ::Type{PersistenceInterval{V}}, matrix::ReductionMatrix{false}, simplex, pivot, cutoff
-) where V<:AbstractVector
-    birth, death = birth_death(matrix, simplex, pivot)
+    ::Type{R}, matrix::ReductionMatrix{false}, column, pivot, cutoff
+) where R<:RepresentativeInterval
+    birth, death = birth_death(matrix, column, pivot)
 
     if death - birth > cutoff
-        return PersistenceInterval(birth, death, move!(matrix.working_boundary))
+        return R(
+            PersistenceInterval(birth, death),
+            simplex(pivot), column, move!(matrix.working_boundary)
+        )
     else
         return nothing
     end
@@ -355,9 +363,16 @@ function compute_intervals!(
 ) where {Co, reps}
     if reps
         representative_type = Co ? simplex_element(matrix) : face_element(matrix)
-        intervals = PersistenceInterval{Vector{representative_type}}[]
+        critical_birth_type = simplex_type(matrix.filtration, dim(matrix))
+        critical_death_type = simplex_type(matrix.filtration, dim(matrix) + 1)
+        intervals = RepresentativeInterval{
+            PersistenceInterval,
+            critical_birth_type,
+            Union{critical_death_type, Nothing},
+            Vector{representative_type},
+        }[]
     else
-        intervals = PersistenceInterval{Nothing}[]
+        intervals = PersistenceInterval[]
     end
     if progress
         progbar = Progress(

src/simplexrecipes.jl

@@ -1,9 +1,3 @@
-# simplex plots ========================================================================== #
-const SxVector{D} = AbstractVector{<:AbstractSimplex{D}}
-const IntervalWithRep{D} = PersistenceInterval{
-    <:AbstractVector{<:Pair{<:AbstractSimplex{D}, <:Any}}
-}
-
 """
     index_data(indices, args...)
 
@@ -35,7 +29,7 @@ end
 
 """
     plottable(sx::AbstractSimplex, args...)
-    plottable(sx::PersistenceInterval, args...)
+    plottable(sx::RepresentativeInterval, args...)
     plottable(sx::Vector{<:AbstractSimplex}, args...)
 
 Turn `sx` into a series that can be plotted. `args...`, should contain the data which tells
@@ -44,23 +38,26 @@ the plot where to place simplex vertices.
 function plottable(sx::AbstractSimplex, args...)
     return plottable([sx], args...)
 end
-function plottable(int::PersistenceInterval, args...)
+function plottable(int::RepresentativeInterval, args...)
     return plottable(simplex.(representative(int)), args...)
 end
-function plottable(int::PersistenceInterval{Nothing}, args...)
-    error("interval has no representative. Run `ripserer` with `representatives=true`")
+function plottable(rep::AbstractVector{<:AbstractChainElement}, args...)
+    return plottable(simplex.(rep), args...)
+end
+function plottable(element::AbstractChainElement, args...)
+    return plottable(simplex(element), args...)
 end
-function plottable(sxs::SxVector{0}, args...)
+function plottable(sxs::AbstractVector{<:AbstractSimplex{0}}, args...)
     indices = only.(vertices.(sxs))
     return index_data(indices, args...), [:seriestype => :scatter], 0
 end
-function plottable(sxs::SxVector{1}, args...)
+function plottable(sxs::AbstractVector{<:AbstractSimplex{1}}, args...)
     indices = mapreduce(vcat, vertices.(sxs)) do (u, v)
         [u, v, 0]
     end
     return index_data(indices, args...), [:seriestype => :path], 1
 end
-function plottable(sxs::SxVector{D}, args...) where D
+function plottable(sxs::AbstractVector{<:AbstractSimplex{D}}, args...) where D
     indices = mapreduce(vcat, vertices.(sxs)) do vs
         idxs = Int[]
         for (u, v, w) in subsets(vs, Val(3))
@@ -75,16 +72,32 @@ end
 function apply_threshold(sx::AbstractSimplex, thresh, thresh_strict)
     return diam(sx) ≤ thresh && diam(sx) < thresh_strict ? sx : nothing
 end
-function apply_threshold(sxs::SxVector, thresh, thresh_strict)
+function apply_threshold(sxs::AbstractVector{<:AbstractSimplex}, thresh, thresh_strict)
     return filter(sx -> diam(sx) ≤ thresh && diam(sx) < thresh_strict, sxs)
 end
-function apply_threshold(int::PersistenceInterval, thresh, thresh_strict)
-    reps = filter(sx -> diam(sx) ≤ thresh && diam(sx) < thresh_strict, representative(int))
-    return PersistenceInterval(birth(int), death(int), reps)
+function apply_threshold(els::AbstractVector{<:AbstractChainElement}, thresh, thresh_strict)
+    return apply_threshold(simplex.(els), thresh, thresh_strict)
+end
+function apply_threshold(int::PersistenceDiagrams.AbstractInterval, thresh, thresh_strict)
+    error("interval has no representative. Run `ripserer` with `representatives=true`")
+end
+function apply_threshold(int::RepresentativeInterval, thresh, thresh_strict)
+    return apply_threshold(representative(int), thresh, thresh_strict)
 end
 
-@recipe function f(sx::Union{AbstractSimplex, SxVector, PersistenceInterval}, args...;
-                   threshold=Inf, threshold_strict=Inf)
+@recipe function f(
+    sx::Union{
+        AbstractSimplex,
+        AbstractChainElement,
+        AbstractVector{<:AbstractSimplex},
+        AbstractVector{<:AbstractChainElement},
+        PersistenceDiagrams.AbstractInterval,
+        RepresentativeInterval,
+    },
+    args...;
+    threshold=Inf,
+    threshold_strict=Inf
+)
     sx = apply_threshold(sx, threshold, threshold_strict)
     isnothing(sx) && return ()
     series, attrs, D = plottable(sx, args...)