API

ripserer(Type{<:AbstractFiltration}, args...; kwargs...)
ripserer(filtration::AbstractFiltration; kwargs...)

Compute the persistent homology of a filtration. The filtration can be given as an AbstractFiltration type, followed by its arguments, or as an initialized object (see examples below). If only data is given, Rips is used by default.

Returns a Vector of PersistenceDiagrams with (dim_max + 1) elements. The diagrams are sorted by dimension; the first element of the result is the 0-dimensional diagram, and the last is the (dim_max)-dimensional diagram.

Keyword 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: use this type of field of coefficients. Defaults to Ripserer.Mod{modulus}.

• threshold: compute persistent homology up to diameter smaller than threshold. Note that this parameter is passed to the filtration constructor. When using low thresholds with Rips filtrations, consider setting sparse=true for optimal performance.

• cutoff: only keep intervals with persistence(interval) > cutoff. Defaults to 0. When cutoff < 0, the result will also contain zero-length intervals.

• reps: if true, attach representative (co)cycles to persistence intervals. Can also be set to collection of integers to only find representatives in specified dimensions, e.g. reps=1:2 will only find representatives in dimensions 1 and 2. This is useful for large filtrations (such as cubical) where calculating zero-dimensional representatives can be very slow. Defaults to false for cohomology and 1:dim_max for homology. Representatives are wrapped in a Chain.

• verbose: If true, show a verbose bar. Defaults to false.

• alg: select the algorithm used in computation. The options are:

  • :cohomology: Default and fastest algorithm. When reps is set, intervals are equipped with representative cocycles.

  • :homology: Significantly slower than :cohomology, but finds representative cycles. Does not find infinite intervals beyond dimension 0.

  • :involuted: Use cohomology result to compute representative cycles. Can be extremely efficient compared to :homology, especially with Rips filtrations. See this paper for more information.

• implicit: If true, an implicit reduction algorithm is used. Defaults to true for :cohomology and :involuted, and false for :homology. implicit=false is not recommended for :cohomology because it disables the emergent pairs optimization.

Other kwargs... are passed to the filtration.

Examples

julia> ts = range(0, 2π; length=20)[1:(end - 1)];

julia> X = [((2 + cos(θ)) * cos(φ), (2 + cos(θ)) * sin(φ), sin(θ)) for θ in ts for φ in ts];

julia> ripserer(X)
2-element Vector{PersistenceDiagrams.PersistenceDiagram}:
 361-element 0-dimensional PersistenceDiagram
 362-element 1-dimensional PersistenceDiagram

julia> ripserer(EdgeCollapsedRips, X; modulus=7, threshold=2)
2-element Vector{PersistenceDiagrams.PersistenceDiagram}:
 361-element 0-dimensional PersistenceDiagram
 362-element 1-dimensional PersistenceDiagram

julia> ripserer(Rips(X; threshold=1); alg=:involuted)
2-element Vector{PersistenceDiagrams.PersistenceDiagram}:
 361-element 0-dimensional PersistenceDiagram
 362-element 1-dimensional PersistenceDiagram

Rips{I, T} <: AbstractRipsFiltration{I, T}

This type represents a filtration of Vietoris-Rips complexes.

Diagonal items in the input matrix are treated as vertex birth times.

Zero values are not allowed due to how sparse matrices work in Julia. If you need zero birth times, try offseting all values by a constant.

Threshold defaults to the radius of the input space. When using low thresholds, consider using the sparse=true keyword argument. It will give the same result, but may be much faster.

Constructors

• Rips(distance_matrix; threshold=nothing)
• Rips(points; metric=Euclidean(1e-12), threshold=nothing)
• Rips{I}(args...): I sets the size of integer used to represent simplices.

Examples

julia> data = [(sin(t), cos(t)) for t in range(0, 2π, length=101)][1:end-1];

julia> ripserer(Rips(data))
2-element Vector{PersistenceDiagrams.PersistenceDiagram}:
 100-element 0-dimensional PersistenceDiagram
 1-element 1-dimensional PersistenceDiagram

julia> ripserer(Rips(data, threshold=1.7))[2]
1-element 1-dimensional PersistenceDiagram:
 [0.0628, ∞)

julia> using Distances

julia> ripserer(Rips(data, metric=Cityblock()))[2]
1-element 1-dimensional PersistenceDiagram:
 [0.0888, 2.0)

Cubical{T, K} <: AbstractFiltration{CartesianIndex{K}, T}

Cubical is used to compute sublevel persistent homology on N-dimensional images, which are of type AbstractArray{T, N}.

This type uses the CubeMap structure to find birth times of cubes (see reference).

Constructor

• Cubical(image::AbstractArray{T, N}, threshold=maximum(image))

Reference

