Note

Go to the end to download the full example code

PyLossless Algorithms

This tutorial explains the calculations that PyLossless performs at each step of the pipeline. We will use example EEG data to demonstrate the calculations.

Note

You can open this notebook in Google Colab!

Notation

Before we begin, we define some notation that will be used throughout the text:

  • We start with a 3D matrix of EEG data, \(X \in \mathbb{R}^{S_\mathcal{G} \times E_\mathcal{G} \times T}\), where \(S_\mathcal{G}\) and \(E_\mathcal{G}\) are the sets of good sensors and epochs, respectively, and \(T\), is the number of samples(i.e. time-points).

  • s, e, and t are sensor, epochs, and samples, respectively.

  • We use superscripts to denote operations across a dimension, and we use subscripts to denote indexing a dimension.

  • We refer to a single sensor \(i\) as \(X\big|_{s=i} \in \mathbb{R}^{E_\mathcal{G} \times T}\), with \(i \in S_\mathcal{G}\).

  • We refer to a single epoch \(j\) as \(X\big|_{e=j} \in \mathbb{R}^{S_\mathcal{G} \times T}\), with \(j \in E_\mathcal{G}\).

  • We denote sensor-specific thresholds for rejecting epochs as \(\tau^e_i \in \mathbb{R}^{S_\mathcal{G}}\)

  • We denote epoch-specific thresholds for rejecting sensors as \(\tau^s_j \in \mathbb{R}^{E_\mathcal{G}}\)

  • We denote quantiles as \(Q\#^{dim}\): i.e. \(Q75^s\) is the 75th quantile along the sensor dimension. The function \(Q75^s(X)\) computes the 75th quantile along the \(s\) dimension of matrix \(X\), resulting in a matrix noted \(X^{Q75^s} \in \mathbb{R}^{E \times T}\).

Throughout the text, we use capital letters for matrices and lowercase letters for scalars. For example, the data point for sensor \(i\), epoch \(j\), and time \(k\) is denoted as \(X\big|_{s=i; e=j; t=k} = x_{ijk} \in \mathbb{R}\), and \(X=\{x_{ijk}\}\).

Imports and data loading

from pathlib import Path

import numpy as np

import mne
from mne.datasets import sample

import pylossless as ll

# Load example mne data
raw = ll.datasets.load_simulated_raw()

Using default location ~/mne_data for sample...
Creating ~/mne_data

  0%|                                              | 0.00/1.65G [00:00<?, ?B/s]
  0%|                                      | 337k/1.65G [00:00<08:10, 3.37MB/s]
  0%|                                     | 1.37M/1.65G [00:00<03:41, 7.45MB/s]
  0%|                                     | 2.50M/1.65G [00:00<02:58, 9.22MB/s]
  0%|                                     | 3.69M/1.65G [00:00<02:53, 9.49MB/s]
  0%|                                     | 4.73M/1.65G [00:00<02:50, 9.65MB/s]
  0%|▏                                    | 5.86M/1.65G [00:00<02:41, 10.2MB/s]
  0%|▏                                    | 7.04M/1.65G [00:00<02:34, 10.7MB/s]
  0%|▏                                    | 8.11M/1.65G [00:01<06:03, 4.52MB/s]
  1%|▏                                    | 8.90M/1.65G [00:01<05:42, 4.80MB/s]
  1%|▏                                    | 9.63M/1.65G [00:01<05:21, 5.12MB/s]
  1%|▏                                    | 10.3M/1.65G [00:01<04:59, 5.49MB/s]
  1%|▎                                    | 11.2M/1.65G [00:01<04:22, 6.25MB/s]
  1%|▎                                    | 12.1M/1.65G [00:01<04:15, 6.43MB/s]
  1%|▎                                    | 13.4M/1.65G [00:01<03:30, 7.78MB/s]
  1%|▎                                    | 14.6M/1.65G [00:02<03:03, 8.91MB/s]
  1%|▎                                    | 15.6M/1.65G [00:02<03:07, 8.71MB/s]
  1%|▍                                    | 16.8M/1.65G [00:02<02:51, 9.52MB/s]
  1%|▍                                    | 21.3M/1.65G [00:02<01:24, 19.3MB/s]
  2%|▌                                    | 26.8M/1.65G [00:02<00:55, 29.3MB/s]
  2%|▋                                    | 32.6M/1.65G [00:02<00:43, 37.5MB/s]
  2%|▊                                    | 38.4M/1.65G [00:02<00:37, 43.5MB/s]
  3%|▉                                    | 44.2M/1.65G [00:02<00:33, 47.9MB/s]
  3%|█                                    | 50.1M/1.65G [00:02<00:31, 51.2MB/s]
  3%|█▎                                   | 56.1M/1.65G [00:02<00:29, 53.6MB/s]
  4%|█▍                                   | 61.5M/1.65G [00:03<00:29, 53.5MB/s]
  4%|█▌                                   | 67.1M/1.65G [00:03<00:29, 54.1MB/s]
  4%|█▋                                   | 72.9M/1.65G [00:03<00:28, 55.2MB/s]
  5%|█▊                                   | 78.9M/1.65G [00:03<00:27, 56.7MB/s]
  5%|█▉                                   | 84.8M/1.65G [00:03<00:27, 57.5MB/s]
  5%|██                                   | 90.7M/1.65G [00:03<00:26, 57.9MB/s]
  6%|██▏                                  | 96.5M/1.65G [00:03<00:34, 45.5MB/s]
  6%|██▎                                   | 101M/1.65G [00:04<01:14, 20.7MB/s]
  6%|██▍                                   | 105M/1.65G [00:04<01:40, 15.5MB/s]
  7%|██▍                                   | 108M/1.65G [00:05<01:55, 13.4MB/s]
  7%|██▌                                   | 110M/1.65G [00:05<02:11, 11.8MB/s]
  7%|██▌                                   | 112M/1.65G [00:05<02:22, 10.8MB/s]
  7%|██▌                                   | 114M/1.65G [00:05<02:23, 10.7MB/s]
  7%|██▋                                   | 115M/1.65G [00:05<02:31, 10.2MB/s]
  7%|██▋                                   | 116M/1.65G [00:06<02:34, 9.95MB/s]
  7%|██▋                                   | 117M/1.65G [00:06<02:35, 9.89MB/s]
  7%|██▋                                   | 118M/1.65G [00:06<02:34, 9.92MB/s]
  7%|██▋                                   | 119M/1.65G [00:06<02:42, 9.41MB/s]
  7%|██▊                                   | 120M/1.65G [00:06<02:53, 8.84MB/s]
  7%|██▊                                   | 121M/1.65G [00:06<03:08, 8.14MB/s]
  7%|██▊                                   | 122M/1.65G [00:06<03:39, 6.98MB/s]
  7%|██▊                                   | 124M/1.65G [00:07<02:54, 8.78MB/s]
  8%|██▊                                   | 125M/1.65G [00:07<02:51, 8.93MB/s]
  8%|██▉                                   | 126M/1.65G [00:07<03:02, 8.36MB/s]
  8%|██▉                                   | 127M/1.65G [00:07<02:57, 8.58MB/s]
  8%|██▉                                   | 128M/1.65G [00:07<02:52, 8.83MB/s]
  8%|██▉                                   | 129M/1.65G [00:07<02:51, 8.88MB/s]
  8%|██▉                                   | 129M/1.65G [00:07<02:55, 8.66MB/s]
  8%|██▉                                   | 130M/1.65G [00:07<03:06, 8.16MB/s]
  8%|███                                   | 132M/1.65G [00:07<02:27, 10.3MB/s]
  8%|███▏                                  | 137M/1.65G [00:08<01:08, 22.0MB/s]
  9%|███▎                                  | 143M/1.65G [00:08<00:48, 31.4MB/s]
  9%|███▍                                  | 149M/1.65G [00:08<00:38, 39.4MB/s]
  9%|███▌                                  | 153M/1.65G [00:08<00:45, 32.7MB/s]
  9%|███▌                                  | 156M/1.65G [00:08<01:18, 19.0MB/s]
 10%|███▋                                  | 159M/1.65G [00:09<01:35, 15.7MB/s]
 10%|███▋                                  | 161M/1.65G [00:09<01:56, 12.8MB/s]
 10%|███▋                                  | 163M/1.65G [00:09<02:05, 11.8MB/s]
 10%|███▊                                  | 167M/1.65G [00:09<01:29, 16.6MB/s]
 10%|███▉                                  | 172M/1.65G [00:09<01:03, 23.2MB/s]
 11%|████                                  | 176M/1.65G [00:09<01:14, 19.8MB/s]
 11%|████                                  | 178M/1.65G [00:10<01:42, 14.4MB/s]
 11%|████▏                                 | 181M/1.65G [00:10<02:00, 12.2MB/s]
 11%|████▏                                 | 182M/1.65G [00:10<02:07, 11.5MB/s]
 11%|████▏                                 | 184M/1.65G [00:11<02:27, 9.97MB/s]
 11%|████▎                                 | 185M/1.65G [00:11<02:36, 9.37MB/s]
 11%|████▎                                 | 186M/1.65G [00:11<02:37, 9.31MB/s]
 11%|████▎                                 | 187M/1.65G [00:11<02:47, 8.74MB/s]
 11%|████▎                                 | 188M/1.65G [00:11<02:57, 8.24MB/s]
 11%|████▎                                 | 189M/1.65G [00:11<03:05, 7.89MB/s]
 11%|████▎                                 | 190M/1.65G [00:11<03:00, 8.12MB/s]
 12%|████▍                                 | 191M/1.65G [00:11<02:51, 8.52MB/s]
 12%|████▍                                 | 192M/1.65G [00:12<02:56, 8.28MB/s]
 12%|████▍                                 | 193M/1.65G [00:12<02:54, 8.38MB/s]
 12%|████▍                                 | 194M/1.65G [00:12<02:48, 8.68MB/s]
 12%|████▍                                 | 195M/1.65G [00:12<02:34, 9.46MB/s]
 12%|████▌                                 | 196M/1.65G [00:12<02:38, 9.19MB/s]
 12%|████▌                                 | 197M/1.65G [00:12<02:44, 8.85MB/s]
 12%|████▌                                 | 198M/1.65G [00:12<03:10, 7.62MB/s]
 12%|████▌                                 | 199M/1.65G [00:12<03:07, 7.76MB/s]
 12%|████▌                                 | 200M/1.65G [00:12<03:15, 7.43MB/s]
 12%|████▌                                 | 200M/1.65G [00:13<03:33, 6.80MB/s]
 12%|████▋                                 | 201M/1.65G [00:13<03:17, 7.35MB/s]
 12%|████▋                                 | 202M/1.65G [00:13<03:14, 7.47MB/s]
 12%|████▋                                 | 203M/1.65G [00:13<03:27, 7.00MB/s]
 12%|████▋                                 | 204M/1.65G [00:13<03:12, 7.51MB/s]
 12%|████▋                                 | 205M/1.65G [00:13<03:19, 7.26MB/s]
 12%|████▋                                 | 206M/1.65G [00:13<03:21, 7.18MB/s]
 12%|████▋                                 | 207M/1.65G [00:13<03:16, 7.35MB/s]
 13%|████▊                                 | 208M/1.65G [00:14<03:05, 7.80MB/s]
 13%|████▊                                 | 211M/1.65G [00:14<01:46, 13.5MB/s]
 13%|████▉                                 | 215M/1.65G [00:14<01:02, 23.0MB/s]
 13%|█████                                 | 221M/1.65G [00:14<00:45, 31.8MB/s]
 14%|█████▏                                | 226M/1.65G [00:14<00:37, 38.3MB/s]
 14%|█████▎                                | 232M/1.65G [00:14<00:32, 44.0MB/s]
 14%|█████▍                                | 238M/1.65G [00:14<00:29, 48.3MB/s]
 15%|█████▌                                | 244M/1.65G [00:14<00:27, 51.4MB/s]
 15%|█████▋                                | 250M/1.65G [00:14<00:26, 53.7MB/s]
 15%|█████▊                                | 255M/1.65G [00:14<00:25, 55.4MB/s]
 16%|██████                                | 261M/1.65G [00:15<00:24, 56.6MB/s]
 16%|██████▏                               | 267M/1.65G [00:15<00:24, 57.3MB/s]
 17%|██████▎                               | 273M/1.65G [00:15<00:23, 57.9MB/s]
 17%|██████▍                               | 279M/1.65G [00:15<00:23, 58.4MB/s]
 17%|██████▌                               | 285M/1.65G [00:15<00:23, 58.7MB/s]
 18%|██████▋                               | 291M/1.65G [00:15<00:23, 58.7MB/s]
 18%|██████▊                               | 297M/1.65G [00:15<00:31, 42.7MB/s]
 18%|██████▉                               | 302M/1.65G [00:16<00:58, 23.2MB/s]
 18%|███████                               | 306M/1.65G [00:16<01:10, 19.1MB/s]
 19%|███████▏                              | 311M/1.65G [00:16<00:56, 23.7MB/s]
 19%|███████▎                              | 316M/1.65G [00:16<00:45, 29.1MB/s]
 19%|███████▍                              | 322M/1.65G [00:16<00:39, 33.9MB/s]
 20%|███████▌                              | 327M/1.65G [00:17<01:06, 19.9MB/s]
 20%|███████▌                              | 330M/1.65G [00:17<01:23, 15.8MB/s]
 20%|███████▋                              | 333M/1.65G [00:18<01:34, 13.9MB/s]
 20%|███████▋                              | 335M/1.65G [00:18<01:51, 11.9MB/s]
 21%|███████▊                              | 339M/1.65G [00:18<01:21, 16.1MB/s]
 21%|███████▉                              | 344M/1.65G [00:18<01:04, 20.3MB/s]
 21%|███████▉                              | 347M/1.65G [00:18<01:08, 19.1MB/s]
 21%|████████                              | 350M/1.65G [00:19<01:27, 15.0MB/s]
 21%|████████▏                             | 354M/1.65G [00:19<01:04, 20.0MB/s]
 22%|████████▎                             | 360M/1.65G [00:19<00:48, 26.6MB/s]
 22%|████████▍                             | 366M/1.65G [00:19<00:38, 33.3MB/s]
 22%|████████▌                             | 372M/1.65G [00:19<00:32, 39.3MB/s]
 23%|████████▋                             | 377M/1.65G [00:19<00:28, 44.0MB/s]
 23%|████████▊                             | 383M/1.65G [00:19<00:26, 48.1MB/s]
 24%|████████▉                             | 389M/1.65G [00:19<00:24, 50.7MB/s]
 24%|█████████                             | 395M/1.65G [00:19<00:23, 53.2MB/s]
 24%|█████████▏                            | 401M/1.65G [00:19<00:22, 55.0MB/s]
 25%|█████████▎                            | 407M/1.65G [00:20<00:24, 51.7MB/s]
 25%|█████████▍                            | 412M/1.65G [00:20<01:01, 20.1MB/s]
 25%|█████████▌                            | 416M/1.65G [00:21<01:27, 14.2MB/s]
 25%|█████████▋                            | 419M/1.65G [00:21<01:46, 11.6MB/s]
 26%|█████████▋                            | 421M/1.65G [00:22<01:57, 10.5MB/s]
 26%|█████████▋                            | 423M/1.65G [00:22<01:53, 10.8MB/s]
 26%|█████████▊                            | 429M/1.65G [00:22<01:15, 16.2MB/s]
 26%|█████████▉                            | 434M/1.65G [00:22<00:54, 22.2MB/s]
 27%|██████████                            | 440M/1.65G [00:22<00:42, 28.7MB/s]
 27%|██████████▎                           | 446M/1.65G [00:22<00:34, 34.9MB/s]
 27%|██████████▎                           | 451M/1.65G [00:22<00:44, 26.7MB/s]
 28%|██████████▍                           | 455M/1.65G [00:23<00:50, 23.8MB/s]
 28%|██████████▌                           | 460M/1.65G [00:23<00:40, 29.2MB/s]
 28%|██████████▋                           | 466M/1.65G [00:23<00:33, 34.9MB/s]
 29%|██████████▊                           | 472M/1.65G [00:23<00:29, 40.4MB/s]
 29%|██████████▉                           | 478M/1.65G [00:23<00:26, 44.9MB/s]
 29%|███████████                           | 484M/1.65G [00:23<00:24, 48.5MB/s]
 30%|███████████▏                          | 489M/1.65G [00:23<00:32, 36.1MB/s]
 30%|███████████▎                          | 493M/1.65G [00:24<00:55, 20.8MB/s]
 30%|███████████▍                          | 497M/1.65G [00:24<01:11, 16.2MB/s]
 30%|███████████▍                          | 499M/1.65G [00:24<01:07, 17.1MB/s]
 30%|███████████▌                          | 504M/1.65G [00:24<00:53, 21.5MB/s]
 31%|███████████▋                          | 509M/1.65G [00:25<00:41, 27.3MB/s]
 31%|███████████▊                          | 513M/1.65G [00:25<01:08, 16.6MB/s]
 31%|███████████▊                          | 516M/1.65G [00:25<01:18, 14.4MB/s]
 31%|███████████▉                          | 519M/1.65G [00:26<01:26, 13.1MB/s]
 31%|███████████▉                          | 520M/1.65G [00:26<01:33, 12.1MB/s]
 32%|████████████                          | 522M/1.65G [00:26<01:38, 11.5MB/s]
 32%|████████████                          | 523M/1.65G [00:26<01:40, 11.2MB/s]
 32%|████████████                          | 525M/1.65G [00:26<01:45, 10.7MB/s]
 32%|████████████                          | 526M/1.65G [00:26<01:54, 9.84MB/s]
 32%|████████████▏                         | 530M/1.65G [00:26<01:15, 15.0MB/s]
 32%|████████████▎                         | 535M/1.65G [00:27<00:47, 23.7MB/s]
 33%|████████████▍                         | 541M/1.65G [00:27<00:35, 31.7MB/s]
 33%|████████████▌                         | 547M/1.65G [00:27<00:28, 38.5MB/s]
 33%|████████████▋                         | 552M/1.65G [00:27<00:25, 43.9MB/s]
 34%|████████████▊                         | 558M/1.65G [00:27<00:22, 48.3MB/s]
 34%|████████████▉                         | 564M/1.65G [00:27<00:21, 51.2MB/s]
 34%|█████████████                         | 570M/1.65G [00:27<00:20, 53.3MB/s]
 35%|█████████████▏                        | 576M/1.65G [00:27<00:19, 55.0MB/s]
 35%|█████████████▍                        | 582M/1.65G [00:27<00:18, 56.4MB/s]
 36%|█████████████▌                        | 588M/1.65G [00:27<00:18, 57.3MB/s]
 36%|█████████████▋                        | 594M/1.65G [00:28<00:18, 58.1MB/s]
 36%|█████████████▊                        | 600M/1.65G [00:28<00:17, 58.7MB/s]
 37%|█████████████▉                        | 606M/1.65G [00:28<00:17, 59.2MB/s]
 37%|██████████████                        | 612M/1.65G [00:28<00:26, 39.9MB/s]
 37%|██████████████▏                       | 617M/1.65G [00:29<00:45, 22.6MB/s]
 38%|██████████████▎                       | 620M/1.65G [00:29<00:59, 17.5MB/s]
 38%|██████████████▎                       | 623M/1.65G [00:29<00:59, 17.4MB/s]
 38%|██████████████▍                       | 629M/1.65G [00:29<00:45, 22.5MB/s]
 38%|██████████████▌                       | 634M/1.65G [00:29<00:36, 28.0MB/s]
 39%|██████████████▋                       | 640M/1.65G [00:29<00:29, 34.1MB/s]
 39%|██████████████▊                       | 646M/1.65G [00:29<00:25, 39.7MB/s]
 39%|██████████████▉                       | 652M/1.65G [00:30<00:22, 44.6MB/s]
 40%|███████████████                       | 657M/1.65G [00:30<00:34, 29.2MB/s]
 40%|███████████████▏                      | 661M/1.65G [00:30<00:53, 18.5MB/s]
 40%|███████████████▎                      | 665M/1.65G [00:31<01:09, 14.2MB/s]
 40%|███████████████▎                      | 667M/1.65G [00:31<01:13, 13.5MB/s]
 41%|███████████████▍                      | 672M/1.65G [00:31<00:53, 18.5MB/s]
 41%|███████████████▌                      | 678M/1.65G [00:31<00:40, 24.1MB/s]
 41%|███████████████▋                      | 683M/1.65G [00:31<00:32, 29.9MB/s]
 42%|███████████████▊                      | 688M/1.65G [00:32<00:38, 25.3MB/s]
 42%|███████████████▉                      | 691M/1.65G [00:32<00:42, 22.5MB/s]
 42%|███████████████▉                      | 695M/1.65G [00:32<00:38, 25.1MB/s]
 42%|████████████████                      | 700M/1.65G [00:32<00:30, 31.2MB/s]
 43%|████████████████▏                     | 706M/1.65G [00:32<00:25, 37.2MB/s]
 43%|████████████████▎                     | 712M/1.65G [00:32<00:22, 42.2MB/s]
 43%|████████████████▍                     | 718M/1.65G [00:32<00:20, 46.2MB/s]
 44%|████████████████▋                     | 723M/1.65G [00:32<00:18, 49.6MB/s]
 44%|████████████████▊                     | 729M/1.65G [00:32<00:17, 52.1MB/s]
 44%|████████████████▉                     | 735M/1.65G [00:33<00:17, 53.8MB/s]
 45%|█████████████████                     | 741M/1.65G [00:33<00:16, 55.1MB/s]
 45%|█████████████████▏                    | 747M/1.65G [00:33<00:16, 56.3MB/s]
 46%|█████████████████▎                    | 753M/1.65G [00:33<00:15, 57.0MB/s]
 46%|█████████████████▍                    | 759M/1.65G [00:33<00:15, 57.6MB/s]
 46%|█████████████████▌                    | 764M/1.65G [00:33<00:15, 57.7MB/s]
 47%|█████████████████▋                    | 770M/1.65G [00:34<00:31, 28.2MB/s]
 47%|█████████████████▊                    | 775M/1.65G [00:34<00:32, 26.7MB/s]
 47%|█████████████████▉                    | 780M/1.65G [00:34<00:28, 30.9MB/s]
 48%|██████████████████                    | 785M/1.65G [00:34<00:24, 36.0MB/s]
 48%|██████████████████▏                   | 790M/1.65G [00:34<00:30, 28.7MB/s]
 48%|██████████████████▎                   | 794M/1.65G [00:35<00:42, 20.3MB/s]
 48%|██████████████████▎                   | 797M/1.65G [00:35<00:55, 15.4MB/s]
 48%|██████████████████▎                   | 799M/1.65G [00:35<01:00, 14.2MB/s]
 48%|██████████████████▍                   | 801M/1.65G [00:35<01:03, 13.4MB/s]
 49%|██████████████████▍                   | 803M/1.65G [00:35<01:06, 12.8MB/s]
 49%|██████████████████▍                   | 804M/1.65G [00:36<01:10, 12.0MB/s]
 49%|██████████████████▌                   | 806M/1.65G [00:36<01:17, 11.0MB/s]
 49%|██████████████████▌                   | 807M/1.65G [00:36<01:21, 10.3MB/s]
 49%|██████████████████▌                   | 808M/1.65G [00:36<01:21, 10.3MB/s]
 49%|██████████████████▋                   | 812M/1.65G [00:36<00:50, 16.6MB/s]
 49%|██████████████████▊                   | 816M/1.65G [00:36<00:37, 22.4MB/s]
 50%|██████████████████▊                   | 820M/1.65G [00:36<00:37, 22.3MB/s]
 50%|██████████████████▉                   | 822M/1.65G [00:37<00:49, 16.7MB/s]
 50%|██████████████████▉                   | 824M/1.65G [00:37<00:58, 14.2MB/s]
 50%|██████████████████▉                   | 826M/1.65G [00:37<01:04, 12.9MB/s]
 50%|███████████████████                   | 827M/1.65G [00:37<01:09, 11.8MB/s]
 50%|███████████████████                   | 829M/1.65G [00:37<01:13, 11.3MB/s]
 50%|███████████████████                   | 830M/1.65G [00:37<01:14, 11.1MB/s]
 50%|███████████████████                   | 831M/1.65G [00:38<01:17, 10.6MB/s]
 50%|███████████████████▏                  | 832M/1.65G [00:38<01:18, 10.5MB/s]
 50%|███████████████████▏                  | 833M/1.65G [00:38<01:21, 10.1MB/s]
 50%|███████████████████▏                  | 834M/1.65G [00:38<01:24, 9.74MB/s]
 51%|███████████████████▏                  | 835M/1.65G [00:38<01:24, 9.69MB/s]
 51%|███████████████████▏                  | 836M/1.65G [00:38<01:25, 9.59MB/s]
 51%|███████████████████▎                  | 837M/1.65G [00:38<01:22, 9.83MB/s]
 51%|███████████████████▎                  | 838M/1.65G [00:38<01:21, 9.97MB/s]
 51%|███████████████████▎                  | 839M/1.65G [00:38<01:23, 9.78MB/s]
 51%|███████████████████▎                  | 840M/1.65G [00:39<01:26, 9.41MB/s]
 51%|███████████████████▍                  | 845M/1.65G [00:39<00:39, 20.2MB/s]
 51%|███████████████████▌                  | 851M/1.65G [00:39<00:26, 30.1MB/s]
 52%|███████████████████▋                  | 857M/1.65G [00:39<00:20, 38.4MB/s]
 52%|███████████████████▊                  | 861M/1.65G [00:39<00:22, 34.7MB/s]
 52%|███████████████████▊                  | 864M/1.65G [00:39<00:40, 19.6MB/s]
 52%|███████████████████▉                  | 867M/1.65G [00:40<00:48, 16.4MB/s]
 53%|███████████████████▉                  | 869M/1.65G [00:40<00:51, 15.1MB/s]
 53%|████████████████████                  | 871M/1.65G [00:40<00:57, 13.6MB/s]
 53%|████████████████████                  | 873M/1.65G [00:40<01:02, 12.4MB/s]
 53%|████████████████████                  | 874M/1.65G [00:40<01:10, 11.1MB/s]
 53%|████████████████████▏                 | 876M/1.65G [00:41<01:14, 10.5MB/s]
 53%|████████████████████▏                 | 877M/1.65G [00:41<01:31, 8.45MB/s]
 53%|████████████████████▏                 | 878M/1.65G [00:41<01:29, 8.61MB/s]
 53%|████████████████████▏                 | 879M/1.65G [00:41<01:28, 8.71MB/s]
 53%|████████████████████▏                 | 880M/1.65G [00:41<01:28, 8.70MB/s]
 53%|████████████████████▎                 | 881M/1.65G [00:41<01:26, 8.97MB/s]
 53%|████████████████████▎                 | 882M/1.65G [00:41<01:21, 9.44MB/s]
 54%|████████████████████▎                 | 884M/1.65G [00:41<00:58, 13.1MB/s]
 54%|████████████████████▍                 | 890M/1.65G [00:42<00:31, 24.1MB/s]
 54%|████████████████████▌                 | 896M/1.65G [00:42<00:22, 33.7MB/s]
 55%|████████████████████▋                 | 902M/1.65G [00:42<00:18, 41.1MB/s]
 55%|████████████████████▊                 | 906M/1.65G [00:42<00:22, 32.5MB/s]
 55%|████████████████████▉                 | 910M/1.65G [00:42<00:31, 23.7MB/s]
 55%|█████████████████████                 | 915M/1.65G [00:42<00:25, 29.5MB/s]
 56%|█████████████████████▏                | 921M/1.65G [00:42<00:20, 36.0MB/s]
 56%|█████████████████████▎                | 927M/1.65G [00:43<00:17, 41.7MB/s]
 56%|█████████████████████▍                | 933M/1.65G [00:43<00:15, 46.5MB/s]
 57%|█████████████████████▌                | 939M/1.65G [00:43<00:14, 49.9MB/s]
 57%|█████████████████████▋                | 944M/1.65G [00:43<00:13, 52.1MB/s]
 58%|█████████████████████▊                | 951M/1.65G [00:43<00:12, 54.5MB/s]
 58%|█████████████████████▉                | 956M/1.65G [00:43<00:15, 43.7MB/s]
 58%|██████████████████████                | 961M/1.65G [00:44<00:29, 23.3MB/s]
 58%|██████████████████████▏               | 965M/1.65G [00:44<00:38, 18.0MB/s]
 59%|██████████████████████▎               | 968M/1.65G [00:44<00:44, 15.4MB/s]
 59%|██████████████████████▎               | 970M/1.65G [00:44<00:48, 14.2MB/s]
 59%|██████████████████████▎               | 972M/1.65G [00:45<00:51, 13.3MB/s]
 59%|██████████████████████▍               | 974M/1.65G [00:45<00:54, 12.6MB/s]
 59%|██████████████████████▍               | 975M/1.65G [00:45<00:56, 12.1MB/s]
 59%|██████████████████████▍               | 977M/1.65G [00:45<00:55, 12.3MB/s]
 59%|██████████████████████▌               | 981M/1.65G [00:45<00:36, 18.3MB/s]
 60%|██████████████████████▋               | 985M/1.65G [00:45<00:29, 22.7MB/s]
 60%|██████████████████████▋               | 988M/1.65G [00:45<00:25, 25.9MB/s]
 60%|██████████████████████▊               | 992M/1.65G [00:45<00:21, 30.2MB/s]
 60%|██████████████████████▉               | 998M/1.65G [00:46<00:17, 36.7MB/s]
 61%|██████████████████████▍              | 1.00G/1.65G [00:46<00:15, 41.2MB/s]
 61%|██████████████████████▌              | 1.01G/1.65G [00:46<00:14, 45.1MB/s]
 61%|██████████████████████▋              | 1.01G/1.65G [00:46<00:14, 45.3MB/s]
 62%|██████████████████████▊              | 1.02G/1.65G [00:46<00:13, 45.4MB/s]
 62%|██████████████████████▉              | 1.02G/1.65G [00:46<00:17, 36.7MB/s]
 62%|██████████████████████▉              | 1.03G/1.65G [00:47<00:29, 21.6MB/s]
 62%|███████████████████████              | 1.03G/1.65G [00:47<00:36, 16.9MB/s]
 62%|███████████████████████              | 1.03G/1.65G [00:47<00:41, 15.1MB/s]
 63%|███████████████████████▏             | 1.03G/1.65G [00:47<00:44, 13.9MB/s]
 63%|███████████████████████▏             | 1.04G/1.65G [00:47<00:47, 13.1MB/s]
 63%|███████████████████████▏             | 1.04G/1.65G [00:48<00:50, 12.3MB/s]
 63%|███████████████████████▎             | 1.04G/1.65G [00:48<00:50, 12.1MB/s]
 63%|███████████████████████▎             | 1.04G/1.65G [00:48<00:51, 11.9MB/s]
 63%|███████████████████████▎             | 1.04G/1.65G [00:48<00:34, 17.4MB/s]
 63%|███████████████████████▍             | 1.05G/1.65G [00:48<00:23, 25.8MB/s]
 64%|███████████████████████▌             | 1.06G/1.65G [00:48<00:17, 34.2MB/s]
 64%|███████████████████████▊             | 1.06G/1.65G [00:48<00:14, 40.8MB/s]
 65%|███████████████████████▉             | 1.07G/1.65G [00:48<00:12, 45.8MB/s]
 65%|████████████████████████             | 1.07G/1.65G [00:48<00:11, 49.7MB/s]
 65%|████████████████████████▏            | 1.08G/1.65G [00:49<00:10, 52.6MB/s]
 66%|████████████████████████▎            | 1.08G/1.65G [00:49<00:10, 54.7MB/s]
 66%|████████████████████████▍            | 1.09G/1.65G [00:49<00:10, 56.1MB/s]
 66%|████████████████████████▌            | 1.10G/1.65G [00:49<00:09, 56.9MB/s]
 67%|████████████████████████▋            | 1.10G/1.65G [00:49<00:09, 57.8MB/s]
 67%|████████████████████████▊            | 1.11G/1.65G [00:49<00:09, 58.3MB/s]
 67%|████████████████████████▉            | 1.11G/1.65G [00:49<00:09, 58.5MB/s]
 68%|█████████████████████████            | 1.12G/1.65G [00:49<00:09, 58.8MB/s]
 68%|█████████████████████████▏           | 1.13G/1.65G [00:49<00:08, 59.0MB/s]
 69%|█████████████████████████▎           | 1.13G/1.65G [00:50<00:20, 25.8MB/s]
 69%|█████████████████████████▍           | 1.14G/1.65G [00:51<00:34, 15.1MB/s]
 69%|█████████████████████████▌           | 1.14G/1.65G [00:51<00:39, 13.0MB/s]
 69%|█████████████████████████▌           | 1.14G/1.65G [00:51<00:41, 12.3MB/s]
 69%|█████████████████████████▋           | 1.15G/1.65G [00:51<00:36, 13.8MB/s]
 70%|█████████████████████████▊           | 1.15G/1.65G [00:51<00:25, 19.4MB/s]
 70%|█████████████████████████▉           | 1.16G/1.65G [00:52<00:19, 25.8MB/s]
 70%|██████████████████████████           | 1.16G/1.65G [00:52<00:15, 32.0MB/s]
 71%|██████████████████████████▏          | 1.17G/1.65G [00:52<00:12, 37.8MB/s]
 71%|██████████████████████████▎          | 1.17G/1.65G [00:52<00:11, 42.8MB/s]
 71%|██████████████████████████▍          | 1.18G/1.65G [00:52<00:10, 47.1MB/s]
 72%|██████████████████████████▌          | 1.19G/1.65G [00:52<00:09, 50.2MB/s]
 72%|██████████████████████████▋          | 1.19G/1.65G [00:52<00:08, 52.9MB/s]
 72%|██████████████████████████▊          | 1.20G/1.65G [00:52<00:08, 54.7MB/s]
 73%|██████████████████████████▉          | 1.20G/1.65G [00:52<00:07, 56.2MB/s]
 73%|███████████████████████████          | 1.21G/1.65G [00:52<00:07, 57.2MB/s]
 74%|███████████████████████████▏         | 1.22G/1.65G [00:53<00:07, 57.9MB/s]
 74%|███████████████████████████▎         | 1.22G/1.65G [00:53<00:07, 58.4MB/s]
 74%|███████████████████████████▍         | 1.23G/1.65G [00:53<00:07, 58.8MB/s]
 75%|███████████████████████████▌         | 1.23G/1.65G [00:53<00:14, 28.2MB/s]
 75%|███████████████████████████▋         | 1.24G/1.65G [00:54<00:21, 19.0MB/s]
 75%|███████████████████████████▊         | 1.24G/1.65G [00:54<00:24, 16.6MB/s]
 75%|███████████████████████████▉         | 1.25G/1.65G [00:54<00:20, 20.3MB/s]
 76%|████████████████████████████         | 1.25G/1.65G [00:54<00:15, 25.5MB/s]
 76%|████████████████████████████▏        | 1.26G/1.65G [00:54<00:12, 31.2MB/s]
 76%|████████████████████████████▎        | 1.26G/1.65G [00:55<00:15, 25.8MB/s]
 77%|████████████████████████████▎        | 1.27G/1.65G [00:55<00:20, 18.7MB/s]
 77%|████████████████████████████▍        | 1.27G/1.65G [00:55<00:24, 15.7MB/s]
 77%|████████████████████████████▍        | 1.27G/1.65G [00:55<00:26, 14.3MB/s]
 77%|████████████████████████████▍        | 1.27G/1.65G [00:56<00:29, 12.9MB/s]
 77%|████████████████████████████▌        | 1.27G/1.65G [00:56<00:30, 12.3MB/s]
 77%|████████████████████████████▌        | 1.28G/1.65G [00:56<00:31, 12.1MB/s]
 77%|████████████████████████████▌        | 1.28G/1.65G [00:56<00:25, 14.5MB/s]
 78%|████████████████████████████▋        | 1.28G/1.65G [00:56<00:16, 22.8MB/s]
 78%|████████████████████████████▊        | 1.29G/1.65G [00:56<00:11, 30.7MB/s]
 78%|████████████████████████████▉        | 1.30G/1.65G [00:56<00:09, 37.6MB/s]
 79%|█████████████████████████████▏       | 1.30G/1.65G [00:56<00:08, 43.2MB/s]
 79%|█████████████████████████████▎       | 1.31G/1.65G [00:57<00:07, 47.4MB/s]
 79%|█████████████████████████████▍       | 1.31G/1.65G [00:57<00:06, 50.7MB/s]
 80%|█████████████████████████████▌       | 1.32G/1.65G [00:57<00:06, 53.2MB/s]
 80%|█████████████████████████████▋       | 1.32G/1.65G [00:57<00:07, 44.1MB/s]
 80%|█████████████████████████████▊       | 1.33G/1.65G [00:57<00:14, 22.6MB/s]
 81%|█████████████████████████████▊       | 1.33G/1.65G [00:58<00:17, 18.2MB/s]
 81%|█████████████████████████████▉       | 1.34G/1.65G [00:58<00:20, 15.4MB/s]
 81%|█████████████████████████████▉       | 1.34G/1.65G [00:58<00:23, 13.5MB/s]
 81%|█████████████████████████████▉       | 1.34G/1.65G [00:58<00:24, 12.6MB/s]
 81%|██████████████████████████████       | 1.34G/1.65G [00:59<00:27, 11.4MB/s]
 81%|██████████████████████████████       | 1.34G/1.65G [00:59<00:29, 10.6MB/s]
 81%|██████████████████████████████       | 1.34G/1.65G [00:59<00:31, 9.77MB/s]
 81%|██████████████████████████████       | 1.35G/1.65G [00:59<00:32, 9.47MB/s]
 81%|██████████████████████████████▏      | 1.35G/1.65G [00:59<00:32, 9.40MB/s]
 82%|██████████████████████████████▏      | 1.35G/1.65G [00:59<00:33, 9.10MB/s]
 82%|██████████████████████████████▏      | 1.35G/1.65G [00:59<00:33, 9.16MB/s]
 82%|██████████████████████████████▏      | 1.35G/1.65G [01:00<00:32, 9.37MB/s]
 82%|██████████████████████████████▏      | 1.35G/1.65G [01:00<00:31, 9.46MB/s]
 82%|██████████████████████████████▎      | 1.35G/1.65G [01:00<00:33, 9.11MB/s]
 82%|██████████████████████████████▎      | 1.35G/1.65G [01:00<00:31, 9.41MB/s]
 82%|██████████████████████████████▎      | 1.35G/1.65G [01:00<00:31, 9.38MB/s]
 82%|██████████████████████████████▎      | 1.35G/1.65G [01:00<00:29, 10.1MB/s]
 82%|██████████████████████████████▎      | 1.36G/1.65G [01:00<00:29, 10.2MB/s]
 82%|██████████████████████████████▎      | 1.36G/1.65G [01:00<00:28, 10.2MB/s]
 82%|██████████████████████████████▍      | 1.36G/1.65G [01:00<00:15, 19.0MB/s]
 83%|██████████████████████████████▌      | 1.37G/1.65G [01:01<00:10, 27.9MB/s]
 83%|██████████████████████████████▋      | 1.37G/1.65G [01:01<00:10, 26.4MB/s]
 83%|██████████████████████████████▋      | 1.37G/1.65G [01:01<00:17, 16.1MB/s]
 83%|██████████████████████████████▊      | 1.37G/1.65G [01:01<00:19, 14.2MB/s]
 83%|██████████████████████████████▊      | 1.38G/1.65G [01:01<00:21, 13.1MB/s]
 83%|██████████████████████████████▊      | 1.38G/1.65G [01:02<00:21, 12.9MB/s]
 84%|██████████████████████████████▉      | 1.38G/1.65G [01:02<00:15, 17.2MB/s]
 84%|███████████████████████████████      | 1.39G/1.65G [01:02<00:10, 25.9MB/s]
 84%|███████████████████████████████▏     | 1.39G/1.65G [01:02<00:07, 34.0MB/s]
 85%|███████████████████████████████▎     | 1.40G/1.65G [01:02<00:06, 39.8MB/s]
 85%|███████████████████████████████▍     | 1.40G/1.65G [01:02<00:05, 45.1MB/s]
 85%|███████████████████████████████▌     | 1.41G/1.65G [01:02<00:04, 48.9MB/s]
 86%|███████████████████████████████▋     | 1.42G/1.65G [01:02<00:04, 51.9MB/s]
 86%|███████████████████████████████▊     | 1.42G/1.65G [01:02<00:04, 54.2MB/s]
 86%|███████████████████████████████▉     | 1.43G/1.65G [01:02<00:04, 55.7MB/s]
 87%|████████████████████████████████     | 1.43G/1.65G [01:03<00:03, 56.1MB/s]
 87%|████████████████████████████████▏    | 1.44G/1.65G [01:03<00:03, 56.5MB/s]
 87%|████████████████████████████████▎    | 1.44G/1.65G [01:03<00:03, 57.5MB/s]
 88%|████████████████████████████████▍    | 1.45G/1.65G [01:03<00:03, 58.3MB/s]
 88%|████████████████████████████████▌    | 1.46G/1.65G [01:03<00:03, 57.8MB/s]
 88%|████████████████████████████████▋    | 1.46G/1.65G [01:03<00:03, 57.7MB/s]
 89%|████████████████████████████████▊    | 1.47G/1.65G [01:04<00:07, 25.3MB/s]
 89%|████████████████████████████████▉    | 1.47G/1.65G [01:04<00:06, 29.0MB/s]
 89%|█████████████████████████████████    | 1.48G/1.65G [01:04<00:05, 34.0MB/s]
 90%|█████████████████████████████████▏   | 1.48G/1.65G [01:04<00:04, 39.2MB/s]
 90%|█████████████████████████████████▎   | 1.49G/1.65G [01:04<00:04, 34.2MB/s]
 90%|█████████████████████████████████▍   | 1.49G/1.65G [01:04<00:04, 36.8MB/s]
 91%|█████████████████████████████████▌   | 1.50G/1.65G [01:04<00:03, 39.6MB/s]
 91%|█████████████████████████████████▋   | 1.50G/1.65G [01:04<00:03, 42.1MB/s]
 91%|█████████████████████████████████▊   | 1.51G/1.65G [01:04<00:03, 45.2MB/s]
 92%|█████████████████████████████████▉   | 1.52G/1.65G [01:05<00:02, 48.8MB/s]
 92%|██████████████████████████████████   | 1.52G/1.65G [01:05<00:02, 51.5MB/s]
 92%|██████████████████████████████████▏  | 1.53G/1.65G [01:05<00:02, 53.8MB/s]
 93%|██████████████████████████████████▎  | 1.53G/1.65G [01:05<00:02, 55.7MB/s]
 93%|██████████████████████████████████▍  | 1.54G/1.65G [01:05<00:01, 57.0MB/s]
 94%|██████████████████████████████████▌  | 1.55G/1.65G [01:05<00:01, 58.0MB/s]
 94%|██████████████████████████████████▋  | 1.55G/1.65G [01:05<00:01, 57.6MB/s]
 94%|██████████████████████████████████▊  | 1.56G/1.65G [01:05<00:01, 57.6MB/s]
 95%|██████████████████████████████████▉  | 1.56G/1.65G [01:05<00:01, 58.3MB/s]
 95%|███████████████████████████████████▏ | 1.57G/1.65G [01:05<00:01, 58.9MB/s]
 95%|███████████████████████████████████▎ | 1.58G/1.65G [01:06<00:02, 34.8MB/s]
 96%|███████████████████████████████████▎ | 1.58G/1.65G [01:06<00:03, 18.8MB/s]
 96%|███████████████████████████████████▍ | 1.58G/1.65G [01:07<00:04, 14.3MB/s]
 96%|███████████████████████████████████▌ | 1.59G/1.65G [01:07<00:04, 14.6MB/s]
 96%|███████████████████████████████████▌ | 1.59G/1.65G [01:07<00:03, 19.6MB/s]
 97%|███████████████████████████████████▊ | 1.60G/1.65G [01:07<00:02, 25.6MB/s]
 97%|███████████████████████████████████▉ | 1.60G/1.65G [01:07<00:01, 31.6MB/s]
 97%|████████████████████████████████████ | 1.61G/1.65G [01:07<00:01, 37.3MB/s]
 98%|████████████████████████████████████▏| 1.61G/1.65G [01:07<00:00, 41.6MB/s]
 98%|████████████████████████████████████▎| 1.62G/1.65G [01:08<00:00, 45.0MB/s]
 98%|████████████████████████████████████▍| 1.63G/1.65G [01:08<00:00, 44.6MB/s]
 99%|████████████████████████████████████▌| 1.63G/1.65G [01:08<00:00, 22.8MB/s]
 99%|████████████████████████████████████▌| 1.63G/1.65G [01:09<00:01, 17.7MB/s]
 99%|████████████████████████████████████▋| 1.64G/1.65G [01:09<00:01, 15.5MB/s]
 99%|████████████████████████████████████▋| 1.64G/1.65G [01:09<00:00, 14.3MB/s]
 99%|████████████████████████████████████▊| 1.64G/1.65G [01:09<00:00, 12.8MB/s]
 99%|████████████████████████████████████▊| 1.64G/1.65G [01:09<00:00, 12.4MB/s]
100%|████████████████████████████████████▊| 1.64G/1.65G [01:10<00:00, 11.4MB/s]
100%|████████████████████████████████████▊| 1.65G/1.65G [01:10<00:00, 11.3MB/s]
100%|████████████████████████████████████▉| 1.65G/1.65G [01:10<00:00, 11.2MB/s]
100%|████████████████████████████████████▉| 1.65G/1.65G [01:10<00:00, 10.3MB/s]
100%|████████████████████████████████████▉| 1.65G/1.65G [01:10<00:00, 9.50MB/s]
100%|████████████████████████████████████▉| 1.65G/1.65G [01:10<00:00, 9.45MB/s]
100%|████████████████████████████████████▉| 1.65G/1.65G [01:10<00:00, 8.86MB/s]
100%|████████████████████████████████████▉| 1.65G/1.65G [01:11<00:00, 8.35MB/s]
  0%|                                              | 0.00/1.65G [00:00<?, ?B/s]
100%|█████████████████████████████████████| 1.65G/1.65G [00:00<00:00, 5.51TB/s]
Attempting to create new mne-python configuration file:
/home/docs/.mne/mne-python.json
Download complete in 01m34s (1576.2 MB)
Opening raw data file /home/docs/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...
    Read a total of 3 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
    Range : 25800 ... 192599 =     42.956 ...   320.670 secs
Ready.
EEG channel type selected for re-referencing
Adding average EEG reference projection.
1 projection items deactivated
Reading forward solution from /home/docs/mne_data/MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif...
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    2 source spaces read
    Desired named matrix (kind = 3523) not available
    Read MEG forward solution (7498 sources, 306 channels, free orientations)
    Desired named matrix (kind = 3523) not available
    Read EEG forward solution (7498 sources, 60 channels, free orientations)
    Forward solutions combined: MEG, EEG
    Source spaces transformed to the forward solution coordinate frame
Setting up raw simulation: 1 position, "cos2" interpolation
Event information stored on channel:              STI 014
    Interval 0.000–2.000 s
Setting up forward solutions
Computing gain matrix for transform #1/1
    Interval 0.000–2.000 s
    Interval 0.000–2.000 s
    Interval 0.000–2.000 s
    Interval 0.000–2.000 s
    Interval 0.000–2.000 s
    Interval 0.000–2.000 s
    Interval 0.000–2.000 s
    Interval 0.000–2.000 s
    Interval 0.000–2.000 s
    10 STC iterations provided
[done]
Adding noise to 366/376 channels (366 channels in cov)
Sphere                : origin at (0.0 0.0 0.0) mm
              radius  : 0.1 mm
Source location file  : dict()
Assuming input in millimeters
Assuming input in MRI coordinates

Positions (in meters) and orientations
1 sources
ecg simulated and trace not stored
Setting up forward solutions
Computing gain matrix for transform #1/1
Sphere                : origin at (0.0 0.0 0.0) mm
              radius  : 0.1 mm
Source location file  : dict()
Assuming input in millimeters
Assuming input in MRI coordinates

Positions (in meters) and orientations
2 sources
blink simulated and trace stored on channel:      EOG 061
Setting up forward solutions
Computing gain matrix for transform #1/1
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>
Adding noise to 60/60 channels (60 channels in cov)
Adding noise to 2/2 channels (2 channels in cov)
Creating RawArray with float64 data, n_channels=60, n_times=12010
    Range : 0 ... 12009 =      0.000 ...    19.995 secs
Ready.

Load a PyLossless configuration file

Let’s load a PyLossless configuration file. This file contains the parameters that will be used for each step of the pipeline. For example, the noisy_channels section contains the parameters for the Flag Noisy Sensors step. We can modify these parameters to change the behavior of the pipeline. For example, we can change the percent of epochs that a sensor must be noisy for it to be flagged via the flag_crit parameter.

config = ll.config.Config()
config.load_default()
config["noisy_channels"]["outliers_kwargs"]["lower"] = 0.25  # lower quantile
config["noisy_channels"]["outliers_kwargs"]["upper"] = 0.75  # upper quantile
config["noisy_channels"]["flag_crit"] = 0.30  # percent of epochs that a sensor must be noisy
config.save("test_config.yaml")

Create a pipeline instance

pipeline = ll.LosslessPipeline("test_config.yaml")
pipeline.raw = raw
raw.plot()
plot 0 implementation
<MNEBrowseFigure size 800x800 with 4 Axes>

Input Data

First, we epoch the data to be used for subsequent steps. Let our 3D matrix below be defined as \(X \in \mathbb{R}^{S \times E \times T}\) where \(X\) is a matrix of real numbers and of dimension \(S\) sensors \(\times$ E\) epochs times T times.

epochs = pipeline.get_epochs()

🧹 Epoching..
Not setting metadata
19 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 19 events and 601 original time points ...
0 bad epochs dropped
🔍 Detecting channels to leave out of reference.
/home/docs/checkouts/readthedocs.org/user_builds/pylossless/checkouts/latest/pylossless/pipeline.py:61: FutureWarning: The current default of copy=False will change to copy=True in 1.7. Set the value of copy explicitly to avoid this warning
  data = epochs.get_data()  # n_epochs, n_channels, n_times
EEG channel type selected for re-referencing
Applying a custom ('EEG',) reference.

Let’s convert our epochs object into a named xarray.DataArray object.

from pylossless.pipeline import epochs_to_xr

#
epochs_xr = epochs_to_xr(epochs, kind="ch")
epochs_xr  # 277 epochs, 50 sensors, 602 samples per epoch

/home/docs/checkouts/readthedocs.org/user_builds/pylossless/checkouts/latest/pylossless/pipeline.py:61: FutureWarning: The current default of copy=False will change to copy=True in 1.7. Set the value of copy explicitly to avoid this warning
  data = epochs.get_data()  # n_epochs, n_channels, n_times

<xarray.DataArray (epoch: 19, ch: 60, time: 601)>
array([[[ 1.13705379e-06, -6.91160248e-06, -1.02693862e-05, ...,
         -1.55051339e-05, -1.96876395e-05, -1.67898363e-05],
        [ 2.76299869e-07,  6.49875684e-06,  9.75350178e-06, ...,
          1.81333284e-05,  2.36450223e-05,  2.09258026e-05],
        [ 7.62474940e-07,  6.90884467e-07, -5.90707015e-08, ...,
          2.80544764e-06,  2.99591046e-06,  1.70102299e-06],
        ...,
        [ 7.72222355e-07,  1.25535897e-06,  1.86463319e-06, ...,
         -7.47594430e-07, -2.15767204e-07, -3.92131192e-07],
        [ 2.40848183e-06,  1.83006184e-06,  2.63754167e-06, ...,
         -6.02942538e-08, -1.44093900e-06, -2.68510990e-06],
        [ 1.86663177e-06,  2.30039526e-06,  2.25657409e-06, ...,
          1.80076778e-06,  8.94905868e-07,  1.58496214e-06]],

       [[-1.67898363e-05, -5.68470102e-06, -4.77038844e-06, ...,
         -9.96482098e-06, -1.06521991e-05, -9.60757659e-06],
        [ 2.09258026e-05,  1.01833300e-05,  8.93191582e-06, ...,
          1.29011842e-05,  1.29407107e-05,  9.83793662e-06],
        [ 1.70102299e-06,  8.17735957e-07, -3.09025237e-07, ...,
         -3.40630115e-06, -2.03742827e-06, -1.44760682e-06],
...
        [-1.75260082e-06, -6.65040505e-07,  5.98434115e-07, ...,
         -3.85315653e-06, -3.48863524e-06, -1.86262811e-06],
        [ 2.59756351e-06,  2.97587549e-06,  2.22933218e-06, ...,
          7.41575231e-07,  9.58793161e-07,  8.96293208e-07],
        [ 2.20825173e-06,  2.63437479e-06,  1.34431406e-06, ...,
          4.59924450e-07,  1.46033214e-07,  1.29851490e-06]],

       [[ 5.47256633e-07, -1.33673297e-05, -8.43730917e-06, ...,
          5.12980366e-06, -5.26759727e-06, -6.84131810e-06],
        [-7.83470169e-07,  1.16691012e-05,  8.05896806e-06, ...,
         -6.88354432e-06,  1.60840737e-06,  3.62822165e-06],
        [ 1.50443446e-06,  1.22282185e-06,  3.32821461e-06, ...,
         -6.88381587e-07,  3.86112624e-07, -8.88390657e-08],
        ...,
        [ 1.02084541e-06,  1.77217630e-06,  2.62686768e-06, ...,
         -1.16207286e-06, -7.86727671e-07,  2.66763846e-08],
        [ 1.80384005e-06,  1.78232559e-06,  4.85693880e-07, ...,
         -9.90702721e-08, -1.53454913e-07,  1.43268509e-06],
        [ 2.49885075e-07, -1.00965428e-06,  1.60455798e-06, ...,
          3.54097415e-06,  3.12948165e-06,  2.87620497e-06]]])
Coordinates:
  * epoch    (epoch) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
  * ch       (ch) <U7 'EEG 001' 'EEG 002' 'EEG 003' ... 'EEG 059' 'EEG 060'
  * time     (time) float64 0.0 0.001665 0.00333 ... 0.9956 0.9973 0.999


Robust Average Reference

Robust Average Reference graphic.

Robust Average Reference. The figure shows the steps for robust average referencing. See the text below for descriptions of mathematical notation.

Before the pipeline can begin, we must average reference the data. This is because the pipeline uses data distributions to identify noisy sensors, and For EEG data that uses an online reference to a single electrode, sensors that are further from the reference will have a higher voltage variance, and the pipeline will be biased to flag these sensors as noisy. The average reference, which subtracts the average signal across sensors from each individual sensor, will ensure an even playing field. Howeer, we dont want to include noisy sensors in the average reference signal. So we will identify noisy sensors and and leave them out of the average reference signal.

sample_std = epochs_xr.std("time")
q25_ch = sample_std.quantile(0.25, dim="ch")
q50_ch = sample_std.median(dim="ch")
q75_ch = sample_std.quantile(0.75, dim="ch")
ch_dist = sample_std - q50_ch  # center the data
ch_dist /= q75_ch - q25_ch  # shape (chans, epoch)

mean_ch_dist = ch_dist.mean(dim="epoch")  # shape (chans)

# find the median and 25 and 75 percentiles
# of the mean of the channel distributions
mdn = np.median(mean_ch_dist)
deviation = np.diff(np.quantile(mean_ch_dist, [0.25, 0.75]))

leave_out = mean_ch_dist.ch[mean_ch_dist > mdn + 6 * deviation].values.tolist()
leave_out

['EEG 001', 'EEG 002', 'EEG 007']

ref_chans = [ch for ch in epochs.pick("eeg").ch_names if ch not in leave_out]
pipeline.raw.set_eeg_reference(ref_channels=ref_chans)

EEG channel type selected for re-referencing
Applying a custom ('EEG',) reference.
General
Measurement date Unknown
Experimenter Unknown
Participant Unknown
Channels
Digitized points 145 points
Good channels 60 EEG
Bad channels None
EOG channels Not available
ECG channels Not available
Data
Sampling frequency 600.61 Hz
Highpass 0.00 Hz
Lowpass 300.31 Hz
Duration 00:00:20 (HH:MM:SS)


Flag Noisy Sensors

Flag Noisy Sensors graphic.

Flag Noisy Sensors. The figure shows the steps for flagging noisy sensors. See the text below for descriptions of mathematical notation.

Since we applied a robust average reference to the raw data, we will need to re-epoch the data:

epochs = pipeline.get_epochs()
epochs_xr = epochs_to_xr(epochs, kind="ch")

# First we take standard deviation of
# :math:`X \in \mathbb{R}^{S \times E \times T}` across the samples dimension :math:`t`
# resulting in a 2D matrix :math:`X^{\sigma_{t}} \in \mathbb{R}^{S \times E}`
trim_ch_sd = epochs_xr.std("time")
trim_ch_sd

🧹 Epoching..
Not setting metadata
19 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 19 events and 601 original time points ...
0 bad epochs dropped
🔍 Detecting channels to leave out of reference.
/home/docs/checkouts/readthedocs.org/user_builds/pylossless/checkouts/latest/pylossless/pipeline.py:61: FutureWarning: The current default of copy=False will change to copy=True in 1.7. Set the value of copy explicitly to avoid this warning
  data = epochs.get_data()  # n_epochs, n_channels, n_times
EEG channel type selected for re-referencing
Applying a custom ('EEG',) reference.
/home/docs/checkouts/readthedocs.org/user_builds/pylossless/checkouts/latest/pylossless/pipeline.py:61: FutureWarning: The current default of copy=False will change to copy=True in 1.7. Set the value of copy explicitly to avoid this warning
  data = epochs.get_data()  # n_epochs, n_channels, n_times

<xarray.DataArray (epoch: 19, ch: 60)>
array([[1.33923477e-05, 1.31993382e-05, 2.35078082e-06, ...,
        2.10354957e-06, 2.26271117e-06, 2.17432334e-06],
       [1.35833591e-05, 1.36148881e-05, 1.79525497e-06, ...,
        1.96148280e-06, 1.84403073e-06, 1.70296604e-06],
       [2.56249005e-05, 2.31628291e-05, 1.83310921e-05, ...,
        1.87991905e-05, 1.89919706e-05, 1.66146307e-05],
       ...,
       [1.27689414e-05, 1.27142496e-05, 2.05817360e-06, ...,
        2.34968279e-06, 1.96616377e-06, 2.25194570e-06],
       [1.19915305e-05, 1.21115757e-05, 1.86364777e-06, ...,
        2.13783198e-06, 1.74270583e-06, 1.85174223e-06],
       [1.88604027e-05, 1.37044171e-05, 6.09928695e-06, ...,
        2.86932172e-06, 2.86073946e-06, 3.40488560e-06]])
Coordinates:
  * epoch    (epoch) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
  * ch       (ch) <U7 'EEG 001' 'EEG 002' 'EEG 003' ... 'EEG 059' 'EEG 060'


a) Take the 50th and 75th quantile across dimension sensor of \(X^{\sigma_{t}}\)

