import numpy as np
import plotly.express as px
import plotly.io as pio
from gtda.homology import VietorisRipsPersistence
from sklearn.cluster import DBSCAN, AgglomerativeClustering, KMeans, SpectralClustering
from sklearn.manifold import Isomap
from dimmer.datasets import cyclooctane, henneberg
from dimmer.gad import GAD
from dimmer.visualization import construct_scatter_with_labels
pio.renderers.default = "plotly_mimetype+notebook"Building Stratifications and Chart Decompositions
In this notebook, we use GAD to construct a stratification of a space and a set of charts. Recall that GAD provides a way of assigning a local dimension estimate to each data point. Looking at how the points in the dataset are grouped as a function of local dimension is essentially the traditional notion of “stratification”. If we push slightly further and understand how each of the equidimensional pieces are incident amongst eachother and their incidence with other dimensional pieces then we arrive at the usual notion of stratification.
We can refine this notion even further to demand that each equidimensional piece of the stratification (a “strata”) be as topologically simple as possible. Further decomposing a stratification this way leads to a “chartification” of a space– a way to understand how the space is built out of charts. This perspective is useful for performing ML on complicated datasets, since ML algorithms (from SVMs to VAEs) tend to perform better on topologically simpler pieces.
Idea: Apply GAD to determine the top dimensional points. Working one dimension at a time decreasing from the maximum dimension, we run various clustering algorithms on the equidimensional pieces to construct each chart within the stratum.
There are many clustering algorithms with many hyperparameters that can be chosen. This notebook explores these choices in a handful of examples.
Decomposing Cyclooctane
We start with the conformation space of cyclooctane. We first attempt to find the top-dimensional strata by finding the high-dimensional points. We use GAD to classify the points and then isolate those which are two dimensional.
cyc_data = cyclooctane()
cyc_gad = GAD(radii=(0.45, 0.7), max_dim=2, n_jobs=3)
cyc_gad_ids = cyc_gad.fit_predict(X=cyc_data)
np.count_nonzero(cyc_gad_ids == 2)4721GAD classifies 4721 points as being 2-dimensional. Since this is the top dimension, we start here and attempt to construct chart(s) for only these points, leaving the remaining points for now.
iso = Isomap(n_components=3)
cyc_iso = iso.fit_transform(X=cyc_data)Spectral Clustering
The number of clusters for spectral clustering is a free parameter of our choice. Regardless, we will see that this method does not cluster well. We instantiate a spectral clusterer and then use it to predict labels on only the 2-dimensional points.
sc = SpectralClustering(n_clusters=8)sc_labels = np.full(cyc_data.shape[0], -1.0)
sc_labels[cyc_gad_ids == 2] = sc.fit_predict(X=cyc_data[cyc_gad_ids == 2])sc_fig = construct_scatter_with_labels(X=cyc_iso, values=sc_labels)
sc_fig.show()The points labelled “-1” are those that are not 2-dimensional so we ignore them (clicking their legend entry once should remove them from the visualization). The remaining points are not topologically simple in any meaningful way. Looking at any individual label (by double-clicking on their legend entry) shows three distinct components.
If you are running this notebook, you’re encouraged to go back and try modifying the n_clusters parameter of sc to see if you can get “charts” which are connected and without loops.
DBSCAN
DBSCAN is a density based clustering scheme that depends heavily on an eps parameter that controls what we mean by ‘dense’. This clustering scheme will be very useful for us!
To determine the density parameter eps, we can investigate the persistence diagram and look for gaps in \(PH_0\).
vr_2 = VietorisRipsPersistence(homology_dimensions=(0,)).fit_transform_plot(
X=[cyc_data[cyc_gad_ids == 2]]
)There are appears to be a gap in \(PH_0\) from about 0.21 to 0.45, so any eps cut-off parameter within that range should work.
db = DBSCAN(eps=0.4)
db_labels = np.full(cyc_data.shape[0], -1)
db_labels[cyc_gad_ids == 2] = db.fit_predict(X=cyc_data[cyc_gad_ids == 2])db_fig = construct_scatter_with_labels(X=cyc_iso, values=db_labels)
db_fig.show()Removing the points labelled -1 again, we see four disjoint components. This looks much better than the previous picture! However the middle band is certainly not contractible and so that piece needs to be decomposed further. There are a variety of methods that can be used to decompose this. Easily enough, one can intersect with half-spaces in this example.
Pushing further, lets stratify the codimension one and two pieces as well. As before, let’s compute the persistence diagram of those points labelled as 1-dimensional.
vr_1 = VietorisRipsPersistence(homology_dimensions=(0,)).fit_transform_plot(
X=[cyc_data[cyc_gad_ids == 1]]
)The same eps parameter of 0.4 seems like it will work well for the codimension one piece as well.
Ultimately, we want to give each point of the original dataset a labelling corresponding to its chart not just its dimension. Therefore, we offset the dimension one labels so they don’t conflict with the dimension 2 labels.
db_labels[cyc_gad_ids == 1] = (
db.fit_predict(X=cyc_data[cyc_gad_ids == 1]) + np.max((db_labels)) + 2
)db_fig = construct_scatter_with_labels(X=cyc_iso, values=db_labels)
db_fig.show()Single linkage clustering
With the same philosophy, we apply single linkage with a similar distance threshold. \(PH_0\) is nothing more than an algebrified version of single linkage, so inferring the parameter from the persistence diagram is even more direct here.
ag = AgglomerativeClustering(linkage="single", distance_threshold=0.4, n_clusters=None)
ag_labels = np.full(cyc_data.shape[0], -1)
ag_labels[cyc_gad_ids == 2] = ag.fit_predict(X=cyc_data[cyc_gad_ids == 2])ag_fig = construct_scatter_with_labels(
X=cyc_iso, values=ag_labels, color_sequence=px.colors.qualitative.Safe
)
ag_fig.show()K-means
km = KMeans(n_clusters=4)
km_labels = np.full(cyc_data.shape[0], -1)
km_labels[cyc_gad_ids == 2] = km.fit_predict(X=cyc_data[cyc_gad_ids == 2])/home/michael-catanzaro/Projects/dimmer/.venv/lib/python3.12/site-packages/sklearn/cluster/_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
super()._check_params_vs_input(X, default_n_init=10)km_fig = construct_scatter_with_labels(
X=cyc_iso, values=km_labels, color_sequence=px.colors.qualitative.Safe
)
km_fig.show()Just as with spectral clustering, k-means clustering does not seem to be appropriate for constructing charts.
We have discovered that both DBSCAN and single-linkage clustering work well for turning stratifications into chart decompositions. Let’s try these same ideas on another example.
Henneberg surface
hen_pts = henneberg()hen_gad = GAD(radii=(1.15, 1.65), max_dim=2, n_jobs=3)
hen_gad_ids = hen_gad.fit_predict(X=hen_pts)Spectral Clustering
We’ll try spectral clustering one more time. Once again, the number of clusters for spectral clustering is a free parameter of our choice.
sc = SpectralClustering(n_clusters=8)sc_labels = np.full(hen_pts.shape[0], -1.0)
sc_labels[hen_gad_ids == 2] = sc.fit_predict(X=hen_pts[hen_gad_ids == 2])sc_fig = construct_scatter_with_labels(
X=hen_pts, values=sc_labels, label="Chart {value}"
)
sc_fig.show()Now Spectral clustering performs extremely well!
DBSCAN
DBSCAN is a density based clustering scheme that depends heavily on an eps parameter that controls what we mean by ‘dense’. To determine this parameter, we can investigate the persistence diagram and look for gaps in \(PH_0\).
vr_2 = VietorisRipsPersistence(homology_dimensions=(0,)).fit_transform_plot(
X=[hen_pts[hen_gad_ids == 2]]
)There are appears to be a gap in \(PH_0\) from 0.85 to 0.92. Let’s try an eps parameter within that range.
db = DBSCAN(eps=0.9)
db_labels = np.full(hen_pts.shape[0], -1)
db_labels[hen_gad_ids == 2] = db.fit_predict(X=hen_pts[hen_gad_ids == 2])db_fig = construct_scatter_with_labels(X=hen_pts, values=db_labels)
db_fig.show()Only one two-dimensional component of the stratification is isolated! We need to modify the eps parameter.
db = DBSCAN(eps=0.7)
db_labels = np.full(hen_pts.shape[0], -1)
db_labels[hen_gad_ids == 2] = db.fit_predict(X=hen_pts[hen_gad_ids == 2])np.unique(db_labels)array([-1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20])db_fig = construct_scatter_with_labels(X=hen_pts, values=db_labels)
db_fig.show()This works pretty well!
Single linkage clustering
With the same philosophy, we apply single linkage with a similar distance threshold. \(PH_0\) is nothing more than an algebrified version of single linkage, so inferring the parameter from the persistence diagram is even more direct here.
ag = AgglomerativeClustering(
linkage="single", distance_threshold=0.5, n_clusters=None
) # (linkage='ward',n_clusters=5)#,
ag_labels = np.full(hen_pts.shape[0], -1)
ag_labels[hen_gad_ids == 2] = ag.fit_predict(X=hen_pts[hen_gad_ids == 2])ag_fig = construct_scatter_with_labels(
X=hen_pts, values=ag_labels, color_sequence=px.colors.qualitative.Safe
)
ag_fig.update_layout(legend={})
ag_fig.show()