Wagner, H., Chen, C., & Vuçini, E. (2012). Efficient computation of persistent homology for cubical data. In Topological methods in data analysis and visualization II (pp. 91-106). Springer, Berlin, Heidelberg.

Example

julia> image = [1 0 0; 0 2 0; 0 0 0];

julia> ripserer(Cubical(image))[1]
1-element 0-dimensional PersistenceDiagram:
 [0.0, ∞)

julia> ripserer(Cubical(image))[2]
1-element 1-dimensional PersistenceDiagram:
 [1.0, 2.0)

Custom{I, T} <: AbstractCustomFiltration{I, T}

Build a custom filtration by specifying simplices and their birth times.

The list of simplices is corrected to form a valid filtration; birth times are corrected so a simplex is never born before its faces and missing simplices are added.

See the examples below for construction. Note how the unlisted 0-simplices were added with birth times equal to the lowest between their cofaces. The order in which simplices are given does not matter.

To create your own types of custom filtrations, subtype AbstractCustomFiltration.

Examples

julia> flt = Custom([(1,) => 0, (4,) => 0, (1, 2) => 1, (1, 3) => 2, (1, 4) => 3, (2, 3) => 4, (2, 4) => 5, (3, 4) => 6, (1, 2, 3) => 7, (1, 2, 4) => 8, (1, 3, 4) => 9]; threshold=8)
Custom{Int64, Int64}(nv=4)

julia> flt[0] # Can be indexed with dimension to list simplices
4-element Vector{Simplex{0, Int64, Int64}}:
 +Simplex{0}((4,), 0)
 +Simplex{0}((2,), 1)
 +Simplex{0}((3,), 2)
 +Simplex{0}((1,), 0)

julia> ripserer(flt)[1]
2-element 0-dimensional PersistenceDiagram:
 [0.0, 3.0)
 [0.0, ∞)

julia> ripserer(flt)[2]
3-element 1-dimensional PersistenceDiagram:
 [5.0, 8.0)
 [4.0, 7.0)
 [6.0, ∞)

Alpha{I, P<:SVector} <: AbstractFiltration{I, Float64}

Alpha filtrations are filtrations of the Delaunay complex.

They have much fewer simplices than Rips, so they are efficient even with large datasets, as long as their dimensionality is low. What "low" means depends on the data, but this is definitely a good choice for 3D or lower. For high dimensional data, filtration construction may take a long time.

Constructors

• Alpha(points; threshold, verbose): points should be a vector of Tuples, SVectors or similar.
• Alpha{I}(args...): I sets the size of integer used to represent simplices. Try using I=Int128 if construction complains about overflow.

Reference

Edelsbrunner, H. (1993, July). The union of balls and its dual shape. In Proceedings of the ninth annual symposium on Computational geometry (pp. 218-231).

Example

julia> data = [(sin(t), cos(t), (t - π)^2) for t in range(0, 2π, length=101)[1:end-1]];

julia> alpha = Alpha(data)
Alpha{Int64, Float64}(nv=100)

julia> rips = Rips(data)
Rips{Int64, Float64}(nv=100, sparse=false)

julia> length(Ripserer.edges(alpha))
197

julia> length(Ripserer.edges(rips))
3613

julia> sort(ripserer(alpha)[2], by=persistence)[end]
[0.375, 2.01) with:
 birth_simplex: Ripserer.Simplex{1, Float64, Int64}
 death_simplex: Ripserer.Simplex{2, Float64, Int64}

julia> sort(ripserer(rips)[2], by=persistence)[end]
[0.375, 2.01) with:
 birth_simplex: Ripserer.Simplex{1, Float64, Int64}
 death_simplex: Ripserer.Simplex{2, Float64, Int64}

EdgeCollapsedRips{I, T} <: AbstractRipsFiltration{I, T}

Perform a sequence of edge collapses on a filtration. This may significantly reduce computation time and does not change the result. The speedup is especially apparent with datasets that have a boundary, and with high-dimensional persistent homology computation.

The drawback is that doing the collapses themselves can be time-consuming. The construction does not require a lot of memory ($\mathcal{O}(n^2)$ for $n$ vertices). This might still be a good choice for large inputs if you are willing to wait but don't have enough memory to compute persistent homology with Rips.

See the reference below for a description of the algorithm.

Constructors

• EdgeCollapsedRips(::AbstractRipsFiltration; verbose=false, threshold=nothing ): Collapse a given filtration. Setting verbose shows a progress bar.

• EdgeCollapsedRips(::EdgeCollapsedRips; verbose=false, threshold=nothing): Allows changing I or threshold without recomputing.

• EdgeCollapsedRips(arg; kwargs...): Use arg and kwargs to construct a Rips filtration and collapse it

• EdgeCollapsedRips{I}(arg; kwargs...): Change the index type used to represent simplices. May be necessary for large inputs and high-dimensional computation.