This operation results in two 1D vectors of size \(E\):

\[X^{{\sigma}_t{Q50^s}} = Q50^s(X^{\sigma_{t}}) \in \mathbb{R}^{E}\]
\[X^{{\sigma}_t{Q75^s}} = Q75^s(X^{\sigma_{t}}) \in \mathbb{R}^{E}\]
q50, q75 = trim_ch_sd.quantile([0.5, 0.75], dim="ch")
q50  # a 1D array of median standard deviation values across channels for each epoch

<xarray.DataArray (epoch: 19)>
array([1.95493588e-06, 1.83232558e-06, 1.86303335e-05, 1.87396007e-06,
       2.49031899e-06, 1.81729805e-06, 2.00626414e-06, 1.82046007e-06,
       3.21349673e-06, 2.02533760e-06, 2.13707921e-06, 1.80897263e-06,
       1.97889010e-06, 1.85854369e-06, 2.00078461e-06, 1.83742319e-06,
       2.00708147e-06, 1.82894494e-06, 2.79858811e-06])
Coordinates:
  * epoch     (epoch) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
    quantile  float64 0.5


b) Define an Upper Quantile Range as \(Q75 - Q50\)

\[UQR^s = X^{{\sigma}_T{Q75}^s} - X^{{\sigma}_T{Q50}^s}\]

