3.1. Tutorial: ECT for Embedded Graphs

This tutorial demonstrates how to use the ect package…

from ect import ECT, EmbeddedGraph

import numpy as np

3.1.1. Basic Usage

First, let’s create a simple graph”””

# Construct an example graph
# Note that this is the same graph that is returned by:
# G = create_example_graph()

G = EmbeddedGraph()

G.add_node("A", [1, 2])
G.add_node("B", [3, 4])
G.add_node("C", [5, 7])
G.add_node("D", [3, 6])
G.add_node("E", [4, 3])
G.add_node("F", [4, 5])

G.add_edge("A", "B")
G.add_edge("B", "C")
G.add_edge("B", "D")
G.add_edge("B", "E")
G.add_edge("C", "D")
G.add_edge("E", "F")

G.plot()

<Axes: >

The embedded graph class inherits from the networkx graph class with the additional attributes coord_matrix and coord_dict.

The coordinates of all vertices can be accessed using the coord_matrix attribute.

G.coord_matrix

array([[1., 2.],
       [3., 4.],
       [5., 7.],
       [3., 6.],
       [4., 3.],
       [4., 5.]])

It’s often useful to center the graph, so you can use the center_coordinates method shift the graph to have the average of the vertex coordinates be 0. Note that this does overwrite the coordinates of the points.

G.center_coordinates(center_type="mean")
print(G.coord_matrix)
G.plot()

[[-2.33333333 -2.5       ]
 [-0.33333333 -0.5       ]
 [ 1.66666667  2.5       ]
 [-0.33333333  1.5       ]
 [ 0.66666667 -1.5       ]
 [ 0.66666667  0.5       ]]

<Axes: >

To get a bounding radius we can use the get_bounding_radius method.

# This is actually getting the radius
r = G.get_bounding_radius()
print(f"The radius of bounding circle centered at the origin is {r}")

# plotting the graph with it's bounding circle of radius r.
G.plot(bounding_circle=True)

The radius of bounding circle centered at the origin is 3.2015621187164243

<Axes: >

We can also rescale our graph to have unit radius using scale_coordinates

G.scale_coordinates(radius=1)
G.plot(bounding_circle=True)

r = G.get_bounding_radius()
print(f"The radius of bounding circle centered at the origin is {r}")

The radius of bounding circle centered at the origin is 0.9362075413977773

myect = ECT(num_dirs=16, num_thresh=20)

# The ECT object will automatically create directions when needed
print(f"Number of directions: {myect.num_dirs}")
print(f"Number of thresholds: {myect.num_thresh}")

Number of directions: 16
Number of thresholds: 20

We can override the bounding radius as follows. Note that some methods will automatically use the bounding radius of the input G if not already set. I’m choosing the radius to be a bit bigger than the bounding radius of G to make some better pictures.

custom_bound_radius = 1.2 * G.get_bounding_radius()
result = myect.calculate(G, override_bound_radius=custom_bound_radius)

print(f"Thresholds chosen are: {myect.thresholds}")

Thresholds chosen are: [-1.12344905 -1.00519125 -0.88693346 -0.76867567 -0.65041787 -0.53216008
 -0.41390228 -0.29564449 -0.17738669 -0.0591289   0.0591289   0.17738669
  0.29564449  0.41390228  0.53216008  0.65041787  0.76867567  0.88693346
  1.00519125  1.12344905]

If we want the Euler characteristic curve for a fixed direction, we use the calculate function with a specific angle. This returns an ECTResult object containing the computed values.

result = myect.calculate(G, theta=np.pi / 2)
print(f"ECT values for direction pi/2: {result[0]}")

ECT values for direction pi/2: [0 0 0 0 1 1 2 2 2 1 1 1 1 1 1 1 0 0 0 0]

To calculate the full ECT, we call the calculate method without specifying theta. The result returns the ECT matrix and associated metadata.

result = myect.calculate(G)

print(f"ECT matrix shape: {result.shape}")
print(f"Number of directions: {myect.num_dirs}")
print(f"Number of thresholds: {myect.num_thresh}")

# We can plot the result matrix
result.plot()

ECT matrix shape: (16, 20)
Number of directions: 16
Number of thresholds: 20

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

## SECT

The Smooth Euler Characteristic Transform (SECT) can be calculated from the ECT. Fix a radius \(R\) bounding the graph. The average ECT in a direction \(\omega\) defined on function values \([-R,R]\) is given by

\[\overline{\text{ECT}_\omega} = \frac{1}{2R} \int_{t = -R}^{R} \chi(g_\omega^{-1}(-\infty,t]) \; dt.\]

