diff --git a/CHANGELOG.md b/CHANGELOG.md index 4198acd0..98341e4d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,7 @@ +# v0.16.3 + +New function: `midlife`. + # v0.16.0 External changes: diff --git a/Project.toml b/Project.toml index 290b8470..ae5c0bf7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Ripserer" uuid = "aa79e827-bd0b-42a8-9f10-2b302677a641" authors = ["mtsch "] -version = "0.16.2" +version = "0.16.3" [deps] Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" @@ -26,7 +26,7 @@ Distances = "0.8, 0.9, 0.10" IterTools = "1" LightGraphs = "1.3.3" MiniQhull = "0.2" -PersistenceDiagrams = "0.8" +PersistenceDiagrams = "^0.8.2" ProgressMeter = "1" RecipesBase = "1" StaticArrays = "0.12, 1" diff --git a/docs/Project.toml b/docs/Project.toml index ad33d460..2b95603f 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" GLMNet = "8d5ece8b-de18-5317-b113-243142960cc6" GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" diff --git a/docs/make.jl b/docs/make.jl index 414dad96..712719a6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,5 +1,6 @@ @info "build started." using Documenter +using Distances using Literate using Ripserer using PersistenceDiagrams @@ -25,9 +26,9 @@ makedocs(; ), pages=[ "Home" => "index.md", - "API" => ["Public" => "api/ripserer.md", "Interfaces" => "api/extensions.md"], + "Usage Guide" => "generated/basics.md", + "API" => "api.md", "Examples" => [ - "generated/basics.md", "generated/stability.md", "generated/cocycles.md", "generated/cubical.md", diff --git a/docs/src/api.md b/docs/src/api.md new file mode 100644 index 00000000..13274269 --- /dev/null +++ b/docs/src/api.md @@ -0,0 +1,227 @@ +# API + +## Ripserer + +```@docs +ripserer +``` + +## Filtrations + +```@docs +Rips +``` + +```@docs +Cubical +``` + +```@docs +Custom +``` + +```@docs +Alpha +``` + +```@docs +EdgeCollapsedRips +``` + +## Persistence Diagrams + +Persistence diagrams live in a separate package, +[PersistenceDiagrams.jl](https://github.com/mtsch/PersistenceDiagrams.jl). The package is +documented in detail [here](https://mtsch.github.io/PersistenceDiagrams.jl/dev/). + +If you are looking for +Wasserstein or bottleneck distances, persistence images, betti curves, landscapes, and +similar, you will need to run `using PersistenceDiagrams`. + +For convenience, the following basic functionality is reexported by Ripserer: + +```@docs +PersistenceDiagrams.PersistenceDiagram +``` + +```@docs +PersistenceDiagrams.PersistenceInterval +``` + +```@docs +birth(::PersistenceDiagrams.PersistenceInterval) +``` + +```@docs +death(::PersistenceDiagrams.PersistenceInterval) +``` + +```@docs +persistence(::PersistenceDiagrams.PersistenceInterval) +``` + +```@docs +midlife +``` + +```@docs +representative(::PersistenceDiagrams.PersistenceInterval) +``` + +```@docs +birth_simplex(::PersistenceDiagrams.PersistenceInterval) +``` + +```@docs +death_simplex(::PersistenceDiagrams.PersistenceInterval) +``` + +```@docs +barcode +``` + +## Simplices and Representatives + +```@docs +Ripserer.Simplex +``` + +```@docs +Ripserer.Cube +``` + +```@docs +Ripserer.dim(::Ripserer.AbstractCell) +``` + +```@docs +Ripserer.birth(::Ripserer.AbstractCell) +``` + +```@docs +Ripserer.index(::Ripserer.AbstractCell) +``` + +```@docs +Ripserer.vertices(::Ripserer.AbstractCell) +``` + +```@docs +Ripserer.Chain +``` + +```@docs +Mod +``` + +## Experimental Features + +```@docs +Ripserer.reconstruct_cycle +``` + +```@docs +Ripserer.Partition +``` + +```@docs +Ripserer.CircularCoordinates +``` + +## Abstract Types and Interfaces + +```@docs +Ripserer.AbstractFiltration +``` + +```@docs +Ripserer.nv(::Ripserer.AbstractFiltration) +``` + +```@docs +Ripserer.births(::Ripserer.AbstractFiltration) +``` + +```@docs +Ripserer.vertices(::Ripserer.AbstractFiltration) +``` + +```@docs +Ripserer.edges(::Ripserer.AbstractFiltration) +``` + +```@docs +Ripserer.simplex_type +``` + +```@docs +Ripserer.simplex +``` + +```@docs +Ripserer.unsafe_simplex +``` + +```@docs +Ripserer.unsafe_cofacet +``` + +```@docs +Ripserer.threshold(::Ripserer.AbstractFiltration) +``` + +```@docs +Ripserer.columns_to_reduce +``` + +```@docs +Ripserer.emergent_pairs +``` + +```@docs +Ripserer.postprocess_diagram +``` + +```@docs +Ripserer.distance_matrix +``` + +```@docs +Ripserer.AbstractRipsFiltration +``` + +```@docs +Ripserer.adjacency_matrix(::Ripserer.AbstractFiltration) +``` + +```@docs +Ripserer.AbstractCustomFiltration +``` + +```@docs +Ripserer.simplex_dicts +``` + +```@docs +Ripserer.AbstractCell +``` + +```@docs +Ripserer.AbstractSimplex +``` + +```@docs +Base.sign(::Ripserer.AbstractCell) +``` + +```@docs +Base.:-(::Ripserer.AbstractCell) +``` + +```@docs +Ripserer.coboundary +``` + +```@docs +Ripserer.boundary +``` diff --git a/docs/src/api/extensions.md b/docs/src/api/extensions.md deleted file mode 100644 index f96d8e2f..00000000 --- a/docs/src/api/extensions.md +++ /dev/null @@ -1,125 +0,0 @@ -# Interfaces - -## Filtration Interface - -```@docs -Ripserer.AbstractFiltration -``` - -```@docs -Ripserer.nv(::Ripserer.AbstractFiltration) -``` - -```@docs -Ripserer.births(::Ripserer.AbstractFiltration) -``` - -```@docs -Ripserer.vertices(::Ripserer.AbstractFiltration) -``` - -```@docs -Ripserer.edges(::Ripserer.AbstractFiltration) -``` - -```@docs -Ripserer.simplex_type -``` - -```@docs -Ripserer.simplex -``` - -```@docs -Ripserer.unsafe_simplex -``` - -```@docs -Ripserer.unsafe_cofacet -``` - -```@docs -Ripserer.threshold(::Ripserer.AbstractFiltration) -``` - -```@docs -Ripserer.columns_to_reduce -``` - -```@docs -Ripserer.emergent_pairs -``` - -```@docs -Ripserer.postprocess_diagram -``` - -```@docs -Ripserer.distance_matrix -``` - -```@docs -Ripserer.AbstractRipsFiltration -``` - -```@docs -Ripserer.adjacency_matrix(::Ripserer.AbstractFiltration) -``` - -```@docs -Ripserer.AbstractCustomFiltration -``` - -```@docs -Ripserer.simplex_dicts -``` - -## Simplex Interface - -```@docs -Ripserer.AbstractCell -``` - -```@docs -Ripserer.AbstractSimplex -``` - -```@docs -Ripserer.dim(::Ripserer.AbstractCell) -``` - -```@docs -Ripserer.birth(::Ripserer.AbstractCell) -``` - -```@docs -Ripserer.index(::Ripserer.AbstractCell) -``` - -```@docs -Ripserer.vertices(::Ripserer.AbstractCell) -``` - -```@docs -Base.sign(::Ripserer.AbstractCell) -``` - -```@docs -Base.:-(::Ripserer.AbstractCell) -``` - -```@docs -Ripserer.coboundary -``` - -```@docs -Ripserer.boundary -``` - -```@docs -Ripserer.Simplex -``` - -```@docs -Ripserer.Cube -``` diff --git a/docs/src/api/ripserer.md b/docs/src/api/ripserer.md deleted file mode 100644 index 477bed4b..00000000 --- a/docs/src/api/ripserer.md +++ /dev/null @@ -1,83 +0,0 @@ -# Public API - -## Ripserer - -```@docs -ripserer -``` - -## Filtrations - -```@docs -Rips -``` - -```@docs -Cubical -``` - -```@docs -Custom -``` - -```@docs -Alpha -``` - -```@docs -EdgeCollapsedRips -``` - -## Fields - -```@docs -Mod -``` - -## Persistence Diagrams - -See also: [PersistenceDiagrams.jl -API](https://mtsch.github.io/PersistenceDiagrams.jl/dev/api/). For convenience, Ripserer -reexports the following: - -```@docs -birth(::PersistenceDiagrams.PersistenceInterval) -``` - -```@docs -death(::PersistenceDiagrams.PersistenceInterval) -``` - -```@docs -persistence(::PersistenceDiagrams.PersistenceInterval) -``` - -```@docs -representative(::PersistenceDiagrams.PersistenceInterval) -``` - -```@docs -birth_simplex(::PersistenceDiagrams.PersistenceInterval) -``` - -```@docs -death_simplex(::PersistenceDiagrams.PersistenceInterval) -``` - -```@docs -barcode -``` - -## Experimental Features - -```@docs -Ripserer.reconstruct_cycle -``` - -```@docs -Ripserer.Partition -``` - -```@docs -Ripserer.CircularCoordinates -``` diff --git a/docs/src/assets/generate-logo.jl b/docs/src/assets/generate-logo.jl new file mode 100644 index 00000000..43d267f1 --- /dev/null +++ b/docs/src/assets/generate-logo.jl @@ -0,0 +1,72 @@ +# This file generates the plot in the logo +using Ripserer +using Plots +gr() +using Random +Random.seed!(1337) + +function circle(n; r=0.7, center, s=0) + points = NTuple{2,Float64}[] + x, y = center + for θ in range(0, 2π; length=n + 1)[circshift(1:n, s)] + push!(points, (x + r * sin(θ), y + r * cos(θ))) + end + return points +end + +data = [ + circle(20; center=(-1, 0), s=0) + circle(20; center=(1, 0), s=0) + circle(20; center=(0, √3), s=0) +] + +res = ripserer(data; reps=true, dim_max=2) + +ripserer_logo = plot(; + aspect_ratio=1, + xlim=(-2.1, 2.1), + ylim=(-1.1, 3.1), + xlabel=:x, + ylabel=:y, + ticks=nothing, + legend=false, + border=:none, +) +for (i, int) in enumerate(sort(res[2]; by=death, rev=true)) + if i == 1 + col = 2 # red + α = 1 + elseif i == 2 + col = 4 # purple + α = 1 + elseif i == 3 + col = 3 # green + α = 1 + else + col = 1 + α = 0.2 + end + + plot!( + ripserer_logo, + int, + data; + alpha=α, + color=col, + threshold_strict=death(int), + linewidth=5, + ) +end +scatter!(ripserer_logo, data; markersize=4, color=1) +display(ripserer_logo) + +diagrams_logo = plot( + res; + title="", + ticks=nothing, + legend=false, + xlabel="", + ylabel="", + markersize=3, + size=(100, 100), +) diff --git a/docs/src/assets/logo-title.png b/docs/src/assets/logo-title.png new file mode 100644 index 00000000..bae1d088 Binary files /dev/null and b/docs/src/assets/logo-title.png differ diff --git a/docs/src/examples/basics.jl b/docs/src/examples/basics.jl index 892dcd2f..17cf04e2 100644 --- a/docs/src/examples/basics.jl +++ b/docs/src/examples/basics.jl @@ -1,18 +1,19 @@ -# # Basics +# # Usage Guide -# In this example, we will present the usage of Ripserer. We start by loading some packages. +# In this example, we will present the basics of using Ripserer. We start by loading some +# packages. -using Ripserer +using Distances using Plots +using Ripserer using Random # hide Random.seed!(1337) # hide gr() # hide nothing # hide -# ## Basic Usage +# ## Using Ripserer With Point Cloud Data -# Let's start with a basic example, points randomly sampled from a noisy circle. -# We start by defining our sampling function. +# 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}[] @@ -23,54 +24,75 @@ function noisy_circle(n; r=1, noise=0.1) return points end -# Next, we sample 100 points from the circle. - circ_100 = noisy_circle(100) scatter(circ_100; aspect_ratio=1, legend=false, title="Noisy Circle") -# To compute the persistent homology, simply run the following. The `dim_max` argument sets -# the maximum dimension persistent homology is computed in. +# !!! tip "Point-like data types" +# Ripserer can interpret various kinds of data as point clouds. The limitation is that +# the data set should be an `AbstractVector` with elements with the following +# properties: +# * all elements are collections numbers; +# * all elements have the same length. +# Examples of element types that work are `Tuple`s, +# [`SVector`](https://github.com/JuliaArrays/StaticArrays.jl)s, and +# [`Point`](https://github.com/JuliaGeometry/GeometryBasics.jl)s. + +# To compute the Vietoris-Rips persistent homology of this data set, run the +# following. + +ripserer(circ_100) + +# 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) -result_circ = ripserer(circ_100; dim_max=3) +# The result can be plotted as a persistence diagram or as a barcode. -# !!! warning "Warning" -# Computing Vietoris-Rips persistent homology in high dimensions for large numbers of -# points is computationally expensive and requires a large amount of memory. Be careful -# or you **will** run out of memory. On an ordinary laptop, you can expect to compute -# one-dimensional persistent homology for datasets of a few thousand points and higher -# (2-3) dimensional persistent homology for datasets of a few hundred points. This, of -# course, depends on the data set itself. +plot(result_rips) +barcode(result_rips) +plot(plot(result_rips), barcode(result_rips)) # hide -# The result can be plotted as a persistence diagram. +# 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 [`PersistenceDiagram`](@ref)s. The +# zero-dimensional diagram is found at index 1. -plot(result_circ) +plot(result_rips[2]) +barcode(result_rips[2:end]; linewidth=2) +plot(plot(result_rips[2]), barcode(result_rips[2:end]; linewidth=2)) # hide -# Or as a barcode. +# Plotting can be further customized using the standard attributes from +# [Plots.jl](http://docs.juliaplots.org/latest/). -barcode(result_circ) +plot(result_rips; markeralpha=1, markershape=:star, color=[:red, :blue, :green, :purple]) -# ``H_1``, ``H_2``, and ``H_3`` in this plot are hard to see because we have too many -# ``H_0`` bars. We can plot only some of the diagrams. +# ## Changing Filtrations -# !!! note "Note" -# `result` is just an array of persistence diagrams, so the zero-dimensional diagram is -# found at index 1. +# By default, calling [`ripserer`](@ref) will compute persistent homology with the +# [`Rips`](@ref) filtration. To use a different filtration, we have two options. -barcode(result_circ[2:end]; linewidth=2) +# The first option is to pass the filtration constructor as the first argument. Any keyword +# arguments the filtration accepts can be passed to [`ripserer`](@ref) and it will be +# forwarded to the constructor. -# We can plot a single diagram in the same manner. +ripserer(EdgeCollapsedRips, circ_100; threshold=1, dim_max=3, metric=Euclidean()) -barcode(result_circ[3]; linewidth=3) +# The second option is to initialize the filtration object first and use that as an argument +# to [`ripserer`](@ref). 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) # ## Distance Matrix Inputs # In the previous example, we got our result by passing a collection of points to -# `ripserer`. Under the hood, the algorithm actually works with distance matrices. Let's -# define a distance matrix of the shortest paths on a [regular -# icosahedron](https://en.wikipedia.org/wiki/Regular_icosahedron) graph. +# [`ripserer`](@ref). Under the hood, [`Rips`](@ref) and [`EdgeCollapsedRips`](@ref) +# actually work with distance matrices. Let's define a distance matrix of the shortest +# paths on a [regular icosahedron](https://en.wikipedia.org/wiki/Regular_icosahedron) graph. # ```@raw html -# +# # ``` icosahedron = [ @@ -89,15 +111,10 @@ icosahedron = [ ] nothing # hide -# To compute the persistent homology, simply feed the distance matrix to `ripserer`. +# To compute the persistent homology, simply feed the distance matrix to [`ripserer`](@ref). result_icosa = ripserer(icosahedron; dim_max=2) -# Because an icosahedron is topologically equivalent to a sphere, we got a single class in -# the second dimension. - -result_icosa[3] - # ## Thresholding # In our next example, we will show how to use thresholding to speed up computation. We @@ -130,14 +147,14 @@ nothing # hide 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. We can extract this interval by doing the -# following. +# 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 = sort(result_cut[2]; by=persistence)[end] +most_persistent = result_cut[2][end] # 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! +# should not change. Let's try it out. @time result_cut_thresh_2 = ripserer(cutout_2000; threshold=2) nothing # hide @@ -149,6 +166,7 @@ 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. +@assert result_cut_thresh_2 == result_cut # hide result_cut_thresh_2 == result_cut # If we pick a threshold that is too low, we still detect the interval, but its death time @@ -159,4 +177,101 @@ nothing # hide # -plot(result_cut_thresh_1; title="Persistence Diagram, threshold=1") +result_cut_thresh_1[2][end] + +# ## Persistence Diagrams + +# The result of a computation is returned as a vector of +# [`PersistenceDiagram`](@ref)s. Let's take a closer look at one of those. + +diagram = result_cut[2] + +# The diagram is a structure that acts as a vector of +# [`PersistenceInterval`](@ref)s. As such, you can use standard Julia +# functions on the diagram. + +# For example, to extract the last three intervals by birth time, you can do +# something like this. + +sort(diagram; by=birth, rev=true)[1:3] + +# To find the [`persistence`](@ref)s of all the intervals, you can use broadcasting. + +persistence.(diagram) + +# Unlike regular vectors, a [`PersistenceDiagram`](@ref) has additional metadata attached +# to it. To see all metadata, use +# [`propertynames`](https://docs.julialang.org/en/v1/base/base/#Base.propertynames). + +propertynames(diagram) + +# You can access the properties with the dot syntax. + +diagram.field + +# The attributes `dim` and `threshold` are given special treatment and can be extracted +# with appropriately named functions. + +dim(diagram), threshold(diagram) + +# Now, let's take a closer look at one of the intervals. + +interval = diagram[end] + +# An interval is very similar to a tuple of two `Float64`s, but also has some metadata +# associated with it. + +interval[1], interval[2] + +# [`birth`](@ref), [`death`](@ref), [`persistence`](@ref), and [`midlife`](@ref) can be used +# to query commonly used values. + +birth(interval), death(interval), persistence(interval), midlife(interval) + +# Accessing metadata works in a similar manner as with diagrams. + +propertynames(interval) + +# + +interval.birth_simplex + +# + +interval.death_simplex + +# ## Simplices + +# In the previous section, we saw each interval has an associated [`birth_simplex`](@ref) +# and [`death_simplex`](@ref). These values are of the type [`Simplex`](@ref). Let's take a +# closer look at simplices. + +simplex = interval.death_simplex + +# [`Simplex`](@ref) is an internal data structure that uses some tricks to increase +# efficiency. For example, if we were to +# [`dump`](https://docs.julialang.org/en/v1/base/io-network/#Base.dump) it, we notice the +# vertices are not actually stored in the simplex itself. + +dump(simplex) + +# To access the vertices, we use [`vertices`](@ref). + +vertices(simplex) + +# Other useful attributes a simplex has are [`index`](@ref), [`dim`](@ref), and +# [`birth`](@ref). + +index(simplex), dim(simplex), birth(simplex) + +# A few additional notes on simplex properties. + +# * A `D`-dimensional simplex is of type `Simplex{D}` and has `D + 1` vertices. +# * [`vertices`](@ref) are always sorted in descending order. +# * [`index`](@ref) and [`dim`](@ref) can be used to uniquely identify a given simplex. +# * [`birth`](@ref) determines when a simplex is added to a filtration. + +# ## Conclusion + +# This concludes the basic usage of Ripserer. For more detailed information, please check +# out the [API](@ref) page, as well as other examples. diff --git a/docs/src/index.md b/docs/src/index.md index 36ee018e..336af570 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -12,10 +12,9 @@ detects the global topological and local geometric structure of data in a noise- stable way. If you are unfamiliar with persistent homology, I recommend reading this [excellent introduction](https://towardsdatascience.com/persistent-homology-with-examples-1974d4b9c3d0). -See the Examples for further info and usage. -The main goal of this project is to provide an easy to use, generic and fast implementation -of persistent homology to the Julia ecosystem. +Please see the [Usage Guide](@ref) for a quick introduction, and the [API](@ref) page for +detailed descriptions of Ripserer's functionality. While this package is fully functional, it is still in development and should not be considered stable. I try to disrupt the public interface as little as possible, but breaking @@ -39,45 +38,47 @@ version of Julia for optimal performance. Ripserer and its companion package [PersistenceDiagrams.jl](https://github.com/mtsch/PersistenceDiagrams.jl) currently support -* Fast Vietoris-Rips and cubical persistent homology computation. -* Representative cocycle and critical simplex computation. +* Fast Vietoris-Rips and cubical, and alpha complex persistent homology computation. +* Representative cocycle, cycle, and critical simplex computation. * Convenient persistence diagram and representative cocycle visualization via - [Plots.jl](https://github.com/JuliaPlots/Plots.jl). [Experimental Makie.jl - support](https://github.com/mtsch/MakieRipserer.jl) is also available. + [Plots.jl](https://github.com/JuliaPlots/Plots.jl). Experimental + [Makie.jl](https://github.com/JuliaPlots/Makie.jl) is also available + [here](https://github.com/mtsch/MakieRipserer.jl). * Bottleneck and Wasserstein matching and distance computation. * Various persistence diagram vectorization functions, implemented with persistence images and persistence curves. * Easy extensibility through a documented API. -* Computing persistent homology and representative cycles (experimental). +* Experimental shortest representative cycle computation. +* Experimental sparse circular coordinate computation. -To access some of the features, you need to use PersistenceDiagrams. +To access some of the features, you need to use the PersistenceDiagrams.jl package. ## Performance -Much like Ripser, Ripserer uses implicit simplicial complex and reduction matrix -representations combined with the clearing optimization and other computational tricks to -achieve its speed. For a more detailed overview of these optimizations, check out [Ulrich -Bauer's article on Ripser](https://arxiv.org/abs/1908.02518). +Much like Ripser, Ripserer uses several computational tricks to achieve its speed. Among +others, these include an implicit simplicial complex representation and the clearing +optimization. For a more detailed overview of these optimizations, check out [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), usually within around 30%. Cubical homology is -up to 3× slower than that of [Cubical Ripser](https://github.com/CubicalRipser/), which uses -a more specialized algorithm. Ripserer is still a good choice for small 3d images and large -2d images. Ripserer's strength performance-wise is very sparse inputs. It also computes -some things Ripser skips, like the critical simplices. See the [Benchmarks](@ref) section for -detailed benchmarks. +[Ripser](https://github.com/Ripser/ripser), usually within around 30%. Ripserer's strength +performance-wise is very sparse inputs, where it can sometimes outperform Ripser. It also +computes some things Ripser skips, like the critical simplices. -## Extensibility +Ripserer's Cubical homology is up to 3× slower than that of [Cubical +Ripser](https://github.com/CubicalRipser/), which uses a more specialized +algorithm. Ripserer is still a good choice for small 3d images and large 2d images. Unlike +Cubical Ripser, it also supports computations on images of dimensions higher than 4. + +See the [Benchmarks](@ref) section for more detailed benchmarks. + +## Extending Ripserer is designed to be easily extended with new simplex or filtration types. See the -[Filtration Interface](@ref) and [Simplex Interface](@ref) API sections for more info. To -see an example of an extension, check out the implementation of cubical homology in -[`src/cubical.jl`](https://github.com/mtsch/Ripserer.jl/blob/master/src/cubical.j). Keep in -mind that this extension customizes almost every aspect of the algorithm and that extensions -can usually be much simpler. +[Abstract Types and Interfaces](@ref) API section for more information. If you have written an extension or are having trouble implementing one, please feel free to -open a pull request or an issue, or contact me directly. +open a pull request or an issue. You may also contact me directly. ## Contributing diff --git a/src/Ripserer.jl b/src/Ripserer.jl index 76793fce..73fced1c 100644 --- a/src/Ripserer.jl +++ b/src/Ripserer.jl @@ -31,7 +31,8 @@ import LightGraphs: vertices, edges, nv, adjacency_matrix # Reexporting basics makes Ripserer usable without having to import another package. import PersistenceDiagrams: birth, threshold, dim -export birth, death, persistence, representative, birth_simplex, death_simplex, barcode +export birth, + death, persistence, midlife, representative, birth_simplex, death_simplex, barcode export Mod, Simplex, diff --git a/src/base/chain.jl b/src/base/chain.jl index d2c82f01..c3eebaf8 100644 --- a/src/base/chain.jl +++ b/src/base/chain.jl @@ -1,13 +1,63 @@ +# A chain should behave exactly like an array of `ChainElements`. Internally, it may store +# its elements in a more efficient way. + +# It can be used as a heap with `heapify!`, `heappush!`, and `heapop!`. When used this way, +# it will only return unique elements, summing elements with the same simplex together. + +# `clean!` can be used to sum duplicates together and sort the chain. """ Chain{F,S,E} <: AbstractVector{ChainElement{S,F}} -A chain should behave exactly like an array of `ChainElements`. Internally, it may store its -elements in a more efficient way. +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`](@ref)). + +Most functions that can be called on `AbstractCell`s can also be called on the elements of a +`Chain`. + +# Examples + +```jldoctest; setup=(using Random; Random.seed!(1337)) +julia> data = [(rand(), rand(), rand()) for _ in 1:100]; + +julia> result = ripserer(data; reps=true, modulus=7); + +julia> chain = result[end][end].representative +6-element Chain{Ripserer.Mod{7},Ripserer.Simplex{1,Float64,Int64}}: + +Simplex{1}((87, 59), 0.23148225999797645) => 6 mod 7 + +Simplex{1}((59, 46), 0.3054281021426286) => 1 mod 7 + +Simplex{1}((87, 14), 0.32453294355760326) => 6 mod 7 + +Simplex{1}((87, 1), 0.34642558062390577) => 6 mod 7 + +Simplex{1}((100, 87), 0.3480194268484163) => 1 mod 7 + +Simplex{1}((79, 1), 0.36519064466525686) => 6 mod 7 + +julia> simplex.(chain) +6-element Array{Simplex{1,Float64,Int64},1}: + +Simplex{1}((87, 59), 0.23148225999797645) + +Simplex{1}((59, 46), 0.3054281021426286) + +Simplex{1}((87, 14), 0.32453294355760326) + +Simplex{1}((87, 1), 0.34642558062390577) + +Simplex{1}((100, 87), 0.3480194268484163) + +Simplex{1}((79, 1), 0.36519064466525686) + +julia> vertices.(chain) +6-element Array{Tuple{Int64,Int64},1}: + (87, 59) + (59, 46) + (87, 14) + (87, 1) + (100, 87) + (79, 1) -It can be used as a heap with `heapify!`, `heappush!`, and `heapop!`. When used this way, it -will only return unique elements, summing elements with the same simplex together. +julia> coefficient.(chain) +6-element Array{Mod{7},1}: + 6 mod 7 + 1 mod 7 + 6 mod 7 + 6 mod 7 + 1 mod 7 + 6 mod 7 -`clean!` can be used to sum duplicates together and sort the chain. +``` """ struct Chain{F,S,E} <: AbstractVector{ChainElement{S,F}} elements::Vector{E} diff --git a/src/base/simplex.jl b/src/base/simplex.jl index cecd7ba7..301a7e6e 100644 --- a/src/base/simplex.jl +++ b/src/base/simplex.jl @@ -13,17 +13,18 @@ time of type `T`. # Examples ```jldoctest -julia> Simplex{2}(2, 1) +julia> sx = Simplex{2}(2, 1) 2-dimensional Simplex(index=2, birth=1): +(4, 2, 1) -julia> Simplex{10}(Int128(-10), 1.0) -10-dimensional Simplex(index=10, birth=1.0): - -(12, 11, 10, 9, 8, 7, 6, 5, 4, 2, 1) +julia> index(sx) +2 -julia> Simplex{2}((5, 2, 1), 1) -2-dimensional Simplex(index=5, birth=1): - +(5, 2, 1) +julia> vertices(sx) +(4, 2, 1) + +julia> birth(sx) +1 ``` """ diff --git a/src/computation/ripserer.jl b/src/computation/ripserer.jl index 2a977a86..78d0cbeb 100644 --- a/src/computation/ripserer.jl +++ b/src/computation/ripserer.jl @@ -42,6 +42,10 @@ Compute the persistent homology of a filtration. The filtration can be given as [`AbstractFiltration`](@ref) 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 [`PersistenceDiagram`](@ref)s 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`. @@ -52,19 +56,19 @@ Compute the persistent homology of a filtration. The filtration can be given as * `field`: use this type of field of coefficients. Defaults to [`Ripserer.Mod`](@ref)`{modulus}`. -* `threshold`: compute persistent homology up to diameter smaller than threshold. This - parameter is only applicable when using distance matrices or points as input. When using - filtrations, threshold must be passed to the filtration constructor. Defaults to the - radius of the input space. When using low thresholds with points or distance matrices, - consider using `sparse=true`. +* `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`. +* `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. + be very slow. Defaults to `false` for cohomology and `1:dim_max` for homology. + Representatives are wrapped in a [`Chain`](@ref). * `verbose`: If `true`, show a verbose bar. Defaults to `false`.