Lecture 2: Simplicial Complexes

Lecture 2: Simplicial Complexes#

In this notebook we will familiarize ourselves a bit more with the concepts from the lecture. We will work with GUDHI, a Python library for topological data analysis (TDA) that lets us handle simplicial complexes.

import itertools, time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Polygon  # for drawing balls and filled simplices
from scipy.spatial import Voronoi, Delaunay, voronoi_plot_2d  # Voronoi diagrams and Delaunay triangulations
import gudhi

plt.rcParams["figure.dpi"] = 110

# Plotting helper used throughout the notebook: takes coordinates and anything that
# yields simplices (a SimplexTree, a list of tuples), and optionally draws a ball of
# radius `r` around each point.
def draw(points, complex_, ax=None, r=None, title=None, ms=30, lw=1.2):
    simplices = [s for s, _ in complex_.get_simplices()] \
        if hasattr(complex_, "get_simplices") else [list(s) for s in complex_]
    points = np.asarray(points, dtype=float)
    if ax is None:
        _, ax = plt.subplots(figsize=(4, 4))
    if r is not None:
        for p in points:
            ax.add_patch(Circle(p, r, color="#C39BD3", alpha=0.3, lw=0))
    for s in simplices:
        if len(s) == 2:
            ax.plot(*points[s].T, color="black", lw=lw, zorder=2)
        elif len(s) == 3:
            ax.add_patch(Polygon(points[s], color="#7FB3D5", alpha=0.6, zorder=1))
    ax.scatter(*points.T, s=ms, color="black", zorder=3)
    ax.set_aspect("equal"); ax.set_xticks([]); ax.set_yticks([])
    ax.set_title(title, fontsize=10)
    return ax

1. Simplicial complexes in GUDHI#

The central object for storing simplicial complexes is the SimplexTree: a simplicial complex where every simplex carries a filtration value (this can be seen as their insertion time, we will talk more about filtrations in upcoming lectures).

Inserting a simplex inserts all of its faces, so the library enforces that a SimplexTree always stores a valid simplicial complex. For example, if we add a filled triangle (a simplex with three elements), GUDHI will add the three vertices and the three edges automatically.

st = gudhi.SimplexTree()
st.insert([0, 1, 2]) # one filled triangle

for s, _ in st.get_simplices(): # the ignored second value is the filtration value
    print(s)

Seven simplices out of one insert: the triangle, its three edges, its three vertices. That is exactly \(K = \{\{v_0\},\{v_1\},\{v_2\},\{v_0v_1\},\{v_1v_2\},\{v_0v_2\},\{v_0v_1v_2\}\}\) from the slides, and you only had to type the top-dimensional piece.

The second value printed for each simplex is its filtration value, i.e. its insertion time. Here it is 0.0 for all of them, since GUDHI created every one of these simplices at once, at instant 0, when we inserted the filled triangle.

SimplexTree.insert expects a sequence of integers: vertices have to be integers, not arbitrary hashables. If you want a geometric complex, keep the coordinates separately, say in a dict mapping each vertex to a point, and only use that dict when you draw the complex or reason about it geometrically. The SimplexTree itself only ever sees the abstract complex: it has no idea where the vertices sit in space, so it will not enforce anything about the embedding, for example that edges do not cross.

coords = {0: (0.0, 0.0), 1: (1.0, 0.0), 2: (0.5, 0.9)} # vertex -> location, kept separately

empty_triangle = gudhi.SimplexTree()
for e in [(0,), (1,), (2,), (0, 1), (1, 2), (0, 2)]: # add every simplex we want, one by one
    empty_triangle.insert(list(e))

print("simplices  ", empty_triangle.num_simplices())
print("dimension  ", empty_triangle.dimension())
print("f-vector   ", empty_triangle.num_simplices_by_dimension())

Alternatively, we could have inserted the filled triangle and then removed only its top simplex, leaving its faces behind: remove_maximal_simplex deletes a simplex that currently has no cofaces, which right after inserting is exactly the simplex we just added.