Examples

julia> using Random; Random.seed!(1337);

julia> data = [tuple(rand(6)...) for _ in 1:100];

julia> rips = Rips(data)
Rips{Int64, Float64}(nv=100, sparse=false)

julia> length(Ripserer.edges(rips))
3934

julia> collapsed = EdgeCollapsedRips(data) # or EdgeCollapsedRips(rips)
EdgeCollapsedRips{Int64, Float64}(nv=100)

julia> length(Ripserer.edges(collapsed))
1324

julia> ripserer(rips) == ripserer(collapsed)
true

julia> ripserer(collapsed; dim_max=4)
5-element Vector{PersistenceDiagrams.PersistenceDiagram}:
 100-element 0-dimensional PersistenceDiagram
 58-element 1-dimensional PersistenceDiagram
 35-element 2-dimensional PersistenceDiagram
 10-element 3-dimensional PersistenceDiagram
 4-element 4-dimensional PersistenceDiagram

Reference

Boissonnat, J. D., & Pritam, S. (2019). Edge Collapse and Persistence of Flag Complexes.

PersistenceDiagram <: AbstractVector{PersistenceInterval}

Type for representing persistence diagrams. Behaves exactly like a vector of PersistenceIntervals, but can have additional metadata attached to it. It supports pretty printing and plotting.

Can be used as a table with any function that uses the Tables.jl interface. Note that using it as a table will only keep interval endpoints and the dim and threshold attributes.

Example

julia> diagram = PersistenceDiagram([(1, 3), (3, 4), (1, Inf)]; dim=1, custom_metadata=:a)
3-element 1-dimensional PersistenceDiagram:
 [1.0, 3.0)
 [3.0, 4.0)
 [1.0, ∞)

julia> diagram[1]
[1.0, 3.0)

julia> sort(diagram; by=persistence, rev=true)
3-element 1-dimensional PersistenceDiagram:
 [1.0, ∞)
 [1.0, 3.0)
 [3.0, 4.0)

julia> propertynames(diagram)
(:intervals, :dim, :custom_metadata)

julia> dim(diagram)
1

julia> diagram.custom_metadata
:a

PersistenceInterval

Type for representing persistence intervals. It behaves exactly like a Tuple{Float64, Float64}, but can have meta data attached to it. The metadata is accessible with getproperty or the dot syntax.

Example

julia> interval = PersistenceInterval(1, Inf; meta1=:a, meta2=:b)
[1.0, ∞) with:
 meta1: Symbol
 meta2: Symbol

julia> birth(interval), death(interval), persistence(interval)
(1.0, Inf, Inf)

julia> isfinite(interval)
false

julia> propertynames(interval)
(:birth, :death, :meta1, :meta2)

julia> interval.meta1
:a

birth(interval)

Get the birth time of interval.

death(interval)

Get the death time of interval.

persistence(interval)

Get the persistence of interval, which is equal to death - birth.

midlife(interval)

Get the midlife of the interval, which is equal to (birth + death) / 2.

representative(interval::PersistenceInterval)

Get the representative (co)cycle attached to interval, if it has one.

birth_simplex(interval::PersistenceInterval)

Get the critical birth simplex of interval, if it has one.

death_simplex(interval::PersistenceInterval)

Get the critical death simplex of interval, if it has one.

barcode(diagram)

Plot the barcode plot of persistence diagram or multiple diagrams in a collection. The infinity keyword argument determines where the infinity line is placed. If unset, the function tries to use threshold(diagram), or guess a good position to place the line at.

Simplex{D, T, I<:Integer} <: AbstractSimplex{D, T, I}

The vanilla simplex type represented by dimension D, an index of type I, and a birth time of type T.

Constructors

• Simplex{D[, T, I]}(index, birth)
• Simplex{D}(vertices, birth): vertices must be sorted descending. This constructor mainly exists for debugging purposes. Using simplex is usually the better option.

Examples

julia> sx = Simplex{2}(2, 1)
2-dimensional Simplex(index=2, birth=1):
  +(4, 2, 1)

julia> index(sx)
2

julia> vertices(sx)
(4, 2, 1)

julia> birth(sx)
1

Cube{D, T, K} <: AbstractCell{D, T, CartesianIndex{K}}

A Cube is similar to a Simplex, but it has 2^D vertices instead of D+1. The vertices are encoded as the position in the CubeMap (see reference in Cubical). A Cube's vertices are of type CartesianIndex{K}.

Example

julia> Cube{1}(CartesianIndex(1, 2), 1.0)
Cube{1, Float64, 2}((1, 2), 1.0)

dim(::AbstractCell)
dim(::Type{<:AbstractCell})

Get the dimension of a cell i.e. the value of D. Can also be called on the type.

Examples

