diff --git a/doc/changes/dev/14210.newfeature.rst b/doc/changes/dev/14210.newfeature.rst new file mode 100644 index 00000000000..a657833cd9c --- /dev/null +++ b/doc/changes/dev/14210.newfeature.rst @@ -0,0 +1 @@ +Allow :class:`mne.Epochs` to hold trials of different duration by passing ``tmin`` and/or ``tmax`` as arrays with one entry per event, with :meth:`mne.Epochs.as_fixed` to obtain a fixed-duration copy spanning their union together with the number of epochs contributing at each time point, and with :meth:`mne.Epochs.plot` browsing them at their own lengths rather than padding, demonstrated in :ref:`tut-variable-duration-epochs`, by `Sina Esmaeili`_. diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 6856ca7e2da..6f48b057789 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -425,6 +425,7 @@ .. _Simon Kornblith: https://simonster.com .. _Simon M. Hofmann: https://github.com/SHEscher .. _Simon-Shlomo Poil: https://github.com/simon-shlomo +.. _Sina Esmaeili: https://github.com/snesmaeili .. _Sondre Foslien: https://github.com/sondrfos .. _Sophie Herbst: https://github.com/SophieHerbst .. _Sourav Singh: https://github.com/souravsingh diff --git a/doc/conf.py b/doc/conf.py index 487e86ab884..c492a91498b 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -779,8 +779,8 @@ def fix_sklearn_inherited_docstrings(app, what, name, obj, options, lines): ( # BaseRaw attributes are documented in Raw "py:obj", "(filename|metadata|proj|times|tmax|tmin|annotations|ch_names" - "|compensation_grade|duration|filenames|first_samp|first_time" - "|last_samp|n_times|proj|times|tmax|tmin)", + "|compensation_grade|duration|durations|filenames|first_samp|first_time" + "|last_samp|n_times|proj|times|tmax|tmin|variable_duration)", ), ] suppress_warnings = [ diff --git a/mne/channels/channels.py b/mne/channels/channels.py index d0e9ae73d6e..7ab98ddd261 100644 --- a/mne/channels/channels.py +++ b/mne/channels/channels.py @@ -636,7 +636,12 @@ def _pick_drop_channels(self, idx, *, verbose=None): else: # All others (Evoked, Epochs, Raw) have chs axis=-2 axis = -2 if hasattr(self, "_data"): # skip non-preloaded Raw - self._data = self._data.take(idx, axis=axis) + if isinstance(self._data, list): + # replacing the contents rather than the attribute keeps `_data` + # an ndarray everywhere else this mixin is used + self._data[:] = [epoch.take(idx, axis=axis) for epoch in self._data] + else: + self._data = self._data.take(idx, axis=axis) else: assert isinstance(self, BaseRaw) and not self.preload diff --git a/mne/epochs.py b/mne/epochs.py index f313a94e2ac..e9b5f0a6931 100644 --- a/mne/epochs.py +++ b/mne/epochs.py @@ -10,7 +10,7 @@ from collections import Counter from collections.abc import Callable, Iterable, Iterator from copy import deepcopy -from functools import partial +from functools import partial, wraps from inspect import getfullargspec from pathlib import Path from typing import TYPE_CHECKING, Literal @@ -97,6 +97,7 @@ _prepare_read_metadata, _prepare_write_metadata, _scale_dataframe_data, + _time_mask, _validate_type, check_fname, check_random_state, @@ -385,6 +386,118 @@ def _handle_event_repeated(events, event_id, event_repeated, selection, drop_log return new_events, event_id, selection, drop_log +def _check_variable_bounds(tmin, tmax, n_events): + """Normalize ``tmin``/``tmax``, which may be given per event.""" + arrays = [np.ndim(bound) > 0 for bound in (tmin, tmax)] + if not any(arrays): + return tmin, tmax, False + + out = [] + for name, bound in (("tmin", tmin), ("tmax", tmax)): + bound = np.atleast_1d(np.asarray(bound, dtype=float)) + if bound.ndim != 1: + raise ValueError(f"{name} must be a scalar or 1D array, got {bound.ndim}D.") + if bound.size == 1: + bound = np.repeat(bound, n_events) + elif bound.size != n_events: + raise ValueError( + f"{name} has {bound.size} entries but there are {n_events} " + "events; per-event bounds must match the number of events." + ) + if not np.all(np.isfinite(bound)): + raise ValueError(f"{name} must be finite.") + out.append(bound) + + tmin, tmax = out + if n_events > 0 and np.all(tmin == tmin[0]) and np.all(tmax == tmax[0]): + # arrays were passed, but every epoch has the same bounds, so this is + # the ordinary fixed-duration case and should behave identically + return float(tmin[0]), float(tmax[0]), False + if np.any(tmin > tmax): + bad = int(np.argmax(tmin > tmax)) + raise ValueError( + f"tmin has to be less than or equal to tmax, but epoch {bad} has " + f"tmin={tmin[bad]:g} > tmax={tmax[bad]:g}." + ) + return tmin, tmax, True + + +def _check_variable_data(data, tmin, tmax, events, sfreq): + """Validate per-epoch data for variable-duration epochs.""" + if isinstance(data, np.ndarray) and data.ndim == 3: + data = list(data) + if not isinstance(data, list | tuple): + raise TypeError( + "With per-event tmin/tmax, data must be a list of arrays of shape " + f"(n_channels, n_times_i), got {type(data).__name__}." + ) + if len(data) != len(events): + raise ValueError( + f"data has {len(data)} epochs but there are {len(events)} events." + ) + n_channels = None + for ii, epoch in enumerate(data): + epoch = np.asarray(epoch) + if epoch.ndim != 2: + raise ValueError( + f"Epoch {ii} has ndim={epoch.ndim}, expected 2 (n_channels, n_times)." + ) + if n_channels is None: + n_channels = epoch.shape[0] + elif epoch.shape[0] != n_channels: + raise ValueError( + f"Epoch {ii} has {epoch.shape[0]} channels but epoch 0 has " + f"{n_channels}. Only the time axis may vary between epochs." + ) + want = round((tmax[ii] - tmin[ii]) * sfreq) + 1 + if epoch.shape[1] != want: + raise ValueError( + f"Epoch {ii} has {epoch.shape[1]} samples but tmin={tmin[ii]:g} " + f"and tmax={tmax[ii]:g} at {sfreq:g} Hz imply {want}." + ) + + +def _check_variable_unsupported(*, baseline, reject_tmin, reject_tmax, decim, preload): + """Reject options this implementation does not yet handle.""" + if baseline is not None: + raise NotImplementedError( + "Baseline correction is not implemented for variable-duration " + "epochs; pass baseline=None and correct afterwards. The baseline " + "window would have to be checked against each epoch separately." + ) + for name, value in (("reject_tmin", reject_tmin), ("reject_tmax", reject_tmax)): + if value is not None: + raise NotImplementedError( + f"{name} is not implemented for variable-duration epochs, " + "because the window is not guaranteed to exist in every epoch." + ) + if decim != 1: + raise NotImplementedError( + "decim is not implemented for variable-duration epochs." + ) + if not preload: + raise NotImplementedError( + "Variable-duration epochs must be preloaded. Reading them lazily " + "from disk requires per-epoch bounds in the I/O layer, which is " + "not part of this change." + ) + + +def _is_variable_duration_data(data): + """Whether ``data`` is a sequence of arrays with differing lengths.""" + if isinstance(data, np.ndarray): + return False + if not isinstance(data, list | tuple) or len(data) == 0: + return False + lengths = set() + for epoch in data: + epoch = np.asarray(epoch) + if epoch.ndim != 2: + return False + lengths.add(epoch.shape[1]) + return len(lengths) > 1 + + @fill_doc class BaseEpochs( ProjMixin, @@ -410,9 +523,11 @@ class BaseEpochs( Parameters ---------- %(info_not_none)s - data : ndarray | None + data : ndarray | list of ndarray | None If ``None``, data will be read from the Raw object. If ndarray, must be - of shape (n_epochs, n_channels, n_times). + of shape (n_epochs, n_channels, n_times). A list of + (n_channels, n_times) arrays, one per epoch, gives epochs of differing + duration. %(events_epochs)s %(event_id)s %(epochs_tmin_tmax)s @@ -459,15 +574,22 @@ class BaseEpochs( used as a constructor for Epochs objects (use instead :class:`mne.Epochs`). """ + # Variable-duration epochs keep one array per epoch in this same slot. The + # declaration describes the rectangular case that every reshaping code path + # below relies on; those paths refuse before they run when durations vary. + _data: np.ndarray | None + _tmin_per_epoch: np.ndarray + _tmax_per_epoch: np.ndarray + @verbose def __init__( self, info: Info, - data: np.ndarray | None, + data: np.ndarray | list[np.ndarray] | None, events: np.ndarray, event_id: int | list[int] | dict | str | list[str] | None = None, - tmin: float = -0.2, - tmax: float = 0.5, + tmin: float | np.ndarray = -0.2, + tmax: float | np.ndarray = 0.5, baseline: tuple[float | None, float | None] | None = (None, 0), raw: "BaseRaw | list | None" = None, picks: str | np.ndarray | slice | None = None, @@ -597,6 +719,13 @@ def __init__( self.metadata = metadata # do not set self.events here, let subclass do it + # Variable-duration epochs: tmin and/or tmax may be given per event. + tmin, tmax, self._variable_duration = _check_variable_bounds( + tmin, + tmax, + len(self.events) if getattr(self, "events", None) is not None else 0, + ) + if (detrend not in [None, 0, 1]) or isinstance(detrend, bool): raise ValueError("detrend must be None, 0, or 1") self.detrend = detrend @@ -614,8 +743,16 @@ def __init__( self.preload = False self._data = None self._do_baseline = True + elif self._variable_duration: + _check_variable_data(data, tmin, tmax, self.events, self.info["sfreq"]) + self.preload = True + self._data = [ # ty: ignore[invalid-assignment] # ragged payload + np.asarray(epoch, dtype=np.float64) for epoch in data + ] + self._do_baseline = False else: assert decim == 1 + assert not isinstance(data, list) # only the branch above takes a list if ( data.ndim != 3 or data.shape[2] != round((tmax - tmin) * self.info["sfreq"]) + 1 @@ -630,13 +767,35 @@ def __init__( self._do_baseline = False self._offset = None - if tmin > tmax: + if not self._variable_duration and tmin > tmax: raise ValueError("tmin has to be less than or equal to tmax") # Handle times sfreq = float(self.info["sfreq"]) - start_idx = int(round(tmin * sfreq)) - self._raw_times = np.arange(start_idx, int(round(tmax * sfreq)) + 1) / sfreq + if self._variable_duration: + # Per-epoch bounds are kept, and `times` spans their union so that + # it agrees with `as_fixed()`. Per-epoch axes are available from + # `get_times()`; see `durations`. + self._tmin_per_epoch = tmin + self._tmax_per_epoch = tmax + start_idx = int(round(tmin.min() * sfreq)) + stop_idx = int(round(tmax.max() * sfreq)) + _check_variable_unsupported( + baseline=baseline, + reject_tmin=reject_tmin, + reject_tmax=reject_tmax, + decim=decim, + preload=self.preload or preload_at_end, + ) + baseline = None + else: + # only the variable-duration path carries per-epoch bounds; every + # reader of these guards on _variable_duration first + self._tmin_per_epoch = None # ty: ignore[invalid-assignment] + self._tmax_per_epoch = None # ty: ignore[invalid-assignment] + start_idx = int(round(tmin * sfreq)) + stop_idx = int(round(tmax * sfreq)) + self._raw_times = np.arange(start_idx, stop_idx + 1) / sfreq self._set_times(self._raw_times) # check reject_tmin and reject_tmax @@ -670,12 +829,22 @@ def __init__( # decimation self._decim = 1 - self.decimate(decim) + if not self._variable_duration: + self.decimate(decim) + else: + # decim != 1 is rejected above, and decimate() is wrapped to raise + # for these, so set the slice it would otherwise have set + self._decim_slice = slice(None, None, None) # baseline correction: replace `None` tuple elements with actual times - self.baseline = _check_baseline( - baseline, times=self.times, sfreq=self.info["sfreq"] - ) + if self._variable_duration: + # baseline is forced to None above, and there is no shared time axis + # to resolve None endpoints against + self.baseline = None + else: + self.baseline = _check_baseline( + baseline, times=self.times, sfreq=self.info["sfreq"] + ) if self.baseline is not None and self.baseline != baseline: logger.info( f"Setting baseline interval to " @@ -743,7 +912,9 @@ def _check_consistency(self): assert len(self.drop_log) >= len(self.events) assert len(self.selection) == sum(len(dl) == 0 for dl in self.drop_log) assert hasattr(self, "_times_readonly") - assert not self.times.flags["WRITEABLE"] + # `times` is deliberately unavailable when durations vary, so check the + # underlying vector rather than going through the property + assert not self._times_readonly.flags["WRITEABLE"] assert isinstance(self.drop_log, tuple) assert all(isinstance(log, tuple) for log in self.drop_log) assert all(isinstance(s, str) for log in self.drop_log for s in log) @@ -773,6 +944,151 @@ def reset_index(self) -> None: self.drop_log = (tuple(),) * len(self.events) self._check_consistency() + # -- variable-duration epochs --------------------------------------- + @property + def variable_duration(self): + """Whether the epochs have per-event ``tmin``/``tmax``.""" + return self._variable_duration + + @property + def times(self): + """The time axis shared by all epochs, in seconds. + + Raises + ------ + RuntimeError + When durations vary, because there is no such axis. Returning the + union would leave ``len(epochs.times) == data.shape[-1]`` false while + looking ordinary, which is how silent errors happen downstream. Use + :meth:`get_times` for one epoch, or :meth:`as_fixed` for the padded + common axis. + """ + if self._variable_duration: + raise RuntimeError( + f"These {len(self.events)} epochs have durations from " + f"{self.durations.min():.3f} to {self.durations.max():.3f} s, so " + "there is no time axis they share. Use get_times(epoch) for one " + "epoch's axis, durations for their lengths, or as_fixed() for " + "the padded common axis." + ) + return super().times + + @property + def tmin(self): + """First time point, or one per epoch when durations vary.""" + if self._variable_duration: + return self._tmin_per_epoch + return self.times[0] + + @property + def tmax(self): + """Last time point, or one per epoch when durations vary.""" + if self._variable_duration: + return self._tmax_per_epoch + return self.times[-1] + + @property + def durations(self): + """Duration of each epoch in seconds.""" + if not self._variable_duration: + return np.full(len(self.events), self.times[-1] - self.times[0]) + return self._tmax_per_epoch - self._tmin_per_epoch + + def get_times(self, epoch=None): + """Return the time vector of one epoch, or of all of them. + + Parameters + ---------- + epoch : int | None + Index of the epoch. If ``None``, return one time vector per epoch. + + Returns + ------- + times : array | list of array + Time vector(s) in seconds. For fixed-duration epochs every entry is + :attr:`times`. + """ + if epoch is None: + return [self.get_times(ii) for ii in range(len(self.events))] + if not self._variable_duration: + return self.times + sfreq = float(self.info["sfreq"]) + start = int(round(self._tmin_per_epoch[epoch] * sfreq)) + stop = int(round(self._tmax_per_epoch[epoch] * sfreq)) + return np.arange(start, stop + 1) / sfreq + + def as_fixed(self, pad_value=np.nan): + """Return fixed-duration epochs spanning the union of all epochs. + + Parameters + ---------- + pad_value : float + Value used where an epoch does not extend across the full window. + + Returns + ------- + epochs : instance of EpochsArray + Fixed-duration epochs from ``min(tmin)`` to ``max(tmax)``. + n_contributing : array, shape (n_times,) + Number of epochs carrying real data at each time point. + + Notes + ----- + Shorter epochs are padded to reach the common window, so the number of + epochs contributing to any reduction varies across the window. That + makes the noise level time-dependent, and the scalar ``nave`` carried by + :class:`mne.Evoked` cannot express it, which matters wherever ``nave`` + scales a noise covariance. ``n_contributing`` is returned rather than + discarded so this stays visible; see + :gh:`14206` for the discussion. + """ + if not self._variable_duration: + return self.copy(), np.full(len(self.times), len(self.events)) + # the union vector, which `times` deliberately refuses to hand out + times = self._times_readonly + n_times = len(times) + data = np.full( + (len(self.events), len(self.ch_names), n_times), pad_value, dtype=float + ) + n_contributing = np.zeros(n_times, int) + sfreq = float(self.info["sfreq"]) + offset = int(round(times[0] * sfreq)) + ragged = self._data + assert ragged is not None # variable-duration epochs are always preloaded + for ii, epoch in enumerate(ragged): + start = int(round(self._tmin_per_epoch[ii] * sfreq)) - offset + stop = start + epoch.shape[1] + data[ii, :, start:stop] = epoch + n_contributing[start:stop] += 1 + out = EpochsArray( + data, + self.info.copy(), + events=self.events.copy(), + event_id=self.event_id, + tmin=float(times[0]), + metadata=self.metadata, + selection=self.selection, + drop_log=self.drop_log, + verbose=False, + ) + return out, n_contributing + + def _get_variable_data(self, *, picks=None, item=None, copy=True): + """Return per-epoch data when durations vary, as one array per epoch.""" + if item is None: + item = slice(None) + sel = np.arange(len(self.events))[item] if not isinstance(item, str) else None + if sel is None: + raise NotImplementedError( + "Selecting variable-duration epochs by condition name is not " + "implemented; index with integers or a slice." + ) + ch_idx = _picks_to_idx(self.info, picks, "all", exclude=()) + ragged = self._data + assert ragged is not None # variable-duration epochs are always preloaded + out = [ragged[ii][ch_idx] for ii in sel] + return [epoch.copy() for epoch in out] if copy else out + def load_data(self) -> Self: """Load the data if not already preloaded. @@ -789,13 +1105,19 @@ def load_data(self) -> Self: """ if self.preload: return self - self._data = self._get_data() + if self._variable_duration: + self._data = self._load_variable_from_raw() + else: + self._data = self._get_data() self.preload = True self._do_baseline = False self._decim_slice = slice(None, None, None) self._decim = 1 - self._raw_times = self.times - assert self._data.shape[-1] == len(self.times) + if not self._variable_duration: + # variable-duration epochs already set _raw_times to the union span + # in __init__, and self.times refuses to answer for them + self._raw_times = self.times + assert self._data.shape[-1] == len(self.times) self._raw = None # shouldn't need it anymore return self @@ -953,14 +1275,59 @@ def _reject_setup(self, reject, flat, *, allow_callable=False): reject_imax = idxs[-1] self._reject_time = slice(reject_imin, reject_imax) + def _load_variable_from_raw(self): + """Read every epoch at its own length into a list, mirroring ``drop_bad``.""" + detrend_picks = self._detrend_picks + drop_log = list(self.drop_log) + good_idx, out = [], [] + for idx, sel in enumerate(self.selection): + epoch_noproj = self._get_epoch_from_raw(idx) + epoch_noproj = self._detrend_offset_decim(epoch_noproj, detrend_picks) + epoch = self._project_epoch(epoch_noproj) + epoch_out = epoch_noproj if self._do_delayed_proj else epoch + is_good, bad_tuple = self._is_good_epoch( + epoch, n_times=self._n_times_per_epoch(idx) + ) + if not is_good: + drop_log[sel] = drop_log[sel] + bad_tuple + continue + good_idx.append(idx) + out.append(epoch_out) + + good_idx = np.asarray(good_idx, dtype=int) + if len(good_idx) != len(self.events): + self._tmin_per_epoch = self._tmin_per_epoch[good_idx] + self._tmax_per_epoch = self._tmax_per_epoch[good_idx] + self.events = self.events[good_idx] + self.selection = self.selection[good_idx] + if self.metadata is not None: + GetEpochsMixin.metadata.fset( + self, self.metadata.iloc[good_idx], verbose=False + ) + self.drop_log = tuple(drop_log) + self._bad_dropped = True + logger.info(f"{len(out)} matching events found after rejection") + return out + + def _n_times_per_epoch(self, idx): + """Return the number of samples in one epoch.""" + if not self._variable_duration: + return len(self.times) + sfreq = float(self.info["sfreq"]) + return ( + int(round((self._tmax_per_epoch[idx] - self._tmin_per_epoch[idx]) * sfreq)) + + 1 + ) + @verbose # verbose is used by mne-realtime - def _is_good_epoch(self, data, verbose=None): + def _is_good_epoch(self, data, verbose=None, *, n_times=None): """Determine if epoch is good.""" if isinstance(data, str): return False, (data,) if data is None: return False, ("NO_DATA",) - n_times = len(self.times) + if n_times is None: + n_times = len(self.times) if data.shape[1] < n_times: # epoch is too short ie at the end of the data return False, ("TOO_SHORT",) @@ -1513,6 +1880,17 @@ def drop_bad( flat = self.flat if any(isinstance(rej, str) and rej != "existing" for rej in (reject, flat)): raise ValueError('reject and flat, if strings, must be "existing"') + if self._variable_duration and (reject or flat): + # the no-arg call has already short-circuited above, so reaching + # here means real thresholds + raise NotImplementedError( + "drop_bad() with reject or flat is not implemented for " + "variable-duration epochs. Amplitude rejection is per-trial " + "and could work here, but the preloaded path needs one shared " + "time axis, which these epochs do not have. Pass reject= to " + "Epochs() at construction instead, which does apply it per " + "epoch. See https://github.com/mne-tools/mne-python/issues/14206." + ) self._reject_setup(reject, flat, allow_callable=True) self._get_data(out=False, verbose=verbose) return self @@ -1873,7 +2251,15 @@ def _get_data( epoch = self._project_epoch(epoch_noproj) epoch_out = epoch_noproj if self._do_delayed_proj else epoch - is_good, bad_tuple = self._is_good_epoch(epoch, verbose=verbose) + is_good, bad_tuple = self._is_good_epoch( + epoch, + verbose=verbose, + n_times=( + self._n_times_per_epoch(idx) + if self._variable_duration + else None + ), + ) if not is_good: assert isinstance(bad_tuple, tuple) assert all(isinstance(x, str) for x in bad_tuple) @@ -1994,9 +2380,12 @@ def get_data( *, copy: bool = True, verbose: bool | str | int | None = None, - ) -> np.ndarray: + ) -> np.ndarray | list[np.ndarray]: """Get all epochs as a 3D array. + Epochs of different durations have no shared time axis, so those are + returned as one array per epoch instead; see :meth:`as_fixed`. + Parameters ---------- %(picks_all)s @@ -2039,10 +2428,20 @@ def get_data( Returns ------- - data : array of shape (n_epochs, n_channels, n_times) + data : array of shape (n_epochs, n_channels, n_times) | list of array The epochs data. Will be a copy when ``copy=True`` and will be a view - when possible when ``copy=False``. + when possible when ``copy=False``. When durations vary, one + ``(n_channels, n_times_i)`` array per epoch. """ + if self._variable_duration: + for name, value in (("units", units), ("tmin", tmin), ("tmax", tmax)): + if value is not None: + raise NotImplementedError( + f"get_data() with {name} is not implemented for " + "variable-duration epochs; it would have to be applied " + "per epoch. Call as_fixed() first to get one array." + ) + return self._get_variable_data(picks=picks, item=item, copy=copy) return self._get_data( picks=picks, item=item, units=units, tmin=tmin, tmax=tmax, copy=copy ) @@ -2158,8 +2557,13 @@ def __repr__(self): """Build string representation.""" s = f"{len(self.events)} events " s += "(all good)" if self._bad_dropped else "(good & bad)" - s += f", {self.tmin:.3f}".rstrip("0").rstrip(".") - s += f" – {self.tmax:.3f}".rstrip("0").rstrip(".") + if self._variable_duration: + s += f", {self.tmin.min():.3f}".rstrip("0").rstrip(".") + s += f" – {self.tmax.max():.3f}".rstrip("0").rstrip(".") + s += f" s, durations {self.durations.min():.3f}–{self.durations.max():.3f}" + else: + s += f", {self.tmin:.3f}".rstrip("0").rstrip(".") + s += f" – {self.tmax:.3f}".rstrip("0").rstrip(".") s += " s (baseline " if self.baseline is None: s += "off" @@ -2250,6 +2654,10 @@ def crop( # XXX this could be made to work on non-preloaded data... _check_preload(self, "Modifying data of epochs") + if self._variable_duration: + self._crop_variable(tmin, tmax, include_tmax) + return self + super().crop(tmin=tmin, tmax=tmax, include_tmax=include_tmax) # Adjust rejection period @@ -2267,6 +2675,98 @@ def crop( self.reject_tmax = self.tmax return self + def _crop_variable(self, tmin, tmax, include_tmax): + """Crop each epoch on its own time axis.""" + for name in ("reject_tmin", "reject_tmax"): + if getattr(self, name, None) is not None: + raise NotImplementedError( + f"{name} is not implemented for variable-duration epochs, " + "because the window is not guaranteed to exist in every " + "epoch." + ) + + sfreq = float(self.info["sfreq"]) + # Work out every selection before changing anything, so a window that + # misses one epoch leaves the object as it was. + masks = list() + clamped_tmin = clamped_tmax = False + for ii in range(len(self.events)): + times = self.get_times(ii) + this_tmin, this_tmax = tmin, tmax + this_include_tmax = include_tmax + if this_tmin is None: + this_tmin = times[0] + elif this_tmin < times[0]: + clamped_tmin = True + this_tmin = times[0] + if this_tmax is None: + this_tmax = times[-1] + elif this_tmax > times[-1]: + clamped_tmax = True + this_tmax = times[-1] + # matches the fixed path: a clamped end keeps its last sample + this_include_tmax = True + # _time_mask raises when the window is inverted, which is what an + # entirely out-of-range request collapses to once tmax is clamped + mask = _time_mask( + times, + this_tmin, + this_tmax, + sfreq=sfreq, + include_tmax=this_include_tmax, + ) + if not mask.any(): + raise ValueError( + f"tmin ({tmin}) and tmax ({tmax}) would leave epoch {ii} " + f"with no samples; it spans {times[0]:g} to {times[-1]:g} s." + ) + masks.append(mask) + + # One warning per bound however many epochs needed clamping + if clamped_tmin: + warn( + "tmin is not in time interval for every epoch. tmin is set to " + "each of those epochs' own first sample." + ) + if clamped_tmax: + warn( + "tmax is not in time interval for every epoch. tmax is set to " + "each of those epochs' own last sample." + ) + + ragged = self._data + assert ragged is not None # variable-duration epochs are always preloaded + starts = np.empty(len(masks)) + stops = np.empty(len(masks)) + for ii, mask in enumerate(masks): + kept = self.get_times(ii)[mask] + ragged[ii] = ragged[ii][..., mask] + starts[ii], stops[ii] = kept[0], kept[-1] + self._tmin_per_epoch = starts + self._tmax_per_epoch = stops + + # Cropping can remove the variation entirely. Compare the axes by their + # sample index and length rather than by float equality: every epoch is + # regularly sampled at one sfreq, so that pair identifies an axis + # exactly and does not depend on how the bound was rounded. + first_idx = np.round(starts * sfreq).astype(int) + lengths = np.array([epoch.shape[-1] for epoch in ragged]) + if ( + len(masks) + and (first_idx == first_idx[0]).all() + and (lengths == lengths[0]).all() + ): + self._data = np.stack(list(ragged)) + self._variable_duration = False + self._tmin_per_epoch = None # ty: ignore[invalid-assignment] + self._tmax_per_epoch = None # ty: ignore[invalid-assignment] + start_idx, stop_idx = first_idx[0], first_idx[0] + lengths[0] - 1 + else: + start_idx = int(round(float(np.min(starts)) * sfreq)) + stop_idx = int(round(float(np.max(stops)) * sfreq)) + self._raw_times = np.arange(start_idx, stop_idx + 1) / sfreq + self._set_times(self._raw_times) + def copy(self) -> Self: """Return copy of Epochs instance. @@ -2374,6 +2874,7 @@ def save( total_size = 0 else: d = self[0].get_data(copy=False) + assert isinstance(d, np.ndarray) # save() refuses ragged epochs # this should be guaranteed by subclasses assert d.dtype in (">f8", "c16", " 0 + and _tmin_arr.size > 0 + and not np.all(_tmin_arr == _tmin_arr.flat[0]) + ) + if not variable and _tmin_arr.ndim > 0: + tmin = float(_tmin_arr.flat[0]) if _tmin_arr.size else 0.0 + if variable: + data = [np.asanyarray(epoch) for epoch in data] + dtype = ( + np.complex128 + if any(np.any(np.iscomplex(epoch)) for epoch in data) + else np.float64 + ) + data = [epoch.astype(dtype, copy=False) for epoch in data] + for ii, epoch in enumerate(data): + if epoch.ndim != 2: + raise ValueError( + f"Epoch {ii} must be 2D (n_channels, n_samples), got " + f"{epoch.ndim}D." + ) + if len(info["ch_names"]) != epoch.shape[0]: + raise ValueError("Info and data must have same number of channels.") + if events is None: + events = _gen_events(len(data)) + info = info.copy() # do not modify original info + tmin = np.broadcast_to(np.asarray(tmin, dtype=float), (len(data),)).copy() + tmax = np.array( + [ + (epoch.shape[1] - 1) / info["sfreq"] + start + for epoch, start in zip(data, tmin) + ] ) + else: + dtype = np.complex128 if np.any(np.iscomplex(data)) else np.float64 + data = np.asanyarray(data, dtype=dtype) + if data.ndim != 3: + raise ValueError( + "Data must be a 3D array of shape (n_epochs, n_channels, n_samples)" + ) - if len(info["ch_names"]) != data.shape[1]: - raise ValueError("Info and data must have same number of channels.") - if events is None: - n_epochs = len(data) - events = _gen_events(n_epochs) - info = info.copy() # do not modify original info - tmax = (data.shape[2] - 1) / info["sfreq"] + tmin + if len(info["ch_names"]) != data.shape[1]: + raise ValueError("Info and data must have same number of channels.") + if events is None: + n_epochs = len(data) + events = _gen_events(n_epochs) + info = info.copy() # do not modify original info + tmax = (data.shape[2] - 1) / info["sfreq"] + tmin super().__init__( info, diff --git a/mne/tests/test_epochs_variable_duration.py b/mne/tests/test_epochs_variable_duration.py new file mode 100644 index 00000000000..347b24b9686 --- /dev/null +++ b/mne/tests/test_epochs_variable_duration.py @@ -0,0 +1,701 @@ +"""Tests for epochs whose trials have different durations.""" + +# Authors: The MNE-Python contributors. +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +import warnings + +import numpy as np +import pytest +from numpy.testing import assert_allclose, assert_array_equal + +from mne import EpochsArray, create_info +from mne.epochs import ( + _VARIABLE_FALLBACK, + _VARIABLE_NEEDS_POLICY, + _VARIABLE_NOT_IMPLEMENTED, +) + +SFREQ = 100.0 +CH_NAMES = ["a", "b", "c"] + + +def _events(n_epochs): + """Return a simple event array.""" + return np.c_[ + np.arange(n_epochs) * 500 + 200, + np.zeros(n_epochs, int), + np.ones(n_epochs, int), + ] + + +def _make(tmin, tmax, seed=0): + """Build variable-duration epochs from per-event bounds.""" + rng = np.random.default_rng(seed) + tmin = np.asarray(tmin, dtype=float) + tmax = np.asarray(tmax, dtype=float) + data = [ + rng.standard_normal((len(CH_NAMES), round((b - a) * SFREQ) + 1)) * 1e-6 + for a, b in zip(tmin, tmax) + ] + info = create_info(CH_NAMES, SFREQ, "eeg") + return EpochsArray( + data, + info, + events=_events(len(data)), + tmin=tmin, + baseline=None, + verbose=False, + ) + + +@pytest.fixture +def variable(): + """Four epochs sharing tmin with different tmax.""" + return _make(np.full(4, -0.2), [0.5, 0.9, 0.7, 0.6]) + + +# -- the scalar path must be untouched ------------------------------------ +def test_scalar_path_unchanged(): + """Test that fixed-duration epochs behave exactly as before.""" + rng = np.random.default_rng(0) + data = rng.standard_normal((5, len(CH_NAMES), 71)) * 1e-6 + info = create_info(CH_NAMES, SFREQ, "eeg") + epochs = EpochsArray(data, info, events=_events(5), tmin=-0.2, verbose=False) + + assert not epochs.variable_duration + assert isinstance(epochs.tmin, float) + assert isinstance(epochs.tmax, float) + assert_allclose(epochs.tmin, -0.2) + assert epochs.get_data().shape == (5, len(CH_NAMES), 71) + assert epochs.average().data.shape == (len(CH_NAMES), 71) + # durations is new but must be defined for the fixed case too + assert_allclose(epochs.durations, np.full(5, 0.7)) + + +def test_equal_bounds_arrays_match_scalar(): + """Test that per-event bounds that happen to be equal match the scalar path.""" + rng = np.random.default_rng(3) + data = rng.standard_normal((4, len(CH_NAMES), 71)) * 1e-6 + info = create_info(CH_NAMES, SFREQ, "eeg") + + scalar = EpochsArray(data, info, events=_events(4), tmin=-0.2, verbose=False) + arrays = EpochsArray( + list(data), info, events=_events(4), tmin=np.full(4, -0.2), verbose=False + ) + + # a list of equal-length arrays is not ragged, so this stays the scalar path + assert not arrays.variable_duration + assert_allclose(arrays.times, scalar.times) + assert_allclose(arrays.get_data(), scalar.get_data()) + + +# -- construction --------------------------------------------------------- +def test_bounds_and_durations(variable): + """Test the per-epoch bounds derived from the data.""" + assert variable.variable_duration + assert_allclose(variable.tmin, np.full(4, -0.2)) + assert_allclose(variable.tmax, [0.5, 0.9, 0.7, 0.6], atol=1e-12) + assert_allclose(variable.durations, [0.7, 1.1, 0.9, 0.8], atol=1e-12) + + +def test_per_epoch_time_vectors(variable): + """Test that each epoch reports its own axis, since none is shared.""" + per_epoch = variable.get_times() + assert len(per_epoch) == 4 + for ii, epoch_times in enumerate(per_epoch): + assert_allclose(epoch_times[0], variable.tmin[ii]) + assert_allclose(epoch_times[-1], variable.tmax[ii], atol=1e-12) + assert len(epoch_times) == variable.get_data()[ii].shape[-1] + + +def test_get_data_returns_one_array_per_epoch(variable): + """Test that data comes back per epoch, since no common length exists.""" + data = variable.get_data() + assert isinstance(data, list) + lengths = [epoch.shape[1] for epoch in data] + assert lengths == [71, 111, 91, 81] + assert all(epoch.shape[0] == len(CH_NAMES) for epoch in data) + + picked = variable.get_data(picks=["a", "c"]) + assert all(epoch.shape[0] == 2 for epoch in picked) + + +def test_get_data_refuses_what_it_cannot_honour(variable): + """Test that units, tmin and tmax raise rather than being ignored.""" + for name, kwargs in ( + ("units", dict(units="uV")), + ("tmin", dict(tmin=0.0)), + ("tmax", dict(tmax=0.4)), + ): + with pytest.raises(NotImplementedError, match=rf"get_data\(\) with {name}"): + variable.get_data(**kwargs) + # all three work on the padded copy the message points at + fixed, _ = variable.as_fixed() + volts = fixed.get_data(tmin=0.0, tmax=0.4) + assert volts.shape == (4, len(CH_NAMES), 40) + micro = fixed.get_data(units="uV", tmin=0.0, tmax=0.4) + assert_allclose(micro, volts * 1e6) + + +def test_only_time_may_vary(): + """Test that a varying channel count is rejected.""" + info = create_info(CH_NAMES, SFREQ, "eeg") + data = [np.zeros((3, 71)), np.zeros((2, 111))] + with pytest.raises(ValueError, match="same number of channels"): + EpochsArray(data, info, events=_events(2), tmin=np.zeros(2), verbose=False) + + +def test_bounds_validation(): + """Test the checks on per-event bounds.""" + from mne.epochs import _check_variable_bounds + + tmin, tmax, variable = _check_variable_bounds(-0.2, 0.5, 3) + assert not variable + assert tmin == -0.2 + + tmin, tmax, variable = _check_variable_bounds(-0.2, np.array([0.5, 0.9, 0.7]), 3) + assert variable + assert_array_equal(tmin, np.full(3, -0.2)) + + with pytest.raises(ValueError, match="entries but there are"): + _check_variable_bounds(np.zeros(2), 0.5, 3) + with pytest.raises(ValueError, match="less than or equal to tmax"): + _check_variable_bounds(np.array([0.5, 0.0]), np.array([0.1, 1.0]), 2) + with pytest.raises(ValueError, match="must be finite"): + _check_variable_bounds(np.array([np.nan, 0.0]), np.array([1.0, 1.0]), 2) + + +@pytest.mark.parametrize( + "kwargs, match", + [ + (dict(baseline=(None, 0)), "Baseline correction is not implemented"), + (dict(reject_tmin=-0.1), "reject_tmin is not implemented"), + (dict(reject_tmax=0.1), "reject_tmax is not implemented"), + ], +) +def test_unsupported_options_raise(kwargs, match): + """Test that options outside this implementation are refused clearly.""" + rng = np.random.default_rng(0) + tmin = np.full(2, -0.2) + data = [rng.standard_normal((len(CH_NAMES), n)) for n in (71, 111)] + info = create_info(CH_NAMES, SFREQ, "eeg") + with pytest.raises(NotImplementedError, match=match): + EpochsArray(data, info, events=_events(2), tmin=tmin, verbose=False, **kwargs) + + +# -- as_fixed ------------------------------------------------------------- +def test_as_fixed_spans_union_and_reports_support(variable): + """Test padding to a common window and the contributor count it implies.""" + fixed, n_contributing = variable.as_fixed() + + assert not fixed.variable_duration + n_union = max(epoch.shape[-1] for epoch in variable.get_data()) + assert fixed.get_data().shape == (4, len(CH_NAMES), n_union) + assert_allclose(fixed.times[0], variable.tmin.min()) + assert_allclose(fixed.times[-1], variable.tmax.max()) + + # every epoch covers the start; only the longest reaches the end + assert n_contributing[0] == 4 + assert n_contributing[-1] == 1 + assert n_contributing.min() < n_contributing.max() + # support never grows once epochs start dropping out + peak = int(np.argmax(n_contributing)) + assert np.all(np.diff(n_contributing[peak:]) <= 0) + + +def test_as_fixed_preserves_data_and_pads_the_rest(variable): + """Test that real samples survive and the remainder is marked.""" + fixed, _ = variable.as_fixed() + dense = fixed.get_data() + for ii, epoch in enumerate(variable.get_data()): + n = epoch.shape[1] + assert_allclose(dense[ii, :, :n], epoch) + assert np.isnan(dense[ii, :, n:]).all() + + +def test_as_fixed_on_fixed_epochs_is_a_copy(): + """Test that as_fixed is defined, and trivial, for fixed-duration epochs.""" + rng = np.random.default_rng(1) + data = rng.standard_normal((3, len(CH_NAMES), 71)) * 1e-6 + info = create_info(CH_NAMES, SFREQ, "eeg") + epochs = EpochsArray(data, info, events=_events(3), tmin=-0.2, verbose=False) + + fixed, n_contributing = epochs.as_fixed() + assert_allclose(fixed.get_data(), epochs.get_data()) + assert_array_equal(n_contributing, np.full(len(epochs.times), 3)) + + +# -- dispatch -------------------------------------------------------------- +@pytest.mark.parametrize("meth", sorted(_VARIABLE_NEEDS_POLICY)) +def test_reductions_ask_for_a_policy(variable, meth): + """Test that combining epochs across a time axis they lack is refused. + + Padding first and reducing afterwards is not a slower answer, it is a + different one: one short epoch turns a whole time point into NaN, and the + scalar ``nave`` keeps reporting the full count. + """ + with pytest.raises(NotImplementedError, match="explicit policy"): + result = getattr(variable, meth)() + list(result) # iter_evoked is a generator + + +def test_policy_message_names_the_varying_count(variable): + """Test that the refusal explains itself rather than just declining.""" + with pytest.raises(NotImplementedError) as excinfo: + variable.average() + message = str(excinfo.value) + assert "varies across the window" in message + assert "as_fixed" in message + + +def test_compute_tfr_does_not_silently_pad(variable): + """Test that the transform is not quietly given padded data. + + Padding before a time-frequency transform is the opposite of the order this + work argues for, which is to transform at native duration and warp the + result. Doing it silently inside ``compute_tfr`` would ship the thing being + argued against. + """ + with pytest.raises(NotImplementedError, match="explicit policy"): + variable.compute_tfr("morlet", freqs=np.arange(10.0, 20.0, 2.0), n_cycles=2) + + +@pytest.mark.parametrize("meth", sorted(_VARIABLE_NOT_IMPLEMENTED)) +def test_per_trial_methods_raise_until_implemented(variable, meth): + """Test that per-trial work refuses rather than running on a padded copy.""" + with pytest.raises(NotImplementedError, match="not implemented"): + getattr(variable, meth)() + + +@pytest.mark.parametrize("meth", sorted(_VARIABLE_FALLBACK)) +def test_display_methods_warn_and_fall_back(variable, meth): + """Test that the remaining inspection method degrades rather than refuses.""" + if meth == "to_data_frame": + pytest.importorskip("pandas") + with pytest.warns(RuntimeWarning, match="ran on as_fixed"): + assert getattr(variable, meth)() is not None + + +def test_plot_is_not_a_fallback(): + """Test that browsing is native, not padded (see mne/viz/tests/test_epochs).""" + assert "plot" not in _VARIABLE_FALLBACK + assert "plot" not in _VARIABLE_NEEDS_POLICY + assert "plot" not in _VARIABLE_NOT_IMPLEMENTED + + +# -- operations that stay native ------------------------------------------- +def test_pick_keeps_durations(variable): + """Test that channel selection leaves the time axis alone.""" + before = variable.durations.copy() + picked = variable.copy().pick(["a", "c"]) + assert picked.ch_names == ["a", "c"] + assert_allclose(picked.durations, before) + for epoch in picked.get_data(): + assert epoch.shape[0] == 2 + + +def test_getitem_keeps_per_epoch_bounds(variable): + """Test that indexing carries the bounds with the epochs.""" + subset = variable[[0, 2]] + assert len(subset) == 2 + assert_allclose(subset.durations, variable.durations[[0, 2]]) + for got, want in zip(subset.get_data(), [variable.get_data()[i] for i in (0, 2)]): + assert_array_equal(got, want) + + +def test_drop_keeps_per_epoch_bounds(variable): + """Test that dropping an epoch drops its bounds too.""" + kept = variable.copy().drop([1]) + assert len(kept) == 3 + assert_allclose(kept.durations, variable.durations[[0, 2, 3]]) + + +def test_shift_time_moves_bounds_not_samples(variable): + """Test that shifting the origin does not resample anything.""" + before_lengths = [epoch.shape[1] for epoch in variable.get_data()] + before_durations = variable.durations.copy() + shifted = variable.copy().shift_time(0.1) + assert_allclose(shifted.tmin, variable.tmin + 0.1) + assert_allclose(shifted.durations, before_durations) + assert [epoch.shape[1] for epoch in shifted.get_data()] == before_lengths + + +# -- the time axis --------------------------------------------------------- +def test_times_refuses_to_invent_a_shared_axis(variable): + """Test that ``times`` raises rather than handing back the union. + + ``len(epochs.times) == data.shape[-1]`` has always held. Returning the union + would leave it false while looking ordinary. + """ + with pytest.raises(RuntimeError, match="no time axis they share"): + variable.times + with pytest.raises(RuntimeError, match="get_times"): + variable.times + + +def test_union_axis_is_reachable_through_as_fixed(variable): + """Test that the padded common axis is available when asked for.""" + fixed, _ = variable.as_fixed() + assert_allclose(fixed.times[0], variable.tmin.min()) + assert_allclose(fixed.times[-1], variable.tmax.max()) + assert len(fixed.times) == fixed.get_data().shape[-1] + + +def test_fixed_epochs_still_have_times(): + """Test that the scalar path is untouched by any of this.""" + rng = np.random.default_rng(2) + info = create_info(CH_NAMES, SFREQ, "eeg") + epochs = EpochsArray( + rng.standard_normal((3, len(CH_NAMES), 71)) * 1e-6, + info, + events=_events(3), + tmin=-0.2, + verbose=False, + ) + assert len(epochs.times) == 71 + assert epochs.average().data.shape == (len(CH_NAMES), 71) + + +def test_nothing_reaches_the_user_as_an_internal_error(variable): + """Test that no public method leaks a NumPy error about lists.""" + import warnings + + names = ( + sorted(_VARIABLE_FALLBACK) + + sorted(_VARIABLE_NEEDS_POLICY) + + sorted(_VARIABLE_NOT_IMPLEMENTED) + ) + for name in names: + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + try: + getattr(variable.copy(), name)() + except (NotImplementedError, RuntimeError): + pass + except TypeError as exc: + assert "argument" in str(exc), f"{name}: {exc}" + except (AttributeError, IndexError) as exc: + raise AssertionError(f"{name} leaked an internal error: {exc}") + except Exception: + pass + + +# -- construction from Raw ------------------------------------------------- +def _raw(n_seconds=30.0, sfreq=SFREQ): + """Return a small continuous recording.""" + from mne.io import RawArray + + rng = np.random.default_rng(3) + info = create_info(CH_NAMES, sfreq, "eeg") + return RawArray( + rng.standard_normal((len(CH_NAMES), int(n_seconds * sfreq))) * 1e-6, + info, + verbose=False, + ) + + +def test_from_raw_preserves_each_slice(): + """Test that epochs read from Raw match the samples they came from.""" + from mne import Epochs + + raw = _raw() + onsets = np.array([100, 700, 1400, 2100]) + events = np.c_[onsets, np.zeros(4, int), np.ones(4, int)] + tmax = np.array([0.5, 1.9, 0.9, 2.4]) + + epochs = Epochs( + raw, + events, + {"a": 1}, + tmin=np.zeros(4), + tmax=tmax, + baseline=None, + preload=True, + verbose=False, + ) + + assert epochs.variable_duration + assert len(epochs) == 4 + # +1 because both endpoints are inclusive + want_n = np.round(tmax * SFREQ).astype(int) + 1 + data = epochs.get_data() + assert [d.shape[-1] for d in data] == list(want_n) + + raw_data = raw.get_data() + for i, start in enumerate(onsets): + assert_array_equal(data[i], raw_data[:, start : start + want_n[i]]) + + +def test_from_raw_load_data_leaves_raw_times_alone(): + """Test that loading from Raw does not ask a ragged object for one axis.""" + from mne import Epochs + + raw = _raw() + events = np.c_[[100, 900], [0, 0], [1, 1]] + # load_data() used to raise here from `self._raw_times = self.times` + epochs = Epochs( + raw, + events, + {"a": 1}, + tmin=np.zeros(2), + tmax=np.array([0.5, 2.0]), + baseline=None, + preload=True, + verbose=False, + ) + assert epochs.preload + # _raw_times spans the union, which is what as_fixed() lays epochs onto + fixed, _ = epochs.as_fixed() + assert len(epochs._raw_times) == fixed.get_data().shape[-1] + + +def test_from_raw_requires_preload(): + """Test that lazy reading is refused with a reason, not an internal error.""" + from mne import Epochs + + raw = _raw() + events = np.c_[[100, 900], [0, 0], [1, 1]] + with pytest.raises(NotImplementedError, match="must be preloaded"): + Epochs( + raw, + events, + {"a": 1}, + tmin=np.zeros(2), + tmax=np.array([0.5, 2.0]), + baseline=None, + preload=False, + verbose=False, + ) + + +def test_from_raw_scalar_bounds_still_scalar(): + """Test that equal per-event bounds do not switch on the ragged path.""" + from mne import Epochs + + raw = _raw() + events = np.c_[[100, 900], [0, 0], [1, 1]] + epochs = Epochs( + raw, + events, + {"a": 1}, + tmin=np.zeros(2), + tmax=np.full(2, 0.5), + baseline=None, + preload=True, + verbose=False, + ) + assert not epochs.variable_duration + assert isinstance(epochs.tmin, float) + assert epochs.get_data().shape == (2, len(CH_NAMES), 51) + + +@pytest.mark.parametrize( + "item", [slice(None, 2), slice(1, None), slice(None, None, 2), slice(None)] +) +def test_getitem_slice_selects_epochs_not_the_slice(variable, item): + """Test that slicing subsets the epochs rather than wrapping the slice.""" + want = np.arange(len(variable))[item] + subset = variable[item] + + assert len(subset) == len(want) + # a slice used to survive into the data list as a single nested element + assert all(isinstance(d, np.ndarray) and d.ndim == 2 for d in subset.get_data()) + assert_allclose(subset.durations, variable.durations[want]) + for got, idx in zip(subset.get_data(), want): + assert_array_equal(got, variable.get_data()[idx]) + + +def test_apply_function_refuses(variable): + """Test that apply_function refuses instead of indexing a list with a tuple.""" + with pytest.raises(NotImplementedError, match="not implemented"): + variable.apply_function(lambda x: x * 2) + + +def test_pick_does_not_reach_back_into_the_parent(variable): + """Test that picking replaces one object's epochs and no other's.""" + before = [epoch.shape for epoch in variable.get_data()] + + # _pick_drop_channels replaces the list contents in place, so anything + # sharing that list would be picked too + variable.copy().pick(["a", "c"]) + assert [epoch.shape for epoch in variable.get_data()] == before + + subset = variable[:2] + assert subset._data is not variable._data + subset.pick(["a"]) + assert [epoch.shape for epoch in variable.get_data()] == before + assert all(epoch.shape[0] == 1 for epoch in subset.get_data()) + + +# -- crop ------------------------------------------------------------------- +def _crop_oracle(epochs, idx, **kwargs): + """Crop epoch ``idx`` as an ordinary one-epoch Epochs, for comparison.""" + data = epochs.get_data()[idx] + tmin = np.atleast_1d(epochs.tmin)[idx] + one = EpochsArray( + data[None], + create_info(list(epochs.ch_names), SFREQ, "eeg"), + tmin=float(tmin), + baseline=None, + verbose=False, + ) + return one.crop(**kwargs) + + +def _assert_crop_matches_mne(epochs, **kwargs): + """Assert every cropped epoch equals ordinary MNE cropping it alone.""" + wanted = [_crop_oracle(epochs, ii, **kwargs) for ii in range(len(epochs))] + got = epochs.copy().crop(**kwargs) + data = got.get_data() + for ii, one in enumerate(wanted): + assert_allclose(data[ii], one.get_data()[0]) + times = got.get_times(ii) if got.variable_duration else got.times + assert_allclose(times, one.times) + return got + + +@pytest.mark.parametrize( + "kwargs", + [ + dict(tmin=0.0, tmax=0.4), # both bounds + dict(tmin=0.1), # only tmin + dict(tmax=0.45), # only tmax + dict(tmin=-0.05, tmax=0.45, include_tmax=True), + dict(tmin=-0.05, tmax=0.45, include_tmax=False), + ], +) +def test_crop_matches_mne_per_epoch(variable, kwargs): + """Test that cropping equals ordinary MNE applied to each epoch alone.""" + _assert_crop_matches_mne(variable, **kwargs) + + +def test_crop_matches_mne_with_unequal_tmin(kwargs=None): + """Test parity when both bounds differ between epochs.""" + epochs = _make([-0.2, -0.35, 0.0, -0.1], [0.5, 0.9, 0.7, 0.6]) + _assert_crop_matches_mne(epochs, tmin=0.05, tmax=0.4) + _assert_crop_matches_mne(epochs, tmax=0.5) + + +def test_crop_keeps_the_object_ragged(variable): + """Test that unequal durations survive a crop that does not equalise them.""" + cropped = variable.copy().crop(tmin=0.0) + assert cropped.variable_duration + assert len(np.unique(cropped.durations)) > 1 + assert isinstance(cropped._data, list) + + +def test_crop_clamps_each_epoch_and_warns_once(variable): + """Test that a bound past some epochs clamps per epoch, warning once.""" + before = variable.durations.copy() + with pytest.warns(RuntimeWarning, match="tmax is not in time interval") as rec: + cropped = variable.copy().crop(tmax=99.0) + assert len(rec) == 1 # not one per epoch + # each epoch kept everything it had + assert_allclose(cropped.durations, before) + assert cropped.variable_duration + + with pytest.warns(RuntimeWarning, match="tmin is not in time interval") as rec: + cropped = variable.copy().crop(tmin=-99.0) + assert len(rec) == 1 + assert_allclose(cropped.durations, before) + + +def test_crop_does_not_warn_when_nothing_is_clamped(variable): + """Test that a window inside every epoch is silent.""" + with warnings.catch_warnings(): + warnings.simplefilter("error") + variable.copy().crop(tmin=0.0, tmax=0.4) + + +def test_crop_clamped_tmax_keeps_the_last_sample(variable): + """Test that clamping tmax includes that epoch's final sample.""" + lengths = [epoch.shape[-1] for epoch in variable.get_data()] + with pytest.warns(RuntimeWarning, match="tmax is not in time interval"): + # include_tmax=False must not drop the endpoint that clamping produced + cropped = variable.copy().crop(tmax=99.0, include_tmax=False) + assert [epoch.shape[-1] for epoch in cropped.get_data()] == lengths + + +def test_crop_outside_every_sample_fails_cleanly(variable): + """Test that a window missing an epoch refuses and changes nothing.""" + before_data = [epoch.copy() for epoch in variable.get_data()] + before_tmin = np.array(variable.tmin) + before_tmax = np.array(variable.tmax) + with pytest.raises(ValueError, match="must be less than or equal to"): + variable.crop(tmin=5.0) + # the failure left the object exactly as it was + for got, want in zip(variable.get_data(), before_data): + assert_array_equal(got, want) + assert_array_equal(np.array(variable.tmin), before_tmin) + assert_array_equal(np.array(variable.tmax), before_tmax) + assert variable.variable_duration + + +def test_crop_bounds_come_from_retained_samples(variable): + """Test that the stored bounds are sample positions, not the request.""" + # only tmin, so the differing ends keep the object ragged + cropped = variable.copy().crop(tmin=0.013) + assert cropped.variable_duration + # the request fell between samples and was snapped to one + assert not np.isclose(np.atleast_1d(cropped.tmin)[0], 0.013) + for ii in range(len(cropped)): + times = cropped.get_times(ii) + assert times[0] == pytest.approx(np.atleast_1d(cropped.tmin)[ii]) + assert times[-1] == pytest.approx(np.atleast_1d(cropped.tmax)[ii]) + # and the axis still describes the block exactly + assert len(times) == cropped.get_data()[ii].shape[-1] + assert_allclose( + cropped.durations, np.atleast_1d(cropped.tmax) - np.atleast_1d(cropped.tmin) + ) + + +def test_crop_that_equalises_axes_returns_fixed_epochs(variable): + """Test that removing the variation gives an ordinary Epochs back.""" + cropped = variable.copy().crop(tmax=0.5) + assert not cropped.variable_duration + assert isinstance(cropped._data, np.ndarray) + assert cropped._tmin_per_epoch is None + assert cropped._tmax_per_epoch is None + assert isinstance(cropped.tmin, float) + assert isinstance(cropped.tmax, float) + # times is answerable again, and agrees with the data + assert len(cropped.times) == cropped.get_data().shape[-1] + # and the reductions come back on their own, without touching the tables + evoked = cropped.average() + assert evoked.nave == len(cropped) + assert_allclose(evoked.times, cropped.times) + + +def test_crop_keeps_epoch_bookkeeping(variable): + """Test that events, metadata and drop_log travel unchanged.""" + pd = pytest.importorskip("pandas") + variable.metadata = pd.DataFrame(dict(kind=list("abcd"))) + events = variable.events.copy() + drop_log = variable.drop_log + selection = variable.selection.copy() + + cropped = variable.copy().crop(tmin=0.0, tmax=0.4) + assert_array_equal(cropped.events, events) + assert cropped.drop_log == drop_log + assert_array_equal(cropped.selection, selection) + assert list(cropped.metadata["kind"]) == list("abcd") + + +def test_crop_does_not_fall_back_to_as_fixed(variable): + """Test that cropping is native, never a padded copy.""" + + def _boom(*args, **kwargs): + raise AssertionError("crop() fell back to as_fixed()") + + variable.as_fixed = _boom + cropped = variable.crop(tmin=0.0) + assert not np.isnan(np.concatenate(cropped.get_data(), axis=-1)).any() + + +def test_crop_refuses_when_rejection_windows_are_set(variable): + """Test that a stray rejection window is refused, not compared to an array.""" + variable.reject_tmin = 0.0 # the constructor forbids this; be defensive + with pytest.raises(NotImplementedError, match="reject_tmin is not implemented"): + variable.crop(tmin=0.0) diff --git a/mne/utils/mixin.py b/mne/utils/mixin.py index 3addd688797..b276f524c5a 100644 --- a/mne/utils/mixin.py +++ b/mne/utils/mixin.py @@ -222,7 +222,11 @@ def _getitem( self._sanity_check_event_id() inst = self.copy() if copy else self if self._data is not None: - np.copyto(inst._data, self._data, casting="no") + if isinstance(self._data, list): + # one array per epoch, so there is nothing to copy into + inst._data = [d.copy() for d in self._data] if copy else self._data + else: + np.copyto(inst._data, self._data, casting="no") del self select = inst._item_to_select(item) @@ -257,9 +261,29 @@ def _getitem( # will reset the index for us GetEpochsMixin.metadata.fset(inst, metadata, verbose=False) if inst.preload and select_data: - # ensure that each Epochs instance owns its own data so we can - # resize later if necessary - inst._data = np.require(inst._data[select], requirements=["O"]) + if isinstance(inst._data, list): + # `select` can still be a slice here, and iterating one yields + # the slice itself rather than the epochs it covers + inst._data = [ + inst._data[ii] for ii in np.arange(len(inst._data))[select] + ] + else: + # ensure that each Epochs instance owns its own data so we can + # resize later if necessary + inst._data = np.require(inst._data[select], requirements=["O"]) + if getattr(inst, "_variable_duration", False): + # an ndarray takes a slice or an index array equally well + inst._tmin_per_epoch = inst._tmin_per_epoch[select] + inst._tmax_per_epoch = inst._tmax_per_epoch[select] + # the union window is defined by those bounds, so it has to be + # re-derived here as `crop` does; otherwise `as_fixed` keeps + # padding out to epochs that are no longer present + if len(inst._tmin_per_epoch): + sfreq = float(inst.info["sfreq"]) + start_idx = int(round(inst._tmin_per_epoch.min() * sfreq)) + stop_idx = int(round(inst._tmax_per_epoch.max() * sfreq)) + inst._raw_times = np.arange(start_idx, stop_idx + 1) / sfreq + inst._set_times(inst._raw_times) if drop_event_id: # update event id to reflect new content of inst inst.event_id = { @@ -767,6 +791,20 @@ def shift_time(self, tshift, relative=True) -> Self: or change the *data* values in any way. """ _check_preload(self, "shift_time") + if getattr(self, "_variable_duration", False): + # each epoch keeps its own length; only the origin moves, so this + # shifts the bounds and leaves the samples alone + if relative: + shift = tshift + else: + shift = tshift - self._tmin_per_epoch.min() + self._tmin_per_epoch = self._tmin_per_epoch + shift + self._tmax_per_epoch = self._tmax_per_epoch + shift + self._raw_times = self._raw_times + shift + self._set_times(self._raw_times) + self._update_first_last() + return self + start = tshift + (self.times[0] if relative else 0.0) new_times = start + np.arange(len(self.times)) / self.info["sfreq"] self._set_times(new_times) diff --git a/mne/viz/_figure.py b/mne/viz/_figure.py index 0579820883b..2cbe146f715 100644 --- a/mne/viz/_figure.py +++ b/mne/viz/_figure.py @@ -44,6 +44,15 @@ def __init__(self, **kwargs): vars(self).update(**kwargs) +def _epoch_window(boundary_times, start_ix, n_epochs): + """Return the start time and duration of a window of whole epochs.""" + n_total = len(boundary_times) - 1 + n_epochs = int(np.clip(n_epochs, 1, n_total)) + start_ix = int(np.clip(start_ix, 0, n_total - n_epochs)) + stop_ix = start_ix + n_epochs + return boundary_times[start_ix], boundary_times[stop_ix] - boundary_times[start_ix] + + class BrowserBase(ABC): """A base class containing for the 2D browser. @@ -78,7 +87,13 @@ def __init__(self, **kwargs): f"Expected an instance of Raw, Epochs, or ICA, got {type(inst)}." ) - if len(inst.times) < 2: + # variable-duration epochs have no one time axis, so count the samples + # the browser will lay end to end instead + if self.mne.instance_type == "epochs": + n_inst_times = self.mne.n_times + else: + n_inst_times = len(inst.times) + if n_inst_times < 2: raise ValueError( "Data from at least two time points are required to open the browser." ) @@ -120,7 +135,12 @@ def __init__(self, **kwargs): self.mne.epoch_traces = list() self.mne.bad_epochs = list() if inst is not None: - self.mne.sampling_period = np.diff(inst.times[:2])[0] / inst.info["sfreq"] + # NB: this is 1 / sfreq**2, not a sampling period; it is only ever + # used as a small nudge before searchsorted on boundary_times, so + # keep the value while deriving it without touching inst.times + # (which variable-duration epochs refuse to provide). + sfreq = inst.info["sfreq"] + self.mne.sampling_period = (1.0 / sfreq) / sfreq # annotations self.mne.annotations = list() self.mne.hscroll_annotations = list() @@ -157,6 +177,23 @@ def __init__(self, **kwargs): self.mne.midpoints = ( np.convolve(self.mne.boundary_times, np.ones(2), mode="valid") / 2 ) + # callers that only ever deal with equal-length epochs (ICA sources) + # do not supply these, so derive them from the boundaries + sfreq = self.mne.info["sfreq"] + if not hasattr(self.mne, "boundary_samples"): + self.mne.boundary_samples = np.round( + np.asarray(self.mne.boundary_times) * sfreq + ).astype(int) + if not hasattr(self.mne, "epoch_tmins"): + n_epochs_total = len(self.mne.boundary_times) - 1 + self.mne.epoch_tmins = np.full( + n_epochs_total, float(self.mne.inst.times[0]) + ) + self.mne.epoch_tmaxs = ( + self.mne.epoch_tmins + + np.diff(self.mne.boundary_times) + - 1.0 / sfreq + ) # initialize picks and projectors self._update_picks() @@ -328,14 +365,35 @@ def _make_butterfly_selections_dict(self): # MANAGE DATA # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # + def _get_epoch_ix_range(self): + """Return the first and last+1 epoch index currently in view. + + Both :meth:`_get_start_stop` and :meth:`_load_data` go through here so + the sample bounds and the concatenated data cannot disagree, which is + what keeps the shape assertions in :meth:`_update_data` meaningful when + epochs differ in duration. + """ + # subtract one sample from tstart before searchsorted, to make sure + # we land on the left side of the boundary time (avoid precision + # errors) + ix_start = int( + np.searchsorted( + self.mne.boundary_times, self.mne.t_start - self.mne.sampling_period + ) + ) + n_total = len(self.mne.boundary_times) - 1 + ix_start = min(ix_start, max(n_total - 1, 0)) + ix_stop = min(ix_start + self.mne.n_epochs, n_total) + return ix_start, ix_stop + def _get_start_stop(self): # update time start_sec = self.mne.t_start - self.mne.first_time if self.mne.is_epochs: - start, stop = np.round( - np.array([start_sec, start_sec + self.mne.duration]) - * self.mne.info["sfreq"] - ).astype(int) + # this agrees with _load_data by construction, not by arithmetic + ix_start, ix_stop = self._get_epoch_ix_range() + start = int(self.mne.boundary_samples[ix_start]) + stop = int(self.mne.boundary_samples[ix_stop]) else: # ensure our end time includes the last sample disp_duration = ( @@ -355,13 +413,7 @@ def _load_data(self, start=None, stop=None): else: return self.mne.inst[:, start:stop] else: - # subtract one sample from tstart before searchsorted, to make sure - # we land on the left side of the boundary time (avoid precision - # errors) - ix_start = np.searchsorted( - self.mne.boundary_times, self.mne.t_start - self.mne.sampling_period - ) - ix_stop = ix_start + self.mne.n_epochs + ix_start, ix_stop = self._get_epoch_ix_range() item = slice(ix_start, ix_stop) data = np.concatenate( self.mne.inst.get_data(item=item, copy=False), axis=-1 @@ -635,7 +687,11 @@ def _create_epoch_histogram(self): """Create peak-to-peak histogram of channel amplitudes.""" epochs = self.mne.inst data = OrderedDict() - ptp = np.ptp(epochs.get_data(copy=False), axis=2) + # per epoch, so that variable-duration epochs (a list of arrays) + # work too; peak-to-peak is a per-trial reduction either way + ptp = np.array( + [np.ptp(epoch, axis=-1) for epoch in epochs.get_data(copy=False)] + ) for ch_type in ("eeg", "mag", "grad"): if ch_type in epochs: data[ch_type] = ptp.T[self.mne.ch_types == ch_type].ravel() @@ -733,6 +789,27 @@ def _load_backend(backend_name): return backend +def _check_variable_duration_backend(): + """Raise unless the active browser backend can draw ragged epochs. + + Matplotlib always can. The Qt backend gained the ability in + mne-qt-browser 0.8, which announces it with a module-level flag, so an + older one declines here rather than drawing a wrong picture from a + boundary model it does not know about. + """ + backend_name = get_browser_backend() + if backend_name == "matplotlib": + return + module = _load_backend(backend_name) + if not getattr(module, "_SUPPORTS_VARIABLE_DURATION", False): + raise NotImplementedError( + f"Browsing variable-duration epochs is not implemented for the " + f"{backend_name} backend of this version, only for matplotlib. " + "Upgrade mne-qt-browser, or select matplotlib with " + 'mne.viz.set_browser_backend("matplotlib").' + ) + + def _get_browser(show, block, **kwargs): """Instantiate a new MNE browse-style figure.""" from .utils import _get_figsize_from_config diff --git a/mne/viz/_mpl_figure.py b/mne/viz/_mpl_figure.py index ff04b3395a6..3ad5989c80d 100644 --- a/mne/viz/_mpl_figure.py +++ b/mne/viz/_mpl_figure.py @@ -58,7 +58,7 @@ from ..defaults import DEFAULTS from ..fixes import _close_event from ..utils import Bunch, _click_ch_name, logger -from ._figure import BrowserBase +from ._figure import BrowserBase, _epoch_window from .utils import ( _BLIT_KWARGS, DraggableLine, @@ -489,7 +489,7 @@ def __init__(self, inst, figsize, ica=None, xlabel="Time (s)", **kwargs): epoch_nums = self.mne.inst.selection for ix, _ in enumerate(epoch_nums): start = self.mne.boundary_times[ix] - width = np.diff(self.mne.boundary_times[:2])[0] + width = self.mne.boundary_times[ix + 1] - start ax_hscroll.add_patch( Rectangle( (start, 0), @@ -791,12 +791,19 @@ def _keypress(self, event): old_t_start = self.mne.t_start direction = 1 if key.endswith("right") else -1 if self.mne.is_epochs: - denom = 1 if key.startswith("shift") else self.mne.n_epochs + # step whole epochs, since they need not share a duration + step = self.mne.n_epochs if key.startswith("shift") else 1 + ix_start, _ = self._get_epoch_ix_range() + self.mne.t_start, self.mne.duration = _epoch_window( + self.mne.boundary_times, + ix_start + direction * step, + self.mne.n_epochs, + ) else: denom = 1 if key.startswith("shift") else 4 - t_max = last_time - self.mne.duration - t_start = self.mne.t_start + direction * self.mne.duration / denom - self.mne.t_start = np.clip(t_start, self.mne.first_time, t_max) + t_max = last_time - self.mne.duration + t_start = self.mne.t_start + direction * self.mne.duration / denom + self.mne.t_start = np.clip(t_start, self.mne.first_time, t_max) if self.mne.t_start != old_t_start: self._update_hscroll() self._redraw(annotations=True, skip_hscroll=True) @@ -828,13 +835,16 @@ def _keypress(self, event): old_dur = self.mne.duration dur_delta = 1 if key == "end" else -1 if self.mne.is_epochs: + ix_start, _ = self._get_epoch_ix_range() # prevent from showing zero epochs, or more epochs than we have - self.mne.n_epochs = np.clip( - self.mne.n_epochs + dur_delta, 1, len(self.mne.inst) + self.mne.n_epochs = int( + np.clip(self.mne.n_epochs + dur_delta, 1, len(self.mne.inst)) + ) + # the epochs added or removed have their own durations + self.mne.t_start, new_dur = _epoch_window( + self.mne.boundary_times, ix_start, self.mne.n_epochs ) - # use the length of one epoch as duration change - min_dur = len(self.mne.inst.times) / self.mne.info["sfreq"] - new_dur = self.mne.duration + dur_delta * min_dur + min_dur = np.diff(self.mne.boundary_times).min() else: # never show fewer than 3 samples min_dur = 3 * np.diff(self.mne.inst.times[:2])[0] @@ -843,8 +853,10 @@ def _keypress(self, event): new_dur = self.mne.duration * dur_delta self.mne.duration = np.clip(new_dur, min_dur, last_time) if self.mne.duration != old_dur: - if self.mne.t_start + self.mne.duration > last_time: - self.mne.t_start = last_time - self.mne.duration + if not self.mne.is_epochs: + if self.mne.t_start + self.mne.duration > last_time: + self.mne.t_start = last_time - self.mne.duration + # (the epochs branch above already clamped t_start to a boundary) self._update_hscroll() self._redraw(annotations=True, skip_hscroll=True) elif key == "?": # help window @@ -998,7 +1010,9 @@ def _mouse_move(self, event): time = np.clip(time, self.mne.first_time, max_time) if self.mne.is_epochs: ix = np.searchsorted(self.mne.boundary_times[1:], time, side="right") - time = self.mne.boundary_times[ix] + time, self.mne.duration = _epoch_window( + self.mne.boundary_times, ix, self.mne.n_epochs + ) if self.mne.t_start != time: self.mne.t_start = time self._update_hscroll() @@ -1917,7 +1931,9 @@ def _check_update_hscroll_clicked(self, event): time = np.clip(time, self.mne.first_time, max_time) if self.mne.is_epochs: ix = np.searchsorted(self.mne.boundary_times[1:], time, side="right") - time = self.mne.boundary_times[ix] + time, self.mne.duration = _epoch_window( + self.mne.boundary_times, ix, self.mne.n_epochs + ) if self.mne.t_start != time: self.mne.t_start = time self._update_hscroll() @@ -2243,8 +2259,10 @@ def _draw_traces(self): # check for bad epochs time_range = (self.mne.times + self.mne.first_time)[[0, -1]] if self.mne.instance_type == "epochs": - epoch_ix = np.searchsorted(self.mne.boundary_times, time_range) - epoch_ix = np.arange(epoch_ix[0], epoch_ix[1]) + # ask the view directly: deriving this from the time range drops + # the last epoch whenever it holds a single sample, because that + # sample's time is its own left boundary + epoch_ix = np.arange(*self._get_epoch_ix_range()) epoch_nums = self.mne.inst.selection[epoch_ix[0] : epoch_ix[-1] + 1] (visible_bad_epoch_ix,) = np.isin(epoch_nums, self.mne.bad_epochs).nonzero() while len(self.mne.epoch_traces): @@ -2252,11 +2270,9 @@ def _draw_traces(self): # handle custom epoch colors (for autoreject integration) if self.mne.epoch_colors is None: # shape: n_traces × RGBA → n_traces × n_epochs × RGBA - custom_colors = np.tile( - ch_colors[:, None, :], (1, self.mne.n_epochs, 1) - ) + custom_colors = np.tile(ch_colors[:, None, :], (1, len(epoch_ix), 1)) else: - custom_colors = np.empty((len(self.mne.picks), self.mne.n_epochs, 4)) + custom_colors = np.empty((len(self.mne.picks), len(epoch_ix), 4)) for ii, _epoch_ix in enumerate(epoch_ix): this_colors = self.mne.epoch_colors[_epoch_ix] custom_colors[:, ii] = to_rgba_array( @@ -2316,7 +2332,11 @@ def _draw_traces(self): _starts = self.mne.boundary_times[epoch_ix][bool_ixs] _stops = self.mne.boundary_times[epoch_ix + 1][bool_ixs] for _start, _stop in zip(_starts, _stops): - _mask = np.logical_and(_start < this_times, this_times <= _stop) + # inclusive at both ends: an epoch owns its own first + # sample, and a one-sample epoch has nothing else + _mask = np.logical_and( + _start <= this_times, this_times <= _stop + ) mask = mask | _mask _times = np.ma.masked_array(this_times, mask=~mask) # always use the existing traces first @@ -2408,26 +2428,51 @@ def _recompute_epochs_vlines(self, xdata): # special case: changed view duration w/ "home" or "end" key # (no click event, hence no xdata) if xdata is None: - xdata = np.array(self.mne.vline.get_segments())[0, 0, 0] - # compute the (continuous) times for the lines on each epoch - epoch_dur = np.diff(self.mne.boundary_times[:2])[0] - rel_time = xdata % epoch_dur - abs_time = self.mne.times[0] - xs = np.arange(self.mne.n_epochs) * epoch_dur + abs_time + rel_time - segs = np.array(self.mne.vline.get_segments()) + segments = self.mne.vline.get_segments() + if not len(segments): # no visible epoch reaches that latency + return None + xdata = np.array(segments)[0, 0, 0] + # Work out which latency relative to its own event was clicked, then + # mark that same latency on every visible epoch. Epochs need not share a + # duration, so an epoch that never reaches this latency gets no line. + sfreq = self.mne.info["sfreq"] + boundary_times = self.mne.boundary_times + clicked_ix = int( + np.clip( + np.searchsorted(boundary_times[1:], xdata, side="right"), + 0, + len(boundary_times) - 2, + ) + ) + # clamp into the clicked epoch: a click in its last half sample would + # otherwise round up to one sample past its end, which no epoch holds, + # and every line would be dropped + n_samp = int( + round((boundary_times[clicked_ix + 1] - boundary_times[clicked_ix]) * sfreq) + ) + offset = int( + np.clip(round((xdata - boundary_times[clicked_ix]) * sfreq), 0, n_samp - 1) + ) + latency = self.mne.epoch_tmins[clicked_ix] + offset / sfreq + ix_start, ix_stop = self._get_epoch_ix_range() + xs = list() + for ix in range(ix_start, ix_stop): + tmin, tmax = self.mne.epoch_tmins[ix], self.mne.epoch_tmaxs[ix] + if tmin - 0.5 / sfreq <= latency <= tmax + 0.5 / sfreq: + xs.append(boundary_times[ix] + (latency - tmin)) + xs = np.array(xs, float) # recreate segs from scratch in case view duration changed # (i.e., handle case when n_segments != n_epochs) segs = np.tile([[0.0], [1.0]], (len(xs), 1, 2)) # y values - segs[..., 0] = np.tile(xs[:, None], 2) # x values + segs[..., 0] = np.tile(xs[:, None], 2) if len(xs) else segs[..., 0] self.mne.vline.set_segments(segs) - return rel_time + return latency def _show_vline(self, xdata): """Show the vertical line(s).""" if self.mne.is_epochs: - # convert xdata to be epoch-relative (for the text) - rel_time = self._recompute_epochs_vlines(xdata) - xdata = rel_time + self.mne.inst.times[0] + # the label shows the latency relative to each epoch's own event + xdata = self._recompute_epochs_vlines(xdata) else: self.mne.vline.set_xdata([xdata]) self.mne.vline_hscroll.set_xdata([xdata]) diff --git a/mne/viz/epochs.py b/mne/viz/epochs.py index 669205eaf4d..9f182b0f64e 100644 --- a/mne/viz/epochs.py +++ b/mne/viz/epochs.py @@ -964,11 +964,24 @@ def plot_epochs( raise ValueError(msg) # handle time dimension + # + # The browser lays the epochs end to end and draws a boundary between them, + # so its x axis is "concatenated seconds". Building that axis from the + # samples each epoch really has keeps it correct when durations vary, and + # reproduces the old uniform grid exactly when they do not. n_epochs = min(n_epochs, len(epochs)) - n_times = len(epochs) * len(epochs.times) - duration = n_epochs * len(epochs.times) / sfreq + lengths = np.array( + [epochs._n_times_per_epoch(ii) for ii in range(len(epochs))], int + ) # NB: this includes start and end of data: - boundary_times = np.arange(len(epochs) + 1) * len(epochs.times) / sfreq + boundary_samples = np.concatenate([[0], np.cumsum(lengths)]) + boundary_times = boundary_samples / sfreq + n_times = int(boundary_samples[-1]) + duration = boundary_times[n_epochs] - boundary_times[0] + # one value per epoch on both paths, so the browser never needs to ask + # whether it is looking at a scalar or an array + epoch_tmins = np.broadcast_to(np.asarray(epochs.tmin, float), (len(epochs),)).copy() + epoch_tmaxs = np.broadcast_to(np.asarray(epochs.tmax, float), (len(epochs),)).copy() # events _validate_type(events, (bool, np.ndarray), "events") @@ -980,14 +993,25 @@ def plot_epochs( events = epochs.events event_nums = events[:, 2] event_samps = events[:, 0] - epoch_n_samps = len(epochs.times) + # first sample of each epoch in the raw recording it was cut from + first_samps = epochs.events[:, 0] - np.round(-epoch_tmins * sfreq).astype(int) # handle overlapping epochs (each event may show up in multiple places) - boundaries = epochs.events[:, [0]] + np.array([-1, 1]) * epochs.time_as_index( - [0, epochs.tmax] - ) - in_bounds = np.logical_and( - boundaries[:, [0]] <= event_samps, event_samps < boundaries[:, [1]] - ) + if epochs._variable_duration: + # each epoch spans the samples it actually holds; both ends are + # inclusive, matching the fixed path's effective behaviour + last_samps = first_samps + lengths - 1 + in_bounds = np.logical_and( + first_samps[:, None] <= event_samps, + event_samps <= last_samps[:, None], + ) + else: + epoch_n_samps = len(epochs.times) + boundaries = epochs.events[:, [0]] + np.array( + [-1, 1] + ) * epochs.time_as_index([0, epochs.tmax]) + in_bounds = np.logical_and( + boundaries[:, [0]] <= event_samps, event_samps < boundaries[:, [1]] + ) event_ixs = [np.nonzero(a)[0] for a in in_bounds.T] warned = False event_times = list() @@ -1000,8 +1024,13 @@ def plot_epochs( "lines may be duplicated in the plot." ) warned = True - offsets = samp - relevant_epoch_events + epochs.time_as_index(0) - this_event_times = (_ixs * epoch_n_samps + offsets) / sfreq + if epochs._variable_duration: + # position inside each containing epoch, then along the strip + offsets = samp - first_samps[_ixs] + this_event_times = boundary_times[_ixs] + offsets / sfreq + else: + offsets = samp - relevant_epoch_events + epochs.time_as_index(0) + this_event_times = (_ixs * epoch_n_samps + offsets) / sfreq event_times.extend(this_event_times) event_numbers.extend([num] * len(_ixs)) event_nums = np.array(event_numbers) @@ -1059,6 +1088,9 @@ def plot_epochs( time_format="float", decim=decim, boundary_times=boundary_times, + boundary_samples=boundary_samples, + epoch_tmins=epoch_tmins, + epoch_tmaxs=epoch_tmaxs, # events event_id_rev=event_id_rev, event_color_dict=event_color_dict, @@ -1100,6 +1132,11 @@ def plot_epochs( figure_class=figure_class, ) + if epochs._variable_duration: + from ._figure import _check_variable_duration_backend + + _check_variable_duration_backend() + fig = _get_browser(show=show, block=block, **params) return fig diff --git a/mne/viz/tests/test_epochs.py b/mne/viz/tests/test_epochs.py index 874efa88e77..dd7800d271a 100644 --- a/mne/viz/tests/test_epochs.py +++ b/mne/viz/tests/test_epochs.py @@ -3,10 +3,12 @@ # Copyright the MNE-Python contributors. import platform +import warnings import matplotlib.pyplot as plt import numpy as np import pytest +from numpy.testing import assert_allclose, assert_array_equal from mne import Epochs, EpochsArray, create_info from mne.datasets import testing @@ -511,3 +513,218 @@ def test_plot_epochs_selection_butterfly(raw, browser_backend): epochs = Epochs(raw, events, tmin=0, tmax=0.5, preload=True, baseline=None) assert len(epochs) == 1 epochs.plot(group_by="selection", butterfly=True) + + +# -- variable-duration epochs ---------------------------------------------- +SFREQ_VAR = 100.0 +LENGTHS_VAR = (100, 250, 75, 180) # deliberately very different + + +def _variable_epochs(tmin=None): + """Return epochs whose trials have deliberately unequal lengths.""" + n = len(LENGTHS_VAR) + info = create_info(["a", "b", "c"], SFREQ_VAR, "eeg") + rng = np.random.default_rng(0) + data = [rng.standard_normal((3, length)) * 1e-6 for length in LENGTHS_VAR] + events = np.column_stack( + [np.arange(n) * 1000 + 500, np.zeros(n, int), np.ones(n, int)] + ) + tmin = np.zeros(n) if tmin is None else np.asarray(tmin, float) + return EpochsArray( + data, + info, + events=events, + tmin=tmin, + event_id={"x": 1}, + baseline=None, + verbose=False, + ) + + +def _boundaries(): + """Return the boundary times the browser should derive.""" + return np.concatenate([[0], np.cumsum(LENGTHS_VAR)]) / SFREQ_VAR + + +@pytest.fixture +def variable_epochs(browser_backend): + """Epochs of unequal duration sharing tmin=0.""" + from mne.viz._figure import _check_variable_duration_backend + + try: # skips qt until an mne-qt-browser that announces support is released + _check_variable_duration_backend() + except NotImplementedError as exc: + pytest.skip(str(exc)) + return _variable_epochs() + + +def test_plot_variable_duration_is_native(variable_epochs, browser_backend): + """Test that browsing ragged epochs neither warns nor pads.""" + + def _boom(*args, **kwargs): + raise AssertionError("plot() fell back to as_fixed() instead of browsing") + + variable_epochs.as_fixed = _boom + with warnings.catch_warnings(): + warnings.simplefilter("error") # any warning fails the test + fig = variable_epochs.plot(n_epochs=2) + assert not np.isnan(fig.mne.data).any() + + +def test_plot_variable_duration_boundaries(variable_epochs, browser_backend): + """Test that the browser lays epochs end to end at their true lengths.""" + fig = variable_epochs.plot(n_epochs=2) + assert_allclose(fig.mne.boundary_times, _boundaries()) + assert_array_equal(fig.mne.boundary_samples, np.r_[0, np.cumsum(LENGTHS_VAR)]) + # the concatenated axis holds every sample once, with nothing invented + assert fig.mne.n_times == sum(LENGTHS_VAR) + assert fig.mne.n_times != len(LENGTHS_VAR) * max(LENGTHS_VAR) + + +def test_plot_variable_duration_window_spans_whole_epochs( + variable_epochs, browser_backend +): + """Test that n_epochs means epochs, not a representative duration.""" + boundaries = _boundaries() + for n_epochs in (1, 2, 3, 4): + fig = variable_epochs.plot(n_epochs=n_epochs) + assert fig.mne.duration == pytest.approx(boundaries[n_epochs]) + assert fig.mne.data.shape[-1] == sum(LENGTHS_VAR[:n_epochs]) + + +def test_plot_variable_duration_data_is_unpadded(variable_epochs, browser_backend): + """Test that a view holds exactly the source samples, in order.""" + fig = variable_epochs.plot(n_epochs=2) + source = variable_epochs.get_data() + for keys in ([], ["right"], ["right", "right"]): + for key in keys: + fig._fake_keypress(key) + start, stop = fig._get_start_stop() + data, times = fig._load_data(start, stop) + ix_start, ix_stop = fig._get_epoch_ix_range() + want = np.concatenate(source[ix_start:ix_stop], axis=-1) + # the raw window is the source samples, in order, with nothing added + assert_array_equal(data, want) + assert not np.isnan(data).any() + # and the sample bounds agree with it, which is what lets the shape + # assertions in _update_data mean something for ragged windows + assert data.shape[-1] == stop - start + assert len(times) == stop - start + + +def test_plot_variable_duration_navigation(variable_epochs, browser_backend): + """Test that arrow keys move by epochs and land on real boundaries.""" + boundaries = _boundaries() + fig = variable_epochs.plot(n_epochs=2) + assert fig.mne.t_start == pytest.approx(boundaries[0]) + + fig._fake_keypress("right") # one epoch + assert fig.mne.t_start == pytest.approx(boundaries[1]) + assert fig.mne.duration == pytest.approx(boundaries[3] - boundaries[1]) + + fig._fake_keypress("right") # clamped: two epochs must stay visible + assert fig.mne.t_start == pytest.approx(boundaries[2]) + assert fig.mne.duration == pytest.approx(boundaries[4] - boundaries[2]) + + fig._fake_keypress("left") + assert fig.mne.t_start == pytest.approx(boundaries[1]) + + # shift moves a whole window + fig._fake_keypress("shift+left") + assert fig.mne.t_start == pytest.approx(boundaries[0]) + fig._fake_keypress("shift+right") + assert fig.mne.t_start == pytest.approx(boundaries[2]) + + +def test_plot_variable_duration_home_end(variable_epochs, browser_backend): + """Test that home/end change the epoch count and recompute the duration.""" + boundaries = _boundaries() + fig = variable_epochs.plot(n_epochs=2) + assert fig.mne.duration == pytest.approx(boundaries[2]) + + fig._fake_keypress("end") # show one more epoch + assert fig.mne.n_epochs == 3 + assert fig.mne.duration == pytest.approx(boundaries[3]) + # the epoch added is 75 samples, not a repeat of the first one + assert fig.mne.duration != pytest.approx(boundaries[2] * 3 / 2) + + fig._fake_keypress("home") + assert fig.mne.n_epochs == 2 + assert fig.mne.duration == pytest.approx(boundaries[2]) + + +def test_plot_variable_duration_hscroll_patches(variable_epochs, browser_backend): + """Test that the scrollbar draws each epoch at its own width.""" + if browser_backend.name != "matplotlib": + pytest.skip("scrollbar patches are matplotlib-specific") + fig = variable_epochs.plot(n_epochs=2) + widths = [p.get_width() for p in fig.mne.ax_hscroll.patches[: len(LENGTHS_VAR)]] + assert_allclose(widths, np.diff(_boundaries())) + assert len(set(np.round(widths, 6))) == len(LENGTHS_VAR) # all different + + +def test_plot_variable_duration_bad_epoch(variable_epochs, browser_backend): + """Test that a click finds the right epoch when the widths differ.""" + if browser_backend.name != "matplotlib": + pytest.skip("epoch marking by click is matplotlib-specific") + boundaries = _boundaries() + fig = variable_epochs.plot(n_epochs=4) + y = fig.mne.traces[0].get_ydata()[0] + # click inside epoch 2, which starts well past twice the first epoch's width + x = (boundaries[2] + boundaries[3]) / 2 + fig._fake_click((x, y), xform="data") + assert list(fig.mne.bad_epochs) == [variable_epochs.selection[2]] + fig._fake_click((x, y), xform="data") # unmark + assert list(fig.mne.bad_epochs) == [] + + +def test_plot_variable_duration_vline_latency(variable_epochs, browser_backend): + """Test that vlines mark a latency, skipping epochs that never reach it.""" + if browser_backend.name != "matplotlib": + pytest.skip("vline segments are matplotlib-specific") + boundaries = _boundaries() + fig = variable_epochs.plot(n_epochs=4) + + # 1.5 s exists only in the 250- and 180-sample epochs + latency = 1.5 + fig._fake_click((boundaries[1] + latency, 0.5), xform="data") + xs = np.sort(np.array(fig.mne.vline.get_segments())[:, 0, 0]) + assert_allclose(xs, [boundaries[1] + latency, boundaries[3] + latency]) + + # 0.5 s exists in every epoch + fig._fake_click((boundaries[0] + 0.5, 0.5), xform="data") + xs = np.sort(np.array(fig.mne.vline.get_segments())[:, 0, 0]) + assert_allclose(xs, boundaries[:4] + 0.5) + + +def test_plot_variable_duration_events(browser_backend): + """Test that events map into the strip using each epoch's own window.""" + if browser_backend.name != "matplotlib": + pytest.skip("event line segments are matplotlib-specific") + # unequal tmin as well as unequal duration + tmin = np.array([0.0, -0.2, 0.0, -0.5]) + epochs = _variable_epochs(tmin=tmin) + boundaries = _boundaries() + first_samps = epochs.events[:, 0] - np.round(-tmin * SFREQ_VAR).astype(int) + + fig = epochs.plot(n_epochs=4, events=True) + lines, _ = _get_event_lines_and_texts(fig) + got = np.sort(np.array(lines)[:, 0, 0]) + # each defining event sits at its own offset inside its own epoch + want = np.sort( + boundaries[: len(LENGTHS_VAR)] + (epochs.events[:, 0] - first_samps) / SFREQ_VAR + ) + assert_allclose(got, want) + + +def test_plot_variable_duration_refuses_old_backends(monkeypatch): + """Test that a backend without the boundary model declines, not fails.""" + import mne.viz._figure + + class _OldBackend: # an mne-qt-browser that predates the boundary model + pass + + monkeypatch.setattr(mne.viz._figure, "get_browser_backend", lambda: "qt") + monkeypatch.setattr(mne.viz._figure, "_load_backend", lambda name: _OldBackend()) + with pytest.raises(NotImplementedError, match="not implemented for the qt"): + _variable_epochs().plot() diff --git a/mne/viz/utils.py b/mne/viz/utils.py index 80a5b146ec8..b36440fe0f9 100644 --- a/mne/viz/utils.py +++ b/mne/viz/utils.py @@ -1401,7 +1401,8 @@ def _compute_scalings(scalings, inst, remove_dc=False, duration=10): data = inst._read_segment(smin, smax) elif isinstance(inst, BaseEpochs): # Load a random subset of epochs up to 100mb in size - n_epochs = 1e8 // (len(inst.ch_names) * len(inst.times) * 8) + longest = max(inst._n_times_per_epoch(ii) for ii in range(len(inst))) + n_epochs = 1e8 // (len(inst.ch_names) * longest * 8) n_epochs = int(np.clip(n_epochs, 1, len(inst))) ixs_epochs = np.random.default_rng(0).choice( len(inst), n_epochs, replace=False @@ -1410,7 +1411,11 @@ def _compute_scalings(scalings, inst, remove_dc=False, duration=10): else: data = inst._data if isinstance(inst, BaseEpochs): - data = inst._data.swapaxes(0, 1).reshape([len(inst.ch_names), -1]) + if isinstance(inst._data, list): + # variable-duration epochs: one array per epoch, already channels-first + data = np.concatenate(inst._data, axis=-1) + else: + data = inst._data.swapaxes(0, 1).reshape([len(inst.ch_names), -1]) # Iterate through ch types and update scaling if ' auto' for key, value in scalings.items(): if key not in ch_types or value != "auto": diff --git a/tutorials/epochs/70_variable_duration_epochs.py b/tutorials/epochs/70_variable_duration_epochs.py new file mode 100644 index 00000000000..52afe8bd2e5 --- /dev/null +++ b/tutorials/epochs/70_variable_duration_epochs.py @@ -0,0 +1,304 @@ +""" +.. _tut-variable-duration-epochs: + +Epochs whose trials have different durations +============================================ + +Most epoching starts from an event and takes the same window around every one of +them, which gives a rectangular ``(n_epochs, n_channels, n_times)`` array and one +time axis shared by every trial. Some experiments do not fit that shape. A gait +cycle, a spoken word, a reaching movement and a sleep stage all last as long as +they last, and the duration is often the thing being studied. + +The usual way to handle this is to pick a fixed window and accept the +consequences: a window long enough for the longest trial pads the short ones, and +a window short enough for the shortest one truncates the rest. This tutorial +shows the other option, keeping each trial at the length it actually had, and +what the resulting object can and cannot do. + +We use the Sleep Physionet data, where the hypnogram annotations mark sleep stage +bouts and each bout carries its own duration. +""" + +# Authors: The MNE-Python contributors. +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +# %% + +import matplotlib.pyplot as plt +import numpy as np + +import mne +from mne.datasets.sleep_physionet.age import fetch_data + +psg_file, hypnogram_file = fetch_data(subjects=[0], recording=[1])[0] + +raw = mne.io.read_raw_edf( + psg_file, + stim_channel=False, + preload=True, + verbose="error", # ignore issues with stored filter settings +) +raw.pick(["EEG Fpz-Cz", "EEG Pz-Oz"]) + +annotations = mne.read_annotations(hypnogram_file) +raw.set_annotations(annotations, emit_warning=False) + +# %% +# The annotations already carry durations +# --------------------------------------- +# +# Each hypnogram entry marks one bout of a sleep stage, and +# :class:`~mne.Annotations` stores its ``duration`` alongside its ``onset``. That +# duration is not a constant. + +stages = { + "Sleep stage 1": 1, + "Sleep stage 2": 2, + "Sleep stage 3": 3, + "Sleep stage 4": 3, # stages 3 and 4 are conventionally merged + "Sleep stage R": 4, +} +event_id = {"N1": 1, "N2": 2, "N3/4": 3, "REM": 4} + +onset = annotations.onset +duration = annotations.duration +description = np.array(annotations.description) + +keep = np.array([desc in stages for desc in description]) +# a five minute cap keeps the padded array in the last section small; nothing +# about the container requires it +keep &= duration <= 300.0 + +onset, duration = onset[keep], duration[keep] +description = description[keep] + +print(f"{len(duration)} bouts, {duration.min():.0f} to {duration.max():.0f} s") +print(f"median {np.median(duration):.0f} s") + +# %% +# Building epochs that keep those durations +# ----------------------------------------- +# +# ``tmin`` and ``tmax`` accept one value per event as well as a single number. +# Here every bout starts at its own onset, so ``tmin`` is zero throughout and +# ``tmax`` is the bout's own length. The last sample is included, which is why +# ``tmax`` is one sample short of the full duration. + +sfreq = raw.info["sfreq"] +events = np.column_stack( + [ + np.round((onset - raw.first_time) * sfreq).astype(int), + np.zeros(len(onset), int), + np.array([stages[desc] for desc in description]), + ] +) + +epochs = mne.Epochs( + raw, + events, + event_id, + tmin=np.zeros(len(events)), + tmax=duration - 1.0 / sfreq, + baseline=None, + preload=True, +) +print(epochs) + +# %% +# The object reports that its trials are not all the same length, and the +# durations it holds are the ones the annotations described. + +print(f"variable_duration: {epochs.variable_duration}") +print(f"durations: {epochs.durations.min():.0f} to {epochs.durations.max():.0f} s") + +# %% +# Because bounds that carry no variation collapse back to a single value, this +# only changes behaviour when the durations really do differ. Passing equal +# bounds gives an ordinary fixed-duration ``Epochs``. + +n2_events = events[events[:, 2] == event_id["N2"]][:5] +fixed_bounds = mne.Epochs( + raw, + n2_events, + {"N2": event_id["N2"]}, + tmin=np.zeros(len(n2_events)), + tmax=np.full(len(n2_events), 29.99), + baseline=None, + preload=True, + verbose=False, +) +print(f"equal bounds -> variable_duration: {fixed_bounds.variable_duration}") + +# %% +# Getting the data out +# -------------------- +# +# There is no rectangular array to return, so :meth:`~mne.Epochs.get_data` gives +# a list with one ``(n_channels, n_times)`` array per epoch. Nothing is padded +# and nothing is cut: each array holds exactly the samples that the bout covered +# in the continuous recording. + +data = epochs.get_data() +print(f"{len(data)} arrays, first four shapes {[d.shape for d in data[:4]]}") + +lengths = np.array([d.shape[-1] for d in data]) +print(f"total samples held: {lengths.sum()}") +print(f"a rectangular array would hold: {lengths.max() * len(lengths)}") + +# %% +# For the same reason there is no single ``times`` attribute. Each epoch has its +# own time axis, which :meth:`~mne.Epochs.get_times` returns. + +for idx in (0, 1): + t = epochs.get_times(idx) + print(f"epoch {idx}: {len(t)} samples, {t[0]:.2f} to {t[-1]:.2f} s") + +# %% +# Asking for ``epochs.times`` raises rather than inventing an axis. Returning the +# longest epoch's axis would make ``len(epochs.times)`` disagree with the data +# for every other epoch while looking perfectly normal. + +try: + epochs.times +except RuntimeError as err: + print(f"RuntimeError: {err}") + +# %% +# Operations that do not touch the time axis +# ------------------------------------------ +# +# Selecting epochs, selecting channels and dropping epochs all work as usual, +# because none of them care how long each trial is. The per-epoch bounds travel +# with the epochs they belong to. + +n2 = epochs["N2"] +print( + f"epochs['N2']: {len(n2)} epochs, " + f"{n2.durations.min():.0f} to {n2.durations.max():.0f} s" +) + +first_ten = epochs[:10] +print(f"epochs[:10]: durations {first_ten.durations.round(0)}") + +one_channel = epochs.copy().pick(["EEG Pz-Oz"]) +print( + f"after pick: {one_channel.ch_names}, durations unchanged: " + f"{np.array_equal(one_channel.durations, epochs.durations)}" +) + +# %% +# Browsing them +# ------------- +# +# :meth:`~mne.Epochs.plot` shows each bout at the length it really has. The +# browser lays the variable-length blocks end to end and rules a line between +# them, so the vertical boundaries are unevenly spaced: a 30 second bout takes a +# fifth of the width of a 150 second one. Nothing is padded or truncated to make +# the picture rectangular, and :meth:`~mne.Epochs.as_fixed` is not involved. +# +# Pick a handful of bouts with genuinely different lengths, taking the first +# occurrence of each distinct duration rather than trusting the first few epochs +# to differ. + +_, first_of_each = np.unique(epochs.durations, return_index=True) +browse_idx = np.sort(first_of_each[:5]) +browse_epochs = epochs[browse_idx] +print(f"browsing durations: {browse_epochs.durations.round(0)} s") + +# the browser's time axis is the real samples, laid end to end +n_browser_samples = sum( + len(browse_epochs.get_times(ii)) for ii in range(len(browse_epochs)) +) +print(f"{n_browser_samples} samples in total, none of them padding") + +# %% +# Browsing variable-duration epochs currently needs the Matplotlib backend; the +# PyQtGraph one does not handle ragged epochs yet. + +with mne.viz.use_browser_backend("matplotlib"): + browse_epochs.plot(n_epochs=len(browse_epochs), picks="eeg") + +# %% +# Operations that need one time axis +# ---------------------------------- +# +# Averaging is the clearest case. :class:`~mne.Evoked` holds one array and one +# ``nave``, and there is no honest way to fill either when the trials stop at +# different times. Rather than pad quietly, the reduction refuses and says what +# it would need. + +try: + epochs.average() +except NotImplementedError as err: + print(f"NotImplementedError: {err}") + +# %% +# Making the padding explicit +# --------------------------- +# +# When a rectangular array is genuinely what you want, +# :meth:`~mne.Epochs.as_fixed` produces one. It returns the padded +# :class:`~mne.EpochsArray` together with the number of epochs contributing at +# each sample, so the cost of the padding is visible rather than implied. + +padded, n_contributing = epochs.as_fixed() +print(f"padded shape: {padded.get_data().shape}") +print( + f"contributing: {n_contributing.max()} at the start, " + f"{n_contributing.min()} at the end" +) + +held = lengths.sum() * len(epochs.ch_names) +allocated = padded.get_data().size +print(f"padding waste: {100 * (1 - held / allocated):.1f}%") + +# %% +# That second return value is the point of the method. Plotted against time it +# shows how quickly the epochs stop contributing, which is exactly the +# information an averaged :class:`~mne.Evoked` cannot carry. + +fig, ax = plt.subplots(figsize=(8, 4), layout="constrained") +times = padded.times +ax.fill_between(times, n_contributing, step="post", alpha=0.25) +ax.plot(times, n_contributing, drawstyle="steps-post") + +half = len(epochs) / 2 +crossing = times[np.argmax(n_contributing < half)] +ax.axhline(half, color="0.4", ls=":", lw=1) +ax.axvline(crossing, color="0.4", ls=":", lw=1) +ax.annotate( + f"half the epochs have ended by {crossing:.0f} s", + xy=(crossing, half), + xytext=(crossing + 20, len(epochs) * 0.7), + arrowprops=dict(arrowstyle="->", color="0.4"), +) + +ax.set( + xlabel="Time (s)", + ylabel="Epochs contributing", + title="How many sleep-stage bouts are still running", + xlim=(0, times[-1]), + ylim=(0, len(epochs) * 1.05), +) + +# %% +# Reading the figure from left to right: every bout contributes at the start, +# and by the end a single long bout is holding up the whole window. An average +# over this padded array would combine all of them at ``t = 0`` and one of them +# at the right-hand edge, while reporting one ``nave`` for the lot. Keeping the +# count alongside the data is what makes that visible. +# +# When fixed windows are the right choice +# --------------------------------------- +# +# None of this argues against fixed-length epochs. Sleep staging is a good +# example of when they are correct: :ref:`tut-sleep-stage-classif` classifies 30 +# second windows, so it passes ``chunk_duration=30.`` to +# :func:`mne.events_from_annotations` and deliberately turns each bout into a +# series of equal windows. That is the right representation when the window is +# the unit of analysis. +# +# Variable-duration epochs are for the other case, when the bout itself is the +# unit and its length is part of what is being measured.