## To install the pre-compiled package TDA for PC and possibly linux.  
## Note a pre-compiled version is not available for the Mac.

install.packages("TDA")       

## To install and compile the package TDA.  
## Note you need to have compilers installed on your computer.

install.packages("TDA", type = "source")  

## load package TDA

library("TDA")

## We will first generate a data set where we know the topology
## circleUnif(30, r = 1) chooses 30 points randomly (with uniform distribution)
## from a circle of radius 1.
## A set of data points is called a point cloud.
## Thus XX is a point cloud.

XX = circleUnif(30, r = 1)
plot (XX, asp = 1) # plot XX with aspect ratio = 1

?ripsDiag
## ripsDiag computes the persistence diagram of the Rips filtration built on top of a point cloud.
## First argument = data set (matrix of data points)
## Second argument = max dimension of the homological features to be computed
## Third argument =  maximum value of the rips filtration (largest threshold).
## The output is a list containing the oobject of class diagram, a P by 3 matrix, 
## where P is the number of cycles in the barcode and 
## the columns correspond to dimension, birth, death of the cycles.
## To access the diagram in the list, use $diagram

maxscale=2.5  # max radius = 2.5
maxdimension=2  # to compute H0, H1, H2
Diag=ripsDiag(XX,maxdimension,maxscale, printProgress=TRUE)$diagram

## Plot barcode

plot(Diag, diagLim=NULL, dimension=NULL, col=NULL, rotated=FALSE, barcode=TRUE, band=NULL, lab.line=2.2, colorBand="pink", colorBorder=NA, add = FALSE)

## Plot persistence diagram

plot(Diag, diagLim=NULL, dimension=NULL, col=NULL, rotated=FALSE, barcode=FALSE, band=NULL, lab.line=2.2, colorBand="pink", colorBorder=NA, add = FALSE)

Diag ##  View birth and death times of each cycle 

## Let's now add noise to our data set:

?runif
noise=cbind(runif(20, -2,2),runif(20, -2,2))
noise

## Add 20 more data points from noise to 30 point data set XX
## for a total of 50 data points in XXnoisy 

XXnoisy=rbind(XX, noise)
plot(XXnoisy, asp = 1)

# Caculate ripsDiag

DiagN=ripsDiag(XXnoisy,maxdimension,maxscale, printProgress=TRUE)$diagram

## Plot barcode

plot(DiagN, diagLim=NULL, dimension=NULL, col=NULL, rotated=FALSE, barcode=TRUE, band=NULL, lab.line=2.2, colorBand="pink", colorBorder=NA, add = FALSE)

## Plot persistence diagram

plot(DiagN, diagLim=NULL, dimension=NULL, col=NULL, rotated=FALSE, barcode=FALSE, band=NULL, lab.line=2.2, colorBand="pink", colorBorder=NA, add = FALSE)

?bottleneck
print( bottleneck(Diag, DiagN, dimension=0) )
? wasserstein
print( wasserstein(Diag, DiagN, p=2, dimension=0) )

print( bottleneck(Diag, DiagN, dimension=1) )
print( wasserstein(Diag, DiagN, p=2, dimension=1) )

print( bottleneck(Diag, DiagN, dimension=2) )
print( wasserstein(Diag, DiagN, p=2, dimension=2) )

##  Add even more noise

noise=cbind(runif(40, -2,2),runif(40, -2,2))

## Add 40 more data points from noise to 30 point data set XX
## for a total of 70 data points in XXnoisy 

XXnoisy2=rbind(XX, noise)
plot(XXnoisy2, asp = 1)

# Caculate ripsDiag

DiagNN=ripsDiag(XXnoisy2,maxdimension,maxscale, printProgress=TRUE)$diagram

## Plot barcode

plot(DiagNN, diagLim=NULL, dimension=NULL, col=NULL, rotated=FALSE, barcode=TRUE, band=NULL, lab.line=2.2, colorBand="pink", colorBorder=NA, add = FALSE)

## Plot persistence diagram

plot(DiagNN, diagLim=NULL, dimension=NULL, col=NULL, rotated=FALSE, barcode=FALSE, band=NULL, lab.line=2.2, colorBand="pink", colorBorder=NA, add = FALSE)

?bottleneck
print( bottleneck(Diag, DiagNN, dimension=0) )
? wasserstein
print( wasserstein(Diag, DiagNN, p=2, dimension=0) )

print( bottleneck(Diag, DiagNN, dimension=1) )
print( wasserstein(Diag, DiagNN, p=2, dimension=1) )

print( bottleneck(Diag, DiagNN, dimension=2) )
print( wasserstein(Diag, DiagNN, p=2, dimension=2) )