import logging
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from sklearn.decomposition import KernelPCA
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
logging.basicConfig(level=logging.INFO, format="%(levelname)s - %(message)s")
logger = logging.getLogger(__name__)
def unit_vectorize(a):
"""
Convert an array with (..., N) angles (in radians) into an array with
(..., N, 2) sine and cosine values for the N angles.
Parameters
----------
a : np.ndarray
Array of angles, in radians.
Returns
-------
v : np.ndarray
Array of shape ``a.shape + (2,)`` with cosine and sine values.
"""
v = np.concatenate([np.cos(a)[..., None], np.sin(a)[..., None]], axis=-1)
return v
def unit_vector_distance(a0, a1, sqrt=True):
"""
Compute the sum of distances between two (..., N) arrays storing the
values of N angles (in radians), using their unit vector representation.
Parameters
----------
a0 : np.ndarray
First array of angles.
a1 : np.ndarray
Second array of angles.
sqrt : bool, optional
Take the square root of the squared distance. The default is True.
Returns
-------
dist : np.ndarray
Distance(s) between a0 and a1.
"""
v0 = unit_vectorize(a0)
v1 = unit_vectorize(a1)
if sqrt:
dist = np.sqrt(np.square(v0 - v1).sum(axis=-1))
else:
dist = np.square(v0 - v1).sum(axis=-1)
dist = dist.sum(axis=-1)
return dist
def unit_vector_kernel(a1, a2, gamma):
"""
Compute a similarity kernel between two arrays of angles (in radians),
based on their unit vector distance.
Parameters
----------
a1 : np.ndarray
First array of angles.
a2 : np.ndarray
Second array of angles.
gamma : float
Kernel coefficient.
Returns
-------
sim : np.ndarray
Similarity value(s).
"""
dist = unit_vector_distance(a1, a2, sqrt=False)
sim = np.exp(-gamma * dist)
return sim
[docs]
def compute_kpca(
df,
feature_cols=None,
circular=False,
n_components=10,
kernel="poly",
gamma=None,
plot=True,
s=5,
color="C0",
context="notebook",
figsize=(10, 5),
edgecolors="black",
ax=None,
):
"""
Reduce the dimensionality of the data using Kernel PCA (KPCA).
Parameters
----------
df : pandas.DataFrame
Data to reduce. If `feature_cols` is given, only those columns are
used to fit KPCA (the returned dataframe keeps ALL of `df`'s
original columns, plus the KPCA components) -- use this whenever
`df` also carries metadata columns (sim_name, replica, frame,
...) that must not be fed to KPCA itself. If `feature_cols` is
None (default, backward-compatible), `df` is used as-is -- every
column must be numeric.
feature_cols : list of str, optional
Column names in `df` to actually reduce. Default None (use all of
`df`).
circular : bool, optional
If True, use a custom similarity kernel suited for angular
features (based on the unit vector representation of the angles,
in radians). This overrides the `kernel` argument. The default is
False.
n_components : int, optional
Number of components to keep. The default is 10.
kernel : str, optional
Kernel used for KPCA, as in the scikit-learn implementation.
Ignored if `circular` is True. The default is "poly".
gamma : float, optional
Kernel coefficient. If None and `circular` is True, it is set to
1 / number of features. The default is None.
plot : bool, optional
Plot the two first components. The default is True.
s : int, optional
Marker size for the scatter plot. The default is 5.
color : str, optional
Marker color for the scatter plot. The default is "C0".
context : str, optional
Seaborn plotting context. The default is "notebook".
figsize : tuple, optional
Figure size. The default is (10, 5).
ax : matplotlib.axes._subplots.AxesSubplot, optional
Axes on which to plot. The default is None.
Returns
-------
new_df : pandas.DataFrame
`df` (all original columns preserved) with the KPCA components
added (columns "KPC1", "KPC2", ...).
fig : matplotlib.figure.Figure or None
Figure containing the plot (if plot=True) -- always retrievable
as an object for further tweaking, e.g. `fig.savefig(...)`.
kpca : sklearn.decomposition.KernelPCA
The fitted KernelPCA object.
Note
----
KPCA computes an N x N kernel matrix (N = number of rows) -- O(N^2)
memory and O(N^3) time for the eigendecomposition. This becomes slow
or memory-heavy well before N reaches tens of thousands of rows;
subsample `df`/`feature_cols` first (e.g. `df.sample(n=2000)`) on
large trajectories.
"""
X = df[feature_cols] if feature_cols is not None else df
if circular:
gamma = 1 / X.shape[1] if gamma is None else gamma
logger.info("Running KPCA with a circular kernel, gamma=%.4f", gamma)
kpca_kernel = lambda a1, a2: unit_vector_kernel(a1, a2, gamma=gamma)
else:
logger.info("Running KPCA with kernel=%s", kernel)
kpca_kernel = kernel
kpca = KernelPCA(
n_components=n_components,
kernel=kpca_kernel,
gamma=gamma,
)
reduced = kpca.fit_transform(X)
columns = [f"KPC{i + 1}" for i in range(reduced.shape[1])]
logger.info("KPCA completed.")
# Positional assignment: reduced preserves X's (and therefore df's)
# row order, so a plain array lines up correctly with new_df's rows.
new_df = df.copy()
new_df[columns] = reduced
fig = None
if plot:
logger.info("Plotting KPCA results...")
sns.set_context(context)
if ax is None:
fig, ax = plt.subplots(figsize=figsize)
else:
fig = ax.figure
ax.scatter(
new_df[columns[0]],
new_df[columns[1]],
s=s,
color=color,
edgecolors=edgecolors,
)
ax.set_title("KPCA")
ax.set_xlabel(columns[0])
ax.set_ylabel(columns[1])
fig.tight_layout()
return new_df, fig, kpca