From bdc2708fe87f830cd6208314bcf551d935146b3a Mon Sep 17 00:00:00 2001 From: mtsch Date: Mon, 15 Jun 2020 20:04:38 +1200 Subject: [PATCH] New Way Of Handling Representatives (#55) * New Way Of Handling Representatives * Fix critical bug in cubical --- Project.toml | 4 +- docs/src/examples/image_sublevel.jl | 8 +- src/Ripserer.jl | 24 ++-- src/chainelement.jl | 4 + src/cubical.jl | 18 ++- src/reductionmatrix.jl | 49 ++++--- src/simplexrecipes.jl | 51 ++++--- src/zerodimensional.jl | 60 +++++--- test/cubical.jl | 44 +++++- test/plotting.jl | 7 +- test/ripserer.jl | 208 ++++++++++++++++++---------- 11 files changed, 318 insertions(+), 159 deletions(-) diff --git a/Project.toml b/Project.toml index 053d72e5..90350a00 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Ripserer" uuid = "aa79e827-bd0b-42a8-9f10-2b302677a641" authors = ["mtsch "] -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" diff --git a/docs/src/examples/image_sublevel.jl b/docs/src/examples/image_sublevel.jl index 6c44bd28..f114c2ab 100644 --- a/docs/src/examples/image_sublevel.jl +++ b/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,7 +44,7 @@ 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) # @@ -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, 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)) diff --git a/src/Ripserer.jl b/src/Ripserer.jl index acbb3913..e422fd8f 100644 --- a/src/Ripserer.jl +++ b/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") diff --git a/src/chainelement.jl b/src/chainelement.jl index 46edbcd2..9d602a9c 100644 --- a/src/chainelement.jl +++ b/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) diff --git a/src/cubical.jl b/src/cubical.jl index 7a624762..dea019ed 100644 --- a/src/cubical.jl +++ b/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 diff --git a/src/reductionmatrix.jl b/src/reductionmatrix.jl index f34b3de4..236d5a16 100644 --- a/src/reductionmatrix.jl +++ b/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( diff --git a/src/simplexrecipes.jl b/src/simplexrecipes.jl index b43b4cb0..7c4e9cdc 100644 --- a/src/simplexrecipes.jl +++ b/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...) diff --git a/src/zerodimensional.jl b/src/zerodimensional.jl index 5d12c73f..bff33c16 100644 --- a/src/zerodimensional.jl +++ b/src/zerodimensional.jl @@ -65,16 +65,31 @@ end birth(dset::DisjointSetsWithBirth, i) = dset.births[i] -""" - zeroth_representative(dset, vertex, reps, CE, V) +function birth_death(dset::DisjointSetsWithBirth, vertex, edge) + return birth(dset, vertex), isnothing(edge) ? Inf : diam(edge) +end -Collect the zero dimensional representative of `vertex`. -""" -function zeroth_representative(filtration, dset, vertex, reps, CE, V) - if reps - return map(find_leaves!(dset, vertex)) do u - CE(V((u,), birth(filtration, u))) +function interval( + ::Type{R}, dset::DisjointSetsWithBirth, filtration, vertex, edge, cutoff +) where {V, R<:RepresentativeInterval{PersistenceInterval, V}} + birth, death = birth_death(dset, vertex, edge) + if death - birth > cutoff + rep = map(find_leaves!(dset, vertex)) do u + V((u,), Ripserer.birth(filtration, u)) end + birth_vertex = findfirst(r -> diam(r) == birth, rep) + return R(PersistenceInterval(birth, death), V((birth_vertex,), birth), edge, rep) + else + return nothing + end +end + +function interval( + ::Type{PersistenceInterval}, dset::DisjointSetsWithBirth, filtration, vertex, edge, cutoff +) + birth, death = birth_death(dset, vertex, edge) + if death - birth > cutoff + return PersistenceInterval(birth, death) else return nothing end @@ -97,9 +112,14 @@ function zeroth_intervals( CE = chain_element_type(V, F) dset = DisjointSetsWithBirth([birth(filtration, v) for v in 1:n_vertices(filtration)]) if reps - intervals = PersistenceInterval{Vector{CE}}[] + intervals = RepresentativeInterval{ + PersistenceInterval, + vertex_type(filtration), + Union{edge_type(filtration), Nothing}, + Vector{CE}, + }[] else - intervals = PersistenceInterval{Nothing}[] + intervals = PersistenceInterval[] end to_skip = edge_type(filtration)[] to_reduce = edge_type(filtration)[] @@ -107,30 +127,28 @@ function zeroth_intervals( if progress progbar = Progress(length(simplices), desc="Computing 0d intervals... ") end - for sx in simplices - u, v = index.(vertices(sx)) + for edge in simplices + u, v = index.(vertices(edge)) i = find_root!(dset, u) j = find_root!(dset, v) if i ≠ j # According to the elder rule, the vertex with the lower birth will fall # into a later interval. - dead = birth(dset, i) > birth(dset, j) ? i : j - if diam(sx) - birth(dset, dead) > cutoff - representative = zeroth_representative(filtration, dset, dead, reps, CE, V) - interval = PersistenceInterval(birth(dset, dead), diam(sx), representative) - push!(intervals, interval) + v = birth(dset, i) > birth(dset, j) ? i : j + int = interval(eltype(intervals), dset, filtration, v, edge, cutoff) + if !isnothing(int) + push!(intervals, int) end union!(dset, i, j) - push!(to_skip, sx) + push!(to_skip, edge) else - push!(to_reduce, sx) + push!(to_reduce, edge) end progress && next!(progbar) end for v in 1:n_vertices(filtration) if find_root!(dset, v) == v - representative = zeroth_representative(filtration, dset, v, reps, CE, V) - push!(intervals, PersistenceInterval(birth(dset, v), Inf, representative)) + push!(intervals, interval(eltype(intervals), dset, filtration, v, nothing, 0)) end end reverse!(to_reduce) diff --git a/test/cubical.jl b/test/cubical.jl index 45f877a0..fe6144a8 100644 --- a/test/cubical.jl +++ b/test/cubical.jl @@ -1,4 +1,6 @@ using Ripserer + +using Compat using StaticArrays using Ripserer: all_equal_in_dim @@ -104,7 +106,47 @@ end @test isempty(cocob) end - @testset "3d" begin + @testset "3d values" begin + data = zeros(2, 2, 2) + data[:, 1, 1] .= 1 + data[:, 2, 1] .= 2 + data[:, 1, 2] .= 3 + data[:, 2, 2] .= 4 + + @testset "Cubelet{1} coboundary" begin + cob = Cubelet{2, Float64, Int}[] + for cube in edges(Cubical(data)) + for c in coboundary(Cubical(data), cube) + @test issubset(cube, c) + @test diam(cube) ≤ diam(c) + end + end + end + @testset "Cubelet{1} coboundary, all=false" begin + cob = Cubelet{2, Float64, Int}[] + for cube in edges(Cubical(data)) + for c in coboundary(Cubical(data), cube, Val(false)) + push!(cob, c) + end + end + @test length(cob) == 6 + @test allunique(cob) + end + @testset "Cubelet{2} coboundary" begin + for cube in [ + Cubelet{2}((1, 2, 3, 4), 2.0), + Cubelet{2}((1, 2, 5, 6), 4.0), + Cubelet{2}((1, 3, 5, 7), 3.0), + Cubelet{2}((2, 4, 6, 8), 4.0), + Cubelet{2}((3, 4, 7, 8), 4.0), + Cubelet{2}((5, 6, 7, 8), 4.0), + ] + coface = only(coboundary(Cubical(data), cube)) + @test abs(coface) == Cubelet{3}((8, 7, 6, 5, 4, 3, 2, 1), 4.0) + end + end + end + @testset "3d counts" begin cob = Cubelet{2, Int, Int}[] flt = Cubical(data3d) cub = Cubelet{1}((14, 13), 1) diff --git a/test/plotting.jl b/test/plotting.jl index 5b6018d3..b7450eb9 100644 --- a/test/plotting.jl +++ b/test/plotting.jl @@ -43,13 +43,14 @@ series(args...; kwargs...) = data = collect(1:100) for dim in 0:3 sx = Simplex{dim}(1, 1) + fsx = Simplex{dim + 1}(1, 1) @test length(series(sx, data)) == 1 @test length(series([sx], data)) == 1 @test length(series([sx], data)) == 1 - @test length(series(PersistenceInterval( - 1.0, 1.0, [ChainElement{typeof(sx), typeof(1//1)}(sx, 1//1)] + @test length(series(RepresentativeInterval( + 1.0, 1.0, sx, fsx, [ChainElement{typeof(sx), typeof(1//1)}(sx, 1//1)] ), data)) == 1 - @test_throws ErrorException series(PersistenceInterval(1.0, 1.0, nothing), data) + @test_throws ErrorException series(PersistenceInterval(1.0, 1.0), data) end end end diff --git a/test/ripserer.jl b/test/ripserer.jl index 127c6835..f7133a9d 100644 --- a/test/ripserer.jl +++ b/test/ripserer.jl @@ -170,55 +170,90 @@ end end end -@testset "reps" begin - _, d1, d2 = ripserer(projective_plane, dim_max=2, reps=true) - - @test simplex.(representative(only(d1))) == [ - Simplex{1}((11, 10), 1), - Simplex{1}((10, 7), 1), - Simplex{1}((10, 6), 1), - Simplex{1}((8, 1), 1), - Simplex{1}((7, 3), 1), - Simplex{1}((7, 1), 1), - Simplex{1}((6, 2), 1), - Simplex{1}((5, 1), 1), - Simplex{1}((2, 1), 1), - ] - @test coefficient.(representative(only(d1))) == fill(Mod{2}(1), 9) - @test simplex.(representative(only(d2))) == [Simplex{2}((6, 2, 1), 1)] - @test coefficient.(representative(only(d2))) == [Mod{2}(1)] -end - -@testset "reps - types" begin - d0, d1, d2, d3 = ripserer(cycle, dim_max=3, reps=true) - @test eltype(d0) <: PersistenceInterval{ - <:Vector{<:PackedElement{Simplex{0, Int, Int}, Mod{2}}}} - @test eltype(d1) <: PersistenceInterval{ - <:Vector{<:PackedElement{Simplex{1, Int, Int}, Mod{2}}}} - @test eltype(d2) <: PersistenceInterval{ - <:Vector{<:PackedElement{Simplex{2, Int, Int}, Mod{2}}}} - @test eltype(d3) <: PersistenceInterval{ - <:Vector{<:PackedElement{Simplex{3, Int, Int}, Mod{2}}}} - - d0, d1, d2, d3 = ripserer(cycle, dim_max=3, reps=true, - field_type=Rational{Int}) - @test eltype(d0) ≡ PersistenceInterval{ - Vector{ChainElement{Simplex{0, Int, Int}, Rational{Int}}}} - @test eltype(d1) ≡ PersistenceInterval{ - Vector{ChainElement{Simplex{1, Int, Int}, Rational{Int}}}} - @test eltype(d2) ≡ PersistenceInterval{ - Vector{ChainElement{Simplex{2, Int, Int}, Rational{Int}}}} - @test eltype(d3) ≡ PersistenceInterval{ - Vector{ChainElement{Simplex{3, Int, Int}, Rational{Int}}}} -end - -@testset "Representative of infinite interval." begin - _, d1 = ripserer( - cycle, dim_max=1, reps=true, threshold=1, field_type=Rational{Int} - ) - @test representative(only(d1)) == ChainElement{ - Simplex{0, Int, Int}, Rational{Int} - }[] +@testset "Representatives" begin + @testset "example from ripser" begin + _, d1, d2 = ripserer(projective_plane, dim_max=2, reps=true) + + @test simplex.(representative(only(d1))) == [ + Simplex{1}((11, 10), 1), + Simplex{1}((10, 7), 1), + Simplex{1}((10, 6), 1), + Simplex{1}((8, 1), 1), + Simplex{1}((7, 3), 1), + Simplex{1}((7, 1), 1), + Simplex{1}((6, 2), 1), + Simplex{1}((5, 1), 1), + Simplex{1}((2, 1), 1), + ] + @test coefficient.(representative(only(d1))) == fill(Mod{2}(1), 9) + @test simplex.(representative(only(d2))) == [Simplex{2}((6, 2, 1), 1)] + @test coefficient.(representative(only(d2))) == [Mod{2}(1)] + end + @testset "types" begin + d0, d1, d2, d3 = ripserer(cycle, dim_max=3, reps=true) + @test eltype(d0) <: RepresentativeInterval{ + PersistenceInterval, + Simplex{0, Int, Int}, + Union{Nothing, Simplex{1, Int, Int}}, + <:Vector{<:PackedElement{Simplex{0, Int, Int}, Mod{2}}}} + @test eltype(d1) <: RepresentativeInterval{ + PersistenceInterval, + Simplex{1, Int, Int}, + Union{Nothing, Simplex{2, Int, Int}}, + <:Vector{<:PackedElement{Simplex{1, Int, Int}, Mod{2}}}} + @test eltype(d2) <: RepresentativeInterval{ + PersistenceInterval, + Simplex{2, Int, Int}, + Union{Nothing, Simplex{3, Int, Int}}, + <:Vector{<:PackedElement{Simplex{2, Int, Int}, Mod{2}}}} + @test eltype(d3) <: RepresentativeInterval{ + PersistenceInterval, + Simplex{3, Int, Int}, + Union{Nothing, Simplex{4, Int, Int}}, + <:Vector{<:PackedElement{Simplex{3, Int, Int}, Mod{2}}}} + + d0, d1, d2, d3 = ripserer(cycle, dim_max=3, reps=true, + field_type=Rational{Int}) + @test eltype(d0) ≡ RepresentativeInterval{ + PersistenceInterval, + Simplex{0, Int, Int}, + Union{Nothing, Simplex{1, Int, Int}}, + Vector{ChainElement{Simplex{0, Int, Int}, Rational{Int}}}} + @test eltype(d1) ≡ RepresentativeInterval{ + PersistenceInterval, + Simplex{1, Int, Int}, + Union{Nothing, Simplex{2, Int, Int}}, + Vector{ChainElement{Simplex{1, Int, Int}, Rational{Int}}}} + @test eltype(d2) ≡ RepresentativeInterval{ + PersistenceInterval, + Simplex{2, Int, Int}, + Union{Nothing, Simplex{3, Int, Int}}, + Vector{ChainElement{Simplex{2, Int, Int}, Rational{Int}}}} + @test eltype(d3) ≡ RepresentativeInterval{ + PersistenceInterval, + Simplex{3, Int, Int}, + Union{Nothing, Simplex{4, Int, Int}}, + Vector{ChainElement{Simplex{3, Int, Int}, Rational{Int}}}} + end + @testset "infinite interval" begin + _, d1 = ripserer( + cycle, dim_max=1, reps=true, threshold=1, field_type=Rational{Int} + ) + @test representative(only(d1)) == ChainElement{ + Simplex{0, Int, Int}, Rational{Int} + }[] + end + @testset "critical simplices" begin + t = torus(100) + result = ripserer(t, reps=true, threshold=0.5) + for diag in result + @test birth.(diag) == diam.(birth_simplex.(diag)) + finite = filter(isfinite, diag) + @test death.(finite) == diam.(death_simplex.(finite)) + infinite = filter(!isfinite, diag) + @test all(isnothing, death_simplex.(infinite)) + end + end end @testset "Sublevel using SparseRips" begin @@ -246,35 +281,66 @@ end @test sort(maxs) == [1.0, 2.0] end -@testset "2D image sublevel using Cubical" begin - data = [0 0 0 0 0; - 0 2 2 2 0; - 0 2 1 2 0; - 0 2 2 2 0; - 0 0 0 0 0] +@testset "Cubical" begin + @testset "2D image" begin + data = [0 0 0 0 0; + 0 2 2 2 0; + 0 2 1 2 0; + 0 2 2 2 0; + 0 0 0 0 0] - d0, d1, d2 = ripserer(Cubical(data), reps=true, dim_max=2) + d0, d1, d2 = ripserer(Cubical(data), reps=true, dim_max=2) - @test d0 == [(0, Inf), (1, 2)] - @test d1 == [(0, 2)] - @test isempty(d2) + @test d0 == [(0, Inf), (1, 2)] + @test d1 == [(0, 2)] + @test isempty(d2) - @test vertices.(representative(d0[1])) == [SVector(i) for i in 1:length(data)] - @test vertices(only(representative(d0[2]))) == SVector(13) + @test vertices.(representative(d0[1])) == [SVector(i) for i in 1:length(data)] + @test vertices(only(representative(d0[2]))) == SVector(13) + end + @testset "3D image" begin + # Cube with hole in the middle. + data = zeros(5, 5, 5) + data[2, 2:4, 2:4] .= 1 + data[3, :, :] .= [0 0 0 0 0; 0 1 1 1 0; 0 1 0 1 0; 0 1 1 1 0; 0 0 0 0 0] + data[4, 2:4, 2:4] .= 1 + + d0, d1, d2 = ripserer(Cubical(data), dim_max=2) + + @test d0 == [(0, 1.0), (0, Inf)] + @test isempty(d1) + @test d2 == [(0, 1)] + end end @testset "Persistent homology." begin - res_hom = ripserer(cycle, cohomology=false, dim_max=3) - res_coh = ripserer(cycle, cohomology=true, dim_max=3) - - @test res_hom == res_coh - - res_hom = ripserer(cycle, cohomology=false, reps=true, dim_max=3) - @test vertices.(simplex.(representative(res_hom[2][1]))) == sort!(vcat( - [SVector(i+1, i) for i in 1:17], [SVector(18, 1)] - )) + @testset "same diagrams" begin + res_hom = ripserer(cycle, cohomology=false, dim_max=3) + res_coh = ripserer(cycle, cohomology=true, dim_max=3) - @test_broken ripserer(cycle, cohomology=false, threshold=2)[2][1] == (1.0, Inf) + @test res_hom == res_coh + end + @testset "same critical simplices" begin + # Add some noise because critical simplices might be different if values are exactly + # the same. + cyc = cycle .+ 0.01 .* rand_dist_matrix(18) + res_hom = ripserer(cyc, cohomology=false, dim_max=3, reps=true) + res_coh = ripserer(cyc, cohomology=true, dim_max=3, reps=true) + + for i in 1:4 + @test birth_simplex.(res_hom[i]) == birth_simplex.(res_coh[i]) + @test death_simplex.(res_hom[i]) == death_simplex.(res_coh[i]) + end + end + @testset "representative cycle" begin + res_hom = ripserer(cycle, cohomology=false, reps=true, dim_max=3) + @test vertices.(simplex.(representative(res_hom[2][1]))) == sort!(vcat( + [SVector(i+1, i) for i in 1:17], [SVector(18, 1)] + )) + end + @testset "infinite intervals" begin + @test_broken ripserer(cycle, cohomology=false, threshold=2)[2][1] == (1.0, Inf) + end end @testset "Only print to stderr and only when progress is enabled." begin