This operation results in a 1D vector of size \(E\).

uqr = q75 - q50
uqr

<xarray.DataArray (epoch: 19)>
array([1.84237618e-07, 5.11518576e-08, 9.44854829e-07, 1.50125438e-07,
       9.62376628e-07, 6.88972902e-08, 1.75000051e-07, 7.36688444e-08,
       1.54999431e-06, 2.54329125e-07, 5.04755307e-07, 7.39668745e-08,
       1.80867488e-07, 7.13429008e-08, 1.86393945e-07, 6.30284767e-08,
       1.10723912e-07, 8.59682952e-08, 1.11922998e-06])
Coordinates:
  * epoch    (epoch) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18


c) Identify outlier Indices \((i, j)\)

We multiply a constant \(k\) by the \(UQR\) to define a measure for the spread of the right tail of the distribution of \(X^{\sigma_{t}}\) values and add it to the median of \(X^{\sigma_{t}}\) to obtain epoch-specific standard deviation threshold for outliers:

\[\tau^s_j = X^{{\sigma}_T{Q50}^S} + UQR^s\times k\]

That is, \(\tau^s_j\) is the epoch-specific threshold for the epoch \(j\)

k = 3
upper_threshold = q50 + q75 * k
upper_threshold  # epoch specific thresholds

