diff --git a/docs/changelog.rst b/docs/changelog.rst index 43c52920..ed0434a4 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -27,6 +27,7 @@ Version 0.7.0 (unreleased) Changelog ~~~~~~~~~ - :meth:`~pyprep.NoisyChannels.find_bad_by_nan_flat` now respects the ``reject_by_annotation`` setting of the :class:`~pyprep.NoisyChannels` instance, consistent with all other detection methods, by `Roy Eric Wieske`_ (:gh:`185`) +- Removed the unmaintained, non-asserting matplotlib-based PREP tests and the ``run_full_prep`` example (along with its committed MATLAB ``.mat`` artifacts); PyPREP's fidelity to MATLAB PREP is validated by ``tests/test_matprep_compare.py`` instead, by `Stefan Appelhoff`_ (:gh:`163`) Bug ~~~ diff --git a/examples/matlab_results/EEG.mat b/examples/matlab_results/EEG.mat deleted file mode 100644 index a7e83c48..00000000 Binary files a/examples/matlab_results/EEG.mat and /dev/null differ diff --git a/examples/matlab_results/EEGNew.mat b/examples/matlab_results/EEGNew.mat deleted file mode 100644 index 1bcd60c1..00000000 Binary files a/examples/matlab_results/EEGNew.mat and /dev/null differ diff --git a/examples/matlab_results/EEG_raw.mat b/examples/matlab_results/EEG_raw.mat deleted file mode 100644 index 7de2db90..00000000 Binary files a/examples/matlab_results/EEG_raw.mat and /dev/null differ diff --git a/examples/matlab_results/EEGinterp.mat b/examples/matlab_results/EEGinterp.mat deleted file mode 100644 index 0f133ea3..00000000 Binary files a/examples/matlab_results/EEGinterp.mat and /dev/null differ diff --git a/examples/matlab_results/EEGref.mat b/examples/matlab_results/EEGref.mat deleted file mode 100644 index 8ba9a293..00000000 Binary files a/examples/matlab_results/EEGref.mat and /dev/null differ diff --git a/examples/run_full_prep.py b/examples/run_full_prep.py deleted file mode 100644 index 03c31bb9..00000000 --- a/examples/run_full_prep.py +++ /dev/null @@ -1,326 +0,0 @@ -""" -================= -Run the full PREP -================= - - -In this example we show how to run PREP with ``pyprep``. We also compare -:class:`~pyprep.PrepPipeline` with PREP's results in Matlab. - -We use sample EEG data from Physionet EEG Motor Movement/Imagery Dataset: -`https://physionet.org/content/eegmmidb/1.0.0/ -`_ - -.. currentmodule:: pyprep -""" # noqa: D205 D400 - -# Authors: The PyPREP developers -# SPDX-License-Identifier: MIT - -############################################################################### -# -# .. warning:: This functionality is work in progress. -# Contributions are welcome! -# - -############################################################################### -# First we import what we need for this example. -from pathlib import Path - -import matplotlib.pyplot as plt -import mne -import numpy as np -import scipy.io as sio - -from pyprep.prep_pipeline import PrepPipeline - -############################################################################### -# Let's download some data for testing. Picking the 1st run of subject 4 here. - -data_paths = mne.datasets.eegbci.load_data(4, 1, update_path=True) - -############################################################################### -# General settings and file paths -mne.set_log_level("WARNING") - -# Raw data -fname_test_file = data_paths[0] - -# mat files for validation -here = Path("__file__").parent.absolute() - -fname_mat1 = here / "matlab_results" / "EEG_raw.mat" -fname_mat2 = here / "matlab_results" / "EEGNew.mat" -fname_mat3 = here / "matlab_results" / "EEG.mat" -fname_mat4 = here / "matlab_results" / "EEGref.mat" -fname_mat5 = here / "matlab_results" / "EEGinterp.mat" - -############################################################################### -# Load data and prepare it -# ------------------------ - -raw = mne.io.read_raw_edf(fname_test_file, preload=True) - -# The eegbci data has non-standard channel names. We need to rename them: -mne.datasets.eegbci.standardize(raw) - -# Add a montage to the data -montage_kind = "standard_1005" -montage = mne.channels.make_standard_montage(montage_kind) - -# Extract some info -sample_rate = raw.info["sfreq"] - -# Make a copy of the data -raw_copy = raw.copy() - -############################################################################### -# Set PREP parameters and run PREP -# -------------------------------- -# -# Notes: We keep all the default parameter settings as described in the PREP, -# though for this case we think that the default fraction of bad time windows -# being 0.01 is too sensitive since it gots only 60 time windows (the EEG data -# is 60s long). As a result this example returns a lot of interpolated channels. - -# Fit prep -prep_params = { - "ref_chs": "eeg", - "reref_chs": "eeg", - "line_freqs": np.arange(60, sample_rate / 2, 60), -} - -prep = PrepPipeline(raw_copy, prep_params, montage) -prep.fit() - -############################################################################### -# Results -# ------- -# -# You can check the detected bad channels in each step of PREP. - -print(f"Bad channels: {prep.interpolated_channels}") -print(f"Bad channels original: {prep.noisy_channels_original['bad_all']}") -print(f"Bad channels after interpolation: {prep.still_noisy_channels}") - -# Matlab's results -# ---------------- -# Bad channels: Fc5, Fc3, Fc1, C3, Cp3, Cp4, Af3, Afz, Af8, F7, F5, F6, F8, -# Ft8, P7, P2 -# Bad channels original: Af3, Af4, Af7, Af8, Fp1, Fp2, Fpz, Ft8 -# Bad channels after interpolation: Cp5, Fp2, Af7, F1 - -############################################################################### -# Validation -# ---------- -# -# To validate each step of pyprep's results, we compare results after each step -# with the results from EEGLAB's PREP. To make it easy to compare, we rescale -# the EEG data to [-1, 1] (divided the data by maximum absolute value) when -# making the plot. - -EEG_raw = raw_copy.get_data(picks="eeg") * 1e6 -EEG_raw_max = np.max(abs(EEG_raw), axis=None) -EEG_raw_matlab = sio.loadmat(fname_mat1) -EEG_raw_matlab = EEG_raw_matlab["save_data"] -EEG_raw_diff = EEG_raw - EEG_raw_matlab -EEG_raw_mse = (EEG_raw_diff / EEG_raw_max**2).mean(axis=None) - -fig, axs = plt.subplots(5, 3, sharex="all", figsize=(16, 12)) -plt.setp(fig, facecolor=[1, 1, 1]) -fig.suptitle("Python versus Matlab PREP results", fontsize=16) - -im = axs[0, 0].imshow( - EEG_raw / EEG_raw_max, - aspect="auto", - extent=[0, (EEG_raw.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[0, 0].set_title("Python", fontsize=14) -axs[0, 1].imshow( - EEG_raw_matlab / EEG_raw_max, - aspect="auto", - extent=[0, (EEG_raw_matlab.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[0, 1].set_title("Matlab", fontsize=14) -axs[0, 2].imshow( - EEG_raw_diff / EEG_raw_max, - aspect="auto", - extent=[0, (EEG_raw_diff.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[0, 2].set_title("Difference", fontsize=14) -axs[0, 1].set_title("Original EEG", fontsize=14) -# axs[0, 0].set_ylabel('Channel Number', fontsize=14) -cb = fig.colorbar(im, ax=axs, fraction=0.05, pad=0.04) -cb.set_label("\u03bcVolt", fontsize=14) - -EEG_new_matlab = sio.loadmat(fname_mat2) -EEG_new_matlab = EEG_new_matlab["save_data"] -EEG_new = prep.EEG_new -EEG_new_max = np.max(abs(EEG_new), axis=None) -EEG_new_diff = EEG_new - EEG_new_matlab -EEG_new_mse = ((EEG_new_diff / EEG_new_max) ** 2).mean(axis=None) -axs[1, 0].imshow( - EEG_new / EEG_new_max, - aspect="auto", - extent=[0, (EEG_new.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[1, 1].imshow( - EEG_new_matlab / EEG_new_max, - aspect="auto", - extent=[0, (EEG_new_matlab.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[1, 2].imshow( - EEG_new_diff / EEG_new_max, - aspect="auto", - extent=[0, (EEG_new_diff.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[1, 1].set_title("High pass filtered EEG", fontsize=14) -# axs[1, 0].set_ylabel('Channel Number', fontsize=14) - -EEG_clean_matlab = sio.loadmat(fname_mat3) -EEG_clean_matlab = EEG_clean_matlab["save_data"] -EEG_clean = prep.EEG -EEG_max = np.max(abs(EEG_clean), axis=None) -EEG_diff = EEG_clean - EEG_clean_matlab -EEG_mse = ((EEG_diff / EEG_max) ** 2).mean(axis=None) -axs[2, 0].imshow( - EEG_clean / EEG_max, - aspect="auto", - extent=[0, (EEG_clean.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[2, 1].imshow( - EEG_clean_matlab / EEG_max, - aspect="auto", - extent=[0, (EEG_clean_matlab.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[2, 2].imshow( - EEG_diff / EEG_max, - aspect="auto", - extent=[0, (EEG_diff.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[2, 1].set_title("Line-noise removed EEG", fontsize=14) -axs[2, 0].set_ylabel("Channel Number", fontsize=14) - -EEG = prep.EEG_before_interpolation -EEG_max = np.max(abs(EEG), axis=None) -EEG_ref_mat = sio.loadmat(fname_mat4) -EEG_ref_matlab = EEG_ref_mat["save_EEG"] -reference_matlab = EEG_ref_mat["save_reference"] -EEG_ref_diff = EEG - EEG_ref_matlab -EEG_ref_mse = ((EEG_ref_diff / EEG_max) ** 2).mean(axis=None) -reference_signal = prep.reference_before_interpolation -reference_max = np.max(abs(reference_signal), axis=None) -reference_diff = reference_signal - reference_matlab -reference_mse = ((reference_diff / reference_max) ** 2).mean(axis=None) -axs[3, 0].imshow( - EEG / EEG_max, - aspect="auto", - extent=[0, (EEG.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[3, 1].imshow( - EEG_ref_matlab / EEG_max, - aspect="auto", - extent=[0, (EEG_ref_matlab.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[3, 2].imshow( - EEG_ref_diff / EEG_max, - aspect="auto", - extent=[0, (EEG_ref_diff.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[3, 1].set_title("Referenced EEG", fontsize=14) -# axs[3, 0].set_ylabel('Channel Number', fontsize=14) - -EEG_final = prep.raw.get_data() * 1e6 -EEG_final_max = np.max(abs(EEG_final), axis=None) -EEG_final_matlab = sio.loadmat(fname_mat5) -EEG_final_matlab = EEG_final_matlab["save_data"] -EEG_final_diff = EEG_final - EEG_final_matlab -EEG_final_mse = ((EEG_final_diff / EEG_final_max) ** 2).mean(axis=None) -axs[4, 0].imshow( - EEG_final / EEG_final_max, - aspect="auto", - extent=[0, (EEG_final.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[4, 1].imshow( - EEG_final_matlab / EEG_final_max, - aspect="auto", - extent=[0, (EEG_final_matlab.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[4, 2].imshow( - EEG_final_diff / EEG_final_max, - aspect="auto", - extent=[0, (EEG_final_diff.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), -) -axs[4, 1].set_title("Interpolated EEG", fontsize=14) -# axs[4, 0].set_ylabel('Channel Number', fontsize=14) -axs[4, 1].set_xlabel("Time(s)", fontsize=14) - -############################################################################### -# Mean square error of each step: - -print(f"Raw data MSE: {EEG_raw_mse}") -print(f"Filtered data MSE: {EEG_new_mse}") -print(f"Line-noise removed data MSE: {EEG_mse}") -print(f"Referenced data MSE: {EEG_ref_mse}") -print(f"Interpolated data MSE: {EEG_final_mse}") - -############################################################################### -# Discussion -# ---------- -# -# It can be seen the results match well on each step except the final step. -# This is due to the difference of find_noisy_channel functions, since the -# channels with relatively large error corresponds to the channels that are only -# interpolated in Python or Matlab. -# -# We think the differences mainly arise from -# -# 1. Difference in bad channels from Ransac criteria, including the random -# number generator -# 2. Difference in some internal functions of Python and Matlab (e.g., filter -# and interpolation function) diff --git a/tests/test_prep_pipeline.py b/tests/test_prep_pipeline.py index b6c3cbb9..bcabd20d 100644 --- a/tests/test_prep_pipeline.py +++ b/tests/test_prep_pipeline.py @@ -3,216 +3,14 @@ # Authors: The PyPREP developers # SPDX-License-Identifier: MIT -import matplotlib.pyplot as plt import numpy as np import pytest -import scipy.io as sio from pyprep.prep_pipeline import PrepPipeline from .conftest import make_random_mne_object -@pytest.mark.usefixtures("raw", "montage") -def test_prep_pipeline(raw, montage): - """Test prep pipeline.""" - eeg_index = np.array( - [idx for idx, typ in enumerate(raw.get_channel_types()) if typ == "eeg"] - ) - raw_copy = raw.copy() - ch_names = raw_copy.info["ch_names"] - ch_names_eeg = list(np.asarray(ch_names)[eeg_index]) - sample_rate = raw_copy.info["sfreq"] - prep_params = { - "ref_chs": ch_names_eeg, - "reref_chs": ch_names_eeg, - "line_freqs": np.arange(60, sample_rate / 2, 60), - } - prep = PrepPipeline(raw_copy, prep_params, montage, random_state=42) - prep.fit() - - EEG_raw = raw_copy.get_data(picks="eeg") * 1e6 - EEG_raw_max = np.max(abs(EEG_raw), axis=None) - EEG_raw_matlab = sio.loadmat("./examples/matlab_results/EEG_raw.mat") - EEG_raw_matlab = EEG_raw_matlab["save_data"] - EEG_raw_diff = EEG_raw - EEG_raw_matlab - # EEG_raw_mse = (EEG_raw_diff / EEG_raw_max ** 2).mean(axis=None) - - fig, axs = plt.subplots(5, 3, sharex="all") - plt.setp(fig, facecolor=[1, 1, 1]) - fig.suptitle("Python versus Matlab PREP results", fontsize=16) - - im = axs[0, 0].imshow( - EEG_raw / EEG_raw_max, - aspect="auto", - extent=[0, (EEG_raw.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - axs[0, 0].set_title("Python", fontsize=14) - axs[0, 1].imshow( - EEG_raw_matlab / EEG_raw_max, - aspect="auto", - extent=[0, (EEG_raw_matlab.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - axs[0, 1].set_title("Matlab", fontsize=14) - axs[0, 2].imshow( - EEG_raw_diff / EEG_raw_max, - aspect="auto", - extent=[0, (EEG_raw_diff.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - axs[0, 2].set_title("Difference", fontsize=14) - # axs[0, 0].set_title('Original EEG', loc='left', fontsize=14) - # axs[0, 0].set_ylabel('Channel Number', fontsize=14) - cb = fig.colorbar(im, ax=axs, fraction=0.05, pad=0.04) - cb.set_label("\u03bcVolt", fontsize=14) - - EEG_new_matlab = sio.loadmat("./examples/matlab_results/EEGNew.mat") - EEG_new_matlab = EEG_new_matlab["save_data"] - EEG_new = prep.EEG_new - EEG_new_max = np.max(abs(EEG_new), axis=None) - EEG_new_diff = EEG_new - EEG_new_matlab - # EEG_new_mse = ((EEG_new_diff / EEG_new_max) ** 2).mean(axis=None) - axs[1, 0].imshow( - EEG_new / EEG_new_max, - aspect="auto", - extent=[0, (EEG_new.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - axs[1, 1].imshow( - EEG_new_matlab / EEG_new_max, - aspect="auto", - extent=[0, (EEG_new_matlab.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - axs[1, 2].imshow( - EEG_new_diff / EEG_new_max, - aspect="auto", - extent=[0, (EEG_new_diff.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - # axs[1, 0].set_title('High pass filter', loc='left', fontsize=14) - # axs[1, 0].set_ylabel('Channel Number', fontsize=14) - - EEG_clean_matlab = sio.loadmat("./examples/matlab_results/EEG.mat") - EEG_clean_matlab = EEG_clean_matlab["save_data"] - EEG_clean = prep.EEG - EEG_max = np.max(abs(EEG_clean), axis=None) - EEG_diff = EEG_clean - EEG_clean_matlab - # EEG_mse = ((EEG_diff / EEG_max) ** 2).mean(axis=None) - axs[2, 0].imshow( - EEG_clean / EEG_max, - aspect="auto", - extent=[0, (EEG_clean.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - axs[2, 1].imshow( - EEG_clean_matlab / EEG_max, - aspect="auto", - extent=[0, (EEG_clean_matlab.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - axs[2, 2].imshow( - EEG_diff / EEG_max, - aspect="auto", - extent=[0, (EEG_diff.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - # axs[2, 0].set_title('Line-noise removal', loc='left', fontsize=14) - axs[2, 0].set_ylabel("Channel Number", fontsize=14) - - EEG = prep.EEG_before_interpolation - EEG_max = np.max(abs(EEG), axis=None) - EEG_ref_mat = sio.loadmat("./examples/matlab_results/EEGref.mat") - EEG_ref_matlab = EEG_ref_mat["save_EEG"] - # reference_matlab = EEG_ref_mat["save_reference"] - EEG_ref_diff = EEG - EEG_ref_matlab - # EEG_ref_mse = ((EEG_ref_diff / EEG_max) ** 2).mean(axis=None) - # reference_signal = prep.reference_before_interpolation - # reference_max = np.max(abs(reference_signal), axis=None) - # reference_diff = reference_signal - reference_matlab - # reference_mse = ((reference_diff / reference_max) ** 2).mean(axis=None) - axs[3, 0].imshow( - EEG / EEG_max, - aspect="auto", - extent=[0, (EEG.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - axs[3, 1].imshow( - EEG_ref_matlab / EEG_max, - aspect="auto", - extent=[0, (EEG_ref_matlab.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - axs[3, 2].imshow( - EEG_ref_diff / EEG_max, - aspect="auto", - extent=[0, (EEG_ref_diff.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - # axs[3, 0].set_title('Referencing', loc='left', fontsize=14) - # axs[3, 0].set_ylabel('Channel Number', fontsize=14) - - EEG_final = prep.raw.get_data() * 1e6 - EEG_final_max = np.max(abs(EEG_final), axis=None) - EEG_final_matlab = sio.loadmat("./examples/matlab_results/EEGinterp.mat") - EEG_final_matlab = EEG_final_matlab["save_data"] - EEG_final_diff = EEG_final - EEG_final_matlab - # EEG_final_mse = ((EEG_final_diff / EEG_final_max) ** 2).mean(axis=None) - axs[4, 0].imshow( - EEG_final / EEG_final_max, - aspect="auto", - extent=[0, (EEG_final.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - axs[4, 1].imshow( - EEG_final_matlab / EEG_final_max, - aspect="auto", - extent=[0, (EEG_final_matlab.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - axs[4, 2].imshow( - EEG_final_diff / EEG_final_max, - aspect="auto", - extent=[0, (EEG_final_diff.shape[1] / sample_rate), 63, 0], - vmin=-1, - vmax=1, - cmap=plt.get_cmap("RdBu"), - ) - # axs[4, 0].set_title('Interpolation', loc='left', fontsize=14) - # axs[4, 0].set_ylabel('Channel Number', fontsize=14) - axs[4, 1].set_xlabel("Time(s)", fontsize=14) - - @pytest.mark.usefixtures("raw", "montage") def test_prep_pipeline_non_eeg(raw, montage): """Test prep pipeline with non eeg channels.""" @@ -280,3 +78,8 @@ def test_prep_pipeline_filter_kwargs(raw, montage): raw_copy, prep_params, montage, random_state=42, filter_kwargs=filter_kwargs ) prep.fit() + + # The input is purely EEG, so there are no non-EEG channels to split off + # and the `raw` property returns the EEG data unchanged. + assert prep.raw_non_eeg is None + assert prep.raw.ch_names == prep.raw_eeg.ch_names