julia> dim(Simplex{2}((3, 2, 1), 3.2))
2

julia> dim(Cube{3, Int, 4})
3

birth(σ::AbstractCell)

Get the birth time of σ, i.e. the time it first appears in the filtration.

Example

julia> birth(Simplex{2}((3, 2, 1), 3.2))
3.2

index(σ::AbstractCell)

Get the combinatorial index of the σ. The index can be any type, but should uniquely identify a cell. It is also used to break ties when comparing simplices with the same birth time.

julia> index(Simplex{2}((3, 2, 1), 3.2))
1

vertices(σ::AbstractCell{D,T,I})::NTuple{length(σ),I}

Get the vertices of σ.

Example

julia> vertices(Simplex{2}((3, 2, 1), 3.2))
(3, 2, 1)

Chain{F,S,E} <: AbstractVector{ChainElement{S,F}}

An internal representation of a chain. Behaves like an array of pairs S => F, where S is the simplex type, and F is the coefficient type (probably a subtype of Mod).

Most functions that can be called on AbstractCells can also be called on the elements of a Chain.

Examples

julia> data = [(rand(), rand(), rand()) for _ in 1:100];

julia> result = ripserer(data; reps=true, modulus=7);

julia> chain = result[end][end].representative
28-element Chain{Ripserer.Mod{7},Ripserer.Simplex{1, Float64, Int64}}:
 +Simplex{1}((68, 54), 0.25178097927369875) => 6 mod 7
 +Simplex{1}((54, 46), 0.2575262844682746) => 1 mod 7
 +Simplex{1}((88, 56), 0.25936896586973557) => 1 mod 7
 +Simplex{1}((79, 22), 0.26690101517910964) => 6 mod 7
 +Simplex{1}((53, 13), 0.27161939814693414) => 6 mod 7
 +Simplex{1}((24, 13), 0.28686687771884056) => 6 mod 7
 +Simplex{1}((56, 54), 0.2869400047523196) => 6 mod 7
 +Simplex{1}((50, 42), 0.29274772369750884) => 6 mod 7
 +Simplex{1}((93, 50), 0.3022995886824861) => 1 mod 7
 +Simplex{1}((73, 54), 0.30771273548368216) => 6 mod 7
 ⋮
 +Simplex{1}((70, 34), 0.3412250751087395) => 6 mod 7
 +Simplex{1}((49, 42), 0.3414959312908324) => 6 mod 7
 +Simplex{1}((49, 7), 0.35568373080560794) => 6 mod 7
 +Simplex{1}((88, 46), 0.35629002982624475) => 1 mod 7
 +Simplex{1}((95, 88), 0.3566986102316938) => 6 mod 7
 +Simplex{1}((95, 54), 0.363029725360828) => 6 mod 7
 +Simplex{1}((50, 13), 0.3648936036309352) => 6 mod 7
 +Simplex{1}((24, 3), 0.3653704317993549) => 6 mod 7
 +Simplex{1}((85, 42), 0.369953146904745) => 6 mod 7

julia> simplex.(chain)
28-element Vector{Simplex{1, Float64, Int64}}:
 +Simplex{1}((68, 54), 0.25178097927369875)
 +Simplex{1}((54, 46), 0.2575262844682746)
 +Simplex{1}((88, 56), 0.25936896586973557)
 +Simplex{1}((79, 22), 0.26690101517910964)
 +Simplex{1}((53, 13), 0.27161939814693414)
 +Simplex{1}((24, 13), 0.28686687771884056)
 +Simplex{1}((56, 54), 0.2869400047523196)
 +Simplex{1}((50, 42), 0.29274772369750884)
 +Simplex{1}((93, 50), 0.3022995886824861)
 +Simplex{1}((73, 54), 0.30771273548368216)
 ⋮
 +Simplex{1}((70, 34), 0.3412250751087395)
 +Simplex{1}((49, 42), 0.3414959312908324)
 +Simplex{1}((49, 7), 0.35568373080560794)
 +Simplex{1}((88, 46), 0.35629002982624475)
 +Simplex{1}((95, 88), 0.3566986102316938)
 +Simplex{1}((95, 54), 0.363029725360828)
 +Simplex{1}((50, 13), 0.3648936036309352)
 +Simplex{1}((24, 3), 0.3653704317993549)
 +Simplex{1}((85, 42), 0.369953146904745)

julia> vertices.(chain)
28-element Vector{Tuple{Int64, Int64}}:
 (68, 54)
 (54, 46)
 (88, 56)
 (79, 22)
 (53, 13)
 (24, 13)
 (56, 54)
 (50, 42)
 (93, 50)
 (73, 54)
 ⋮
 (70, 34)
 (49, 42)
 (49, 7)
 (88, 46)
 (95, 88)
 (95, 54)
 (50, 13)
 (24, 3)
 (85, 42)