To really make the point that SimplexTree does not care about geometry, let’s break the coordinates first: collapse all three vertices onto the same point. As a geometric simplicial complex this would not be a triangle at all; the convex hull of three coincident points is just a dot. But the SimplexTree only stores the combinatorics, so it builds the exact same empty triangle regardless.

coords = {0: (0.5, 0.5), 1: (0.5, 0.5), 2: (0.5, 0.5)}   # degenerate: all three vertices coincide

also_empty_triangle = gudhi.SimplexTree()
also_empty_triangle.insert([0, 1, 2])
also_empty_triangle.remove_maximal_simplex([0, 1, 2])

print("simplices  ", also_empty_triangle.num_simplices())
print("dimension  ", also_empty_triangle.dimension())
print("f-vector   ", also_empty_triangle.num_simplices_by_dimension())

print("same complex as before:",
      sorted(s for s, _ in also_empty_triangle.get_simplices())
      == sorted(s for s, _ in empty_triangle.get_simplices()))

We just used num_simplices, dimension, and num_simplices_by_dimension to query a simplicial complex about its properties. SimplexTree has many more. See the full method list in the documentation.

How many simplices do we need to be able to have an \(N-1\)-dimensional simplex?#

Small exercise. If you insert a single simplex with \(N\) elements (an \((N-1)\)-simplex), how many simplices end up in the SimplexTree? Face containment means every non-empty subset of those \(N\) vertices has to be there too. Write num_simplices_in_simplex(N) returning the formula, then check it against GUDHI below.

def num_simplices_in_simplex(N):
    '''Number of simplices you get from inserting a single N-element simplex.'''
    return 0

for N in range(1, 8):
    st = gudhi.SimplexTree()
    st.insert(list(range(N)))
    got, want = st.num_simplices(), num_simplices_in_simplex(N)
    print(f"{'ok  ' if got == want else 'FAIL'} N={N} (a {N - 1}-simplex): "
          f"gudhi={got}  yours={want}")

Faces, boundary, interior, skeleton#

A few more SimplexTree queries, on an easy complex: two triangles glued along a shared edge.

  • get_boundaries gives the facets of a simplex, its proper codimension-1 faces;

  • get_skeleton(k) gives the \(k\)-skeleton, every simplex of dimension \(\le k\);

  • get_cofaces(simplex, 1) tells you what sits one dimension up. The shared edge has two cofaces, the outer edges has only one.

K = gudhi.SimplexTree()
for t in [(0, 1, 2), (1, 2, 3)]: # two filled triangles sharing the edge {1, 2}
    K.insert(list(t))

print("all simplices     ", [s for s, _ in K.get_simplices()])
print("facets of {0,1,2} ", [s for s, _ in K.get_boundaries([0, 1, 2])])
print("1-skeleton        ", [s for s, _ in K.get_skeleton(1)])
print("cofaces of {1,2}  ", [s for s, _ in K.get_cofaces([1, 2], 1)], "-> interior edge")
print("cofaces of {0,1}  ", [s for s, _ in K.get_cofaces([0, 1], 1)], "-> boundary edge")

2. Point clouds, distances, balls#

A point cloud is a finite collection of points in a metric space \((M, d)\). This section has a few helper functions to build the point clouds we will use for the rest of the notebook.

def circle_cloud(n=24, radius=1.0, noise=0.07, seed=0):
    g = np.random.default_rng(seed)
    t = g.uniform(0, 2 * np.pi, n)
    rad = radius + noise * g.standard_normal(n)
    return np.column_stack([rad * np.cos(t), rad * np.sin(t)])

g = np.random.default_rng(0)

# six small clouds with different shapes; the later comparisons run on all of
# them. "grid" is cocircular in groups of four, the degenerate case that makes
# Delaunay implementations sweat.
CLOUDS = {
    "circle":      circle_cloud(n=18, radius=1.0, noise=0.05, seed=1),
    "big circle":  circle_cloud(n=24, radius=3.2, noise=0.22, seed=7),
    "two circles": np.vstack([circle_cloud(9, 0.6, 0.04, seed=2),
                              circle_cloud(9, 0.6, 0.04, seed=3) + [1.7, 0]]),
    "uniform":     g.uniform(0, 1, (18, 2)),
    "grid":        np.array([[i, j] for i in range(4) for j in range(4)]) / 3,
    "clusters":    np.vstack([g.normal(c, 0.12, (6, 2))
                              for c in [(0, 0), (1, 0.2), (0.4, 1)]]),
}

