GAD Usage

In this notebook, we highlight some basic functionality of GAD.

The GAD constructor takes seven optional parameters:

radii : a tuple of the inner and outer radius of each annular neighborhood.
neighbors : a tuple of the inner and outer number of neighbors of each annular neighborhood.
max_dim : the maximum intrinsic dimension to check; defaults to 1.
n_jobs : the number of jobs to use for the computation; defaults to 1. This permits the persistent homology computation of each annular neighborhood to be parallelized.
collapse_edges : a boolean flag determining if the collapse_edges algorithm should be run in giotto-tda; defaults to False if max_dim < 4, otherwise defaults to True.
threshold : an optional float used to determine when a persistent feature is relevant. The default behavior is the difference radii[1] - radii[0] if radii is specified or computes this difference if neighbors is specified.
outlier_label : a number to assign to any outlier, defaults to -1.0.

Notes:

The first two parameters radii and neighbors determine the size of the local neighborhood and have the greatest impact on the results.

⚠️ Only one of radii or neighbors should be specified.

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.

tl;dr: All classes in dimmer are sklearn transformers. To compute all the local neighborhoods, use the fit method. To perform the dimension estimation, use predict. To perform stratification estimation, use predict_labels.

Example: Rectangle

The first example consists of a rectangle in the plane.

import numpy as np
import plotly.io as pio

from dimmer.gad import GAD
from dimmer.visualization import (
    construct_scatter_with_labels,
    plot_annular_neighborhood,
)

pio.renderers.default = "plotly_mimetype+notebook"

rng = np.random.default_rng(0)
x = rng.random(1000)
y = 0.5 * rng.random(1000)
data_rect = np.array(list(zip(x, y)))

rect_fig = construct_scatter_with_labels(X=data_rect)
rect_fig.update_scenes(xaxis_visible=False, yaxis_visible=False)
rect_fig.update_layout(yaxis=dict(scaleanchor="x"))
rect_fig.show()

We instantiate a GAD instance and perform the local neighborhood computations in fit.

gad_l = GAD(radii=(0.08, 0.12), max_dim=2, n_jobs=3)
gad_l.fit(X=data_rect)

GAD instance with radii=(0.08, 0.12), 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 iterative algorithm is implemented in the predict method:

id_labels = gad_l.predict(X=data_rect)

np.unique(id_labels).shape

(3,)

fig = construct_scatter_with_labels(
    X=data_rect, values=id_labels, label="Dimension {value} points"
)
fig.update_scenes(xaxis_visible=False, yaxis_visible=False)
fig.update_layout(yaxis=dict(scaleanchor="x"))

fig.show()

The plotting methods are explained in more detail in the visualization notebook.

Both computations can be performed with one call using the fit_predict method.

gad_ids = GAD(radii=(0.08, 0.12), max_dim=2, n_jobs=3).fit_predict(X=data_rect)

(gad_ids == id_labels).all()

True

The algorithm computes the (two-dimensional) interior of the rectangle, its (one-dimensional) boundary, and the (“zero-dimensional”) corners.

The subset parameter can be passed to fit to localize the computation of the annular neighborhoods of only those specified by subset.

