Extracting \(\mu\)-wave from the somato-sensory dataset

This example illustrates how to learn rank-1 atoms 1 on the multivariate somato-sensorymotor dataset from mne. The displayed results highlight the presence of \(\mu\)-waves located in the SI cortex.

1

Dupré La Tour, T., Moreau, T., Jas, M., & Gramfort, A. (2018). Multivariate Convolutional Sparse Coding for Electromagnetic Brain Signals. Advances in Neural Information Processing Systems (NIPS).

# Authors: Thomas Moreau <thomas.moreau@inria.fr>
#          Mainak Jas <mainak.jas@telecom-paristech.fr>
#          Tom Dupre La Tour <tom.duprelatour@telecom-paristech.fr>
#          Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>
#
# License: BSD (3-clause)

Let us first define the parameters of our model.

sfreq = 150.

# Define the shape of the dictionary
n_atoms = 25
n_times_atom = int(round(sfreq * 1.0))  # 1000. ms

Next, we define the parameters for multivariate CSC

from alphacsc import BatchCDL
cdl = BatchCDL(
    # Shape of the dictionary
    n_atoms=n_atoms,
    n_times_atom=n_times_atom,
    # Request a rank1 dictionary with unit norm temporal and spatial maps
    rank1=True, uv_constraint='separate',
    # Initialize the dictionary with random chunk from the data
    D_init='chunk',
    # rescale the regularization parameter to be 20% of lambda_max
    lmbd_max="scaled", reg=.2,
    # Number of iteration for the alternate minimization and cvg threshold
    n_iter=100, eps=1e-4,
    # solver for the z-step
    solver_z="lgcd", solver_z_kwargs={'tol': 1e-2, 'max_iter': 1000},
    # solver for the d-step
    solver_d='alternate_adaptive', solver_d_kwargs={'max_iter': 300},
    # Technical parameters
    verbose=1, random_state=0, n_jobs=6)

Here, we load the data from the somato-sensory dataset and preprocess them in epochs. The epochs are selected around the stim, starting 2 seconds before and finishing 4 seconds after.

from alphacsc.datasets.somato import load_data
t_lim = (-2, 4)
X, info = load_data(epoch=t_lim, sfreq=sfreq)

Traceback (most recent call last):
  File "/home/sed-sac/hgozukan/dev/alphacsc/examples/multicsc/plot_somato_mu_waves.py", line 63, in <module>
    X, info = load_data(epoch=t_lim, sfreq=sfreq)
  File "/scratch/hgozukan/miniconda3/lib/python3.8/site-packages/joblib/memory.py", line 594, in __call__
    return self._cached_call(args, kwargs)[0]
  File "/scratch/hgozukan/miniconda3/lib/python3.8/site-packages/joblib/memory.py", line 537, in _cached_call
    out, metadata = self.call(*args, **kwargs)
  File "/scratch/hgozukan/miniconda3/lib/python3.8/site-packages/joblib/memory.py", line 779, in call
    output = self.func(*args, **kwargs)
  File "/home/sed-sac/hgozukan/dev/alphacsc/alphacsc/datasets/somato.py", line 50, in load_data
    data_path = mne.datasets.somato.data_path()
  File "</scratch/hgozukan/miniconda3/lib/python3.8/site-packages/mne/externals/decorator.py:decorator-gen-354>", line 2, in data_path
  File "/scratch/hgozukan/miniconda3/lib/python3.8/site-packages/mne/utils/_logging.py", line 90, in wrapper
    return function(*args, **kwargs)
  File "/scratch/hgozukan/miniconda3/lib/python3.8/site-packages/mne/datasets/somato/somato.py", line 19, in data_path
    return _data_path(path=path, force_update=force_update,
  File "/scratch/hgozukan/miniconda3/lib/python3.8/site-packages/mne/datasets/utils.py", line 390, in _data_path
    remove_archive, full = _download(path, u, an, h)
  File "/scratch/hgozukan/miniconda3/lib/python3.8/site-packages/mne/datasets/utils.py", line 447, in _download
    _fetch_file(url, full_name, print_destination=False,
  File "</scratch/hgozukan/miniconda3/lib/python3.8/site-packages/mne/externals/decorator.py:decorator-gen-3>", line 2, in _fetch_file
  File "/scratch/hgozukan/miniconda3/lib/python3.8/site-packages/mne/utils/_logging.py", line 90, in wrapper
    return function(*args, **kwargs)
  File "/scratch/hgozukan/miniconda3/lib/python3.8/site-packages/mne/utils/fetching.py", line 162, in _fetch_file
    _get_http(url, temp_file_name, initial_size, file_size, timeout,
  File "/scratch/hgozukan/miniconda3/lib/python3.8/site-packages/mne/utils/fetching.py", line 60, in _get_http
    with open(temp_file_name, mode) as local_file:
FileNotFoundError: [Errno 2] No such file or directory: '/home/sed-sac/hgozukan/mne_data/MNE-somato-data.tar.gz.part'

Fit the model and learn rank1 atoms

cdl.fit(X)

Display the 4-th atom, which displays a \(\mu\)-waveform in its temporal pattern.

import mne
import numpy as np
import matplotlib.pyplot as plt

i_atom = 4
n_plots = 3
figsize = (n_plots * 3.5, 5)
fig, axes = plt.subplots(1, n_plots, figsize=figsize, squeeze=False)

# Plot the spatial map of the learn atom using mne topomap
ax = axes[0, 0]
u_hat = cdl.u_hat_[i_atom]
mne.viz.plot_topomap(u_hat, info, axes=ax, show=False)
ax.set(title='Learned spatial pattern')

# Plot the temporal pattern of the learn atom
ax = axes[0, 1]
v_hat = cdl.v_hat_[i_atom]
t = np.arange(v_hat.size) / sfreq
ax.plot(t, v_hat)
ax.set(xlabel='Time (sec)', title='Learned temporal waveform')
ax.grid(True)

# Plot the psd of the time atom
ax = axes[0, 2]
psd = np.abs(np.fft.rfft(v_hat)) ** 2
frequencies = np.linspace(0, sfreq / 2.0, len(psd))
ax.semilogy(frequencies, psd)
ax.set(xlabel='Frequencies (Hz)', title='Power Spectral Density')
ax.grid(True)
ax.set_xlim(0, 30)

plt.tight_layout()
plt.show()