small = circle_cloud(n=12, radius=1.0, noise=0.05, seed=5)   # used in sections 5 and 8

Let’s visualize the clouds we’ve just created.

fig, axes = plt.subplots(1, len(CLOUDS), figsize=(3 * len(CLOUDS), 3.2))
for ax, (name, X) in zip(axes, CLOUDS.items()):
    ax.scatter(*X.T, s=18, color="black")
    ax.set_title(f"{name} (n={len(X)})", fontsize=9)
    ax.set_aspect("equal"); ax.set_xticks([]); ax.set_yticks([])
plt.tight_layout(); plt.show()

3. Čech Complex#

Given a finite collection of sets, the nerve records which subfamilies have a common point: $\(\mathrm{Nrv}(\mathcal{U}) = \Big\{X \subseteq \mathcal{U} \ \Big|\ \bigcap_{U \in X} U \neq \emptyset\Big\}.\)$ This is automatically a simplicial complex: intersecting fewer sets can only make the intersection bigger, so faces come for free.

The Čech complex is the nerve of the cover by balls, $\(\check{C}^{\,r}(P) = \Big\{\alpha \subseteq P \ \Big|\ \bigcap_{x \in \alpha} B(x,r) \neq \emptyset\Big\} = \mathrm{Nrv}\big(\{B(x,r)\}_{x \in P}\big),\)$ and by the Nerve Lemma (since balls are convex every intersection of them is contractible) it has the homotopy type of the union of balls.

def enclosing_radius(pts):
    """Radius of the smallest ball containing 2 or 3 points in the plane."""
    if len(pts) == 2: # two points: the ball is centered at the midpoint of the segment connecting them
        return np.linalg.norm(pts[0] - pts[1]) / 2
    a, b, c = pts
    sides = np.array([np.linalg.norm(b - c), np.linalg.norm(a - c),
                      np.linalg.norm(a - b)])
    if sides.max() ** 2 >= (sides ** 2).sum() - sides.max() ** 2:  # the thre points make an obtuse triangle: the longest edge decides
        return sides.max() / 2
    # otherwise, the points make an acute triangle, and the circumradius is the right answer
    # https://en.wikipedia.org/wiki/Circumcircle 
    area = abs((b - a)[0] * (c - a)[1] - (b - a)[1] * (c - a)[0]) / 2
    return sides.prod() / (4 * area)

two = np.array([[0.0, 0.0], [1.0, 0.0]])
far = np.array([[0.0, 0.0], [10.0, 0.0], [0.0, 10.0]])
unit_tri = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3) / 2]], dtype=float)

Exercise. Write cech(points, r, max_dim=2) returning a SimplexTree whose vertices are the row indices of points. Straight from the definition: for every subset of size up to max_dim + 1, decide whether the balls have a common point, and insert the ones that do. enclosing_radius handles 2 and 3 points, which is why max_dim stops at 2.

def my_cech(points, r, max_dim=2):
    '''Cech complex C^r(points) as a SimplexTree.'''
    # tip: use itertools.combinations to generate all subsets of size 1, 2, ..., max_dim + 1
    # tip: use the enclosing_radius function to check if a subset of points is contained in a ball of radius r
    raise NotImplementedError

Five checks with answers you can do on paper: two points span an edge exactly when \(r\) reaches half their distance, an equilateral triangle of side 1 fills in at its circumradius \(1/\sqrt{3} \approx 0.577\), and points far apart stay isolated.

