From 21a185fe8c680da6aced6e6ba06309a140082842 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 4 Feb 2025 19:30:44 +0000 Subject: [PATCH 01/13] [pre-commit.ci] pre-commit autoupdate (#13099) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- .pre-commit-config.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index fb5a6bd4247..5abf98c9945 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,7 +1,7 @@ repos: # Ruff mne - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.9.3 + rev: v0.9.4 hooks: - id: ruff name: ruff lint mne @@ -23,7 +23,7 @@ repos: # Codespell - repo: https://github.com/codespell-project/codespell - rev: v2.4.0 + rev: v2.4.1 hooks: - id: codespell additional_dependencies: @@ -82,7 +82,7 @@ repos: # zizmor - repo: https://github.com/woodruffw/zizmor-pre-commit - rev: v1.2.2 + rev: v1.3.0 hooks: - id: zizmor From afd69afe9ed9b07d48e860537d37fc349def6635 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Wed, 5 Feb 2025 20:21:44 +0200 Subject: [PATCH 02/13] Take units (m or mm) into account when showing fieldmaps on top of brains (#13101) Co-authored-by: Eric Larson --- azure-pipelines.yml | 1 + doc/changes/devel/13101.bugfix.rst | 1 + mne/conftest.py | 10 +++++++--- mne/viz/evoked_field.py | 6 +++++- mne/viz/tests/test_3d.py | 18 +++++++++++++++--- 5 files changed, 29 insertions(+), 7 deletions(-) create mode 100644 doc/changes/devel/13101.bugfix.rst diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 7149edac50b..9c42b3286c1 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -89,6 +89,7 @@ stages: DISPLAY: ':99' OPENBLAS_NUM_THREADS: '1' MNE_TEST_ALLOW_SKIP: '^.*(PySide6 causes segfaults).*$' + MNE_BROWSER_PRECOMPUTE: 'false' steps: - bash: | set -e diff --git a/doc/changes/devel/13101.bugfix.rst b/doc/changes/devel/13101.bugfix.rst new file mode 100644 index 00000000000..d24e55b5056 --- /dev/null +++ b/doc/changes/devel/13101.bugfix.rst @@ -0,0 +1 @@ +Take units (m or mm) into account when drawing :func:`~mne.viz.plot_evoked_field` on top of :class:`~mne.viz.Brain`, by `Marijn van Vliet`_. diff --git a/mne/conftest.py b/mne/conftest.py index 2a73c7a1b8e..fabf9f93e5a 100644 --- a/mne/conftest.py +++ b/mne/conftest.py @@ -6,6 +6,7 @@ import inspect import os import os.path as op +import platform import re import shutil import sys @@ -190,6 +191,8 @@ def pytest_configure(config: pytest.Config): ignore:Starting field name with a underscore.*: # joblib ignore:process .* is multi-threaded, use of fork/exec.*:DeprecationWarning + # sklearn + ignore:Python binding for RankQuantileOptions.*:RuntimeWarning """ # noqa: E501 for warning_line in warning_lines.split("\n"): warning_line = warning_line.strip() @@ -285,9 +288,10 @@ def __init__(self, exception_handler=None, signals=None): @pytest.fixture(scope="session") def azure_windows(): """Determine if running on Azure Windows.""" - return os.getenv( - "AZURE_CI_WINDOWS", "false" - ).lower() == "true" and sys.platform.startswith("win") + return ( + os.getenv("AZURE_CI_WINDOWS", "false").lower() == "true" + and platform.system() == "Windows" + ) @pytest.fixture(scope="function") diff --git a/mne/viz/evoked_field.py b/mne/viz/evoked_field.py index cf5a9996216..839259ee117 100644 --- a/mne/viz/evoked_field.py +++ b/mne/viz/evoked_field.py @@ -7,6 +7,7 @@ # License: BSD-3-Clause # Copyright the MNE-Python contributors. +from copy import deepcopy from functools import partial import numpy as np @@ -185,6 +186,7 @@ def __init__( if isinstance(fig, Brain): self._renderer = fig._renderer self._in_brain_figure = True + self._units = fig._units if _get_3d_backend() == "notebook": raise NotImplementedError( "Plotting on top of an existing Brain figure " @@ -195,6 +197,7 @@ def __init__( fig, bgcolor=(0.0, 0.0, 0.0), size=(600, 600) ) self._in_brain_figure = False + self._units = "m" self.plotter = self._renderer.plotter self.interaction = interaction @@ -277,7 +280,8 @@ def _prepare_surf_map(self, surf_map, color, alpha): # Make a solid surface surf = surf_map["surf"] - if self._in_brain_figure: + if self._units == "mm": + surf = deepcopy(surf) surf["rr"] *= 1000 map_vmax = self._vmax.get(surf_map["kind"]) if map_vmax is None: diff --git a/mne/viz/tests/test_3d.py b/mne/viz/tests/test_3d.py index 34022d59768..97b10621108 100644 --- a/mne/viz/tests/test_3d.py +++ b/mne/viz/tests/test_3d.py @@ -191,10 +191,21 @@ def test_plot_evoked_field(renderer): ch_type=t, ) evoked.plot_field(maps, time=0.1, n_contours=n_contours) + renderer.backend._close_all() - # Test plotting inside an existing Brain figure. - brain = Brain("fsaverage", "lh", "inflated", subjects_dir=subjects_dir) - fig = evoked.plot_field(maps, time=0.1, fig=brain) + # Test plotting inside an existing Brain figure. Check that units are taken into + # account. + for units in ["mm", "m"]: + brain = Brain( + "fsaverage", "lh", "inflated", units=units, subjects_dir=subjects_dir + ) + fig = evoked.plot_field(maps, time=0.1, fig=brain) + assert brain._units == fig._units + scale = 1000 if units == "mm" else 1 + assert ( + fig._surf_maps[0]["surf"]["rr"][0, 0] == scale * maps[0]["surf"]["rr"][0, 0] + ) + renderer.backend._close_all() # Test some methods fig = evoked.plot_field(maps, time_viewer=True) @@ -214,6 +225,7 @@ def test_plot_evoked_field(renderer): fig = evoked.plot_field(maps, time_viewer=False) assert isinstance(fig, Figure3D) + renderer.backend._close_all() @testing.requires_testing_data From 47b036953128b1fdc8e7d9631a89db9b296dce73 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Fri, 7 Feb 2025 11:25:03 -0500 Subject: [PATCH 03/13] MAINT: Use statsmodels pre and fix CircleCI (#13106) --- doc/sphinxext/related_software.py | 4 ++++ tools/install_pre_requirements.sh | 8 +------- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/doc/sphinxext/related_software.py b/doc/sphinxext/related_software.py index 2548725390a..a739e545f8f 100644 --- a/doc/sphinxext/related_software.py +++ b/doc/sphinxext/related_software.py @@ -85,6 +85,10 @@ "Home-page": "https://eeg-positions.readthedocs.io", "Summary": "Compute and plot standard EEG electrode positions.", }, + "mne-faster": { + "Home-page": "https://github.com/wmvanvliet/mne-faster", + "Summary": "MNE-FASTER: automatic bad channel/epoch/component detection.", # noqa: E501 + }, "mne-features": { "Home-page": "https://mne.tools/mne-features", "Summary": "MNE-Features software for extracting features from multivariate time series", # noqa: E501 diff --git a/tools/install_pre_requirements.sh b/tools/install_pre_requirements.sh index 01bfc086fe5..c717b1b477b 100755 --- a/tools/install_pre_requirements.sh +++ b/tools/install_pre_requirements.sh @@ -27,15 +27,9 @@ python -m pip install $STD_ARGS --only-binary ":all:" --default-timeout=60 \ # statsmodels requires formulaic@main so we need to use --extra-index-url echo "statsmodels" -# https://github.com/statsmodels/statsmodels/issues/9501 -if [ "$PLATFORM" == "Windows" ]; then - STATS_VER=0.14.4 -else - STATS_VER=0.15.0.dev0 -fi python -m pip install $STD_ARGS --only-binary ":all:" \ --extra-index-url "https://pypi.anaconda.org/scientific-python-nightly-wheels/simple" \ - "statsmodels>=$STATS_VER" + "statsmodels>=0.15.0.dev0" # No Numba because it forces an old NumPy version From e4cc4e27106774455347b5d95d8b7b58953af10b Mon Sep 17 00:00:00 2001 From: Clemens Brunner Date: Fri, 7 Feb 2025 23:15:51 +0100 Subject: [PATCH 04/13] Skip first "New Segment" BrainVision marker (#13100) --- doc/changes/devel/13100.bugfix.rst | 1 + mne/io/brainvision/brainvision.py | 308 +++++++++---------- mne/io/brainvision/tests/test_brainvision.py | 6 - 3 files changed, 155 insertions(+), 160 deletions(-) create mode 100644 doc/changes/devel/13100.bugfix.rst diff --git a/doc/changes/devel/13100.bugfix.rst b/doc/changes/devel/13100.bugfix.rst new file mode 100644 index 00000000000..edd4c75264d --- /dev/null +++ b/doc/changes/devel/13100.bugfix.rst @@ -0,0 +1 @@ +Do not convert the first "New Segment" marker in a BrainVision file to an annotation, as it only contains the recording date (which is already available in ``info["meas_date"]``), by `Clemens Brunner`_. \ No newline at end of file diff --git a/mne/io/brainvision/brainvision.py b/mne/io/brainvision/brainvision.py index 07be6ab388a..f9dfffccd55 100644 --- a/mne/io/brainvision/brainvision.py +++ b/mne/io/brainvision/brainvision.py @@ -32,18 +32,18 @@ class RawBrainVision(BaseRaw): ---------- vhdr_fname : path-like Path to the EEG header file. - eog : list or tuple - Names of channels or list of indices that should be designated - EOG channels. Values should correspond to the header file. - Default is ``('HEOGL', 'HEOGR', 'VEOGb')``. - misc : list or tuple of str | ``'auto'`` - Names of channels or list of indices that should be designated - MISC channels. Values should correspond to the electrodes - in the header file. If ``'auto'``, units in header file are used for - inferring misc channels. Default is ``'auto'``. + eog : list of (int | str) | tuple of (int | str) + Names of channels or list of indices that should be designated EOG channels. + Values should correspond to the header file. Default is ``('HEOGL', 'HEOGR', + 'VEOGb')``. + misc : list of (int | str) | tuple of (int | str) | ``'auto'`` + Names of channels or list of indices that should be designated MISC channels. + Values should correspond to the electrodes in the header file. If ``'auto'``, + units in header file are used for inferring misc channels. Default is + ``'auto'``. scale : float - The scaling factor for EEG data. Unless specified otherwise by - header file, units are in microvolts. Default scale factor is 1. + The scaling factor for EEG data. Unless specified otherwise by header file, + units are in microvolts. Default scale factor is 1. ignore_marker_types : bool If ``True``, ignore marker types and only use marker descriptions. Default is ``False``. @@ -64,10 +64,9 @@ class RawBrainVision(BaseRaw): Notes ----- If the BrainVision header file contains impedance measurements, these may be - accessed using ``raw.impedances`` after reading using this function. However, - this attribute will NOT be available after a save and re-load of the data. - That is, it is only available when reading data directly from the BrainVision - header file. + accessed using ``raw.impedances`` after reading using this function. However, this + attribute will NOT be available after a save and re-load of the data. That is, it is + only available when reading data directly from the BrainVision header file. BrainVision markers consist of a type and a description (in addition to other fields like onset and duration). In contrast, annotations in MNE only have a description. @@ -75,6 +74,10 @@ class RawBrainVision(BaseRaw): converted to an annotation "Stimulus/S 1" by default. If you want to ignore the type and instead only use the description, set ``ignore_marker_types=True``, which will convert the same marker to an annotation "S 1". + + The first marker in a BrainVision file is usually a "New Segment" marker, which + contains the recording time. This time is stored in the ``info['meas_date']`` + attribute of the returned object and is not converted to an annotation. """ _extra_attributes = ("impedances",) @@ -110,7 +113,7 @@ def __init__( if isinstance(fmt, dict): # ASCII, this will be slow :( if order == "F": # multiplexed, channels in columns n_skip = 0 - for ii in range(int(fmt["skiplines"])): + for _ in range(int(fmt["skiplines"])): n_skip += len(f.readline()) offsets = np.cumsum([n_skip] + [len(line) for line in f]) n_samples = len(offsets) - 1 @@ -183,9 +186,8 @@ def _read_segment_file(self, data, idx, fi, start, stop, cals, mult): line = fid.readline().decode("ASCII") line = line.strip() - # Not sure why we special-handle the "," character here, - # but let's just keep this for historical and backward- - # compat reasons + # Not sure why we special-handle the "," character here, but let's + # just keep this for historical and backward- compat reasons if ( isinstance(fmt, dict) and "decimalsymbol" in fmt @@ -224,69 +226,66 @@ def _read_segments_c(raw, data, idx, fi, start, stop, cals, mult): _mult_cal_one(data, block, idx, cals, mult) -def _read_mrk(fname, ignore_marker_types=False): +def _read_mrk(fname): """Read annotations from a vmrk/amrk file. Parameters ---------- fname : str vmrk/amrk file to be read. - ignore_marker_types : bool - If True, ignore marker types and only use marker descriptions. Default is False. Returns ------- - onset : array, shape (n_annots,) + onset : list of float The onsets in seconds. - duration : array, shape (n_annots,) + duration : list of float The onsets in seconds. - description : array, shape (n_annots,) - The description of each annotation. + type_ : list of str + The marker types. + description : list of str + The marker descriptions. date_str : str - The recording time as a string. Defaults to empty string if no - recording time is found. + The recording time. Defaults to empty string if no recording time is found. """ # read marker file with open(fname, "rb") as fid: txt = fid.read() - # we don't actually need to know the coding for the header line. - # the characters in it all belong to ASCII and are thus the - # same in Latin-1 and UTF-8 + # we don't actually need to know the encoding for the header line. the characters in + # it all belong to ASCII and are thus the same in Latin-1 and UTF-8 header = txt.decode("ascii", "ignore").split("\n")[0].strip() _check_bv_version(header, "marker") - # although the markers themselves are guaranteed to be ASCII (they - # consist of numbers and a few reserved words), we should still - # decode the file properly here because other (currently unused) - # blocks, such as that the filename are specifying are not - # guaranteed to be ASCII. + # although the markers themselves are guaranteed to be ASCII (they consist of + # numbers and a few reserved words), we should still decode the file properly here + # because other (currently unused) blocks are not guaranteed to be ASCII try: - # if there is an explicit codepage set, use it - # we pretend like it's ascii when searching for the codepage + # if there is an explicit codepage set, use it; we pretend like it's ASCII when + # searching for the codepage cp_setting = re.search( "Codepage=(.+)", txt.decode("ascii", "ignore"), re.IGNORECASE & re.MULTILINE ) codepage = "utf-8" if cp_setting: codepage = cp_setting.group(1).strip() - # BrainAmp Recorder also uses ANSI codepage - # an ANSI codepage raises a LookupError exception - # python recognize ANSI decoding as cp1252 + # BrainAmp Recorder also uses ANSI codepage; an ANSI codepage raises a + # LookupError exception; Python recognize ANSI decoding as cp1252 if codepage == "ANSI": codepage = "cp1252" txt = txt.decode(codepage) except UnicodeDecodeError: - # if UTF-8 (new standard) or explicit codepage setting fails, - # fallback to Latin-1, which is Windows default and implicit - # standard in older recordings + # if UTF-8 (new standard) or explicit codepage setting fails, fallback to + # Latin-1, which is Windows default and implicit standard in older recordings txt = txt.decode("latin-1") # extract Marker Infos block + onset, duration, type_, description = [], [], [], [] + date_str = "" + m = re.search(r"\[Marker Infos\]", txt, re.IGNORECASE) if not m: - return np.array(list()), np.array(list()), np.array(list()), "" + return onset, duration, type_, description, date_str mk_txt = txt[m.end() :] m = re.search(r"^\[.*\]$", mk_txt) @@ -295,29 +294,25 @@ def _read_mrk(fname, ignore_marker_types=False): # extract event information items = re.findall(r"^Mk\d+=(.*)", mk_txt, re.MULTILINE) - onset, duration, description = list(), list(), list() - date_str = "" for info in items: info_data = info.split(",") mtype, mdesc, this_onset, this_duration = info_data[:4] - # commas in mtype and mdesc are handled as "\1". convert back to comma + # commas in mtype and mdesc are handled as "\1", convert back to comma mtype = mtype.replace(r"\1", ",") mdesc = mdesc.replace(r"\1", ",") if date_str == "" and len(info_data) == 5 and mtype == "New Segment": - # to handle the origin of time and handle the presence of multiple - # New Segment annotations. We only keep the first one that is - # different from an empty string for date_str. + # to handle the origin of time and handle the presence of multiple New + # Segment annotations, we only keep the first one that is different from an + # empty string for date_str date_str = info_data[-1] this_duration = int(this_duration) if this_duration.isdigit() else 0 duration.append(this_duration) onset.append(int(this_onset) - 1) # BV is 1-indexed, not 0-indexed - if not ignore_marker_types: - description.append(mtype + "/" + mdesc) - else: - description.append(mdesc) + type_.append(mtype) + description.append(mdesc) - return np.array(onset), np.array(duration), np.array(description), date_str + return onset, duration, type_, description, date_str def _read_annotations_brainvision(fname, sfreq="auto", ignore_marker_types=False): @@ -344,9 +339,7 @@ def _read_annotations_brainvision(fname, sfreq="auto", ignore_marker_types=False annotations : instance of Annotations The annotations present in the file. """ - onset, duration, description, date_str = _read_mrk( - fname, ignore_marker_types=ignore_marker_types - ) + onset, duration, type_, description, date_str = _read_mrk(fname) orig_time = _str_to_meas_date(date_str) if sfreq == "auto": @@ -358,8 +351,17 @@ def _read_annotations_brainvision(fname, sfreq="auto", ignore_marker_types=False _, _, _, info = _aux_hdr_info(hdr_fname) sfreq = info["sfreq"] + # skip the first "New Segment" marker (as it only contains the recording time) + if len(type_) > 0 and type_[0] == "New Segment": + onset = onset[1:] + duration = duration[1:] + type_ = type_[1:] + description = description[1:] + onset = np.array(onset, dtype=float) / sfreq duration = np.array(duration, dtype=float) / sfreq + if not ignore_marker_types: + description = [f"{t}/{d}" for t, d in zip(type_, description)] annotations = Annotations( onset=onset, duration=duration, description=description, orig_time=orig_time ) @@ -372,8 +374,8 @@ def _check_bv_version(header, kind): "MNE-Python currently only supports %s versions 1.0 and 2.0, got unparsable " "%r. Contact MNE-Python developers for support." ) - # optional space, optional Core or V-Amp, optional Exchange, - # Version/Header, optional comma, 1/2 + # optional space, optional Core or V-Amp, optional Exchange, Version/Header, + # optional comma, 1/2 _data_re = r"Brain ?Vision( Core| V-Amp)? Data( Exchange)? %s File,? Version %s\.0" assert kind in ("header", "marker") @@ -416,8 +418,8 @@ def _str_to_meas_date(date_str): if date_str in ["", "0", "00000000000000000000"]: return None - # these calculations are in naive time but should be okay since - # they are relative (subtraction below) + # these calculations are in naive time but should be okay since they are relative + # (subtraction below) try: meas_date = datetime.strptime(date_str, "%Y%m%d%H%M%S%f") except ValueError as e: @@ -436,16 +438,15 @@ def _aux_hdr_info(hdr_fname): # extract the first section to resemble a cfg header = f.readline() codepage = "utf-8" - # we don't actually need to know the coding for the header line. - # the characters in it all belong to ASCII and are thus the - # same in Latin-1 and UTF-8 + # we don't actually need to know the coding for the header line; the characters + # in it all belong to ASCII and are thus the same in Latin-1 and UTF-8 header = header.decode("ascii", "ignore").strip() _check_bv_version(header, "header") settings = f.read() try: # if there is an explicit codepage set, use it - # we pretend like it's ascii when searching for the codepage + # we pretend like it's ASCII when searching for the codepage cp_setting = re.search( "Codepage=(.+)", settings.decode("ascii", "ignore"), @@ -453,16 +454,15 @@ def _aux_hdr_info(hdr_fname): ) if cp_setting: codepage = cp_setting.group(1).strip() - # BrainAmp Recorder also uses ANSI codepage - # an ANSI codepage raises a LookupError exception - # python recognize ANSI decoding as cp1252 + # BrainAmp Recorder also uses ANSI codepage; an ANSI codepage raises a + # LookupError exception; Python recognize ANSI decoding as cp1252 if codepage == "ANSI": codepage = "cp1252" settings = settings.decode(codepage) except UnicodeDecodeError: - # if UTF-8 (new standard) or explicit codepage setting fails, - # fallback to Latin-1, which is Windows default and implicit - # standard in older recordings + # if UTF-8 (new standard) or explicit codepage setting fails, fallback to + # Latin-1, which is Windows default and implicit standard in older + # recordings settings = settings.decode("latin-1") if settings.find("[Comment]") != -1: @@ -499,13 +499,12 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): Names of channels that should be designated EOG channels. Names should correspond to the header file. misc : list or tuple of str | 'auto' - Names of channels or list of indices that should be designated - MISC channels. Values should correspond to the electrodes in the - header file. If 'auto', units in header file are used for inferring - misc channels. Default is ``'auto'``. + Names of channels or list of indices that should be designated MISC channels. + Values should correspond to the electrodes in the header file. If 'auto', units + in header file are used for inferring misc channels. Default is ``'auto'``. scale : float - The scaling factor for EEG data. Unless specified otherwise by - header file, units are in microvolts. Default scale factor is 1. + The scaling factor for EEG data. Unless specified otherwise by header file, + units are in microvolts. Default scale factor is 1. Returns ------- @@ -523,16 +522,16 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): montage : DigMontage Coordinates of the channels, if present in the header file. orig_units : dict - Dictionary mapping channel names to their units as specified in - the header file. Example: {'FC1': 'nV'} + Dictionary mapping channel names to their units as specified in the header file. + Example: {'FC1': 'nV'} """ scale = float(scale) ext = op.splitext(hdr_fname)[-1] ahdr_format = ext == ".ahdr" if ext not in (".vhdr", ".ahdr"): raise OSError( - "The header file must be given to read the data, " - f"not a file with extension '{ext}'." + "The header file must be given to read the data, not a file with extension " + f"'{ext}'." ) settings, cfg, cinfostr, info = _aux_hdr_info(hdr_fname) @@ -552,9 +551,8 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): else: if order == "C": # channels in rows raise NotImplementedError( - "BrainVision files with ASCII data in " - "vectorized order (i.e. channels in rows" - ") are not supported yet." + "BrainVision files with ASCII data in vectorized order (i.e. channels " + "in rows) are not supported yet." ) fmt = {key: cfg.get("ASCII Infos", key) for key in cfg.options("ASCII Infos")} @@ -590,8 +588,8 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): n_samples = cfg.getint(cinfostr, "DataPoints") except configparser.NoOptionError: warn( - "No info on DataPoints found. Inferring number of " - "samples from the data file size." + "No info on DataPoints found. Inferring number of samples from the " + "data file size." ) with open(data_fname, "rb") as fid: fid.seek(0, 2) @@ -644,12 +642,11 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): misc = list(misc_chs.keys()) if misc == "auto" else misc - # create montage: 'Coordinates' section in VHDR/AHDR file corresponds to - # "BVEF" BrainVision Electrode File. The data are based on BrainVision - # Analyzer coordinate system: Defined between standard electrode positions: - # X-axis from T7 to T8, Y-axis from Oz to Fpz, Z-axis orthogonal from - # XY-plane through Cz, fit to a sphere if idealized (when radius=1), - # specified in mm + # create montage: 'Coordinates' section in VHDR/AHDR file corresponds to "BVEF" + # BrainVision Electrode File. The data are based on BrainVision Analyzer coordinate + # system: Defined between standard electrode positions: X-axis from T7 to T8, Y-axis + # from Oz to Fpz, Z-axis orthogonal from XY-plane through Cz, fit to a sphere if + # idealized (when radius=1), specified in mm montage = None if cfg.has_section("Coordinates"): montage_pos = list() @@ -671,8 +668,8 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): if (pos == 0).all() and ch_name not in list(eog) + misc: to_misc.append(ch_name) montage_pos.append(pos) - # Make a montage, normalizing from BrainVision units "mm" to "m", the - # unit used for montages in MNE + # Make a montage, normalizing from BrainVision units "mm" to "m", the unit used + # for montages in MNE montage_pos = np.array(montage_pos) / 1e3 montage = make_dig_montage( ch_pos=dict(zip(montage_names, montage_pos)), coord_frame="head" @@ -688,8 +685,8 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): if np.isnan(cals).any(): raise RuntimeError("Missing channel units") - # Attempts to extract filtering info from header. If not found, both are - # set to zero. + # Attempts to extract filtering info from header. If not found, both are set to + # zero. settings = settings.splitlines() idx = None @@ -703,9 +700,9 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): else: idx = None - # If software filters are active, then they override the hardware setup - # But we still want to be able to double check the channel names - # for alignment purposes, we keep track of the hardware setting idx + # If software filters are active, then they override the hardware setup; we still + # want to be able to double check the channel names for alignment purposes, we keep + # track of the hardware setting idx idx_amp = idx filter_list_has_ch_name = True @@ -716,8 +713,8 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): hp_col, lp_col = 1, 2 filter_list_has_ch_name = False warn( - "Online software filter detected. Using software " - "filter settings and ignoring hardware values" + "Online software filter detected. Using software filter settings " + "and ignoring hardware values" ) break else: @@ -727,19 +724,18 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): lowpass = [] highpass = [] - # for newer BV files, the unit is specified for every channel - # separated by a single space, while for older files, the unit is - # specified in the column headers + # for newer BV files, the unit is specified for every channel separated by a + # single space, while for older files, the unit is specified in the column + # headers divider = r"\s+" if "Resolution / Unit" in settings[idx]: shift = 1 # shift for unit else: shift = 0 - # Extract filter units and convert from seconds to Hz if necessary. - # this cannot be done as post-processing as the inverse t-f - # relationship means that the min/max comparisons don't make sense - # unless we know the units. + # Extract filter units and convert from seconds to Hz if necessary. this cannot + # be done as post-processing as the inverse t-f relationship means that the + # min/max comparisons don't make sense unless we know the units. # # For reasoning about the s to Hz conversion, see this reference: # `Ebersole, J. S., & Pedley, T. A. (Eds.). (2003). @@ -787,20 +783,20 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): else: heterogeneous_hp_filter = True if hp_s: - # We convert channels with disabled filters to having - # highpass relaxed / no filters + # We convert channels with disabled filters to having highpass relaxed / + # no filters highpass = [ float(filt) if filt not in ("NaN", "Off", "DC") else np.inf for filt in highpass ] info["highpass"] = np.max(np.array(highpass, dtype=np.float64)) - # Conveniently enough 1 / np.inf = 0.0, so this works for - # DC / no highpass filter + # Conveniently enough 1 / np.inf = 0.0, so this works for DC / no + # highpass filter # filter time constant t [secs] to Hz conversion: 1/2*pi*t info["highpass"] = 1.0 / (2 * np.pi * info["highpass"]) - # not exactly the cleanest use of FP, but this makes us - # more conservative in *not* warning. + # not exactly the cleanest use of FP, but this makes us more + # conservative in *not* warning. if info["highpass"] == 0.0 and len(set(highpass)) == 1: # not actually heterogeneous in effect # ... just heterogeneously disabled @@ -818,9 +814,8 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): if heterogeneous_hp_filter: warn( - "Channels contain different highpass filters. " - f"Lowest (weakest) filter setting ({info['highpass']:0.2f} Hz) " - "will be stored." + "Channels contain different highpass filters. Lowest (weakest) " + f"filter setting ({info['highpass']:0.2f} Hz) will be stored." ) if len(lowpass) == 0: @@ -837,8 +832,8 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): else: heterogeneous_lp_filter = True if lp_s: - # We convert channels with disabled filters to having - # infinitely relaxed / no filters + # We convert channels with disabled filters to having infinitely relaxed + # / no filters lowpass = [ float(filt) if filt not in ("NaN", "Off", "0") else 0.0 for filt in lowpass @@ -850,19 +845,19 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): except ZeroDivisionError: if len(set(lowpass)) == 1: - # No lowpass actually set for the weakest setting - # so we set lowpass to the Nyquist frequency + # No lowpass actually set for the weakest setting so we set + # lowpass to the Nyquist frequency info["lowpass"] = info["sfreq"] / 2.0 # not actually heterogeneous in effect # ... just heterogeneously disabled heterogeneous_lp_filter = False else: - # no lowpass filter is the weakest filter, - # but it wasn't the only filter + # no lowpass filter is the weakest filter, but it wasn't the + # only filter pass else: - # We convert channels with disabled filters to having - # infinitely relaxed / no filters + # We convert channels with disabled filters to having infinitely relaxed + # / no filters lowpass = [ float(filt) if filt not in ("NaN", "Off", "0") else np.inf for filt in lowpass @@ -870,8 +865,8 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): info["lowpass"] = np.max(np.array(lowpass, dtype=np.float64)) if np.isinf(info["lowpass"]): - # No lowpass actually set for the weakest setting - # so we set lowpass to the Nyquist frequency + # No lowpass actually set for the weakest setting so we set lowpass + # to the Nyquist frequency info["lowpass"] = info["sfreq"] / 2.0 if len(set(lowpass)) == 1: # not actually heterogeneous in effect @@ -879,10 +874,10 @@ def _get_hdr_info(hdr_fname, eog, misc, scale): heterogeneous_lp_filter = False if heterogeneous_lp_filter: - # this isn't clean FP, but then again, we only want to provide - # the Nyquist hint when the lowpass filter was actually - # calculated from dividing the sampling frequency by 2, so the - # exact/direct comparison (instead of tolerance) makes sense + # this isn't clean FP, but then again, we only want to provide the + # Nyquist hint when the lowpass filter was actually calculated from + # dividing the sampling frequency by 2, so the exact/direct comparison + # (instead of tolerance) makes sense if info["lowpass"] == info["sfreq"] / 2.0: nyquist = ", Nyquist limit" else: @@ -953,29 +948,31 @@ def read_raw_brainvision( ---------- vhdr_fname : path-like Path to the EEG header file. - eog : list or tuple of str - Names of channels or list of indices that should be designated - EOG channels. Values should correspond to the header file - Default is ``('HEOGL', 'HEOGR', 'VEOGb')``. - misc : list or tuple of str | ``'auto'`` - Names of channels or list of indices that should be designated - MISC channels. Values should correspond to the electrodes in the - header file. If ``'auto'``, units in header file are used for inferring - misc channels. Default is ``'auto'``. + eog : list of (int | str) | tuple of (int | str) + Names of channels or list of indices that should be designated EOG channels. + Values should correspond to the header file Default is ``('HEOGL', 'HEOGR', + 'VEOGb')``. + misc : list of (int | str) | tuple of (int | str) | ``'auto'`` + Names of channels or list of indices that should be designated MISC channels. + Values should correspond to the electrodes in the header file. If ``'auto'``, + units in header file are used for inferring misc channels. Default is + ``'auto'``. scale : float - The scaling factor for EEG data. Unless specified otherwise by - header file, units are in microvolts. Default scale factor is 1. + The scaling factor for EEG data. Unless specified otherwise by header file, + units are in microvolts. Default scale factor is 1. ignore_marker_types : bool If ``True``, ignore marker types and only use marker descriptions. Default is ``False``. + + .. versionadded:: 1.8 %(preload)s %(verbose)s Returns ------- raw : instance of RawBrainVision - A Raw object containing BrainVision data. - See :class:`mne.io.Raw` for documentation of attributes and methods. + A Raw object containing BrainVision data. See :class:`mne.io.Raw` for + documentation of attributes and methods. See Also -------- @@ -984,10 +981,9 @@ def read_raw_brainvision( Notes ----- If the BrainVision header file contains impedance measurements, these may be - accessed using ``raw.impedances`` after reading using this function. However, - this attribute will NOT be available after a save and re-load of the data. - That is, it is only available when reading data directly from the BrainVision - header file. + accessed using ``raw.impedances`` after reading using this function. However, this + attribute will NOT be available after a save and re-load of the data. That is, it is + only available when reading data directly from the BrainVision header file. BrainVision markers consist of a type and a description (in addition to other fields like onset and duration). In contrast, annotations in MNE only have a description. @@ -995,6 +991,10 @@ def read_raw_brainvision( converted to an annotation "Stimulus/S 1" by default. If you want to ignore the type and instead only use the description, set ``ignore_marker_types=True``, which will convert the same marker to an annotation "S 1". + + The first marker in a BrainVision file is usually a "New Segment" marker, which + contains the recording time. This time is stored in the ``info['meas_date']`` + attribute of the returned object and is not converted to an annotation. """ return RawBrainVision( vhdr_fname=vhdr_fname, @@ -1070,8 +1070,8 @@ def _parse_impedance(settings, recording_date=None): impedance_unit = impedance_setting[1].lstrip("[").rstrip("]") impedance_time = None - # If we have a recording date, we can update it with the time of - # impedance measurement + # If we have a recording date, we can update it with the time of impedance + # measurement if recording_date is not None: meas_time = [int(i) for i in impedance_setting[3].split(":")] impedance_time = recording_date.replace( @@ -1081,8 +1081,8 @@ def _parse_impedance(settings, recording_date=None): microsecond=0, ) for setting in settings[idx + 1 :]: - # Parse channel impedances until we find a line that doesn't start - # with a channel name and optional +/- polarity for passive elecs + # Parse channel impedances until we find a line that doesn't start with a + # channel name and optional +/- polarity for passive elecs match = re.match(r"[ a-zA-Z0-9_+-]+:", setting) if match: channel_name = match.group().rstrip(":") diff --git a/mne/io/brainvision/tests/test_brainvision.py b/mne/io/brainvision/tests/test_brainvision.py index 704af064b48..8366f15dc3a 100644 --- a/mne/io/brainvision/tests/test_brainvision.py +++ b/mne/io/brainvision/tests/test_brainvision.py @@ -665,7 +665,6 @@ def test_ignore_marker_types(): # default behavior (do not ignore marker types) raw = read_raw_brainvision(vhdr_path) expected_descriptions = [ - "New Segment/", "Stimulus/S253", "Stimulus/S255", "Event/254", @@ -685,7 +684,6 @@ def test_ignore_marker_types(): # ignore marker types raw = read_raw_brainvision(vhdr_path, ignore_marker_types=True) expected_descriptions = [ - "", "S253", "S255", "254", @@ -720,7 +718,6 @@ def test_read_vhdr_annotations_and_events(tmp_path): expected_orig_time = _stamp_to_dt((1384359243, 794232)) expected_onset_latency = np.array( [ - 0, 486.0, 496.0, 1769.0, @@ -738,7 +735,6 @@ def test_read_vhdr_annotations_and_events(tmp_path): ] ) expected_annot_description = [ - "New Segment/", "Stimulus/S253", "Stimulus/S255", "Event/254", @@ -760,7 +756,6 @@ def test_read_vhdr_annotations_and_events(tmp_path): expected_onset_latency, np.zeros_like(expected_onset_latency), [ - 99999, 253, 255, 254, @@ -782,7 +777,6 @@ def test_read_vhdr_annotations_and_events(tmp_path): .T ) expected_event_id = { - "New Segment/": 99999, "Stimulus/S253": 253, "Stimulus/S255": 255, "Event/254": 254, From 8eaa521270c17096586e45211f4df8a9f6bb57be Mon Sep 17 00:00:00 2001 From: "Thomas S. Binns" Date: Mon, 10 Feb 2025 15:12:39 +0000 Subject: [PATCH 05/13] Add support for n-dimensional arrays in `_tfr_from_mt` (#13104) Co-authored-by: Eric Larson --- mne/time_frequency/tfr.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/mne/time_frequency/tfr.py b/mne/time_frequency/tfr.py index f4a01e87895..0c8bb0f4fb0 100644 --- a/mne/time_frequency/tfr.py +++ b/mne/time_frequency/tfr.py @@ -4291,19 +4291,20 @@ def _tfr_from_mt(x_mt, weights): Parameters ---------- - x_mt : array, shape (n_channels, n_tapers, n_freqs, n_times) + x_mt : array, shape (..., n_tapers, n_freqs, n_times) The complex-valued multitaper coefficients. weights : array, shape (n_tapers, n_freqs) The weights to use to combine the tapered estimates. Returns ------- - tfr : array, shape (n_channels, n_freqs, n_times) + tfr : array, shape (..., n_freqs, n_times) The time-frequency power estimates. """ - weights = weights[np.newaxis, :, :, np.newaxis] # add singleton channel & time dims + # add singleton dim for time and any dims preceding the tapers + weights = weights[..., np.newaxis] tfr = weights * x_mt tfr *= tfr.conj() - tfr = tfr.real.sum(axis=1) - tfr *= 2 / (weights * weights.conj()).real.sum(axis=1) + tfr = tfr.real.sum(axis=-3) + tfr *= 2 / (weights * weights.conj()).real.sum(axis=-3) return tfr From 64ed25561841e6d91413564222e498c44244ad94 Mon Sep 17 00:00:00 2001 From: Stefan Appelhoff Date: Mon, 10 Feb 2025 17:38:04 +0100 Subject: [PATCH 06/13] add overwrite and verbose params to info.save (#13107) Co-authored-by: Eric Larson --- doc/changes/devel/13107.newfeature.rst | 1 + mne/_fiff/meas_info.py | 13 +++++++++++-- mne/_fiff/tests/test_meas_info.py | 2 +- 3 files changed, 13 insertions(+), 3 deletions(-) create mode 100644 doc/changes/devel/13107.newfeature.rst diff --git a/doc/changes/devel/13107.newfeature.rst b/doc/changes/devel/13107.newfeature.rst new file mode 100644 index 00000000000..a19381fbdb2 --- /dev/null +++ b/doc/changes/devel/13107.newfeature.rst @@ -0,0 +1 @@ +The :meth:`mne.Info.save` method now has an ``overwrite`` and a ``verbose`` parameter, by `Stefan Appelhoff`_. diff --git a/mne/_fiff/meas_info.py b/mne/_fiff/meas_info.py index 51612824a6a..28f0629c323 100644 --- a/mne/_fiff/meas_info.py +++ b/mne/_fiff/meas_info.py @@ -1935,15 +1935,24 @@ def _repr_html_(self): info_template = _get_html_template("repr", "info.html.jinja") return info_template.render(info=self) - def save(self, fname): + @verbose + def save(self, fname, *, overwrite=False, verbose=None): """Write measurement info in fif file. Parameters ---------- fname : path-like The name of the file. Should end by ``'-info.fif'``. + %(overwrite)s + + .. versionadded:: 1.10 + %(verbose)s + + See Also + -------- + mne.io.write_info """ - write_info(fname, self) + write_info(fname, self, overwrite=overwrite) def _simplify_info(info, *, keep=()): diff --git a/mne/_fiff/tests/test_meas_info.py b/mne/_fiff/tests/test_meas_info.py index a38ecaade50..9c0639a830c 100644 --- a/mne/_fiff/tests/test_meas_info.py +++ b/mne/_fiff/tests/test_meas_info.py @@ -975,7 +975,7 @@ def test_field_round_trip(tmp_path): meas_date=_stamp_to_dt((1, 2)), ) fname = tmp_path / "temp-info.fif" - write_info(fname, info) + info.save(fname) info_read = read_info(fname) assert_object_equal(info, info_read) with pytest.raises(TypeError, match="datetime"): From 7b316cb500a5039a9e79ada02f117f650801d17b Mon Sep 17 00:00:00 2001 From: Antoine Collas Date: Mon, 10 Feb 2025 20:02:10 +0100 Subject: [PATCH 07/13] ENH: add interpolate_to method (#13044) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: autofix-ci[bot] <114827586+autofix-ci[bot]@users.noreply.github.com> Co-authored-by: Eric Larson --- doc/changes/devel/13044.newfeature.rst | 1 + doc/changes/names.inc | 1 + doc/references.bib | 8 ++ examples/preprocessing/interpolate_to.py | 81 ++++++++++++ mne/channels/channels.py | 158 ++++++++++++++++++++++- mne/channels/tests/test_interpolation.py | 84 +++++++++++- 6 files changed, 331 insertions(+), 2 deletions(-) create mode 100644 doc/changes/devel/13044.newfeature.rst create mode 100644 examples/preprocessing/interpolate_to.py diff --git a/doc/changes/devel/13044.newfeature.rst b/doc/changes/devel/13044.newfeature.rst new file mode 100644 index 00000000000..9633aba66b9 --- /dev/null +++ b/doc/changes/devel/13044.newfeature.rst @@ -0,0 +1 @@ +Add :meth:`mne.Evoked.interpolate_to` to allow interpolating EEG data to other montages, by :newcontrib:`Antoine Collas`. diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 18753a0c872..282fa8341a0 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -24,6 +24,7 @@ .. _Anna Padee: https://github.com/apadee/ .. _Annalisa Pascarella: https://www.iac.cnr.it/personale/annalisa-pascarella .. _Anne-Sophie Dubarry: https://github.com/annesodub +.. _Antoine Collas: https://www.antoinecollas.fr .. _Antoine Gauthier: https://github.com/Okamille .. _Antti Rantala: https://github.com/Odingod .. _Apoorva Karekal: https://github.com/apoorva6262 diff --git a/doc/references.bib b/doc/references.bib index e2578ed18f2..f0addb5f3b2 100644 --- a/doc/references.bib +++ b/doc/references.bib @@ -2514,3 +2514,11 @@ @article{OyamaEtAl2015 year = {2015}, pages = {24--36}, } + +@inproceedings{MellotEtAl2024, + title = {Physics-informed and Unsupervised Riemannian Domain Adaptation for Machine Learning on Heterogeneous EEG Datasets}, + author = {Mellot, Apolline and Collas, Antoine and Chevallier, Sylvain and Engemann, Denis and Gramfort, Alexandre}, + booktitle = {Proceedings of the 32nd European Signal Processing Conference (EUSIPCO)}, + year = {2024}, + address = {Lyon, France} +} diff --git a/examples/preprocessing/interpolate_to.py b/examples/preprocessing/interpolate_to.py new file mode 100644 index 00000000000..b97a7251cbb --- /dev/null +++ b/examples/preprocessing/interpolate_to.py @@ -0,0 +1,81 @@ +""" +.. _ex-interpolate-to-any-montage: + +====================================================== +Interpolate EEG data to any montage +====================================================== + +This example demonstrates how to interpolate EEG channels to match a given montage. +This can be useful for standardizing +EEG channel layouts across different datasets (see :footcite:`MellotEtAl2024`). + +- Using the field interpolation for EEG data. +- Using the target montage "biosemi16". + +In this example, the data from the original EEG channels will be +interpolated onto the positions defined by the "biosemi16" montage. +""" + +# Authors: Antoine Collas +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +import matplotlib.pyplot as plt + +import mne +from mne.channels import make_standard_montage +from mne.datasets import sample + +print(__doc__) +ylim = (-10, 10) + +# %% +# Load EEG data +data_path = sample.data_path() +eeg_file_path = data_path / "MEG" / "sample" / "sample_audvis-ave.fif" +evoked = mne.read_evokeds(eeg_file_path, condition="Left Auditory", baseline=(None, 0)) + +# Select only EEG channels +evoked.pick("eeg") + +# Plot the original EEG layout +evoked.plot(exclude=[], picks="eeg", ylim=dict(eeg=ylim)) + +# %% +# Define the target montage +standard_montage = make_standard_montage("biosemi16") + +# %% +# Use interpolate_to to project EEG data to the standard montage +evoked_interpolated_spline = evoked.copy().interpolate_to( + standard_montage, method="spline" +) + +# Plot the interpolated EEG layout +evoked_interpolated_spline.plot(exclude=[], picks="eeg", ylim=dict(eeg=ylim)) + +# %% +# Use interpolate_to to project EEG data to the standard montage +evoked_interpolated_mne = evoked.copy().interpolate_to(standard_montage, method="MNE") + +# Plot the interpolated EEG layout +evoked_interpolated_mne.plot(exclude=[], picks="eeg", ylim=dict(eeg=ylim)) + +# %% +# Comparing before and after interpolation +fig, axs = plt.subplots(3, 1, figsize=(8, 6), constrained_layout=True) +evoked.plot(exclude=[], picks="eeg", axes=axs[0], show=False, ylim=dict(eeg=ylim)) +axs[0].set_title("Original EEG Layout") +evoked_interpolated_spline.plot( + exclude=[], picks="eeg", axes=axs[1], show=False, ylim=dict(eeg=ylim) +) +axs[1].set_title("Interpolated to Standard 1020 Montage using spline interpolation") +evoked_interpolated_mne.plot( + exclude=[], picks="eeg", axes=axs[2], show=False, ylim=dict(eeg=ylim) +) +axs[2].set_title("Interpolated to Standard 1020 Montage using MNE interpolation") + +# %% +# References +# ---------- +# .. footbibliography:: diff --git a/mne/channels/channels.py b/mne/channels/channels.py index bf9e58f2819..d0e57eecb5f 100644 --- a/mne/channels/channels.py +++ b/mne/channels/channels.py @@ -41,7 +41,7 @@ pick_info, pick_types, ) -from .._fiff.proj import setup_proj +from .._fiff.proj import _has_eeg_average_ref_proj, setup_proj from .._fiff.reference import add_reference_channels, set_eeg_reference from .._fiff.tag import _rename_list from ..bem import _check_origin @@ -960,6 +960,162 @@ def interpolate_bads( return self + def interpolate_to(self, sensors, origin="auto", method="spline", reg=0.0): + """Interpolate EEG data onto a new montage. + + .. warning:: + Be careful, only EEG channels are interpolated. Other channel types are + not interpolated. + + Parameters + ---------- + sensors : DigMontage + The target montage containing channel positions to interpolate onto. + origin : array-like, shape (3,) | str + Origin of the sphere in the head coordinate frame and in meters. + Can be ``'auto'`` (default), which means a head-digitization-based + origin fit. + method : str + Method to use for EEG channels. + Supported methods are 'spline' (default) and 'MNE'. + reg : float + The regularization parameter for the interpolation method + (only used when the method is 'spline'). + + Returns + ------- + inst : instance of Raw, Epochs, or Evoked + The instance with updated channel locations and data. + + Notes + ----- + This method is useful for standardizing EEG layouts across datasets. + However, some attributes may be lost after interpolation. + + .. versionadded:: 1.10.0 + """ + from ..epochs import BaseEpochs, EpochsArray + from ..evoked import Evoked, EvokedArray + from ..forward._field_interpolation import _map_meg_or_eeg_channels + from ..io import RawArray + from ..io.base import BaseRaw + from .interpolation import _make_interpolation_matrix + from .montage import DigMontage + + # Check that the method option is valid. + _check_option("method", method, ["spline", "MNE"]) + _validate_type(sensors, DigMontage, "sensors") + + # Get target positions from the montage + ch_pos = sensors.get_positions().get("ch_pos", {}) + target_ch_names = list(ch_pos.keys()) + if not target_ch_names: + raise ValueError( + "The provided sensors configuration has no channel positions." + ) + + # Get original channel order + orig_names = self.info["ch_names"] + + # Identify EEG channel + picks_good_eeg = pick_types(self.info, meg=False, eeg=True, exclude="bads") + if len(picks_good_eeg) == 0: + raise ValueError("No good EEG channels available for interpolation.") + # Also get the full list of EEG channel indices (including bad channels) + picks_remove_eeg = pick_types(self.info, meg=False, eeg=True, exclude=[]) + eeg_names_orig = [orig_names[i] for i in picks_remove_eeg] + + # Identify non-EEG channels in original order + non_eeg_names_ordered = [ch for ch in orig_names if ch not in eeg_names_orig] + + # Create destination info for new EEG channels + sfreq = self.info["sfreq"] + info_interp = create_info( + ch_names=target_ch_names, + sfreq=sfreq, + ch_types=["eeg"] * len(target_ch_names), + ) + info_interp.set_montage(sensors) + info_interp["bads"] = [ch for ch in self.info["bads"] if ch in target_ch_names] + # Do not assign "projs" directly. + + # Compute the interpolation mapping + if method == "spline": + origin_val = _check_origin(origin, self.info) + pos_from = self.info._get_channel_positions(picks_good_eeg) - origin_val + pos_to = np.stack(list(ch_pos.values()), axis=0) + + def _check_pos_sphere(pos): + d = np.linalg.norm(pos, axis=-1) + d_norm = np.mean(d / np.mean(d)) + if np.abs(1.0 - d_norm) > 0.1: + warn("Your spherical fit is poor; interpolation may be inaccurate.") + + _check_pos_sphere(pos_from) + _check_pos_sphere(pos_to) + mapping = _make_interpolation_matrix(pos_from, pos_to, alpha=reg) + + else: + assert method == "MNE" + info_eeg = pick_info(self.info, picks_good_eeg) + # If the original info has an average EEG reference projector but + # the destination info does not, + # update info_interp via a temporary RawArray. + if _has_eeg_average_ref_proj(self.info) and not _has_eeg_average_ref_proj( + info_interp + ): + # Create dummy data: shape (n_channels, 1) + temp_data = np.zeros((len(info_interp["ch_names"]), 1)) + temp_raw = RawArray(temp_data, info_interp, first_samp=0) + # Using the public API, add an average reference projector. + temp_raw.set_eeg_reference( + ref_channels="average", projection=True, verbose=False + ) + # Extract the updated info. + info_interp = temp_raw.info + mapping = _map_meg_or_eeg_channels( + info_eeg, info_interp, mode="accurate", origin=origin + ) + + # Interpolate EEG data + data_good = self.get_data(picks=picks_good_eeg) + data_interp = mapping @ data_good + + # Create a new instance for the interpolated EEG channels + # TODO: Creating a new instance leads to a loss of information. + # We should consider updating the existing instance in the future + # by 1) drop channels, 2) add channels, 3) re-order channels. + if isinstance(self, BaseRaw): + inst_interp = RawArray(data_interp, info_interp, first_samp=self.first_samp) + elif isinstance(self, BaseEpochs): + inst_interp = EpochsArray(data_interp, info_interp) + else: + assert isinstance(self, Evoked) + inst_interp = EvokedArray(data_interp, info_interp) + + # Merge only if non-EEG channels exist + if not non_eeg_names_ordered: + return inst_interp + + inst_non_eeg = self.copy().pick(non_eeg_names_ordered).load_data() + inst_out = inst_non_eeg.add_channels([inst_interp], force_update_info=True) + + # Reorder channels + # Insert the entire new EEG block at the position of the first EEG channel. + orig_names_arr = np.array(orig_names) + mask_eeg = np.isin(orig_names_arr, eeg_names_orig) + if mask_eeg.any(): + first_eeg_index = np.where(mask_eeg)[0][0] + pre = orig_names_arr[:first_eeg_index] + new_eeg = np.array(info_interp["ch_names"]) + post = orig_names_arr[first_eeg_index:] + post = post[~np.isin(orig_names_arr[first_eeg_index:], eeg_names_orig)] + new_order = np.concatenate((pre, new_eeg, post)).tolist() + else: + new_order = orig_names + inst_out.reorder_channels(new_order) + return inst_out + @verbose def rename_channels(info, mapping, allow_duplicates=False, *, verbose=None): diff --git a/mne/channels/tests/test_interpolation.py b/mne/channels/tests/test_interpolation.py index a881a41edcc..62c7d79e3eb 100644 --- a/mne/channels/tests/test_interpolation.py +++ b/mne/channels/tests/test_interpolation.py @@ -12,7 +12,7 @@ from mne import Epochs, pick_channels, pick_types, read_events from mne._fiff.constants import FIFF from mne._fiff.proj import _has_eeg_average_ref_proj -from mne.channels import make_dig_montage +from mne.channels import make_dig_montage, make_standard_montage from mne.channels.interpolation import _make_interpolation_matrix from mne.datasets import testing from mne.io import RawArray, read_raw_ctf, read_raw_fif, read_raw_nirx @@ -439,3 +439,85 @@ def test_method_str(): raw.interpolate_bads(method="spline") raw.pick("eeg", exclude=()) raw.interpolate_bads(method="spline") + + +@pytest.mark.parametrize("montage_name", ["biosemi16", "standard_1020"]) +@pytest.mark.parametrize("method", ["spline", "MNE"]) +@pytest.mark.parametrize("data_type", ["raw", "epochs", "evoked"]) +def test_interpolate_to_eeg(montage_name, method, data_type): + """Test the interpolate_to method for EEG for raw, epochs, and evoked.""" + # Load EEG data + raw, epochs_eeg = _load_data("eeg") + epochs_eeg = epochs_eeg.copy() + + # Load data for raw + raw.load_data() + + # Create a target montage + montage = make_standard_montage(montage_name) + + # Prepare data to interpolate to + if data_type == "raw": + inst = raw.copy() + elif data_type == "epochs": + inst = epochs_eeg.copy() + elif data_type == "evoked": + inst = epochs_eeg.average() + shape = list(inst._data.shape) + orig_total = len(inst.info["ch_names"]) + n_eeg_orig = len(pick_types(inst.info, eeg=True)) + + # Assert first and last channels are not EEG + if data_type == "raw": + ch_types = inst.get_channel_types() + assert ch_types[0] != "eeg" + assert ch_types[-1] != "eeg" + + # Record the names and data of the first and last channels. + if data_type == "raw": + first_name = inst.info["ch_names"][0] + last_name = inst.info["ch_names"][-1] + data_first = inst._data[..., 0, :].copy() + data_last = inst._data[..., -1, :].copy() + + # Interpolate the EEG channels. + inst_interp = inst.copy().interpolate_to(montage, method=method) + + # Check that the new channel names include the montage channels. + assert set(montage.ch_names).issubset(set(inst_interp.info["ch_names"])) + # Check that the overall channel order is changed. + assert inst.info["ch_names"] != inst_interp.info["ch_names"] + + # Check that the data shape is as expected. + new_nchan_expected = orig_total - n_eeg_orig + len(montage.ch_names) + expected_shape = (new_nchan_expected, shape[-1]) + if len(shape) == 3: + expected_shape = (shape[0],) + expected_shape + assert inst_interp._data.shape == expected_shape + + # Verify that the first and last channels retain their positions. + if data_type == "raw": + assert inst_interp.info["ch_names"][0] == first_name + assert inst_interp.info["ch_names"][-1] == last_name + + # Verify that the data for the first and last channels is unchanged. + if data_type == "raw": + np.testing.assert_allclose( + inst_interp._data[..., 0, :], + data_first, + err_msg="Data for the first non-EEG channel has changed.", + ) + np.testing.assert_allclose( + inst_interp._data[..., -1, :], + data_last, + err_msg="Data for the last non-EEG channel has changed.", + ) + + # Validate that bad channels are carried over. + # Mark the first non eeg channel as bad + all_ch = inst_interp.info["ch_names"] + eeg_ch = [all_ch[i] for i in pick_types(inst_interp.info, eeg=True)] + bads = [ch for ch in all_ch if ch not in eeg_ch][:1] + inst.info["bads"] = bads + inst_interp = inst.copy().interpolate_to(montage, method=method) + assert inst_interp.info["bads"] == bads From cf5ef5fea4412dcc4b8f7c79c14b0faf9342efb1 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 11 Feb 2025 00:00:28 +0000 Subject: [PATCH 08/13] [pre-commit.ci] pre-commit autoupdate (#13110) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- .pre-commit-config.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 5abf98c9945..41210aac357 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,7 +1,7 @@ repos: # Ruff mne - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.9.4 + rev: v0.9.6 hooks: - id: ruff name: ruff lint mne @@ -82,7 +82,7 @@ repos: # zizmor - repo: https://github.com/woodruffw/zizmor-pre-commit - rev: v1.3.0 + rev: v1.3.1 hooks: - id: zizmor From 77195a2a5ff42fa0e193f889e712b279b6c65fd6 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 12 Feb 2025 09:10:27 -0500 Subject: [PATCH 09/13] BUG: Fix bug with reading dig strings (#13083) Co-authored-by: autofix-ci[bot] <114827586+autofix-ci[bot]@users.noreply.github.com> --- doc/changes/devel/13083.bugfix.rst | 1 + mne/_fiff/_digitization.py | 3 +++ mne/_fiff/meas_info.py | 26 ++++---------------------- mne/_fiff/open.py | 2 +- mne/_fiff/tag.py | 24 +++++++++++++++++------- mne/_fiff/tests/test_meas_info.py | 11 ++++++++++- 6 files changed, 36 insertions(+), 31 deletions(-) create mode 100644 doc/changes/devel/13083.bugfix.rst diff --git a/doc/changes/devel/13083.bugfix.rst b/doc/changes/devel/13083.bugfix.rst new file mode 100644 index 00000000000..20ffba40936 --- /dev/null +++ b/doc/changes/devel/13083.bugfix.rst @@ -0,0 +1 @@ +Fix bug with reading digitization points from digitization strings with newer MEGIN systems, by `Eric Larson`_. \ No newline at end of file diff --git a/mne/_fiff/_digitization.py b/mne/_fiff/_digitization.py index eb8b6bc396a..9e74dcf1c96 100644 --- a/mne/_fiff/_digitization.py +++ b/mne/_fiff/_digitization.py @@ -181,6 +181,9 @@ def _read_dig_fif(fid, meas_info, *, return_ch_names=False): if kind == FIFF.FIFF_DIG_POINT: tag = read_tag(fid, pos) dig.append(tag.data) + elif kind == FIFF.FIFF_DIG_STRING: + tag = read_tag(fid, pos) + dig.extend(tag.data) elif kind == FIFF.FIFF_MNE_COORD_FRAME: tag = read_tag(fid, pos) coord_frame = _coord_frame_named.get(int(tag.data.item())) diff --git a/mne/_fiff/meas_info.py b/mne/_fiff/meas_info.py index 28f0629c323..76d78782ac8 100644 --- a/mne/_fiff/meas_info.py +++ b/mne/_fiff/meas_info.py @@ -46,7 +46,7 @@ write_dig, ) from .compensator import get_current_comp -from .constants import FIFF, _ch_unit_mul_named, _coord_frame_named +from .constants import FIFF, _ch_unit_mul_named from .ctf_comp import _read_ctf_comp, write_ctf_comp from .open import fiff_open from .pick import ( @@ -1970,7 +1970,7 @@ def _simplify_info(info, *, keep=()): @verbose -def read_fiducials(fname, verbose=None): +def read_fiducials(fname, *, verbose=None): """Read fiducials from a fiff file. Parameters @@ -1990,26 +1990,8 @@ def read_fiducials(fname, verbose=None): fname = _check_fname(fname=fname, overwrite="read", must_exist=True) fid, tree, _ = fiff_open(fname) with fid: - isotrak = dir_tree_find(tree, FIFF.FIFFB_ISOTRAK) - isotrak = isotrak[0] - pts = [] - coord_frame = FIFF.FIFFV_COORD_HEAD - for k in range(isotrak["nent"]): - kind = isotrak["directory"][k].kind - pos = isotrak["directory"][k].pos - if kind == FIFF.FIFF_DIG_POINT: - tag = read_tag(fid, pos) - pts.append(DigPoint(tag.data)) - elif kind == FIFF.FIFF_MNE_COORD_FRAME: - tag = read_tag(fid, pos) - coord_frame = tag.data[0] - coord_frame = _coord_frame_named.get(coord_frame, coord_frame) - - # coord_frame is not stored in the tag - for pt in pts: - pt["coord_frame"] = coord_frame - - return pts, coord_frame + pts = _read_dig_fif(fid, tree) + return pts, pts[0]["coord_frame"] @verbose diff --git a/mne/_fiff/open.py b/mne/_fiff/open.py index cf50cf96353..1d99bd8ddc2 100644 --- a/mne/_fiff/open.py +++ b/mne/_fiff/open.py @@ -348,7 +348,7 @@ def _show_tree( postpend += " ... list len=" + str(len(tag.data)) elif issparse(tag.data): postpend += ( - f" ... sparse ({tag.data.getformat()}) shape=" + f" ... sparse ({tag.data.__class__.__name__}) shape=" f"{tag.data.shape}" ) else: diff --git a/mne/_fiff/tag.py b/mne/_fiff/tag.py index 3fd36454d58..6a4636ef264 100644 --- a/mne/_fiff/tag.py +++ b/mne/_fiff/tag.py @@ -270,19 +270,26 @@ def _read_id_struct(fid, tag, shape, rlims): ) -def _read_dig_point_struct(fid, tag, shape, rlims): +def _read_dig_point_struct(fid, tag, shape, rlims, *, string=False): """Read dig point struct tag.""" kind = int(np.frombuffer(fid.read(4), dtype=">i4").item()) kind = _dig_kind_named.get(kind, kind) ident = int(np.frombuffer(fid.read(4), dtype=">i4").item()) if kind == FIFF.FIFFV_POINT_CARDINAL: ident = _dig_cardinal_named.get(ident, ident) - return dict( - kind=kind, - ident=ident, - r=np.frombuffer(fid.read(12), dtype=">f4"), - coord_frame=FIFF.FIFFV_COORD_UNKNOWN, - ) + n = 1 if not string else int(np.frombuffer(fid.read(4), dtype=">i4").item()) + out = [ + dict( + kind=kind, + ident=ident, + r=np.frombuffer(fid.read(12), dtype=">f4"), + coord_frame=FIFF.FIFFV_COORD_UNKNOWN, + ) + for _ in range(n) + ] + if not string: + out = out[0] + return out def _read_coord_trans_struct(fid, tag, shape, rlims): @@ -378,11 +385,13 @@ def _read_julian(fid, tag, shape, rlims): FIFF.FIFFT_COMPLEX_DOUBLE: _read_complex_double, FIFF.FIFFT_ID_STRUCT: _read_id_struct, FIFF.FIFFT_DIG_POINT_STRUCT: _read_dig_point_struct, + FIFF.FIFFT_DIG_STRING_STRUCT: partial(_read_dig_point_struct, string=True), FIFF.FIFFT_COORD_TRANS_STRUCT: _read_coord_trans_struct, FIFF.FIFFT_CH_INFO_STRUCT: _read_ch_info_struct, FIFF.FIFFT_OLD_PACK: _read_old_pack, FIFF.FIFFT_DIR_ENTRY_STRUCT: _read_dir_entry_struct, FIFF.FIFFT_JULIAN: _read_julian, + FIFF.FIFFT_VOID: lambda fid, tag, shape, rlims: None, } _call_dict_names = { FIFF.FIFFT_STRING: "str", @@ -390,6 +399,7 @@ def _read_julian(fid, tag, shape, rlims): FIFF.FIFFT_COMPLEX_DOUBLE: "c16", FIFF.FIFFT_ID_STRUCT: "ids", FIFF.FIFFT_DIG_POINT_STRUCT: "dps", + FIFF.FIFFT_DIG_STRING_STRUCT: "dss", FIFF.FIFFT_COORD_TRANS_STRUCT: "cts", FIFF.FIFFT_CH_INFO_STRUCT: "cis", FIFF.FIFFT_OLD_PACK: "op_", diff --git a/mne/_fiff/tests/test_meas_info.py b/mne/_fiff/tests/test_meas_info.py index 9c0639a830c..9ca83697df0 100644 --- a/mne/_fiff/tests/test_meas_info.py +++ b/mne/_fiff/tests/test_meas_info.py @@ -29,7 +29,7 @@ write_cov, write_forward_solution, ) -from mne._fiff import meas_info +from mne._fiff import meas_info, tag from mne._fiff._digitization import DigPoint, _make_dig_points from mne._fiff.constants import FIFF from mne._fiff.meas_info import ( @@ -1268,3 +1268,12 @@ def test_get_montage(): assert len(raw.info.get_montage().ch_names) == len(ch_names) raw.info["bads"] = [ch_names[0]] assert len(raw.info.get_montage().ch_names) == len(ch_names) + + +def test_tag_consistency(): + """Test that structures for tag reading are consistent.""" + call_set = set(tag._call_dict) + call_names = set(tag._call_dict_names) + assert call_set == call_names, "Mismatch between _call_dict and _call_dict_names" + # TODO: This was inspired by FIFF_DIG_STRING gh-13083, we should ideally add a test + # that those dig points can actually be read in correctly at some point. From b4a7c64e5b36747cfd3621e7a043ddd839db4e96 Mon Sep 17 00:00:00 2001 From: GreasyCat <131315c@gmail.com> Date: Thu, 13 Feb 2025 15:50:13 -0500 Subject: [PATCH 10/13] BUG: Fix GDF Out of Bound Error (#13113) Co-authored-by: Scott Huberty <52462026+scott-huberty@users.noreply.github.com> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- doc/changes/devel/13113.bugfix.rst | 1 + doc/changes/names.inc | 1 + mne/io/edf/edf.py | 4 +--- 3 files changed, 3 insertions(+), 3 deletions(-) create mode 100644 doc/changes/devel/13113.bugfix.rst diff --git a/doc/changes/devel/13113.bugfix.rst b/doc/changes/devel/13113.bugfix.rst new file mode 100644 index 00000000000..78589822574 --- /dev/null +++ b/doc/changes/devel/13113.bugfix.rst @@ -0,0 +1 @@ +Fix bug in :func:`mne.io.read_raw_gdf`, by :newcontrib:`Rongfei Jin`. \ No newline at end of file diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 282fa8341a0..a87adaf3047 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -257,6 +257,7 @@ .. _Romain Derollepot: https://github.com/rderollepot .. _Romain Trachel: https://fr.linkedin.com/in/trachelr .. _Roman Goj: https://romanmne.blogspot.co.uk +.. _Rongfei Jin: https://github.com/greasycat .. _Ross Maddox: https://medicine.umich.edu/dept/khri/ross-maddox-phd .. _Rotem Falach: https://github.com/Falach .. _Roy Eric Wieske: https://github.com/Randomidous diff --git a/mne/io/edf/edf.py b/mne/io/edf/edf.py index 09ac24f753e..763ef4f91eb 100644 --- a/mne/io/edf/edf.py +++ b/mne/io/edf/edf.py @@ -1453,9 +1453,7 @@ def _read_gdf_header(fname, exclude, include=None): n_events = np.fromfile(fid, UINT32, 1)[0] else: ne = np.fromfile(fid, UINT8, 3) - n_events = ne[0] - for i in range(1, len(ne)): - n_events = n_events + int(ne[i]) * 2 ** (i * 8) + n_events = sum(int(ne[i]) << (i * 8) for i in range(len(ne))) event_sr = np.fromfile(fid, FLOAT32, 1)[0] pos = np.fromfile(fid, UINT32, n_events) - 1 # 1-based inds From 4c98f8ae5524fea24968fb1aa2bf2028886d5934 Mon Sep 17 00:00:00 2001 From: Stefan Appelhoff Date: Mon, 17 Feb 2025 16:17:04 +0100 Subject: [PATCH 11/13] DOC: add link to governance people so that it is findable. (#13116) --- doc/development/governance.rst | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/development/governance.rst b/doc/development/governance.rst index c66d1489211..77d01b66c85 100644 --- a/doc/development/governance.rst +++ b/doc/development/governance.rst @@ -46,7 +46,9 @@ Governance Model Leadership Roles ^^^^^^^^^^^^^^^^ -The MNE-Python leadership structure shall consist of: +The MNE-Python leadership structure shall consist of the following groups. +A list of the current members of the respective groups is maintained at the +page :ref:`governance-people`. Maintainer Team --------------- From a79809d018d4e261b39e5eed27402528befc9ac8 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Mon, 24 Feb 2025 14:25:16 -0500 Subject: [PATCH 12/13] MAINT: Fix CIs (#13125) --- azure-pipelines.yml | 2 +- tools/install_pre_requirements.sh | 5 ++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 9c42b3286c1..8be7d507375 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -244,7 +244,7 @@ stages: PYTHONIOENCODING: 'utf-8' AZURE_CI_WINDOWS: 'true' PYTHON_ARCH: 'x64' - timeoutInMinutes: 80 + timeoutInMinutes: 90 strategy: maxParallel: 4 matrix: diff --git a/tools/install_pre_requirements.sh b/tools/install_pre_requirements.sh index c717b1b477b..44fb77e665f 100755 --- a/tools/install_pre_requirements.sh +++ b/tools/install_pre_requirements.sh @@ -23,7 +23,10 @@ python -m pip install $STD_ARGS --only-binary ":all:" --default-timeout=60 \ --index-url "https://pypi.anaconda.org/scientific-python-nightly-wheels/simple" \ "numpy>=2.1.0.dev0" "scikit-learn>=1.6.dev0" "scipy>=1.15.0.dev0" \ "pandas>=3.0.0.dev0" \ - "h5py>=3.12.1" "dipy>=1.10.0.dev0" "pyarrow>=19.0.0.dev0" "tables>=3.10.2.dev0" + "h5py>=3.12.1" "dipy>=1.10.0.dev0" +# TODO: should have above: "pyarrow>=20.0.0.dev0" "tables>=3.10.2.dev0" +# https://github.com/apache/arrow/issues/40216#issuecomment-2678899168 +# https://github.com/PyTables/PyTables/issues/1115 # statsmodels requires formulaic@main so we need to use --extra-index-url echo "statsmodels" From c570cfccee22f9869ae1f40c0ce5e8046433f083 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 24 Feb 2025 21:48:31 +0000 Subject: [PATCH 13/13] [pre-commit.ci] pre-commit autoupdate (#13126) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 41210aac357..32c8098aedb 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,7 +1,7 @@ repos: # Ruff mne - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.9.6 + rev: v0.9.7 hooks: - id: ruff name: ruff lint mne