subset = rng.choice(len(data_rect), len(data_rect) // 10, replace=False)

gad_ids_with_subset = GAD(radii=(0.08, 0.12), max_dim=2, n_jobs=3).fit_predict(
    X=data_rect, subset=subset
)

fig = construct_scatter_with_labels(
    X=data_rect, values=gad_ids_with_subset, label="Dimension {value} points"
)
fig.update_scenes(xaxis_visible=False, yaxis_visible=False)
fig.update_layout(yaxis=dict(scaleanchor="x"))

fig.show()

We see the subset should be sampled respecting the stratification since the maximal stratum dominates the sampling.

Detecting stratifications

As mentioned above, GAD is a stratification finding tool, not a dimension estimation tool. Although the previous rectangle example highlights the dimension estimation aspect, it can be used as a stratification finding tool by using different methods.

In particular, the predict_labels method will return integer labels for each point. The output should be interpretted as follows:

self.outlier_label (defaults to -1) means an outlier,
1 means a boundary point,
2 means a manifold point,
3 means a singular point.

We can see this with the following example on the Henneberg surface.

from dimmer.datasets import henneberg

henpoints = henneberg()

hen_gad = GAD(radii=(1.5, 2), max_dim=2, n_jobs=3)
hen_gad.fit(X=henpoints)

hen_gad_ids = hen_gad.predict_labels(X=henpoints, k=1)

hen_fig = construct_scatter_with_labels(X=henpoints, values=hen_gad_ids)

We update the title and names of traces to their proper names.

names = ["Boundary point", "Manifold point", "Singular point"]

for idx, trace in enumerate(hen_fig.select_traces()):
    hen_fig.update_traces(name=f"{names[idx]}", selector=dict(name=f"Label {idx+1}"))
hen_fig.update_layout(title_text="A stratification of the Henneberg surface")
hen_fig.show()

Constructing and Visualizing local neighborhoods

The GAD class tracks the neighborhoods used in the computation. This makes it easy to plot them to see whats going on. Here we run through a simple example from the rectangle example.

rng = np.random.default_rng(0)
x = rng.random(5000)
y = 0.5 * rng.random(5000)
data_rect = np.array(list(zip(x, y)))

rect_fig = construct_scatter_with_labels(X=data_rect, values=[0 for _ in data_rect])
rect_fig.show()

gad_rect = GAD(radii=(0.08, 0.12), max_dim=2, n_jobs=3)
gad_rect.fit(X=data_rect)
gad_ids = gad_rect.predict(X=data_rect)

Comparing manifold vs boundary vs corner point neighborhoods

First, lets determine a manifold index.

mfld_index = np.where(gad_ids == 2)[0][0]
print(mfld_index)

The attribute nbhd_dict is a dictionary whose keys are indices of X and whose values are the original indices of those points in the annular shell around the key. The outliers of X are removed from this dictionary, so calling their index can raise a KeyError. The updated neighborhoods are stored in nbhd_dict_update_, which can be useful for debugging.

These dictionaries can be manually queried as shown below.

Alternatively, the method plot_annular_neighborhood will plot these neighborhoods quickly. We first show what the full annular neighborhood looks like using nbhd_dict.

mfld_fig = plot_annular_neighborhood(
    gad=gad_rect, X=data_rect, idx=mfld_index, values=gad_ids
)

mfld_fig.update_scenes(xaxis_visible=False, yaxis_visible=False)
mfld_fig.update_layout(yaxis=dict(scaleanchor="x"))
mfld_fig.update_layout(
    {"plot_bgcolor": "rgba(0,0,0,0)", "paper_bgcolor": "rgba(0,0,0,0)"}
)
mfld_fig.update_layout(title_text="Annular neighborhood around a manifold point")
mfld_fig.show()

Let’s visualize a boundary point next.

bdy_index = np.where(gad_ids == 1)[0][0]
print(bdy_index)

bdy_fig = plot_annular_neighborhood(
    gad=gad_rect, X=data_rect, idx=bdy_index, values=gad_ids
)
bdy_fig.update_scenes(xaxis_visible=False, yaxis_visible=False)
bdy_fig.update_layout(yaxis=dict(scaleanchor="x"))
bdy_fig.update_layout(
    {
        "plot_bgcolor": "rgba(0,0,0,0)",
        "paper_bgcolor": "rgba(0,0,0,0)",
        "title_text": "Annular neighborhood around a boundary point.",
    }
)
bdy_fig.show()

Finally, let’s see what the local neighborhood of a corner point looks like.

corner_idx = np.where(gad_ids == 0)[0][5]
print(corner_idx)

corner_fig = plot_annular_neighborhood(
    gad=gad_rect, X=data_rect, idx=corner_idx, values=gad_ids
)
corner_fig.update_scenes(xaxis_visible=False, yaxis_visible=False)
corner_fig.update_layout(yaxis=dict(scaleanchor="x"))
corner_fig.update_layout(
    {
        "plot_bgcolor": "rgba(0,0,0,0)",
        "paper_bgcolor": "rgba(0,0,0,0)",
        "title_text": "Annular neighborhood around a corner point.",
    }
)
corner_fig.show()

It is difficult to distinguish between the boundary and the corner cases by looking at the original neighborhoods. In other words, both look contractible from a topological perspective. Since GAD is an iterative algorithm, we can see the difference if we look at their updated neighborhoods. We plot this by setting the original_nbhd flag of plot_annular_neighborhood to False:

point_idx = bdy_index  # boundary point
updated_bdy_fig = plot_annular_neighborhood(
    gad=gad_rect, X=data_rect, idx=point_idx, original_nbhd=False
)
updated_bdy_fig.update_scenes(xaxis_visible=False, yaxis_visible=False)
updated_bdy_fig.update_layout(yaxis=dict(scaleanchor="x"))
updated_bdy_fig.update_layout(
    {
        "plot_bgcolor": "rgba(0,0,0,0)",
        "paper_bgcolor": "rgba(0,0,0,0)",
        "title_text": "Updated annular neighborhood around a boundary point.",
    }
)
updated_bdy_fig.show()

point_idx = corner_idx  # corner point
updated_corner_fig = plot_annular_neighborhood(
    gad=gad_rect, X=data_rect, idx=point_idx, original_nbhd=False
)
updated_corner_fig.update_scenes(xaxis_visible=False, yaxis_visible=False)
updated_corner_fig.update_layout(yaxis=dict(scaleanchor="x"))
updated_corner_fig.update_layout(
    {
        "plot_bgcolor": "rgba(0,0,0,0)",
        "paper_bgcolor": "rgba(0,0,0,0)",
        "title_text": "Updated annular neighborhood around a corner point.",
    }
)
updated_corner_fig.show()

The updated neighborhood around a boundary part is disjoint, whereas around a corner point, it is connected. This distinction allows GAD to differentiate the two, even though their original neighborhoods are both topologically simple.