julia> coefficient.(chain)
28-element Vector{Mod{7}}:
 6 mod 7
 1 mod 7
 1 mod 7
 6 mod 7
 6 mod 7
 6 mod 7
 6 mod 7
 6 mod 7
 1 mod 7
 6 mod 7
       ⋮
 6 mod 7
 6 mod 7
 6 mod 7
 1 mod 7
 6 mod 7
 6 mod 7
 6 mod 7
 6 mod 7
 6 mod 7

Mod{M} <: Integer

Mod{M} is the default field used by Ripserer. It is a representation of a finite field $\mathbb{Z}_M$, integers modulo small, prime M. Supports field arithmetic and can be converted to integer with Int.

Its values are not comparable on purpose.

Example

julia> Mod{3}(5)
2 mod 3

julia> Mod{3}(5) + 1
0 mod 3

RipsPersistentHomology

Compute Vietoris-Rips persistent homology and convert the results to a table of continuous values.

!!! warning Warning Computing Vietoris-Rips persistent homology may be CPU and memory intensive even for modestly-sized data sets. Consider using AlphaPersistentHomology for low-dimensional data.

Hyperparameters

• vectorizer = PersistenceImageVectorizer(): AbstractVectorizer used to transform persistence diagrams to continuous vectors. See the PersistenceDiagrams.jl Documentation for more information.

• dim_max = 1: compute persistent homology up to this dimension.

• modulus = 2: compute persistent homology with coefficients in the prime field of integers mod modulus.

• threshold = nothing: compute persistent homology up to diameter smaller than threshold.

• cutoff = 0: only keep intervals with persistence(interval) > cutoff.

• sparse = false: use a sparse Rips filtration.

• collapse = false: use the EdgeCollapsedRips filtration. This is usually a good choice for computations with a high dim_max.

See also

• ripserer
• Rips
• EdgeCollapsedRips

AlphaPersistentHomology

Compute alpha complex persistent homology and convert the results to a table of continuous values.

!!! warning Warning Using high-dimensional data with this model may be computationaly expensive. Consider using RipsPersistentHomology.

Hyperparameters

• vectorizer = PersistenceImageVectorizer(): AbstractVectorizer used to transform persistence diagrams to continuous vectors. See the PersistenceDiagrams.jl Documentation for more information.

• dim_max = 1: compute persistent homology up to this dimension.

• modulus = 2: compute persistent homology with coefficients in the prime field of integers mod modulus.

• threshold = nothing: compute persistent homology up to diameter smaller than threshold.

• cutoff = 0: only keep intervals with persistence(interval) > cutoff.

See also

• ripserer
• Alpha

CubicalPersistentHomology

Compute cubical persistent homology and convert the results to a table of continuous values.

Hyperparameters

• vectorizer = PersistenceImageVectorizer(): AbstractVectorizer used to transform persistence diagrams to continuous vectors. See the PersistenceDiagrams.jl Documentation for more information.

• dim_max = 1: compute persistent homology up to this dimension.

• threshold = nothing: compute persistent homology up to diameter smaller than threshold.

• cutoff = 0: only keep intervals with persistence(interval) > cutoff.

• negate = false: negate the image before computation.

See also

• ripserer
• Cubical

reconstruct_cycle(filtration, interval[, t]; distances=distance_matrix(filtration))

Reconstruct the shortest representative cycle for the first homology group of given interval. The optional argument t sets the time at which the cycle is to be computed. It defaults to interval birth time, which gives a cycle similar to a representative cycle computed from homology. In general, higher times will yield nicer cycles. t can be a simplex or a number.

The optional distances keyword argument can be used to change the distance matrix used for determining edge lengths.

This method uses the representative cocycle to compute the cycle. As such, the interval must be computed with the default cohomology algorithm and must include a representative. To get such an interval, run ripserer with the keyword argument reps=true or reps=1.

Note that this method only works in the first dimension, as it is based on finding shortest paths in a graph.

module Partition

This submodule contains the following partition functions to be used with circular coordinates.

• Partition.linear(r, d) = max(r - d, 0.0)
• Partition.quadratic(r, d) = (max(r - d, 0.0))^2
• Partition.exponential(r, d) = r - d > 0 ? exp(r^2/(d^2 - r^2)) : 0.0

CircularCoordinates

This struct implements Perea's sparse circular coordinates (see reference below).

The idea behind this method is that you pick a subset of your data set (also called landmarks), compute persistent cohomology on that subset, and use the result to construct circular coordinates on the whole data set.

For this to work correctly, you need the following.

• A set of landmarks that cover the points well. The default minmax sampling method is

usually a good choice here.

• The persistence diagram on the landmarks must have an interval that is persistent enough.

