From 2821954cf91360c042fb625cdedf4283c0b938e9 Mon Sep 17 00:00:00 2001 From: gcattan Date: Thu, 22 Jun 2023 12:22:36 +0200 Subject: [PATCH] Resting state with dataset and example (#400) * In some places, the virtual reality dataset code was wrong. * fix: PC data not downloading. fix: inversion 12 blocks of 5 repetitions * push example from Pedro * fix error with datframe initialization * [pre-commit.ci] auto fixes from pre-commit.com hooks * add whats new * add test * [pre-commit.ci] auto fixes from pre-commit.com hooks * fix pytest/unittest * [pre-commit.ci] auto fixes from pre-commit.com hooks * replace logging by warnings library * move docstring to the top * [pre-commit.ci] auto fixes from pre-commit.com hooks * test completed * [pre-commit.ci] auto fixes from pre-commit.com hooks * leftover * typo >< * Update examples/vr_pc_p300_different_epoch_size.py Co-authored-by: Sylvain Chevallier * rename into plot_vr_pc_p300_different_epoch_size.py * - Add figure plot - add comments * [pre-commit.ci] auto fixes from pre-commit.com hooks * Update plot_vr_pc_p300_different_epoch_size.py * [pre-commit.ci] auto fixes from pre-commit.com hooks * Create resting_state.py * push resting state * add dataset * push example * couple of bug fixes * add a condition to p300 to ignore Target/NonTarget check Fix loading of the mat file * working example * improve doc * [pre-commit.ci] auto fixes from pre-commit.com hooks * Update whats_new.rst * Update phmd_ml.py * Update plot_phmd_ml_spectrum.py flake8 * complete documentation * improve lisibility * push test * [pre-commit.ci] auto fixes from pre-commit.com hooks * fix tests * event_list missing in initialization. Correct typo. * [pre-commit.ci] auto fixes from pre-commit.com hooks * fix typo * [pre-commit.ci] auto fixes from pre-commit.com hooks * Applying and improving small details inside the tutorial * [pre-commit.ci] auto fixes from pre-commit.com hooks --------- Co-authored-by: Gregoire Cattan Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Bru Co-authored-by: Sylvain Chevallier --- docs/source/dataset_summary.rst | 13 ++++ docs/source/datasets.rst | 11 +++ docs/source/whats_new.rst | 1 + examples/plot_phmd_ml_spectrum.py | 74 ++++++++++++++++++ moabb/datasets/__init__.py | 1 + moabb/datasets/phmd_ml.py | 124 ++++++++++++++++++++++++++++++ moabb/paradigms/__init__.py | 1 + moabb/paradigms/p300.py | 22 ++++-- moabb/paradigms/resting_state.py | 81 +++++++++++++++++++ moabb/tests/paradigms.py | 48 ++++++++++++ 10 files changed, 371 insertions(+), 5 deletions(-) create mode 100644 examples/plot_phmd_ml_spectrum.py create mode 100644 moabb/datasets/phmd_ml.py create mode 100644 moabb/paradigms/resting_state.py diff --git a/docs/source/dataset_summary.rst b/docs/source/dataset_summary.rst index 0d12df98d..2f02e1276 100644 --- a/docs/source/dataset_summary.rst +++ b/docs/source/dataset_summary.rst @@ -78,6 +78,19 @@ SSVEP Wang2016,34,62,40,6,5s,250Hz,1 +Resting States +====================== + +Include neuro experiments where the participant is not actively doing something. +For example, recoding the EEG of a subject while s/he is having the eye closed or opened +is a resting state experiment. + +.. csv-table:: + :header: Dataset, #Subj, #Chan, #Classes, #Blocks / class, Trials length, Sampling rate, #Sessions + :class: sortable + + HeadMountedDisplay,12,16,2,10,60s,512Hz,1 + Submit a new dataset ~~~~~~~~~~~~~~~~~~~~ diff --git a/docs/source/datasets.rst b/docs/source/datasets.rst index c47f7cd97..345496196 100644 --- a/docs/source/datasets.rst +++ b/docs/source/datasets.rst @@ -75,6 +75,17 @@ SSVEP Datasets Lee2019_SSVEP +---------------------- +Resting State Datasets +---------------------- + +.. autosummary:: + :toctree: generated/ + :template: class.rst + + HeadMountedDisplay + + ------------ Base & Utils ------------ diff --git a/docs/source/whats_new.rst b/docs/source/whats_new.rst index a86b42834..762dcdaa7 100644 --- a/docs/source/whats_new.rst +++ b/docs/source/whats_new.rst @@ -25,6 +25,7 @@ Enhancements - Adding second deployment of the documentation (:gh:`374` by `Bruno Aristimunha`_) - Adding Parallel evaluation for :func:`moabb.evaluations.WithinSessionEvaluation` , :func:`moabb.evaluations.CrossSessionEvaluation` (:gh:`364` by `Bruno Aristimunha`_) - Add example with VirtualReality BrainInvaders dataset (:gh:`393` by `Gregoire Cattan`_ and `Pedro L. C. Rodrigues`_) +- Add resting state paradigm with dataset and example (:gh:`400` by `Gregoire Cattan`_ and `Pedro L. C. Rodrigues`_) Bugs ~~~~ diff --git a/examples/plot_phmd_ml_spectrum.py b/examples/plot_phmd_ml_spectrum.py new file mode 100644 index 000000000..722f7c91a --- /dev/null +++ b/examples/plot_phmd_ml_spectrum.py @@ -0,0 +1,74 @@ +""" +================================ +Spectral analysis of the trials +================================ + +This example demonstrates how to perform spectral +analysis on epochs extracted from a specific subject +within the :class:`moabb.datasets.HeadMountedDisplay` dataset. + +""" + +# Authors: Pedro Rodrigues +# Modified by: Gregoire Cattan +# License: BSD (3-clause) + +import warnings + +import matplotlib.pyplot as plt +import numpy as np +from scipy.signal import welch + +from moabb.datasets import HeadMountedDisplay +from moabb.paradigms import RestingStateToP300Adapter + + +warnings.filterwarnings("ignore") + +############################################################################### +# Initialization +# --------------- +# +# 1) Specify the channel and subject to compute the power spectrum. +# 2) Create an instance of the :class:`moabb.datasets.HeadMountedDisplay` dataset. +# 3) Create an instance of the :class:`moabb.paradigms.RestingStateToP300Adapter` paradigm. +# By default, the data is filtered between 1-35 Hz, +# and epochs are extracted from 10 to 50 seconds after event tagging. + +# Select channel and subject for the remaining of the example. +channel = "Cz" +subject = 1 + +dataset = HeadMountedDisplay() +events = ["on", "off"] +paradigm = RestingStateToP300Adapter(events=events, channels=[channel]) + + +############################################################################### +# Estimate Power Spectral Density +# --------------- +# 1) Obtain the epochs for the specified subject. +# 2) Use Welch's method to estimate the power spectral density. + +X, y, _ = paradigm.get_data(dataset, [subject]) +f, S = welch(X, axis=-1, nperseg=1024, fs=paradigm.resample) + +############################################################################### +# Display of the data +# --------------- +# +# Plot the averaged Power Spectral Density (PSD) for each label condition, +# using the selected channel specified at the beginning of the script. + +fig, ax = plt.subplots(facecolor="white", figsize=(8.2, 5.1)) +for condition in events: + mean_power = np.mean(S[y == condition], axis=0).flatten() + ax.plot(f, 10 * np.log10(mean_power), label=condition) + +ax.set_xlim(paradigm.fmin, paradigm.fmax) +ax.set_ylim(100, 135) +ax.set_ylabel("Spectrum Magnitude (dB)", fontsize=14) +ax.set_xlabel("Frequency (Hz)", fontsize=14) +ax.set_title("PSD for Channel " + channel, fontsize=16) +ax.legend() +fig.show() diff --git a/moabb/datasets/__init__.py b/moabb/datasets/__init__.py index 62ebc2ee4..d1ad513aa 100644 --- a/moabb/datasets/__init__.py +++ b/moabb/datasets/__init__.py @@ -35,6 +35,7 @@ from .Lee2019 import Lee2019_ERP, Lee2019_MI, Lee2019_SSVEP from .mpi_mi import MunichMI from .neiry import DemonsP300 +from .phmd_ml import HeadMountedDisplay from .physionet_mi import PhysionetMI from .schirrmeister2017 import Schirrmeister2017 from .sosulski2019 import Sosulski2019 diff --git a/moabb/datasets/phmd_ml.py b/moabb/datasets/phmd_ml.py new file mode 100644 index 000000000..b7dc6a4f2 --- /dev/null +++ b/moabb/datasets/phmd_ml.py @@ -0,0 +1,124 @@ +import os + +import mne +import numpy as np +from scipy.io import loadmat + +from . import download as dl +from .base import BaseDataset + + +HEADMOUNTED_URL = "https://zenodo.org/record/2617085/files/" + + +class HeadMountedDisplay(BaseDataset): + """ + Passive Head Mounted Display with Music Listening dataset. + + .. admonition:: Dataset summary + + + ================= ======= ======= ========== ================= ============ =============== =========== + Name #Subj #Chan #Classes #Blocks/class Trials len Sampling rate #Sessions + ================== ======= ======= ========== ================= ============ =============== =========== + HeadMountedDisplay 12 16 2 10 60s 512Hz 1 + ================== ======= ======= ========== ================= ============ =============== =========== + + We describe the experimental procedures for a dataset that we have made publicly available + at https://doi.org/10.5281/zenodo.2617084 in mat (Mathworks, Natick, USA) and csv formats. + This dataset contains electroencephalographic recordings of 12 subjects listening to music + with and without a passive head-mounted display, that is, a head-mounted display which does + not include any electronics at the exception of a smartphone. The electroencephalographic + headset consisted of 16 electrodes. Data were recorded during a pilot experiment taking + place in the GIPSA-lab, Grenoble, France, in 2017 (Cattan and al, 2018). + The ID of this dataset is PHMDML.EEG.2017-GIPSA. + + **full description of the experiment** + https://hal.archives-ouvertes.fr/hal-02085118 + + **Link to the data** + https://doi.org/10.5281/zenodo.2617084 + + **Authors** + Principal Investigator: Eng. Grégoire Cattan + Technical Supervisors: Eng. Pedro L. C. Rodrigues + Scientific Supervisor: Dr. Marco Congedo + + **ID of the dataset** + PHMDML.EEG.2017-GIPSA + + Notes + ----- + + .. versionadded:: 0.6.0 + + References + ---------- + + .. [1] G. Cattan, P. L. Coelho Rodrigues, and M. Congedo, + ‘Passive Head-Mounted Display Music-Listening EEG dataset’, + Gipsa-Lab ; IHMTEK, Research Report 2, Mar. 2019. doi: 10.5281/zenodo.2617084. + + + """ + + def __init__(self): + super().__init__( + subjects=list(range(1, 12 + 1)), + sessions_per_subject=1, + events=dict(on=1, off=2), + code="PHMD-ML", + interval=[0, 1], + paradigm="rstate", + doi="https://doi.org/10.5281/zenodo.2617084 ", + ) + self._chnames = [ + "Fp1", + "Fp2", + "Fc5", + "Fz", + "Fc6", + "T7", + "Cz", + "T8", + "P7", + "P3", + "Pz", + "P4", + "P8", + "O1", + "Oz", + "O2", + "stim", + ] + self._chtypes = ["eeg"] * 16 + ["stim"] + + def _get_single_subject_data(self, subject): + """return data for a single subject""" + + filepath = self.data_path(subject)[0] + data = loadmat(os.path.join(filepath, os.listdir(filepath)[0])) + + first_channel = 1 + last_channel = 17 + S = data["data"][:, first_channel:last_channel] + stim = data["data"][:, -1] + + X = np.concatenate([S, stim[:, None]], axis=1).T + + info = mne.create_info( + ch_names=self._chnames, sfreq=512, ch_types=self._chtypes, verbose=False + ) + raw = mne.io.RawArray(data=X, info=info, verbose=False) + return {"session_0": {"run_0": raw}} + + def data_path( + self, subject, path=None, force_update=False, update_path=None, verbose=None + ): + if subject not in self.subject_list: + raise (ValueError("Invalid subject number")) + + url = "{:s}subject_{:02d}.mat".format(HEADMOUNTED_URL, subject) + file_path = dl.data_path(url, "HEADMOUNTED") + + return [file_path] diff --git a/moabb/paradigms/__init__.py b/moabb/paradigms/__init__.py index 055e278a3..a1dac30d8 100644 --- a/moabb/paradigms/__init__.py +++ b/moabb/paradigms/__init__.py @@ -9,4 +9,5 @@ # flake8: noqa from moabb.paradigms.p300 import * +from moabb.paradigms.resting_state import * from moabb.paradigms.ssvep import * diff --git a/moabb/paradigms/p300.py b/moabb/paradigms/p300.py index b3b8def71..09650ad61 100644 --- a/moabb/paradigms/p300.py +++ b/moabb/paradigms/p300.py @@ -168,11 +168,15 @@ def process_raw( # noqa: C901 # pick events, based on event_id try: - if type(event_id["Target"]) is list and type(event_id["NonTarget"]) == list: - event_id_new = dict(Target=1, NonTarget=0) - events = mne.merge_events(events, event_id["Target"], 1) - events = mne.merge_events(events, event_id["NonTarget"], 0) - event_id = event_id_new + if "Target" in event_id and "NonTarget" in event_id: + if ( + type(event_id["Target"]) is list + and type(event_id["NonTarget"]) == list + ): + event_id_new = dict(Target=1, NonTarget=0) + events = mne.merge_events(events, event_id["Target"], 1) + events = mne.merge_events(events, event_id["NonTarget"], 0) + event_id = event_id_new events = mne.pick_events(events, include=list(event_id.values())) except RuntimeError: # skip raw if no event found @@ -317,6 +321,14 @@ def __init__(self, fmin=1, fmax=24, **kwargs): raise (ValueError("P300 does not take argument filters")) super().__init__(filters=[[fmin, fmax]], **kwargs) + @property + def fmax(self): + return self.filters[0][1] + + @property + def fmin(self): + return self.filters[0][0] + class P300(SinglePass): """P300 for Target/NonTarget classification diff --git a/moabb/paradigms/resting_state.py b/moabb/paradigms/resting_state.py new file mode 100644 index 000000000..f41ab4c67 --- /dev/null +++ b/moabb/paradigms/resting_state.py @@ -0,0 +1,81 @@ +"""Resting state Paradigms + +Regroups paradigms for experience where we record the EEG +and the participant is not doing an active task, such +as focusing, counting or speaking. + +Typically, a open/close eye experiment, where we +record the EEG of a subject while he is having the eye open or close +is a resting state experiment. +""" + +from moabb.paradigms.p300 import SinglePass + + +class RestingStateToP300Adapter(SinglePass): + """Adapter to the P300 paradigm for resting state experiments. + It implements a SinglePass processing as for P300, except that: + - the name of the event is free (it is not enforced to Target/NonTarget as for P300) + - the default values are different. In particular, the length of the epochs is larger. + + Parameters + ---------- + fmin: float (default 1) + cutoff frequency (Hz) for the high pass filter + + fmax: float (default 35) + cutoff frequency (Hz) for the low pass filter + + events: List of str | None (default None) + event to use for epoching. If None, default to all events defined in + the dataset. + + tmin: float (default 10s) + Start time (in second) of the epoch, relative to the dataset specific + task interval e.g. tmin = 1 would mean the epoch will start 1 second + after the beginning of the task as defined by the dataset. + + tmax: float | None, (default 50s) + End time (in second) of the epoch, relative to the beginning of the + dataset specific task interval. tmax = 5 would mean the epoch will end + 5 second after the beginning of the task as defined in the dataset. If + None, use the dataset value. + + resample: float | None (default 128) + If not None, resample the eeg data with the sampling rate provided. + + baseline: None | tuple of length 2 + The time interval to consider as “baseline” when applying baseline + correction. If None, do not apply baseline correction. + If a tuple (a, b), the interval is between a and b (in seconds), + including the endpoints. + Correction is applied by computing the mean of the baseline period + and subtracting it from the data (see mne.Epochs) + + channels: list of str | None (default None) + list of channel to select. If None, use all EEG channels available in + the dataset. + """ + + def __init__(self, fmin=1, fmax=35, tmin=10, tmax=50, resample=128, **kwargs): + super().__init__( + fmin=fmin, fmax=fmax, tmin=tmin, tmax=tmax, resample=resample, **kwargs + ) + + def used_events(self, dataset): + return {ev: dataset.event_id[ev] for ev in self.events} + + def is_valid(self, dataset): + ret = True + if not (dataset.paradigm == "rstate"): + ret = False + + if self.events: + if not set(self.events) <= set(dataset.event_id.keys()): + ret = False + + return ret + + @property + def scoring(self): + return "roc_auc" diff --git a/moabb/tests/paradigms.py b/moabb/tests/paradigms.py index 4917aa3c5..f69f4072f 100644 --- a/moabb/tests/paradigms.py +++ b/moabb/tests/paradigms.py @@ -16,6 +16,7 @@ FilterBankMotorImagery, FilterBankSSVEP, LeftRightImagery, + RestingStateToP300Adapter, ) @@ -329,6 +330,53 @@ def test_P300_paradigm(self): self.assertIsInstance(epochs, BaseEpochs) +class Test_RestingState(unittest.TestCase): + def test_RestingState_paradigm(self): + event_list = ["Open", "Close"] + paradigm = RestingStateToP300Adapter(events=event_list) + dataset = FakeDataset(paradigm="rstate", event_list=event_list) + X, labels, metadata = paradigm.get_data(dataset, subjects=[1]) + + # we should have all the same length + self.assertEqual(len(X), len(labels), len(metadata)) + # X must be a 3D Array + self.assertEqual(len(X.shape), 3) + # labels must contain 2 values (Open/Close) + self.assertEqual(len(np.unique(labels)), 2) + # metadata must have subjets, sessions, runs + self.assertTrue("subject" in metadata.columns) + self.assertTrue("session" in metadata.columns) + self.assertTrue("run" in metadata.columns) + # we should have only one subject in the metadata + self.assertEqual(np.unique(metadata.subject), 1) + # we should have two sessions in the metadata + self.assertEqual(len(np.unique(metadata.session)), 2) + # should return epochs + epochs, _, _ = paradigm.get_data(dataset, subjects=[1], return_epochs=True) + self.assertIsInstance(epochs, BaseEpochs) + # should return raws + raws, _, _ = paradigm.get_data(dataset, subjects=[1], return_raws=True) + for raw in raws: + self.assertIsInstance(raw, BaseRaw) + # should raise error + self.assertRaises( + ValueError, + paradigm.get_data, + dataset, + subjects=[1], + return_epochs=True, + return_raws=True, + ) + + def test_RestingState_default_values(self): + paradigm = RestingStateToP300Adapter() + assert paradigm.tmin == 10 + assert paradigm.tmax == 50 + assert paradigm.fmin == 1 + assert paradigm.fmax == 35 + assert paradigm.resample == 128 + + class Test_SSVEP(unittest.TestCase): def test_BaseSSVEP_paradigm(self): paradigm = BaseSSVEP(n_classes=None)