<xarray.DataArray (epoch: 19)>
array([8.37245637e-06, 7.48275789e-06, 7.73558986e-05, 7.94621660e-06,
       1.28484059e-05, 7.47588408e-06, 8.55005672e-06, 7.50284682e-06,
       1.75039698e-05, 8.86433779e-06, 1.00625828e-05, 7.45779116e-06,
       8.45816287e-06, 7.64820345e-06, 8.56232029e-06, 7.53877820e-06,
       8.36049760e-06, 7.57368467e-06, 1.45520424e-05])
Coordinates:
  * epoch    (epoch) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18


Now, we compare our 2D standard deviation matrix to the threshold vector of \(\tau^e_j\):

\[X^{\sigma_{t}} \big|_{e=j} > \tau^s_j\]

resulting in the indicator matrix \(C \in \{0, 1\}^{S \times E}=\{c_{ij}\}\):

\[\begin{split}c_{ij} = \begin{cases} 0 & \text{if } x^{\sigma_{t}}_{ij} < \tau^s_j \\ 1 & \text{if } x^{\sigma_{t}}_{ij} \geq \tau^s_j \end{cases}\end{split}\]

Each element of this matrix indicates whether sensor \(i\) is an outlier at an epoch \(j\).

outlier_mask = trim_ch_sd > upper_threshold
outlier_mask

