Lecture 12: Computing Rips Complex Persistence for a Point Cloud#
Installing packages#
This tutorial will be based on the following packages.
Also note that some figures in this notbook are done in plotly. However, these figures do not render on the website. Be sure to download the notebook to see them.
Uncomment and run the commands below if these packages are not already installed in your system.
# !pip install teaspoon
# !pip install scikit-tda
import ripser # This is in the scikit-tda package
from persim import plot_diagrams # Also in scikit-tda
import teaspoon.MakeData.PointCloud as gPC
import teaspoon.TDA.Draw as Draw
# Standard imports
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy.spatial import distance_matrix
import plotly.express as px
import plotly.graph_objects as go
%matplotlib inline
Next, we’re going to compute persistence for a point cloud using the Rips filtration.
This section requires ripser (from scikit-tda) and teaspoon installed. Ripser is a package that is very optimized for point cloud Rips persistence.
Rips Definition#
The setup for persistence in this section is as follows.
We are given a point cloud \(P = \{x_0, \ldots,x_{N-1}\} \subseteq \mathbb{R}^d\)
The Vietoris-Rips complex is defined by \(VR(P, r) = \{ \sigma \subseteq P \mid \|x_i-x_j\| \leq r \; \forall x_i,x_j \in \sigma\} \)
Note that I’m not using the \(2r\) from the book since this is what usually shows up in code.
Persistence module: \(H_p(VR(P,0)) \to H_p(VR(P,r_1)) \to \cdots H_p(VR(P,r_N))\)
Goal: Compute persistence diagram of this module and draw persistence diagram
Tiny example#
Let’s compute the Rips complex for a tiny example: a point cloud with just 6 points. I’ve got them labeled with their order in the matrix \(P\).
P = np.array([[0,0],[1,0],[0,2],[1,1.5],[0.5,0.5],[2,2], [2,0.5], [2.5,1]])
P
array([[0. , 0. ],
[1. , 0. ],
[0. , 2. ],
[1. , 1.5],
[0.5, 0.5],
[2. , 2. ],
[2. , 0.5],
[2.5, 1. ]])
fig, ax = plt.subplots(figsize=(6,6))
ax.scatter(P[:,0], P[:,1], color='tab:purple', edgecolors = 'k', s=80, zorder=3)
labels = [f"{i}" for i in range(len(P))]
for i, txt in enumerate(labels):
ax.text(P[i,0], P[i,1]+0.03, txt, ha='left', va='bottom', fontsize=9, zorder=4)
ax.set_aspect('equal', 'box')
ax.grid(True, linestyle='--', linewidth=0.5, color='lightgray', alpha=0.9)
plt.show()
Remember, for any fixed radius \(r\), the Rips complex has all edges for pairs of points \(x\) and \(y\) that are within distance \(r\) of each other; and then all higher dimensional simplicies that are available.
D = distance_matrix(P, P)
plt.matshow(D)
plt.colorbar();
- What triangles are in VR(P,1.3)? Tetrahedra?
- What triangles are in VR(P,1.6)? Tetrahedra?
r = 1.3
fig, ax = plt.subplots(figsize=(6,6))
for i in range(len(P)):
for j in range(i+1, len(P)):
if D[i, j] < r:
ax.plot([P[i,0], P[j,0]], [P[i,1], P[j,1]], color='gray', linewidth=2)
ax.scatter(P[:,0], P[:,1], color='tab:purple',
edgecolors='k', s=80, zorder=3)
labels = [f"{i}" for i in range(len(P))]
for i, txt in enumerate(labels):
ax.text(P[i,0], P[i,1]+0.03, txt, ha='center', va='bottom', fontsize=9, zorder=4)
ax.set_title(f'Points and edges with distance < {r}')
ax.set_aspect('equal', 'box')
ax.grid(True, linestyle='--', linewidth=0.5, color='lightgray', alpha=0.9)
plt.show()
Answer commented out in this cell.
The ripser package., which you installed with the scikit-tda package, computes Rips persistence for point clouds. To do this, we need only run the following command and pass in the point cloud P.
The ripser package computes persistence of the point cloud for the Rips (aka Vietoris-Rips) complex. Again, see the
output = ripser.ripser(P)
print(f"Ripser's output is of type {type(output)}\nwith keys {output.keys()}")
Ripser's output is of type <class 'dict'>
with keys dict_keys(['dgms', 'cocycles', 'num_edges', 'dperm2all', 'idx_perm', 'r_cover'])
We care about the dgms portion in this case. It is a list of arrays, the ith entry is the i-dimensional persistence diagram.
output['dgms']
[array([[0. , 0.70710677],
[0. , 0.70710677],
[0. , 0.70710677],
[0. , 1.11803401],
[0. , 1.11803401],
[0. , 1.11803401],
[0. , 1.11803401],
[0. , inf]]),
array([[1.41421354, 1.5 ],
[1.11803401, 1.5 ]])]
print('0-dimensional diagram points:')
print(output['dgms'][0])
0-dimensional diagram points:
[[0. 0.70710677]
[0. 0.70710677]
[0. 0.70710677]
[0. 1.11803401]
[0. 1.11803401]
[0. 1.11803401]
[0. 1.11803401]
[0. inf]]
print('1-dimensional diagram points:')
print(output['dgms'][1])
1-dimensional diagram points:
[[1.41421354 1.5 ]
[1.11803401 1.5 ]]
Then we can also use scikit-tda’s plotting to visualize this easily.
plot_diagrams(output['dgms'])
ax = plt.gca()
ax.grid(True, linestyle='--', linewidth=1.2, color='lightgray', alpha=1)
Rips for point clouds sampled from spaces#
Annulus#
Let’s do a more interesting example by computing the persistent homology of a point cloud from an annulus.
Image credit: Wikipedia
Answer commented out in this cell.
The teaspoon has code for creating random point clouds for testing purposes. In particular, we have the Annulus command which generates a point cloud from an annulus. See the documentation for specific usage.
We first generate a point cloud with 200 points from an annulus which has inside radius \(r=2\) and outside radius \(R=3\).
P = gPC.Annulus(N=200, r = 2, R = 3)
Draw.drawPtCloud(P)
ax = plt.gca()
ax.axis('square')
ax.grid(True, linestyle='--', linewidth=0.5, color='lightgray', alpha=0.9)
P
array([[ 1.76243294, 0.97205737],
[-1.26672625, 2.27607755],
[-2.49430639, -1.2800484 ],
[ 2.83770111, 0.09580523],
[ 2.8644357 , -0.38387615],
[ 2.03038535, -1.36911699],
[ 0.49754451, 2.58180015],
[ 0.38025519, -2.94497835],
[ 1.15954566, 2.52816422],
[-2.04693913, 0.78140259],
[ 2.22777347, -0.47923324],
[ 2.3500861 , -1.58805883],
[ 2.05806357, -1.84835211],
[-2.04777456, -1.30399017],
[-2.15252848, -2.06196926],
[-1.04559809, 2.71517326],
[-1.21111187, -2.39947752],
[-2.019565 , -1.68303387],
[-1.56022083, 1.50040389],
[ 1.07717745, 2.4392088 ],
[-0.21604177, -2.40344244],
[-0.47848901, -2.66142265],
[-0.02776209, 2.6492904 ],
[-2.90595973, 0.3210419 ],
[-1.46636893, -1.66905686],
[ 0.90625111, -1.86529861],
[-2.49340269, 0.28944743],
[ 2.84324216, -0.59631852],
[ 1.48364243, -1.69611387],
[ 0.04626871, 2.90373007],
[-2.4766083 , 1.02454059],
[-0.38980522, 2.24549147],
[ 0.00402259, 2.66762951],
[ 0.17186335, -2.15167745],
[-2.66230308, 0.22585939],
[-0.47781738, -2.11175657],
[ 1.41093249, -1.53042717],
[ 1.01179998, -2.05682477],
[ 1.42645998, 2.33057493],
[ 2.29014624, 0.3152334 ],
[ 2.5776112 , 0.36558618],
[ 1.4686417 , 2.56228504],
[ 0.24250334, 2.64648712],
[-2.82779065, -0.45644002],
[ 0.28064451, 2.33564368],
[ 1.10257405, -2.54001123],
[-1.40149213, 1.92612691],
[ 2.63405962, 0.56504005],
[ 2.49247177, 1.27475393],
[-1.77978761, -1.65593865],
[-2.26959587, -1.33390187],
[-0.24062323, 2.51630889],
[-0.32997001, -2.10442957],
[-0.84953516, 1.93031812],
[ 2.82574208, -0.26507861],
[-2.64782397, 1.3907923 ],
[ 0.26333632, -2.33336768],
[ 0.67432154, 2.62716382],
[-1.3457415 , 2.11482066],
[ 1.43857738, -1.87998323],
[-2.86631595, -0.40006947],
[-1.5003611 , 2.52356556],
[ 2.90734303, -0.48698858],
[ 1.68941934, 1.65186164],
[-1.54174209, -2.15688489],
[-1.41709997, -2.42830245],
[ 1.99317487, 1.98483806],
[-2.66731067, -0.07399672],
[ 1.91541029, -1.35987517],
[ 2.50392383, -0.38185748],
[-0.88954589, -2.09392074],
[-2.32586324, 1.80816576],
[-1.94724498, -1.68415527],
[-2.84503877, 0.2552908 ],
[-1.79423258, -2.40310682],
[ 0.61664103, -2.10470217],
[-1.44291411, 1.93911541],
[-0.87326655, -2.5690005 ],
[ 2.85553963, 0.16195922],
[-2.61946997, -0.47751869],
[-1.4358022 , 1.96515746],
[-1.20460551, 1.87762869],
[-2.18348589, -0.83296522],
[-1.19533713, -1.6618921 ],
[ 1.88798519, -0.88440115],
[-2.72743752, 0.26046041],
[ 0.23400844, 2.39902814],
[-2.18041528, -1.43837798],
[-1.63630215, -2.39523752],
[-0.83112571, -2.43628816],
[ 0.95895232, 2.83049933],
[ 1.97590378, -1.37541697],
[-2.33015657, -1.05825042],
[-0.20774524, 2.42530855],
[-2.02718163, -1.12793753],
[-0.05857569, 2.21441586],
[ 2.16633875, -1.33391594],
[-1.94888242, 1.67517294],
[ 2.07111946, 0.932817 ],
[ 2.4371451 , -0.35106948],
[-2.80329537, -0.76442792],
[ 1.65465296, 1.46682304],
[-1.70520708, 2.41620266],
[ 2.80476886, 0.15360835],
[ 2.12875352, 0.44796196],
[ 0.96486068, 2.14147839],
[-2.36238195, -1.51427336],
[-0.71645192, -2.87532347],
[ 0.0158555 , 2.57710989],
[ 0.97630578, -2.60382147],
[-2.62569779, -0.43108784],
[-2.13419516, 1.43410999],
[-2.27350559, -0.26897064],
[-2.16441397, 0.5993368 ],
[-0.10822852, 2.3328303 ],
[ 2.7037723 , 0.25824985],
[ 1.5747935 , -2.3235471 ],
[ 0.21105743, -2.64143135],
[ 2.06815706, -0.18587762],
[ 2.34829223, 0.06454898],
[ 0.30283948, -2.93627354],
[-0.54181773, 2.00818029],
[-0.7529509 , -2.40343136],
[-2.27147549, 0.02276391],
[-0.71165603, -1.91132106],
[ 1.11334786, -2.60200712],
[-1.73247483, 2.19956029],
[-2.26551455, -1.72560151],
[ 0.62891638, 2.57056352],
[-1.11792625, 1.7065487 ],
[-2.50377833, -1.58485172],
[-2.27641468, -1.65666021],
[ 0.89925124, 2.08224407],
[ 0.41555918, 2.88312545],
[-1.37190278, 2.65186908],
[ 2.02073919, -0.79184719],
[-0.08950867, -2.69726008],
[ 0.12077471, -2.18466736],
[-2.38395535, 0.88890668],
[-1.05014165, 2.06255931],
[-2.17876932, -0.30665543],
[ 0.34025338, 2.32935051],
[ 1.95847036, -0.72964278],
[-1.03474818, -1.8522477 ],
[ 2.32881473, -1.73204748],
[-2.48555012, 1.28906456],
[-2.14248618, 0.12928078],
[ 1.69529413, -1.90684823],
[ 0.63826844, 2.7701254 ],
[-2.45792227, -1.65797642],
[ 2.49550758, 1.04625264],
[-0.08812661, -2.03297074],
[-0.59062634, 2.238792 ],
[ 1.93137007, -2.16907727],
[-1.59772231, -2.19894248],
[-2.18214411, -0.11614279],
[-1.2096695 , 2.65871866],
[-1.74232942, -2.41589824],
[-2.84789412, 0.7343372 ],
[-0.85278805, -1.81225075],
[ 2.35148865, -0.80933228],
[ 2.69364347, -0.99492707],
[-2.40856622, -0.7598865 ],
[-2.95438461, -0.22558623],
[-2.97472376, -0.33144545],
[-0.7561349 , -1.85876501],
[ 0.73413466, -2.00195865],
[ 0.33191277, -2.38660804],
[-1.28798254, 2.44645597],
[-0.28334635, 2.08174109],
[-2.40363714, 1.43536019],
[ 2.702349 , 0.27246188],
[ 0.76580191, 1.87271538],
[-0.56342167, 1.99228791],
[-2.77065137, -0.33580568],
[ 1.00944695, 1.91093312],
[-1.15419965, -2.68149142],
[-2.23569652, -0.14361017],
[-2.20793834, -0.75465916],
[-2.43470183, -0.7332519 ],
[ 2.11520409, -0.4686775 ],
[ 1.94898385, -0.8118823 ],
[-1.39482747, 2.27775507],
[ 1.05736827, -1.93849771],
[-2.04764553, -1.53894872],
[-2.78919145, 1.04843516],
[-2.38144292, -1.03320686],
[-2.44730761, -0.98052759],
[-0.03204719, 2.73581902],
[ 0.25340191, 2.12459508],
[-0.08103629, 2.21583356],
[ 1.99593032, -1.15377649],
[ 2.03815488, -0.64748671],
[ 1.57778489, -2.43308223],
[-1.68532796, -2.20144187],
[ 0.60468723, -2.0147109 ],
[-2.28426076, 1.1401741 ],
[-1.2484405 , -2.35509733],
[-1.66227652, -1.37794664],
[ 2.71413302, 1.1687207 ]])
Then computing persistence is one line. Note that ripser defaults to only computing up to 1-dimensional persistence unless you tell it otherwise. For fun, I’m telling it to go up to 2-dimensional persistence.
output = ripser.ripser(P,maxdim=2)
dgms = output['dgms']
print(f"2-dimensional diagram points:")
print(dgms[2])
2-dimensional diagram points:
[[3.70902872 3.72065973]
[3.64157963 3.69366431]
[3.63115883 3.91910648]]
Then we can plot the resulting diagram.
plot_diagrams(dgms)
ax = plt.gca()
ax.grid(True, linestyle='--', linewidth=0.5, color='lightgray', alpha=0.9)
Sphere#
The following code will generate a point cloud drawn from a sphere. Note this is done with plotly, so will not render on the website. Download the jupyter notebook to see it.
P_sphere = gPC.Sphere(300)
# Create a 3d plot of the sphere using plotly
fig = px.scatter_3d(x=P_sphere[:,0], y=P_sphere[:,1], z=P_sphere[:,2])
fig.update_traces(marker=dict(size=3))
fig.show()
- What is the homology of the sphere?
- Compute the persistence diagram for the sphere point cloud. Does this fit with what you said above?
- Increase the number of points in the sphere to 400 and compute its persistence diagram. Do the same for 100 points and 500 points (although you might need to be patient, my laptop took about 30 seconds for that last one). What changes in the diagram??
# Your code here
Torus#
Here’s a point cloud sampled from a torus, with inner radius \(r=1\) and outer radius \(R=2\). Again note that the figure is done with plotly, so will not work on the website.
Image: Wikipedia
P_torus = gPC.Torus(200)
fig = px.scatter_3d(x=P_torus[:,0], y=P_torus[:,1], z=P_torus[:,2])
fig.update_traces(marker=dict(size=3))
fig.show()
- What is the homology of the torus?
- Compute the persistence diagram for the torus point cloud. Does this fit with what you said above? Hint: If it doesn't, try increasing the number of points in the point cloud.
- What is different between this diagram and the sphere diagram?
# Put your code here.
Properties of Rips Persistence#
Running times#
One important feature to notice is that you should only be computing persistence up to the dimension you care about since the running time gets way worse with more simplices.
P_torus = gPC.Torus(300)
# Note: If you want to see a major difference in running times, increase the number of points here.
%%time
# Here the code is computing up to 1 dimensional homology, and jupyter will print out the running time.
output = ripser.ripser(P_torus, maxdim=1)
CPU times: user 52.7 ms, sys: 6.46 ms, total: 59.1 ms
Wall time: 57.8 ms
%%time
# Here the code is computing up to 2 dimensional homology, and jupyter will print out the running time.
output = ripser.ripser(P_torus, maxdim=2)
CPU times: user 764 ms, sys: 29.8 ms, total: 793 ms
Wall time: 794 ms
%%time
# Here the code is computing up to 3 dimensional homology, and jupyter will print out the running time.
# Warning: this may take a while!
output = ripser.ripser(P_torus, maxdim=3)
CPU times: user 30.7 s, sys: 2.08 s, total: 32.7 s
Wall time: 33.2 s
Special case of 0-dimensional Rips persistence#
Note that in this setting, all vertices enter at \(r=0\), so all connected components are born at 0. This means the only thing interesting about 0-dimensional Rips persistence is the distribution of death times.
P_clusters = gPC.Clusters(N = 300, centers = np.array([[0,0], [4,4], [-1,6]]), sd = .3)
Draw.drawPtCloud(P_clusters)
ax = plt.gca()
ax.grid(True, linestyle='--', linewidth=0.5, color='lightgray', alpha=0.9)
Dgms_clusters = ripser.ripser(P_clusters, maxdim = 1)['dgms']
plot_diagrams(Dgms_clusters[0])
ax = plt.gca()
ax.grid(True, linestyle='--', linewidth=0.5, color='lightgray', alpha=0.9)
One option for visualization is to draw a histogram of the death times. In this case, I draw it horizontally since that matches with the death axis in the persistence diagrams.
Dgm0 = Dgms_clusters[0][:,1] # Just keep the death times
Dgm0 = Dgm0[~np.isinf(Dgm0)] # Drop the infinite bar because can't plot that
Dgm0
plt.hist(Dgm0,log=True,orientation = 'horizontal',bins = 20)
plt.xlabel('Number of points in the diagram')
plt.ylabel('Death time')
Text(0, 0.5, 'Death time')
Persistence for similar point clouds#
If we generate a bunch of random point clouds drawn from the same distribution, the persistence diagrams won’t be exactly the same, but they will have similar structure.
Here’s two point clouds, each pulled from an annulus with \(N = 200, r = 1, R = 2\) and draw their one-dimensional persistence diagrams on top of each other.
P1 = gPC.Annulus(200,1,2)
P2 = gPC.Annulus(200,1,2)
import matplotlib.pyplot as plt
plt.figure(figsize=(6,6))
plt.scatter(P1[:,0], P1[:,1], marker='*', color='purple', s=60, label='Annulus 1')
plt.scatter(P2[:,0], P2[:,1], marker='o', color='orange', s=40, label='Annulus 2')
plt.gca().set_aspect('equal', 'box')
plt.legend()
<matplotlib.legend.Legend at 0x16d606210>
We compute the persistence for the two point clouds first.
Dgms1 = ripser.ripser(P1)['dgms']
Dgms2 = ripser.ripser(P2)['dgms']
Then we can take a look at the two diagrams. Here’s the entirety of both diagrams in both dimension 0 and 1 side by side.
fig, ax = plt.subplots(1,2,figsize = (8,4))
plot_diagrams(Dgms1, ax = ax[0])
ax[0].set_title('Annulus 1')
plot_diagrams(Dgms2, ax = ax[1])
ax[1].set_title('Annulus 2')
Text(0.5, 1.0, 'Annulus 2')
In order to see the comparison better, let’s draw the two 1-dimensional diagrams on top of each other.
fig, ax = plt.subplots(1,2)
ax[0].scatter(P1[:,0], P1[:,1], marker = '*', color = 'purple', label = 'P1')
ax[0].scatter(P2[:,0], P2[:,1], color = 'orange', label = 'P2')
ax[0].legend()
ax[0].set_aspect('equal', 'box')
ax[0].set_title('Point Clouds')
ax[1].plot([0,2],[0,2], color='gray', linestyle='--')
ax[1].scatter(Dgms1[1][:,0], Dgms1[1][:,1], marker = '*', color = 'purple', label = 'P1')
ax[1].scatter(Dgms2[1][:,0], Dgms2[1][:,1], color = 'orange', label = 'P2')
ax[1].axis('square')
ax[1].set_title('1-D Persistence Diagrams')
ax[1].legend();
- What do you notice about the two persistence diagrams?
- Overlay the 1-dimensional torus diagram from earlier. What do you see? Which are more similar: the two annuli diagrams, or an annulus and the torus?
# Your code here
Changing the radius on the annulus#
Below is code to compute the persistence of the annulus with the default parameters.
P = gPC.Annulus(500,r = 1, R = 2)
Dgms = ripser.ripser(P)['dgms']
fig, axs = plt.subplots(1,2, figsize = (10,5))
plt.sca(axs[0])
Draw.drawPtCloud(P)
plt.sca(axs[1])
plot_diagrams(Dgms)
- What changes in the persistence diagram when you change the inside radius r?
- What changes in the persistence diagram when you change the outside radius R
The effect of outliers#
Consider the following point clouds \(P\) and \(Q\), where \(Q\) is just \(P\) with an extra point at the origin.
P = gPC.Annulus(300,r = 1, R = 2)
Q = np.concatenate([P,[[0,0]]])
DgmsP = ripser.ripser(P)['dgms']
DgmsQ = ripser.ripser(Q)['dgms']
fig, axs = plt.subplots(2,2, figsize = (8,8))
plt.sca(axs[0,0])
Draw.drawPtCloud(Q)
plt.axis('equal')
plt.sca(axs[0,1])
plot_diagrams(DgmsQ)
plt.sca(axs[1,0])
Draw.drawPtCloud(P)
plt.axis('equal')
plt.sca(axs[1,1])
plot_diagrams(DgmsP)
- What changes in the persistence diagram when you add the point?