def check_cech(fn):
    cases = [
        ("two points, r just under d/2", two, 0.49,
         [(0,), (1,)]),
        ("two points, r just over d/2",  two, 0.51,
         [(0,), (0, 1), (1,)]),
        ("triangle, r under 1/sqrt(3)",  unit_tri, 0.55,
         [(0,), (0, 1), (0, 2), (1,), (1, 2), (2,)]),
        ("triangle, r over 1/sqrt(3)",   unit_tri, 0.58,
         [(0,), (0, 1), (0, 1, 2), (0, 2), (1,), (1, 2), (2,)]),
        ("far apart, r = 1",             far, 1.00,
         [(0,), (1,), (2,)]),
    ]
    for name, pts, r, want in cases:
        tree = fn(pts, r)   # bind before iterating get_simplices(), or GUDHI can segfault on the temporary
        got = sorted(tuple(s) for s, _ in tree.get_simplices())
        print(f"{'ok  ' if got == want else 'FAIL'} {name:30s} simplices = {got}")


check_cech(cech)
# check_cech(my_cech)  # uncomment to test your implementation

And on the point cloud, as \(r\) grows.

POINT_CLOUDS = {"two points": two, "unit triangle": unit_tri, "far apart": far,
                 "small (circle)": small, **CLOUDS}

cloud = POINT_CLOUDS["big circle"]   # experiment: try any other key above

fig, axes = plt.subplots(1, 4, figsize=(14, 3.8))
for ax, r in zip(axes, [0.4, 0.7, 1.0, 1.5]):
    st = cech(cloud, r)
    draw(cloud, st, ax=ax, r=r, ms=18, lw=0.9,
         title=f"r = {r}\n num simplices by dimension: {st.num_simplices_by_dimension()}")
plt.tight_layout(); plt.show()

Remember: the Čech complex is an abstract simplicial complex. Drawing it by putting each vertex at its own point, as above, does not generally give a geometric simplicial complex.

4. Vietoris–Rips#

\[\mathrm{VR}^r(P) = \{\alpha \subseteq P \mid d(x_i,x_j) \le 2r \ \text{ for all } x_i,x_j \in \alpha\}\]

Čech asks the points to be jointly close; Vietoris-Rips only asks them to be pairwise close. This key difference makes Rips much cheaper: we don’t need the underlying geometry/coordinates, just the distance matrix (the distance between any two points in the cloud).

GUDHI builds it directly with RipsComplex, which takes max_edge_length which acts as the diameter, while the slides use a radius. So max_edge_length = 2 * r; getting this wrong shifts every picture by a factor of two.

We also cap the dimension with max_dim: without it, the complex would keep growing to every clique in the graph, up to \(2^n - 1\) simplices for \(n\) points.

def vietoris_rips(points, r, max_dim=2):
    return gudhi.RipsComplex(points=points, max_edge_length=2 * r).create_simplex_tree(max_dimension=max_dim)

The equilateral triangle#

Let’s look at an example where the two are actually different. Three points at mutual distance 1. Rips fills the triangle as soon as all three edges exist, at \(2r = 1\); Čech waits until the three balls actually share a point, at \(r = 1/\sqrt{3}\). In between, Rips claims a filled triangle while the union of balls still has a curved gap.

for r in [0.45, 0.5, 0.55, 0.6]:
    print(f"r={r}   Rips num simplices by dimension: {vietoris_rips(unit_tri, r).num_simplices_by_dimension()};"
          f"   Cech num simplices by dimension: {cech(unit_tri, r).num_simplices_by_dimension()}.")

# plot the Vietoris-Rips and Cech complexes for the unit triangle at a specific radius
r = 0.52 # play around with this value to see how the two complexes differ
fig, axes = plt.subplots(1, 2, figsize=(8, 3.8))
draw(unit_tri, vietoris_rips(unit_tri, r), ax=axes[0], r=r, title=f"VR, r={r}")
draw(unit_tri, cech(unit_tri, r), ax=axes[1], r=r, title=f"Cech, r={r}")
plt.tight_layout(); plt.show()

The Rips–Čech lemma#

So how wrong can Rips be? As we saw in the lecture we know that the follwoing chain of inclusions hold:

\[\check{C}^{\,r}(P) \ \subseteq\ \mathrm{VR}^r(P) \ \subseteq\ \check{C}^{\,2r}(P).\]