<xarray.DataArray (epoch: 19, ch: 60)>
array([[ True,  True, False, ..., False, False, False],
       [ True,  True, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [ True,  True, False, ..., False, False, False],
       [ True,  True, False, ..., False, False, False],
       [ True, False, False, ..., False, False, False]])
Coordinates:
  * epoch    (epoch) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
  * ch       (ch) <U7 'EEG 001' 'EEG 002' 'EEG 003' ... 'EEG 059' 'EEG 060'


d) Identify noisy sensors part 1

To identify outlier sensors, we average across the epoch dimension of our indicator matrix \(C\) and obtain \(C^{\mu_e} \in \mathbb{R}^{S_\mathcal{G}}\), which is a vector of fractional numbers \(c^{\mu_e}_i\) representing the percentage of epochs for which that sensor is an outlier.

percent_outliers = outlier_mask.astype(float).mean("epoch")
percent_outliers  # percent of epochs that sensor is an outlier

<xarray.DataArray (ch: 60)>
array([0.94736842, 0.84210526, 0.        , 0.        , 0.        ,
       0.        , 0.26315789, 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ])
Coordinates:
  * ch       (ch) <U7 'EEG 001' 'EEG 002' 'EEG 003' ... 'EEG 059' 'EEG 060'


e) Identify noisy sensors part 2