The circular coordinates returned are on the interval [0, 1). If you are looking for angles, make sure to multiply them by 2π.

If you try to compute circular coordinates for a point that is not near any landmark, the result will be missing.

Constructor

CircularCoordinates([::AbstractFiltration, ]points, landmarks; kwargs...)

Arguments

• points: a vector of points. The eltype of this vector can be Tuples, SVectors, or similar.

• landmarks: can be an integer, a vector of indices, or a vector of points. If set to an integer, it sets the number of landmarks to choose with maxmin sampling. If you are looking for non-sparse circular coordinates, use landmarks=eachindex(points).

• out_dim: number of most persistent persistence intervals to use to compute coordinates. If less than out_dim suitable intervals are found, the construction will show a warning and construct CircularCoordinates with a lower out_dim. This warning can be suppressed by passing warn=false.

• partition: a function that defines the partition of unity used when determining the coordinates. Should take two arguments, r, the radius of the balls, and d, the distance from the landmark. The partition function should only have support on the ball around the landmark (in other words, it should evaluate to 0 when d ≥ r). If that is not the case, the circular coordinates are not well-defined and the results may not always make sense. See Partition for a small collection of predefined functions.

• metric: a metric from Distances.jl. Can only be set if the filtration argument is Rips.

• kwargs...: additional keyword arguments passed to ripserer. Note that modulus is set to a random prime between 7 and 79 by default.

Example

julia> data = [(sin(t), cos(t)) for t in range(0, 2π, length=101)[1:end-1]];

julia> cc = CircularCoordinates(data, 1:10:100)
CircularCoordinates(
  out_dim=1,
  radius=(0.9510565162951535,),
  n_landmarks=10,
  partition=linear,
  metric=Distances.Euclidean(0.0),
)

julia> summary(cc(data))
"100×1 Matrix{Union{Missing, Float64}}"

julia> summary(cc(data, 1))
"100-element Vector{Union{Missing, Float64}}"

Reference

Perea, J. A. (2020, June). Sparse circular coordinates via principal Z-bundles. In Topological Data Analysis: The Abel Symposium 2018 (Vol. 15, p. 435). Springer Nature.

AbstractFiltration{I,T}

A filtration is used to find the edges in filtration and to create simplices. An AbstractFiltration{I,T}'s simplex type is expected to return simplices of type <:AbstractCell{_,T,I}.

An AbstractFiltration constructor must accept the verbose keyword argument.

Interface

• nv(::AbstractFiltration)
• births(::AbstractFiltration)
• vertices(::AbstractFiltration)
• edges(::AbstractFiltration)
• simplex_type(::Type{AbstractFiltration}, dim)
• simplex(::AbstractFiltration, ::Val{dim}, vertices)
• unsafe_simplex(::AbstractFiltration, ::Val{dim}, vertices)
• unsafe_cofacet(::AbstractFiltration, simplex, vertices, vertex[, edges])
• threshold(::AbstractFiltration)
• columns_to_reduce(::AbstractFiltration, ::Any)
• emergent_pairs(::AbstractFiltration)
• postprocess_diagram(::AbstractFiltration, ::Any)

nv(::AbstractFiltration)

Return the number of vertices in filtration.

Example

julia> Ripserer.nv(Rips([1 1; 1 1]))
2

births(::AbstractFiltration)

Get the birth times of vertices in filtration. Defaults to all births being 0. Must return array of the same shape as the filtration's vertices.

Examples

julia> flt = Rips([1 1 2; 1 0 1; 2 1 0]);

julia> Ripserer.births(flt)
3-element view(::Vector{Int64}, 1:4:9) with eltype Int64:
 1
 0
 0

vertices(::AbstractFiltration)

Return the vertices in filtration. Defaults to 1:n. The shape and eltype of the result can be anything as long as result[result[i]] == result[i] holds.

Example

julia> vertices(Rips([0 1 1; 1 0 1; 1 1 0]))
Base.OneTo(3)

julia> vertices(Cubical([0 1 1; 1 0 1; 1 1 0]))
3×3 CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)  CartesianIndex(1, 3)
 CartesianIndex(2, 1)  CartesianIndex(2, 2)  CartesianIndex(2, 3)
 CartesianIndex(3, 1)  CartesianIndex(3, 2)  CartesianIndex(3, 3)

edges(::AbstractFiltration)

Get edges (1-simplices) in filtration. Edges should be of type simplex_type(filtration, 1).

Example

julia> Ripserer.edges(Rips([0 2 1; 2 0 1; 1 1 0], threshold=2))
3-element Array{Simplex{1,Int64,Int64},1}:
 +Simplex{1}([2, 1], 2)
 +Simplex{1}([3, 1], 1)
 +Simplex{1}([3, 2], 1)

simplex_type(::Type{<:AbstractFiltration}, D)
simplex_type(::AbstractFiltration, D)

