Example Usage

[2]:

import chalc.chromatic as chromatic
from chalc.sixpack import SubChromaticInclusion, SixPack
from chalc.plotting import plot_sixpack, plot_diagram, draw_filtration, animate_filtration
import numpy as np
import matplotlib.pyplot as plt

Computing chromatic filtrations

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

[4]:

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()
_images/example_5_0.svg

To quantify the spatial relationship between the orange and blue points, we need a filtration that contains blue and orange subfiltrations. The Čech and Vietoris–Rips filtrations satisfy this requirement, but they each have an enormous number of simplices (\(\displaystyle \binom{200}{2} = 19900\) edges and \(\displaystyle \binom{200}{3} = 1313400\) 2-simplices). We can compute the chromatic alpha/Delaunay–Čech/Delaunay–Rips complex filt of the point cloud using the functions chromatic.alpha, chromatic.delaunay_cech, and chromatic.delaunay_rips respectively. filt has far fewer simplices than the Čech and Vietoris–Rips filtrations, and contains subfiltrations spanned by the orange and blue points. This is the primary benefit of chromatic Delaunay filtrations.

[5]:

chromatic_alpha_filt  = chromatic.alpha(points, colours)
print(f'\'chromatic_alpha_filt\' has {len(chromatic_alpha_filt.simplices[1])} 1-simplices and {len(chromatic_alpha_filt.simplices[2])} 2-simplices.')

'chromatic_alpha_filt' has 954 1-simplices and 1312 2-simplices.

chromatic_alpha_filt is an object of type Filtration, and simplices in the filtration are objects of type Simplex. Both of these have convenience methods:

[6]:

# __len__()
print(f'\'chromatic_alpha_filt\' has {len(chromatic_alpha_filt)} simplices and {chromatic_alpha_filt.num_vertices} vertices.')
print(f"Dimension of the filtration: {chromatic_alpha_filt.dimension}.")

facets = []
# __iter__()
for simplex in chromatic_alpha_filt:
    # Do something with the simplex
    if simplex.dimension == 2:
        print(f"Found a simplex of dimension {simplex.dimension} with vertices {simplex.vertices} and colours {simplex.colours}")
        # Get a list of handles to the facets of the simplex
        facets = simplex.facets
        break

if facets:
    verts = facets[0].vertices
    # __contains__()
    if verts in chromatic_alpha_filt:
        print(f"{verts} is a simplex in the filtration.")
    else:
        print(f"{verts} is not a simplex in the filtration.")

'chromatic_alpha_filt' has 3023 simplices and 200 vertices.
Dimension of the filtration: 3.
Found a simplex of dimension 2 with vertices [30, 33, 87] and colours [0]
[33, 87] is a simplex in the filtration.

We can visualize the filtration at any given time using plotting.draw_filtration. For example, at time 0.3 we have:

[7]:

fig, ax = plt.subplots(1, 3, figsize=(12, 4))
draw_filtration(chromatic_alpha_filt, points, time=0.3, include_colours=[0], ax=ax[0])
ax[0].set_title('Colour 0')
draw_filtration(chromatic_alpha_filt, points, time=0.3, include_colours=[1], ax=ax[1])
ax[1].set_title('Colour 1')
draw_filtration(chromatic_alpha_filt, 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()
_images/example_11_0.svg

Computing 6-packs of persistence diagrams

Next we quantify how the orange and blue points are related to each other spatially using persistent homology. To do this we need more than just the persistent homology of the filtration as a whole. We should also consider the persistent homology of the individual subfiltrations, as well the kernel, cokernel, and image of the map of persistence modules induced by the inclusion of the blue subfiltration into the orange subfiltration (or the other way around). Since our persistent homology modules arise from actual simplicial filtrations, we can also compute the relative homology persistence modules induced by these inclusions. The kernel, cokernel, image, domain, codomain, and relative diagram are bundled together into a “six-pack” of persistence diagrams. For this example, we will consider the inclusion of the blue points (colour 0) into all points by using the SubChromaticInclusion class.

[8]:

# Using the previously computed filtration:
dgms_alpha = SubChromaticInclusion(chromatic_alpha_filt, {0}).sixpack()
# Using the Delaunay--Čech and Delaunay--Rips filtrations
chromatic_delcech_filt = chromatic.delaunay_cech(points, colours)
chromatic_delrips_filt = chromatic.delaunay_rips(points, colours)
dgms_delcech = SubChromaticInclusion(chromatic_delcech_filt, {0}).sixpack()
dgms_delrips = SubChromaticInclusion(chromatic_delrips_filt, {0}).sixpack()

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.

[9]:

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.

[10]:

print(dgms_alpha["ker"])

Paired: frozenset({(466, 470), (682, 692), (707, 709), (439, 663), (737, 750), (283, 342), (669, 673), (582, 660), (494, 495), (228, 235), (456, 480), (396, 594), (651, 686), (804, 806), (281, 320), (560, 644), (262, 380), (339, 587), (848, 849), (349, 384), (525, 604), (771, 772), (616, 621), (373, 393), (592, 828), (799, 800), (657, 658), (854, 855), (364, 457), (667, 684), (500, 543), (252, 441), (362, 417), (561, 562), (357, 555), (698, 711), (419, 432), (293, 390), (590, 625), (749, 786), (426, 505), (352, 540), (433, 516), (361, 414), (337, 420), (648, 654), (614, 634)})
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.

[11]:

# __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:

[12]:

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.03060548]
 [0.         0.01173629]
 ...
 [0.         0.0164209 ]
 [0.         0.04228335]
 [0.                inf]]
