Example Usage¶
[2]:
from chalc.morphisms import SixPackMorphisms
import matplotlib.pyplot as plt
import numpy as np
np.set_printoptions(threshold=10)
from itertools import chain, combinations
from dataclasses import fields
from chalc import chromatic
from chalc.plotting import (
animate_filtration,
draw_filtration,
plot_persistence_diagram,
plot_sixpack,
)
from chalc.morphisms import (
KChromaticGluing,
SixPack,
SubChromaticInclusion,
SixPackRepresentatives,
)
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()
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 using plotting.draw_filtration.
A selector decides which simplices to draw. For example, at time 0.3 we can include one colour or both:
[7]:
fig, ax = plt.subplots(1, 3, figsize=(12, 4), layout="compressed")
draw_filtration(
chromatic_alpha_filt,
points,
lambda simplex, _: simplex.filtration_value <= 0.3 and set(simplex.colours) <= {0},
ax=ax[0],
)
ax[0].set_title("Colour 0")
draw_filtration(
chromatic_alpha_filt,
points,
lambda simplex, _: simplex.filtration_value <= 0.3 and set(simplex.colours) <= {1},
ax=ax[1],
)
ax[1].set_title("Colour 1")
draw_filtration(
chromatic_alpha_filt,
points,
lambda simplex, _: simplex.filtration_value <= 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.show()
Computing 6-packs of persistence modules¶
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 a single-colour subfiltration into the whole filtration.
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 persistent homology modules are bundled together into a “6-pack” of persistence modules.
For this example, we will consider the inclusion of the blue points (colour 0) into all points by using the SubChromaticInclusion class.
A 6-pack can be extracted from any instance of this class by calling its method sixpack.
Pass the argument max_homological_degree to limit homological degrees for which persistence modules are computed.
[8]:
# Using the previously computed filtration:
sixpack_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)
sixpack_delcech = SubChromaticInclusion(chromatic_delcech_filt, {0}).sixpack()
sixpack_delrips = SubChromaticInclusion(chromatic_delrips_filt, {0}).sixpack()
sixpack_delrips_max_deg_1 = SubChromaticInclusion(chromatic_delrips_filt, {0}).sixpack(
max_homological_degree=1
)
Each 6-pack is an object of type SixPack containing the kernel, cokernel, domain, codomain, image, and relative persistent homology modules.
These modules are instance variables of type AugmentedPersistenceModule.
A convenience method max_non_empty_degree gives the maximum degree of any interval in any persistence module in the 6-pack.
[9]:
print(f"The persistence module names are: {SixPack.module_names()}")
for module_name, module in sixpack_alpha:
print(f"{module_name}: {len(module)} features")
print(f"The cokernel has {len(sixpack_alpha.cokernel)} features.")
assert sixpack_delrips_max_deg_1.max_nonempty_degree() <= 1
("The persistence module names are: ('kernel', 'domain', 'image', 'codomain', "
"'cokernel', 'relative')")
'kernel: 47 features'
'domain: 155 features'
'image: 141 features'
'codomain: 372 features'
'cokernel: 234 features'
'relative: 274 features'
'The cokernel has 234 features.'
Each interval in each of the modules is represented as a pairing of a birth index to a death index, or to None if the feature persists indefinitely.
These indices refer to column indices in the boundary matrix of the underlying codomain filtration (or the mapping cylinder for gluing morphisms).
A kernel feature has homological degree one less than its birth-column degree, while other features have the degree of their birth column.
For more details on the boundary matrix ordering and format, see boundary_matrix.
[10]:
print(f"The domain has a feature born at index 20: {20 in sixpack_alpha.pairings.domain}")
kernel_pairings = sixpack_alpha.pairings.kernel
for birth_index, death_index in kernel_pairings.items():
if death_index is None:
# This feature has infinite death time.
...
else:
# This feature dies at death_index.
...
print(repr(kernel_pairings))
'The domain has a feature born at index 20: True'
'{228: 235, 254: 521, 264: 412, 287: 328, ...}'
Extracting persistence modules from the 6-pack¶
An AugmentedPersistenceModule is a persistence module of graded \(\mathbb{F}_2\) vector spaces, annotated with its interval decomposition, and a choice of generator-relation pairings that yield this interval decomposition.
The generators and relations are strictly ordered by entrance time (ties resolved arbitrarily), so that the pairings are stored as pairings of integer indices and can be accessed through pairings.
These pairing indices are NOT equal to the pairing indices from the parent 6-pack, since they are re-indexed to remove gaps in the ordering.
They are, however, compatible with the keys of the representatives and morphisms in the parent 6-pack (see the next section).
Use the pairing indices from the parent 6-pack only if you need to index into the boundary matrix of the filtration that generated the 6-pack;
otherwise, use the module pairing indices since they are independent of the ambient boundary matrix coordinates.
The attribute degrees maps each interval (identified by its birth index) to its degrees.
Use intervals to extract the interval decomposition as a matrix of birth/death times, optionally specifying degrees in which to extract intervals.
A single integer degree returns an \(m \times 2\) matrix, while a collection of degrees returns a dictionary of matrices keyed by the requested degrees.
If there are no intervals in the specified degree, an empty matrix with shape (0, 2) is returned.
Omitting degrees returns a dictionary of all degrees containing any intervals.
[11]:
kernel_module = sixpack_alpha.kernel
print("The kernel is non-empty." if kernel_module else "The kernel is empty.")
print(repr(kernel_module.pairings))
print(repr(kernel_module.degrees))
assert kernel_module.degrees.keys() == kernel_module.pairings.keys()
np.set_printoptions(threshold=10)
domain_module = sixpack_alpha.domain
# The domain diagram in degree zero.
print(domain_module.intervals(0))
# The kernel diagram in degrees zero and one.
print(kernel_module.intervals([0, 1]))
# The kernel diagram in all inhabited degrees.
print(kernel_module.intervals())
'The kernel is non-empty.'
'mappingproxy(...: 20, 15: 24})'
'mappingproxy(... 6: 0, 15: 0})'
array([[0. , inf],
[0. , 0.0645652 ],
[0. , 0.17278701],
...,
[0. , 0.02176687],
[0. , 0.04707739],
[0. , 0.03121891]], shape=(100, 2))
{0: array([[0.022234 , 0.02443656],
[0.03181886, 0.06730145],
[0.03412368, 0.05823877],
...,
[0.08579831, 0.10430597],
[0.09189952, 0.17278701],
[0.09594719, 0.10077087]], shape=(25, 2)),
1: array([[0.0566658 , 0.06077279],
[0.07147837, 0.07166417],
[0.07494828, 0.07500632],
...,
[0.14259426, 0.14317964],
[0.22895679, 0.23057114],
[0.23772323, 0.78357409]], shape=(22, 2))}
{0: array([[0.022234 , 0.02443656],
[0.03181886, 0.06730145],
[0.03412368, 0.05823877],
...,
[0.08579831, 0.10430597],
[0.09189952, 0.17278701],
[0.09594719, 0.10077087]], shape=(25, 2)),
1: array([[0.0566658 , 0.06077279],
[0.07147837, 0.07166417],
[0.07494828, 0.07500632],
...,
[0.14259426, 0.14317964],
[0.22895679, 0.23057114],
[0.23772323, 0.78357409]], shape=(22, 2))}
A module also exposes the sorted entrance_times and the degrees of its generators and relations.
[12]:
print(kernel_module.entrance_times)
print(kernel_module.nonempty_degrees())
print(f"Largest inhabited degree in the 6-pack: {sixpack_alpha.max_nonempty_degree()}")
array([0.022234 , 0.02443656, 0.03181886, ..., 0.23057114, 0.23772323,
0.78357409], shape=(94,))
array([0, 1])
'Largest inhabited degree in the 6-pack: 2'
Plotting 6-packs of persistence diagrams¶
Use plot_sixpack to plot all six diagrams and plot_persistence_diagram to plot an augmented persistence module.
A persistence threshold is applied to the original intervals to discard low-persistence features, and the displayed coordinates are truncated to fit in the plot.
[13]:
fig1, ax1 = plot_sixpack(sixpack_alpha, threshold=1e-3, truncations=0.5)
fig1.suptitle("Chromatic Alpha")
fig1.set_figwidth(15)
fig1.set_figheight(9)
plt.show()
Specify homological degrees with a single integer or a collection of integers.
If degrees is omitted, all inhabited degrees are plotted.
[14]:
fig, ax = plt.subplots(1, 3)
plot_persistence_diagram(
sixpack_alpha.relative, degrees={0, 2}, truncation=0.5, ax=ax[0], threshold=1e-3
)
ax[0].set_title("Chromatic Alpha")
plot_persistence_diagram(
sixpack_delcech.relative, degrees=1, truncation=0.5, ax=ax[1], threshold=1e-3
)
ax[1].set_title("Chromatic Delaunay-Čech")
plot_persistence_diagram(sixpack_delrips.relative, truncation=0.5, ax=ax[2], threshold=1e-3)
ax[2].set_title("Chromatic Delaunay-Rips")
fig.suptitle("Individual relative diagrams in specific degrees")
fig.set_figwidth(15)
plt.show()
Customize the plotting across all degrees by passing keyword arguments to the scatter plot.
Use kwargs_per_dim for a dictionary of scatter options keyed by homological degree.
Options for a degree take precedence over global options.
[15]:
fig, ax = plt.subplots(1, 3)
plot_persistence_diagram(
sixpack_alpha.relative,
degrees={0, 2},
truncation=0.5,
ax=ax[0],
threshold=1e-3,
kwargs_per_dim={
0: {"marker": "o", "color": "black"},
2: {"marker": "^", "color": "orange", "s": 100},
},
)
ax[0].set_title("Chromatic Alpha")
plot_persistence_diagram(
sixpack_delcech.relative, degrees=1, truncation=0.5, ax=ax[1], threshold=1e-3, color="k"
)
ax[1].set_title("Chromatic Delaunay-Čech")
plot_persistence_diagram(
sixpack_delrips.relative,
truncation=0.5,
ax=ax[2],
threshold=1e-3,
kwargs_per_dim={1: {"marker": "o", "s": 40}},
marker="X",
s=100,
)
ax[2].set_title("Chromatic Delaunay-Rips")
fig.suptitle("Individual relative diagrams in specific degrees")
fig.set_figwidth(15)
plt.show()
Saving diagrams to file¶
We can save the diagrams to a HDF5 file or load a diagram from a HDF5 file.
[16]:
import h5py
from pathlib import Path
from tempfile import TemporaryDirectory
with TemporaryDirectory() as directory:
path = Path(directory) / "diagrams.h5"
with h5py.File(path, "w") as file:
sixpack_alpha.to_hdf5(file)
with h5py.File(path, "r") as file:
sixpack_alpha_from_file = SixPack.from_hdf5(file)
assert sixpack_alpha_from_file == sixpack_alpha
Representative cycles and persistent morphisms¶
We can extract representative cycles for each interval in each persistence module in the 6-pack, and the natural maps between the persistence modules in the 6-pack.
To do this, pass the argument get_representatives=True when calling sixpack.
[17]:
sixpack_alpha = SubChromaticInclusion(chromatic_alpha_filt, {0}).sixpack(get_representatives=True)
representatives = sixpack_alpha.representatives
assert isinstance(representatives, SixPackRepresentatives)
The representatives are a dataclass mapping names of persistence modules in the 6-pack to a collection of representative cycles. Each collection is a dictionary mapping intervals in the module (identified by their birth index in the module) to a one-dimensional NumPy array of integer indices of columns in the boundary matrix of the filtration.
[18]:
print(list(field.name for field in fields(representatives)))
assert representatives.domain.keys() == sixpack_alpha.domain.pairings.keys()
print(repr(representatives.domain))
['kernel', 'domain', 'image', 'codomain', 'cokernel', 'relative']
'{0: array([0]), 1: array([0, 1]), 2: array([0, 2]), 3: array([1, 3]), ...}'
The morphisms are a dataclass mapping the name of the morphism to a collection of interval mappings. For example, to get the map from the domain to codomain:
[19]:
morphisms = sixpack_alpha.morphisms
assert isinstance(morphisms, SixPackMorphisms)
print(list(field.name for field in fields(morphisms)))
assert morphisms.domain_to_codomain.keys() == sixpack_alpha.domain.pairings.keys()
codomain_intervals_hit = {
codomain_interval
for image_of_domain_interval in morphisms.domain_to_codomain.values()
for codomain_interval in image_of_domain_interval
}
assert codomain_intervals_hit.issubset(sixpack_alpha.codomain.pairings.keys())
print(repr(morphisms.domain_to_codomain))
['kernel_to_domain',
'domain_to_image',
'image_to_codomain',
'domain_to_codomain',
'codomain_to_cokernel']
'{0: array([0]), 1: array([1]), 2: array([1, 2]), 3: array([2, 3]), ...}'
Plotting representative cycles¶
The selector passed to draw_filtration receives a simplex handle and its index in the filtration’s boundary matrix. It returns True for each simplex to draw. We can use this to plot representative cycles for each diagram in the 6-pack.
For example, when we include the blue points into the union of blue and orange points, there is a persistent cycle in the kernel that corresponds to the cycle formed by the blue points. The code below uses the representative blue cycle returned by :meth:sixpack(get_representatives=True) <chalc.morphisms.FiltrationMorphism.sixpack> to visualize this blue cycle.
The representative dictionary uses generator indices as keys. Its values are arrays of boundary-matrix indices. We choose the degree-one representative with the longest kernel interval.
[20]:
sixpack_alpha = SubChromaticInclusion(chromatic_alpha_filt, {0}).sixpack(get_representatives=True)
assert sixpack_alpha.representatives is not None
boundary_matrix = chromatic_alpha_filt.boundary_matrix()
kernel = sixpack_alpha.kernel
longest_interval_len = 0.0
longest_interval_birth_index = -1
for i in kernel.pairings:
if kernel.degrees[i] == 1:
birth, death = kernel.interval(i)
if death is None:
longest_interval_birth_index = i
break
if death - birth > longest_interval_len:
longest_interval_len = death - birth
longest_interval_birth_index = i
interval_birth, interval_death = kernel.interval(longest_interval_birth_index)
representative_indices = sixpack_alpha.representatives.kernel[longest_interval_birth_index]
fig, ax = plt.subplots(figsize=(6, 6))
# Draw the filtration in the background
draw_filtration(
chromatic_alpha_filt,
points,
select=lambda simplex, _: simplex.filtration_value <= interval_birth,
ax=ax,
edges_kwargs={"alpha": 0.15},
faces_kwargs={"alpha": 0.10},
)
# Overlay the representatives
draw_filtration(
chromatic_alpha_filt,
points,
select=lambda _, index: index in representative_indices,
ax=ax,
edges_kwargs={"alpha": 1.0, "linewidth": 1.8, "zorder": 4},
)
ax.set_title(
f"SubChromaticInclusion: kernel cycle\n"
f"Kernel interval [{interval_birth:.3f}, {interval_death:.3f})"
)
plt.show()
Animating the filtration¶
We can visualize the 2-skeleton of a chromatic filtration for points in 2D:
[21]:
animation = animate_filtration(
chromatic_alpha_filt,
points,
filtration_times=np.linspace(0, 1.0, 45).tolist(),
animation_length=5,
)
animation
[21]:
Other mapping methods¶
In all of the previous examples, 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(filt, [[0, 1], [1, 3], [2]]) above.
We could also consider other mapping methods:
The \(k\)-chromatic subfiltration of
filtis the subfiltration spanned by all simplices having at most \(k\) colours. We can specify the inclusion of the \(k\)-chromatic subfiltration intofiltby using the classKChromaticInclusion.We could choose a family of subchromatic filtrations (each of which given as in the SubChromaticInclusion example above), and consider the union of their inclusions into
filt(the sub-chromatic gluing map). To do this we can useSubChromaticGluing.The subchromatic gluing map over the family of subfiltrations spanned by the various \(k\)-subsets of colours is specified with
KChromaticGluing.
Below we show an example using the KChromaticGluing by considering the union of three non-overlapping arcs that make up a circle.
These arcs are taken two at a time and glued together in the 2-chromatic gluing map, as shown below.
[22]:
# Generate three semi-circles
rng = np.random.default_rng(41)
num_points_arc = 40
mean = [0, 0]
cov = np.eye(2) * 0.005
arcs = tuple(
np.array(
[
[np.sin((t + i) * 2 * np.pi / 3), np.cos((t + i) * 2 * np.pi / 3)]
for t in rng.random(num_points_arc)
]
).T
+ rng.multivariate_normal(mean, cov, num_points_arc).T
for i in range(3)
)
points = np.concatenate(tuple(arcs), axis=1)
colours = np.array([i for i in chain.from_iterable([j] * num_points_arc for j in range(3))])
# Plot the point cloud
for i in range(3):
plt.scatter(arcs[i][0, :], arcs[i][1, :], s=10)
plt.gca().set_aspect("equal")
plt.gca().set_title("Point cloud")
# Compute and draw the bichromatic filtrations
fig, axes = plt.subplots(1, 3, figsize=(8, 18), constrained_layout=True)
cmap = np.asarray(plt.rcParams["axes.prop_cycle"].by_key()["color"])
kchromatic_example_filts = {}
for i, colour_pair in enumerate(combinations([0, 1, 2], 2)):
colour_pair = frozenset(colour_pair)
mask = [colour in colour_pair for colour in colours]
x = points[:, mask]
c = colours[mask]
filt = chromatic.delaunay_cech(x, c)
kchromatic_example_filts[colour_pair] = filt
draw_filtration(
filtration=filt,
points=x,
select=lambda simplex, _: simplex.filtration_value < 0.15,
ax=axes[i],
plot_colours=np.array([cmap[i] for i in colour_pair]),
)
axes[i].set_aspect("equal")
# Compute and draw the total filtration
total_filt = chromatic.delaunay_cech(points, colours)
ax = draw_filtration(
total_filt,
points,
select=lambda simplex, _: simplex.filtration_value < 0.15,
)
ax.set_title("Total filtration")
# Plot the 6-pack
gluing_morphism = KChromaticGluing(total_filt, k=2)
sixpack = gluing_morphism.sixpack(get_representatives=True)
fig, ax = plot_sixpack(sixpack)
fig.suptitle("6-pack from 2-chromatic gluing")
plt.show()
The the long-lived feature in the cokernel indicates the presence of a cycle in the total filtration that does not exist in any of the bichromatic subfiltrations. We can visualize a representative cycle for this feature, as in the code below. The required mapping is stored in the metadata variable in the 6-pack.
[23]:
assert sixpack.representatives is not None
assert sixpack.metadata is not None
cokernel = sixpack.cokernel
longest_interval_len = 0.0
longest_interval_birth_index = -1
for i in cokernel.pairings:
if cokernel.degrees[i] == 1:
birth, death = cokernel.interval(i)
if death is None:
longest_interval_birth_index = i
break
if death - birth > longest_interval_len:
longest_interval_len = death - birth
longest_interval_birth_index = i
interval_birth, interval_death = cokernel.interval(longest_interval_birth_index)
representative_indices = sixpack.representatives.cokernel[longest_interval_birth_index]
fig, ax = plt.subplots(figsize=(6, 6))
# Draw the filtration in the background
draw_filtration(
total_filt,
points,
select=lambda simplex, _: simplex.filtration_value <= interval_birth,
ax=ax,
edges_kwargs={"alpha": 0.15},
faces_kwargs={"alpha": 0.10},
)
# Overlay the representatives
draw_filtration(
total_filt,
points,
select=lambda _, index: index in representative_indices,
ax=ax,
edges_kwargs={"alpha": 1.0, "linewidth": 1.8, "zorder": 4},
)
ax.set_title(
f"KChromaticGluing: Cokernel cycle\n"
f"Cokernel interval [{interval_birth:.3f}, {interval_death:.3f})"
)
plt.show()