Next, we define a threshold \(\tau^{p}\) (\(p\) for percentile; default, .20) as a cutoff point for determining if a sensor should be marked artifactual. The sensor \(i\) is flagged as noisy if \(c^{\mu_e}_i > \tau^{p}\). That is, if the sensor is an outlier for more than \(\tau^{p}\) percent of the epochs, it is flagged as noisy.

p_threshold = config["noisy_channels"]["flag_crit"]  # 0.3, or 30%
noisy_chs = percent_outliers[percent_outliers > p_threshold].coords.to_index().values
noisy_chs

array(['EEG 001', 'EEG 002'], dtype=object)

f) Add the noisy sensors to the pipeline flags

Let’s add the noisy sensors to the pipeline flags.

pipeline.flags["ch"].add_flag_cat(kind="noisy", bad_ch_names=noisy_chs)
pipeline.raw.info["bads"].extend(pipeline.flags["ch"]["noisy"].tolist())
pipeline.flags["ch"]

Flagged channels: |
  Noisy: None
  Bridged: None
  Uncorrelated: None
  Rank: None

Flag Noisy Epochs

This step closely resembles the Flag Noisy Sensors step. For the sake of brevity we will be more concise in the documentation.

https://raw.githubusercontent.com/scott-huberty/wip_pipeline-figures/main/Flag_noisy_epochs.png

Flag Noisy Epochs. The figure shows the steps for flagging noisy epochs. See the text below for descriptions of mathematical notation.

a) Take standard deviation across the samples dimension \(t\)

Take a moment below to notice that the sensors flagged in the prior setp are not in epochs_xr below:

epochs = pipeline.get_epochs()
# Let's make our epochs array into a named Array
epochs_xr = epochs_to_xr(epochs, kind="ch")
trim_ch_sd = epochs_xr.std("time")
trim_ch_sd.coords["ch"]

🧹 Epoching..
Not setting metadata
19 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 19 events and 601 original time points ...
0 bad epochs dropped
🔍 Detecting channels to leave out of reference.
/home/docs/checkouts/readthedocs.org/user_builds/pylossless/checkouts/latest/pylossless/pipeline.py:61: FutureWarning: The current default of copy=False will change to copy=True in 1.7. Set the value of copy explicitly to avoid this warning
  data = epochs.get_data()  # n_epochs, n_channels, n_times
EEG channel type selected for re-referencing
Applying a custom ('EEG',) reference.
/home/docs/checkouts/readthedocs.org/user_builds/pylossless/checkouts/latest/pylossless/pipeline.py:61: FutureWarning: The current default of copy=False will change to copy=True in 1.7. Set the value of copy explicitly to avoid this warning
  data = epochs.get_data()  # n_epochs, n_channels, n_times

<xarray.DataArray 'ch' (ch: 58)>
array(['EEG 003', 'EEG 004', 'EEG 005', 'EEG 006', 'EEG 007', 'EEG 008',
       'EEG 009', 'EEG 010', 'EEG 011', 'EEG 012', 'EEG 013', 'EEG 014',
       'EEG 015', 'EEG 016', 'EEG 017', 'EEG 018', 'EEG 019', 'EEG 020',
       'EEG 021', 'EEG 022', 'EEG 023', 'EEG 024', 'EEG 025', 'EEG 026',
       'EEG 027', 'EEG 028', 'EEG 029', 'EEG 030', 'EEG 031', 'EEG 032',
       'EEG 033', 'EEG 034', 'EEG 035', 'EEG 036', 'EEG 037', 'EEG 038',
       'EEG 039', 'EEG 040', 'EEG 041', 'EEG 042', 'EEG 043', 'EEG 044',
       'EEG 045', 'EEG 046', 'EEG 047', 'EEG 048', 'EEG 049', 'EEG 050',
       'EEG 051', 'EEG 052', 'EEG 053', 'EEG 054', 'EEG 055', 'EEG 056',
       'EEG 057', 'EEG 058', 'EEG 059', 'EEG 060'], dtype='<U7')
Coordinates:
  * ch       (ch) <U7 'EEG 003' 'EEG 004' 'EEG 005' ... 'EEG 059' 'EEG 060'


b) Compute 50th and 75th quantile across epochs and the UQR

Like before, We Take the median and 70th quantile, but now we operate across epochs, resulting in two 1D vector’s of size n_good_sensors \(S_\mathcal{G}\)

\[X^{{\sigma}_t{Q50^e}} = Q50^e(X^{\sigma_{t}}) \in \mathbb{R}^{S_\mathcal{G}}\]
\[X^{{\sigma}_t{Q75^e}} = Q75^e(X^{\sigma_{t}}) \in \mathbb{R}^{S_\mathcal{G}}\]
\[UQR^e = (X^{{\sigma}_T{Q75}^e} - X^{{\sigma}_T{Q50}^e})\]

c) Define sensor-specific thresholds for rejecting epochs \(\tau^e_i\)

Our sensor-specifc threshold for rejecting epochs is defined by:

\[\tau^e_i = X^{{\sigma}_T{Q50}^e} + UQR^e\times k\]
q50, q75 = trim_ch_sd.quantile([0.5, 0.75], dim="epoch")
uqr_epoch = q75 - q50
uqr_epoch

<xarray.DataArray (ch: 58)>
array([1.00013545e-06, 3.45363819e-06, 3.82001346e-06, 1.76506725e-06,
       8.56916740e-06, 1.10041292e-07, 7.89424895e-07, 8.35588450e-07,
       3.59619817e-07, 2.24753021e-07, 4.16920007e-07, 5.41729070e-07,
       1.35679512e-06, 3.68594924e-07, 4.89241699e-07, 3.26336145e-07,
       8.90045179e-08, 2.74341470e-07, 5.06464178e-07, 1.59053581e-07,
       4.20889632e-08, 2.61037838e-07, 2.00681999e-07, 3.81847616e-07,
       9.94571229e-08, 1.48083330e-07, 1.22486569e-06, 1.24387913e-07,
       1.69031663e-07, 1.34891549e-07, 7.19681135e-08, 1.68176799e-07,
       1.34279990e-07, 3.95788264e-07, 2.04583376e-07, 1.10848156e-07,
       8.64738563e-08, 1.11935198e-07, 1.22813229e-07, 2.07262840e-07,
       4.30558412e-07, 2.89290954e-07, 1.48163635e-07, 1.93594851e-07,
       2.24343781e-07, 1.25047727e-07, 1.96501569e-07, 1.39680962e-07,
       2.17582054e-07, 1.19969136e-07, 1.33108394e-07, 1.33108394e-07,
       2.18810351e-07, 2.50953169e-07, 8.38153870e-08, 1.57694187e-07,
       2.88693153e-07, 2.90642956e-07])
