Note
Go to the end to download the full example code.
Detect bad sensors using RANSAC#
This example demonstrates how to use RANSAC [1] from the PREP pipeline to
detect bad sensors and repair them. Note that this implementation in
autoreject [2] is an extension of the original implementation and
works for MEG sensors as well.
References#
# Author: Mainak Jas <mainak.jas@telecom-paristech.fr>
# License: BSD-3-Clause
For the purposes of this example, we shall use the MNE sample dataset. Therefore, let us make some MNE related imports.
import mne
from mne import io
from mne import Epochs
from mne.datasets import sample
Let us now read in the raw fif file for MNE sample dataset.
data_path = sample.data_path()
meg_path = data_path / 'MEG' / 'sample'
raw_fname = meg_path / 'sample_audvis_filt-0-40_raw.fif'
raw = io.read_raw_fif(raw_fname, preload=True)
Opening raw data file /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis_filt-0-40_raw.fif...
Read a total of 4 projection items:
PCA-v1 (1 x 102) idle
PCA-v2 (1 x 102) idle
PCA-v3 (1 x 102) idle
Average EEG reference (1 x 60) idle
Range : 6450 ... 48149 = 42.956 ... 320.665 secs
Ready.
Reading 0 ... 41699 = 0.000 ... 277.709 secs...
We can then read in the events
event_fname = meg_path / 'sample_audvis_filt-0-40_raw-eve.fif'
event_id = {'Auditory/Left': 1}
tmin, tmax = -0.2, 0.5
events = mne.read_events(event_fname)
And pick MEG channels for repairing. Currently, autoreject can repair
only one channel type at a time.
raw.info['bads'] = []
Now, we can create epochs. The reject params will be set to None
because we do not want epochs to be dropped when instantiating
mne.Epochs.
raw.del_proj() # remove proj, don't proj while interpolating
epochs = Epochs(raw, events, event_id, tmin, tmax,
baseline=(None, 0), reject=None,
verbose=False, detrend=0, preload=True)
picks = mne.pick_types(epochs.info, meg='grad', eeg=False,
stim=False, eog=False,
include=[], exclude=[])
We import Ransac and run the familiar fit_transform method.
0%| | interpolating channels : 0/50 [00:00<?, ?it/s]
2%|▏ | interpolating channels : 1/50 [00:01<01:22, 1.69s/it]
4%|▍ | interpolating channels : 2/50 [00:03<01:21, 1.71s/it]
6%|▌ | interpolating channels : 3/50 [00:05<01:20, 1.71s/it]
8%|▊ | interpolating channels : 4/50 [00:06<01:18, 1.71s/it]
10%|█ | interpolating channels : 5/50 [00:08<01:16, 1.70s/it]
12%|█▏ | interpolating channels : 6/50 [00:10<01:14, 1.70s/it]
14%|█▍ | interpolating channels : 7/50 [00:11<01:13, 1.70s/it]
16%|█▌ | interpolating channels : 8/50 [00:13<01:11, 1.71s/it]
18%|█▊ | interpolating channels : 9/50 [00:15<01:09, 1.71s/it]
20%|██ | interpolating channels : 10/50 [00:17<01:08, 1.71s/it]
22%|██▏ | interpolating channels : 11/50 [00:18<01:06, 1.70s/it]
24%|██▍ | interpolating channels : 12/50 [00:20<01:04, 1.71s/it]
26%|██▌ | interpolating channels : 13/50 [00:22<01:03, 1.71s/it]
28%|██▊ | interpolating channels : 14/50 [00:23<01:01, 1.71s/it]
30%|███ | interpolating channels : 15/50 [00:25<01:00, 1.72s/it]
32%|███▏ | interpolating channels : 16/50 [00:27<00:58, 1.72s/it]
34%|███▍ | interpolating channels : 17/50 [00:29<00:56, 1.72s/it]
36%|███▌ | interpolating channels : 18/50 [00:30<00:55, 1.73s/it]
38%|███▊ | interpolating channels : 19/50 [00:32<00:53, 1.73s/it]
40%|████ | interpolating channels : 20/50 [00:34<00:51, 1.73s/it]
42%|████▏ | interpolating channels : 21/50 [00:36<00:50, 1.74s/it]
44%|████▍ | interpolating channels : 22/50 [00:38<00:48, 1.74s/it]
46%|████▌ | interpolating channels : 23/50 [00:39<00:47, 1.74s/it]
48%|████▊ | interpolating channels : 24/50 [00:41<00:45, 1.74s/it]
50%|█████ | interpolating channels : 25/50 [00:43<00:43, 1.75s/it]
52%|█████▏ | interpolating channels : 26/50 [00:45<00:42, 1.75s/it]
54%|█████▍ | interpolating channels : 27/50 [00:47<00:40, 1.76s/it]
56%|█████▌ | interpolating channels : 28/50 [00:48<00:38, 1.76s/it]
58%|█████▊ | interpolating channels : 29/50 [00:50<00:37, 1.77s/it]
60%|██████ | interpolating channels : 30/50 [00:52<00:35, 1.77s/it]
62%|██████▏ | interpolating channels : 31/50 [00:54<00:33, 1.78s/it]
64%|██████▍ | interpolating channels : 32/50 [00:56<00:31, 1.78s/it]
66%|██████▌ | interpolating channels : 33/50 [00:57<00:30, 1.78s/it]
68%|██████▊ | interpolating channels : 34/50 [00:59<00:28, 1.78s/it]
70%|███████ | interpolating channels : 35/50 [01:01<00:26, 1.78s/it]
72%|███████▏ | interpolating channels : 36/50 [01:03<00:24, 1.78s/it]
74%|███████▍ | interpolating channels : 37/50 [01:05<00:23, 1.79s/it]
76%|███████▌ | interpolating channels : 38/50 [01:07<00:21, 1.79s/it]
78%|███████▊ | interpolating channels : 39/50 [01:08<00:19, 1.79s/it]
80%|████████ | interpolating channels : 40/50 [01:10<00:17, 1.79s/it]
82%|████████▏ | interpolating channels : 41/50 [01:12<00:16, 1.79s/it]
84%|████████▍ | interpolating channels : 42/50 [01:14<00:14, 1.78s/it]
86%|████████▌ | interpolating channels : 43/50 [01:15<00:12, 1.79s/it]
88%|████████▊ | interpolating channels : 44/50 [01:17<00:10, 1.79s/it]
90%|█████████ | interpolating channels : 45/50 [01:19<00:08, 1.79s/it]
92%|█████████▏| interpolating channels : 46/50 [01:21<00:07, 1.79s/it]
94%|█████████▍| interpolating channels : 47/50 [01:23<00:05, 1.79s/it]
96%|█████████▌| interpolating channels : 48/50 [01:24<00:03, 1.79s/it]
98%|█████████▊| interpolating channels : 49/50 [01:26<00:01, 1.79s/it]
100%|██████████| interpolating channels : 50/50 [01:28<00:00, 1.79s/it]
100%|██████████| interpolating channels : 50/50 [01:28<00:00, 1.77s/it]
[Parallel(n_jobs=1)]: Done 1 tasks | elapsed: 1.5min
[Parallel(n_jobs=1)]: Done 1 tasks | elapsed: 1.5min
0%| | Iterating epochs : 0/72 [00:00<?, ?it/s]
1%|▏ | Iterating epochs : 1/72 [00:00<00:06, 10.93it/s]
3%|▎ | Iterating epochs : 2/72 [00:00<00:04, 15.08it/s]
4%|▍ | Iterating epochs : 3/72 [00:00<00:03, 17.31it/s]
6%|▌ | Iterating epochs : 4/72 [00:00<00:03, 18.35it/s]
7%|▋ | Iterating epochs : 5/72 [00:00<00:03, 19.30it/s]
8%|▊ | Iterating epochs : 6/72 [00:00<00:03, 20.00it/s]
10%|▉ | Iterating epochs : 7/72 [00:00<00:03, 20.54it/s]
11%|█ | Iterating epochs : 8/72 [00:00<00:03, 20.96it/s]
12%|█▎ | Iterating epochs : 9/72 [00:00<00:02, 21.29it/s]
14%|█▍ | Iterating epochs : 10/72 [00:00<00:02, 21.57it/s]
15%|█▌ | Iterating epochs : 11/72 [00:00<00:02, 21.80it/s]
17%|█▋ | Iterating epochs : 12/72 [00:00<00:02, 22.00it/s]
18%|█▊ | Iterating epochs : 13/72 [00:00<00:02, 22.17it/s]
19%|█▉ | Iterating epochs : 14/72 [00:00<00:02, 22.31it/s]
21%|██ | Iterating epochs : 15/72 [00:00<00:02, 22.45it/s]
22%|██▏ | Iterating epochs : 16/72 [00:00<00:02, 22.59it/s]
24%|██▎ | Iterating epochs : 17/72 [00:00<00:02, 22.69it/s]
25%|██▌ | Iterating epochs : 18/72 [00:00<00:02, 22.78it/s]
26%|██▋ | Iterating epochs : 19/72 [00:00<00:02, 22.86it/s]
28%|██▊ | Iterating epochs : 20/72 [00:00<00:02, 22.93it/s]
29%|██▉ | Iterating epochs : 21/72 [00:00<00:02, 23.00it/s]
31%|███ | Iterating epochs : 22/72 [00:00<00:02, 23.07it/s]
32%|███▏ | Iterating epochs : 23/72 [00:01<00:02, 23.12it/s]
33%|███▎ | Iterating epochs : 24/72 [00:01<00:02, 23.16it/s]
35%|███▍ | Iterating epochs : 25/72 [00:01<00:02, 23.18it/s]
36%|███▌ | Iterating epochs : 26/72 [00:01<00:01, 23.21it/s]
38%|███▊ | Iterating epochs : 27/72 [00:01<00:01, 23.25it/s]
39%|███▉ | Iterating epochs : 28/72 [00:01<00:01, 23.28it/s]
40%|████ | Iterating epochs : 29/72 [00:01<00:01, 23.30it/s]
42%|████▏ | Iterating epochs : 30/72 [00:01<00:01, 23.33it/s]
43%|████▎ | Iterating epochs : 31/72 [00:01<00:01, 23.38it/s]
44%|████▍ | Iterating epochs : 32/72 [00:01<00:01, 23.41it/s]
46%|████▌ | Iterating epochs : 33/72 [00:01<00:01, 23.44it/s]
47%|████▋ | Iterating epochs : 34/72 [00:01<00:01, 23.45it/s]
49%|████▊ | Iterating epochs : 35/72 [00:01<00:01, 23.48it/s]
50%|█████ | Iterating epochs : 36/72 [00:01<00:01, 23.51it/s]
51%|█████▏ | Iterating epochs : 37/72 [00:01<00:01, 23.53it/s]
53%|█████▎ | Iterating epochs : 38/72 [00:01<00:01, 23.54it/s]
54%|█████▍ | Iterating epochs : 39/72 [00:01<00:01, 23.55it/s]
56%|█████▌ | Iterating epochs : 40/72 [00:01<00:01, 23.56it/s]
57%|█████▋ | Iterating epochs : 41/72 [00:01<00:01, 23.57it/s]
58%|█████▊ | Iterating epochs : 42/72 [00:01<00:01, 23.58it/s]
60%|█████▉ | Iterating epochs : 43/72 [00:01<00:01, 23.59it/s]
61%|██████ | Iterating epochs : 44/72 [00:01<00:01, 23.59it/s]
62%|██████▎ | Iterating epochs : 45/72 [00:01<00:01, 23.61it/s]
64%|██████▍ | Iterating epochs : 46/72 [00:01<00:01, 23.63it/s]
65%|██████▌ | Iterating epochs : 47/72 [00:02<00:01, 23.65it/s]
67%|██████▋ | Iterating epochs : 48/72 [00:02<00:01, 23.66it/s]
68%|██████▊ | Iterating epochs : 49/72 [00:02<00:00, 23.66it/s]
69%|██████▉ | Iterating epochs : 50/72 [00:02<00:00, 23.67it/s]
71%|███████ | Iterating epochs : 51/72 [00:02<00:00, 23.68it/s]
72%|███████▏ | Iterating epochs : 52/72 [00:02<00:00, 23.70it/s]
74%|███████▎ | Iterating epochs : 53/72 [00:02<00:00, 23.70it/s]
75%|███████▌ | Iterating epochs : 54/72 [00:02<00:00, 23.71it/s]
76%|███████▋ | Iterating epochs : 55/72 [00:02<00:00, 23.71it/s]
78%|███████▊ | Iterating epochs : 56/72 [00:02<00:00, 23.70it/s]
79%|███████▉ | Iterating epochs : 57/72 [00:02<00:00, 23.71it/s]
81%|████████ | Iterating epochs : 58/72 [00:02<00:00, 23.72it/s]
82%|████████▏ | Iterating epochs : 59/72 [00:02<00:00, 23.73it/s]
83%|████████▎ | Iterating epochs : 60/72 [00:02<00:00, 23.73it/s]
85%|████████▍ | Iterating epochs : 61/72 [00:02<00:00, 23.73it/s]
86%|████████▌ | Iterating epochs : 62/72 [00:02<00:00, 23.74it/s]
88%|████████▊ | Iterating epochs : 63/72 [00:02<00:00, 23.72it/s]
89%|████████▉ | Iterating epochs : 64/72 [00:02<00:00, 23.70it/s]
90%|█████████ | Iterating epochs : 65/72 [00:02<00:00, 23.70it/s]
92%|█████████▏| Iterating epochs : 66/72 [00:02<00:00, 23.70it/s]
93%|█████████▎| Iterating epochs : 67/72 [00:02<00:00, 23.70it/s]
94%|█████████▍| Iterating epochs : 68/72 [00:02<00:00, 23.69it/s]
96%|█████████▌| Iterating epochs : 69/72 [00:02<00:00, 23.69it/s]
97%|█████████▋| Iterating epochs : 70/72 [00:02<00:00, 23.68it/s]
99%|█████████▊| Iterating epochs : 71/72 [00:03<00:00, 23.68it/s]
100%|██████████| Iterating epochs : 72/72 [00:03<00:00, 23.68it/s]
100%|██████████| Iterating epochs : 72/72 [00:03<00:00, 23.35it/s]
[Parallel(n_jobs=1)]: Done 1 tasks | elapsed: 3.1s
[Parallel(n_jobs=1)]: Done 1 tasks | elapsed: 3.1s
[Done]
Setting channel interpolation method to {'eeg': 'spline', 'meg': 'MNE'}.
Interpolating bad channels.
Automatic origin fit: head of radius 91.0 mm
Computing dot products for 292 MEG channels...
Computing cross products for 292 → 14 MEG channels...
Preparing the mapping matrix...
Truncating at 80/292 components to omit less than 0.0001 (9.4e-05)
We can also get the list of bad channels computed by Ransac.
print('\n'.join(ransac.bad_chs_))
MEG 0722
MEG 0813
MEG 0812
MEG 0822
MEG 0912
MEG 0942
MEG 0943
MEG 1013
MEG 1023
MEG 1022
MEG 1032
MEG 1422
MEG 1633
MEG 2443
Then we compute the evoked before and after interpolation.
We will manually mark the bad channels just for plotting.
evoked.info['bads'] = ['MEG 2443']
evoked_clean.info['bads'] = ['MEG 2443']
Let us plot the results.
from autoreject.utils import set_matplotlib_defaults # noqa
import matplotlib.pyplot as plt # noqa
set_matplotlib_defaults(plt)
fig, axes = plt.subplots(2, 1, figsize=(6, 6))
for ax in axes:
ax.tick_params(axis='x', which='both', bottom='off', top='off')
ax.tick_params(axis='y', which='both', left='off', right='off')
ylim = dict(grad=(-170, 200))
evoked.pick_types(meg='grad', exclude=[])
evoked.plot(exclude=[], axes=axes[0], ylim=ylim, show=False)
axes[0].set_title('Before RANSAC')
evoked_clean.pick_types(meg='grad', exclude=[])
evoked_clean.plot(exclude=[], axes=axes[1], ylim=ylim)
axes[1].set_title('After RANSAC')
fig.tight_layout()
NOTE: pick_types() is a legacy function. New code should use inst.pick(...).
NOTE: pick_types() is a legacy function. New code should use inst.pick(...).
/home/circleci/project/examples/plot_ransac.py:115: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
fig.tight_layout()
To top things up, we can also visualize the bad sensors for each trial using a heatmap.
ch_names = [epochs.ch_names[ii] for ii in ransac.picks][7::10]
fig, ax = plt.subplots(1, 1, figsize=(12, 6))
ax.imshow(ransac.bad_log, cmap='Reds',
interpolation='nearest')
ax.grid(False)
ax.set_xlabel('Sensors')
ax.set_ylabel('Trials')
plt.setp(ax, xticks=range(7, len(ransac.picks), 10),
xticklabels=ch_names)
plt.setp(ax.get_yticklabels(), rotation=0)
plt.setp(ax.get_xticklabels(), rotation=90)
ax.tick_params(axis=u'both', which=u'both', length=0)
fig.tight_layout(rect=[None, None, None, 1.1])
plt.show()
Total running time of the script: (1 minutes 38.987 seconds)