Below we test it empirically

lemma_cloud = np.random.default_rng(0).uniform(0, 1, (5, 2))

r = 0.14 # play around with this value to see how the three complexes differ

# build the complexes
C_tree, VR_tree, C2_tree = cech(lemma_cloud, r), vietoris_rips(lemma_cloud, r), cech(lemma_cloud, 2 * r)

# print the simplices
print("Cech(r)   ", sorted(tuple(sorted(s)) for s, _ in C_tree.get_simplices()))
print("VR(r)     ", sorted(tuple(sorted(s)) for s, _ in VR_tree.get_simplices()))
print("Cech(2r)  ", sorted(tuple(sorted(s)) for s, _ in C2_tree.get_simplices()))

# visually compare the three complexes
fig, axes = plt.subplots(1, 3, figsize=(11, 3.8))
draw(lemma_cloud, C_tree, ax=axes[0], r=r, title=f"Cech(r), r={r}")
draw(lemma_cloud, VR_tree, ax=axes[1], r=r, title=f"VR(r), r={r}")
draw(lemma_cloud, C2_tree, ax=axes[2], r=2 * r, title=f"Cech(2r), r={2 * r}")
plt.tight_layout(); plt.show()

5. Voronoi, Delaunay, Alpha#

The Voronoi cell of \(u\) is the region of the plane closer to \(u\) than to any other sample point. The Delaunay complex is the nerve of the Voronoi diagram, and it has no parameter: it is fixed by the point cloud.

Q = circle_cloud(n=18, radius=2.0, noise=0.25, seed=3)
# Q = unit_tri # a much simpler example, feel free to try other point clouds from the POINT_CLOUDS dictionary above
vor, dela = Voronoi(Q), Delaunay(Q)

fig, axes = plt.subplots(1, 2, figsize=(9, 4.4))
voronoi_plot_2d(vor, ax=axes[0], show_vertices=False, line_colors="#2E86C1", point_size=12)
axes[0].set_title("Voronoi diagram")

draw(Q, [tuple(t) for t in dela.simplices], ax=axes[1], ms=12, title="Delaunay")
axes[1].triplot(*Q.T, dela.simplices, color="black", lw=0.9, zorder=2)   # crisp black edges over the fill
voronoi_plot_2d(vor, ax=axes[1], show_vertices=False, line_colors="#2E86C1", point_size=12)

for ax in axes:
    ax.set_aspect("equal"); ax.set_xticks([]); ax.set_yticks([])
    ax.set_xlim(-3.5, 3.5); ax.set_ylim(-3.5, 3.5)
plt.tight_layout(); plt.show()

The alpha complex shrinks each ball to its own Voronoi cell, \(D_p^{\alpha} = B(p,r) \cap V_p\), and takes the nerve of that. It is a subcomplex of Delaunay that grows with \(r\) (for large enough \(r\) it becomes Delaunay), and, like Čech but unlike Rips, it has the homotopy type of the union of balls.

AlphaComplex returns squared filtration values by default; output_squared_values=False gives plain radii, matching the \(r\) on the slides. GUDHI may reorder the points, so read the coordinates back with get_point.

ac = gudhi.AlphaComplex(points=Q)
alpha = ac.create_simplex_tree(output_squared_values=False)
A = np.array([ac.get_point(i) for i in range(len(Q))])

fig, axes = plt.subplots(1, 4, figsize=(14, 3.8))
for ax, r in zip(axes, [0.35, 0.55, 0.8, 1.2]):
    st = alpha.copy()
    st.prune_above_filtration(r)
    draw(A, st, ax=ax, r=r, ms=18, lw=0.9,
         title=f"Alpha, r={r}\n{st.num_simplices_by_dimension()}")
plt.tight_layout(); plt.show()

At a shared radius, Rips, Čech, and Alpha can have very different sizes. Since \(\mathrm{Alpha}(r) \subseteq \check{C}^{\,r} \subseteq \mathrm{VR}^r\), Rips is never smaller: the extra simplices are triangles Rips invents, points pairwise within \(2r\) whose balls share no common point.

