Note
Click here to download the full example code
Extracting artifact and evoked response atoms from the sample dataset#
This example illustrates how to learn rank-1 1 atoms on the multivariate
sample dataset from mne. We display a selection of atoms, featuring
heartbeat and eyeblink artifacts, two atoms of evoked responses, and a
non-sinusoidal oscillation.
- 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.
# sampling frequency. The signal will be resampled to match this.
sfreq = 150.
# Define the shape of the dictionary
n_atoms = 40
n_times_atom = int(round(sfreq * 1.0)) # 1000. ms
# Regularization parameter which control sparsity
reg = 0.1
# number of processors for parallel computing
n_jobs = 5
# To accelerate the run time of this example, we split the signal in n_slits.
# The number of splits should actually be the smallest possible to avoid
# introducing border artifacts in the learned atoms and it should be not much
# larger than n_jobs.
n_splits = 10
Next, we define the parameters for multivariate CSC
from alphacsc import GreedyCDL
cdl = GreedyCDL(
# 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',
# apply a temporal window reparametrization
window=True,
# at the end, refit the activations with fixed support and no reg to unbias
unbiased_z_hat=True,
# Initialize the dictionary with random chunk from the data
D_init='chunk',
# rescale the regularization parameter to be a percentage of lambda_max
lmbd_max="scaled",
reg=reg,
# 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-3,
'max_iter': 100000},
# solver for the d-step
solver_d='alternate_adaptive',
solver_d_kwargs={'max_iter': 300},
# sort atoms by explained variances
sort_atoms=True,
# Technical parameters
verbose=1,
random_state=0,
n_jobs=n_jobs)
Load the sample data from MNE-python and select the gradiometer channels. The MNE sample data contains MEG recordings of a subject with visual and auditory stimuli. We load the data using utilities from MNE-python as a Raw object and select the gradiometers from the signal.
import os
import mne
import numpy as np
print("Loading the data...", end='', flush=True)
data_path = mne.datasets.sample.data_path()
subjects_dir = os.path.join(data_path, "subjects")
data_dir = os.path.join(data_path, 'MEG', 'sample')
file_name = os.path.join(data_dir, 'sample_audvis_raw.fif')
raw = mne.io.read_raw_fif(file_name, preload=True, verbose=False)
raw.pick_types(meg='grad', eeg=False, eog=False, stim=True)
print('done')
Out:
Loading the data...Using default location ~/mne_data for sample...
0%| | 0.00/1.65G [00:00<?, ?B/s]
0%| | 4.95M/1.65G [00:00<00:33, 49.5MB/s]
1%|▎ | 11.5M/1.65G [00:00<00:27, 58.9MB/s]
1%|▍ | 18.3M/1.65G [00:00<00:25, 62.9MB/s]
2%|▌ | 25.1M/1.65G [00:00<00:25, 64.9MB/s]
2%|▋ | 31.9M/1.65G [00:00<00:24, 66.2MB/s]
2%|▊ | 38.8M/1.65G [00:00<00:24, 67.0MB/s]
3%|█ | 45.7M/1.65G [00:00<00:23, 67.6MB/s]
3%|█▏ | 52.5M/1.65G [00:00<00:23, 67.9MB/s]
4%|█▎ | 59.4M/1.65G [00:00<00:23, 68.2MB/s]
4%|█▍ | 66.3M/1.65G [00:01<00:23, 68.3MB/s]
4%|█▋ | 73.2M/1.65G [00:01<00:23, 68.5MB/s]
5%|█▊ | 80.0M/1.65G [00:01<00:22, 68.5MB/s]
5%|█▉ | 86.9M/1.65G [00:01<00:22, 68.5MB/s]
6%|██ | 93.8M/1.65G [00:01<00:22, 68.6MB/s]
6%|██▎ | 101M/1.65G [00:01<00:22, 68.5MB/s]
7%|██▍ | 107M/1.65G [00:01<00:22, 68.5MB/s]
7%|██▋ | 114M/1.65G [00:01<00:22, 68.5MB/s]
7%|██▊ | 121M/1.65G [00:01<00:22, 68.3MB/s]
8%|██▉ | 128M/1.65G [00:01<00:22, 68.2MB/s]
8%|███ | 135M/1.65G [00:02<00:22, 68.3MB/s]
9%|███▎ | 142M/1.65G [00:02<00:22, 68.2MB/s]
9%|███▍ | 149M/1.65G [00:02<00:22, 68.3MB/s]
9%|███▌ | 155M/1.65G [00:02<00:21, 68.4MB/s]
10%|███▋ | 162M/1.65G [00:02<00:21, 68.5MB/s]
10%|███▉ | 169M/1.65G [00:02<00:21, 68.7MB/s]
11%|████ | 176M/1.65G [00:02<00:21, 68.6MB/s]
11%|████▏ | 183M/1.65G [00:02<00:21, 68.7MB/s]
11%|████▎ | 190M/1.65G [00:02<00:21, 68.7MB/s]
12%|████▌ | 197M/1.65G [00:02<00:21, 68.7MB/s]
12%|████▋ | 204M/1.65G [00:03<00:21, 68.5MB/s]
13%|████▊ | 210M/1.65G [00:03<00:21, 68.6MB/s]
13%|████▉ | 217M/1.65G [00:03<00:20, 68.6MB/s]
14%|█████▏ | 224M/1.65G [00:03<00:20, 68.7MB/s]
14%|█████▎ | 231M/1.65G [00:03<00:20, 68.8MB/s]
14%|█████▍ | 238M/1.65G [00:03<00:20, 68.7MB/s]
15%|█████▋ | 245M/1.65G [00:03<00:20, 68.7MB/s]
15%|█████▊ | 252M/1.65G [00:03<00:20, 68.9MB/s]
16%|█████▉ | 259M/1.65G [00:03<00:20, 68.9MB/s]
16%|██████ | 266M/1.65G [00:03<00:20, 68.8MB/s]
16%|██████▎ | 272M/1.65G [00:04<00:20, 68.7MB/s]
17%|██████▍ | 279M/1.65G [00:04<00:20, 68.6MB/s]
17%|██████▌ | 286M/1.65G [00:04<00:19, 68.6MB/s]
18%|██████▋ | 293M/1.65G [00:04<00:19, 68.6MB/s]
18%|██████▉ | 300M/1.65G [00:04<00:19, 68.5MB/s]
19%|███████ | 307M/1.65G [00:04<00:19, 68.7MB/s]
19%|███████▏ | 314M/1.65G [00:04<00:19, 68.7MB/s]
19%|███████▎ | 321M/1.65G [00:04<00:19, 68.6MB/s]
20%|███████▌ | 327M/1.65G [00:04<00:19, 68.7MB/s]
20%|███████▋ | 334M/1.65G [00:04<00:19, 68.7MB/s]
21%|███████▊ | 341M/1.65G [00:05<00:19, 68.8MB/s]
21%|████████ | 348M/1.65G [00:05<00:22, 57.3MB/s]
21%|████████▏ | 354M/1.65G [00:05<00:23, 56.3MB/s]
22%|████████▎ | 361M/1.65G [00:05<00:21, 59.4MB/s]
22%|████████▍ | 368M/1.65G [00:05<00:20, 62.0MB/s]
23%|████████▌ | 375M/1.65G [00:05<00:20, 63.8MB/s]
23%|████████▊ | 382M/1.65G [00:05<00:19, 65.1MB/s]
23%|████████▉ | 388M/1.65G [00:05<00:19, 66.0MB/s]
24%|█████████ | 395M/1.65G [00:05<00:18, 66.8MB/s]
24%|█████████▏ | 402M/1.65G [00:05<00:18, 67.4MB/s]
25%|█████████▍ | 409M/1.65G [00:06<00:18, 67.2MB/s]
25%|█████████▌ | 416M/1.65G [00:06<00:18, 67.4MB/s]
26%|█████████▋ | 423M/1.65G [00:06<00:18, 67.7MB/s]
26%|█████████▊ | 429M/1.65G [00:06<00:18, 67.8MB/s]
26%|██████████ | 436M/1.65G [00:06<00:17, 68.0MB/s]
27%|██████████▏ | 443M/1.65G [00:06<00:17, 68.1MB/s]
27%|██████████▎ | 450M/1.65G [00:06<00:17, 68.1MB/s]
28%|██████████▌ | 457M/1.65G [00:06<00:17, 68.2MB/s]
28%|██████████▋ | 464M/1.65G [00:06<00:17, 68.3MB/s]
28%|██████████▊ | 470M/1.65G [00:06<00:17, 68.4MB/s]
29%|██████████▉ | 477M/1.65G [00:07<00:17, 68.5MB/s]
29%|███████████▏ | 484M/1.65G [00:07<00:17, 68.6MB/s]
30%|███████████▎ | 491M/1.65G [00:07<00:16, 68.6MB/s]
30%|███████████▍ | 498M/1.65G [00:07<00:16, 68.6MB/s]
31%|███████████▌ | 505M/1.65G [00:07<00:16, 68.6MB/s]
31%|███████████▊ | 512M/1.65G [00:07<00:16, 68.6MB/s]
31%|███████████▉ | 519M/1.65G [00:07<00:16, 68.6MB/s]
32%|████████████ | 525M/1.65G [00:07<00:16, 68.4MB/s]
32%|████████████▏ | 532M/1.65G [00:07<00:16, 68.3MB/s]
33%|████████████▍ | 539M/1.65G [00:07<00:16, 68.3MB/s]
33%|████████████▌ | 546M/1.65G [00:08<00:16, 68.5MB/s]
33%|████████████▋ | 553M/1.65G [00:08<00:16, 68.5MB/s]
34%|████████████▊ | 560M/1.65G [00:08<00:15, 68.6MB/s]
34%|█████████████ | 567M/1.65G [00:08<00:15, 68.5MB/s]
35%|█████████████▏ | 573M/1.65G [00:08<00:15, 68.6MB/s]
35%|█████████████▎ | 580M/1.65G [00:08<00:15, 68.6MB/s]
36%|█████████████▍ | 587M/1.65G [00:08<00:15, 68.5MB/s]
36%|█████████████▋ | 594M/1.65G [00:08<00:15, 68.6MB/s]
36%|█████████████▊ | 601M/1.65G [00:08<00:15, 68.8MB/s]
37%|█████████████▉ | 608M/1.65G [00:08<00:15, 68.8MB/s]
37%|██████████████▏ | 615M/1.65G [00:09<00:15, 68.8MB/s]
38%|██████████████▎ | 622M/1.65G [00:09<00:14, 68.8MB/s]
38%|██████████████▍ | 628M/1.65G [00:09<00:14, 68.8MB/s]
38%|██████████████▌ | 635M/1.65G [00:09<00:14, 68.8MB/s]
39%|██████████████▊ | 642M/1.65G [00:09<00:14, 68.7MB/s]
39%|██████████████▉ | 649M/1.65G [00:09<00:14, 68.7MB/s]
40%|███████████████ | 656M/1.65G [00:09<00:14, 68.7MB/s]
40%|███████████████▏ | 663M/1.65G [00:09<00:14, 68.5MB/s]
41%|███████████████▍ | 670M/1.65G [00:09<00:14, 68.4MB/s]
41%|███████████████▌ | 677M/1.65G [00:09<00:14, 68.4MB/s]
41%|███████████████▋ | 683M/1.65G [00:10<00:14, 68.3MB/s]
42%|███████████████▊ | 690M/1.65G [00:10<00:14, 68.4MB/s]
42%|████████████████ | 697M/1.65G [00:10<00:13, 68.4MB/s]
43%|████████████████▏ | 704M/1.65G [00:10<00:13, 68.5MB/s]
43%|████████████████▎ | 711M/1.65G [00:10<00:13, 68.5MB/s]
43%|████████████████▌ | 718M/1.65G [00:10<00:13, 68.5MB/s]
44%|████████████████▋ | 725M/1.65G [00:10<00:13, 68.5MB/s]
44%|████████████████▊ | 731M/1.65G [00:10<00:13, 68.5MB/s]
45%|████████████████▉ | 738M/1.65G [00:10<00:16, 55.1MB/s]
45%|█████████████████▏ | 745M/1.65G [00:11<00:15, 57.8MB/s]
45%|█████████████████▎ | 752M/1.65G [00:11<00:14, 60.5MB/s]
46%|█████████████████▍ | 759M/1.65G [00:11<00:14, 62.8MB/s]
46%|█████████████████▌ | 765M/1.65G [00:11<00:13, 64.4MB/s]
47%|█████████████████▊ | 772M/1.65G [00:11<00:13, 65.7MB/s]
47%|█████████████████▉ | 779M/1.65G [00:11<00:13, 66.5MB/s]
48%|██████████████████ | 786M/1.65G [00:11<00:12, 67.1MB/s]
48%|██████████████████▏ | 793M/1.65G [00:11<00:12, 67.7MB/s]
48%|██████████████████▍ | 800M/1.65G [00:11<00:12, 67.9MB/s]
49%|██████████████████▌ | 807M/1.65G [00:11<00:12, 68.1MB/s]
49%|██████████████████▋ | 813M/1.65G [00:12<00:12, 68.3MB/s]
50%|██████████████████▊ | 820M/1.65G [00:12<00:12, 68.5MB/s]
50%|███████████████████ | 827M/1.65G [00:12<00:12, 68.6MB/s]
50%|███████████████████▏ | 834M/1.65G [00:12<00:11, 68.6MB/s]
51%|███████████████████▎ | 841M/1.65G [00:12<00:11, 68.5MB/s]
51%|███████████████████▍ | 848M/1.65G [00:12<00:11, 68.5MB/s]
52%|███████████████████▋ | 855M/1.65G [00:12<00:11, 68.1MB/s]
52%|███████████████████▊ | 862M/1.65G [00:12<00:11, 68.2MB/s]
53%|███████████████████▉ | 868M/1.65G [00:12<00:11, 68.2MB/s]
53%|████████████████████ | 875M/1.65G [00:12<00:11, 68.3MB/s]
53%|████████████████████▎ | 882M/1.65G [00:13<00:11, 68.3MB/s]
54%|████████████████████▍ | 889M/1.65G [00:13<00:11, 68.5MB/s]
54%|████████████████████▌ | 896M/1.65G [00:13<00:11, 68.3MB/s]
55%|████████████████████▊ | 903M/1.65G [00:13<00:11, 68.1MB/s]
55%|████████████████████▉ | 909M/1.65G [00:13<00:10, 68.1MB/s]
55%|█████████████████████ | 916M/1.65G [00:13<00:10, 68.1MB/s]
56%|█████████████████████▏ | 923M/1.65G [00:13<00:10, 68.0MB/s]
56%|█████████████████████▍ | 930M/1.65G [00:13<00:10, 67.8MB/s]
57%|█████████████████████▌ | 937M/1.65G [00:13<00:10, 67.8MB/s]
57%|█████████████████████▋ | 944M/1.65G [00:13<00:10, 68.1MB/s]
57%|█████████████████████▊ | 950M/1.65G [00:14<00:10, 68.1MB/s]
58%|██████████████████████ | 957M/1.65G [00:14<00:10, 68.1MB/s]
58%|██████████████████████▏ | 964M/1.65G [00:14<00:10, 68.1MB/s]
59%|██████████████████████▎ | 971M/1.65G [00:14<00:10, 68.1MB/s]
59%|██████████████████████▍ | 978M/1.65G [00:14<00:09, 68.1MB/s]
60%|██████████████████████▋ | 984M/1.65G [00:14<00:09, 68.0MB/s]
60%|██████████████████████▊ | 991M/1.65G [00:14<00:09, 68.2MB/s]
60%|██████████████████████▉ | 998M/1.65G [00:14<00:09, 68.4MB/s]
61%|██████████████████████▍ | 1.01G/1.65G [00:14<00:09, 68.4MB/s]
61%|██████████████████████▋ | 1.01G/1.65G [00:14<00:09, 68.5MB/s]
62%|██████████████████████▊ | 1.02G/1.65G [00:15<00:09, 68.4MB/s]
62%|██████████████████████▉ | 1.03G/1.65G [00:15<00:09, 68.5MB/s]
62%|███████████████████████ | 1.03G/1.65G [00:15<00:09, 68.5MB/s]
63%|███████████████████████▎ | 1.04G/1.65G [00:15<00:08, 68.6MB/s]
63%|███████████████████████▍ | 1.05G/1.65G [00:15<00:08, 68.5MB/s]
64%|███████████████████████▌ | 1.05G/1.65G [00:15<00:08, 68.5MB/s]
64%|███████████████████████▋ | 1.06G/1.65G [00:15<00:08, 68.5MB/s]
65%|███████████████████████▉ | 1.07G/1.65G [00:15<00:08, 68.4MB/s]
65%|████████████████████████ | 1.07G/1.65G [00:15<00:08, 68.5MB/s]
65%|████████████████████████▏ | 1.08G/1.65G [00:15<00:08, 68.3MB/s]
66%|████████████████████████▎ | 1.09G/1.65G [00:16<00:08, 68.0MB/s]
66%|████████████████████████▍ | 1.09G/1.65G [00:16<00:08, 67.8MB/s]
67%|████████████████████████▋ | 1.10G/1.65G [00:16<00:08, 67.6MB/s]
67%|████████████████████████▊ | 1.11G/1.65G [00:16<00:08, 67.6MB/s]
67%|████████████████████████▉ | 1.11G/1.65G [00:16<00:07, 67.5MB/s]
68%|█████████████████████████ | 1.12G/1.65G [00:16<00:07, 67.5MB/s]
68%|█████████████████████████▎ | 1.13G/1.65G [00:16<00:07, 67.4MB/s]
69%|█████████████████████████▍ | 1.13G/1.65G [00:16<00:07, 67.4MB/s]
69%|█████████████████████████▌ | 1.14G/1.65G [00:16<00:09, 55.7MB/s]
69%|█████████████████████████▋ | 1.15G/1.65G [00:17<00:09, 55.0MB/s]
70%|█████████████████████████▊ | 1.15G/1.65G [00:17<00:08, 58.2MB/s]
70%|█████████████████████████▉ | 1.16G/1.65G [00:17<00:08, 60.8MB/s]
71%|██████████████████████████▏ | 1.17G/1.65G [00:17<00:07, 62.6MB/s]
71%|██████████████████████████▎ | 1.17G/1.65G [00:17<00:07, 63.9MB/s]
71%|██████████████████████████▍ | 1.18G/1.65G [00:17<00:07, 64.9MB/s]
72%|██████████████████████████▌ | 1.19G/1.65G [00:17<00:07, 65.7MB/s]
72%|██████████████████████████▋ | 1.19G/1.65G [00:17<00:06, 66.2MB/s]
73%|██████████████████████████▉ | 1.20G/1.65G [00:17<00:06, 66.4MB/s]
73%|███████████████████████████ | 1.21G/1.65G [00:17<00:06, 66.7MB/s]
73%|███████████████████████████▏ | 1.21G/1.65G [00:18<00:06, 66.8MB/s]
74%|███████████████████████████▎ | 1.22G/1.65G [00:18<00:06, 66.9MB/s]
74%|███████████████████████████▍ | 1.23G/1.65G [00:18<00:06, 67.1MB/s]
75%|███████████████████████████▋ | 1.23G/1.65G [00:18<00:06, 67.1MB/s]
75%|███████████████████████████▊ | 1.24G/1.65G [00:18<00:06, 67.2MB/s]
76%|███████████████████████████▉ | 1.25G/1.65G [00:18<00:06, 67.2MB/s]
76%|████████████████████████████ | 1.26G/1.65G [00:18<00:05, 67.4MB/s]
76%|████████████████████████████▏ | 1.26G/1.65G [00:18<00:05, 67.5MB/s]
77%|████████████████████████████▍ | 1.27G/1.65G [00:18<00:05, 67.5MB/s]
77%|████████████████████████████▌ | 1.28G/1.65G [00:18<00:05, 67.5MB/s]
78%|████████████████████████████▋ | 1.28G/1.65G [00:19<00:05, 67.4MB/s]
78%|████████████████████████████▊ | 1.29G/1.65G [00:19<00:05, 67.4MB/s]
78%|█████████████████████████████ | 1.30G/1.65G [00:19<00:05, 67.3MB/s]
79%|█████████████████████████████▏ | 1.30G/1.65G [00:19<00:05, 67.1MB/s]
79%|█████████████████████████████▎ | 1.31G/1.65G [00:19<00:05, 67.2MB/s]
80%|█████████████████████████████▍ | 1.32G/1.65G [00:19<00:05, 67.2MB/s]
80%|█████████████████████████████▌ | 1.32G/1.65G [00:19<00:04, 67.2MB/s]
80%|█████████████████████████████▊ | 1.33G/1.65G [00:19<00:04, 67.2MB/s]
81%|█████████████████████████████▉ | 1.34G/1.65G [00:19<00:04, 67.2MB/s]
81%|██████████████████████████████ | 1.34G/1.65G [00:19<00:04, 67.2MB/s]
82%|██████████████████████████████▏ | 1.35G/1.65G [00:20<00:04, 67.3MB/s]
82%|██████████████████████████████▎ | 1.36G/1.65G [00:20<00:04, 67.4MB/s]
82%|██████████████████████████████▌ | 1.36G/1.65G [00:20<00:04, 67.4MB/s]
83%|██████████████████████████████▋ | 1.37G/1.65G [00:20<00:04, 67.4MB/s]
83%|██████████████████████████████▊ | 1.38G/1.65G [00:20<00:04, 67.2MB/s]
84%|██████████████████████████████▉ | 1.38G/1.65G [00:20<00:04, 67.1MB/s]
84%|███████████████████████████████ | 1.39G/1.65G [00:20<00:03, 67.2MB/s]
85%|███████████████████████████████▎ | 1.40G/1.65G [00:20<00:03, 67.3MB/s]
85%|███████████████████████████████▍ | 1.40G/1.65G [00:20<00:03, 67.5MB/s]
85%|███████████████████████████████▌ | 1.41G/1.65G [00:20<00:03, 67.4MB/s]
86%|███████████████████████████████▋ | 1.42G/1.65G [00:21<00:03, 67.6MB/s]
86%|███████████████████████████████▊ | 1.42G/1.65G [00:21<00:03, 67.5MB/s]
87%|████████████████████████████████ | 1.43G/1.65G [00:21<00:03, 67.7MB/s]
87%|████████████████████████████████▏ | 1.44G/1.65G [00:21<00:03, 67.8MB/s]
87%|████████████████████████████████▎ | 1.44G/1.65G [00:21<00:03, 67.8MB/s]
88%|████████████████████████████████▍ | 1.45G/1.65G [00:21<00:03, 67.1MB/s]
88%|████████████████████████████████▋ | 1.46G/1.65G [00:21<00:02, 67.1MB/s]
89%|████████████████████████████████▊ | 1.46G/1.65G [00:21<00:02, 67.2MB/s]
89%|████████████████████████████████▉ | 1.47G/1.65G [00:21<00:02, 67.4MB/s]
89%|█████████████████████████████████ | 1.48G/1.65G [00:21<00:02, 67.5MB/s]
90%|█████████████████████████████████▏ | 1.48G/1.65G [00:22<00:02, 67.4MB/s]
90%|█████████████████████████████████▍ | 1.49G/1.65G [00:22<00:02, 67.5MB/s]
91%|█████████████████████████████████▌ | 1.50G/1.65G [00:22<00:02, 67.6MB/s]
91%|█████████████████████████████████▋ | 1.50G/1.65G [00:22<00:02, 67.6MB/s]
91%|█████████████████████████████████▊ | 1.51G/1.65G [00:22<00:02, 67.7MB/s]
92%|█████████████████████████████████▉ | 1.52G/1.65G [00:22<00:01, 67.6MB/s]
92%|██████████████████████████████████▏ | 1.53G/1.65G [00:22<00:01, 67.4MB/s]
93%|██████████████████████████████████▎ | 1.53G/1.65G [00:22<00:01, 67.0MB/s]
93%|██████████████████████████████████▍ | 1.54G/1.65G [00:22<00:01, 67.0MB/s]
94%|██████████████████████████████████▌ | 1.55G/1.65G [00:22<00:01, 67.3MB/s]
94%|██████████████████████████████████▊ | 1.55G/1.65G [00:23<00:01, 67.3MB/s]
94%|██████████████████████████████████▉ | 1.56G/1.65G [00:23<00:01, 67.5MB/s]
95%|███████████████████████████████████ | 1.57G/1.65G [00:23<00:01, 67.5MB/s]
95%|███████████████████████████████████▏ | 1.57G/1.65G [00:23<00:01, 67.6MB/s]
96%|███████████████████████████████████▎ | 1.58G/1.65G [00:23<00:01, 67.7MB/s]
96%|███████████████████████████████████▌ | 1.59G/1.65G [00:23<00:00, 67.8MB/s]
96%|███████████████████████████████████▋ | 1.59G/1.65G [00:23<00:00, 67.9MB/s]
97%|███████████████████████████████████▊ | 1.60G/1.65G [00:23<00:00, 67.7MB/s]
97%|███████████████████████████████████▉ | 1.61G/1.65G [00:23<00:00, 67.6MB/s]
98%|████████████████████████████████████ | 1.61G/1.65G [00:23<00:00, 67.6MB/s]
98%|████████████████████████████████████▎| 1.62G/1.65G [00:24<00:00, 67.5MB/s]
98%|████████████████████████████████████▍| 1.63G/1.65G [00:24<00:00, 67.5MB/s]
99%|████████████████████████████████████▌| 1.63G/1.65G [00:24<00:00, 67.6MB/s]
99%|████████████████████████████████████▋| 1.64G/1.65G [00:24<00:00, 67.6MB/s]
100%|████████████████████████████████████▉| 1.65G/1.65G [00:24<00:00, 67.6MB/s]
0%| | 0.00/1.65G [00:00<?, ?B/s]
100%|█████████████████████████████████████| 1.65G/1.65G [00:00<00:00, 3.09TB/s]
Removing projector <Projection | PCA-v1, active : False, n_channels : 102>
Removing projector <Projection | PCA-v2, active : False, n_channels : 102>
Removing projector <Projection | PCA-v3, active : False, n_channels : 102>
done
Then, we remove the powerline artifacts and high-pass filter to remove the drift which can impact the CSC technique. The signal is also resampled to 150 Hz to reduce the computationnal burden.
print("Preprocessing the data...", end='', flush=True)
raw.notch_filter(np.arange(60, 181, 60), n_jobs=n_jobs, verbose=False)
raw.filter(2, None, n_jobs=n_jobs, verbose=False)
raw = raw.resample(sfreq, npad='auto', n_jobs=n_jobs, verbose=False)
print('done')
Out:
Preprocessing the data.../github/workspace/examples/multicsc/plot_sample_evoked_response.py:115: RuntimeWarning: Resampling of the stim channels caused event information to become unreliable. Consider finding events on the original data and passing the event matrix as a parameter.
raw = raw.resample(sfreq, npad='auto', n_jobs=n_jobs, verbose=False)
done
Load the data as an array and split it in chunks to allow parallel processing during the model fit. Each split is considered as independent. To reduce the impact of border artifacts, we use apply_window=True which scales down the border of each split with a tukey window.
Fit the model and learn rank1 atoms
cdl.fit(X_split)
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(
................................................+
......
[GreedyCDL] Converged after 56 iteration, (dz, du) = 9.879e-05, 9.739e-05
[GreedyCDL] Fit in 4937.0s
<alphacsc.convolutional_dictionary_learning.GreedyCDL object at 0x7f587feca340>
Then we call the transform method, which returns the sparse codes associated with X, without changing the dictionary learned during the fit. Note that we transform on the unsplit data so that the sparse codes reflect the original data and not the windowed data.
Out:
Refitting the activation to avoid amplitude bias...done
Display a selection of atoms#
We recognize a heartbeat artifact, an eyeblink artifact, two atoms of evoked responses, and a non-sinusoidal oscillation.
import matplotlib.pyplot as plt
# preselected atoms of interest
plotted_atoms = [0, 1, 2, 6, 4]
n_plots = 3 # number of plots by atom
n_columns = min(6, len(plotted_atoms))
split = int(np.ceil(len(plotted_atoms) / n_columns))
figsize = (4 * n_columns, 3 * n_plots * split)
fig, axes = plt.subplots(n_plots * split, n_columns, figsize=figsize)
for ii, kk in enumerate(plotted_atoms):
# Select the axes to display the current atom
print("\rDisplaying {}-th atom".format(kk), end='', flush=True)
i_row, i_col = ii // n_columns, ii % n_columns
it_axes = iter(axes[i_row * n_plots:(i_row + 1) * n_plots, i_col])
# Select the current atom
u_k = cdl.u_hat_[kk]
v_k = cdl.v_hat_[kk]
# Plot the spatial map of the atom using mne topomap
ax = next(it_axes)
mne.viz.plot_topomap(u_k, info, axes=ax, show=False)
ax.set(title="Spatial pattern %d" % (kk, ))
# Plot the temporal pattern of the atom
ax = next(it_axes)
t = np.arange(n_times_atom) / sfreq
ax.plot(t, v_k)
ax.set_xlim(0, n_times_atom / sfreq)
ax.set(xlabel='Time (sec)', title="Temporal pattern %d" % kk)
# Plot the power spectral density (PSD)
ax = next(it_axes)
psd = np.abs(np.fft.rfft(v_k, n=256)) ** 2
frequencies = np.linspace(0, sfreq / 2.0, len(psd))
ax.semilogy(frequencies, psd, label='PSD', color='k')
ax.set(xlabel='Frequencies (Hz)', title="Power spectral density %d" % kk)
ax.grid(True)
ax.set_xlim(0, 30)
ax.set_ylim(1e-4, 1e2)
ax.legend()
print("\rDisplayed {} atoms".format(len(plotted_atoms)).rjust(40))
fig.tight_layout()
Out:
Displaying 0-th atom
Displaying 1-th atom
Displaying 2-th atom
Displaying 6-th atom
Displaying 4-th atom
Displayed 5 atoms
Display the evoked reconstructed envelope#
The MNE sample data contains data for auditory (event_id=1 and 2) and visual stimuli (event_id=3 and 4). We extract the events now so that we can later identify the atoms related to different events. Note that the convolutional sparse coding method does not need to know the events for learning atoms.
event_id = [1, 2, 3, 4]
events = mne.find_events(raw, stim_channel='STI 014')
events = mne.pick_events(events, include=event_id)
events[:, 0] -= raw.first_samp
Out:
319 events found
Event IDs: [ 1 2 3 4 5 32]
For each atom (columns), and for each event (rows), we compute the envelope of the reconstructed signal, align it with respect to the event onsets, and take the average. For some atoms, the activations are correlated with the events, leading to a large evoked envelope. The gray area corresponds to not statistically significant values, computing with sampling.
from alphacsc.utils.signal import fast_hilbert
from alphacsc.viz.epoch import plot_evoked_surrogates
from alphacsc.utils.convolution import construct_X_multi
# time window around the events. Note that for the sample datasets, the time
# inter-event is around 0.5s
t_lim = (-0.1, 0.5)
n_plots = len(event_id)
n_columns = min(6, len(plotted_atoms))
split = int(np.ceil(len(plotted_atoms) / n_columns))
figsize = (4 * n_columns, 3 * n_plots * split)
fig, axes = plt.subplots(n_plots * split, n_columns, figsize=figsize)
for ii, kk in enumerate(plotted_atoms):
# Select the axes to display the current atom
print("\rDisplaying {}-th atom envelope".format(kk), end='', flush=True)
i_row, i_col = ii // n_columns, ii % n_columns
it_axes = iter(axes[i_row * n_plots:(i_row + 1) * n_plots, i_col])
# Select the current atom
v_k = cdl.v_hat_[kk]
v_k_1 = np.r_[[1], v_k][None]
z_k = z_hat[:, kk:kk + 1]
X_k = construct_X_multi(z_k, v_k_1, n_channels=1)[0, 0]
# compute the 'envelope' of the reconstructed signal X_k
correlation = np.abs(fast_hilbert(X_k))
# loop over all events IDs
for this_event_id in event_id:
this_events = events[events[:, 2] == this_event_id]
# plotting function
ax = next(it_axes)
this_info = info.copy()
event_info = dict(event_id = this_event_id, events=events)
this_info['temp'] = event_info
plot_evoked_surrogates(correlation, info=this_info, t_lim=t_lim, ax=ax,
n_jobs=n_jobs, label='event %d' % this_event_id)
ax.set(xlabel='Time (sec)', title="Evoked envelope %d" % kk)
print("\rDisplayed {} atoms".format(len(plotted_atoms)).rjust(40))
fig.tight_layout()
Out:
Displaying 0-th atom envelopeUsing data from preloaded Raw for 72 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 73 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 73 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 70 events and 91 original time points ...
0 bad epochs dropped
Displaying 1-th atom envelopeUsing data from preloaded Raw for 72 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 73 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 73 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 70 events and 91 original time points ...
0 bad epochs dropped
Displaying 2-th atom envelopeUsing data from preloaded Raw for 72 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 73 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 73 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 70 events and 91 original time points ...
0 bad epochs dropped
Displaying 6-th atom envelopeUsing data from preloaded Raw for 72 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 73 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 73 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 70 events and 91 original time points ...
0 bad epochs dropped
Displaying 4-th atom envelopeUsing data from preloaded Raw for 72 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 73 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 73 events and 91 original time points ...
0 bad epochs dropped
Using data from preloaded Raw for 70 events and 91 original time points ...
0 bad epochs dropped
Displayed 5 atoms
Display the equivalent dipole for a learned topomap#
Finally, let us fit a dipole to one of the atoms. To fit a dipole, we need the following:
BEM solution: Obtained by running the cortical reconstruction pipeline of Freesurfer and describes the conductivity of different tissues in the head.
Trans: An affine transformation matrix needed to bring the data from sensor space to head space. This is usually done by coregistering the fiducials with the MRI.
Noise covariance matrix: To whiten the data so that the assumption of Gaussian noise model with identity covariance matrix is satisfied.
We recommend users to consult the MNE documentation for further information.
subjects_dir = os.path.join(data_path, 'subjects')
fname_bem = os.path.join(subjects_dir, 'sample', 'bem',
'sample-5120-bem-sol.fif')
fname_trans = os.path.join(data_path, 'MEG', 'sample',
'sample_audvis_raw-trans.fif')
fname_cov = os.path.join(data_path, 'MEG', 'sample', 'sample_audvis-cov.fif')
Let us construct an evoked object for MNE with the spatial pattern of the atoms.
evoked = mne.EvokedArray(cdl.u_hat_.T, info)
Fit a dipole to each of the atoms.
dip = mne.fit_dipole(evoked, fname_cov, fname_bem, fname_trans,
n_jobs=n_jobs, verbose=False)[0]
Plot the dipole fit from the 3rd atom, linked to mu-wave and display the goodness of fit.
atom_dipole_idx = 4
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(10, 4))
# Display the dipole fit
ax = fig.add_subplot(1, 3, 1, projection='3d')
dip.plot_locations(fname_trans, 'sample', subjects_dir, idx=atom_dipole_idx,
ax=ax)
ax.set_title('Atom #{} (GOF {:.2f}%)'.format(atom_dipole_idx,
dip.gof[atom_dipole_idx]))
# Plot the spatial map
ax = fig.add_subplot(1, 3, 2)
mne.viz.plot_topomap(cdl.u_hat_[atom_dipole_idx], info, axes=ax)
# Plot the temporal atom
ax = fig.add_subplot(1, 3, 3)
t = np.arange(n_times_atom) / sfreq
ax.plot(t, cdl.v_hat_[atom_dipole_idx])
ax.set_xlim(0, n_times_atom / sfreq)
ax.set(xlabel='Time (sec)', title="Temporal pattern {}"
.format(atom_dipole_idx))
fig.suptitle('')
fig.tight_layout()
Total running time of the script: ( 91 minutes 23.318 seconds)