3.2. Tutorial: ECT for CW complexes

This tutorial walks you through how to build a CW complex with the EmbeddedCW class, and then use the ECT class to compute the Euler characteristic transform

from ect import ECT, EmbeddedCW
from ect.utils.examples import create_example_cw
import numpy as np

The CW complex is the same as the EmbeddedGraph class with that additional ability to add faces. Faces are added by passing in a list of vertices. Note that we are generally assuming that these vertices follow around an empty region (as in, no other vertex is in the interior) in the graph bounded by the vertices, and further that all edges are already included in the graph. However the class does not yet check for this so you need to be careful!

K = EmbeddedCW()

# Add vertices with coordinates
K.add_node("A", [0, 0])
K.add_node("B", [1, 0])
K.add_node("C", [1, 1])
K.add_node("D", [0, 1])

# Add edges to form a square
K.add_edges_from([("A", "B"), ("B", "C"), ("C", "D"), ("D", "A")])

# Add the square face
K.add_face(["A", "B", "C", "D"])

K.center_coordinates()
K.plot()

<Axes: >

Just to have something a bit more interesting, let’s make a more complicated example that’s built into the class.

K = create_example_cw()
K.plot(bounding_circle=True)

<Axes: >

As with the EmbeddedGraph class, we can initialize the ECT class by deciding how many directions and how many thresholds to use.

ect = ECT(num_dirs=100, num_thresh=80)

Then we can compute the ECC for a single direction. In this case, the \(x\)-axis will be computed for the num_thresh=80 stopping points in the interval \([-1.2r,1.2r]\) where \(r\) is the minimum bounding radius for the input complex.

override_bound_radius = 1.2 * K.get_bounding_radius()
result = ect.calculate(K, theta=0, override_bound_radius=override_bound_radius)
result.plot();

/Users/lizliz/Library/CloudStorage/Dropbox/Math/Code/ect/src/ect/ect_graph.py:211: NumbaPerformanceWarning: 
The keyword argument 'parallel=True' was specified but no transformation for parallel execution was possible.

To find out why, try turning on parallel diagnostics, see https://numba.readthedocs.io/en/stable/user/parallel.html#diagnostics for help.

File "src/ect/ect_graph.py", line 216:
    @njit(parallel=True, fastmath=True)
    def shape_descriptor(simplex_counts_list):
    ^

  result[i, j] = shape_descriptor(simplex_counts_list)
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

But of course it’s easier to see this in a plot. This command calculates the ECC and immediately plots it.

Similarly, we can compute the ECT and return the matrix. We make sure to internally set the bounding radius to use to control the \(y\) axis of the plot.

result = ect.calculate(K, override_bound_radius=override_bound_radius)
result.plot()

<Axes: xlabel='Direction $\\omega$ (radians)', ylabel='Threshold $a$'>

We can also look at the Smooth ECT:

# Calculate SECT and plot
smooth = result.smooth()
smooth.plot()

<Axes: xlabel='Direction $\\omega$ (radians)', ylabel='Threshold $a$'>

We can also compute the ECT in 3D.

import numpy as np

vertices = [
    (letter, coordinates) for letter, coordinates in zip("abcde", np.random.randn(5, 3))
]
edges = [("a", "b"), ("a", "c"), ("a", "d"), ("b", "c"), ("b", "d"), ("c", "d")]
faces = [
    ("a", "b", "c"),
    ("a", "b", "d"),
    ("a", "c", "d"),
    ("b", "c", "d"),
    ("a", "b", "c", "d"),
]
K = EmbeddedCW()
K.add_nodes_from(vertices)
K.add_edges_from(edges)

K.add_faces_from(faces)
K.plot(bounding_circle=True)

<Axes3D: xlabel='X', ylabel='Y', zlabel='Z'>

ect = ECT(num_dirs=100, num_thresh=80)
result = ect.calculate(K)
result.plot()

<Axes: xlabel='Direction Index', ylabel='Threshold $a$'>