Example Usage

import chalc as ch
import numpy as np
import matplotlib.pyplot as plt

For our data we sample 100 points on a circle with some noise and 100 points from inside the unit disk.

rng = np.random.default_rng(40)
num_points_circle = 100
num_points_disk = 100
mean = [0, 0]
cov = np.eye(2)*0.01
circle = np.array([[np.sin(2*np.pi*t), np.cos(2*np.pi*t)] for t in rng.random(num_points_circle)]).T +\
    rng.multivariate_normal(mean, cov, num_points_circle).T # points as columns
disk = np.sqrt(rng.random(num_points_disk)) * np.array([[np.sin(2*np.pi*t), np.cos(2*np.pi*t)] for t in rng.random(num_points_disk)]).T
points = np.concatenate((circle, disk), axis=1)
plt.scatter(circle[0, :], circle[1, :], s=10)
plt.scatter(disk[0, :], disk[1, :], s=10)
plt.gca().set_aspect('equal')
colours = [0]*num_points_circle + [1]*num_points_disk
plt.title('Point cloud')
plt.show()

We can compute the chromatic alpha/Delaunay–Cech/Delaunay–Rips complex $K$ of the point cloud using the functions chromatic.alpha, chromatic.delcech, and chromatic.delrips respectively. In each of these cases, \(K\) has far fewer simplices than either the Čech or Vietoris–Rips complex, which each have \(\displaystyle \binom{200}{2} = 19900\) edges and \(\displaystyle \binom{200}{3} = 1313400\) 2-simplices.

K, _ = ch.chromatic.alpha(points, colours)
printmd(f'$K$ has {len(K.simplices[1])} 1-simplices and {len(K.simplices[2])} 2-simplices')

\(K\) has 954 1-simplices and 1312 2-simplices

\(K\) is a FilteredComplex object, and simplices in the filtration are objects of type Simplex (although you should not generally need to interact with the API of these objects). We can visualise the filtration at any given time using plotting.draw_filtration. For example, at time 0.3 we have:

fig, ax = plt.subplots(1, 3, figsize=(12, 4))
ch.plotting.draw_filtration(K, points, time=0.3, include_colours=[0], ax=ax[0])
ax[0].set_title('Colour 0')
ch.plotting.draw_filtration(K, points, time=0.3, include_colours=[1], ax=ax[1])
ax[1].set_title('Colour 1')
ch.plotting.draw_filtration(K, points, time=0.3, ax=ax[2])
ax[2].set_title('Both colours')
for i in range(3):
    ax[i].set_xlim(-1.2, 1.2)
    ax[i].set_ylim(-1.2, 1.2)
fig.suptitle('Complexes at $t=0.3$')
plt.tight_layout()
plt.show()

The functions sixpack.from_filtration and sixpack.compute compute the 6-pack of persistent homology diagrams of the chromatic alpha/Delaunay–Čech complex/Delaunay–Rips complexes.

# using the previously computed filtration
dgms_alpha = ch.sixpack.from_filtration(K, dom=[0])
# or directly from the point cloud
dgms_delcech = ch.sixpack.compute(points, colours, dom=[0], method="delcech")
dgms_delrips = ch.sixpack.compute(points, colours, dom=[0], method="delrips")

Each 6-pack is an object of type DiagramEnsemble that has the kernel, cokernel, domain, codomain, and relative persistent homology diagrams. The class DiagramEnsemble implements the abstract base class Mapping from the Python standard library and supports dictionary-like syntax. In particular you can access the diagrams by indexing into the 6-pack.

print(f"The diagram names are: {list(dgms_alpha.keys())}")

# iterate over the diagrams like in a dictionary
for diagram_name, diagram in dgms_alpha:
    print(diagram_name)
    # do something with the diagram

# access individual diagrams by indexing
ker = dgms_alpha['ker']

The diagram names are: ['ker', 'cok', 'dom', 'cod', 'im', 'rel']
ker
cok
dom
cod
im
rel

Every feature in the diagram is represented by either a pair of simplices \((\sigma, \tau)\) if the feature has a finite birth and death time, or an unpaired simplex \(\sigma\) if the feature does not die. For a \(p\)-dimensional feature, \(\dim(\sigma) = p+1\) if the diagram is the kernel, and \(\dim(\sigma) = p\) otherwise.

print(dgms_alpha['ker'])

Paired: {(1187, 1200), (358, 409), (1445, 1453), (584, 883), (1697, 1701), (1995, 1999), (228, 235), (333, 581), (361, 456), (454, 471), (1648, 1650), (1486, 1499), (1447, 1495), (1461, 1494), (1565, 1600), (429, 502), (1510, 1514), (606, 626), (1360, 1458), (1248, 1252), (1393, 1403), (1267, 1268), (387, 587), (487, 531), (506, 596), (1454, 1457), (292, 383), (348, 528), (1329, 1331), (421, 491), (252, 438), (1592, 1671), (280, 317), (1366, 1408), (359, 412), (1464, 1470), (2013, 2433), (345, 376), (282, 338), (413, 428), (544, 646), (352, 540), (262, 372), (1692, 1693), (436, 658), (1505, 1516), (331, 416)}
Unpaired: set()