Coordinates:
  * ch       (ch) <U7 'EEG 003' 'EEG 004' 'EEG 005' ... 'EEG 059' 'EEG 060'


k = 8
upper_threshold = q50 + uqr_epoch * k
upper_threshold

<xarray.DataArray (ch: 58)>
array([1.00592572e-05, 2.94869852e-05, 3.26392012e-05, 1.62648688e-05,
       7.05506615e-05, 2.80169357e-06, 8.22389450e-06, 8.77924053e-06,
       5.06337810e-06, 4.00799945e-06, 5.57097025e-06, 6.43864045e-06,
       1.27783481e-05, 4.79807706e-06, 6.42977147e-06, 4.53279203e-06,
       2.59074383e-06, 4.23301768e-06, 6.01481975e-06, 3.19295037e-06,
       2.23379226e-06, 3.92688771e-06, 3.98410797e-06, 4.98664216e-06,
       2.65475237e-06, 3.09882747e-06, 1.19776328e-05, 2.97287669e-06,
       3.30703112e-06, 2.99518479e-06, 2.52960465e-06, 3.19560132e-06,
       3.02933035e-06, 5.35446348e-06, 3.67814536e-06, 2.74318580e-06,
       2.53668876e-06, 2.71152057e-06, 2.85830228e-06, 3.55914984e-06,
       5.32412901e-06, 4.35545101e-06, 3.13076799e-06, 3.44582207e-06,
       3.68865538e-06, 2.83537563e-06, 3.51300214e-06, 2.91043574e-06,
       3.64020116e-06, 2.93078184e-06, 2.97544570e-06, 2.97544570e-06,
       3.63553730e-06, 3.98879565e-06, 2.75154221e-06, 3.35120085e-06,
       4.29810525e-06, 4.49946698e-06])
Coordinates:
  * ch        (ch) <U7 'EEG 003' 'EEG 004' 'EEG 005' ... 'EEG 059' 'EEG 060'
    quantile  float64 0.5


d) Identify Outlier indices

The indicator matrix is defined by:

\[\begin{split}c_{ij} = \begin{cases} 0 & \text{if } x^{\sigma_{t}}_{ij} < \tau^e_i \\ 1 & \text{if } x^{\sigma_{t}}_{ij} \geq \tau^e_i \end{cases}\end{split}\]

To identify outlier epochs, we average across the sensor dimension of our indicator matrix \(C\) and obtain \(C^{\mu_s} \in \mathbb{R}^{E_\mathcal{G}}\), which is a vector of numbers \(c^{\mu_s}_j\) representing the percentage of sensors for which that epoch is an outlier.

outlier_mask = trim_ch_sd > upper_threshold
outlier_mask

<xarray.DataArray (epoch: 19, ch: 58)>
array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [ True, False, False, ...,  True,  True,  True],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]])
Coordinates:
  * epoch     (epoch) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
  * ch        (ch) <U7 'EEG 003' 'EEG 004' 'EEG 005' ... 'EEG 059' 'EEG 060'
    quantile  float64 0.5


percent_outliers = outlier_mask.astype(float).mean("ch")
percent_outliers

<xarray.DataArray (epoch: 19)>
array([0.        , 0.        , 0.93103448, 0.        , 0.01724138,
       0.        , 0.        , 0.        , 0.13793103, 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.03448276])
Coordinates:
  * epoch     (epoch) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
    quantile  float64 0.5


e) Identify noisy epochs

Next, we define a fractional threshold \(\tau^{p}\) as a cutoff point for determining if a epoch should be marked artifactual. The epoch \(j\) is flagged as noisy if \(c^{\mu_s}_j > \tau^{p}\).

bad_epochs = percent_outliers[percent_outliers > p_threshold].coords.to_index().values
bad_epochs

array([2])

f) Add the noisy epochs to the pipeline flags

Let’s add the outlier epochs to our flags These will be added directly as mne.Annotations to the raw data.

pipeline.flags["epoch"].add_flag_cat(
    kind="noisy", bad_epoch_inds=bad_epochs, epochs=epochs
)
pipeline.raw.annotations.description

📋 LOSSLESS: 1.0 second(s) flagged as BAD_LL_noisy

array(['BAD_LL_noisy'], dtype='<U12')

pipeline.raw.plot()
plot 0 implementation
<MNEBrowseFigure size 800x800 with 4 Axes>

Filtering

After flagging noisy sensors and epochs, we filter the data. By default, The pipeline uses a 1-100Hz bandpass filter. This is because 1), ICA decompositions are more stable when low frequency drifts are removed, and 2) the ICLabel classifier is trained on data that has been filtered between 1-100Hz. A notch filter can also be optionally specified.

pipeline.config["filtering"]["notch_filter_args"]["freqs"] = [50]
pipeline.filter()

LOSSLESS: 👇 filter.
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 1 - 1e+02 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass 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)
- Upper passband edge: 100.00 Hz
- Upper transition bandwidth: 25.00 Hz (-6 dB cutoff frequency: 112.50 Hz)
- Filter length: 1983 samples (3.302 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
Filtering raw data in 1 contiguous segment
Setting up band-stop filter from 49 - 51 Hz

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 passband edge: 49.38
- Lower transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 49.12 Hz)
- Upper passband edge: 50.62 Hz
- Upper transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 50.88 Hz)
- Filter length: 3965 samples (6.602 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
Filtering raw data in 1 contiguous segment
Setting up band-stop filter from 49 - 51 Hz

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 passband edge: 49.38
- Lower transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 49.12 Hz)
- Upper passband edge: 50.62 Hz
- Upper transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 50.88 Hz)
- Filter length: 3965 samples (6.602 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
LOSSLESS: 🏁 Finished filter after 0.21 seconds.

Find Nearest Neighbours & return Maximum Correlation

Nearest Neighbors graphic.

Nearest Neighbors. The figure shows the steps for finding nearest neighbors. See the text below for descriptions of mathematical notation.

Whereas Flag Noisy Sensors and Flag Noisy Epochs operated on a 2D matrix of standard deviation values, The next few steps will operate on correlation coefficients. Here we describe the procedure for defining the 2D matrix of correlation coefficients.

from pylossless.pipeline import chan_neighbour_r

Notice that our flagged epochs are dropped.

epochs = pipeline.get_epochs()

🧹 Epoching..
Not setting metadata
19 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 19 events and 601 original time points ...
2 bad epochs dropped
🔍 Detecting channels to leave out of reference.
/home/docs/checkouts/readthedocs.org/user_builds/pylossless/checkouts/latest/pylossless/pipeline.py:61: FutureWarning: The current default of copy=False will change to copy=True in 1.7. Set the value of copy explicitly to avoid this warning
  data = epochs.get_data()  # n_epochs, n_channels, n_times
EEG channel type selected for re-referencing
Applying a custom ('EEG',) reference.

a) Calculate Correlation Coefficients between each Sensor and its neighboring eighbors

  • For each good sensor $i$ in \(S_{\mathcal{G}}\), we select its \(N\) nearest neighbors. I.e. the \(N\) sensors that are closest to it.

  • We call the sensor \(i\) the origin, and its nearest neighbors :math`hat{s_l}` with \(l \in \{1, 2, \ldots, N\}\)

  • Then, for each epoch \(j\), we calculate the correlation coefficient \(\rho^t_{(i,\hat{s_l}),j}\) between origin sensor \(i\) and each neighbor \(\hat{s_l}\) across dimension \(t\) (samples), returning a 3D matrice of correlation coefficients:

\[\mathrm{P}^t = \{\rho^t_{(i, \hat{s_l}),j}\} \in \mathbb{R}^{S_G \times E_G \times n}\]

Finally, we select the maximum correlation coefficient across the neighbor dimension \(n\):

\[\mathrm{P}^{t,{\text{max}}^n}= \max\limits_{\hat{s_l}} \rho^t_{(i, \hat{s_l}),j}\]

Returning a 2D matrix where each value at \((i, j)\) is the maximum correlation coefficient between sensor \(i\) and its \(N\) nearest neighbors, at each epoch \(j\)

data_r_ch = chan_neighbour_r(epochs, nneigbr=3, method="max")
# maximum correlation out of correlations between ch and its 3 neighbors
data_r_ch

  0%|          | 0/58 [00:00<?, ?it/s]
 10%|█         | 6/58 [00:00<00:00, 55.84it/s]
 21%|██        | 12/58 [00:00<00:00, 56.39it/s]
 31%|███       | 18/58 [00:00<00:00, 56.50it/s]
 41%|████▏     | 24/58 [00:00<00:00, 56.48it/s]
 52%|█████▏    | 30/58 [00:00<00:00, 56.37it/s]
 62%|██████▏   | 36/58 [00:00<00:00, 56.14it/s]
 72%|███████▏  | 42/58 [00:00<00:00, 56.23it/s]
 83%|████████▎ | 48/58 [00:00<00:00, 56.14it/s]
 93%|█████████▎| 54/58 [00:00<00:00, 56.20it/s]
100%|██████████| 58/58 [00:01<00:00, 56.26it/s]

<xarray.DataArray (ch: 58, epoch: 17)>
array([[0.47509132, 0.10390411, 0.92440533, 0.0673749 , 0.35361236,
        0.09878501, 0.95152367, 0.79810564, 0.81160935, 0.08853303,
        0.317296  , 0.05409002, 0.33585425, 0.12059436, 0.42590711,
        0.22011646, 0.9408646 ],
       [0.12825352, 0.29361574, 0.95814263, 0.21848927, 0.08475667,
        0.11232181, 0.983097  , 0.87469634, 0.88550819, 0.10098821,
        0.18006057, 0.06655078, 0.14338218, 0.06734977, 0.11634325,
        0.37225122, 0.97075647],
       [0.25862341, 0.30038307, 0.95814263, 0.21848927, 0.15158083,
        0.19218838, 0.983097  , 0.87469634, 0.88550819, 0.10630883,
        0.37466424, 0.08551102, 0.14815096, 0.13974137, 0.26023929,
        0.37225122, 0.97075647],
       [0.48762523, 0.1364003 , 0.93883571, 0.1565395 , 0.42390116,
        0.10945235, 0.97434588, 0.87228409, 0.87775807, 0.10253312,
        0.34017858, 0.1408269 , 0.49947826, 0.07968345, 0.42590711,
        0.27671893, 0.95981436],
       [0.42200668, 0.1364003 , 0.93883571, 0.24258576, 0.34117969,
        0.12457922, 0.97434588, 0.87228409, 0.87775807, 0.10253312,
        0.20519497, 0.1971861 , 0.37873406, 0.10518128, 0.36585942,
        0.28780298, 0.95981436],
...
       [0.29042428, 0.15838339, 0.48891537, 0.09079786, 0.31102356,
        0.18043789, 0.72488604, 0.10977672, 0.2813894 , 0.15409647,
        0.25672442, 0.15930577, 0.23447126, 0.11240078, 0.24234849,
        0.09767763, 0.56241932],
       [0.25967832, 0.17170418, 0.55231119, 0.10258568, 0.3008919 ,
        0.17661823, 0.73512567, 0.19009366, 0.39053985, 0.14610338,
        0.29190444, 0.11591035, 0.164513  , 0.21104717, 0.26288218,
        0.04975481, 0.59656811],
       [0.25967832, 0.21588847, 0.60288648, 0.05461549, 0.2594587 ,
        0.06640184, 0.73512567, 0.13003028, 0.4694078 , 0.1431741 ,
        0.36572649, 0.10992836, 0.24173394, 0.16121281, 0.41740997,
        0.12890009, 0.70291817],
       [0.2857795 , 0.15838339, 0.54993913, 0.21882875, 0.31102356,
        0.05738706, 0.72488604, 0.13026371, 0.44592835, 0.1551928 ,
        0.40007881, 0.0823491 , 0.26668939, 0.06750884, 0.26337664,
        0.12890009, 0.70291817],
       [0.30896966, 0.21588847, 0.6081601 , 0.19998163, 0.31506985,
        0.06640184, 0.75291729, 0.20817757, 0.4694078 , 0.1431741 ,
        0.35385313, 0.10992836, 0.36842147, 0.13241676, 0.41740997,
        0.13631191, 0.70638937]])
