Run RANSAC

In this example we show how to run the RANSAC of pyprep.

# Authors: Yorguin Mantilla <yjmantilla@gmail.com>
#
# License: MIT
# Based On: use_noisy_module.py

First we import what we need for this example.

import numpy as np
import mne
from scipy import signal as signal
from time import perf_counter

from pyprep.find_noisy_channels import NoisyChannels

Now let’s make some arbitrary MNE raw object for demonstration purposes. We will think of good channels as sine waves and bad channels correlated with each other as sawtooths. The RANSAC will be biased towards sines in its prediction (they are the majority) so it will identify the sawtooths as bad. We will need to set a montage because the RANSAC needs to interpolate.

sfreq = 1000.0

# We need a montage, because RANSAC uses spherical splines for interpolation
montage = mne.channels.make_standard_montage("standard_1020")

ch_names = montage.ch_names

n_chans = len(ch_names)

info = mne.create_info(ch_names=ch_names, sfreq=sfreq, ch_types=["eeg"] * n_chans)

time = np.arange(0, 30, 1.0 / sfreq)  # 30 seconds of recording
n_bad_chans = 3

rng = np.random.default_rng(42)
bad_channels = rng.choice(np.arange(n_chans), n_bad_chans, replace=False)
bad_channels = [int(i) for i in bad_channels]
bad_ch_names = [ch_names[i] for i in bad_channels]

# The frequency components to use in the signal for good and bad channels
freq_good = 20
freq_bad = 20

# Generate the data
X = [
    signal.sawtooth(2 * np.pi * freq_bad * time)
    if i in bad_channels
    else np.sin(2 * np.pi * freq_good * time)
    for i in range(n_chans)
]
# Scale the signal amplitude and add noise.
X = 2e-5 * np.array(X) + 1e-5 * rng.random((n_chans, time.shape[0]))

raw = mne.io.RawArray(X, info)

raw.set_montage(montage, verbose=False)

Creating RawArray with float64 data, n_channels=94, n_times=30000
    Range : 0 ... 29999 =      0.000 ...    29.999 secs
Ready.

General

Measurement date	Unknown
Experimenter	Unknown
Participant	Unknown

Channels

Digitized points	97 points
Good channels	94 EEG
Bad channels	None
EOG channels	Not available
ECG channels	Not available

Data

Sampling frequency	1000.00 Hz
Highpass	0.00 Hz
Lowpass	500.00 Hz
Duration	00:00:30 (HH:MM:SS)

Assign the mne object to the NoisyChannels class. The resulting object will be the place where all following methods are performed.

nd = NoisyChannels(raw, random_state=1337)
nd2 = NoisyChannels(raw, random_state=1337)

NOTE: pick_types() is a legacy function. New code should use inst.pick(...).
Setting up high-pass filter at 1 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: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 3301 samples (3.301 s)

NOTE: pick_types() is a legacy function. New code should use inst.pick(...).
Setting up high-pass filter at 1 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: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 3301 samples (3.301 s)

Find all bad channels using channel-wise RANSAC and print a summary

start_time = perf_counter()
nd.find_bad_by_ransac(channel_wise=True)
print("--- %s seconds ---" % (perf_counter() - start_time))

# Repeat channel-wise RANSAC using a single channel at a time. This is slower
# but needs less memory.
start_time = perf_counter()
nd2.find_bad_by_ransac(channel_wise=True, max_chunk_size=1)
print("--- %s seconds ---" % (perf_counter() - start_time))

Executing RANSAC
This may take a while, so be patient...
Finding optimal chunk size : 94
Total # of chunks: 1
Current chunk:
1

RANSAC done!
--- 6.430755437000698 seconds ---
Executing RANSAC
This may take a while, so be patient...
Finding optimal chunk size : 1
Total # of chunks: 94
Current chunk:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

RANSAC done!
--- 12.828276427000674 seconds ---

Now the bad channels are saved in bads and we can continue processing our raw object. For more information, we can access attributes of the nd instance:

# Check channels that go bad together by correlation (RANSAC)
print(nd.bad_by_ransac)
assert set(bad_ch_names) == set(nd.bad_by_ransac)

# Check that the channel wise RANSAC yields identical results
print(nd2.bad_by_ransac)
assert set(bad_ch_names) == set(nd2.bad_by_ransac)

['AFz', 'P3', 'PO5']
['AFz', 'P3', 'PO5']