<jupyter_start><jupyter_text>Rips Complex TutorialThis is a quick Rips Filtration tutorial used to illustrate options provided in BATS.py.<jupyter_code>import bats
import numpy as np
import matplotlib.pyplot as plt
import scipy.spatial.distance as distance
import bats
import time
# first, generate a circle
n = 100
X = np.random.normal(size=(n,2))
X = X / np.linalg.norm(X, axis=1).reshape(-1,1)
X = X + np.random.normal(size=(n,2), scale = 0.1 )
fig = plt.scatter(X[:,0], X[:,1])
fig.axes.set_aspect('equal')
plt.savefig('figures/RipsEx_data.png')<jupyter_output><empty_output><jupyter_text>Rips filtrations are commonly used in conjunction with persistent homology to create features for finite dimensional metric spaces (point clouds). Given a metric space $(X, d)$, a Rips complex consists of simplices with a maximum pairwise distance between vertices is less than some threshold $r$:$$X_r = \{(x_0,\dots,x_k) \mid x_i\in X, d(x_i,x_j) \le r\}.$$A Rips filtration is a filtration of Rips complexes $X_r \subseteq X_s$ if $r \le s$.<jupyter_code># compute pairwise distances
D = distance.squareform(distance.pdist(X))

# Rips complex for full metric space
# ie., threshold r = infinity now
F = bats.LightRipsFiltration(bats.Matrix(D), np.inf, 2)

# compute with F2 coefficents
t0 = time.monotonic()
R = bats.reduce(F, bats.F2())
t1 = time.monotonic()
print("time of compute persistent homology: {} sec.".format(t1-t0))<jupyter_output>time of compute persistent homology: 1.251744269 sec.<jupyter_text>Now you are able to see persistent diagrams.<jupyter_code># find persistence pairs at each dimension
ps = []
for d in range(R.maxdim()):
    ps.extend(R.persistence_pairs(d))

# Draw persistent diagram
# 'tmax' is the axis maixmum value
fig, ax = bats.persistence_diagram(ps, tmax = 2.0)<jupyter_output><empty_output><jupyter_text>Efficient Computation The number of simplices in Rips filtrations quickly grows with the size of the data set, and much effort has gone into developing efficient algorithms for computing persistent homology of Rips filtrations.   Construction optimizationThe first method that has been applied in several high-performance packages for Rips computations is to stop a filtration at the enclosing ***radius*** of the metric space, at which point the complex becomes contractible, which can reduce the total number of simplices in the filtration considerably without changing persistent homology.<jupyter_code># Two ways to find Enclosing Radius
r_enc = np.min(np.max(D, axis=0))
print("enclosing radius = {}".format(r_enc))

r_enc = bats.enclosing_radius(bats.Matrix(D))
print("enclosing radius = {}".format(r_enc))
# Rips complex up to enclosing radius
t0 = time.monotonic()
F_enc = bats.LightRipsFiltration(bats.Matrix(D), r_enc, 2)
t1 = time.monotonic()
print("construction time: {} sec.".format(t1-t0))

# compute with F2 coefficents
t0 = time.monotonic()
R_enc = bats.reduce(F_enc, bats.F2())
t1 = time.monotonic()
print("reduction time: {} sec.".format(t1-t0))<jupyter_output>construction time: 0.02580976899999987 sec.
reduction time: 0.45710929600000005 sec.<jupyter_text>You can see the obvious improvement with about 2x speedup.  Algorithm optimization There are also many efficent algorithms provided in BATS.py:- [**Reduction/Compression**](https://arxiv.org/abs/1303.0477): two options without basis returned.- [**Cohomology**](https://arxiv.org/abs/1107.5665): faster on some filtration.- [**Update Persistence**](https://arxiv.org/abs/2108.05022): suitable when there are several similar datasets needed to be computed PH (e.g., optimization on persistence penalty).  - **Extra Reduction**: peform extra reduction to eliminate nonzeros even after pivot has been found- **Combination of the above options**<jupyter_code>def time_BATS_flags(X, flags=(bats.standard_reduction_flag(), bats.compute_basis_flag())):
    
    t0 = time.monotonic()
    D = distance.squareform(distance.pdist(X))
    r_enc = bats.enclosing_radius(bats.Matrix(D))
    F_enc = bats.LightRipsFiltration(bats.Matrix(D), r_enc, 2)
    R = bats.reduce(F_enc, bats.F2(), *flags)
    t1 = time.monotonic()
    print("compute PH needs: {} sec.".format(t1 - t0))
flags = [
    (bats.standard_reduction_flag(), bats.compute_basis_flag()),
    (bats.standard_reduction_flag(),),
    (bats.standard_reduction_flag(), bats.clearing_flag()),
    (bats.standard_reduction_flag(), bats.compression_flag()),
    (bats.extra_reduction_flag(), bats.compute_basis_flag()),
    (bats.extra_reduction_flag(),),
    (bats.extra_reduction_flag(), bats.clearing_flag()),
    (bats.extra_reduction_flag(), bats.compression_flag()),
]
labels = [
    "standard w/ basis",
    "standard w/ no basis",
    "standard w/ clearing",
    "standard w/ compression",
    "extra w/ basis",
    "extra w/ no basis",
    "extra w/ clearing",
    "extra w/ compression"
]

for flag, label in zip(flags, labels):
    print("\n{}:".format(label))
    time_BATS_flags(X, flag)
# add some noise to the original datasets 
# to create a similar datasets to show the performance of update persistence
X2 = X + np.random.normal(size=(n,2), scale = 0.001)
# PH computation on X1
D = distance.squareform(distance.pdist(X))
r_enc = bats.enclosing_radius(bats.Matrix(D))
F_X = bats.LightRipsFiltration(bats.Matrix(D), r_enc, 2)
R = bats.reduce(F_X, bats.F2())

# PH computation on X2 by update persistence on X1
t0 = time.monotonic()
D2 = distance.squareform(distance.pdist(X2))
r_enc2 = bats.enclosing_radius(bats.Matrix(D2))
F_Y = bats.LightRipsFiltration(bats.Matrix(D2), r_enc2, 2) # generate a RipsFiltration
UI = bats.UpdateInfoLightFiltration(F_X, F_Y) # find updating information
R.update_filtration_general(UI)
t1 = time.monotonic()
print("compute PH by updating persistenc needs: {} sec.".format(t1 - t0))<jupyter_output>compute PH by updating persistenc needs: 0.1257570130000003 sec.
StarcoderdataJupyterStructuredCleanDedup