Note that a more detailed tutorial on the R TDA package is available here.

Before we begin, load the TDA package. (If the package is not installed, you can install it using RStudio, or type install.packages("TDA").)

library("TDA")

Vietoris-Rips Persistence

The ripsDiag function computes the persistence diagram of a Vietoris-Rips complex. It returns a list containing a diagram object, which can be plotted either as a persistence diagram or as a barcode. The ripsDiag function can use any of the libraries GUDHI, Dionysus, or PHAT.

points1 <- read.csv("~/Dropbox/data/point cloud scripts/points1.csv")
plot(points1, xlab="", ylab="", main="points1")

Diag1 <- ripsDiag(X = points1, maxdimension = 1, maxscale = 4, library="GUDHI", location = FALSE, printProgress = FALSE)
plot(Diag1[["diagram"]], main="persistence diagram from points1")

plot(Diag1[["diagram"]], barcode=TRUE, main="barcode from points1")

Superlevelset Persistence

We first compute a density estimator for the point cloud. We will use the Gaussian Kernel Density Estimator (KDE), computed over a grid of points in the plane. We first set up the grid, then compute the density estimator using the kde function provided by the TDA package. We can also plot the density estimator using the perspective plot function persp. Note that the persp function requires many parameters to set the viewing direction and shading.

Xlim <- c(-3,3)
Ylim <- c(-3,3)
by <- 0.1
Xseq <-seq(Xlim[1], Xlim[2], by=by)
Yseq <-seq(Ylim[1], Ylim[2], by=by)
Grid <- expand.grid(Xseq, Yseq)
points1kde <- kde(X=points1, Grid = Grid, h = 0.3)
persp(Xseq, Yseq, matrix(points1kde, ncol = length(Yseq), nrow = length(Xseq)), xlab = "",ylab = "", zlab = "", theta = -20, phi = 35, ltheta = 50, col = 2, d = 0.5, main="KDE points1")

We now compute persistent homology of the superlevel sets of the density estimator using the gridDiag function. The gridDiag function can use either of the libraries Dionysus or PHAT.

Diag1kde <- gridDiag(X = points1, FUN = kde, h = 0.3, lim = cbind(Xlim, Ylim), by = by, sublevel = FALSE, library = "PHAT", printProgress = TRUE)
## # Generated complex of size: 21841 
## # Persistence timer: Elapsed time [ 0.016496 ] seconds
plot(Diag1kde[["diagram"]], main="KDE diagram")

Bottleneck Distance

We now measure the similarity between two persistence diagrams by computing the bottleneck distance. The bottleneck distance is the shortest distance b for which there exists a perfect matching between points of the two diagrams (points may be matched to the diagonal) such that the distance between any pair of matched points is at most b.

First, we generate a persistence diagram from points (exactly) on a circle of radius 2. The TDA package provides a function circleUnif to sample points on a circle.

circle <- circleUnif(200, 2)
plot(circle)

Diag2 <- ripsDiag(circle, 1, 4)
plot(Diag2[["diagram"]])

Notice that the diagram above is similar to the first persistence diagram that we computed (Diag1), but with less noise near the diagonal.

We now compute the bottleneck distance between the Diag1 and Diag2, using the bottleneck command provided by the TDA package.

print(bottleneck(Diag1[["diagram"]], Diag2[["diagram"]]))
## [1] 0.9459715

The bottleneck command also allows us to compute the bottleneck distance using only the features of a particular dimension.

print(bottleneck(Diag1[["diagram"]], Diag2[["diagram"]], dimension=0))
## [1] 0.2720294
print(bottleneck(Diag1[["diagram"]], Diag2[["diagram"]], dimension=1))
## [1] 0.9459715

Your Turn

Now experiment with the data sets available on the conference website, or work through the R-TDA package tutorial. You might also want to refer to the paper Introduction to the R package TDA.