Return the D-dimensional simplex type in the filtration. Only the method for the type needs to be overloaded.

Examples

julia> Ripserer.simplex_type(Rips{Int,Float64}, 1)
Simplex{1, Float64, Int64}

julia> Ripserer.simplex_type(Cubical{2,Float16}, 2)
Cube{2, Float16, 2}

 simplex(::AbstractFiltration, ::Val{D}, vertices)

Return D-simplex constructed from vertices. Return nothing if simplex is not in filtration. This function is safe to call with vertices that are out of order. Default implementation sorts vertices and calls unsafe_simplex.

Example

julia> simplex(Rips([0 2 1; 2 0 1; 1 1 0], threshold=2), Val(1), (1, 2))
1-dimensional Simplex(index=1, birth=2):
  +[2, 1]

unsafe_simplex(::AbstractFiltration, ::Val{D}, vertices)

Return D-simplex constructed from vertices. Return nothing if simplex is not in filtration. The unsafe in the name implies that it's up to the caller to ensure vertices are sorted and unique.

unsafe_cofacet(filtration, simplex, cofacet_vertices, v[, edges])
unsafe_cofacet(::Type{S}, filtration, simplex, cofacet_vertices, v[, edges])

Return cofacet of simplex with vertices equal to cofacet_vertices. v is the vertex that was added to construct the cofacet. In the case of sparse rips filtrations, an additional argument edges is used. edges is a vector that contains the weights on edges connecting the new vertex to old vertices. S is the simplex type which can be used for dispatch.

The unsafe in the name implies that it's up to the caller to ensure vertices are sorted and unique.

Default implementation uses unsafe_simplex.

threshold(::AbstractFiltration)

Get the threshold of filtration. This is the maximum diameter a simplex in the filtration can have. Defaults to Inf.

Examples

julia> threshold(Rips([0 2 1; 2 0 1; 1 1 0]))
1

julia> threshold(Cubical([1 1 2; 3 2 1; 0 0 0]))
3

columns_to_reduce(::AbstractFilration, prev_column_itr)

List all columns to reduce in next dimension, possibly computing it from previous columns. Default implementation uses coboundary with the all cofacets parameter set to Val(false).

Example

julia> flt = Rips([0 1 1; 1 0 1; 1 1 0]);

julia> Ripserer.columns_to_reduce(flt, Ripserer.edges(flt)) |> collect
1-element Vector{Simplex{2, Int64, Int64}}:
 +Simplex{2}((3, 2, 1), 1)

emergent_pairs(::AbstractFiltration)

Return true if the emergent pairs optimization is to be performed. Default to returning true. Should be set to false for a filtration type that is unable to produce (co)boundary simplices in the correct order.

postprocess_diagram(::AbstractFiltration, diagram)

This function is called on each resulting persistence diagram after all intervals have been computed. Defaults to not doing anything.

distance_matrix(::AbstractFiltration)

Return a matrix with distances between vertices of the filtration. These distances are used to determine edge length when finding shortest represenatative cycle in reconstruct_cycle.

Defaults to all distances being 1.

Examples

julia> flt = Rips([1 1 2; 1 0 1; 2 1 0]);

julia> Ripserer.distance_matrix(flt)
3×3 Matrix{Int64}:
 1  1  2
 1  0  1
 2  1  0

julia> flt = Custom([(1,2,3,4) => 10.0, (1,2) => 3.0, (1,3) => 4.0, (2,3) => 5.0]);

julia> Ripserer.distance_matrix(flt)
4×4 Ripserer.DefaultDist{Float64}:
 0.0  1.0  1.0  1.0
 1.0  0.0  1.0  1.0
 1.0  1.0  0.0  1.0
 1.0  1.0  1.0  0.0

AbstractRipsFiltration{I<:Signed, T} <: AbstractFiltration{I, T}

An abstract Vietoris-Rips filtration. Its subtypes can only overload adjacency_matrix and get default implementations for the rest of the filtration interface.

Example

julia> struct MyRips <: Ripserer.AbstractRipsFiltration{Int, Float16} end

julia> Ripserer.adjacency_matrix(::MyRips) = [0 1 1; 1 0 1; 1 1 0]

julia> ripserer(MyRips())
2-element Vector{PersistenceDiagrams.PersistenceDiagram}:
 3-element 0-dimensional PersistenceDiagram
 0-element 1-dimensional PersistenceDiagram

adjacency_matrix(::AbstractFiltration)

Return the adjacency matrix. For sparse filtrations, this should return a SparseMatrixCSC.

Examples

julia> Ripserer.adjacency_matrix(Rips([0 2 1; 2 0 1; 1 1 0]))
3×3 Matrix{Int64}:
 0  2  1
 2  0  1
 1  1  0

