Usage Guide

In this example, we will present the basics of using Ripserer. We start by loading some packages.

using Distances
using Plots
using Ripserer

Using Ripserer With Point Cloud Data

Let's start with generating some points, randomly sampled from a noisy circle.

function noisy_circle(n; r=1, noise=0.1)
    points = NTuple{2,Float64}[]
    for _ in 1:n
        θ = 2π * rand()
        push!(points, (r * sin(θ) + noise * rand(), r * cos(θ) + noise * rand()))
    end
    return points
end

circ_100 = noisy_circle(100)
scatter(circ_100; aspect_ratio=1, legend=false, title="Noisy Circle")

To compute the Vietoris-Rips persistent homology of this data set, run the following.

ripserer(circ_100)

2-element Vector{PersistenceDiagram}:
 100-element 0-dimensional PersistenceDiagram
 2-element 1-dimensional PersistenceDiagram

You can use the dim_max argument to set the maximum dimension persistent homology is computed in.

result_rips = ripserer(circ_100; dim_max=3)

4-element Vector{PersistenceDiagram}:
 100-element 0-dimensional PersistenceDiagram
 2-element 1-dimensional PersistenceDiagram
 8-element 2-dimensional PersistenceDiagram
 1-element 3-dimensional PersistenceDiagram

The result can be plotted as a persistence diagram or as a barcode.

plot(result_rips)
barcode(result_rips)

We can also plot a single diagram or a subset of all diagrams in the same manner. Keep in mind that the result is just a vector of PersistenceDiagrams. The zero-dimensional diagram is found at index 1.

plot(result_rips[2])
barcode(result_rips[2:end]; linewidth=2)

Plotting can be further customized using the standard attributes from Plots.jl.

plot(result_rips; markeralpha=1, markershape=:star, color=[:red, :blue, :green, :purple])

Changing Filtrations

By default, calling ripserer will compute persistent homology with the Rips filtration. To use a different filtration, we have two options.

The first option is to pass the filtration constructor as the first argument. Any keyword arguments the filtration accepts can be passed to ripserer and it will be forwarded to the constructor.

ripserer(EdgeCollapsedRips, circ_100; threshold=1, dim_max=3, metric=Euclidean())

4-element Vector{PersistenceDiagram}:
 100-element 0-dimensional PersistenceDiagram
 2-element 1-dimensional PersistenceDiagram
 0-element 2-dimensional PersistenceDiagram
 0-element 3-dimensional PersistenceDiagram

The second option is to initialize the filtration object first and use that as an argument to ripserer. This can be useful in cases where constructing the filtration takes a long time.

collapsed_rips = EdgeCollapsedRips(circ_100; threshold=1, metric=Euclidean())
ripserer(collapsed_rips; dim_max=3)

4-element Vector{PersistenceDiagram}:
 100-element 0-dimensional PersistenceDiagram
 2-element 1-dimensional PersistenceDiagram
 0-element 2-dimensional PersistenceDiagram
 0-element 3-dimensional PersistenceDiagram

Distance Matrix Inputs

In the previous example, we got our result by passing a collection of points to ripserer. Under the hood, Rips and EdgeCollapsedRips actually work with distance matrices. Let's define a distance matrix of the shortest paths on a regular icosahedron graph.

icosahedron = [
    0 1 2 2 1 2 1 1 2 2 1 3
    1 0 3 2 1 1 2 1 2 1 2 2
    2 3 0 1 2 2 1 2 1 2 1 1
    2 2 1 0 3 2 1 1 2 1 2 1
    1 1 2 3 0 1 2 2 1 2 1 2
    2 1 2 2 1 0 3 2 1 1 2 1
    1 2 1 1 2 3 0 1 2 2 1 2
    1 1 2 1 2 2 1 0 3 1 2 2
    2 2 1 2 1 1 2 3 0 2 1 1
    2 1 2 1 2 1 2 1 2 0 3 1
    1 2 1 2 1 2 1 2 1 3 0 2
    3 2 1 1 2 1 2 2 1 1 2 0
]

To compute the persistent homology, simply feed the distance matrix to ripserer.

result_icosa = ripserer(icosahedron; dim_max=2)

3-element Vector{PersistenceDiagram}:
 12-element 0-dimensional PersistenceDiagram
 0-element 1-dimensional PersistenceDiagram
 1-element 2-dimensional PersistenceDiagram

Thresholding

In our next example, we will show how to use thresholding to speed up computation. We start by defining a sampling function that generates $n$ points from the square $[-4,4]\times[-4,4]$ with a circular hole of radius 1 in the middle.

function cutout(n)
    points = NTuple{2,Float64}[]
    while length(points) < n
        x, y = (8rand() - 4, 8rand() - 4)
        if x^2 + y^2 > 1
            push!(points, (x, y))
        end
    end
    return points
end

cutout (generic function with 1 method)

We sample 2000 points from this space.

cutout_2000 = cutout(2000)
scatter(cutout_2000; markersize=1, aspect_ratio=1, legend=false, title="Cutout")

We calculate the persistent homology and time the calculation.

@time result_cut = ripserer(cutout_2000)

  4.914381 seconds (11.30 k allocations: 630.610 MiB, 2.70% gc time)

plot(result_cut)

Notice that while there are many 1-dimensional classes, one of them stands out. This class represents the hole in the middle of our square. Since the intervals are sorted by persistence, we know the last interval in the diagram will be the most persistent.

most_persistent = result_cut[2][end]

[0.233, 1.86) with:
 birth_simplex: Simplex{1, Float64, Int64}
 death_simplex: Simplex{2, Float64, Int64}

Notice the death time of this interval is around 1.83 and that no intervals occur after 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_thresh_2 = ripserer(cutout_2000; threshold=2)

  1.796538 seconds (94.98 k allocations: 134.316 MiB, 15.59% gc time, 3.42% compilation time)

plot(result_cut_thresh_2; title="Persistence Diagram, threshold=2")