Note
Click here to download the full example code
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.
Out:
Using default location ~/mne_data for somato...
0%| | 0.00/611M [00:00<?, ?B/s]
1%|▎ | 5.71M/611M [00:00<00:10, 57.1MB/s]
2%|▋ | 11.4M/611M [00:00<00:12, 46.9MB/s]
3%|█ | 16.2M/611M [00:00<00:13, 44.8MB/s]
3%|█▎ | 20.8M/611M [00:00<00:15, 37.3MB/s]
4%|█▌ | 24.6M/611M [00:00<00:16, 36.5MB/s]
5%|█▊ | 29.6M/611M [00:00<00:14, 40.5MB/s]
6%|██▏ | 35.7M/611M [00:00<00:12, 46.5MB/s]
7%|██▌ | 41.2M/611M [00:00<00:11, 48.9MB/s]
8%|██▉ | 47.9M/611M [00:01<00:10, 54.1MB/s]
9%|███▍ | 54.5M/611M [00:01<00:09, 57.8MB/s]
10%|███▊ | 61.2M/611M [00:01<00:09, 60.7MB/s]
11%|████▏ | 68.0M/611M [00:01<00:08, 62.6MB/s]
12%|████▋ | 74.8M/611M [00:01<00:08, 64.4MB/s]
13%|█████ | 81.8M/611M [00:01<00:07, 66.1MB/s]
15%|█████▌ | 88.9M/611M [00:01<00:07, 67.3MB/s]
16%|█████▉ | 95.9M/611M [00:01<00:07, 68.1MB/s]
17%|██████▌ | 103M/611M [00:01<00:07, 68.7MB/s]
18%|███████ | 110M/611M [00:01<00:07, 69.1MB/s]
19%|███████▍ | 117M/611M [00:02<00:07, 69.5MB/s]
20%|███████▉ | 124M/611M [00:02<00:06, 69.7MB/s]
21%|████████▎ | 131M/611M [00:02<00:06, 69.8MB/s]
23%|████████▊ | 138M/611M [00:02<00:06, 70.0MB/s]
24%|█████████▎ | 145M/611M [00:02<00:06, 70.1MB/s]
25%|█████████▋ | 152M/611M [00:02<00:06, 70.4MB/s]
26%|██████████▏ | 159M/611M [00:02<00:06, 70.5MB/s]
27%|██████████▌ | 166M/611M [00:02<00:06, 70.5MB/s]
28%|███████████ | 173M/611M [00:02<00:06, 70.6MB/s]
30%|███████████▌ | 180M/611M [00:02<00:06, 70.6MB/s]
31%|███████████▉ | 188M/611M [00:03<00:05, 70.7MB/s]
32%|████████████▍ | 195M/611M [00:03<00:05, 70.8MB/s]
33%|████████████▉ | 202M/611M [00:03<00:05, 70.9MB/s]
34%|█████████████▎ | 209M/611M [00:03<00:05, 70.9MB/s]
35%|█████████████▊ | 216M/611M [00:03<00:05, 70.9MB/s]
37%|██████████████▏ | 223M/611M [00:03<00:05, 70.9MB/s]
38%|██████████████▋ | 230M/611M [00:03<00:05, 70.8MB/s]
39%|███████████████▏ | 237M/611M [00:03<00:05, 70.8MB/s]
40%|███████████████▌ | 244M/611M [00:03<00:05, 70.8MB/s]
41%|████████████████ | 251M/611M [00:03<00:05, 70.9MB/s]
42%|████████████████▌ | 258M/611M [00:04<00:04, 70.9MB/s]
44%|████████████████▉ | 266M/611M [00:04<00:04, 71.0MB/s]
45%|█████████████████▍ | 273M/611M [00:04<00:04, 70.9MB/s]
46%|█████████████████▊ | 280M/611M [00:04<00:04, 70.8MB/s]
47%|██████████████████▎ | 287M/611M [00:04<00:04, 70.9MB/s]
48%|██████████████████▊ | 294M/611M [00:04<00:04, 70.9MB/s]
49%|███████████████████▏ | 301M/611M [00:04<00:04, 70.9MB/s]
50%|███████████████████▋ | 308M/611M [00:04<00:04, 71.0MB/s]
52%|████████████████████▏ | 315M/611M [00:04<00:04, 71.0MB/s]
53%|████████████████████▌ | 322M/611M [00:04<00:04, 71.0MB/s]
54%|█████████████████████ | 330M/611M [00:05<00:03, 71.0MB/s]
55%|█████████████████████▌ | 337M/611M [00:05<00:03, 71.1MB/s]
56%|█████████████████████▉ | 344M/611M [00:05<00:03, 71.0MB/s]
57%|██████████████████████▍ | 351M/611M [00:05<00:03, 71.1MB/s]
59%|██████████████████████▊ | 358M/611M [00:05<00:03, 71.0MB/s]
60%|███████████████████████▎ | 365M/611M [00:05<00:03, 71.0MB/s]
61%|███████████████████████▊ | 372M/611M [00:05<00:03, 71.1MB/s]
62%|████████████████████████▏ | 379M/611M [00:05<00:03, 58.4MB/s]
63%|████████████████████████▋ | 386M/611M [00:05<00:03, 58.0MB/s]
64%|█████████████████████████ | 393M/611M [00:06<00:03, 61.3MB/s]
65%|█████████████████████████▌ | 400M/611M [00:06<00:03, 63.8MB/s]
67%|█████████████████████████▉ | 407M/611M [00:06<00:03, 65.7MB/s]
68%|██████████████████████████▍ | 414M/611M [00:06<00:02, 67.2MB/s]
69%|██████████████████████████▉ | 421M/611M [00:06<00:02, 68.2MB/s]
70%|███████████████████████████▎ | 428M/611M [00:06<00:02, 68.8MB/s]
71%|███████████████████████████▊ | 435M/611M [00:06<00:02, 69.2MB/s]
72%|████████████████████████████▏ | 442M/611M [00:06<00:02, 69.7MB/s]
74%|████████████████████████████▋ | 449M/611M [00:06<00:02, 69.9MB/s]
75%|█████████████████████████████▏ | 456M/611M [00:06<00:02, 70.2MB/s]
76%|█████████████████████████████▌ | 463M/611M [00:07<00:02, 70.3MB/s]
77%|██████████████████████████████ | 470M/611M [00:07<00:01, 70.4MB/s]
78%|██████████████████████████████▍ | 477M/611M [00:07<00:01, 70.7MB/s]
79%|██████████████████████████████▉ | 484M/611M [00:07<00:01, 70.7MB/s]
81%|███████████████████████████████▍ | 492M/611M [00:07<00:01, 70.7MB/s]
82%|███████████████████████████████▊ | 499M/611M [00:07<00:01, 70.7MB/s]
83%|████████████████████████████████▎ | 506M/611M [00:07<00:01, 70.9MB/s]
84%|████████████████████████████████▊ | 513M/611M [00:07<00:01, 71.0MB/s]
85%|█████████████████████████████████▏ | 520M/611M [00:07<00:01, 70.9MB/s]
86%|█████████████████████████████████▋ | 527M/611M [00:07<00:01, 70.9MB/s]
87%|██████████████████████████████████ | 534M/611M [00:08<00:01, 70.8MB/s]
89%|██████████████████████████████████▌ | 541M/611M [00:08<00:00, 70.8MB/s]
90%|███████████████████████████████████ | 548M/611M [00:08<00:00, 70.8MB/s]
91%|███████████████████████████████████▍ | 555M/611M [00:08<00:00, 70.8MB/s]
92%|███████████████████████████████████▉ | 562M/611M [00:08<00:00, 70.4MB/s]
93%|████████████████████████████████████▍ | 570M/611M [00:08<00:00, 70.5MB/s]
94%|████████████████████████████████████▊ | 577M/611M [00:08<00:00, 70.6MB/s]
96%|█████████████████████████████████████▎ | 584M/611M [00:08<00:00, 70.9MB/s]
97%|█████████████████████████████████████▋ | 591M/611M [00:08<00:00, 70.7MB/s]
98%|██████████████████████████████████████▏| 598M/611M [00:08<00:00, 70.7MB/s]
99%|██████████████████████████████████████▋| 605M/611M [00:09<00:00, 70.8MB/s]
0%| | 0.00/611M [00:00<?, ?B/s]
100%|███████████████████████████████████████| 611M/611M [00:00<00:00, 1.22TB/s]
Attempting to create new mne-python configuration file:
/github/home/.mne/mne-python.json
Opening raw data file /github/home/mne_data/MNE-somato-data/sub-01/meg/sub-01_task-somato_meg.fif...
Range : 237600 ... 506999 = 791.189 ... 1688.266 secs
Ready.
Reading 0 ... 269399 = 0.000 ... 897.077 secs...
Setting up band-stop filter
FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandstop filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower transition bandwidth: 0.50 Hz
- Upper transition bandwidth: 0.50 Hz
- Filter length: 1983 samples (6.603 sec)
Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 2 Hz
FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 2.00
- Lower transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 1.00 Hz)
- Filter length: 497 samples (1.655 sec)
111 events found
Event IDs: [1]
Not setting metadata
111 matching events found
Setting baseline interval to [-3.9992341833870637, 0.0] sec
Applying baseline correction (mode: mean)
0 projection items activated
Using data from preloaded Raw for 111 events and 1202 original time points ...
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
5 bad epochs dropped
Computing rank from data with rank=None
Using tolerance 7.5e-09 (2.2e-16 eps * 204 dim * 1.7e+05 max singular value)
Estimated rank (grad): 204
GRAD: rank 204 computed from 204 data channels with 0 projectors
/github/workspace/alphacsc/datasets/mne_data.py:94: RuntimeWarning: Something went wrong in the data-driven estimation of the data rank as it exceeds the theoretical rank from the info (204 > 64). Consider setting rank to "auto" or setting it explicitly as an integer.
cov = mne.compute_covariance(epochs_cov)
Reducing data rank from 204 -> 204
Estimating covariance using EMPIRICAL
Done.
Number of samples used : 127412
[done]
Not setting metadata
111 matching events found
Setting baseline interval to [-2.001282051803185, 0.0] sec
Applying baseline correction (mode: mean)
0 projection items activated
Using data from preloaded Raw for 111 events and 1803 original time points ...
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
8 bad epochs dropped
Fit the model and learn rank1 atoms
cdl.fit(X)
Out:
./github/workspace/alphacsc/utils/optim.py:136: DeprecationWarning: Please use `scalar_search_armijo` from the `scipy.optimize` namespace, the `scipy.optimize.linesearch` namespace is deprecated.
step_size, obj_uv = optimize.linesearch.scalar_search_armijo(
.............
[BatchCDL] Converged after 14 iteration, (dz, du) = 8.570e-05, 8.726e-05
[BatchCDL] Fit in 873.8s
<alphacsc.convolutional_dictionary_learning.BatchCDL object at 0x7f587ff16ca0>
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 * 5, 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()
Total running time of the script: ( 15 minutes 10.788 seconds)