Intersecter Usage

We show how Intersecter works in the context of several examples, starting with the conformation space of Cyclooctane.

import numpy as np
import plotly.io as pio
from sklearn.manifold import Isomap

from dimmer.datasets import cyclooctane
from dimmer.intersecter import Intersecter
from dimmer.visualization import construct_scatter_with_labels

pio.renderers.default = "plotly_mimetype+notebook"

Example: Cyclooctane

The first example is the conformation space of cyclooctane. It can be loaded from the helper function dimmer.datasets.cyclooctane. It consists of 6040 points in 24-dimensional Euclidean space.

cyclopoints = cyclooctane()

print(cyclopoints.shape)

(6040, 24)

inter = Intersecter(radii=(0.4, 0.7), max_dim=2, n_jobs=3)

The Intersecter constructor takes seven (optional) parameters:

radii : a tuple of the inner and outer radii of each annular neighborhood,
neighbors : a tuple of the number of inner and outer neighbors,
max_dim : the maximum intrinsic dimension to check, and
n_jobs : the number of jobs to use for the computation.
collapse_edges : a boolean flag determining if the collapse_edges algorithm should be run in giotto-tda.
threshold : optional float used to determine when a persistent feature is relevant. The default behavior is the difference outer_rad - inner_rad.
outlier_label : the number to assign to every outlier.

The computation can be performed on a full data set X, or on a subset of it, specified by subset, a list of indices corresponding to X. Thus, subset is a subset of range(len(X)). If subset is not specified, the computation is performed on all of X.

The radii parameters can be floats or array-like objects, specifying the radii to use at each point of subset.

To compute all the local neighborhoods, call the fit method on the dataset.

inter.fit(X=cyclopoints)

Intersecter instance with radii=(0.4, 0.7), neighbors=None, and max_dim=2.

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

The intersection algorithm is implemented in the predict method.

strat_indices = inter.predict(X=cyclopoints)

The output is an with 0’s for manifold points and 1’s for intersection points.

np.count_nonzero(strat_indices)

The example of cyclooctane is nice because Isomap provides a convenient projection into 3-dimensions where we can see the \(S^2\) intersecting the Klein bottle in some compressed fashion.

iso = Isomap(n_components=3)
cyc_iso = iso.fit_transform(X=cyclopoints)

Let’s plot the isomap projection and color each point with its index as determined by Intersecter.

cyc_fig = construct_scatter_with_labels(X=cyc_iso, values=strat_indices)
cyc_fig.show()

Note that in the middle of this figure, we see two cones coming together, like the middle of an hour-glass. None of these points have been labelled as intersection points. This ‘intersection’ is an artifact of the isomap projection–it does not occur in the ambient 24-dimensional space. Contrast this with the two parallel circles, which do come from genuine intersections of the data set.

Example: Line and plane in \(\mathbb{R}^3\)

The second example is a line intersecting a plane in \(\mathbb{R}^3\). We construct the union of these two objects and run Intersecter on them.

rng = np.random.default_rng(42)

line = [(2 * rng.random() - 1, 0, 0) for _ in range(800)]
plane = [(0, 2 * rng.random() - 1, 2 * rng.random() - 1) for _ in range(700)]
data = np.concatenate((line, plane), axis=0)

Now data consists of points sampled from \([-1,1) \times \{(0,0)\}\) (the line) and \(\{0\} \times [-1,1)^2\) (the plane).

inter = Intersecter(radii=(0.3, 0.5), max_dim=2, n_jobs=3)

inter_ids = inter.fit_predict(X=data)

fig = construct_scatter_with_labels(X=data, values=inter_ids)
fig.show()

Points on the plane (except for a handful near the corners) and points on the line far from the intersection have appropriately connected annular neighborhoods, and therefore Intersecter labels them as manifold points. Near the intersection, points have neighborhoods which look roughly like the union two line segments (from the line) and a circle (from the plane). This has non-trivial \(H_0\) and \(H_1\), so these points are classified as intersection points. The number of intersection points is determined by the radii parameter, so this can be tuned to meet domain specific needs.

Example: Two planes in \(\mathbb{R}^3\)

Finally, we have two planes intersecting in \(\mathbb{R}^3\).

plane1 = [
    (
        2 * rng.random() - 1,
        2 * rng.random() - 1,
        0,
    )  # 5000
    for _ in range(700)
]
plane2 = [
    (0,) + (2 * rng.random() - 1, 2 * rng.random() - 1) for _ in range(700)
]  # 5000

data2 = np.concatenate((plane1, plane2), axis=0)

We instantiate Intersecter with the neighbors parameter instead of radii. This specifies the exact number of points to be used in each annular neighborhood. In the following, the annular neighborhood of each point consist of its 40th to 120th nearest neighbors. The radial threshold for determining relevance is computed by taking distance to the 120th nearest neighbor and subtracting the distance to the 40th nearest neighbor.

inter_planes = Intersecter(neighbors=(30, 80), max_dim=2, n_jobs=5)
inter_planes_ids = inter_planes.fit_predict(X=data2)

fig = construct_scatter_with_labels(X=data2, values=inter_planes_ids)
fig.show()

Once again, points near the intersection of the two planes are determined to be intersection points, while those far away are labelled as non-intersection (i.e. manifold) points.