def simplices(st):
    return {frozenset(s) for s, _ in st.get_simplices()}


ac_small = gudhi.AlphaComplex(points=small)
alpha_small = ac_small.create_simplex_tree(output_squared_values=False)
S = np.array([ac_small.get_point(i) for i in range(len(small))])

r = 0.92
R_tree, C_tree = vietoris_rips(small, r), cech(small, r)
A_tree = alpha_small.copy()
A_tree.prune_above_filtration(r)

extra = [sorted(t) for t in simplices(R_tree) - simplices(C_tree) if len(t) == 3]
print(f"Rips {R_tree.num_simplices()}   Cech {C_tree.num_simplices()}   "
      f"Alpha {A_tree.num_simplices()}   extra triangles in Rips: {len(extra)}")

fig, axes = plt.subplots(1, 3, figsize=(11, 3.8))
for ax, name, pts, st in [(axes[0], "Rips", small, R_tree),
                          (axes[1], "Cech", small, C_tree),
                          (axes[2], "Alpha", S, A_tree)]:
    draw(pts, st, ax=ax, r=r, ms=22, lw=0.9, title=f"{name}, {st.num_simplices()} simplices")
for t in extra:
    axes[0].add_patch(Polygon(small[t], color="crimson", alpha=0.55, zorder=2))
plt.tight_layout(); plt.show()

Two more claims from the slides: \(\mathrm{Alpha}(r) \subseteq \mathrm{Delaunay}\), with equality once \(r\) is large enough, and \(\mathrm{Alpha}(r) \subseteq \check{C}^{\,r}(P)\). Both are set inclusions, so we can check them directly.

ac_ = gudhi.AlphaComplex(points=small)
tree = ac_.create_simplex_tree(output_squared_values=False)
coords_ = np.array([ac_.get_point(i) for i in range(len(small))])
delaunay = simplices(tree)

for r in [0.2, 0.5, 0.8, 1.5, 10.0]:
    a = tree.copy()
    a.prune_above_filtration(r)
    c = cech(coords_, r)
    print(f"r={r:5.2f}   Alpha < Delaunay: {simplices(a) <= delaunay}   "
          f"Alpha < Cech: {simplices(a) <= simplices(c)}   "
          f"Alpha == Delaunay: {simplices(a) == delaunay}")

Alpha sits between Delaunay and Čech, and is the smaller complex nearly everywhere: that is the practical argument for it. Once you can measure homotopy type (next lecture), the question becomes whether the smaller complex has lost anything.

6. What it costs#

With \(n\) points, the \(k\)-skeleton of a Rips complex can have \(\sum_{j \le k} \binom{n}{j+1}\) simplices, and the full complex \(2^n - 1\). A Delaunay or Alpha complex in the plane has \(O(n)\) simplices: a triangulation of \(n\) planar points has at most \(3n-6\) edges and \(2n-5\) triangles, however the points are arranged.

Fix the radius, grow the cloud, and watch. These are GUDHI’s implementations; the by-hand versions above enumerate subsets and would be far slower.

def timed(build, pts):
    t0 = time.perf_counter()
    st = build(pts)
    return st.num_simplices(), time.perf_counter() - t0


rows = []
for n in [50, 100, 200, 400, 800]:
    pts = circle_cloud(n=n, radius=1.0, noise=0.08, seed=n)
    s2, t2 = timed(lambda p: vietoris_rips(p, 0.2, 2), pts)
    sa, ta = timed(lambda p: gudhi.AlphaComplex(points=p).create_simplex_tree(), pts)
    if n <= 400:                       # past this the 3-skeleton eats the RAM
        s3, t3 = timed(lambda p: vietoris_rips(p, 0.2, 3), pts)
    else:
        s3 = t3 = np.nan
    rows.append([n, s2, t2, s3, t3, sa, ta])
    d3 = "skipped" if np.isnan(s3) else f"{s3:,.0f} ({t3*1e3:.1f} ms)"
    print(f"n={n:4d}   rips dim2 {s2:>9,.0f} ({t2*1e3:5.1f} ms)   "
          f"rips dim3 {d3:>22}   alpha {sa:>6,.0f} ({ta*1e3:4.1f} ms)")