julia> Ripserer.adjacency_matrix(Rips([0 10 2; 10 0 1; 2 1 0]; sparse=true))
3×3 SparseArrays.SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 ⋅  ⋅  2
 ⋅  ⋅  1
 2  1  ⋅

julia> Ripserer.adjacency_matrix(Custom([(2, 1) => 1, (5, 1) => 2, (3, 4) => 3]))
5×5 SparseArrays.SparseMatrixCSC{Bool, Int64} with 6 stored entries:
 ⋅  1  ⋅  ⋅  1
 1  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  1  ⋅
 ⋅  ⋅  1  ⋅  ⋅
 1  ⋅  ⋅  ⋅  ⋅

abstract type AbstractCustomFiltration{I, T} <: AbstractFiltration{I, T}

This abstract type is for filtrations that have all simplices stored in Dicts. The dicts should be accessible by the function simplex_dicts and should be a vector of Dict{I, T}. A custom filtration should also have adjacency_matrix defined. This matrix is only used as an adjacency matrix. Its values are ignored.

simplex_dicts(::AbstractCustomFiltration)

Get the dictionaries used to get simplex birth times. Should return a Vector of Dict{I, T} that maps a simplex index to its birth time. The first element of this Vector corresponds to vertices, second to 1-simplices etc.

AbstractCell{D, T, I}

An abstract type for representing simplices and other types of cells, such as cubes. A cell is determined by the following.

• Its dimension encoded in the type parameter D. Can be accessed with dim

• The index of type I is used to differentiate between different cells. In general, two cells of the same type (and hence dimension) and the same index are considered to be equal.

• The birth of type T determines when the cell enters the filtration. Note that two cells with the same index should also have the same birth.

• vertices should return the cell's vertices as a tuple. For 0-cells (vertices), the vertices are also used to index into a filtration's vertices.

• The sign determens its orientation. Note that cell == -cell.

Interface

• birth(::AbstractCell)::T
• index(::AbstractCell)::I
• vertices(σ::AbstractCell)::NTuple{length(σ),I}
• length(::Type{AbstractCell})
• sign(::AbstractCell)
• Base.:-(::AbstractCell)
• coboundary(::Any, ::AbstractCell)
• boundary(::Any, ::AbstractCell)

AbstractSimplex{D, T, I<:Integer} <: AbstractCell{I}

An abstract type for representing simplices. A simplex is an AbstractCell with an integer index. A D-simplex has D + 1 vertices.

If a type implements the interface below, default implementations of boundary, coboundary, and vertices are provided.

Interface

• AbstractSimplex{0}(::I, ::T)
• AbstractSimplex{1}(::I, ::T)
• birth(::AbstractSimplex)
• index(::AbstractSimplex)
• sign(::AbstractSimplex)
• Base.:-(::AbstractSimplex)
• unsafe_simplex
• unsafe_cofacet

sign(σ::AbstractCell)

Get the orientation of σ. Should return -1 or 1.

Examples

julia> sign(Simplex{2}((3, 2, 1), 3.2))
1

julia> sign(-Simplex{2}((3, 2, 1), 3.2))
-1

-(σ::AbstractCell)

Reverse the cell orientation.

Example

julia> -Simplex{2}((3, 2, 1), 3.2)
2-dimensional Simplex(index=1, birth=3.2):
  -(3, 2, 1)

coboundary(filtration, simplex[, Val{all_cofacets}])

Iterate over the coboundary of simplex by decreasing index. Use the filtration to determine the diameters and validity of cofacets.

If all_cofacets is false, only return cofaces with vertices added to the beginning of vertex list. The method with all_cofacets only has to be implemented if the filtration does not overload columns_to_reduce.

Comes with a default implementation.

Examples

filtration = Rips([0 1 1 1; 1 0 1 1; 1 1 0 1; 1 1 1 0])

for c in Ripserer.coboundary(filtration, Simplex{1}(2, 1))
    println(c)
end

# output

+Simplex{2}((4, 3, 1), 1)
-Simplex{2}((3, 2, 1), 1)

for c in Ripserer.coboundary(filtration, Simplex{1}(2, 1), Val(false))
    println(c)
end

# output

+Simplex{2}((4, 3, 1), 1)

boundary(filtration, simplex[, Val{all_cofacets}])

Iterate over the boundary of simplex by increasing index. Use the filtration to determine the diameters and validity of cofacets.

Comes with a default implementation.

Example

filtration = Rips([0 1 1 1; 1 0 1 1; 1 1 0 1; 1 1 1 0])

for f in Ripserer.boundary(filtration, Simplex{2}(2, 1))
    println(f)
end

# output

+Simplex{1}((2, 1), 1)
-Simplex{1}((4, 1), 1)
+Simplex{1}((4, 2), 1)