---
title: "Calculating Persistent Homology with a Vietoris-Rips Complex"
author: "Raoul R. Wadhwa, Matthew Piekenbrock, Jacob G. Scott"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Calculating Persistent Homology with a Vietoris-Rips Complex}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(ripserr)
```

## Sample dataset

In this vignette, we will generate a point cloud using a sample of 2-dimensional points on the unit circle's circumference; this will be stored in a variable named `circle_df`.

```{r load-data}
# create reproducible dataset
set.seed(42)
angles <- runif(25, 0, 2 * pi)
circle_df <- data.frame(x = cos(angles),
                        y = sin(angles))

# take a peek at first 6 rows
head(circle_df)
```

Above, each of the 100 rows represents a single point, with each of the 2 columns representing a Cartesian coordinate for a single dimension.
Column `x` contains the x-coordinates of the 100 points and column `y` contains the respective y-coordinates.
To confirm that the points in `circle_df` do lie on the circumference of a circle, we can quickly create a scatterplot.

```{r plot-circle, fig.width = 4, fig.height = 4.5}
# scatterplot of circle2d
plot(circle_df, xlab = "x", ylab = "y", main = "2-d circle point cloud")
```

## Calculating persistent homology

Given that the points in `circle_df` are uniformly distributed across the circumference of a circle without any error or noise, we expect a single prominent 1-cycle to be present in its persistent homology.
The [Ripser](https://github.com/Ripser/ripser) C++ library is wrapped by R using [Rcpp](https://github.com/RcppCore/Rcpp), and performs calculations on a Vietoris-Rips complex created with the input point cloud.
These calculations result in a data frame that contains all the necessary information to characterize the persistence of homological features within `circle_df`, and can be performed with a single line of R code using ripserr.

```{r calc-hom}
# calculate persistent homology
circ_phom <- vietoris_rips(circle_df)

# print first 6 features (ordered by dimension and birth)
head(circ_phom)

# print last 6 features (ordered by dimension and birth)
tail(circ_phom)
```

Each row in the homology matrix returned by the `vietoris_rips` function (variable named `circ_phom`) represents a single feature (cycle).
The homology matrix has 3 columns in the following order:

1. **dimension:** if 0, represents a 0-cycle; if 1, represents a 1-cycle; and so on.
1. **birth:** the radius of the Vietoris-Rips complex at which this feature was first detected.
1. **death:** the radius of the Vietoris-Rips complex at which this feature was last detected.

Persistence of a feature is generally defined as the length of the interval of the radius within which the feature exists.
This can be calculated as the numerical difference between the second (birth) and third (death) columns of the homology matrix.
Confirmed in the output of the `head` and `tail` functions above, the homology matrix is ordered by dimension, with the birth column used to compare features of the same dimension.
As expected for `circle_df`, the homology matrix contains a single prominent 1-cycle (last line of `tail`'s output).
Although we suspect the feature to be a persistent 1-cycle, comparison with the other features in the homology matrix is required to confirm that it is sufficiently persistent.
This task is done far more easily with an appropriate visualization (e.g. topological barcode or persistence diagram) than by eyeballing the contents of `circ_phom`.
The [ggtda](https://github.com/tdaverse/ggtda) and [TDAstats](https://CRAN.R-project.org/package=TDAstats) R packages could be useful to create such visualizations.

## Dataset formats

Although we used a data frame to store the example point cloud, a matrix or dist object would work equally well.
Let's test this out by creating the matrix and dist equivalents of `circle_df` and checking if the resulting `vietoris_rips` outputs are equal using `base::all.equal` to compare data frames.

```{r formats}
# matrix format
circ_mat <- as.matrix(circle_df)
head(circ_mat)

# dist object
circ_dist <- dist(circ_mat)
head(circ_dist)

# calculate persistent homology of each
mat_phom <- vietoris_rips(circ_mat)
dist_phom<- vietoris_rips(circ_dist)

# compare equality
all.equal(circ_phom, mat_phom)
all.equal(circ_phom, dist_phom)
```

## Coefficients in various prime fields

Ripser can calculate persistent homology with coefficients in various prime fields ($\mathbb{Z}/p\mathbb{Z}$).
ripserr implements this functionality via the `p` parameter in `vietoris_rips`.

```{r prime-fields}
# prime field Z/2Z (default)
circ_phom2 <- vietoris_rips(circle_df)

# prime field Z/5Z
circ_phom5 <- vietoris_rips(circle_df, p = 5L)

# confirm equal outputs
all.equal(circ_phom2, circ_phom5)
```