Each individual diagram is an object of type SimplexPairings, which supports set-like syntax and implements the abstract base class Collection from the Python standard library.

# check if the diagram is empty
print("The kernel is non-empty." if dgms_alpha['ker'] else "The kernel is empty.")
# or you can use bool(dgms_alpha.ker)

# __len__()
# get the number of features in the diagram
print(f"The cokernel has {len(dgms_alpha['cok'])} features.")

# __contains__()
# check if a feature exists in a diagram
print("The domain contains a feature represented by the simplex pair (20, 40): "
      f"{(20, 40) in dgms_alpha['dom']}")

# __iter__()
# iterate over the paired and unpaired simplices together
for feature in dgms_alpha['ker']:
    if isinstance(feature, tuple):
        sigma, tau = feature
        # do something with a pair of simplices
        ...
    else:
        sigma = feature
        # do something with an unpaired simplex
        ...

# the paired and unpaired simplices can be considered separately
print("The \"paired\" property of the kernel has type: "
      f"{type(dgms_alpha['ker'].paired)}")

print("The \"unpaired\" property of the kernel has type: "
      f"{type(dgms_alpha['ker'].unpaired)}")

for sigma, tau in dgms_alpha['ker'].paired:
    # do something with paired simplices
    ...

for sigma in dgms_alpha['dom'].unpaired:
    # do something with unpaired simplices
    ...

The kernel is non-empty.
The cokernel has 234 features.
The domain contains a feature represented by the simplex pair (20, 40): False
The "paired" property of the kernel has type: <class 'set'>
The "unpaired" property of the kernel has type: <class 'set'>

A convenience function get_matrix is provided to extract the diagrams as a matrix of birth/death times:

np.set_printoptions(threshold=10)
print(dgms_alpha.get_matrix('dom', 0)) # get the domain in dimension zero
# dgms_alpha.get('ker', [0, 1]) to get the kernel in dimensions zero and one
# dgms_alpha.get('ker') to get the kernel in all dimensions from zero to max(dgms_alpha.dimensions)

[[0.         0.03060548]
 [0.         0.01367401]
 [0.         0.09263013]
 ...
 [0.         0.0164209 ]
 [0.         0.03560387]
 [0.                inf]]

You can also access the entrance_times and dimensions of simplices in the 6-pack diagrams.

print(dgms_alpha.entrance_times)
print(dgms_alpha.dimensions)

[ 0.          0.          0.         ...  3.94302018  5.36996584
 16.97606576]
[0 0 0 ... 3 3 3]

The convenience functions plotting.plot_sixpack and plotting.plot_diagram are provided for plotting either a 6-pack or a specific diagram from a 6-pack.

fig1, ax1 = ch.plotting.plot_sixpack(dgms_alpha, tolerance = 1e-3)
fig1.suptitle('Chromatic Alpha')
fig1.set_figwidth(15)
fig1.set_figheight(9)
fig1.subplots_adjust(top=0.92)
plt.show()
fig2, ax2 = plt.subplots(1, 3)
# You can specify the dimensions of features to include in the plot.
ch.plotting.plot_diagram(dgms_alpha, 'rel', dimensions = {0, 2}, truncation = 0.3, ax = ax2[0], tolerance = 1e-3)
ax2[0].set_title('Chromatic Alpha')
# You can also specify a single dimension as an integer.
ch.plotting.plot_diagram(dgms_delcech, 'rel', dimensions = 1, truncation = 0.3, ax = ax2[1], tolerance = 1e-3)
ax2[1].set_title('Chromatic Delaunay-Cech')
# If no dimensions are specified, all features will be included.
ch.plotting.plot_diagram(dgms_delrips, 'rel', truncation = 0.3, ax = ax2[2], tolerance = 1e-3)
ax2[2].set_title('Chromatic Delaunay-Rips')
fig2.suptitle('Individual relative diagrams')
fig2.set_figwidth(15)
plt.show()

We can save the diagrams to a HDF5 file or load a diagram from a HDF5 file. The HDF5 format is

with h5py.File('test.h5', 'w') as f:
    dgms_alpha.save(f)

with h5py.File('test.h5', 'r') as f:
    dgms_alpha_from_file = ch.sixpack.DiagramEnsemble.from_file(f)

print(dgms_alpha_from_file == dgms_alpha)

True

We can visualise the 2-skeleton of the filtration for points in 2D:

animation = ch.plotting.animate_filtration(
    K, points, filtration_times=list(np.linspace(0, 1.0, 45)), animation_length=5)
animation

Once Loop Reflect