rows = np.array(rows, dtype=float)

n = rows[:, 0]
fig, axes = plt.subplots(1, 2, figsize=(11, 4))

axes[0].loglog(n, rows[:, 1], "o-", label="Rips 2-skeleton")
axes[0].loglog(n, rows[:, 3], "s-", label="Rips 3-skeleton")
axes[0].loglog(n, rows[:, 5], "^-", label="Alpha")
axes[0].loglog(n, rows[0, 1] * (n / n[0]) ** 3, "k:", lw=1, label="$n^3$")
axes[0].loglog(n, rows[0, 5] * (n / n[0]), "k--", lw=1, label="$n$")
axes[0].set_xticks(n)
axes[0].set_xticklabels([f"{int(v)}" for v in n])
axes[0].xaxis.set_minor_formatter(plt.NullFormatter())
axes[0].set_xlabel("points"); axes[0].set_ylabel("simplices")
axes[0].set_title("log-log")
axes[0].legend(fontsize=8); axes[0].grid(alpha=.3)

axes[1].plot(n, rows[:, 1], "o-", label="Rips 2-skeleton")
axes[1].plot(n, rows[:, 3], "s-", label="Rips 3-skeleton")
axes[1].plot(n, rows[:, 5], "^-", label="Alpha")
axes[1].set_xlabel("points"); axes[1].set_ylabel("simplices")
axes[1].set_title("linear, same data")
axes[1].legend(fontsize=8); axes[1].grid(alpha=.3)

plt.tight_layout(); plt.show()

7. Other complexes, if you need them#

The lecture covered Čech, Rips, Delaunay and Alpha. GUDHI ships several more, mostly to dodge the cost problem above or to handle data that is not a point cloud.

what it is for

docs

Witness

a few landmarks carry the complex, the rest only vote

manual

Cubical

images and grids, where the natural cells are pixels

manual

Delaunay–Čech

Čech filtration values on the Delaunay complex

manual

Weighted Rips / DTM-Rips

Rips that is robust to outliers

manual

Tangential

manifold data in high ambient dimension

manual

Cover complexes / Mapper

graphs and visual summaries from a cover

manual

8. Optional exercise: The Nerve of a cover#

The nerve is not only for balls: take any cover (a list of sets, each one containing whatever you like: letters, words, numbers) and build its nerve. The nerve itself is an abstract simplicial complex that does not care what is inside the sets.

Write cover_nerve(cover, max_dim=2). The vertices of the result are positions in the cover list, one per set.

def cover_nerve(cover, max_dim=2):
    '''Nerve of a list of sets. Vertices are indices into `cover`.'''
    raise NotImplementedError

COVER_CASES = {
    "three sets, pairwise only": (
        [{"cat", "dog"}, {"dog", "bird"}, {"bird", "cat"}],
        [(0,), (0, 1), (0, 2), (1,), (1, 2), (2,)]),
    "three sets with a common point": (
        [{"cat", "dog", "bird"}, {"dog", "bird", "fish"}, {"bird", "fish", "frog"}],
        [(0,), (0, 1), (0, 1, 2), (0, 2), (1,), (1, 2), (2,)]),
    "disjoint": (
        [{"a"}, {"b"}, {"c"}],
        [(0,), (1,), (2,)]),
    "one set inside another": (
        [{"x", "y", "z"}, {"y"}],
        [(0,), (0, 1), (1,)]),
}


def check_cover_nerve(fn):
    for name, (cover, expected) in COVER_CASES.items():
        st = fn(cover)
        got = sorted(tuple(s) for s, _ in st.get_simplices())
        ok = got == sorted(expected)
        print(f"{'ok  ' if ok else 'FAIL'} {name}")
        if not ok:
            print("     got     ", got)
            print("     expected", sorted(expected))


check_cover_nerve(cover_nerve)

Where to read more#