{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"id": "5fb0b379",
"metadata": {},
"source": [
"# 5. Scattering coefficients reconstruction\n",
"\n",
"After retrieving the most relevant features with an ICA from the scattering coefficients, we can now check what these features are actually telling us."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "07a95975",
"metadata": {},
"outputs": [],
"source": [
"import pickle\n",
"\n",
"import numpy as np\n",
"import matplotlib.pylab as plt"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "49a58946",
"metadata": {},
"outputs": [],
"source": [
"%config InlineBackend.figure_format = \"svg\""
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "1f86a5c2",
"metadata": {},
"source": [
"## Load the data\n",
"\n",
"For this notebook we need to load the scattering network model, the dimensionality reduction model and the independent components."
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "f07c3d55",
"metadata": {},
"outputs": [],
"source": [
"# Load the scattering network\n",
"network = pickle.load(open(\"../example/scattering_network.pickle\", \"rb\"))\n",
"\n",
"# Load the dimensionality reduction model\n",
"dimension_model = pickle.load(open(\"../example/dimension_model.pickle\", \"rb\"))\n",
"\n",
"# Load data from file\n",
"with np.load(\n",
" \"../example/independent_components.npz\", allow_pickle=True\n",
") as data:\n",
" features = data[\"features\"]\n",
" times = data[\"times\"]"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "07564103",
"metadata": {},
"source": [
"## Reconstruct scattering coefficients from one feature\n",
"\n",
"FastICA finds a set of independent components $\\mathbf{s}$ by multiplying an unmixing matrix $\\mathbf{W}$ to the scattering coefficients stored in the matrix $\\mathbf{x}$: \n",
"\n",
"$$\\mathbf{s} = \\mathbf{W} \\mathbf{x}$$\n",
"\n",
"The FastICA model stores the unmixing matrix $\\mathbf{W}$ in the attribute `.components_`. We can extract the weights of the unmixing matrix and visualize it to understand better what each independent component means, since it connects the scattering coefficients to the independent components.\n",
"\n",
"To give an example on how to use the information of the unmixing matrix, we choose one single independent component and extract the unmixing weights for that component. The feature we have chosen represents the evolution of the tremor signal over the day. Note that the sign of a feature is arbitrary—a negative or positive change in the feature amplitude means only a relative change. Without the unmixing weights we cannot tell if the feature amplitude change is related to an amplitude increase or decrease in the scattering coefficient matrix."
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "9e0f691c",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Pick a feature\n",
"feature_id = 2\n",
"feature = features[:, feature_id]\n",
"\n",
"# Figure creation\n",
"fig = plt.figure(figsize=(10, 3))\n",
"ax = plt.axes()\n",
"\n",
"# Plot the weights\n",
"ax.plot(times, feature)\n",
"ax.set_ylabel(f\"Amplitude of feature {feature_id}\")\n",
"\n",
"# Show\n",
"plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "2a2d6062",
"metadata": {},
"source": [
"## Visualize feature direction in scattering space\n",
"\n",
"The top row shows the unmixing weights for the first-order scattering coefficients, which all show negative values with a minimum at a first-order frequency $f_1$ of 2 Hz. The bottom row shows the unmixing weights for the second-order coefficients, which show both positive and negative values."
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "1af97c15",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Etract weights from the dimensionality reduction model\n",
"weights = dimension_model.components_[feature_id]\n",
"vmax = np.abs(weights).max()\n",
"\n",
"# Scattering coefficients shape and frequencies\n",
"n_cha = 3\n",
"n_order_1 = network.banks[0].octaves * network.banks[0].resolution\n",
"n_order_2 = network.banks[1].octaves * network.banks[1].resolution\n",
"f_1 = network.banks[0].centers\n",
"f_2 = network.banks[1].centers\n",
"\n",
"# Extract and reshape weights\n",
"order_1 = weights[: n_cha * n_order_1].reshape(n_cha, n_order_1)\n",
"order_2 = weights[n_cha * n_order_1 :].reshape(n_cha, n_order_1, n_order_2)\n",
"\n",
"# Show weights\n",
"fig, ax = plt.subplots(nrows=2, ncols=3, sharex=True, sharey=\"row\")\n",
"image_kw = dict(vmin=-vmax, vmax=vmax, rasterized=True, cmap=\"PuOr\")\n",
"for id, channel in enumerate(\"ENZ\"):\n",
"\n",
" # Show\n",
" ax[0, id].plot(f_1, order_1[id], label=channel)\n",
" mappable = ax[1, id].pcolormesh(f_1, f_2, order_2[id].T, **image_kw)\n",
"\n",
" # Labels\n",
" ax[0, id].set_title(channel)\n",
" ax[1, id].set_xlabel(\"$f_1$ (Hz)\")\n",
"\n",
"# Labels\n",
"ax[0, 0].set_ylabel(\"Unmixing weights\")\n",
"ax[1, 0].set_ylabel(\"$f_2$ (Hz)\")\n",
"ax[1, 0].set_xscale(\"log\")\n",
"ax[1, 0].set_yscale(\"log\")\n",
"\n",
"# Colorbar\n",
"colorbar = fig.colorbar(mappable, orientation=\"horizontal\", ax=ax, shrink=0.3)\n",
"colorbar.set_label(\"Unmixing weights\")\n",
"\n",
"# Show\n",
"plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "aa5a52f2",
"metadata": {},
"source": [
"## Interpret partial reconstruction\n",
"\n",
"We here reconstruct the scattering coefficients from a selected independent component. We only show the first order coefficients, since the second order coefficients are hard to map as a function of time. \n",
"\n",
"The reconstruction shows clearly the evolving tremor signal, starting between 9:00 and 12:00. The tremor signal is characterized by an increased amplitude peaking around 2 Hz. The negative values in the unmixing matrix were applied to the scattering coefficients matrix to retrieve feature 2, tracking the tremor onset with an amplitude decrease in the feature. If feature 2 would track the tremor signal with an amplitude increase, the weights in the unmixing matrix would have had the opposite sign."
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "61278419",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Filter out latent space\n",
"features_filtered = np.zeros(features.shape)\n",
"features_filtered[:, feature_id] = feature\n",
"\n",
"# Extract all scattering coefficients\n",
"reconstructed = dimension_model.inverse_transform(features_filtered)\n",
"reconstructed_order_1 = reconstructed[:, : n_cha * n_order_1].reshape(\n",
" -1, n_cha, n_order_1\n",
")\n",
"vmin = reconstructed_order_1.min()\n",
"vmax = reconstructed_order_1.max()\n",
"\n",
"# Plot\n",
"fig, ax = plt.subplots(\n",
" nrows=3, sharex=True, sharey=\"row\", constrained_layout=True\n",
")\n",
"\n",
"# Plot\n",
"for id, channel in enumerate(\"ENZ\"):\n",
" data = reconstructed_order_1[:, id, :].squeeze().T\n",
" mappable = ax[id].pcolormesh(\n",
" times, f_1, data, rasterized=True, vmin=vmin, vmax=vmax\n",
" )\n",
" ax[id].set_ylabel(\"$f_1$ (Hz)\")\n",
" ax[id].set_yscale(\"log\")\n",
"\n",
"# Colorbar\n",
"colorbar = fig.colorbar(mappable, orientation=\"vertical\", ax=ax, shrink=0.5)\n",
"colorbar.set_label(\"Scattering coefficients\")\n",
"fig.autofmt_xdate()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": ".venv",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.7"
}
},
"nbformat": 4,
"nbformat_minor": 5
}