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–Čech/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.

# max_num_threads=0 gives automatic parallelism
K, numerical_instability_flag = ch.chromatic.alpha(points, colours, max_num_threads=0)
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 visualize 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(filtration=K, mapping_method=ch.sixpack.SubChromaticInclusion({0}))
# Or directly from the point cloud:
dgms_delcech = ch.sixpack.compute(points, colours, mapping_method=ch.sixpack.SubChromaticInclusion({0}), filtration_algorithm="delcech")
dgms_delrips = ch.sixpack.compute(points, colours, mapping_method=ch.sixpack.SubChromaticInclusion({0}), filtration_algorithm="delrips")

In code above, SubChromaticInclusion corresponds to the inclusion of some subset (or subsets) of colours into the whole filtration. In this example, we are considering the map of filtrations induced by the inclusion of colour 0 into the union of colours {0, 1}. If we had more colours, say {0, 1, 2, 3}, we might also consider the inclusion of the sub-filtration spanned by simplices whose colours lie in any subset among (say) {{0, 1}, {1, 3}, {2}}. This would be specified by setting mapping_method=SubChromaticInclusion([[0, 1], [1, 3], [2]]) above. You can also consider other mapping methods; see KChromaticInclusion and KChromaticQuotient.

Each 6-pack is an object of type SixPack that has the kernel, cokernel, domain, codomain, and relative persistent homology diagrams. The class SixPack 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 in dgms_alpha:
    print(diagram_name)
    # do something with dgms_alpha[diagram_name]

# 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: frozenset({(379, 737), (367, 418), (289, 358), (996, 1010), (494, 623), (228, 235), (2118, 2957), (1498, 1504), (562, 568), (390, 477), (1080, 1110), (507, 644), (479, 506), (618, 707), (1035, 1047), (828, 919), (392, 553), (1122, 1143), (742, 986), (550, 584), (349, 484), (1403, 1406), (1295, 1446), (653, 862), (1133, 1139), (519, 1023), (254, 521), (299, 428), (2084, 2091), (264, 412), (884, 938), (401, 433), (814, 1020), (832, 1754), (372, 702), (743, 746), (438, 842), (886, 907), (1245, 1310), (287, 328), (1013, 1018), (1489, 1490), (1029, 1096), (351, 825), (999, 1098), (608, 609), (389, 470)})
Unpaired: frozenset()

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.

# __bool__()
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.
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 'frozenset'>
The "unpaired" property of the kernel has type: <class 'frozenset'>

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

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

[[0.         0.01367401]
 [0.         0.01173629]
 [0.         0.02611904]
 ...
 [0.         0.04149252]
 [0.         0.06114399]
 [0.                inf]]
Paired: frozenset({(379, 737), (367, 418), (289, 358), (996, 1010), (494, 623), (228, 235), (2118, 2957), (1498, 1504), (562, 568), (390, 477), (1080, 1110), (507, 644), (479, 506), (618, 707), (1035, 1047), (828, 919), (392, 553), (1122, 1143), (742, 986), (550, 584), (349, 484), (1403, 1406), (1295, 1446), (653, 862), (1133, 1139), (519, 1023), (254, 521), (299, 428), (2084, 2091), (264, 412), (884, 938), (401, 433), (814, 1020), (832, 1754), (372, 702), (743, 746), (438, 842), (886, 907), (1245, 1310), (287, 328), (1013, 1018), (1489, 1490), (1029, 1096), (351, 825), (999, 1098), (608, 609), (389, 470)})
Unpaired: frozenset()
Paired: frozenset({(379, 737), (367, 418), (289, 358), (996, 1010), (494, 623), (228, 235), (2118, 2957), (1498, 1504), (562, 568), (390, 477), (1080, 1110), (507, 644), (479, 506), (618, 707), (1035, 1047), (828, 919), (392, 553), (1122, 1143), (742, 986), (550, 584), (349, 484), (1403, 1406), (1295, 1446), (653, 862), (1133, 1139), (519, 1023), (254, 521), (299, 428), (2084, 2091), (264, 412), (884, 938), (401, 433), (814, 1020), (832, 1754), (372, 702), (743, 746), (438, 842), (886, 907), (1245, 1310), (287, 328), (1013, 1018), (1489, 1490), (1029, 1096), (351, 825), (999, 1098), (608, 609), (389, 470)})
Unpaired: frozenset()

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.         ... 16.97606576 16.97606576
 16.97606576]
[0 0 0 ... 2 2 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, threshold = 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], threshold = 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], threshold = 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], threshold = 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.SixPack.from_file(f)

print(dgms_alpha_from_file == dgms_alpha)

True

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

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

Once Loop Reflect