diff --git a/Project.toml b/Project.toml index 6e062118..daefa7fb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Ripserer" uuid = "aa79e827-bd0b-42a8-9f10-2b302677a641" authors = ["mtsch "] -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"] diff --git a/docs/src/examples/basics.jl b/docs/src/examples/basics.jl index f4f3f7b7..7c1da353 100644 --- a/docs/src/examples/basics.jl +++ b/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") diff --git a/docs/src/examples/cocycles.jl b/docs/src/examples/cocycles.jl index fd1406f8..a2e11cac 100644 --- a/docs/src/examples/cocycles.jl +++ b/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", diff --git a/docs/src/examples/image_sublevel.jl b/docs/src/examples/image_sublevel.jl index ed6187fe..6c44bd28 100644 --- a/docs/src/examples/image_sublevel.jl +++ b/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)) diff --git a/docs/src/examples/time_series_sublevel.jl b/docs/src/examples/time_series_sublevel.jl index 13fbb9e6..cf5b3e0b 100644 --- a/docs/src/examples/time_series_sublevel.jl +++ b/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 diff --git a/docs/src/index.md b/docs/src/index.md index cc04c1ea..acc346cb 100644 --- a/docs/src/index.md +++ b/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. diff --git a/src/Ripserer.jl b/src/Ripserer.jl index 8b4d188e..fe34e79f 100644 --- a/src/Ripserer.jl +++ b/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 diff --git a/src/abstractsimplex.jl b/src/abstractsimplex.jl index f8226edf..d30bbd21 100644 --- a/src/abstractsimplex.jl +++ b/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}) diff --git a/src/cubical.jl b/src/cubical.jl index 65906bdb..091a169e 100644 --- a/src/cubical.jl +++ b/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 diff --git a/src/reduction.jl b/src/reduction.jl deleted file mode 100644 index 7e19b0b2..00000000 --- a/src/reduction.jl +++ /dev/null @@ -1,363 +0,0 @@ -""" - ReductionState - -This structure represents the reduction matrix in the current dimension. A new one is -created for every dimension. - -# Fields: - -* `filtration`: the filtration we are analyzing. -* `reduction_matrix`: the reduction matrix. Each column of the matrix records the operations - that were performed when reducing the column. -* `working_column`: the current working column, the column we are currently reducing. -* `reduction_entries`: this is where we record which simplices we added to the working - column. -""" -struct ReductionState{ - Field, - S<:AbstractSimplex, SE<:AbstractChainElement{S, Field}, - C<:AbstractSimplex, CE<:AbstractChainElement{C, Field}, - F, -} - filtration ::F - reduction_matrix ::ReductionMatrix{C, SE} - working_column ::Column{CE} - reduction_entries ::Vector{SE} - - function ReductionState{Field, S}(filtration::F) where {Field, S<:AbstractSimplex, F} - SE = chain_element_type(S, Field) - C = coface_type(S) - CE = chain_element_type(C, Field) - working_column = Column{CE}() - reduction_entries = SE[] - - new{Field, S, SE, C, CE, F}( - filtration, ReductionMatrix{C, SE}(), working_column, reduction_entries - ) - end -end - -simplex_type(rs::ReductionState{<:Any, S}) where S = S -coface_element(rs::ReductionState{<:Any, <:Any, <:Any, <:Any, CE}) where CE = CE - - -""" - initialize!(rs::ReductionState, column_simplex) - -Initialize the columns in `rs`. Empty both `rs.working_column` and `rs.reduction_entries`, -then push the coboundary of `column_simplex` to `rs.working_column`. -""" -function initialize!(rs::ReductionState, column_simplex::AbstractSimplex) - empty!(rs.working_column) - empty!(rs.reduction_entries) - - for coface in coboundary(rs.filtration, column_simplex) - if diam(coface) == diam(column_simplex) && !has_column(rs.reduction_matrix, coface) - empty!(rs.working_column) - return coface_element(rs)(coface) - end - nonheap_push!(rs.working_column, coface) - end - repair!(rs.working_column) - return first(rs.working_column) # There are no duplicates on the heap at this point. -end - -""" - add!(rs::ReductionState, current_pivot) - -Add column in `rs.reduction_matrix` indexed by `current_pivot` and multiplied by the correct -factor to `rs.working_column`. Also record the addition in `rs.reduction_entries`. -""" -function add!(rs::ReductionState, current_pivot) - λ = -coefficient(current_pivot) - for element in rs.reduction_matrix[current_pivot] - push!(rs.reduction_entries, λ * element) - for coface in coboundary(rs.filtration, simplex(element)) - push!(rs.working_column, coface_element(rs)(coface, λ * coefficient(element))) - end - end - return pivot(rs.working_column) -end - -""" - reduce_working_column!(rs::ReductionState, column_simplex, cutoff, reps) - -Reduce the working column by adding other columns to it until it has the lowest pivot or is -reduced. Record it in the reduction matrix and return the persistence intervals with -persistence larger than `cutoff`. If `reps` is `true`, add representative cocycles to the -interval. -""" -function reduce_working_column!( - rs::ReductionState{F, S}, column_simplex, cutoff, reps -) where {F, S} - current_pivot = initialize!(rs, column_simplex) - - while !isnothing(current_pivot) && has_column(rs.reduction_matrix, current_pivot) - current_pivot = add!(rs, current_pivot) - end - if isnothing(current_pivot) - death = Inf - else - insert_column!(rs.reduction_matrix, current_pivot) - push!(rs.reduction_entries, chain_element_type(S, F)(column_simplex)) - append_unique_times!( - rs.reduction_matrix, rs.reduction_entries, inv(coefficient(current_pivot)) - ) - death = diam(simplex(current_pivot)) - end - birth = diam(column_simplex) - - if reps && !isfinite(death) - return PersistenceInterval(birth, death, chain_element_type(S, F)[]) - elseif reps && death - birth > cutoff - representative = collect(rs.reduction_matrix[current_pivot]) - return PersistenceInterval(birth, death, representative) - elseif !isfinite(death) || death - birth > cutoff - return PersistenceInterval(birth, death) - else - return nothing - end -end - -""" - compute_pairs!(rs::ReductionState, columns, cutoff, reps, progress) - -Compute persistence intervals by reducing `columns`, a collection of simplices. Return -`PersistenceDiagram`. Only keep intervals with persistence larger than `cutoff` and return -representative cocycles if `reps` is `true`. If `progress` is set, show a progress bar. -""" -function compute_intervals!( - rs::ReductionState{F, S}, columns, cutoff, reps, progress -) where {S, F} - T = dist_type(rs.filtration) - if reps - intervals = PersistenceInterval{Vector{chain_element_type(S, F)}}[] - else - intervals = PersistenceInterval{Nothing}[] - end - if progress - progbar = Progress(length(columns), desc="Computing $(dim(S))d intervals... ") - end - for column in columns - interval = reduce_working_column!(rs, column, cutoff, reps) - if !isnothing(interval) - push!(intervals, interval) - end - progress && next!(progbar) - end - return sort!( - PersistenceDiagram(dim(eltype(columns)), intervals, threshold(rs.filtration)) - ) -end - -""" - assemble_columns!(rs::ReductionState, unreduced, reduced, progress) - -Assemble columns that need to be reduced in the next dimension. Apply clearing optimization. -""" -function assemble_columns!( - rs::ReductionState{<:Any, S}, unreduced, reduced, progress -) where S - C = coface_type(S) - new_unreduced = C[] - new_reduced = C[] - - if progress - progbar = Progress( - length(unreduced) + length(reduced), desc="Assembling columns... " - ) - end - for cols in (unreduced, reduced) - for simplex in cols - for coface in coboundary(rs.filtration, simplex, Val(false)) - if !has_column(rs.reduction_matrix, coface) - push!(new_unreduced, abs(coface)) - else - push!(new_reduced, abs(coface)) - end - end - progress && next!(progbar) - end - end - sort!(new_unreduced, rev=true) - progress && printstyled("Assembled $(length(new_unreduced)) columns.\n", color=:green) - return new_unreduced, new_reduced -end - -""" - nth_intervals(filtration, columns, simplices, cutoff, field_type, ::Val{reps}; next=true) - -Compute the ``n``-th intervals of persistent cohomology. The ``n`` is determined from the -`eltype` of `columns`. If `assemble` is `true`, assemble columns for the next dimension. - -Only keep intervals with persistence larger than `cutoff`. Compute homology with -coefficients in `field_type`. If `reps` is `true`, compute representative cocycles. Show a -progress bar if `progress` is set. -""" -function nth_intervals( - filtration, unreduced, reduced, cutoff, field_type, reps, assemble, progress -) - rs = ReductionState{field_type, eltype(unreduced)}(filtration) - sizehint!(rs.reduction_matrix, length(unreduced)) - intervals = compute_intervals!(rs, unreduced, cutoff, reps, progress) - if assemble - return (intervals, assemble_columns!(rs, unreduced, reduced, progress)...) - else - return (intervals, nothing, nothing) - end -end - -""" - zeroth_representative(dset, vertex, reps, CE, V) - -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))) - end - else - return nothing - end -end - -""" - zeroth_intervals(filtration, cutoff, field_type, reps, progress) - -Compute 0-dimensional persistent homology using Kruskal's Algorithm. - -Only keep intervals with desired birth/death `cutoff`. Compute homology with coefficients in -`field_type`. If `reps` is `true`, compute representative cocycles. Show a progress bar if -`progress` is set. -""" -function zeroth_intervals(filtration, cutoff, field_type, reps, progress) - T = dist_type(filtration) - V = vertex_type(filtration) - CE = chain_element_type(V, field_type) - dset = DisjointSetsWithBirth([birth(filtration, v) for v in 1:n_vertices(filtration)]) - if reps - intervals = PersistenceInterval{Vector{CE}}[] - else - intervals = PersistenceInterval{Nothing}[] - end - reduced = edge_type(filtration)[] - unreduced = edge_type(filtration)[] - simplices = edges(filtration) - if progress - progbar = Progress(length(simplices), desc="Computing 0d intervals... ") - end - for sx in simplices - u, v = vertices(sx) - 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) - end - union!(dset, i, j) - push!(reduced, sx) - else - push!(unreduced, sx) - 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)) - end - end - reverse!(unreduced) - progress && printstyled("Assembled $(length(unreduced)) columns.\n", color=:green) - return ( - sort!(PersistenceDiagram(0, intervals, threshold(filtration))), unreduced, reduced, - ) -end - -""" - ripserer(dists::AbstractMatrix; kwargs...) - ripserer(points; metric=Euclidean(), births, kwargs...) - ripserer(filtration::AbstractFiltration; kwargs...) - -Compute the persistent homology of metric space represented by `dists`, `points` and -`metric` or an `::AbstractFiltration`. - -If using points, `points` must be an array of bitstypes, such as `NTuple`s or `SVectors`. - -# Keyoword Arguments - -* `dim_max`: compute persistent homology up to this dimension. Defaults to `1`. -* `modulus`: compute persistent homology with coefficients in the prime field of integers - mod `modulus`. Defaults to `2`. -* `field_type`: use this type of field of coefficients. Defaults to `Mod{modulus}`. -* `threshold`: compute persistent homology up to diameter smaller than threshold. - For non-sparse Rips filtrations, it defaults to radius of input space. -* `cutoff`: only keep intervals with `persistence(interval) > cutoff`. Defaults to `0`. -* `representatives`: if `true`, return representative cocycles along with persistence - intervals. Defaults to `false`. -* `progress`: If `true`, show a progress bar. Defaults to `false`. -* `metric`: when calculating persistent homology from points, any metric from - [`Distances.jl`](https://github.com/JuliaStats/Distances.jl) can be used. Defaults to - `Euclidean()`. -* `births`: when calculating persistent homology from points, births can be used to add - birth times to vertices. Defaults to all births equal to `0`. -""" -function ripserer( - dists::AbstractMatrix; - dim_max=1, - sparse=false, - cutoff=0, - representatives=false, - modulus=2, - field_type=Mod{modulus}, - progress=false, - kwargs..., # kwargs for filtration -) - if issparse(dists) - filtration = SparseRips(dists; kwargs...) - else - filtration = Rips(dists; kwargs...) - end - return ripserer( - filtration; - dim_max=dim_max, - representatives=representatives, - cutoff=cutoff, - field_type=field_type, - progress=progress, - ) -end - -function ripserer(points; metric=Euclidean(), births=nothing, kwargs...) - dists = distances(metric, points, births) - return ripserer(dists; kwargs...) -end - -function ripserer( - filtration::AbstractFiltration; - dim_max=1, representatives=false, cutoff=0, field_type=Mod{2}, progress=false -) - return ripserer(filtration, cutoff, field_type, Val(dim_max), representatives, progress) -end - -function ripserer( - filtration, cutoff, field_type, ::Val{dim_max}, reps, progress -) where dim_max - res = PersistenceDiagram[] - res_0, cols, sxs = zeroth_intervals(filtration, cutoff, field_type, reps, progress) - push!(res, res_0) - for dim in 1:dim_max - res_n, cols, sxs = nth_intervals( - filtration, cols, sxs, cutoff, field_type, reps, dim ≠ dim_max, progress) - push!(res, res_n) - # TODO: A lot of things get allocated and then destroyed in nth_intervals, so it - # might make sense to do a GC.gc() here. - end - return res -end diff --git a/src/reduction_structs.jl b/src/reduction_structs.jl deleted file mode 100644 index 712f551d..00000000 --- a/src/reduction_structs.jl +++ /dev/null @@ -1,250 +0,0 @@ -""" - ReductionMatrix{C, SE} - -A representation of the reduction matrix. It is indexed by simplices of type `C` (or chain -element types) and holds chain elements of one dimension lower. Supports the following -operations. - -* `insert_column!(::ReductionMatrix, i::K)`: add a new column with column index `i`. -* `has_column(::ReductionMatrix, i::K)`: return `true` if the matrix has any entries in the - `i`-th column. -* `push!(::ReductionMatrix, val::V)`: push `val` to *the last* column that was added to the - matrix. -""" -struct ReductionMatrix{C<:AbstractSimplex, SE<:AbstractChainElement} - column_index ::Dict{C, Int} - colptr ::Vector{Int} - nzval ::Vector{SE} - - function ReductionMatrix{C, SE}() where {C, SE} - return new{C, SE}(Dict{C, Int}(), Int[1], SE[]) - end -end - -function has_column(rm::ReductionMatrix{C}, sx::C) where C - return haskey(rm.column_index, abs(sx)) -end -function has_column(rm::ReductionMatrix{C}, ce::AbstractChainElement{C}) where C - return haskey(rm.column_index, simplex(ce)) -end - -function insert_column!(rm::ReductionMatrix{C}, ce::AbstractChainElement{C}) where C - return insert_column!(rm, simplex(ce)) -end -function insert_column!(rm::ReductionMatrix{C}, sx::C) where C - rm.column_index[abs(sx)] = length(rm.colptr) - push!(rm.colptr, rm.colptr[end]) - return rm -end - -function Base.push!(rm::ReductionMatrix, value) - push!(rm.nzval, value) - rm.colptr[end] += 1 - return value -end - -""" - append_unique_times!(rm::ReductionMatrix, values, factor) - -Append `values`, sorted with duplicates added together and multiplied by `factor` to `rm`. -""" -function append_unique_times!(rm::ReductionMatrix, values, factor) - sort!(values, alg=QuickSort) - prev = values[1] - for i in 2:length(values) - @inbounds current = values[i] - if current == prev - prev += current - else - !iszero(prev) && push!(rm, prev * factor) - prev = current - end - end - !iszero(prev) && push!(rm, prev * factor) -end - -function Base.sizehint!(rm::ReductionMatrix, n) - sizehint!(rm.column_index, n) - sizehint!(rm.colptr, n) - sizehint!(rm.nzval, n) -end - -Base.eltype(::Type{<:ReductionMatrix{<:Any, SE}}) where SE = SE -Base.length(rm::ReductionMatrix) = length(rm.colptr) - 1 -Base.lastindex(rm::ReductionMatrix) = length(rm.colptr) - 1 -function Base.getindex(rm::ReductionMatrix{C}, ce::AbstractChainElement{C}) where C - return RMColumnIterator{typeof(rm)}(rm, rm.column_index[simplex(ce)]) -end -function Base.getindex(rm::ReductionMatrix{C}, sx::C) where C - return RMColumnIterator{typeof(rm)}(rm, rm.column_index[sx]) -end - -""" - RMColumnIterator{R} - -An iterator over a column of a `ReductionMatrix`. -""" -struct RMColumnIterator{R<:ReductionMatrix} - rm ::R - idx ::Int -end - -Base.IteratorSize(::Type{RMColumnIterator}) = Base.HasLength() -Base.IteratorEltype(::Type{RMColumnIterator}) = Base.HasEltype() -Base.eltype(::Type{<:RMColumnIterator{R}}) where R = eltype(R) -Base.length(ci::RMColumnIterator) = ci.rm.colptr[ci.idx + 1] - ci.rm.colptr[ci.idx] - -function Base.iterate(ci::RMColumnIterator, i=1) - colptr = ci.rm.colptr - index = i + colptr[ci.idx] - 1 - if index ≥ colptr[ci.idx + 1] - return nothing - else - return (ci.rm.nzval[index], i + 1) - end -end - - -# columns ================================================================================ # -""" - Column{CE<:AbstractChainElement} - -Wrapper around `BinaryMinHeap{CE}`. Acts like a heap of chain elements, where elements with -the same simplex are summed together and elements with coefficient value `0` are ignored. -Support `push!`ing chain elements or simplices. Unlike a regular heap, it returns nothing -when `pop!` is used on an empty column. -""" -struct Column{CE<:AbstractChainElement} - heap::Vector{CE} - - Column{CE}() where CE = new{CE}(CE[]) -end - -Base.empty!(col::Column) = empty!(col.heap) -Base.isempty(col::Column) = isempty(col.heap) - -function Base.sizehint!(col::Column, size) - sizehint!(col.heap, size) - return col -end - -""" - pop_pivot!(column::Column) - -Pop the pivot from `column`. If there are multiple simplices with the same index on the top -of the column, sum them together. If they sum to 0, pop the next column. Return -`nothing` if column is empty or if all pivots summed to 0. -""" -function pop_pivot!(column::Column) - isempty(column) && return nothing - heap = column.heap - - pivot = heappop!(heap) - while !isempty(heap) - if iszero(pivot) - pivot = heappop!(heap) - elseif first(heap) == pivot - pivot += heappop!(heap) - else - break - end - end - return iszero(pivot) ? nothing : pivot -end - -""" - pivot(column) - -Return the pivot of the column - the element with the lowest diameter. -""" -function pivot(column::Column) - pivot = pop_pivot!(column) - if !isnothing(pivot) - heappush!(column.heap, pivot) - end - return pivot -end - -Base.push!(column::Column{CE}, simplex) where CE = push!(column, CE(simplex)) -function Base.push!(column::Column{CE}, element::CE) where CE - heap = column.heap - @inbounds if !isempty(heap) && heap[1] == element - heap[1] += element - else - heappush!(heap, element) - end - return column -end - -nonheap_push!(column::Column{CE}, simplex) where CE = push!(column.heap, CE(simplex)) -repair!(column::Column) = heapify!(column.heap) -Base.first(column::Column) = isempty(column) ? nothing : first(column.heap) - -# disjointset with birth ================================================================= # -""" - DisjointSetsWithBirth{T} - -Almost identical to `DataStructures.IntDisjointSets`, but keeps track of vertex birth times. -Has no `num_groups` method. -""" -struct DisjointSetsWithBirth{T} - parents ::Vector{Int} - ranks ::Vector{Int} - births ::Vector{T} - - function DisjointSetsWithBirth(births::AbstractVector{T}) where T - n = length(births) - return new{T}(collect(1:n), fill(0, n), copy(births)) - end -end - -Base.length(s::DisjointSetsWithBirth) = - length(s.parents) - -function DataStructures.find_root!(s::DisjointSetsWithBirth, x) - parents = s.parents - p = parents[x] - @inbounds if parents[p] != p - parents[x] = p = find_root!(s, p) - end - return p -end - -""" - find_leaves!(s::DisjointSetsWithBirth, x) - -Find all leaves below `x`, i.e. vertices that have `x` as root. -""" -function find_leaves!(s::DisjointSetsWithBirth, x) - leaves = Int[] - for i in 1:length(s) - find_root!(s, i) == x && push!(leaves, i) - end - return leaves -end - -function Base.union!(s::DisjointSetsWithBirth, x, y) - parents = s.parents - xroot = find_root!(s, x) - yroot = find_root!(s, y) - xroot != yroot ? root_union!(s, xroot, yroot) : xroot -end - -function DataStructures.root_union!(s::DisjointSetsWithBirth, x, y) - parents = s.parents - rks = s.ranks - births = s.births - @inbounds xrank = rks[x] - @inbounds yrank = rks[y] - - if xrank < yrank - x, y = y, x - elseif xrank == yrank - rks[x] += 1 - end - @inbounds parents[y] = x - @inbounds births[x] = min(births[x], births[y]) - return x -end - -birth(dset::DisjointSetsWithBirth, i) = dset.births[i] diff --git a/src/reductionmatrix.jl b/src/reductionmatrix.jl new file mode 100644 index 00000000..c8a5c083 --- /dev/null +++ b/src/reductionmatrix.jl @@ -0,0 +1,414 @@ +""" + ReducedMatrix{S<:AbstractSimplex, E<:AbstractChainElement} + +Representation of the reduced part of the matrix. Is indexed by `S`, indexing iterates +values of type `E`. `is_homology` controls the ordering of values. + +Changes to the matrix are recorded in a buffer with `record!`. Once a column is ready to be +added, use `commit!` to move the changes from the buffer to the matrix and create a +column. Changes can also be discarded with `discard!`. +""" +struct ReducedMatrix{S<:AbstractSimplex, E<:AbstractChainElement, O<:Base.Ordering} + column_index::Dict{S, Int} + indices::Vector{Int} + values::Vector{E} + buffer::Vector{E} + ordering::O + + function ReducedMatrix{S, E}(ordering::O) where {S, E, O} + return new{S, E, O}(Dict{S, Int}(), Int[1], E[], E[], ordering) + end +end + +function Base.sizehint!(matrix::ReducedMatrix, size) + sizehint!(matrix.column_index, size) + sizehint!(matrix.indices, size + 1) + sizehint!(matrix.values, size) + return matrix +end + +""" + record!(matrix::ReducedMatrix, element) + +Record the operation that was performed in the buffer. +""" +function record!(matrix::ReducedMatrix, element::AbstractChainElement) + push!(matrix.buffer, element) + return element +end + +function record!(matrix::ReducedMatrix, elements, factor) + i = length(matrix.buffer) + resize!(matrix.buffer, length(matrix.buffer) + length(elements)) + @inbounds for element in elements + i += 1 + matrix.buffer[i] = element * factor + end + return elements +end + +""" + commit!(matrix::ReducedMatrix, simplex, factor) + +Commit the changes that were `record!`ed, creating a new column indexed by `simplex`. This +sorts the buffer and adds duplicates together. All entries are multiplied by `factor`. +""" +function commit!(matrix::ReducedMatrix, simplex, factor) + @assert sign(simplex) == 1 + isempty(matrix.buffer) && return matrix + + sort!(matrix.buffer, alg=QuickSort, order=matrix.ordering) + i = 0 + @inbounds prev = matrix.buffer[1] + @inbounds for j in 2:length(matrix.buffer) + current = matrix.buffer[j] + if current == prev + prev += current + else + if !iszero(prev) + i += 1 + matrix.buffer[i] = prev * factor + end + prev = current + end + end + if !iszero(prev) + i += 1 + @inbounds matrix.buffer[i] = prev * factor + end + if i > 0 + resize!(matrix.buffer, i) + append!(matrix.values, matrix.buffer) + push!(matrix.indices, matrix.indices[end] + i) + matrix.column_index[simplex] = length(matrix.indices) - 1 + end + + empty!(matrix.buffer) + return matrix +end + +""" + discard!(matrix::ReducedMatrix) + +Undo recorded changes. +""" +function discard!(matrix::ReducedMatrix) + empty!(matrix.buffer) + return matrix +end + +Base.eltype(::Type{<:ReducedMatrix{<:Any, E}}) where E = E +Base.length(matrix::ReducedMatrix) = length(matrix.indices) - 1 + +Base.haskey(matrix::ReducedMatrix, simplex) = haskey(matrix.column_index, abs(simplex)) + +function Base.getindex(matrix::ReducedMatrix{S}, element::AbstractChainElement{S}) where S + return matrix[simplex(element)] +end +function Base.getindex(matrix::ReducedMatrix{S}, simplex::S) where S + index = get(matrix.column_index, abs(simplex), 0) + if index == 0 + from = 0 + to = -1 + else + from = matrix.indices[index] + to = matrix.indices[index + 1] - 1 + end + + # Constructing directly is equivalent to Base.unsafe_view + return SubArray(matrix.values, (from:to,)) +end + +# ======================================================================================== # +struct WorkingBoundary{E<:AbstractChainElement, O<:Base.Ordering} + heap::Vector{E} + ordering::O + + function WorkingBoundary{E}(ordering::O) where {E, O} + new{E, O}(E[], ordering) + end +end + +Base.empty!(col::WorkingBoundary) = empty!(col.heap) +Base.isempty(col::WorkingBoundary) = isempty(col.heap) + +function Base.sizehint!(col::WorkingBoundary, size) + sizehint!(col.heap, size) + return col +end + +function _pop_pivot!(column::WorkingBoundary) + isempty(column) && return nothing + heap = column.heap + + pivot = heappop!(heap, column.ordering) + while !isempty(heap) + if iszero(pivot) + pivot = heappop!(heap, column.ordering) + elseif first(heap) == pivot + pivot += heappop!(heap, column.ordering) + else + break + end + end + return iszero(pivot) ? nothing : pivot +end + +""" + get_pivot!(column) + +Return the pivot of the column - the element with the lowest diameter. Duplicates are summed +together, zero elements are discarded. Return `nothing` if there is no pivot. +""" +function get_pivot!(column::WorkingBoundary) + pivot = _pop_pivot!(column) + if !isnothing(pivot) + heappush!(column.heap, pivot, column.ordering) + end + return pivot +end + +function Base.push!(column::WorkingBoundary{E}, simplex::AbstractSimplex) where E + push!(column, E(simplex)) +end +function Base.push!(column::WorkingBoundary{E}, element::E) where E + heap = column.heap + @inbounds if !isempty(heap) && heap[1] == element + heap[1] += element + if iszero(heap[1]) + heappop!(heap, column.ordering) + end + else + heappush!(heap, element, column.ordering) + end + return column +end + +function nonheap_push!(column::WorkingBoundary{E}, simplex::AbstractSimplex) where E + push!(column.heap, E(simplex)) +end +function nonheap_push!(column::WorkingBoundary{E}, element::E) where E + push!(column.heap, element) +end + +repair!(column::WorkingBoundary) = heapify!(column.heap, column.ordering) +Base.first(column::WorkingBoundary) = first(column.heap) + +function move!(column::WorkingBoundary{E}) where E + dst = E[] + while (pivot = _pop_pivot!(column)) ≠ nothing + push!(dst, pivot) + end + return dst +end + +# ======================================================================================== # +struct ReductionMatrix{ + Co, Field, Filtration, Simplex, SimplexElem, Face, FaceElem, O<:Base.Ordering +} + filtration::Filtration + reduced::ReducedMatrix{Face, SimplexElem, O} + working_boundary::WorkingBoundary{FaceElem, O} + columns_to_reduce::Vector{Simplex} + columns_to_skip::Vector{Simplex} +end + +function ReductionMatrix{Co, Field}( + filtration::Filtration, columns_to_reduce, columns_to_skip +) where {Co, Field, Filtration} + + Simplex = eltype(columns_to_reduce) + Face = Co ? coface_type(Simplex) : face_type(Simplex) + ordering = Co ? Base.Order.Forward : Base.Order.Reverse + O = typeof(ordering) + SimplexElem = chain_element_type(Simplex, Field) + FaceElem = chain_element_type(Face, Field) + + reduced = ReducedMatrix{Face, SimplexElem}(ordering) + sizehint!(reduced, length(columns_to_reduce)) + working_boundary = WorkingBoundary{FaceElem}(ordering) + + return ReductionMatrix{Co, Field, Filtration, Simplex, SimplexElem, Face, FaceElem, O}( + filtration, + reduced, + working_boundary, + columns_to_reduce, + columns_to_skip + ) +end + +""" + co_boundary(matrix, simplex) + +Iterate over the (co)boundary of the `simplex`. Chooses between the boundary and the +coboundary based on the matrix type. +""" +function co_boundary(matrix::ReductionMatrix{true}, simplex::AbstractSimplex) + return coboundary(matrix.filtration, simplex) +end +function co_boundary(matrix::ReductionMatrix{false}, simplex::AbstractSimplex) + return boundary(matrix.filtration, simplex) +end + +simplex_type(::ReductionMatrix{<:Any, <:Any, <:Any, S}) where S = S +simplex_element(::ReductionMatrix{<:Any, <:Any, <:Any, <:Any, E}) where E = E +face_element(::ReductionMatrix{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, E}) where E = E +dim(::ReductionMatrix{true, <:Any, <:Any, S}) where S = dim(S) +dim(::ReductionMatrix{false, <:Any, <:Any, S}) where S = dim(S) - 1 + +function initialize_boundary!(matrix::ReductionMatrix{Co}, column_simplex) where Co + # TODO: can emergent pairs be safely be enabled for all kinds of complexes? + # An argument that disables it could be added. + empty!(matrix.working_boundary) + for face in co_boundary(matrix, column_simplex) + # Checking this on every face helps if more than one face has the same diameter. + if Co && diam(face) == diam(column_simplex) && !haskey(matrix.reduced, face) + empty!(matrix.working_boundary) + return face_element(matrix)(face) + end + nonheap_push!(matrix.working_boundary, face) + end + if isempty(matrix.working_boundary) + return nothing + else + repair!(matrix.working_boundary) + return first(matrix.working_boundary) + end +end + +function add!(matrix::ReductionMatrix, column, factor) + record!(matrix.reduced, column, factor) + for element in column + for face in co_boundary(matrix, simplex(element)) + push!( + matrix.working_boundary, + face_element(matrix)(face, coefficient(element) * factor) + ) + end + end + return matrix +end + +function reduce_column!(matrix::ReductionMatrix, column_simplex, cutoff, ::Type{I}) where I + pivot = initialize_boundary!(matrix, column_simplex) + + while !isnothing(pivot) + column = matrix.reduced[pivot] + isempty(column) && break + + add!(matrix, column, -coefficient(pivot)) + pivot = get_pivot!(matrix.working_boundary) + end + if isnothing(pivot) + discard!(matrix.reduced) + else + record!(matrix.reduced, simplex_element(matrix)(column_simplex)) + commit!(matrix.reduced, simplex(pivot), inv(coefficient(pivot))) + end + + return interval(I, matrix, column_simplex, pivot, cutoff) +end + +function birth_death(matrix::ReductionMatrix{true}, simplex, pivot) + return diam(simplex), isnothing(pivot) ? Inf : diam(pivot) +end +function birth_death(matrix::ReductionMatrix{false}, simplex, pivot) + return isnothing(pivot) ? Inf : diam(pivot), diam(simplex) +end + +# Interval with no representative, (co)homology. +function interval( + ::Type{PersistenceInterval{Nothing}}, matrix, simplex, pivot, cutoff +) + birth, death = birth_death(matrix, simplex, 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) + + if !isfinite(death) + return PersistenceInterval(birth, death, V()) + elseif death - birth > cutoff + return PersistenceInterval(birth, death, 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) + + if death - birth > cutoff + return PersistenceInterval(birth, death, move!(matrix.working_boundary)) + else + return nothing + end +end + +function compute_intervals!( + matrix::ReductionMatrix{Co}, cutoff, progress, ::Val{reps} +) where {Co, reps} + if reps + representative_type = Co ? simplex_element(matrix) : face_element(matrix) + intervals = PersistenceInterval{Vector{representative_type}}[] + else + intervals = PersistenceInterval{Nothing}[] + end + if progress + progbar = Progress( + length(matrix.columns_to_reduce), desc="Computing $(dim(matrix))d intervals... " + ) + end + for column in matrix.columns_to_reduce + interval = reduce_column!(matrix, column, cutoff, eltype(intervals)) + if !isnothing(interval) + push!(intervals, interval) + end + progress && next!(progbar) + end + return sort!( + PersistenceDiagram(dim(matrix), intervals, threshold(matrix.filtration)) + ) +end + +simplex_name(::Type{<:Simplex{2}}) = "triangles" +simplex_name(::Type{<:Simplex{3}}) = "tetrahedra" +simplex_name(::Type{<:AbstractSimplex{D}}) where D = "$D-simplices" + +function next_matrix(matrix::ReductionMatrix{Co, Field}, progress) where {Co, Field} + C = coface_type(simplex_type(matrix)) + new_to_reduce = C[] + new_to_skip = C[] + Co && sizehint!(new_to_skip, length(matrix.reduced)) + + if progress + progbar = Progress( + length(matrix.columns_to_reduce) + length(matrix.columns_to_skip), + desc="Assembling... " + ) + end + for simplex in Iterators.flatten((matrix.columns_to_reduce, matrix.columns_to_skip)) + for coface in coboundary(matrix.filtration, simplex, Val(false)) + # Clearing optimization only enabled for cohomology. + if Co && haskey(matrix.reduced, coface) + push!(new_to_skip, abs(coface)) + else + push!(new_to_reduce, abs(coface)) + end + end + progress && next!(progbar) + end + progress && printstyled( + stderr, "Assembled $(length(new_to_reduce)) $(simplex_name(C)). Sorting... ", + color=:green + ) + sort!(new_to_reduce, rev=Co) + progress && printstyled(stderr, "done.\n", color=:green) + + return ReductionMatrix{Co, Field}(matrix.filtration, new_to_reduce, new_to_skip) +end diff --git a/src/ripserer.jl b/src/ripserer.jl new file mode 100644 index 00000000..7ae4f945 --- /dev/null +++ b/src/ripserer.jl @@ -0,0 +1,130 @@ +""" + ripserer(dists::AbstractMatrix; kwargs...) + ripserer(points; metric=Euclidean(), births, kwargs...) + ripserer(filtration::AbstractFiltration; kwargs...) + +Compute the persistent homology of metric space represented by `dists`, `points` and +`metric` or an `::AbstractFiltration`. + +If using points, `points` must be an array of bitstypes, such as `NTuple`s or `SVectors`. + +# Keyoword Arguments + +* `dim_max`: compute persistent homology up to this dimension. Defaults to `1`. +* `modulus`: compute persistent homology with coefficients in the prime field of integers + mod `modulus`. Defaults to `2`. +* `field_type`: use this type of field of coefficients. Defaults to `Mod{modulus}`. +* `threshold`: compute persistent homology up to diameter smaller than threshold. + For non-sparse Rips filtrations, it defaults to radius of input space. +* `cutoff`: only keep intervals with `persistence(interval) > cutoff`. Defaults to `0`. +* `reps`: if `true`, return representative cocycles along with persistence intervals. + Defaults to `false`. +* `progress`: If `true`, show a progress bar. Defaults to `false`. +* `metric`: when calculating persistent homology from points, any metric from + [`Distances.jl`](https://github.com/JuliaStats/Distances.jl) can be used. Defaults to + `Euclidean()`. +* `births`: when calculating persistent homology from points, births can be used to add + birth times to vertices. Defaults to all births equal to `0`. +""" +function ripserer( + dists::AbstractMatrix; + dim_max=1, + sparse=false, + cutoff=0, + reps=false, + modulus=2, + field_type=Mod{modulus}, + progress=false, + cohomology=true, + kwargs..., # kwargs for filtration +) + if issparse(dists) + filtration = SparseRips(dists; kwargs...) + else + filtration = Rips(dists; kwargs...) + end + return ripserer( + filtration; + dim_max=dim_max, + reps=reps, + cutoff=cutoff, + field_type=field_type, + progress=progress, + cohomology=cohomology, + ) +end + +function ripserer(points; metric=Euclidean(), births=nothing, kwargs...) + dists = distances(metric, points, births) + return ripserer(dists; kwargs...) +end + +function ripserer( + filtration::AbstractFiltration; + dim_max=1, reps=false, cutoff=0, field_type=Mod{2}, progress=false, cohomology=true +) + if cohomology + return _cohomology( + filtration, cutoff, progress, field_type, Val(dim_max), Val(reps) + ) + else + return _homology( + filtration, cutoff, progress, field_type, Val(dim_max), Val(reps) + ) + end +end + +function _cohomology( + filtration, cutoff, progress, ::Type{F}, ::Val{dim_max}, ::Val{reps} +) where {F, dim_max, reps} + result = PersistenceDiagram[] + zeroth, to_reduce, to_skip = zeroth_intervals( + filtration, cutoff, progress, F, Val(reps) + ) + push!(result, zeroth) + if dim_max == 0 + return result + else + matrix = ReductionMatrix{true, F}(filtration, to_reduce, to_skip) + for dim in 1:dim_max + push!(result, compute_intervals!(matrix, cutoff, progress, Val(reps))) + + if dim < dim_max + matrix = next_matrix(matrix, progress) + end + end + + return result + end +end + +# Homology is still experimental and does not always work correctly. It does not yet output +# infinite intervals. +function _homology( + filtration, cutoff, progress, ::Type{F}, ::Val{dim_max}, ::Val{reps} +) where {F, dim_max, reps} + result = PersistenceDiagram[] + zeroth, to_reduce, to_skip = zeroth_intervals( + filtration, cutoff, progress, F, Val(reps) + ) + push!(result, zeroth) + + # We want to start with triangles. Starting at the highest dimension and going down, as + # dictated by the twist algorithm might not be worth it. The highest dimension may be + # sparse and twist is not supposed to bring a huge improvement to regular persistent + # homology. + # Constructing a matrix and throwing it away has some overhead, but is nothing compared + # to the slowness of homology. + matrix = next_matrix( + ReductionMatrix{false, F}(filtration, to_reduce, to_skip), progress + ) + for dim in 1:dim_max + push!(result, compute_intervals!(matrix, cutoff, progress, Val(reps))) + + if dim < dim_max + matrix = next_matrix(matrix, progress) + end + end + + return result +end diff --git a/src/simplex.jl b/src/simplex.jl index 1eb8389b..82252f1a 100644 --- a/src/simplex.jl +++ b/src/simplex.jl @@ -168,27 +168,7 @@ where ``i_k`` are the simplex vertex indices. return expr end -# coboundaries =========================================================================== # -""" - IndexedCobounary{all_cofaces, D, F, S<:IndexedSimplex{D-1}} - -Iterator that evaluates the coboundary of a `(D - 1)`-dimensional indexed simplex. Uses the -filtration to determine which simplices are valid cofaces in the filtration and to determine -their diameter. If the type parameter `all_cofaces` is `true`, return all cofaces, otherwise -only return cofaces where the new vertex index is larger than all vertex indices in -`simplex`. `all_cofaces=false` is used with during the call to [`assemble_columns!`](@ref) -to generate the list of all simplices. - -# Fields - -* `filtration ::F` -* `simplex ::S` -* `vertices ::NTuple{D, Int}` - -# Constructor - - coboundary(filtration, simplex[, Val(false)]) -""" +# (co)boundaries ========================================================================= # struct IndexedCobounary{all_cofaces, D, F, S<:IndexedSimplex} filtration ::F simplex ::S @@ -231,6 +211,38 @@ function Base.iterate( end end +struct IndexedBoundary{D, F, S<:IndexedSimplex} + filtration::F + simplex::S + vertices::NTuple{D, Int} + + function IndexedBoundary(filtration::F, simplex::S) where {D, F, S<:IndexedSimplex{D}} + return new{D + 1, F, S}(filtration, simplex, vertices(simplex)) + end +end + +function boundary(filtration, simplex::IndexedSimplex) + return IndexedBoundary(filtration, simplex) +end + +function Base.iterate(bi::IndexedBoundary{D}, k=1) where D + diameter = missing + face_vertices = TupleTools.unsafe_tail(bi.vertices) + while ismissing(diameter) && k ≤ D + face_vertices = TupleTools.deleteat(bi.vertices, k) + diameter = diam(bi.filtration, face_vertices) + k += 1 + end + if !ismissing(diameter) + sign = ifelse(iseven(k), 1, -1) + new_index = index(face_vertices) * sign + + return face_type(bi.simplex)(new_index, diameter), k + else + return nothing + end +end + """ Simplex{D, T, I} <: IndexedSimplex{D, T, I} @@ -297,3 +309,4 @@ end index(sx::Simplex) = sx.index diam(sx::Simplex) = sx.diam coface_type(::Type{<:Simplex{D, T, I}}) where {D, T, I} = Simplex{D+1, T, I} +face_type(::Type{<:Simplex{D, T, I}}) where {D, T, I} = Simplex{D-1, T, I} diff --git a/src/zerodimensional.jl b/src/zerodimensional.jl new file mode 100644 index 00000000..fecc827b --- /dev/null +++ b/src/zerodimensional.jl @@ -0,0 +1,141 @@ +""" + DisjointSetsWithBirth{T} + +Almost identical to `DataStructures.IntDisjointSets`, but keeps track of vertex birth times. +Has no `num_groups` method. +""" +struct DisjointSetsWithBirth{T} + parents ::Vector{Int} + ranks ::Vector{Int} + births ::Vector{T} + + function DisjointSetsWithBirth(births::AbstractVector{T}) where T + n = length(births) + return new{T}(collect(1:n), fill(0, n), copy(births)) + end +end + +Base.length(s::DisjointSetsWithBirth) = length(s.parents) + +function DataStructures.find_root!(s::DisjointSetsWithBirth, x) + parents = s.parents + p = parents[x] + @inbounds if parents[p] != p + parents[x] = p = find_root!(s, p) + end + return p +end + +""" + find_leaves!(s::DisjointSetsWithBirth, x) + +Find all leaves below `x`, i.e. vertices that have `x` as root. +""" +function find_leaves!(s::DisjointSetsWithBirth, x) + leaves = Int[] + for i in 1:length(s) + find_root!(s, i) == x && push!(leaves, i) + end + return leaves +end + +function Base.union!(s::DisjointSetsWithBirth, x, y) + parents = s.parents + xroot = find_root!(s, x) + yroot = find_root!(s, y) + xroot != yroot ? root_union!(s, xroot, yroot) : xroot +end + +function DataStructures.root_union!(s::DisjointSetsWithBirth, x, y) + parents = s.parents + rks = s.ranks + births = s.births + @inbounds xrank = rks[x] + @inbounds yrank = rks[y] + + if xrank < yrank + x, y = y, x + elseif xrank == yrank + rks[x] += 1 + end + @inbounds parents[y] = x + @inbounds births[x] = min(births[x], births[y]) + return x +end + +birth(dset::DisjointSetsWithBirth, i) = dset.births[i] + +""" + zeroth_representative(dset, vertex, reps, CE, V) + +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))) + end + else + return nothing + end +end + +""" + zeroth_intervals(filtration, cutoff, progress, field_type, Val(reps)) + +Compute 0-dimensional persistent homology using Kruskal's Algorithm. + +Only keep intervals with desired birth/death `cutoff`. Compute homology with coefficients in +`field_type`. If `reps` is `true`, compute representative cocycles. Show a progress bar if +`progress` is set. +""" +function zeroth_intervals( + filtration, cutoff, progress, ::Type{F}, ::Val{reps} +) where {F, reps} + T = dist_type(filtration) + V = vertex_type(filtration) + CE = chain_element_type(V, F) + dset = DisjointSetsWithBirth([birth(filtration, v) for v in 1:n_vertices(filtration)]) + if reps + intervals = PersistenceInterval{Vector{CE}}[] + else + intervals = PersistenceInterval{Nothing}[] + end + to_skip = edge_type(filtration)[] + to_reduce = edge_type(filtration)[] + simplices = edges(filtration) + if progress + progbar = Progress(length(simplices), desc="Computing 0d intervals... ") + end + for sx in simplices + u, v = vertices(sx) + 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) + end + union!(dset, i, j) + push!(to_skip, sx) + else + push!(to_reduce, sx) + 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)) + end + end + reverse!(to_reduce) + progress && printstyled(stderr, "Assembled $(length(to_reduce)) edges.\n", color=:green) + return ( + sort!(PersistenceDiagram(0, intervals, threshold(filtration))), to_reduce, to_skip, + ) +end diff --git a/test/cubical.jl b/test/cubical.jl index a210e194..3c27092a 100644 --- a/test/cubical.jl +++ b/test/cubical.jl @@ -118,3 +118,42 @@ end @test length(cocob) == 1 end end + +@testset "boundary" begin + @testset "1d" begin + bnd = Cubelet{0, Float64, Int}[] + flt = Cubical(data3d) + cub = Cubelet{1}((3, 2), 1.0) + for f in boundary(flt, cub) + push!(bnd, f) + end + @test bnd == [Cubelet{0}((3,), 1.0), -Cubelet{0}((2,), 1.0)] + end + + @testset "2d" begin + bnd = Cubelet{1, Int, Int}[] + flt = Cubical(data2d) + cub = Cubelet{2}((10, 9, 3, 2), 2) + for f in boundary(flt, cub) + push!(bnd, f) + end + @test bnd == [Cubelet{1}((10, 3), 2), -Cubelet{1}((9, 2), 2), + Cubelet{1}((10, 9), 1), -Cubelet{1}((3, 2), 2)] + end + + @testset "3d" begin + bnd = Cubelet{2, Float64, Int}[] + flt = Cubical(data3d) + cub = Cubelet{3}((112, 111, 102, 101, 12, 11, 2, 1), 1.0) + for f in boundary(flt, cub) + push!(bnd, f) + end + @test bnd == [Cubelet{2}((112, 102, 12, 2), 1.0), + -Cubelet{2}((111, 101, 11, 1), 1.0), + Cubelet{2}((112, 111, 12, 11), 1.0), + -Cubelet{2}((102, 101, 2, 1), 1.0), + Cubelet{2}((112, 111, 102, 101), 1.0), + -Cubelet{2}((12, 11, 2, 1), 1.0), + ] + end +end diff --git a/test/reduction.jl b/test/reduction.jl deleted file mode 100644 index 64aafddc..00000000 --- a/test/reduction.jl +++ /dev/null @@ -1,266 +0,0 @@ -using Ripserer -using Ripserer: zeroth_intervals, ChainElement, PackedElement - -using Compat - -include("data.jl") - -@testset "ripserer" begin - @testset "full matrix, no threshold" begin - @testset "icosahedron" begin - res = ripserer(icosahedron, dim_max=2) - @test res[1] == [fill(PersistenceInterval(0.0, 1.0), 11); - PersistenceInterval(0.0, Inf)] - @test isempty(res[2]) - @test res[3] == [PersistenceInterval(1.0, 2.0)] - end - @testset "torus 16" begin - d0, d1, d2 = ripserer(torus(16), dim_max=2) - - @test length(d0) == 16 - - @test all(x -> birth(x) ≈ 0.5, d1) - @test count(x -> death(x) ≈ 1, d1) == 2 - @test count(x -> isapprox(death(x), 0.71, atol=0.1), d1) == 15 - - @test death(only(d2)) == 1 - end - @testset "torus 100" begin - d0, d1 = ripserer(torus(100), dim_max=1) - - @test length(d0) == 100 - - deaths = sort(death.(d1)) - @test deaths[end] ≈ 0.8 - @test deaths[end-1] ≈ 0.8 - @test deaths[end-2] < 0.5 - end - @testset "cycle" begin - d0, d1, d2, d3, d4 = ripserer(cycle, dim_max=4) - @test d0 == [fill(PersistenceInterval(0, 1), size(cycle, 1) - 1); - PersistenceInterval(0, Inf)] - @test d1 == [PersistenceInterval(1, 6)] - @test d2 == fill(PersistenceInterval(6, 7), 5) - @test d3 == [PersistenceInterval(7, 8)] - @test d4 == [] - - d0_7, d1_7, d2_7, d3_7, d4_7 = ripserer(cycle, dim_max=4, modulus=7) - @test all(d0 .== d0_7) - @test all(d1 .== d1_7) - @test all(d2 .== d2_7) - @test all(d3 .== d3_7) - @test all(d4 .== d4_7) - - d0r, d1r, d2r, d3r, d4r = ripserer(cycle, dim_max=4, field_type=Rational{Int}) - @test all(d0 .== d0r) - @test all(d1 .== d1r) - @test all(d2 .== d2r) - @test all(d3 .== d3r) - @test all(d4 .== d4r) - end - @testset "projective plane (modulus)" begin - _, d1_2, d2_2 = ripserer(projective_plane, dim_max=2) - _, d1_3, d2_3 = ripserer(projective_plane, dim_max=2, modulus=3) - @test d1_2 == [PersistenceInterval(1, 2)] - @test d2_2 == [PersistenceInterval(1, 2)] - @test isempty(d1_3) - @test isempty(d2_3) - end - end - - @testset "full matrix, with threshold" begin - @testset "icosahedron, high threshold" begin - res = ripserer(icosahedron, threshold=2, dim_max=2) - @test res[1] == [fill(PersistenceInterval(0.0, 1.0), 11); - PersistenceInterval(0.0, Inf)] - @test isempty(res[2]) - @test res[3] == [PersistenceInterval(1.0, 2.0)] - end - @testset "icosahedron, med threshold" begin - res = ripserer(icosahedron, dim_max=2, threshold=1) - @test res[1] == [fill(PersistenceInterval(0.0, 1.0), 11); - PersistenceInterval(0.0, Inf)] - @test isempty(res[2]) - @test res[3] == [PersistenceInterval(1.0, Inf)] - end - @testset "icosahedron, low threshold" begin - res = ripserer(icosahedron, dim_max=2, threshold=0.5) - @test res[1] == fill(PersistenceInterval(0.0, Inf), 12) - @test isempty(res[2]) - @test isempty(res[3]) - end - @testset "torus 16, high threshold" begin - d0, d1, d2 = ripserer(torus(16), dim_max=2, threshold=2) - - @test length(d0) == 16 - - @test all(x -> birth(x) ≈ 0.5, d1) - @test count(x -> death(x) ≈ 1, d1) == 2 - @test count(x -> isapprox(death(x), 0.71, atol=0.1), d1) == 15 - - @test death(only(d2)) == 1 - end - @testset "torus 16, med threshold" begin - d0, d1, d2 = ripserer(torus(16), dim_max=2, threshold=0.9) - - @test length(d0) == 16 - - @test all(x -> birth(x) ≈ 0.5, d1) - @test count(x -> death(x) == Inf, d1) == 2 - @test count(x -> isapprox(death(x), 0.71, atol=0.1), d1) == 15 - - @test last(only(d2)) == Inf - end - @testset "torus 16, low threshold" begin - d0, d1, d2 = ripserer(torus(16), dim_max=2, threshold=0.5) - - @test length(d0) == 16 - - @test all(x -> birth(x) ≈ 0.5, d1) - @test all(x -> death(x) == Inf, d1) - - @test isempty(d2) - end - @testset "projective plane (modulus), med threshold" begin - _, d1_2, d2_2 = ripserer(projective_plane, - dim_max=2, threshold=1) - _, d1_3, d2_3 = ripserer(projective_plane, - dim_max=2, modulus=3, threshold=1) - @test d1_2 == [PersistenceInterval(1, Inf)] - @test d2_2 == [PersistenceInterval(1, Inf)] - @test isempty(d1_3) - @test isempty(d2_3) - end - end - - @testset "sparse matrix" begin - @testset "icosahedron" begin - flt = SparseRips(icosahedron, threshold=2) - res = ripserer(flt, dim_max=2) - @test res[1] == [fill(PersistenceInterval(0.0, 1.0), 11); - PersistenceInterval(0.0, Inf)] - @test isempty(res[2]) - @test res[3] == [PersistenceInterval(1.0, 2.0)] - end - @testset "torus 16" begin - dists = sparse(torus(16)) - SparseArrays.fkeep!(dists, (_, _, v) -> v ≤ 1) - - d0, d1, d2 = ripserer(dists, dim_max=2) - - @test length(d0) == 16 - - @test all(x -> birth(x) ≈ 0.5, d1) - @test count(x -> death(x) ≈ 1, d1) == 2 - @test count(x -> isapprox(death(x), 0.71, atol=0.1), d1) == 15 - - @test last(only(d2)) == 1 - end - @testset "projective plane (modulus), med threshold" begin - dists = sparse(projective_plane) - SparseArrays.fkeep!(dists, (_, _, v) -> v ≤ 2) - - _, d1_2, d2_2 = ripserer(dists, dim_max=2, threshold=1) - _, d1_3, d2_3 = ripserer(dists, dim_max=2, modulus=3, threshold=1) - @test d1_2 == [PersistenceInterval(1, Inf)] - @test d2_2 == [PersistenceInterval(1, Inf)] - @test isempty(d1_3) - @test isempty(d2_3) - end - end - - @testset "representatives" begin - _, d1, d2 = ripserer(projective_plane, dim_max=2, representatives=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 "representatives - types" begin - d0, d1, d2, d3 = ripserer(cycle, dim_max=3, representatives=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, representatives=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 "representatives - thresh" begin - _, d1 = ripserer( - cycle, dim_max=1, representatives=true, threshold=1, field_type=Rational{Int} - ) - @test representative(only(d1)) == ChainElement{ - Simplex{0, Int, Int}, Rational{Int} - }[] - end - - @testset "lower star w/Rips" begin - data = [range(0, 1, length=5); - range(1, 0.5, length=5)[2:end]; - range(0.5, 2, length=4)[2:end]; - range(2, -1, length=4)[2:end]] - - # Create distance matrix from data, where neighboring points are connected by edges - # and the edge weights are equal to the max of both vertex births. - n = length(data) - dists = spzeros(n, n) - for i in 1:n - dists[i, i] = data[i] - end - for i in 1:n-1 - j = i + 1 - dists[i, j] = dists[j, i] = max(dists[i, i], dists[j, j]) - end - # 0-dimensional persistence should find values of minima and maxima of our data. - res = first(ripserer(dists, dim_max=0)) - mins = birth.(res) - maxs = death.(filter(isfinite, res)) - @test sort(mins) == [-1.0, 0.0, 0.0, 0.5] - @test sort(maxs) == [1.0, 2.0] - end - - @testset "image lower star" 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, d3, d4 = ripserer(Cubical(data), representatives=true, dim_max=4) - - @test d0 == [(0, Inf), (1, 2)] - @test d1 == [(0, 2)] - @test isempty(d2) - @test isempty(d3) - @test isempty(d4) - - @test vertices.(representative(d0[1])) == [(i,) for i in 1:length(data)] - @test vertices(only(representative(d0[2]))) == (13,) - end -end diff --git a/test/reduction_structs.jl b/test/reduction_structs.jl deleted file mode 100644 index 2f31cfab..00000000 --- a/test/reduction_structs.jl +++ /dev/null @@ -1,128 +0,0 @@ -using Ripserer -using Ripserer: ChainElement, coefficient, - ReductionMatrix, insert_column!, has_column, append_unique_times!, - Column, pop_pivot!, pivot, - DisjointSetsWithBirth - -using Compat - -@testset "ReductionMatrix" begin - @testset "push!, insert_column!, getindex, iterate" begin - for T in (Mod{3}, Rational{Int}) - CE = ChainElement{Simplex{3, Int, Int}, T} - SE = ChainElement{Simplex{2, Int, Int}, T} - rm = @inferred ReductionMatrix{Simplex{3, Int, Int}, SE}() - @test length(rm) == 0 - @test eltype(rm) ≡ SE - - insert_column!(rm, -Simplex{3}(3, 1)) - push!(rm, SE(Simplex{2}(1, 1), 1)) - push!(rm, SE(Simplex{2}(2, 1), 1)) - push!(rm, SE(Simplex{2}(3, 1), 1)) - push!(rm, SE(Simplex{2}(4, 1), 1)) - - insert_column!(rm, CE(Simplex{3}(10, 1), 1)) - push!(rm, SE(Simplex{2}(5, 2), 2)) - push!(rm, SE(Simplex{2}(6, 2), 2)) - push!(rm, SE(Simplex{2}(7, 2), 2)) - - insert_column!(rm, Simplex{3}(1, 1)) - - insert_column!(rm, Simplex{3}(15, 1)) - push!(rm, SE(Simplex{2}(8, 3), -1)) - - @test has_column(rm, Simplex{3}(3, 1)) - @test collect(rm[Simplex{3}(3, 1)]) == SE.(Simplex{2}.([1, 2, 3, 4], 1), 1) - - @test has_column(rm, Simplex{3}(10, 1)) - @test length(rm[Simplex{3}(10, 1)]) == 3 - @test all(ce -> coefficient(ce) == 2, rm[Simplex{3}(10, 1)]) - - @test has_column(rm, -Simplex{3}(1, 1)) - @test length(rm[Simplex{3}(1, 1)]) == 0 - @test eltype(collect(rm[Simplex{3}(1, 1)])) ≡ SE - - @test has_column(rm, CE(Simplex{3}(15, 1), 0)) - @test length(rm[CE(Simplex{3}(15, 1))]) == 1 - @test first(rm[CE(Simplex{3}(15, 1))]) == SE(Simplex{2}(8, 3), -1) - - @test !has_column(rm, Simplex{3}(2, 1)) - @test !has_column(rm, CE(Simplex{3}(100, 1))) - end - end - @testset "append_unique_times!" begin - CE = ChainElement{Simplex{3, Int, Int}, Mod{3}} - SE = ChainElement{Simplex{2, Int, Int}, Mod{3}} - rm = @inferred ReductionMatrix{Simplex{3, Int, Int}, SE}() - - insert_column!(rm, Simplex{3}(3, 1)) - append_unique_times!( - rm, SE.(Simplex{2}.([1, 4, 6, 4, 3, 1, 1, 4, 6, 5], 1)), Mod{3}(1) - ) - - result = collect(rm[Simplex{3}(3, 1)]) - @test issorted(result) - @test allunique(result) - @test length(result) == 3 - end -end - -@testset "Column" begin - @testset "single element" begin - CE = ChainElement{Simplex{1, Float64, Int}, Mod{2}} - col = Column{CE}() - push!(col, Simplex{1}(1, 2.0)) - push!(col, Simplex{1}(1, 2.0)) - push!(col, Simplex{1}(1, 2.0)) - push!(col, Simplex{1}(1, 2.0)) - push!(col, Simplex{1}(1, 2.0)) - - @test pop_pivot!(col) == CE(Simplex{1}(1, 2.0), 1) - @test isnothing(pivot(col)) - @test isempty(col) - - CE = ChainElement{Simplex{2, Float64, Int}, Mod{3}} - col = Column{CE}() - push!(col, Simplex{2}(1, 2.0)) - push!(col, Simplex{2}(1, 2.0)) - push!(col, Simplex{2}(1, 2.0)) - - @test isnothing(pivot(col)) - @test isnothing(pop_pivot!(col)) - @test isempty(col) - - CE = ChainElement{Simplex{2, Int, Int}, Rational{Int}} - col = Column{CE}() - push!(col, Simplex{2}(1, 2)) - push!(col, Simplex{2}(2, 3)) - push!(col, Simplex{2}(-1, 2)) - push!(col, Simplex{2}(-2, 3)) - push!(col, -Simplex{2}(1, 2)) - - @test coefficient(pivot(col)) == -1 - @test pop_pivot!(col) == CE(Simplex{2}(1, 2), -1) - @test isnothing(pivot(col)) - @test isnothing(pop_pivot!(col)) - @test isempty(col) - end - @testset "multiple" begin - CE = ChainElement{Simplex{3, Float64, Int}, Mod{5}} - col = Column{CE}() - push!(col, CE(Simplex{3}(2, 1.0), 3)) - push!(col, CE(Simplex{3}(3, 2.0), 4)) - push!(col, CE(Simplex{3}(2, 1.0), 2)) - push!(col, CE(Simplex{3}(1, 3.0), 2)) - push!(col, CE(Simplex{3}(3, 2.0), 1)) - push!(col, CE(Simplex{3}(4, 4.0), 4)) - push!(col, CE(Simplex{3}(4, 4.0), 4)) - push!(col, CE(Simplex{3}(4, 4.0), 4)) - push!(col, CE(Simplex{3}(4, 5.0), 4)) - push!(col, CE(Simplex{3}(4, 5.0), 1)) - - @test pop_pivot!(col) == CE(Simplex{3}(1, 3.0), 2) - @test pivot(col) == CE(Simplex{3}(4, 4.0), 2) - @test pop_pivot!(col) == CE(Simplex{3}(4, 4.0), 2) - @test isnothing(pop_pivot!(col)) - @test isnothing(pop_pivot!(col)) - end -end diff --git a/test/reductionmatrix.jl b/test/reductionmatrix.jl new file mode 100644 index 00000000..c935354f --- /dev/null +++ b/test/reductionmatrix.jl @@ -0,0 +1,275 @@ +using Ripserer + +using Random + +using Ripserer: chain_element_type, coefficient + +using Ripserer: ReducedMatrix, record!, commit!, discard! +using Ripserer: WorkingBoundary, nonheap_push!, get_pivot!, repair! +using Ripserer: ReductionMatrix, simplex_type, simplex_element, face_element, next_matrix + +@testset "ReducedMatrix" begin + for S in (Simplex{2, Int, Int}, Cubelet{1, Int, Int}), T in (Mod{3}, Rational{Int}) + C = coface_type(S) + CE = chain_element_type(C, T) + SE = chain_element_type(S, T) + + columns = C.([3, 10, 6, 8], 1) + colelems = CE.(columns) + + fwd = Base.Order.Forward + rev = Base.Order.Reverse + + @testset "ReducedMatrix with simplex type $S and field type $T" begin + @testset "a fresh ReducedMatrix is empty even if you commit nothing." begin + matrix = ReducedMatrix{C, SE}(fwd) + commit!(matrix, columns[1], T(2)) + + @test length(matrix) == 0 + + for col in Iterators.flatten((columns, colelems)) + @test isempty(matrix[col]) + end + end + + @testset "is still empty after recording values" begin + matrix = ReducedMatrix{C, SE}(rev) + vals = [ + SE(S(1, 1)), + SE(S(2, 1)), + SE(S(3, 1)), + SE(S(4, 1)), + ] + for v in vals + record!(matrix, v) + end + + @test length(matrix) == 0 + + for col in Iterators.flatten((columns, colelems)) + @test isempty(matrix[col]) + end + end + + @testset "commiting adds a column and changes length" begin + matrix = ReducedMatrix{C, SE}(fwd) + vals = [ + SE(S(1, 1)), + SE(S(2, 1)), + SE(S(3, 1)), + SE(S(4, 1)), + SE(S(4, 1)), + SE(S(-1, 1)), + ] + for v in vals + record!(matrix, v) + end + + commit!(matrix, columns[1], T(2)) + commit!(matrix, columns[2], T(2)) + + @test length(matrix) == 1 + @test index.(simplex.(collect(matrix[columns[1]]))) == [4, 3, 2] + @test index.(simplex.(collect(matrix[colelems[1]]))) == [4, 3, 2] + + @test coefficient.(collect(matrix[columns[1]])) == [T(4), T(2), T(2)] + @test coefficient.(collect(matrix[colelems[1]])) == [T(4), T(2), T(2)] + + for col in Iterators.flatten((columns[2:end], colelems[2:end])) + @test isempty(matrix[col]) + end + end + + @testset "discard! undos changes" begin + matrix = ReducedMatrix{C, SE}(rev) + vals = [ + SE(S(1, 1)), + SE(S(2, 1)), + SE(S(3, 1)), + SE(S(4, 1)), + SE(S(4, 1)), + SE(S(-1, 1)), + ] + for v in vals + record!(matrix, v) + end + discard!(matrix) + commit!(matrix, columns[2], T(2)) + + for v in vals + record!(matrix, v) + end + commit!(matrix, columns[1], T(2)) + + @test length(matrix) == 1 + @test index.(simplex.(collect(matrix[columns[1]]))) == [2, 3, 4] + @test index.(simplex.(collect(matrix[colelems[1]]))) == [2, 3, 4] + + @test coefficient.(collect(matrix[columns[1]])) == [T(2), T(2), T(4)] + @test coefficient.(collect(matrix[colelems[1]])) == [T(2), T(2), T(4)] + + for col in Iterators.flatten((columns[2:end], colelems[2:end])) + @test isempty(matrix[col]) + end + end + + @testset "committing values that sum to 0 does not create a column" begin + matrix = ReducedMatrix{C, SE}(fwd) + vals = [ + SE(S(1, 1)), + SE(S(-2, 1)), + SE(S(4, 1)), + SE(S(-4, 1)), + SE(S(2, 1)), + SE(S(-1, 1)), + ] + for v in vals + record!(matrix, v) + end + commit!(matrix, columns[1], T(2)) + + @test length(matrix) == 0 + for col in Iterators.flatten((columns, colelems)) + @test isempty(matrix[col]) + end + end + + @testset "committing multiple times creates multiple columns" begin + matrix = ReducedMatrix{C, SE}(rev) + vals_1 = [ + SE(S(1, 1)), + SE(S(-4, 1)), + SE(S(2, 1)), + SE(S(3, 1)), + SE(S(4, 1)), + ] + for v in vals_1 + record!(matrix, v) + end + commit!(matrix, columns[1], T(2)) + + vals_1 = [ + SE(S(4, 1)), + SE(S(1, 1)), + SE(S(5, 1)), + SE(S(6, 1)), + SE(S(-1, 1)), + ] + for v in vals_1 + record!(matrix, v) + end + commit!(matrix, columns[2], T(-1)) + + @test length(matrix) == 2 + + @test index.(simplex.(collect(matrix[columns[1]]))) == [1, 2, 3] + @test index.(simplex.(collect(matrix[columns[2]]))) == [4, 5, 6] + + for col in Iterators.flatten((columns[3:end], colelems[3:end])) + @test isempty(matrix[col]) + end + end + end + end +end + +@testset "WorkingBoundary" begin + for S in (Simplex{2, Int, Int}, Cubelet{1, Int, Int}), T in (Mod{3}, Rational{Int}) + SE = chain_element_type(S, T) + + elements = SE.(S.([1, -7, 2, 3, 4, 7, 5, 6, -1], [1, 7, 1, 1, 4, 7, 5, 6, 1])) + + fwd = Base.Order.Forward + rev = Base.Order.Reverse + + @testset "a fresh WorkingBoundary is empty and its pivot is nothing" begin + working_boundary = WorkingBoundary{SE}(fwd) + + @test isempty(working_boundary) + @test get_pivot!(working_boundary) ≡ nothing + @test_throws BoundsError first(working_boundary) + end + + @testset "pushing elements and getting pivot finds the lowest simplex (cohomology)" begin + working_boundary = WorkingBoundary{SE}(fwd) + for e in elements + push!(working_boundary, e) + end + + @test get_pivot!(working_boundary) == SE(S(3, 1)) + end + + @testset "the same happens for nonheap_push! with repair! (homology)" begin + working_boundary = WorkingBoundary{SE}(rev) + for e in elements + nonheap_push!(working_boundary, e) + end + repair!(working_boundary) + + @test get_pivot!(working_boundary) == SE(S(6, 6)) + end + + @testset "getting pivot and adding its inverse to the boundary repeatedly" begin + working_boundary = WorkingBoundary{SE}(fwd) + for e in elements[1:3] + nonheap_push!(working_boundary, e) + end + repair!(working_boundary) + for e in elements[4:end] + push!(working_boundary, e) + end + + result = [] + while (p = get_pivot!(working_boundary)) ≢ nothing + push!(working_boundary, -p) + push!(result, p) + end + @test length(result) == 5 + @test allunique(result) + @test issorted(result) + end + end +end + +@testset "ReducedMatrix" begin + for Co in (true, false), + S in (Simplex{2, Int, Int}, Cubelet{1, Int, Int}), + T in (Mod{3}, Rational{Int}) + + columns_1 = (S.([1, 2, 3, 4, 5], [1, 2, 3, 4, 5])) + columns_2 = (S.([6, 7, 8, 9, 10], [6, 7, 8, 9, 10])) + data = [0 1 1 1 1; 1 0 1 1 1; 1 1 0 1 1; 1 1 1 0 1; 1 1 1 1 0] + flt = S <: Simplex ? Rips(data) : Cubical(data) + co_str = Co ? "co" : "" + + SE = chain_element_type(S, T) + F = Co ? coface_type(S) : face_type(S) + FE = chain_element_type(F, T) + + @testset "$(co_str)homology, $S, $T" begin + @testset "construction inferrence, parameters" begin + @test begin + @inferred ReductionMatrix{Co, T}(flt, columns_1, columns_2) + true + end + + matrix = ReductionMatrix{Co, T}(flt, columns_1, columns_2) + + @test simplex_type(matrix) == S + @test simplex_element(matrix) == SE + @test face_element(matrix) == FE + @test dim(matrix) == (Co ? dim(S) : dim(S) - 1) + end + + @testset "next matrix inferrence" begin + matrix = ReductionMatrix{Co, T}(flt, columns_1, columns_2) + @test begin + @inferred next_matrix(matrix, false) + true + end + next = next_matrix(matrix, false) + @test dim(next) == dim(matrix) + 1 + end + end + end +end diff --git a/test/ripserer.jl b/test/ripserer.jl new file mode 100644 index 00000000..c7f4288b --- /dev/null +++ b/test/ripserer.jl @@ -0,0 +1,289 @@ +using Ripserer +using Ripserer: zeroth_intervals, ChainElement, PackedElement + +using Compat +using Suppressor + +include("data.jl") + +@testset "full matrix, no threshold" begin + @testset "icosahedron" begin + res = ripserer(icosahedron, dim_max=2) + @test res[1] == [fill(PersistenceInterval(0.0, 1.0), 11); + PersistenceInterval(0.0, Inf)] + @test isempty(res[2]) + @test res[3] == [PersistenceInterval(1.0, 2.0)] + end + @testset "torus 16" begin + d0, d1, d2 = ripserer(torus(16), dim_max=2) + + @test length(d0) == 16 + + @test all(x -> birth(x) ≈ 0.5, d1) + @test count(x -> death(x) ≈ 1, d1) == 2 + @test count(x -> isapprox(death(x), 0.71, atol=0.1), d1) == 15 + + @test death(only(d2)) == 1 + end + @testset "torus 100" begin + d0, d1 = ripserer(torus(100), dim_max=1) + + @test length(d0) == 100 + + deaths = sort(death.(d1)) + @test deaths[end] ≈ 0.8 + @test deaths[end-1] ≈ 0.8 + @test deaths[end-2] < 0.5 + end + @testset "cycle" begin + d0, d1, d2, d3, d4 = ripserer(cycle, dim_max=4) + @test d0 == [fill(PersistenceInterval(0, 1), size(cycle, 1) - 1); + PersistenceInterval(0, Inf)] + @test d1 == [PersistenceInterval(1, 6)] + @test d2 == fill(PersistenceInterval(6, 7), 5) + @test d3 == [PersistenceInterval(7, 8)] + @test d4 == [] + + d0_7, d1_7, d2_7, d3_7, d4_7 = ripserer(cycle, dim_max=4, modulus=7) + @test all(d0 .== d0_7) + @test all(d1 .== d1_7) + @test all(d2 .== d2_7) + @test all(d3 .== d3_7) + @test all(d4 .== d4_7) + + d0r, d1r, d2r, d3r, d4r = ripserer(cycle, dim_max=4, field_type=Rational{Int}) + @test all(d0 .== d0r) + @test all(d1 .== d1r) + @test all(d2 .== d2r) + @test all(d3 .== d3r) + @test all(d4 .== d4r) + end + @testset "projective plane (modulus)" begin + _, d1_2, d2_2 = ripserer(projective_plane, dim_max=2) + _, d1_3, d2_3 = ripserer(projective_plane, dim_max=2, modulus=3) + @test d1_2 == [PersistenceInterval(1, 2)] + @test d2_2 == [PersistenceInterval(1, 2)] + @test isempty(d1_3) + @test isempty(d2_3) + end +end + +@testset "full matrix, with threshold" begin + @testset "icosahedron, high threshold" begin + res = ripserer(icosahedron, threshold=2, dim_max=2) + @test res[1] == [fill(PersistenceInterval(0.0, 1.0), 11); + PersistenceInterval(0.0, Inf)] + @test isempty(res[2]) + @test res[3] == [PersistenceInterval(1.0, 2.0)] + end + @testset "icosahedron, med threshold" begin + res = ripserer(icosahedron, dim_max=2, threshold=1) + @test res[1] == [fill(PersistenceInterval(0.0, 1.0), 11); + PersistenceInterval(0.0, Inf)] + @test isempty(res[2]) + @test res[3] == [PersistenceInterval(1.0, Inf)] + end + @testset "icosahedron, low threshold" begin + res = ripserer(icosahedron, dim_max=2, threshold=0.5) + @test res[1] == fill(PersistenceInterval(0.0, Inf), 12) + @test isempty(res[2]) + @test isempty(res[3]) + end + @testset "torus 16, high threshold" begin + d0, d1, d2 = ripserer(torus(16), dim_max=2, threshold=2) + + @test length(d0) == 16 + + @test all(x -> birth(x) ≈ 0.5, d1) + @test count(x -> death(x) ≈ 1, d1) == 2 + @test count(x -> isapprox(death(x), 0.71, atol=0.1), d1) == 15 + + @test death(only(d2)) == 1 + end + @testset "torus 16, med threshold" begin + d0, d1, d2 = ripserer(torus(16), dim_max=2, threshold=0.9) + + @test length(d0) == 16 + + @test all(x -> birth(x) ≈ 0.5, d1) + @test count(x -> death(x) == Inf, d1) == 2 + @test count(x -> isapprox(death(x), 0.71, atol=0.1), d1) == 15 + + @test last(only(d2)) == Inf + end + @testset "torus 16, low threshold" begin + d0, d1, d2 = ripserer(torus(16), dim_max=2, threshold=0.5) + + @test length(d0) == 16 + + @test all(x -> birth(x) ≈ 0.5, d1) + @test all(x -> death(x) == Inf, d1) + + @test isempty(d2) + end + @testset "projective plane (modulus), med threshold" begin + _, d1_2, d2_2 = ripserer(projective_plane, + dim_max=2, threshold=1) + _, d1_3, d2_3 = ripserer(projective_plane, + dim_max=2, modulus=3, threshold=1) + @test d1_2 == [PersistenceInterval(1, Inf)] + @test d2_2 == [PersistenceInterval(1, Inf)] + @test isempty(d1_3) + @test isempty(d2_3) + end +end + +@testset "sparse matrix" begin + @testset "icosahedron" begin + flt = SparseRips(icosahedron, threshold=2) + res = ripserer(flt, dim_max=2) + @test res[1] == [fill(PersistenceInterval(0.0, 1.0), 11); + PersistenceInterval(0.0, Inf)] + @test isempty(res[2]) + @test res[3] == [PersistenceInterval(1.0, 2.0)] + end + @testset "torus 16" begin + dists = sparse(torus(16)) + SparseArrays.fkeep!(dists, (_, _, v) -> v ≤ 1) + + d0, d1, d2 = ripserer(dists, dim_max=2) + + @test length(d0) == 16 + + @test all(x -> birth(x) ≈ 0.5, d1) + @test count(x -> death(x) ≈ 1, d1) == 2 + @test count(x -> isapprox(death(x), 0.71, atol=0.1), d1) == 15 + + @test last(only(d2)) == 1 + end + @testset "projective plane (modulus), med threshold" begin + dists = sparse(projective_plane) + SparseArrays.fkeep!(dists, (_, _, v) -> v ≤ 2) + + _, d1_2, d2_2 = ripserer(dists, dim_max=2, threshold=1) + _, d1_3, d2_3 = ripserer(dists, dim_max=2, modulus=3, threshold=1) + @test d1_2 == [PersistenceInterval(1, Inf)] + @test d2_2 == [PersistenceInterval(1, Inf)] + @test isempty(d1_3) + @test isempty(d2_3) + 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} + }[] +end + +@testset "Sublevel using SparseRips" begin + data = [range(0, 1, length=5); + range(1, 0.5, length=5)[2:end]; + range(0.5, 2, length=4)[2:end]; + range(2, -1, length=4)[2:end]] + + # Create distance matrix from data, where neighboring points are connected by edges + # and the edge weights are equal to the max of both vertex births. + n = length(data) + dists = spzeros(n, n) + for i in 1:n + dists[i, i] = data[i] + end + for i in 1:n-1 + j = i + 1 + dists[i, j] = dists[j, i] = max(dists[i, i], dists[j, j]) + end + # 0-dimensional persistence should find values of minima and maxima of our data. + res = first(ripserer(dists, dim_max=0)) + mins = birth.(res) + maxs = death.(filter(isfinite, res)) + @test sort(mins) == [-1.0, 0.0, 0.0, 0.5] + @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] + + d0, d1, d2, d3, d4 = ripserer(Cubical(data), reps=true, dim_max=4) + + @test d0 == [(0, Inf), (1, 2)] + @test d1 == [(0, 2)] + @test isempty(d2) + @test isempty(d3) + @test isempty(d4) + + @test vertices.(representative(d0[1])) == [(i,) for i in 1:length(data)] + @test vertices(only(representative(d0[2]))) == (13,) +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!([ + [(i+1, i) for i in 1:17]; (18, 1) + ]) + + @test_broken ripserer(cycle, cohomology=false, threshold=2)[2][1] == (1.0, Inf) +end + +@testset "Only print to stderr and only when progress is enabled." begin + @suppress begin + @test (@capture_out ripserer(torus(100), dim_max=2)) == "" + @test (@capture_out ripserer(torus(100), dim_max=2, progress=true)) == "" + + @test (@capture_err ripserer(torus(100), dim_max=2)) == "" + @test (@capture_err ripserer(torus(100), dim_max=2, progress=true)) != "" + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 73968791..dbbb3764 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,11 +21,14 @@ end @safetestset "chainelement" begin include("chainelement.jl") end -@safetestset "reduction_structs" begin - include("reduction_structs.jl") +@safetestset "zerodimensional" begin + include("zerodimensional.jl") end -@safetestset "reduction" begin - include("reduction.jl") +@safetestset "reductionmatrix" begin + include("reductionmatrix.jl") +end +@safetestset "ripserer" begin + include("ripserer.jl") end @safetestset "plotting" begin include("plotting.jl") diff --git a/test/simplex.jl b/test/simplex.jl index aa2f0e0b..20edfeba 100644 --- a/test/simplex.jl +++ b/test/simplex.jl @@ -110,4 +110,37 @@ end end end end + + @testset "boundary" begin + @testset "number of faces" begin + for dim in 1:10 + simplex = Simplex{dim}(10, 1) + simplex_vxs = vertices(simplex) + bnd = face_type(simplex)[] + for face in boundary(FakeFiltration(), simplex) + push!(bnd, face) + end + @test length(bnd) == dim + 1 + @test all(isequal(1), diam.(bnd)) + @test all(issubset(vertices(f), simplex_vxs) for f in bnd) + end + end + @testset "type stability" begin + for dim in 1:10 + simplex = Simplex{dim}(10, 1) + @test begin @inferred boundary(FakeFiltration(), simplex); true end + biter = boundary(FakeFiltration(), simplex) + @static if VERSION ≥ v"1.1.0" + @test begin @inferred Union{ + Nothing, + Tuple{face_type(simplex), Int}, + } iterate(biter); true end + @test begin @inferred Union{ + Nothing, + Tuple{face_type(simplex), Int}, + } iterate(biter, 1); true end + end + end + end + end end diff --git a/test/zerodimensional.jl b/test/zerodimensional.jl new file mode 100644 index 00000000..38a6798b --- /dev/null +++ b/test/zerodimensional.jl @@ -0,0 +1,35 @@ +using Ripserer + +using DataStructures +using Ripserer: DisjointSetsWithBirth, find_leaves! + +s = DisjointSetsWithBirth(collect(1:10)) + +@testset "basic tests" begin + @test length(s) == 10 + for i = 1:10 + @test find_root!(s, i) == i + @test find_leaves!(s, i) == [i] + @test birth(s, i) == i + end + @test_throws BoundsError find_root!(s, 11) +end + +@testset "union!, find_root!, find_leaves!" begin + union!(s, 2, 3) + @test find_root!(s, 3) == 2 + @test find_leaves!(s, 2) == [2, 3] + + @test union!(s, 8, 7) == 8 + @test union!(s, 5, 6) == 5 + @test find_leaves!(s, 5) == [5, 6] + @test union!(s, 8, 5) == 8 + @test find_root!(s, 6) == 8 + @test find_leaves!(s, 8) == [5, 6, 7, 8] + union!(s, 2, 6) + @test find_root!(s, 2) == 8 + root1 = find_root!(s, 6) + root2 = find_root!(s, 2) + @test root_union!(s, Int(root1), Int(root2)) == 8 + @test union!(s, 5, 6) == 8 +end