Paired: frozenset({(466, 470), (682, 692), (707, 709), (439, 663), (737, 750), (283, 342), (669, 673), (582, 660), (494, 495), (228, 235), (456, 480), (396, 594), (651, 686), (804, 806), (281, 320), (560, 644), (262, 380), (339, 587), (848, 849), (349, 384), (525, 604), (771, 772), (616, 621), (373, 393), (592, 828), (799, 800), (657, 658), (854, 855), (364, 457), (667, 684), (500, 543), (252, 441), (362, 417), (561, 562), (357, 555), (698, 711), (419, 432), (293, 390), (590, 625), (749, 786), (426, 505), (352, 540), (433, 516), (361, 414), (337, 420), (648, 654), (614, 634)})
Unpaired: frozenset()
Paired: frozenset({(466, 470), (682, 692), (707, 709), (439, 663), (737, 750), (283, 342), (669, 673), (582, 660), (494, 495), (228, 235), (456, 480), (396, 594), (651, 686), (804, 806), (281, 320), (560, 644), (262, 380), (339, 587), (848, 849), (349, 384), (525, 604), (771, 772), (616, 621), (373, 393), (592, 828), (799, 800), (657, 658), (854, 855), (364, 457), (667, 684), (500, 543), (252, 441), (362, 417), (561, 562), (357, 555), (698, 711), (419, 432), (293, 390), (590, 625), (749, 786), (426, 505), (352, 540), (433, 516), (361, 414), (337, 420), (648, 654), (614, 634)})
Unpaired: frozenset()

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

[13]:

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

[0.         0.         0.         ... 0.23772323 0.78357409 0.78357409]
[0 0 0 ... 2 2 3]

Plotting 6-packs of persistence diagrams

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.

[14]:

fig1, ax1 = plot_sixpack(dgms_alpha, threshold = 1e-3, truncations=0.5)
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.
plot_diagram(dgms_alpha, 'rel', dimensions = {0, 2}, truncation = 0.5, ax = ax2[0], threshold = 1e-3)
ax2[0].set_title('Chromatic Alpha')
# You can also specify a single dimension as an integer.
plot_diagram(dgms_delcech, 'rel', dimensions = 1, truncation = 0.5, ax = ax2[1], threshold = 1e-3)
ax2[1].set_title('Chromatic Delaunay-Čech')
# If no dimensions are specified, all features will be included.
plot_diagram(dgms_delrips, 'rel', truncation = 0.5, ax = ax2[2], threshold = 1e-3)
ax2[2].set_title('Chromatic Delaunay-Rips')
fig2.suptitle('Individual relative diagrams in specific dimensions')
fig2.set_figwidth(15)
plt.show()
_images/example_27_0.svg
_images/example_27_1.svg

Saving diagrams to file

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

[15]:

import h5py

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

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

print(dgms_alpha_from_file == dgms_alpha)

True

Other mapping methods

In the example given previous, SubChromaticInclusion was used to specify the inclusion of some subset of colours into the whole filtration. We can even specify multiple subsets of colours to include. For example, 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 SubChromaticInclusion([[0, 1], [1, 3], [2]]) above. We could also consider other mapping methods:

  1. The \(k\)-chromatic subfiltration of filt is the subfiltration spanned by all simplices having at most \(k\) colours. We can specify the inclusion of the \(k\)-chromatic subfiltration into filt by using the class KChromaticInclusion.

  2. We could choose a family of subchromatic filtrations (each of which given as in the example above), and consider the union of their inclusions into filt (the sub-chromatic quotient map). To do this we can use SubChromaticQuotient.

  3. The subchromatic quotient map over the family of subfiltrations spanned by the various \(k\)-subsets of colours is specified with KChromaticQuotient.

Animating the filtration

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

[16]:

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

[16]:
_images/example_33_1.svg