computePI: A Vector Summary of the Persistence Surface

In TDAvec: Vector Summaries of Persistence Diagrams

A Vector Summary of the Persistence Surface

Description

For a given persistence diagram D=\{(b_i,p_i)\}_{i=1}^N, computePI() computes the persistence image (PI) - a vector summary of the persistence surface:

ρ(x,y)=∑_{i=1}^N f(b_i,p_i)φ_{(b_i,p_i)}(x,y),

where φ_{(b_i,p_i)}(x,y) is the Gaussian distribution with mean (b_i,p_i) and covariance matrix σ^2 I_{2\times 2} and

f(b,p) = w(p)=≤ft\{ \begin{array}{ll} 0 & \quad y≤q 0 \\ p/p_{max} & \quad 0<p<p_{max}\\ 1& \quad y≥q p_{max} \end{array} \right.

is the weighting function with p_{max} being the maximum persistence value among all persistence diagrams considered in the experiment. Points of D with infinite persistence value are ignored

Usage

computePI(D, homDim, xSeq, ySeq, sigma)

Arguments

D	matrix with three columns containing the dimension, birth and persistence values respectively
homDim	homological dimension (0 for H_0, 1 for H_1, etc.)
xSeq	numeric vector of increasing x (birth) values used for vectorization
ySeq	numeric vector of increasing y (persistence) values used for vectorization
sigma	standard deviation of the Gaussian

Value

A numeric vector whose elements are the average values of the persistence surface computed over each cell of the two-dimensional grid constructred from xSeq=\{x_1,x_2,…,x_n\} and ySeq=\{y_1,y_2,…,y_m\}:

\Big(\frac{1}{Δ x_1Δ y_1}\int_{[x_1,x_2]\times [y_1,y_2]}ρ(x,y)dA,…,\frac{1}{Δ x_{n-1}Δ y_{m-1}}\int_{[x_{n-1},x_n]\times [y_{m-1},y_m]}ρ(x,y)dA\Big),

where dA=dxdy, Δ x_k=x_{k+1}-x_k and Δ y_j=y_{j+1}-y_j

Author(s)

Umar Islambekov

References

1. Adams, H., Emerson, T., Kirby, M., Neville, R., Peterson, C., Shipman, P., ... & Ziegelmeier, L. (2017). Persistence images: A stable vector representation of persistent homology. Journal of Machine Learning Research, 18.

Examples

N <- 100 
set.seed(123)
# sample N points uniformly from unit circle and add Gaussian noise
X <- TDA::circleUnif(N,r=1) + rnorm(2*N,mean = 0,sd = 0.2)

# compute a persistence diagram using the Rips filtration built on top of X
D <- TDA::ripsDiag(X,maxdimension = 1,maxscale = 2)$diagram 

# switch from the birth-death to the birth-persistence coordinates
D[,3] <- D[,3] - D[,2] 
colnames(D)[3] <- "Persistence"

resB <- 5 # resolution (or grid size) along the birth axis
resP <- 5 # resolution (or grid size) along the persistence axis 

# compute PI for homological dimension H_0
minPH0 <- min(D[D[,1]==0,3]); maxPH0 <- max(D[D[,1]==0,3])
ySeqH0 <- seq(minPH0,maxPH0,length.out=resP+1)
sigma <- 0.5*(maxPH0-minPH0)/resP 
computePI(D,homDim=0,xSeq=NA,ySeqH0,sigma) 

# compute PI for homological dimension H_1
minBH1 <- min(D[D[,1]==1,2]); maxBH1 <- max(D[D[,1]==1,2])
minPH1 <- min(D[D[,1]==1,3]); maxPH1 <- max(D[D[,1]==1,3])
xSeqH1 <- seq(minBH1,maxBH1,length.out=resB+1)
ySeqH1 <- seq(minPH1,maxPH1,length.out=resP+1)
sigma <- 0.5*(maxPH1-minPH1)/resP
computePI(D,homDim=1,xSeqH1,ySeqH1,sigma) 

TDAvec documentation built on Oct. 31, 2022, 5:06 p.m.