Coordinates:
  * epoch    (epoch) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
  * ch       (ch) <U7 'EEG 003' 'EEG 004' 'EEG 005' ... 'EEG 059' 'EEG 060'


This matrix \(\mathrm{P}^{t,{\text{max}}^n}\) will be used in the steps below.

Flag Bridged Sensors

Flag Bridged Sensors graphic.

Flag Bridged Sensors. The figure shows the steps for flagging bridged sensors. See the text below for descriptions of mathematical notation.

a) Calculate the 50th, 75th quantile and IQR across epochs

\[IQR^e = \mathrm{P}^{t,{\text{max}}^nQ75^e} - \mathrm{P}^{t,{\text{max}}^nQ25^e}\]

For each sensor, divide the median across epochs by the IQR across epochs. Bridged channels should have a high median correlation but a low IQR of the correlation. We call this measure the bridge-indicator.

\[\mathcal{B}_s = \frac{\mathrm{P}^{t,{\text{max}}^nQ50^e}}{IQR^e}\]

b) Define a bridging threshold

Now, take the 25th, 50th, and 75th quantile of \(\mathcal{B}_s\) across sensors, And calculate the \(IQR^s\). A channel \(i\) is bridged if

\[\mathcal{B}_i > B^{Q50^s} +k \times IQR^s\]
import scipy
from functools import partial

msr = data_r_ch.median("epoch") / data_r_ch.reduce(scipy.stats.iqr, dim="epoch")
# msr is a 1D vector of size n_sensors
config_trim = 40
config_bridge_z = 6
#
trim = config_trim
if trim >= 1:
    trim /= 100
trim /= 2  # .20 and will be used as (.20, .20)
#
trim_mean = partial(scipy.stats.mstats.trimmed_mean, limits=(trim, trim))
trim_std = partial(scipy.stats.mstats.trimmed_std, limits=(trim, trim))
#
z_val = config_bridge_z  # 6
mask = msr > msr.reduce(trim_mean, dim="ch") + z_val * msr.reduce(
    trim_std, dim="ch"
)  # bridged chans
#
bridged_ch_names = data_r_ch.ch.values[mask]
bridged_ch_names

array(['EEG 053', 'EEG 054'], dtype='<U7')

Let’s add the outlier channels to our flags

bad_chs = bridged_ch_names
pipeline.flags["ch"].add_flag_cat(kind="bridged", bad_ch_names=bad_chs)
pipeline.flags["ch"]

Flagged channels: |
  Noisy: None
  Bridged: None
  Uncorrelated: None
  Rank: None

Identify the Rank Channel

Because the pipeline uses an average reference before the ICA decomposition, it is necessary to account for rank deficiency (i.e., every sensor in the montage is linearly dependent on the other channels due to the common average reference). To account for this, the pipeline flags the sensor (out of the remaining good sensors) with the highest median of the max correlation coefficient with its neighbors (across epochs):

\[\begin{equation} i = \text{arg}\max\limits_i \rho_{i}^{t,{\text{max}}^n,median^j} \end{equation}\]

This sensor has the least unique time-series out of the remaining set of good sensors \(S_\mathcal{G}\) and is flagged by the pipeline as ”rank”. Note that this sensor is not flagged because it contains artifact, but only because one of the remaining sensors needs to be removed to address rank deficiency before ICA decomposition is performed. By choosing this sensor, we are likely to lose little information because of its high correlation with its neighbors. This sensor can be reintroduced after the ICA has been applied for artifact corrections.

good_chs = [
    ch for ch in data_r_ch.ch.values if ch not in pipeline.flags["ch"].get_flagged()
]
data_r_ch_good = data_r_ch.sel(ch=good_chs)

flag_ch = [str(data_r_ch_good.median("epoch").idxmax(dim="ch").to_numpy())]
pipeline.flags["ch"].add_flag_cat(kind="rank", bad_ch_names=flag_ch)
pipeline.flags["ch"]

Flagged channels: |
  Noisy: None
  Bridged: None
  Uncorrelated: None
  Rank: ['EEG 006']

Flag low correlation Epochs

This step is designed to identify time periods in which many sensors are uncorrelated with neighboring sensors. It is similar to the Flag Noisy Sensors step,

Again we calculate the 25th and 50th quantile of \(\mathrm{P}^{t,{\text{max}}^n}\), across the epochs dimension, and calculate the lower quantile range \(LQR^s\). This results in vectors \(\mathrm{P}^{t,{\text{max}}^nQ25^e}\) and \(\mathrm{P}^{t,{\text{max}}^nQ50^e}\) of size \(S_\mathcal{G}\). As for previous steps, we define sensor-specific thresholds for flagging epochs:

\[\begin{equation} \tau^e = \mathrm{P}^{t,{\text{max}}^nQ50^e} - LQR^e\times k \end{equation}\]

And the corresponding indicator matrix:

\[\begin{split}\begin{equation} c_{ij} = \begin{cases} 1 & \text{if } \rho^{t,{\text{max}}^n}_{ij} < \tau^e_i \\ 0 & \text{if } \rho^{t,{\text{max}}^n}_{ij} \geq \tau^e_i \end{cases} \end{equation}\end{split}\]

We average the indicator matrix across sensors and obtain a vector \(C^{\mu_s}\) that we use to flag uncorrelated epochs using the following criterion:

\[c^{\mu_e}_i > \tau^{p}.\]

Step a

q25, q50 = data_r_ch.quantile([0.25, 0.5], dim="epoch")
#
# Define the LQR
lqr = q50 - q25
#
# define a threshold
k = 3
lower_threshold = q50 - lqr * k
#
outlier_mask = data_r_ch < lower_threshold
#
percent_outliers = outlier_mask.astype(float).mean("ch")
#
p_threshold = 0.2
bad_epochs = percent_outliers[percent_outliers > p_threshold].coords.to_index().values
#
# Add the outlier epochs to our flags
pipeline.flags["epoch"].add_flag_cat(
    kind="uncorrelated", bad_epoch_inds=bad_epochs, epochs=epochs
)
pipeline.raw.annotations.description

array(['BAD_LL_noisy'], dtype='<U12')

in this case, no epochs were flagged as uncorrelated.

Flag low correlation Sensors

Flag Uncorrelated Sensors graphic.

Flag Uncorrelated Sensors. The figure shows the steps for flagging uncorrelated sensors. See the text below for descriptions of mathematical notation.

This step is designed to identify sensors that have an unusually low correlation with neighboring sensors. The operations involved by this step are similar to those of the Flag Noisy Sensors step, except we use maximal nearest neighbor correlations instead of dispersion and the left instead of the right tail of the distribution to set the threshold for outliers.

a) Take lower quantile range and defined sensor-specific thresholds

We get the indicator matrix as described previously, using

\[\tau^e_i = \mathrm{P}^{t,{\text{max}}^nQ50^e} - LQR^e\times k\]

and

\[\begin{split}c_{ij} = \begin{cases} 1 & \text{if } \rho^{t,{\text{max}}^n}_{ij} < \tau^e_i \\ 0 & \text{if } \rho^{t,{\text{max}}^n}_{ij} \geq \tau^e_i \end{cases}\end{split}\]

b) Identify uncorrelated sensors

We define a threshold as we did in the previous step and flag uncorrelated epochs \(j\) if \(c^{\mu_s}_j > \tau^{p}\).

q25, q50 = data_r_ch.quantile([0.25, 0.5], dim="ch")
#
# Define LQR
lqr = q50 - q25
#
# define a threshold
k = 3
lower_threshold = q50 - lqr * k
#
# Identify correlations less than the threshold
outlier_mask = data_r_ch < lower_threshold
percent_outliers = outlier_mask.astype(float).mean("epoch")
#
p_threshold = 0.2
bad_chs = percent_outliers[percent_outliers > p_threshold].coords.to_index().values
#
# Add the outlier channels to our flags
pipeline.flags["ch"].add_flag_cat(kind="uncorrelated", bad_ch_names=bad_chs)
pipeline.flags["ch"]

Flagged channels: |
  Noisy: None
  Bridged: None
  Uncorrelated: None
  Rank: ['EEG 006']

In this case, no sensors were flagged as uncorrelated.

Run Initial ICA

The pipeline by runs ICA two times. The first ICA is only used to identify noisy periods in its IC activation time-series. For this reason, the pipeline uses the FastICA algorithm for speed.

pipeline.run_ica("run1")

LOSSLESS: 👇 run_ica.
🧹 Epoching..
Not setting metadata
19 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 19 events and 601 original time points ...
2 bad epochs dropped
🔍 Detecting channels to leave out of reference.
/home/docs/checkouts/readthedocs.org/user_builds/pylossless/checkouts/latest/pylossless/pipeline.py:61: FutureWarning: The current default of copy=False will change to copy=True in 1.7. Set the value of copy explicitly to avoid this warning
  data = epochs.get_data()  # n_epochs, n_channels, n_times
EEG channel type selected for re-referencing
Applying a custom ('EEG',) reference.
Fitting ICA to data using 55 channels (please be patient, this may take a while)
Selecting by non-zero PCA components: 54 components
/home/docs/checkouts/readthedocs.org/user_builds/pylossless/envs/latest/lib/python3.10/site-packages/sklearn/decomposition/_fastica.py:128: ConvergenceWarning: FastICA did not converge. Consider increasing tolerance or the maximum number of iterations.
  warnings.warn(
Fitting ICA took 20.5s.
LOSSLESS: 🏁 Finished run_ica after 20.55 seconds.

Flag Noisy IC Activation time-periods

This step follows the same procedure as the Flag Noisy Sensors step, except that the data is now the IC activation time-series. thus we start with a 3D matrix \(X_{ica} \in \mathbb{R}^{I_\mathcal{G} \times E_\mathcal{G} \times T}\) of IC time-courses rather than scalp EEG data and where \(I\) is the set of independent components.

pipeline.flag_noisy_ics()

LOSSLESS: 👇 flag_noisy_ics.
🧹 Epoching..
Not setting metadata
19 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 19 events and 601 original time points ...
2 bad epochs dropped
🔍 Detecting channels to leave out of reference.
/home/docs/checkouts/readthedocs.org/user_builds/pylossless/checkouts/latest/pylossless/pipeline.py:61: FutureWarning: The current default of copy=False will change to copy=True in 1.7. Set the value of copy explicitly to avoid this warning
  data = epochs.get_data()  # n_epochs, n_channels, n_times
EEG channel type selected for re-referencing
Applying a custom ('EEG',) reference.
/home/docs/checkouts/readthedocs.org/user_builds/pylossless/checkouts/latest/pylossless/pipeline.py:64: FutureWarning: The current default of copy=False will change to copy=True in 1.7. Set the value of copy explicitly to avoid this warning
  data = ica.get_sources(epochs).get_data()
LOSSLESS: 🏁 Finished flag_noisy_ics after 0.19 seconds.

pipeline.raw.annotations.description

array(['BAD_LL_noisy'], dtype='<U12')

Run Final ICA

Now The pipeline runs the final ICA decomposition, this time using the extended Infomax algorithm. Note that any sensors or time-periods that have been flagged up to this point will not be passed into the ICA decomposition. For the sake of time, we will not run the second ICA here, as there are no more pipeline calculations.

Run ICLabel Classifier

The pipeline will run the ICLabel classifier on the final ICA, which will produce a label for each IC, one of "brain", "muscle", "eog" (eye), "ecg" (heart), line_noise, or "channel_noise".

Conclusion

And that’s all! See the other pylossless tutorials for brief examples on running the pipeline on your own data, and rejecting the flagged data.

Total running time of the script: (2 minutes 1.299 seconds)

Gallery generated by Sphinx-Gallery