Then the SECT is defined by

$$

$$

The SECT can be computed from the ECT result:

sect = result.smooth()

sect.plot()

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

## ECT for higher dimensions

The ECT class can also be used for higher dimensional graphs.

# generate a 3d graph
list_3d = [
    ("A", [0, 0, 0]),
    ("B", [1, 0, 0]),
    ("C", [0, 1, 0]),
    ("D", [1, 1, 0]),
    ("E", [0, 0, 1]),
    ("F", [1, 0, 1]),
    ("G", [0, 1, 1]),
    ("H", [1, 1, 1]),
]
graph_3d = EmbeddedGraph()
graph_3d.add_nodes_from(list_3d)

# add edges
graph_3d.add_edge("A", "B")
graph_3d.add_edge("B", "C")
graph_3d.add_edge("C", "D")
graph_3d.add_edge("D", "E")
graph_3d.add_edge("E", "F")
graph_3d.add_edge("F", "G")
graph_3d.add_edge("G", "H")
graph_3d.add_edge("H", "A")

# plot the graph
graph_3d.plot()

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

lets center the graph

graph_3d.center_coordinates(center_type="mean")
graph_3d.plot()

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

Now we can compute the ECT of the 3d graph.

ect_3d = ECT(num_dirs=16, num_thresh=20)
result_3d = ect_3d.calculate(graph_3d)
result_3d.plot()

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

Note that the each of the directions are appended in a list for the ECT result, so we won’t see the same periodic behavior as in the 2d case.

ECT results inherit from ndarrays but they store the associated directions and thresholds.

result_3d.directions.vectors

array([[ 0.48609512, -0.71244342, -0.50609872],
       [ 0.00364975, -0.82190818,  0.56960831],
       [-0.77710663, -0.23260984,  0.5848059 ],
       [ 0.75362655, -0.61431121, -0.23381352],
       [-0.18133219, -0.98270096, -0.03764916],
       [-0.92728701,  0.26549447,  0.26391569],
       [-0.81397529,  0.14802792,  0.56172232],
       [ 0.54073592, -0.29337918,  0.78837385],
       [ 0.48363696,  0.74722911, -0.45578937],
       [ 0.99631892,  0.01735462,  0.08394891],
       [ 0.74179171, -0.62822234, -0.23469503],
       [ 0.09906268, -0.01998319, -0.99488052],
       [ 0.57301878,  0.71110718, -0.40740159],
       [ 0.93645869, -0.11437729, -0.33160663],
       [-0.19170715,  0.53609831, -0.82209913],
       [ 0.69801901,  0.38101835, -0.6062957 ]])

result_3d.thresholds

array([-0.8660254 , -0.77486483, -0.68370427, -0.5925437 , -0.50138313,
       -0.41022256, -0.31906199, -0.22790142, -0.13674085, -0.04558028,
        0.04558028,  0.13674085,  0.22790142,  0.31906199,  0.41022256,
        0.50138313,  0.5925437 ,  0.68370427,  0.77486483,  0.8660254 ])

We can also define custom directions and thresholds for the ECT in case we need finer control. We use random sampling from the unit sphere for the directions and cosine sampling for the thresholds.

from ect import Directions

directions = Directions.random(16, dim=3)
thresholds = np.cos(np.linspace(0, np.pi, 20))
ect_3d = ECT(directions=directions, thresholds=thresholds)
result_3d = ect_3d.calculate(graph_3d)
result_3d.plot()

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

# another example sampling from the half sphere
sample_size = 100
theta = np.linspace(0, np.pi / 2, sample_size)  # Only go up to pi/2 for half sphere
phi = np.linspace(0, 2 * np.pi, sample_size)
theta, phi = np.meshgrid(theta, phi)

# Flatten the meshgrid arrays and create vectors
half_sphere_vectors = np.column_stack(
    [
        np.sin(theta.flatten()) * np.cos(phi.flatten()),
        np.sin(theta.flatten()) * np.sin(phi.flatten()),
        np.cos(theta.flatten()),
    ]
)

# Normalize the vectors
half_sphere_vectors = half_sphere_vectors / np.linalg.norm(
    half_sphere_vectors, axis=1, keepdims=True
)

directions = Directions.from_vectors(half_sphere_vectors)
print(f"Number of direction vectors: {len(directions)}")
ect_3d = ECT(directions=directions, num_thresh=20)  # Reduced number of thresholds
result_3d = ect_3d.calculate(graph_3d)
result_3d.plot()

Number of direction vectors: 10000

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