epochs.py

"""Tools for working with epoched data."""

# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.

import json
import operator
import os.path as op
from collections import Counter
from copy import deepcopy
from functools import partial
from inspect import getfullargspec
from pathlib import Path

import numpy as np
from scipy.interpolate import interp1d

from ._fiff.constants import FIFF
from ._fiff.meas_info import (
    ContainsMixin,
    SetChannelsMixin,
    _ensure_infos_match,
    read_meas_info,
    write_meas_info,
)
from ._fiff.open import _get_next_fname, fiff_open
from ._fiff.pick import (
    _DATA_CH_TYPES_SPLIT,
    _pick_data_channels,
    _picks_to_idx,
    channel_indices_by_type,
    channel_type,
    pick_channels,
    pick_info,
)
from ._fiff.proj import ProjMixin, setup_proj
from ._fiff.tag import _read_tag_header, read_tag
from ._fiff.tree import dir_tree_find
from ._fiff.utils import _make_split_fnames
from ._fiff.write import (
    _NEXT_FILE_BUFFER,
    INT32_MAX,
    _get_split_size,
    end_block,
    start_and_end_file,
    start_block,
    write_complex_double_matrix,
    write_complex_float_matrix,
    write_double_matrix,
    write_float,
    write_float_matrix,
    write_id,
    write_int,
    write_string,
)
from .annotations import (
    EpochAnnotationsMixin,
    _read_annotations_fif,
    _write_annotations,
    events_from_annotations,
)
from .baseline import _check_baseline, _log_rescale, rescale
from .bem import _check_origin
from .channels.channels import InterpolationMixin, ReferenceMixin, UpdateChannelsMixin
from .event import _read_events_fif, make_fixed_length_events, match_event_names
from .evoked import EvokedArray
from .filter import FilterMixin, _check_fun, detrend
from .fixes import rng_uniform
from .html_templates import _get_html_template
from .parallel import parallel_func
from .time_frequency.spectrum import EpochsSpectrum, SpectrumMixin, _validate_method
from .time_frequency.tfr import AverageTFR, EpochsTFR
from .utils import (
    ExtendedTimeMixin,
    GetEpochsMixin,
    SizeMixin,
    _build_data_frame,
    _check_combine,
    _check_event_id,
    _check_fname,
    _check_option,
    _check_pandas_index_arguments,
    _check_pandas_installed,
    _check_preload,
    _check_time_format,
    _convert_times,
    _ensure_events,
    _gen_events,
    _on_missing,
    _path_like,
    _pl,
    _prepare_read_metadata,
    _prepare_write_metadata,
    _scale_dataframe_data,
    _validate_type,
    check_fname,
    check_random_state,
    copy_function_doc_to_method_doc,
    logger,
    object_size,
    repr_html,
    sizeof_fmt,
    verbose,
    warn,
)
from .utils.docs import fill_doc
from .viz import plot_drop_log, plot_epochs, plot_epochs_image, plot_topo_image_epochs


def _pack_reject_params(epochs):
    reject_params = dict()
    for key in ("reject", "flat", "reject_tmin", "reject_tmax"):
        val = getattr(epochs, key, None)
        if val is not None:
            reject_params[key] = val
    return reject_params


def _save_split(epochs, split_fnames, part_idx, n_parts, fmt, overwrite):
    """Split epochs.

    Anything new added to this function also needs to be added to
    BaseEpochs.save to account for new file sizes.
    """
    # insert index in filename
    this_fname = split_fnames[part_idx]
    _check_fname(this_fname, overwrite=overwrite)

    next_fname, next_idx = None, None
    if part_idx < n_parts - 1:
        next_idx = part_idx + 1
        next_fname = split_fnames[next_idx]

    with start_and_end_file(this_fname) as fid:
        _save_part(fid, epochs, fmt, n_parts, next_fname, next_idx)


def _save_part(fid, epochs, fmt, n_parts, next_fname, next_idx):
    info = epochs.info
    meas_id = info["meas_id"]

    start_block(fid, FIFF.FIFFB_MEAS)
    write_id(fid, FIFF.FIFF_BLOCK_ID)
    if info["meas_id"] is not None:
        write_id(fid, FIFF.FIFF_PARENT_BLOCK_ID, info["meas_id"])

    # Write measurement info
    write_meas_info(fid, info)

    # One or more evoked data sets
    start_block(fid, FIFF.FIFFB_PROCESSED_DATA)
    start_block(fid, FIFF.FIFFB_MNE_EPOCHS)

    # write events out after getting data to ensure bad events are dropped
    data = epochs.get_data(copy=False)

    _check_option("fmt", fmt, ["single", "double"])

    if np.iscomplexobj(data):
        if fmt == "single":
            write_function = write_complex_float_matrix
        elif fmt == "double":
            write_function = write_complex_double_matrix
    else:
        if fmt == "single":
            write_function = write_float_matrix
        elif fmt == "double":
            write_function = write_double_matrix

    # Epoch annotations are written if there are any
    annotations = getattr(epochs, "annotations", [])
    if annotations is not None and len(annotations):
        _write_annotations(fid, annotations)

    # write Epoch event windows
    start_block(fid, FIFF.FIFFB_MNE_EVENTS)
    write_int(fid, FIFF.FIFF_MNE_EVENT_LIST, epochs.events.T)
    write_string(fid, FIFF.FIFF_DESCRIPTION, _event_id_string(epochs.event_id))
    end_block(fid, FIFF.FIFFB_MNE_EVENTS)

    # Metadata
    if epochs.metadata is not None:
        start_block(fid, FIFF.FIFFB_MNE_METADATA)
        metadata = _prepare_write_metadata(epochs.metadata)
        write_string(fid, FIFF.FIFF_DESCRIPTION, metadata)
        end_block(fid, FIFF.FIFFB_MNE_METADATA)

    # First and last sample
    first = int(round(epochs.tmin * info["sfreq"]))  # round just to be safe
    last = first + len(epochs.times) - 1
    write_int(fid, FIFF.FIFF_FIRST_SAMPLE, first)
    write_int(fid, FIFF.FIFF_LAST_SAMPLE, last)

    # write raw original sampling rate
    write_float(fid, FIFF.FIFF_MNE_EPOCHS_RAW_SFREQ, epochs._raw_sfreq)

    # save baseline
    if epochs.baseline is not None:
        bmin, bmax = epochs.baseline
        write_float(fid, FIFF.FIFF_MNE_BASELINE_MIN, bmin)
        write_float(fid, FIFF.FIFF_MNE_BASELINE_MAX, bmax)

    # The epochs itself
    decal = np.empty(info["nchan"])
    for k in range(info["nchan"]):
        decal[k] = 1.0 / (info["chs"][k]["cal"] * info["chs"][k].get("scale", 1.0))

    data *= decal[np.newaxis, :, np.newaxis]

    write_function(fid, FIFF.FIFF_EPOCH, data)

    # undo modifications to data
    data /= decal[np.newaxis, :, np.newaxis]

    write_string(fid, FIFF.FIFF_MNE_EPOCHS_DROP_LOG, json.dumps(epochs.drop_log))

    reject_params = _pack_reject_params(epochs)
    if reject_params:
        write_string(fid, FIFF.FIFF_MNE_EPOCHS_REJECT_FLAT, json.dumps(reject_params))

    write_int(fid, FIFF.FIFF_MNE_EPOCHS_SELECTION, epochs.selection)

    # And now write the next file info in case epochs are split on disk
    if next_fname is not None and n_parts > 1:
        start_block(fid, FIFF.FIFFB_REF)
        write_int(fid, FIFF.FIFF_REF_ROLE, FIFF.FIFFV_ROLE_NEXT_FILE)
        write_string(fid, FIFF.FIFF_REF_FILE_NAME, op.basename(next_fname))
        if meas_id is not None:
            write_id(fid, FIFF.FIFF_REF_FILE_ID, meas_id)
        write_int(fid, FIFF.FIFF_REF_FILE_NUM, next_idx)
        end_block(fid, FIFF.FIFFB_REF)

    end_block(fid, FIFF.FIFFB_MNE_EPOCHS)
    end_block(fid, FIFF.FIFFB_PROCESSED_DATA)
    end_block(fid, FIFF.FIFFB_MEAS)


def _event_id_string(event_id):
    return ";".join([k + ":" + str(v) for k, v in event_id.items()])


def _merge_events(events, event_id, selection):
    """Merge repeated events."""
    event_id = event_id.copy()
    new_events = events.copy()
    event_idxs_to_delete = list()
    unique_events, counts = np.unique(events[:, 0], return_counts=True)
    for ev in unique_events[counts > 1]:
        # indices at which the non-unique events happened
        idxs = (events[:, 0] == ev).nonzero()[0]

        # Figure out new value for events[:, 1]. Set to 0, if mixed vals exist
        unique_priors = np.unique(events[idxs, 1])
        new_prior = unique_priors[0] if len(unique_priors) == 1 else 0

        # If duplicate time samples have same event val, "merge" == "drop"
        # and no new event_id key will be created
        ev_vals = np.unique(events[idxs, 2])
        if len(ev_vals) <= 1:
            new_event_val = ev_vals[0]

        # Else, make a new event_id for the merged event
        else:
            # Find all event_id keys involved in duplicated events. These
            # keys will be merged to become a new entry in "event_id"
            event_id_keys = list(event_id.keys())
            event_id_vals = list(event_id.values())
            new_key_comps = [
                event_id_keys[event_id_vals.index(value)] for value in ev_vals
            ]

            # Check if we already have an entry for merged keys of duplicate
            # events ... if yes, reuse it
            for key in event_id:
                if set(key.split("/")) == set(new_key_comps):
                    new_event_val = event_id[key]
                    break

            # Else, find an unused value for the new key and make an entry into
            # the event_id dict
            else:
                ev_vals = np.unique(
                    np.concatenate(
                        (list(event_id.values()), events[:, 1:].flatten()), axis=0
                    )
                )
                if ev_vals[0] > 1:
                    new_event_val = 1
                else:
                    diffs = np.diff(ev_vals)
                    idx = np.where(diffs > 1)[0]
                    idx = -1 if len(idx) == 0 else idx[0]
                    new_event_val = ev_vals[idx] + 1

                new_event_id_key = "/".join(sorted(new_key_comps))
                event_id[new_event_id_key] = int(new_event_val)

        # Replace duplicate event times with merged event and remember which
        # duplicate indices to delete later
        new_events[idxs[0], 1] = new_prior
        new_events[idxs[0], 2] = new_event_val
        event_idxs_to_delete.extend(idxs[1:])

    # Delete duplicate event idxs
    new_events = np.delete(new_events, event_idxs_to_delete, 0)
    new_selection = np.delete(selection, event_idxs_to_delete, 0)

    return new_events, event_id, new_selection


def _handle_event_repeated(events, event_id, event_repeated, selection, drop_log):
    """Handle repeated events.

    Note that drop_log will be modified inplace
    """
    assert len(events) == len(selection)
    selection = np.asarray(selection)

    unique_events, u_ev_idxs = np.unique(events[:, 0], return_index=True)

    # Return early if no duplicates
    if len(unique_events) == len(events):
        return events, event_id, selection, drop_log

    # Else, we have duplicates. Triage ...
    _check_option("event_repeated", event_repeated, ["error", "drop", "merge"])
    drop_log = list(drop_log)
    if event_repeated == "error":
        raise RuntimeError(
            "Event time samples were not unique. Consider "
            'setting the `event_repeated` parameter."'
        )

    elif event_repeated == "drop":
        logger.info(
            "Multiple event values for single event times found. "
            "Keeping the first occurrence and dropping all others."
        )
        new_events = events[u_ev_idxs]
        new_selection = selection[u_ev_idxs]
        drop_ev_idxs = np.setdiff1d(selection, new_selection)
        for idx in drop_ev_idxs:
            drop_log[idx] = drop_log[idx] + ("DROP DUPLICATE",)
        selection = new_selection
    elif event_repeated == "merge":
        logger.info(
            "Multiple event values for single event times found. "
            "Creating new event value to reflect simultaneous events."
        )
        new_events, event_id, new_selection = _merge_events(events, event_id, selection)
        drop_ev_idxs = np.setdiff1d(selection, new_selection)
        for idx in drop_ev_idxs:
            drop_log[idx] = drop_log[idx] + ("MERGE DUPLICATE",)
        selection = new_selection
    drop_log = tuple(drop_log)

    # Remove obsolete kv-pairs from event_id after handling
    keys = new_events[:, 1:].flatten()
    event_id = {k: v for k, v in event_id.items() if v in keys}

    return new_events, event_id, selection, drop_log


@fill_doc
class BaseEpochs(
    ProjMixin,
    ContainsMixin,
    UpdateChannelsMixin,
    ReferenceMixin,
    SetChannelsMixin,
    InterpolationMixin,
    FilterMixin,
    ExtendedTimeMixin,
    SizeMixin,
    GetEpochsMixin,
    EpochAnnotationsMixin,
    SpectrumMixin,
):
    """Abstract base class for `~mne.Epochs`-type classes.

    .. note::
        This class should not be instantiated directly via
        ``mne.BaseEpochs(...)``. Instead, use one of the functions listed in
        the See Also section below.

    Parameters
    ----------
    %(info_not_none)s
    data : ndarray | None
        If ``None``, data will be read from the Raw object. If ndarray, must be
        of shape (n_epochs, n_channels, n_times).
    %(events_epochs)s
    %(event_id)s
    %(epochs_tmin_tmax)s
    %(baseline_epochs)s
        Defaults to ``(None, 0)``, i.e. beginning of the the data until
        time point zero.
    %(raw_epochs)s
    %(picks_all)s
    %(reject_epochs)s
    %(flat)s
    %(decim)s
    %(epochs_reject_tmin_tmax)s
    %(detrend_epochs)s
    %(proj_epochs)s
    %(on_missing_epochs)s
    preload_at_end : bool
        %(epochs_preload)s
    %(selection)s

        .. versionadded:: 0.16
    %(drop_log)s
    filename : Path | None
        The filename (if the epochs are read from disk).
    %(metadata_epochs)s

        .. versionadded:: 0.16
    %(event_repeated_epochs)s
    %(raw_sfreq)s
    annotations : instance of mne.Annotations | None
        Annotations to set.
    %(verbose)s

    See Also
    --------
    Epochs
    EpochsArray
    make_fixed_length_epochs

    Notes
    -----
    The ``BaseEpochs`` class is public to allow for stable type-checking in
    user code (i.e., ``isinstance(my_epochs, BaseEpochs)``) but should not be
    used as a constructor for Epochs objects (use instead :class:`mne.Epochs`).
    """

    @verbose
    def __init__(
        self,
        info,
        data,
        events,
        event_id=None,
        tmin=-0.2,
        tmax=0.5,
        baseline=(None, 0),
        raw=None,
        picks=None,
        reject=None,
        flat=None,
        decim=1,
        reject_tmin=None,
        reject_tmax=None,
        detrend=None,
        proj=True,
        on_missing="raise",
        preload_at_end=False,
        selection=None,
        drop_log=None,
        filename=None,
        metadata=None,
        event_repeated="error",
        *,
        raw_sfreq=None,
        annotations=None,
        verbose=None,
    ):
        if events is not None:  # RtEpochs can have events=None
            events = _ensure_events(events)
            # Allow reading empty epochs (ToDo: Maybe not anymore in the future)
            if len(events) == 0:
                self._allow_empty = True
                selection = None
            else:
                self._allow_empty = False
                events_max = events.max()
                if events_max > INT32_MAX:
                    raise ValueError(
                        f"events array values must not exceed {INT32_MAX}, "
                        f"got {events_max}"
                    )
        event_id = _check_event_id(event_id, events)
        self.event_id = event_id
        del event_id

        if events is not None:  # RtEpochs can have events=None
            for key, val in self.event_id.items():
                if val not in events[:, 2]:
                    msg = f"No matching events found for {key} (event id {val})"
                    _on_missing(on_missing, msg)

            # ensure metadata matches original events size
            self.selection = np.arange(len(events))
            self.events = events

            # same as self.metadata = metadata, but suppress log in favor
            # of logging below (after setting self.selection)
            GetEpochsMixin.metadata.fset(self, metadata, verbose=False)
            del events

            values = list(self.event_id.values())
            selected = np.where(np.isin(self.events[:, 2], values))[0]
            if selection is None:
                selection = selected
            else:
                selection = np.array(selection, int)
            if selection.shape != (len(selected),):
                raise ValueError(
                    f"selection must be shape {selected.shape} got shape "
                    f"{selection.shape}"
                )
            self.selection = selection
            if drop_log is None:
                self.drop_log = tuple(
                    () if k in self.selection else ("IGNORED",)
                    for k in range(max(len(self.events), max(self.selection) + 1))
                )
            else:
                self.drop_log = drop_log

            self.events = self.events[selected]

            (
                self.events,
                self.event_id,
                self.selection,
                self.drop_log,
            ) = _handle_event_repeated(
                self.events,
                self.event_id,
                event_repeated,
                self.selection,
                self.drop_log,
            )

            # then subselect
            sub = np.where(np.isin(selection, self.selection))[0]
            if isinstance(metadata, list):
                metadata = [metadata[s] for s in sub]
            elif metadata is not None:
                metadata = metadata.iloc[sub]

            # Remove temporarily set metadata from above, and set
            # again to get the correct log ("adding metadata", instead of
            # "replacing existing metadata")
            GetEpochsMixin.metadata.fset(self, None, verbose=False)
            self.metadata = metadata
            del metadata

            n_events = len(self.events)
            if n_events > 1:
                if np.diff(self.events.astype(np.int64)[:, 0]).min() <= 0:
                    warn(
                        "The events passed to the Epochs constructor are not "
                        "chronologically ordered.",
                        RuntimeWarning,
                    )

            if n_events > 0:
                logger.info(f"{n_events} matching events found")
            else:
                # Allow reading empty epochs (ToDo: Maybe not anymore in the future)
                if not self._allow_empty:
                    raise ValueError("No desired events found.")
        else:
            self.drop_log = tuple()
            self.selection = np.array([], int)
            self.metadata = metadata
            # do not set self.events here, let subclass do it

        if (detrend not in [None, 0, 1]) or isinstance(detrend, bool):
            raise ValueError("detrend must be None, 0, or 1")
        self.detrend = detrend

        self._raw = raw
        info._check_consistency()
        self.picks = _picks_to_idx(
            info, picks, none="all", exclude=(), allow_empty=False
        )
        self.info = pick_info(info, self.picks)
        del info
        self._current = 0

        if data is None:
            self.preload = False
            self._data = None
            self._do_baseline = True
        else:
            assert decim == 1
            if (
                data.ndim != 3
                or data.shape[2] != round((tmax - tmin) * self.info["sfreq"]) + 1
            ):
                raise RuntimeError("bad data shape")
            if data.shape[0] != len(self.events):
                raise ValueError(
                    "The number of epochs and the number of events must match"
                )
            self.preload = True
            self._data = data
            self._do_baseline = False
        self._offset = None

        if 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
        self._set_times(self._raw_times)

        # check reject_tmin and reject_tmax
        if reject_tmin is not None:
            if np.isclose(reject_tmin, tmin):
                # adjust for potential small deviations due to sampling freq
                reject_tmin = self.tmin
            elif reject_tmin < tmin:
                raise ValueError(
                    f"reject_tmin needs to be None or >= tmin (got {reject_tmin})"
                )

        if reject_tmax is not None:
            if np.isclose(reject_tmax, tmax):
                # adjust for potential small deviations due to sampling freq
                reject_tmax = self.tmax
            elif reject_tmax > tmax:
                raise ValueError(
                    f"reject_tmax needs to be None or <= tmax (got {reject_tmax})"
                )

        if (reject_tmin is not None) and (reject_tmax is not None):
            if reject_tmin >= reject_tmax:
                raise ValueError(
                    f"reject_tmin ({reject_tmin}) needs to be "
                    f" < reject_tmax ({reject_tmax})"
                )

        self.reject_tmin = reject_tmin
        self.reject_tmax = reject_tmax

        # decimation
        self._decim = 1
        self.decimate(decim)

        # baseline correction: replace `None` tuple elements  with actual times
        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 "
                f"[{self.baseline[0]}, {self.baseline[1]}] s"
            )

        logger.info(_log_rescale(self.baseline))

        # setup epoch rejection
        self.reject = None
        self.flat = None
        self._reject_setup(reject, flat)

        # do the rest
        valid_proj = [True, "delayed", False]
        if proj not in valid_proj:
            raise ValueError(f'"proj" must be one of {valid_proj}, not {proj}')
        if proj == "delayed":
            self._do_delayed_proj = True
            logger.info("Entering delayed SSP mode.")
        else:
            self._do_delayed_proj = False
        activate = False if self._do_delayed_proj else proj
        self._projector, self.info = setup_proj(self.info, False, activate=activate)
        if preload_at_end:
            assert self._data is None
            assert self.preload is False
            self.load_data()  # this will do the projection
        elif proj is True and self._projector is not None and data is not None:
            # let's make sure we project if data was provided and proj
            # requested
            # we could do this with np.einsum, but iteration should be
            # more memory safe in most instances
            for ii, epoch in enumerate(self._data):
                self._data[ii] = np.dot(self._projector, epoch)
        self.filename = filename if filename is not None else filename
        if raw_sfreq is None:
            raw_sfreq = self.info["sfreq"]
        self._raw_sfreq = raw_sfreq
        self._check_consistency()
        self.set_annotations(annotations, on_missing="ignore")

    def _check_consistency(self):
        """Check invariants of epochs object."""
        if hasattr(self, "events"):
            assert len(self.selection) == len(self.events)
            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"]
        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)

    def reset_drop_log_selection(self):
        """Reset the drop_log and selection entries.

        This method will simplify ``self.drop_log`` and ``self.selection``
        so that they are meaningless (tuple of empty tuples and increasing
        integers, respectively). This can be useful when concatenating
        many Epochs instances, as ``drop_log`` can accumulate many entries
        which can become problematic when saving.
        """
        self.selection = np.arange(len(self.events))
        self.drop_log = (tuple(),) * len(self.events)
        self._check_consistency()

    def load_data(self):
        """Load the data if not already preloaded.

        Returns
        -------
        epochs : instance of Epochs
            The epochs object.

        Notes
        -----
        This function operates in-place.

        .. versionadded:: 0.10.0
        """
        if self.preload:
            return self
        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)
        self._raw = None  # shouldn't need it anymore
        return self

    @verbose
    def apply_baseline(self, baseline=(None, 0), *, verbose=None):
        """Baseline correct epochs.

        Parameters
        ----------
        %(baseline_epochs)s
            Defaults to ``(None, 0)``, i.e. beginning of the the data until
            time point zero.
        %(verbose)s

        Returns
        -------
        epochs : instance of Epochs
            The baseline-corrected Epochs object.

        Notes
        -----
        Baseline correction can be done multiple times, but can never be
        reverted once the data has been loaded.

        .. versionadded:: 0.10.0
        """
        baseline = _check_baseline(baseline, times=self.times, sfreq=self.info["sfreq"])

        if self.preload:
            if self.baseline is not None and baseline is None:
                raise RuntimeError(
                    "You cannot remove baseline correction "
                    "from preloaded data once it has been "
                    "applied."
                )
            self._do_baseline = True
            picks = self._detrend_picks
            rescale(self._data, self.times, baseline, copy=False, picks=picks)
            self._do_baseline = False
        else:  # logging happens in "rescale" in "if" branch
            logger.info(_log_rescale(baseline))
            # For EpochsArray and Epochs, this is already True:
            # assert self._do_baseline is True
            # ... but for EpochsFIF it's not, so let's set it explicitly
            self._do_baseline = True
        self.baseline = baseline
        return self

    def _reject_setup(self, reject, flat, *, allow_callable=False):
        """Set self._reject_time and self._channel_type_idx."""
        idx = channel_indices_by_type(self.info)
        reject = deepcopy(reject) if reject is not None else dict()
        flat = deepcopy(flat) if flat is not None else dict()
        for rej, kind in zip((reject, flat), ("reject", "flat")):
            _validate_type(rej, dict, kind)
            bads = set(rej.keys()) - set(idx.keys())
            if len(bads) > 0:
                raise KeyError(f"Unknown channel types found in {kind}: {bads}")

        for key in idx.keys():
            # don't throw an error if rejection/flat would do nothing
            if len(idx[key]) == 0 and (
                np.isfinite(reject.get(key, np.inf)) or flat.get(key, -1) >= 0
            ):
                # This is where we could eventually add e.g.
                # self.allow_missing_reject_keys check to allow users to
                # provide keys that don't exist in data
                raise ValueError(
                    f"No {key.upper()} channel found. Cannot reject based on "
                    f"{key.upper()}."
                )

            # check for invalid values
            for rej, kind in zip((reject, flat), ("Rejection", "Flat")):
                for key, val in rej.items():
                    name = f"{kind} dict value for {key}"
                    if callable(val) and allow_callable:
                        continue
                    extra_str = ""
                    if allow_callable:
                        extra_str = "or callable"
                    _validate_type(val, "numeric", name, extra=extra_str)
                    if val is None or val < 0:
                        raise ValueError(
                            f"If using numerical {name} criteria, the value "
                            f"must be >= 0, not {repr(val)}"
                        )

        # now check to see if our rejection and flat are getting more
        # restrictive
        old_reject = self.reject if self.reject is not None else dict()
        old_flat = self.flat if self.flat is not None else dict()
        bad_msg = (
            '{kind}["{key}"] == {new} {op} {old} (old value), new '
            "{kind} values must be at least as stringent as "
            "previous ones"
        )

        # copy thresholds for channel types that were used previously, but not
        # passed this time
        for key in set(old_reject) - set(reject):
            reject[key] = old_reject[key]
        # make sure new thresholds are at least as stringent as the old ones
        for key in reject:
            # Skip this check if old_reject and reject are callables
            if callable(reject[key]) and allow_callable:
                continue
            if key in old_reject and reject[key] > old_reject[key]:
                raise ValueError(
                    bad_msg.format(
                        kind="reject",
                        key=key,
                        new=reject[key],
                        old=old_reject[key],
                        op=">",
                    )
                )

        # same for flat thresholds
        for key in set(old_flat) - set(flat):
            flat[key] = old_flat[key]
        for key in flat:
            if callable(flat[key]) and allow_callable:
                continue
            if key in old_flat and flat[key] < old_flat[key]:
                raise ValueError(
                    bad_msg.format(
                        kind="flat", key=key, new=flat[key], old=old_flat[key], op="<"
                    )
                )

        # after validation, set parameters
        self._bad_dropped = False
        self._channel_type_idx = idx
        self.reject = reject if len(reject) > 0 else None
        self.flat = flat if len(flat) > 0 else None

        if (self.reject_tmin is None) and (self.reject_tmax is None):
            self._reject_time = None
        else:
            if self.reject_tmin is None:
                reject_imin = None
            else:
                idxs = np.nonzero(self.times >= self.reject_tmin)[0]
                reject_imin = idxs[0]
            if self.reject_tmax is None:
                reject_imax = None
            else:
                idxs = np.nonzero(self.times <= self.reject_tmax)[0]
                reject_imax = idxs[-1]
            self._reject_time = slice(reject_imin, reject_imax)

    @verbose  # verbose is used by mne-realtime
    def _is_good_epoch(self, data, verbose=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 data.shape[1] < n_times:
            # epoch is too short ie at the end of the data
            return False, ("TOO_SHORT",)
        if self.reject is None and self.flat is None:
            return True, None
        else:
            if self._reject_time is not None:
                data = data[:, self._reject_time]

            return _is_good(
                data,
                self.ch_names,
                self._channel_type_idx,
                self.reject,
                self.flat,
                full_report=True,
                ignore_chs=self.info["bads"],
            )

    @verbose
    def _detrend_offset_decim(self, epoch, picks, verbose=None):
        """Aux Function: detrend, baseline correct, offset, decim.

        Note: operates inplace
        """
        if (epoch is None) or isinstance(epoch, str):
            return epoch

        # Detrend
        if self.detrend is not None:
            # We explicitly detrend just data channels (not EMG, ECG, EOG which
            # are processed by baseline correction)
            use_picks = _pick_data_channels(self.info, exclude=())
            epoch[use_picks] = detrend(epoch[use_picks], self.detrend, axis=1)

        # Baseline correct
        if self._do_baseline:
            rescale(
                epoch,
                self._raw_times,
                self.baseline,
                picks=picks,
                copy=False,
                verbose=False,
            )

        # Decimate if necessary (i.e., epoch not preloaded)
        epoch = epoch[:, self._decim_slice]

        # handle offset
        if self._offset is not None:
            epoch += self._offset

        return epoch

    def iter_evoked(self, copy=False):
        """Iterate over epochs as a sequence of Evoked objects.

        The Evoked objects yielded will each contain a single epoch (i.e., no
        averaging is performed).

        This method resets the object iteration state to the first epoch.

        Parameters
        ----------
        copy : bool
            If False copies of data and measurement info will be omitted
            to save time.
        """
        self.__iter__()

        while True:
            try:
                out = self.__next__(True)
            except StopIteration:
                break
            data, event_id = out
            tmin = self.times[0]
            info = self.info
            if copy:
                info = deepcopy(self.info)
                data = data.copy()

            yield EvokedArray(data, info, tmin, comment=str(event_id))

    def subtract_evoked(self, evoked=None):
        """Subtract an evoked response from each epoch.

        Can be used to exclude the evoked response when analyzing induced
        activity, see e.g. [1]_.

        Parameters
        ----------
        evoked : instance of Evoked | None
            The evoked response to subtract. If None, the evoked response
            is computed from Epochs itself.

        Returns
        -------
        self : instance of Epochs
            The modified instance (instance is also modified inplace).

        References
        ----------
        .. [1] David et al. "Mechanisms of evoked and induced responses in
               MEG/EEG", NeuroImage, vol. 31, no. 4, pp. 1580-1591, July 2006.
        """
        logger.info("Subtracting Evoked from Epochs")
        if evoked is None:
            picks = _pick_data_channels(self.info, exclude=[])
            evoked = self.average(picks)

        # find the indices of the channels to use
        picks = pick_channels(evoked.ch_names, include=self.ch_names, ordered=False)

        # make sure the omitted channels are not data channels
        if len(picks) < len(self.ch_names):
            sel_ch = [evoked.ch_names[ii] for ii in picks]
            diff_ch = list(set(self.ch_names).difference(sel_ch))
            diff_idx = [self.ch_names.index(ch) for ch in diff_ch]
            diff_types = [channel_type(self.info, idx) for idx in diff_idx]
            bad_idx = [
                diff_types.index(t) for t in diff_types if t in _DATA_CH_TYPES_SPLIT
            ]
            if len(bad_idx) > 0:
                bad_str = ", ".join([diff_ch[ii] for ii in bad_idx])
                raise ValueError(
                    "The following data channels are missing "
                    f"in the evoked response: {bad_str}"
                )
            logger.info(
                "    The following channels are not included in the subtraction: "
                + ", ".join(diff_ch)
            )

        # make sure the times match
        if (
            len(self.times) != len(evoked.times)
            or np.max(np.abs(self.times - evoked.times)) >= 1e-7
        ):
            raise ValueError(
                "Epochs and Evoked object do not contain the same time points."
            )

        # handle SSPs
        if not self.proj and evoked.proj:
            warn("Evoked has SSP applied while Epochs has not.")
        if self.proj and not evoked.proj:
            evoked = evoked.copy().apply_proj()

        # find the indices of the channels to use in Epochs
        ep_picks = [self.ch_names.index(evoked.ch_names[ii]) for ii in picks]

        # do the subtraction
        if self.preload:
            self._data[:, ep_picks, :] -= evoked.data[picks][None, :, :]
        else:
            if self._offset is None:
                self._offset = np.zeros(
                    (len(self.ch_names), len(self.times)), dtype=np.float64
                )
            self._offset[ep_picks] -= evoked.data[picks]
        logger.info("[done]")

        return self

    @fill_doc
    def average(self, picks=None, method="mean", by_event_type=False):
        """Compute an average over epochs.

        Parameters
        ----------
        %(picks_all_data)s
        method : str | callable
            How to combine the data. If "mean"/"median", the mean/median
            are returned.
            Otherwise, must be a callable which, when passed an array of shape
            (n_epochs, n_channels, n_time) returns an array of shape
            (n_channels, n_time).
            Note that due to file type limitations, the kind for all
            these will be "average".
        %(by_event_type)s

        Returns
        -------
        %(evoked_by_event_type_returns)s

        Notes
        -----
        Computes an average of all epochs in the instance, even if
        they correspond to different conditions. To average by condition,
        do ``epochs[condition].average()`` for each condition separately.

        When picks is None and epochs contain only ICA channels, no channels
        are selected, resulting in an error. This is because ICA channels
        are not considered data channels (they are of misc type) and only data
        channels are selected when picks is None.

        The ``method`` parameter allows e.g. robust averaging.
        For example, one could do:

            >>> from scipy.stats import trim_mean  # doctest:+SKIP
            >>> trim = lambda x: trim_mean(x, 0.1, axis=0)  # doctest:+SKIP
            >>> epochs.average(method=trim)  # doctest:+SKIP

        This would compute the trimmed mean.
        """
        self._handle_empty("raise", "average")
        if by_event_type:
            evokeds = list()
            for event_type in self.event_id.keys():
                ev = self[event_type]._compute_aggregate(picks=picks, mode=method)
                ev.comment = event_type
                evokeds.append(ev)
        else:
            evokeds = self._compute_aggregate(picks=picks, mode=method)
        return evokeds

    @fill_doc
    def standard_error(self, picks=None, by_event_type=False):
        """Compute standard error over epochs.

        Parameters
        ----------
        %(picks_all_data)s
        %(by_event_type)s

        Returns
        -------
        %(std_err_by_event_type_returns)s
        """
        return self.average(picks=picks, method="std", by_event_type=by_event_type)

    def _compute_aggregate(self, picks, mode="mean"):
        """Compute the mean, median, or std over epochs and return Evoked."""
        # if instance contains ICA channels they won't be included unless picks
        # is specified
        if picks is None:
            check_ICA = [x.startswith("ICA") for x in self.ch_names]
            if np.all(check_ICA):
                raise TypeError(
                    "picks must be specified (i.e. not None) for ICA channel data"
                )
            elif np.any(check_ICA):
                warn(
                    "ICA channels will not be included unless explicitly "
                    "selected in picks"
                )

        n_channels = len(self.ch_names)
        n_times = len(self.times)

        if self.preload:
            n_events = len(self.events)
            fun = _check_combine(mode, valid=("mean", "median", "std"))
            data = fun(self._data)
            assert len(self.events) == len(self._data)
            if data.shape != self._data.shape[1:]:
                raise RuntimeError(
                    f"You passed a function that resulted n data of shape "
                    f"{data.shape}, but it should be {self._data.shape[1:]}."
                )
        else:
            if mode not in {"mean", "std"}:
                raise ValueError(
                    "If data are not preloaded, can only compute "
                    "mean or standard deviation."
                )
            data = np.zeros((n_channels, n_times))
            n_events = 0
            for e in self:
                if np.iscomplexobj(e):
                    data = data.astype(np.complex128)
                data += e
                n_events += 1

            if n_events > 0:
                data /= n_events
            else:
                data.fill(np.nan)

            # convert to stderr if requested, could do in one pass but do in
            # two (slower) in case there are large numbers
            if mode == "std":
                data_mean = data.copy()
                data.fill(0.0)
                for e in self:
                    data += (e - data_mean) ** 2
                data = np.sqrt(data / n_events)

        if mode == "std":
            kind = "standard_error"
            data /= np.sqrt(n_events)
        else:
            kind = "average"

        return self._evoked_from_epoch_data(
            data, self.info, picks, n_events, kind, self._name
        )

    @property
    def _name(self):
        """Give a nice string representation based on event ids."""
        return self._get_name()

    def _get_name(self, count="frac", ms="×", sep="+"):
        """Generate human-readable name for epochs and evokeds from event_id.

        Parameters
        ----------
        count : 'frac' | 'total'
            Whether to include the fraction or total number of epochs that each
            event type contributes to the number of all epochs.
            Ignored if only one event type is present.
        ms : str | None
            The multiplication sign to use. Pass ``None`` to omit the sign.
            Ignored if only one event type is present.
        sep : str
            How to separate the different events names. Ignored if only one
            event type is present.
        """
        _check_option("count", value=count, allowed_values=["frac", "total"])

        if len(self.event_id) == 1:
            comment = next(iter(self.event_id.keys()))
        else:
            counter = Counter(self.events[:, 2])
            comments = list()

            # Take care of padding
            if ms is None:
                ms = " "
            else:
                ms = f" {ms} "

            for event_name, event_code in self.event_id.items():
                if count == "frac":
                    frac = float(counter[event_code]) / len(self.events)
                    comment = f"{frac:.2f}{ms}{event_name}"
                else:  # 'total'
                    comment = f"{counter[event_code]}{ms}{event_name}"
                comments.append(comment)

            comment = f" {sep} ".join(comments)
        return comment

    def _evoked_from_epoch_data(self, data, info, picks, n_events, kind, comment):
        """Create an evoked object from epoch data."""
        info = deepcopy(info)
        # don't apply baseline correction; we'll set evoked.baseline manually
        evoked = EvokedArray(
            data,
            info,
            tmin=self.times[0],
            comment=comment,
            nave=n_events,
            kind=kind,
            baseline=None,
        )
        evoked.baseline = self.baseline

        # the above constructor doesn't recreate the times object precisely
        # due to numerical precision issues
        evoked._set_times(self.times.copy())

        # pick channels
        picks = _picks_to_idx(self.info, picks, "data_or_ica", ())
        ch_names = [evoked.ch_names[p] for p in picks]
        evoked.pick(ch_names)

        if len(evoked.info["ch_names"]) == 0:
            raise ValueError("No data channel found when averaging.")

        if evoked.nave < 1:
            warn("evoked object is empty (based on less than 1 epoch)")

        return evoked

    @property
    def ch_names(self):
        """Channel names."""
        return self.info["ch_names"]

    @copy_function_doc_to_method_doc(plot_epochs)
    def plot(
        self,
        picks=None,
        scalings=None,
        n_epochs=20,
        n_channels=20,
        title=None,
        events=False,
        event_color=None,
        order=None,
        show=True,
        block=False,
        decim="auto",
        noise_cov=None,
        butterfly=False,
        show_scrollbars=True,
        show_scalebars=True,
        epoch_colors=None,
        event_id=None,
        group_by="type",
        precompute=None,
        use_opengl=None,
        *,
        theme=None,
        overview_mode=None,
        splash=True,
    ):
        return plot_epochs(
            self,
            picks=picks,
            scalings=scalings,
            n_epochs=n_epochs,
            n_channels=n_channels,
            title=title,
            events=events,
            event_color=event_color,
            order=order,
            show=show,
            block=block,
            decim=decim,
            noise_cov=noise_cov,
            butterfly=butterfly,
            show_scrollbars=show_scrollbars,
            show_scalebars=show_scalebars,
            epoch_colors=epoch_colors,
            event_id=event_id,
            group_by=group_by,
            precompute=precompute,
            use_opengl=use_opengl,
            theme=theme,
            overview_mode=overview_mode,
            splash=splash,
        )

    @copy_function_doc_to_method_doc(plot_topo_image_epochs)
    def plot_topo_image(
        self,
        layout=None,
        sigma=0.0,
        vmin=None,
        vmax=None,
        colorbar=None,
        order=None,
        cmap="RdBu_r",
        layout_scale=0.95,
        title=None,
        scalings=None,
        border="none",
        fig_facecolor="k",
        fig_background=None,
        font_color="w",
        select=False,
        show=True,
    ):
        return plot_topo_image_epochs(
            self,
            layout=layout,
            sigma=sigma,
            vmin=vmin,
            vmax=vmax,
            colorbar=colorbar,
            order=order,
            cmap=cmap,
            layout_scale=layout_scale,
            title=title,
            scalings=scalings,
            border=border,
            fig_facecolor=fig_facecolor,
            fig_background=fig_background,
            font_color=font_color,
            select=select,
            show=show,
        )

    @verbose
    def drop_bad(self, reject="existing", flat="existing", verbose=None):
        """Drop bad epochs without retaining the epochs data.

        Should be used before slicing operations.

        .. warning:: This operation is slow since all epochs have to be read
                     from disk. To avoid reading epochs from disk multiple
                     times, use :meth:`mne.Epochs.load_data()`.

        .. note:: To constrain the time period used for estimation of signal
                  quality, set ``epochs.reject_tmin`` and
                  ``epochs.reject_tmax``, respectively.

        Parameters
        ----------
        %(reject_drop_bad)s
        %(flat_drop_bad)s
        %(verbose)s

        Returns
        -------
        epochs : instance of Epochs
            The epochs with bad epochs dropped. Operates in-place.

        Notes
        -----
        Dropping bad epochs can be done multiple times with different
        ``reject`` and ``flat`` parameters. However, once an epoch is
        dropped, it is dropped forever, so if more lenient thresholds may
        subsequently be applied, :meth:`epochs.copy <mne.Epochs.copy>` should be
        used.
        """
        if reject == "existing":
            if flat == "existing" and self._bad_dropped:
                return
            reject = self.reject
        if flat == "existing":
            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"')
        self._reject_setup(reject, flat, allow_callable=True)
        self._get_data(out=False, verbose=verbose)
        return self

    def drop_log_stats(self, ignore=("IGNORED",)):
        """Compute the channel stats based on a drop_log from Epochs.

        Parameters
        ----------
        ignore : list
            The drop reasons to ignore.

        Returns
        -------
        perc : float
            Total percentage of epochs dropped.

        See Also
        --------
        plot_drop_log
        """
        return _drop_log_stats(self.drop_log, ignore)

    @copy_function_doc_to_method_doc(plot_drop_log)
    def plot_drop_log(
        self,
        threshold=0,
        n_max_plot=20,
        subject=None,
        color=(0.9, 0.9, 0.9),
        width=0.8,
        ignore=("IGNORED",),
        show=True,
    ):
        if not self._bad_dropped:
            raise ValueError(
                "You cannot use plot_drop_log since bad "
                "epochs have not yet been dropped. "
                "Use epochs.drop_bad()."
            )
        return plot_drop_log(
            self.drop_log,
            threshold,
            n_max_plot,
            subject,
            color=color,
            width=width,
            ignore=ignore,
            show=show,
        )

    @copy_function_doc_to_method_doc(plot_epochs_image)
    def plot_image(
        self,
        picks=None,
        sigma=0.0,
        vmin=None,
        vmax=None,
        colorbar=True,
        order=None,
        show=True,
        units=None,
        scalings=None,
        cmap=None,
        fig=None,
        axes=None,
        overlay_times=None,
        combine=None,
        group_by=None,
        evoked=True,
        ts_args=None,
        title=None,
        clear=False,
    ):
        return plot_epochs_image(
            self,
            picks=picks,
            sigma=sigma,
            vmin=vmin,
            vmax=vmax,
            colorbar=colorbar,
            order=order,
            show=show,
            units=units,
            scalings=scalings,
            cmap=cmap,
            fig=fig,
            axes=axes,
            overlay_times=overlay_times,
            combine=combine,
            group_by=group_by,
            evoked=evoked,
            ts_args=ts_args,
            title=title,
            clear=clear,
        )

    @verbose
    def drop(self, indices, reason="USER", verbose=None):
        """Drop epochs based on indices or boolean mask.

        .. note:: The indices refer to the current set of undropped epochs
                  rather than the complete set of dropped and undropped epochs.
                  They are therefore not necessarily consistent with any
                  external indices (e.g., behavioral logs). To drop epochs
                  based on external criteria, do not use the ``preload=True``
                  flag when constructing an Epochs object, and call this
                  method before calling the :meth:`mne.Epochs.drop_bad` or
                  :meth:`mne.Epochs.load_data` methods.

        Parameters
        ----------
        indices : array of int or bool
            Set epochs to remove by specifying indices to remove or a boolean
            mask to apply (where True values get removed). Events are
            correspondingly modified.
        reason : list | tuple | str
            Reason(s) for dropping the epochs ('ECG', 'timeout', 'blink' etc).
            Reason(s) are applied to all indices specified.
            Default: 'USER'.
        %(verbose)s

        Returns
        -------
        epochs : instance of Epochs
            The epochs with indices dropped. Operates in-place.
        """
        indices = np.atleast_1d(indices)

        if indices.ndim > 1:
            raise TypeError("indices must be a scalar or a 1-d array")
        # Check if indices and reasons are of the same length
        # if using collection to drop epochs

        if indices.dtype == np.dtype(bool):
            indices = np.where(indices)[0]
        try_idx = np.where(indices < 0, indices + len(self.events), indices)

        out_of_bounds = (try_idx < 0) | (try_idx >= len(self.events))
        if out_of_bounds.any():
            first = indices[out_of_bounds][0]
            raise IndexError(f"Epoch index {first} is out of bounds")
        keep = np.setdiff1d(np.arange(len(self.events)), try_idx)
        self._getitem(keep, reason, copy=False, drop_event_id=False)
        count = len(try_idx)
        logger.info(
            "Dropped %d epoch%s: %s",
            count,
            _pl(count),
            ", ".join(map(str, np.sort(try_idx))),
        )

        return self

    def _get_epoch_from_raw(self, idx, verbose=None):
        """Get a given epoch from disk."""
        raise NotImplementedError

    def _project_epoch(self, epoch):
        """Process a raw epoch based on the delayed param."""
        # whenever requested, the first epoch is being projected.
        if (epoch is None) or isinstance(epoch, str):
            # can happen if t < 0 or reject based on annotations
            return epoch
        proj = self._do_delayed_proj or self.proj
        if self._projector is not None and proj is True:
            epoch = np.dot(self._projector, epoch)
        return epoch

    def _handle_empty(self, on_empty, meth):
        if len(self.events) == 0:
            msg = (
                f"epochs.{meth}() can't run because this Epochs-object is empty. "
                f"You might want to check Epochs.drop_log or Epochs.plot_drop_log()"
                f" to see why epochs were dropped."
            )
            _on_missing(on_empty, msg, error_klass=RuntimeError)

    @verbose
    def _get_data(
        self,
        out=True,
        picks=None,
        item=None,
        *,
        units=None,
        tmin=None,
        tmax=None,
        copy=False,
        on_empty="warn",
        verbose=None,
    ):
        """Load all data, dropping bad epochs along the way.

        Parameters
        ----------
        out : bool
            Return the data. Setting this to False is used to reject bad
            epochs without caching all the data, which saves memory.
        %(picks_all)s
        item : slice | array-like | str | list | None
            See docstring of get_data method.
        %(units)s
        tmin : int | float | None
            Start time of data to get in seconds.
        tmax : int | float | None
            End time of data to get in seconds.
        %(verbose)s
        """
        from .io.base import _get_ch_factors

        if copy is not None:
            _validate_type(copy, bool, "copy")

        # Handle empty epochs
        self._handle_empty(on_empty, "_get_data")
        # if called with 'out=False', the call came from 'drop_bad()'
        # if no reasons to drop, just declare epochs as good and return
        if not out:
            # make sure first and last epoch not out of bounds of raw
            in_bounds = self.preload or (
                self._get_epoch_from_raw(idx=0) is not None
                and self._get_epoch_from_raw(idx=-1) is not None
            )
            # might be BaseEpochs or Epochs, only the latter has the attribute
            reject_by_annotation = getattr(self, "reject_by_annotation", False)
            if (
                self.reject is None
                and self.flat is None
                and in_bounds
                and self._reject_time is None
                and not reject_by_annotation
            ):
                logger.debug("_get_data is a noop, returning")
                self._bad_dropped = True
                return None
        start, stop = self._handle_tmin_tmax(tmin, tmax)

        if item is None:
            item = slice(None)
        elif not self._bad_dropped:
            raise ValueError(
                "item must be None in epochs.get_data() unless bads have been "
                "dropped. Consider using epochs.drop_bad()."
            )
        select = self._item_to_select(item)  # indices or slice
        use_idx = np.arange(len(self.events))[select]
        n_events = len(use_idx)
        # in case there are no good events
        if self.preload:
            # we will store our result in our existing array
            data = self._data
        else:
            # we start out with an empty array, allocate only if necessary
            data = np.empty((0, len(self.info["ch_names"]), len(self.times)))
            msg = (
                f"for {n_events} events and {len(self._raw_times)} original time points"
            )
            if self._decim > 1:
                msg += " (prior to decimation)"
            if getattr(self._raw, "preload", False):
                logger.info(f"Using data from preloaded Raw {msg} ...")
            else:
                logger.info(f"Loading data {msg} ...")

        orig_picks = picks
        if orig_picks is None:
            picks = _picks_to_idx(self.info, picks, "all", exclude=())
        else:
            picks = _picks_to_idx(self.info, picks)

        # handle units param only if we are going to return data (out==True)
        if (units is not None) and out:
            ch_factors = _get_ch_factors(self, units, picks)
        else:
            ch_factors = None

        if self._bad_dropped:
            if not out:
                return
            if self.preload:
                return self._data_sel_copy_scale(
                    data,
                    select=select,
                    orig_picks=orig_picks,
                    picks=picks,
                    ch_factors=ch_factors,
                    start=start,
                    stop=stop,
                    copy=copy,
                )

            # we need to load from disk, drop, and return data
            detrend_picks = self._detrend_picks
            for ii, idx in enumerate(use_idx):
                # faster to pre-allocate memory here
                epoch_noproj = self._get_epoch_from_raw(idx)
                epoch_noproj = self._detrend_offset_decim(epoch_noproj, detrend_picks)
                if self._do_delayed_proj:
                    epoch_out = epoch_noproj
                else:
                    epoch_out = self._project_epoch(epoch_noproj)
                if ii == 0:
                    data = np.empty(
                        (n_events, len(self.ch_names), len(self.times)),
                        dtype=epoch_out.dtype,
                    )
                data[ii] = epoch_out
        else:
            # bads need to be dropped, this might occur after a preload
            # e.g., when calling drop_bad w/new params
            good_idx = []
            n_out = 0
            drop_log = list(self.drop_log)
            assert n_events == len(self.selection)
            if not self.preload:
                detrend_picks = self._detrend_picks
            for idx, sel in enumerate(self.selection):
                if self.preload:  # from memory
                    if self._do_delayed_proj:
                        epoch_noproj = self._data[idx]
                        epoch = self._project_epoch(epoch_noproj)
                    else:
                        epoch_noproj = None
                        epoch = self._data[idx]
                else:  # from disk
                    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, verbose=verbose)
                if not is_good:
                    assert isinstance(bad_tuple, tuple)
                    assert all(isinstance(x, str) for x in bad_tuple)
                    drop_log[sel] = drop_log[sel] + bad_tuple
                    continue
                good_idx.append(idx)

                # store the epoch if there is a reason to (output or update)
                if out or self.preload:
                    # faster to pre-allocate, then trim as necessary
                    if n_out == 0 and not self.preload:
                        data = np.empty(
                            (n_events, epoch_out.shape[0], epoch_out.shape[1]),
                            dtype=epoch_out.dtype,
                            order="C",
                        )
                    data[n_out] = epoch_out
                    n_out += 1
            self.drop_log = tuple(drop_log)
            del drop_log

            self._bad_dropped = True
            n_bads_dropped = n_events - len(good_idx)
            logger.info(f"{n_bads_dropped} bad epochs dropped")

            if n_bads_dropped == n_events:
                warn(
                    "All epochs were dropped!\n"
                    "You might need to alter reject/flat-criteria "
                    "or drop bad channels to avoid this. "
                    "You can use Epochs.plot_drop_log() to see which "
                    "channels are responsible for the dropping of epochs."
                )

            # adjust the data size if there is a reason to (output or update)
            if out or self.preload:
                if data.flags["OWNDATA"] and data.flags["C_CONTIGUOUS"]:
                    data.resize((n_out,) + data.shape[1:], refcheck=False)
                else:
                    data = data[:n_out]
                    if self.preload:
                        self._data = data

            # Now update our properties (excepd data, which is already fixed)
            self._getitem(
                good_idx, None, copy=False, drop_event_id=False, select_data=False
            )

        if not out:
            return
        return self._data_sel_copy_scale(
            data,
            select=slice(None),
            orig_picks=orig_picks,
            picks=picks,
            ch_factors=ch_factors,
            start=start,
            stop=stop,
            copy=copy,
        )

    def _data_sel_copy_scale(
        self, data, *, select, orig_picks, picks, ch_factors, start, stop, copy
    ):
        # data arg starts out as self._data when data is preloaded
        data_is_self_data = bool(self.preload)
        logger.debug(f"Data is self data: {data_is_self_data}")
        # only two types of epoch subselection allowed
        assert isinstance(select, slice | np.ndarray), type(select)
        if not isinstance(select, slice):
            logger.debug("  Copying, fancy indexed epochs")
            data_is_self_data = False  # copy (fancy indexing)
        elif select != slice(None):
            logger.debug("  Slicing epochs")
        if orig_picks is not None:
            logger.debug("  Copying, fancy indexed picks")
            assert isinstance(picks, np.ndarray), type(picks)
            data_is_self_data = False  # copy (fancy indexing)
        else:
            picks = slice(None)
        if not all(isinstance(x, slice) and x == slice(None) for x in (select, picks)):
            data = data[select][:, picks]
        del picks
        if start != 0 or stop != self.times.size:
            logger.debug("  Slicing time")
            data = data[..., start:stop]  # view (slice)
        if ch_factors is not None:
            if data_is_self_data:
                logger.debug("  Copying, scale factors applied")
                data = data.copy()
                data_is_self_data = False
            data *= ch_factors[:, np.newaxis]
        if not data_is_self_data:
            return data
        if copy:
            logger.debug("  Copying, copy=True")
            data = data.copy()
        return data

    @property
    def _detrend_picks(self):
        if self._do_baseline:
            return _pick_data_channels(
                self.info, with_ref_meg=True, with_aux=True, exclude=()
            )
        else:
            return []

    @verbose
    def get_data(
        self,
        picks=None,
        item=None,
        units=None,
        tmin=None,
        tmax=None,
        *,
        copy=True,
        verbose=None,
    ):
        """Get all epochs as a 3D array.

        Parameters
        ----------
        %(picks_all)s
        item : slice | array-like | str | list | None
            The items to get. See :meth:`mne.Epochs.__getitem__` for
            a description of valid options. This can be substantially faster
            for obtaining an ndarray than :meth:`~mne.Epochs.__getitem__`
            for repeated access on large Epochs objects.
            None (default) is an alias for ``slice(None)``.

            .. versionadded:: 0.20
        %(units)s

            .. versionadded:: 0.24
        tmin : int | float | None
            Start time of data to get in seconds.

            .. versionadded:: 0.24.0
        tmax : int | float | None
            End time of data to get in seconds.

            .. versionadded:: 0.24.0
        copy : bool
            Whether to return a copy of the object's data, or (if possible) a view.
            See :ref:`the NumPy docs <numpy:basics.copies-and-views>` for an
            explanation. Default is ``False`` in 1.6 but will change to ``True`` in 1.7,
            set it explicitly to avoid a warning in some cases. A view is only possible
            when ``item is None``, ``picks is None``, ``units is None``, and data are
            preloaded.

            .. warning::
               Using ``copy=False`` and then modifying the returned ``data`` will in
               turn modify the Epochs object. Use with caution!

            .. versionchanged:: 1.7
               The default changed from ``False`` to ``True``.

            .. versionadded:: 1.6
        %(verbose)s

        Returns
        -------
        data : array of shape (n_epochs, n_channels, n_times)
            The epochs data. Will be a copy when ``copy=True`` and will be a view
            when possible when ``copy=False``.
        """
        return self._get_data(
            picks=picks, item=item, units=units, tmin=tmin, tmax=tmax, copy=copy
        )

    @verbose
    def apply_function(
        self,
        fun,
        picks=None,
        dtype=None,
        n_jobs=None,
        channel_wise=True,
        verbose=None,
        **kwargs,
    ):
        """Apply a function to a subset of channels.

        %(applyfun_summary_epochs)s

        Parameters
        ----------
        %(fun_applyfun)s
        %(picks_all_data_noref)s
        %(dtype_applyfun)s
        %(n_jobs)s Ignored if ``channel_wise=False`` as the workload
            is split across channels.
        %(channel_wise_applyfun_epo)s
        %(verbose)s
        %(kwargs_fun)s

        Returns
        -------
        self : instance of Epochs
            The epochs object with transformed data.
        """
        _check_preload(self, "epochs.apply_function")
        picks = _picks_to_idx(self.info, picks, exclude=(), with_ref_meg=False)

        if not callable(fun):
            raise ValueError("fun needs to be a function")

        data_in = self._data
        if dtype is not None and dtype != self._data.dtype:
            self._data = self._data.astype(dtype)

        args = getfullargspec(fun).args + getfullargspec(fun).kwonlyargs
        if channel_wise is False:
            if ("ch_idx" in args) or ("ch_name" in args):
                raise ValueError(
                    "apply_function cannot access ch_idx or ch_name "
                    "when channel_wise=False"
                )
        if "ch_idx" in args:
            logger.info("apply_function requested to access ch_idx")
        if "ch_name" in args:
            logger.info("apply_function requested to access ch_name")

        if channel_wise:
            parallel, p_fun, n_jobs = parallel_func(_check_fun, n_jobs)
            if n_jobs == 1:
                _fun = partial(_check_fun, fun)
                # modify data inplace to save memory
                for ch_idx in picks:
                    if "ch_idx" in args:
                        kwargs.update(ch_idx=ch_idx)
                    if "ch_name" in args:
                        kwargs.update(ch_name=self.info["ch_names"][ch_idx])
                    self._data[:, ch_idx, :] = np.apply_along_axis(
                        _fun, -1, data_in[:, ch_idx, :], **kwargs
                    )
            else:
                # use parallel function
                _fun = partial(np.apply_along_axis, fun, -1)
                data_picks_new = parallel(
                    p_fun(
                        _fun,
                        data_in[:, ch_idx, :],
                        **kwargs,
                        **{
                            k: v
                            for k, v in [
                                ("ch_name", self.info["ch_names"][ch_idx]),
                                ("ch_idx", ch_idx),
                            ]
                            if k in args
                        },
                    )
                    for ch_idx in picks
                )
                for run_idx, ch_idx in enumerate(picks):
                    self._data[:, ch_idx, :] = data_picks_new[run_idx]
        else:
            self._data = _check_fun(fun, data_in, **kwargs)

        return self

    @property
    def filename(self) -> Path | None:
        """The filename if the epochs are loaded from disk.

        :type: :class:`pathlib.Path` | ``None``
        """
        return self._filename

    @filename.setter
    def filename(self, value):
        if value is not None:
            value = _check_fname(value, overwrite="read", must_exist=True)
        self._filename = value

    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(".")
        s += " s (baseline "
        if self.baseline is None:
            s += "off"
        else:
            s += f"{self.baseline[0]:.3f}".rstrip("0").rstrip(".")
            s += f" – {self.baseline[1]:.3f}".rstrip("0").rstrip(".")
            s += " s"
            if self.baseline != _check_baseline(
                self.baseline,
                times=self.times,
                sfreq=self.info["sfreq"],
                on_baseline_outside_data="adjust",
            ):
                s += " (baseline period was cropped after baseline correction)"

        s += f"), ~{sizeof_fmt(self._size)}"
        s += f", data{'' if self.preload else ' not'} loaded"
        s += ", with metadata" if self.metadata is not None else ""
        max_events = 10
        counts = [
            f"{k!r}: {sum(self.events[:, 2] == v)}"
            for k, v in list(self.event_id.items())[:max_events]
        ]
        if len(self.event_id) > 0:
            s += "," + "\n ".join([""] + counts)
        if len(self.event_id) > max_events:
            not_shown_events = len(self.event_id) - max_events
            s += f"\n and {not_shown_events} more events ..."
        class_name = self.__class__.__name__
        class_name = "Epochs" if class_name == "BaseEpochs" else class_name
        return f"<{class_name} | {s}>"

    @repr_html
    def _repr_html_(self):
        if isinstance(self.event_id, dict):
            event_strings = []
            for k, v in sorted(self.event_id.items()):
                n_events = sum(self.events[:, 2] == v)
                event_strings.append(f"{k}: {n_events}")
        elif isinstance(self.event_id, list):
            event_strings = []
            for k in self.event_id:
                n_events = sum(self.events[:, 2] == k)
                event_strings.append(f"{k}: {n_events}")
        elif isinstance(self.event_id, int):
            n_events = len(self.events[:, 2])
            event_strings = [f"{self.event_id}: {n_events}"]
        else:
            event_strings = None

        t = _get_html_template("repr", "epochs.html.jinja")
        t = t.render(
            inst=self,
            filenames=(
                [Path(self.filename).name]
                if getattr(self, "filename", None) is not None
                else None
            ),
            event_counts=event_strings,
        )
        return t

    @verbose
    def crop(self, tmin=None, tmax=None, include_tmax=True, verbose=None):
        """Crop a time interval from the epochs.

        Parameters
        ----------
        tmin : float | None
            Start time of selection in seconds.
        tmax : float | None
            End time of selection in seconds.
        %(include_tmax)s
        %(verbose)s

        Returns
        -------
        epochs : instance of Epochs
            The cropped epochs object, modified in-place.

        Notes
        -----
        %(notes_tmax_included_by_default)s
        """
        # XXX this could be made to work on non-preloaded data...
        _check_preload(self, "Modifying data of epochs")

        super().crop(tmin=tmin, tmax=tmax, include_tmax=include_tmax)

        # Adjust rejection period
        if self.reject_tmin is not None and self.reject_tmin < self.tmin:
            logger.info(
                f"reject_tmin is not in epochs time interval. "
                f"Setting reject_tmin to epochs.tmin ({self.tmin} s)"
            )
            self.reject_tmin = self.tmin
        if self.reject_tmax is not None and self.reject_tmax > self.tmax:
            logger.info(
                f"reject_tmax is not in epochs time interval. "
                f"Setting reject_tmax to epochs.tmax ({self.tmax} s)"
            )
            self.reject_tmax = self.tmax
        return self

    def copy(self):
        """Return copy of Epochs instance.

        Returns
        -------
        epochs : instance of Epochs
            A copy of the object.
        """
        return deepcopy(self)

    def __deepcopy__(self, memodict):
        """Make a deepcopy."""
        cls = self.__class__
        result = cls.__new__(cls)
        for k, v in self.__dict__.items():
            # drop_log is immutable and _raw is private (and problematic to
            # deepcopy)
            if k in ("drop_log", "_raw", "_times_readonly"):
                memodict[id(v)] = v
            else:
                v = deepcopy(v, memodict)
            result.__dict__[k] = v
        return result

    @verbose
    def save(
        self,
        fname,
        split_size="2GB",
        fmt="single",
        overwrite=False,
        split_naming="neuromag",
        verbose=None,
    ):
        """Save epochs in a fif file.

        Parameters
        ----------
        fname : path-like
            The name of the file, which should end with ``-epo.fif`` or
            ``-epo.fif.gz``.
        split_size : str | int
            Large raw files are automatically split into multiple pieces. This
            parameter specifies the maximum size of each piece. If the
            parameter is an integer, it specifies the size in Bytes. It is
            also possible to pass a human-readable string, e.g., 100MB.
            Note: Due to FIFF file limitations, the maximum split size is 2GB.

            .. versionadded:: 0.10.0
        fmt : str
            Format to save data. Valid options are 'double' or
            'single' for 64- or 32-bit float, or for 128- or
            64-bit complex numbers respectively. Note: Data are processed with
            double precision. Choosing single-precision, the saved data
            will slightly differ due to the reduction in precision.

            .. versionadded:: 0.17
        %(overwrite)s
            To overwrite original file (the same one that was loaded),
            data must be preloaded upon reading. This defaults to True in 0.18
            but will change to False in 0.19.

            .. versionadded:: 0.18
        %(split_naming)s

            .. versionadded:: 0.24
        %(verbose)s

        Returns
        -------
        fnames : List of path-like
            List of path-like objects containing the path to each file split.
            .. versionadded:: 1.9

        Notes
        -----
        Bad epochs will be dropped before saving the epochs to disk.
        """
        check_fname(
            fname, "epochs", ("-epo.fif", "-epo.fif.gz", "_epo.fif", "_epo.fif.gz")
        )

        # check for file existence and expand `~` if present
        fname = str(
            _check_fname(
                fname=fname,
                overwrite=overwrite,
                check_bids_split=True,
                name="fname",
            )
        )

        split_size_bytes = _get_split_size(split_size)

        _check_option("fmt", fmt, ["single", "double"])

        # to know the length accurately. The get_data() call would drop
        # bad epochs anyway
        self.drop_bad()
        # total_size tracks sizes that get split
        # over_size tracks overhead (tags, things that get written to each)
        if len(self) == 0:
            warn("Saving epochs with no data")
            total_size = 0
        else:
            d = self[0].get_data(copy=False)
            # this should be guaranteed by subclasses
            assert d.dtype in (">f8", "<f8", ">c16", "<c16")
            total_size = d.nbytes * len(self)
        self._check_consistency()
        over_size = 0
        if fmt == "single":
            total_size //= 2  # 64bit data converted to 32bit before writing.
        over_size += 32  # FIF tags
        # Account for all the other things we write, too
        # 1. meas_id block plus main epochs block
        over_size += 132
        # 2. measurement info (likely slight overestimate, but okay)
        over_size += object_size(self.info) + 16 * len(self.info)
        # 3. events and event_id in its own block
        total_size += self.events.size * 4
        over_size += len(_event_id_string(self.event_id)) + 72
        # 4. Metadata in a block of its own
        if self.metadata is not None:
            total_size += len(_prepare_write_metadata(self.metadata))
        over_size += 56
        # 5. first sample, last sample, baseline
        over_size += 40 * (self.baseline is not None) + 40
        # 6. drop log: gets written to each, with IGNORE for ones that are
        #    not part of it. So make a fake one with all having entries.
        drop_size = len(json.dumps(self.drop_log)) + 16
        drop_size += 8 * (len(self.selection) - 1)  # worst case: all but one
        over_size += drop_size
        # 7. reject params
        reject_params = _pack_reject_params(self)
        if reject_params:
            over_size += len(json.dumps(reject_params)) + 16
        # 8. selection
        total_size += self.selection.size * 4
        over_size += 16
        # 9. end of file tags
        over_size += _NEXT_FILE_BUFFER
        logger.debug(f"    Overhead size:   {str(over_size).rjust(15)}")
        logger.debug(f"    Splittable size: {str(total_size).rjust(15)}")
        logger.debug(f"    Split size:      {str(split_size_bytes).rjust(15)}")
        # need at least one per
        n_epochs = len(self)
        n_per = total_size // n_epochs if n_epochs else 0
        min_size = n_per + over_size
        if split_size_bytes < min_size:
            raise ValueError(
                f"The split size {split_size} is too small to safely write "
                "the epochs contents, minimum split size is "
                f"{sizeof_fmt(min_size)} ({min_size} bytes)"
            )

        # This is like max(int(ceil(total_size / split_size)), 1) but cleaner
        n_parts = max((total_size - 1) // (split_size_bytes - over_size) + 1, 1)
        assert n_parts >= 1, n_parts
        if n_parts > 1:
            logger.info(f"Splitting into {n_parts} parts")
            if n_parts > 100:  # This must be an error
                raise ValueError(
                    f"Split size {split_size} would result in writing {n_parts} files"
                )

        if len(self.drop_log) > 100000:
            warn(
                f"epochs.drop_log contains {len(self.drop_log)} entries "
                f"which will incur up to a {sizeof_fmt(drop_size)} writing "
                f"overhead (per split file), consider using "
                f"epochs.reset_drop_log_selection() prior to writing"
            )

        epoch_idxs = np.array_split(np.arange(n_epochs), n_parts)

        _check_option("split_naming", split_naming, ("neuromag", "bids"))
        split_fnames = _make_split_fnames(fname, n_parts, split_naming)
        for part_idx, epoch_idx in enumerate(epoch_idxs):
            this_epochs = self[epoch_idx] if n_parts > 1 else self
            # avoid missing event_ids in splits
            this_epochs.event_id = self.event_id

            _save_split(this_epochs, split_fnames, part_idx, n_parts, fmt, overwrite)
        return split_fnames

    @verbose
    def export(self, fname, fmt="auto", *, overwrite=False, verbose=None):
        """Export Epochs to external formats.

        %(export_fmt_support_epochs)s

        %(export_warning)s

        Parameters
        ----------
        %(fname_export_params)s
        %(export_fmt_params_epochs)s
        %(overwrite)s

            .. versionadded:: 0.24.1
        %(verbose)s

        Notes
        -----
        .. versionadded:: 0.24

        %(export_warning_note_epochs)s
        %(export_eeglab_note)s
        """
        from .export import export_epochs

        export_epochs(fname, self, fmt, overwrite=overwrite, verbose=verbose)

    @fill_doc
    def equalize_event_counts(
        self, event_ids=None, method="mintime", *, random_state=None
    ):
        """Equalize the number of trials in each condition.

        It tries to make the remaining epochs occurring as close as possible in
        time. This method works based on the idea that if there happened to be
        some time-varying (like on the scale of minutes) noise characteristics
        during a recording, they could be compensated for (to some extent) in
        the equalization process. This method thus seeks to reduce any of
        those effects by minimizing the differences in the times of the events
        within a `~mne.Epochs` instance. For example, if one event type
        occurred at time points ``[1, 2, 3, 4, 120, 121]`` and the another one
        at ``[3.5, 4.5, 120.5, 121.5]``, this method would remove the events at
        times ``[1, 2]`` for the first event type – and not the events at times
        ``[120, 121]``.

        Parameters
        ----------
        event_ids : None | list | dict
            The event types to equalize.

            If ``None`` (default), equalize the counts of **all** event types
            present in the `~mne.Epochs` instance.

            If a list, each element can either be a string (event name) or a
            list of strings. In the case where one of the entries is a list of
            strings, event types in that list will be grouped together before
            equalizing trial counts across conditions.

            If a dictionary, the keys are considered as the event names whose
            counts to equalize, i.e., passing ``dict(A=1, B=2)`` will have the
            same effect as passing ``['A', 'B']``. This is useful if you intend
            to pass an ``event_id`` dictionary that was used when creating
            `~mne.Epochs`.

            In the case where partial matching is used (using ``/`` in
            the event names), the event types will be matched according to the
            provided tags, that is, processing works as if the ``event_ids``
            matched by the provided tags had been supplied instead.
            The ``event_ids`` must identify non-overlapping subsets of the
            epochs.
        %(equalize_events_method)s
        %(random_state)s Used only if ``method='random'``.

        Returns
        -------
        epochs : instance of Epochs
            The modified instance. It is modified in-place.
        indices : array of int
            Indices from the original events list that were dropped.

        Notes
        -----
        For example (if ``epochs.event_id`` was ``{'Left': 1, 'Right': 2,
        'Nonspatial':3}``:

            epochs.equalize_event_counts([['Left', 'Right'], 'Nonspatial'])

        would equalize the number of trials in the ``'Nonspatial'`` condition
        with the total number of trials in the ``'Left'`` and ``'Right'``
        conditions combined.

        If multiple indices are provided (e.g. ``'Left'`` and ``'Right'`` in
        the example above), it is not guaranteed that after equalization the
        conditions will contribute equally. E.g., it is possible to end up
        with 70 ``'Nonspatial'`` epochs, 69 ``'Left'`` and 1 ``'Right'``.

        .. versionchanged:: 0.23
            Default to equalizing all events in the passed instance if no
            event names were specified explicitly.
        """
        from collections.abc import Iterable

        _validate_type(
            event_ids,
            types=(Iterable, None),
            item_name="event_ids",
            type_name="list-like or None",
        )
        if isinstance(event_ids, str):
            raise TypeError(
                f"event_ids must be list-like or None, but "
                f"received a string: {event_ids}"
            )

        if event_ids is None:
            event_ids = list(self.event_id)
        elif not event_ids:
            raise ValueError("event_ids must have at least one element")

        if not self._bad_dropped:
            self.drop_bad()
        # figure out how to equalize
        eq_inds = list()

        # deal with hierarchical tags
        ids = self.event_id
        orig_ids = list(event_ids)
        tagging = False
        if "/" in "".join(ids):
            # make string inputs a list of length 1
            event_ids = [[x] if isinstance(x, str) else x for x in event_ids]
            for ids_ in event_ids:  # check if tagging is attempted
                if any([id_ not in ids for id_ in ids_]):
                    tagging = True
            # 1. treat everything that's not in event_id as a tag
            # 2a. for tags, find all the event_ids matched by the tags
            # 2b. for non-tag ids, just pass them directly
            # 3. do this for every input
            event_ids = [
                [
                    k for k in ids if all(tag in k.split("/") for tag in id_)
                ]  # ids matching all tags
                if all(id__ not in ids for id__ in id_)
                else id_  # straight pass for non-tag inputs
                for id_ in event_ids
            ]
            for ii, id_ in enumerate(event_ids):
                if len(id_) == 0:
                    raise KeyError(
                        f"{orig_ids[ii]} not found in the epoch object's event_id."
                    )
                elif len({sub_id in ids for sub_id in id_}) != 1:
                    err = (
                        "Don't mix hierarchical and regular event_ids"
                        f" like in '{', '.join(id_)}'."
                    )
                    raise ValueError(err)

            # raise for non-orthogonal tags
            if tagging is True:
                events_ = [set(self[x].events[:, 0]) for x in event_ids]
                doubles = events_[0].intersection(events_[1])
                if len(doubles):
                    raise ValueError(
                        "The two sets of epochs are "
                        "overlapping. Provide an "
                        "orthogonal selection."
                    )

        for eq in event_ids:
            eq_inds.append(self._keys_to_idx(eq))

        sample_nums = [self.events[e, 0] for e in eq_inds]
        indices = _get_drop_indices(sample_nums, method, random_state)
        # need to re-index indices
        indices = np.concatenate([e[idx] for e, idx in zip(eq_inds, indices)])
        self.drop(indices, reason="EQUALIZED_COUNT")
        # actually remove the indices
        return self, indices

    @verbose
    def compute_psd(
        self,
        method="multitaper",
        fmin=0,
        fmax=np.inf,
        tmin=None,
        tmax=None,
        picks=None,
        proj=False,
        remove_dc=True,
        exclude=(),
        *,
        n_jobs=1,
        verbose=None,
        **method_kw,
    ):
        """Perform spectral analysis on sensor data.

        Parameters
        ----------
        %(method_psd)s
            Default is ``'multitaper'``.
        %(fmin_fmax_psd)s
        %(tmin_tmax_psd)s
        %(picks_good_data_noref)s
        %(proj_psd)s
        %(remove_dc)s
        %(exclude_psd)s
        %(n_jobs)s
        %(verbose)s
        %(method_kw_psd)s

        Returns
        -------
        spectrum : instance of EpochsSpectrum
            The spectral representation of each epoch.

        Notes
        -----
        .. versionadded:: 1.2

        References
        ----------
        .. footbibliography::
        """
        method = _validate_method(method, type(self).__name__)
        self._set_legacy_nfft_default(tmin, tmax, method, method_kw)

        return EpochsSpectrum(
            self,
            method=method,
            fmin=fmin,
            fmax=fmax,
            tmin=tmin,
            tmax=tmax,
            picks=picks,
            exclude=exclude,
            proj=proj,
            remove_dc=remove_dc,
            n_jobs=n_jobs,
            verbose=verbose,
            **method_kw,
        )

    @verbose
    def compute_tfr(
        self,
        method,
        freqs,
        *,
        tmin=None,
        tmax=None,
        picks=None,
        proj=False,
        output="power",
        average=False,
        return_itc=False,
        decim=1,
        n_jobs=None,
        verbose=None,
        **method_kw,
    ):
        """Compute a time-frequency representation of epoched data.

        Parameters
        ----------
        %(method_tfr_epochs)s
        %(freqs_tfr_epochs)s
        %(tmin_tmax_psd)s
        %(picks_good_data_noref)s
        %(proj_psd)s
        %(output_compute_tfr)s
        average : bool
            Whether to return average power across epochs (instead of single-trial
            power). ``average=True`` is not compatible with ``output="complex"`` or
            ``output="phase"``. Ignored if ``method="stockwell"`` (Stockwell method
            *requires* averaging). Default is ``False``.
        return_itc : bool
            Whether to return inter-trial coherence (ITC) as well as power estimates.
            If ``True`` then must specify ``average=True`` (or ``method="stockwell",
            average="auto"``). Default is ``False``.
        %(decim_tfr)s
        %(n_jobs)s
        %(verbose)s
        %(method_kw_epochs_tfr)s

        Returns
        -------
        tfr : instance of EpochsTFR or AverageTFR
            The time-frequency-resolved power estimates.
        itc : instance of AverageTFR
            The inter-trial coherence (ITC). Only returned if ``return_itc=True``.

        Notes
        -----
        If ``average=True`` (or ``method="stockwell", average="auto"``) the result will
        be an :class:`~mne.time_frequency.AverageTFR` instead of an
        :class:`~mne.time_frequency.EpochsTFR`.

        .. versionadded:: 1.7

        References
        ----------
        .. footbibliography::
        """
        if method == "stockwell" and not average:  # stockwell method *must* average
            logger.info(
                'Requested `method="stockwell"` so ignoring parameter `average=False`.'
            )
            average = True
        if average:
            # augment `output` value for use by tfr_array_* functions
            _check_option("output", output, ("power",), extra=" when average=True")
            method_kw["output"] = "avg_power_itc" if return_itc else "avg_power"
        else:
            msg = (
                "compute_tfr() got incompatible parameters `average=False` and `{}` "
                "({} requires averaging over epochs)."
            )
            if return_itc:
                raise ValueError(msg.format("return_itc=True", "computing ITC"))
            if method == "stockwell":
                raise ValueError(msg.format('method="stockwell"', "Stockwell method"))
            # `average` and `return_itc` both False, so "phase" and "complex" are OK
            _check_option("output", output, ("power", "phase", "complex"))
            method_kw["output"] = output

        if method == "stockwell":
            method_kw["return_itc"] = return_itc
            method_kw.pop("output")
            if isinstance(freqs, str):
                _check_option("freqs", freqs, "auto")
            else:
                _validate_type(freqs, "array-like")
                _check_option(
                    "freqs", np.array(freqs).shape, ((2,),), extra=" (wrong shape)."
                )
        if average:
            out = AverageTFR(
                inst=self,
                method=method,
                freqs=freqs,
                tmin=tmin,
                tmax=tmax,
                picks=picks,
                proj=proj,
                decim=decim,
                n_jobs=n_jobs,
                verbose=verbose,
                **method_kw,
            )
            # tfr_array_stockwell always returns ITC (but sometimes it's None)
            if hasattr(out, "_itc"):
                if out._itc is not None:
                    state = out.__getstate__()
                    state["data"] = out._itc
                    state["data_type"] = "Inter-trial coherence"
                    itc = AverageTFR(inst=state)
                    del out._itc
                    return out, itc
                del out._itc
            return out
        # now handle average=False
        return EpochsTFR(
            inst=self,
            method=method,
            freqs=freqs,
            tmin=tmin,
            tmax=tmax,
            picks=picks,
            proj=proj,
            decim=decim,
            n_jobs=n_jobs,
            verbose=verbose,
            **method_kw,
        )

    @verbose
    def plot_psd(
        self,
        fmin=0,
        fmax=np.inf,
        tmin=None,
        tmax=None,
        picks=None,
        proj=False,
        *,
        method="auto",
        average=False,
        dB=True,
        estimate="power",
        xscale="linear",
        area_mode="std",
        area_alpha=0.33,
        color="black",
        line_alpha=None,
        spatial_colors=True,
        sphere=None,
        exclude="bads",
        ax=None,
        show=True,
        n_jobs=1,
        verbose=None,
        **method_kw,
    ):
        """%(plot_psd_doc)s.

        Parameters
        ----------
        %(fmin_fmax_psd)s
        %(tmin_tmax_psd)s
        %(picks_good_data_noref)s
        %(proj_psd)s
        %(method_plot_psd_auto)s
        %(average_plot_psd)s
        %(dB_plot_psd)s
        %(estimate_plot_psd)s
        %(xscale_plot_psd)s
        %(area_mode_plot_psd)s
        %(area_alpha_plot_psd)s
        %(color_plot_psd)s
        %(line_alpha_plot_psd)s
        %(spatial_colors_psd)s
        %(sphere_topomap_auto)s

            .. versionadded:: 0.22.0
        exclude : list of str | 'bads'
            Channels names to exclude from being shown. If 'bads', the bad
            channels are excluded. Pass an empty list to plot all channels
            (including channels marked "bad", if any).

            .. versionadded:: 0.24.0
        %(ax_plot_psd)s
        %(show)s
        %(n_jobs)s
        %(verbose)s
        %(method_kw_psd)s

        Returns
        -------
        fig : instance of Figure
            Figure with frequency spectra of the data channels.

        Notes
        -----
        %(notes_plot_psd_meth)s
        """
        return super().plot_psd(
            fmin=fmin,
            fmax=fmax,
            tmin=tmin,
            tmax=tmax,
            picks=picks,
            proj=proj,
            reject_by_annotation=False,
            method=method,
            average=average,
            dB=dB,
            estimate=estimate,
            xscale=xscale,
            area_mode=area_mode,
            area_alpha=area_alpha,
            color=color,
            line_alpha=line_alpha,
            spatial_colors=spatial_colors,
            sphere=sphere,
            exclude=exclude,
            ax=ax,
            show=show,
            n_jobs=n_jobs,
            verbose=verbose,
            **method_kw,
        )

    @verbose
    def to_data_frame(
        self,
        picks=None,
        index=None,
        scalings=None,
        copy=True,
        long_format=False,
        time_format=None,
        *,
        verbose=None,
    ):
        """Export data in tabular structure as a pandas DataFrame.

        Channels are converted to columns in the DataFrame. By default,
        additional columns "time", "epoch" (epoch number), and "condition"
        (epoch event description) are added, unless ``index`` is not ``None``
        (in which case the columns specified in ``index`` will be used to form
        the DataFrame's index instead).

        Parameters
        ----------
        %(picks_all)s
        %(index_df_epo)s
            Valid string values are 'time', 'epoch', and 'condition'.
            Defaults to ``None``.
        %(scalings_df)s
        %(copy_df)s
        %(long_format_df_epo)s
        %(time_format_df)s

            .. versionadded:: 0.20
        %(verbose)s

        Returns
        -------
        %(df_return)s
        """
        # check pandas once here, instead of in each private utils function
        pd = _check_pandas_installed()  # noqa
        # arg checking
        valid_index_args = ["time", "epoch", "condition"]
        valid_time_formats = ["ms", "timedelta"]
        index = _check_pandas_index_arguments(index, valid_index_args)
        time_format = _check_time_format(time_format, valid_time_formats)
        # get data
        picks = _picks_to_idx(self.info, picks, "all", exclude=())
        data = self._get_data(on_empty="raise")[:, picks, :]
        times = self.times
        n_epochs, n_picks, n_times = data.shape
        data = np.hstack(data).T  # (time*epochs) x signals
        if copy:
            data = data.copy()
        data = _scale_dataframe_data(self, data, picks, scalings)
        # prepare extra columns / multiindex
        mindex = list()
        times = np.tile(times, n_epochs)
        times = _convert_times(times, time_format, self.info["meas_date"])
        mindex.append(("time", times))
        rev_event_id = {v: k for k, v in self.event_id.items()}
        conditions = [rev_event_id[k] for k in self.events[:, 2]]
        mindex.append(("condition", np.repeat(conditions, n_times)))
        mindex.append(("epoch", np.repeat(self.selection, n_times)))
        assert all(len(mdx) == len(mindex[0]) for mdx in mindex)
        # build DataFrame
        df = _build_data_frame(
            self,
            data,
            picks,
            long_format,
            mindex,
            index,
            default_index=["condition", "epoch", "time"],
        )
        return df

    def as_type(self, ch_type="grad", mode="fast"):
        """Compute virtual epochs using interpolated fields.

        .. Warning:: Using virtual epochs to compute inverse can yield
            unexpected results. The virtual channels have ``'_v'`` appended
            at the end of the names to emphasize that the data contained in
            them are interpolated.

        Parameters
        ----------
        ch_type : str
            The destination channel type. It can be 'mag' or 'grad'.
        mode : str
            Either ``'accurate'`` or ``'fast'``, determines the quality of the
            Legendre polynomial expansion used. ``'fast'`` should be sufficient
            for most applications.

        Returns
        -------
        epochs : instance of mne.EpochsArray
            The transformed epochs object containing only virtual channels.

        Notes
        -----
        This method returns a copy and does not modify the data it
        operates on. It also returns an EpochsArray instance.

        .. versionadded:: 0.20.0
        """
        from .forward import _as_meg_type_inst

        self._handle_empty("raise", "as_type")
        return _as_meg_type_inst(self, ch_type=ch_type, mode=mode)


def _drop_log_stats(drop_log, ignore=("IGNORED",)):
    """Compute drop log stats.

    Parameters
    ----------
    drop_log : list of list
        Epoch drop log from Epochs.drop_log.
    ignore : list
        The drop reasons to ignore.

    Returns
    -------
    perc : float
        Total percentage of epochs dropped.
    """
    if (
        not isinstance(drop_log, tuple)
        or not all(isinstance(d, tuple) for d in drop_log)
        or not all(isinstance(s, str) for d in drop_log for s in d)
    ):
        raise TypeError("drop_log must be a tuple of tuple of str")
    perc = 100 * np.mean(
        [len(d) > 0 for d in drop_log if not any(r in ignore for r in d)]
    )
    return perc


def make_metadata(
    events,
    event_id,
    tmin,
    tmax,
    sfreq,
    row_events=None,
    keep_first=None,
    keep_last=None,
):
    """Automatically generate metadata for use with `mne.Epochs` from events.

    This function mimics the epoching process (it constructs time windows
    around time-locked "events of interest") and collates information about
    any other events that occurred within those time windows. The information
    is returned as a :class:`pandas.DataFrame`, suitable for use as
    `~mne.Epochs` metadata: one row per time-locked event, and columns
    indicating presence or absence and latency of each ancillary event type.

    The function will also return a new ``events`` array and ``event_id``
    dictionary that correspond to the generated metadata, which together can then be
    readily fed into `~mne.Epochs`.

    Parameters
    ----------
    events : array, shape (m, 3)
        The :term:`events array <events>`. By default, the returned metadata
        :class:`~pandas.DataFrame` will have as many rows as the events array.
        To create rows for only a subset of events, pass the ``row_events``
        parameter.
    event_id : dict
        A mapping from event names (keys) to event IDs (values). The event
        names will be incorporated as columns of the returned metadata
        :class:`~pandas.DataFrame`.
    tmin, tmax : float | str | list of str | None
        If float, start and end of the time interval for metadata generation in seconds,
        relative to the time-locked event of the respective time window (the "row
        events").

        .. note::
           If you are planning to attach the generated metadata to
           `~mne.Epochs` and intend to include only events that fall inside
           your epoch's time interval, pass the same ``tmin`` and ``tmax``
           values here as you use for your epochs.

        If ``None``, the time window used for metadata generation is bounded by the
        ``row_events``. This is can be particularly practical if trial duration varies
        greatly, but each trial starts with a known event (e.g., a visual cue or
        fixation).

        .. note::
           If ``tmin=None``, the first time window for metadata generation starts with
           the first row event. If ``tmax=None``, the last time window for metadata
           generation ends with the last event in ``events``.

        If a string or a list of strings, the events bounding the metadata around each
        "row event". For ``tmin``, the events are assumed to occur **before** the row
        event, and for ``tmax``, the events are assumed to occur **after** – unless
        ``tmin`` or ``tmax`` are equal to a row event, in which case the row event
        serves as the bound.

        .. versionchanged:: 1.6.0
           Added support for ``None``.

        .. versionadded:: 1.7.0
           Added support for strings.
    sfreq : float
        The sampling frequency of the data from which the events array was
        extracted.
    row_events : list of str | str | None
        Event types around which to create the time windows. For each of these
        time-locked events, we will create a **row** in the returned metadata
        :class:`pandas.DataFrame`. If provided, the string(s) must be keys of
        ``event_id``. If ``None`` (default), rows are created for **all** event types
        present in ``event_id``.
    keep_first : str | list of str | None
        Specify subsets of :term:`hierarchical event descriptors` (HEDs,
        inspired by :footcite:`BigdelyShamloEtAl2013`) matching events of which
        the **first occurrence** within each time window shall be stored in
        addition to the original events.

        .. note::
           There is currently no way to retain **all** occurrences of a
           repeated event. The ``keep_first`` parameter can be used to specify
           subsets of HEDs, effectively creating a new event type that is the
           union of all events types described by the matching HED pattern.
           Only the very first event of this set will be kept.

        For example, you might have two response events types,
        ``response/left`` and ``response/right``; and in trials with both
        responses occurring, you want to keep only the first response. In this
        case, you can pass ``keep_first='response'``. This will add two new
        columns to the metadata: ``response``, indicating at what **time** the
        event  occurred, relative to the time-locked event; and
        ``first_response``, stating which **type** (``'left'`` or ``'right'``)
        of event occurred.
        To match specific subsets of HEDs describing different sets of events,
        pass a list of these subsets, e.g.
        ``keep_first=['response', 'stimulus']``. If ``None`` (default), no
        event aggregation will take place and no new columns will be created.

        .. note::
           By default, this function will always retain  the first instance
           of any event in each time window. For example, if a time window
           contains two ``'response'`` events, the generated ``response``
           column will automatically refer to the first of the two events. In
           this specific case, it is therefore **not** necessary to make use of
           the ``keep_first`` parameter – unless you need to differentiate
           between two types of responses, like in the example above.

    keep_last : list of str | None
        Same as ``keep_first``, but for keeping only the **last**  occurrence
        of matching events. The column indicating the **type** of an event
        ``myevent`` will be named ``last_myevent``.

    Returns
    -------
    metadata : pandas.DataFrame
        Metadata for each row event, with the following columns:

        - ``event_name``, with strings indicating the name of the time-locked
          event ("row event") for that specific time window

        - one column per event type in ``event_id``, with the same name; floats
          indicating the latency of the event in seconds, relative to the
          time-locked event

        - if applicable, additional columns named after the ``keep_first`` and
          ``keep_last`` event types; floats indicating the latency  of the
          event in seconds, relative to the time-locked event

        - if applicable, additional columns ``first_{event_type}`` and
          ``last_{event_type}`` for ``keep_first`` and ``keep_last`` event
          types, respetively; the values will be strings indicating which event
          types were matched by the provided HED patterns

    events : array, shape (n, 3)
        The events corresponding to the generated metadata, i.e. one
        time-locked event per row.
    event_id : dict
        The event dictionary corresponding to the new events array. This will
        be identical to the input dictionary unless ``row_events`` is supplied,
        in which case it will only contain the events provided there.

    Notes
    -----
    The time window used for metadata generation need not correspond to the
    time window used to create the `~mne.Epochs`, to which the metadata will
    be attached; it may well be much shorter or longer, or not overlap at all,
    if desired. This can be useful, for example, to include events that
    occurred before or after an epoch, e.g. during the inter-trial interval.
    If either ``tmin``, ``tmax``, or both are ``None``, or a string referring e.g. to a
    response event, the time window will typically vary, too.

    .. versionadded:: 0.23

    References
    ----------
    .. footbibliography::
    """
    pd = _check_pandas_installed()

    _validate_type(events, types=("array-like",), item_name="events")
    _validate_type(event_id, types=(dict,), item_name="event_id")
    _validate_type(sfreq, types=("numeric",), item_name="sfreq")
    _validate_type(tmin, types=("numeric", str, "array-like", None), item_name="tmin")
    _validate_type(tmax, types=("numeric", str, "array-like", None), item_name="tmax")
    _validate_type(row_events, types=(None, str, "array-like"), item_name="row_events")
    _validate_type(keep_first, types=(None, str, "array-like"), item_name="keep_first")
    _validate_type(keep_last, types=(None, str, "array-like"), item_name="keep_last")

    if not event_id:
        raise ValueError("event_id dictionary must contain at least one entry")

    def _ensure_list(x):
        if x is None:
            return []
        elif isinstance(x, str):
            return [x]
        else:
            return list(x)

    row_events = _ensure_list(row_events)
    keep_first = _ensure_list(keep_first)
    keep_last = _ensure_list(keep_last)

    # Turn tmin, tmax into a list if they're strings or arrays of strings
    try:
        _validate_type(tmin, types=(str, "array-like"), item_name="tmin")
        tmin = _ensure_list(tmin)
    except TypeError:
        pass

    try:
        _validate_type(tmax, types=(str, "array-like"), item_name="tmax")
        tmax = _ensure_list(tmax)
    except TypeError:
        pass

    keep_first_and_last = set(keep_first) & set(keep_last)
    if keep_first_and_last:
        raise ValueError(
            f"The event names in keep_first and keep_last must "
            f"be mutually exclusive. Specified in both: "
            f"{', '.join(sorted(keep_first_and_last))}"
        )
    del keep_first_and_last

    for param_name, values in dict(keep_first=keep_first, keep_last=keep_last).items():
        for first_last_event_name in values:
            try:
                match_event_names(event_id, [first_last_event_name])
            except KeyError:
                raise ValueError(
                    f'Event "{first_last_event_name}", specified in '
                    f"{param_name}, cannot be found in event_id dictionary"
                )

    # If tmin, tmax are strings, ensure these event names are present in event_id
    def _diff_input_strings_vs_event_id(input_strings, input_name, event_id):
        event_name_diff = sorted(set(input_strings) - set(event_id.keys()))
        if event_name_diff:
            raise ValueError(
                f"Present in {input_name}, but missing from event_id: "
                f"{', '.join(event_name_diff)}"
            )

    _diff_input_strings_vs_event_id(
        input_strings=row_events, input_name="row_events", event_id=event_id
    )
    if isinstance(tmin, list):
        _diff_input_strings_vs_event_id(
            input_strings=tmin, input_name="tmin", event_id=event_id
        )
    if isinstance(tmax, list):
        _diff_input_strings_vs_event_id(
            input_strings=tmax, input_name="tmax", event_id=event_id
        )

    # First and last sample of each epoch, relative to the time-locked event
    # This follows the approach taken in mne.Epochs
    # For strings and None, we don't know the start and stop samples in advance as the
    # time window can vary.
    if isinstance(tmin, type(None) | list):
        start_sample = None
    else:
        start_sample = int(round(tmin * sfreq))

    if isinstance(tmax, type(None) | list):
        stop_sample = None
    else:
        stop_sample = int(round(tmax * sfreq)) + 1

    # Make indexing easier
    # We create the DataFrame before subsetting the events so we end up with
    # indices corresponding to the original event indices. Not used for now,
    # but might come in handy sometime later
    events_df = pd.DataFrame(events, columns=("sample", "prev_id", "id"))
    id_to_name_map = {v: k for k, v in event_id.items()}

    # Only keep events that are of interest
    events = events[np.isin(events[:, 2], list(event_id.values()))]
    events_df = events_df.loc[events_df["id"].isin(event_id.values()), :]

    # Prepare & condition the metadata DataFrame

    # Avoid column name duplications if the exact same event name appears in
    # event_id.keys() and keep_first / keep_last simultaneously
    keep_first_cols = [col for col in keep_first if col not in event_id]
    keep_last_cols = [col for col in keep_last if col not in event_id]
    first_cols = [f"first_{col}" for col in keep_first_cols]
    last_cols = [f"last_{col}" for col in keep_last_cols]

    columns = [
        "event_name",
        *event_id.keys(),
        *keep_first_cols,
        *keep_last_cols,
        *first_cols,
        *last_cols,
    ]

    data = np.empty((len(events_df), len(columns)), float)
    metadata = pd.DataFrame(data=data, columns=columns, index=events_df.index)

    # Event names
    metadata["event_name"] = ""

    # Event times
    start_idx = 1
    stop_idx = start_idx + len(event_id.keys()) + len(keep_first_cols + keep_last_cols)
    metadata.iloc[:, start_idx:stop_idx] = np.nan

    # keep_first and keep_last names
    start_idx = stop_idx
    metadata[columns[start_idx:]] = None

    # We're all set, let's iterate over all events and fill in in the
    # respective cells in the metadata. We will subset this to include only
    # `row_events` later
    for row_event in events_df.itertuples(name="RowEvent"):
        row_idx = row_event.Index
        metadata.loc[row_idx, "event_name"] = id_to_name_map[row_event.id]

        # Determine which events fall into the current time window
        if start_sample is None and isinstance(tmin, list):
            # Lower bound is the the current or the closest previpus event with a name
            # in "tmin"; if there is no such event (e.g., beginning of the recording is
            # being approached), the upper lower becomes the last event in the
            # recording.
            prev_matching_events = events_df.loc[
                (events_df["sample"] <= row_event.sample)
                & (events_df["id"].isin([event_id[name] for name in tmin])),
                :,
            ]
            if prev_matching_events.size == 0:
                # No earlier matching event. Use the current one as the beginning of the
                # time window. This may occur at the beginning of a recording.
                window_start_sample = row_event.sample
            else:
                # At least one earlier matching event. Use the closest one.
                window_start_sample = prev_matching_events.iloc[-1]["sample"]
        elif start_sample is None:
            # Lower bound is the current event.
            window_start_sample = row_event.sample
        else:
            # Lower bound is determined by tmin.
            window_start_sample = row_event.sample + start_sample

        if stop_sample is None and isinstance(tmax, list):
            # Upper bound is the the current or the closest following event with a name
            # in "tmax"; if there is no such event (e.g., end of the recording is being
            # approached), the upper bound becomes the last event in the recording.
            next_matching_events = events_df.loc[
                (events_df["sample"] >= row_event.sample)
                & (events_df["id"].isin([event_id[name] for name in tmax])),
                :,
            ]
            if next_matching_events.size == 0:
                # No matching event after the current one; use the end of the recording
                # as upper bound. This may occur at the end of a recording.
                window_stop_sample = events_df["sample"].iloc[-1]
            else:
                # At least one matching later event. Use the closest one..
                window_stop_sample = next_matching_events.iloc[0]["sample"]
        elif stop_sample is None:
            # Upper bound: next event of the same type, or the last event (of
            # any type) if no later event of the same type can be found.
            next_events = events_df.loc[
                (events_df["sample"] > row_event.sample),
                :,
            ]
            if next_events.size == 0:
                # We've reached the last event in the recording.
                window_stop_sample = row_event.sample
            elif next_events.loc[next_events["id"] == row_event.id, :].size > 0:
                # There's still an event of the same type appearing after the
                # current event. Stop one sample short, we don't want to include that
                # last event here, but in the next iteration.
                window_stop_sample = (
                    next_events.loc[next_events["id"] == row_event.id, :].iloc[0][
                        "sample"
                    ]
                    - 1
                )
            else:
                # There are still events after the current one, but not of the
                # same type.
                window_stop_sample = next_events.iloc[-1]["sample"]
        else:
            # Upper bound is determined by tmax.
            window_stop_sample = row_event.sample + stop_sample

        events_in_window = events_df.loc[
            (events_df["sample"] >= window_start_sample)
            & (events_df["sample"] <= window_stop_sample),
            :,
        ]

        assert not events_in_window.empty

        # Store the metadata
        for event in events_in_window.itertuples(name="Event"):
            event_sample = event.sample - row_event.sample
            event_time = event_sample / sfreq
            event_time = 0 if np.isclose(event_time, 0) else event_time
            event_name = id_to_name_map[event.id]

            if not np.isnan(metadata.loc[row_idx, event_name]):
                # Event already exists in current time window!
                assert metadata.loc[row_idx, event_name] <= event_time

                if event_name not in keep_last:
                    continue

            metadata.loc[row_idx, event_name] = event_time

            # Handle keep_first and keep_last event aggregation
            for event_group_name in keep_first + keep_last:
                if event_name not in match_event_names(event_id, [event_group_name]):
                    continue

                if event_group_name in keep_first:
                    first_last_col = f"first_{event_group_name}"
                else:
                    first_last_col = f"last_{event_group_name}"

                old_time = metadata.loc[row_idx, event_group_name]
                if not np.isnan(old_time):
                    if (event_group_name in keep_first and old_time <= event_time) or (
                        event_group_name in keep_last and old_time >= event_time
                    ):
                        continue

                if event_group_name not in event_id:
                    # This is an HED. Strip redundant information from the
                    # event name
                    name = (
                        event_name.replace(event_group_name, "")
                        .replace("//", "/")
                        .strip("/")
                    )
                    metadata.loc[row_idx, first_last_col] = name
                    del name

                metadata.loc[row_idx, event_group_name] = event_time

    # Only keep rows of interest
    if row_events:
        event_id_timelocked = {
            name: val for name, val in event_id.items() if name in row_events
        }
        events = events[np.isin(events[:, 2], list(event_id_timelocked.values()))]
        metadata = metadata.loc[metadata["event_name"].isin(event_id_timelocked)]
        assert len(events) == len(metadata)
        event_id = event_id_timelocked

    return metadata, events, event_id


def _events_from_annotations(raw, events, event_id, annotations, on_missing):
    """Generate events and event_ids from annotations."""
    events, event_id_tmp = events_from_annotations(raw)
    if events.size == 0:
        raise RuntimeError(
            "No usable annotations found in the raw object. "
            "Either `events` must be provided or the raw "
            "object must have annotations to construct epochs"
        )
    if any(raw.annotations.duration > 0):
        logger.info(
            "Ignoring annotation durations and creating fixed-duration epochs "
            "around annotation onsets."
        )
    if event_id is None:
        event_id = event_id_tmp
    # if event_id is the names of events, map to events integers
    if isinstance(event_id, str):
        event_id = [event_id]
    if isinstance(event_id, list | tuple | set):
        if not set(event_id).issubset(set(event_id_tmp)):
            msg = (
                "No matching annotations found for event_id(s) "
                f"{set(event_id) - set(event_id_tmp)}"
            )
            _on_missing(on_missing, msg)
        # remove extras if on_missing not error
        event_id = set(event_id) & set(event_id_tmp)
        event_id = {my_id: event_id_tmp[my_id] for my_id in event_id}
        # remove any non-selected annotations
        annotations.delete(~np.isin(raw.annotations.description, list(event_id)))
    return events, event_id, annotations


@fill_doc
class Epochs(BaseEpochs):
    """Epochs extracted from a Raw instance.

    Parameters
    ----------
    %(raw_epochs)s

        .. note::
            If ``raw`` contains annotations, ``Epochs`` can be constructed around
            ``raw.annotations.onset``, but note that the durations of the annotations
            are ignored in this case.
    %(events_epochs)s

        .. versionchanged:: 1.7
            Allow ``events=None`` to use ``raw.annotations.onset`` as the source of
            epoch times.
    %(event_id)s
    %(epochs_tmin_tmax)s
    %(baseline_epochs)s
        Defaults to ``(None, 0)``, i.e. beginning of the the data until
        time point zero.
    %(picks_all)s
    preload : bool
        %(epochs_preload)s
    %(reject_epochs)s
    %(flat)s
    %(proj_epochs)s
    %(decim)s
    %(epochs_reject_tmin_tmax)s
    %(detrend_epochs)s
    %(on_missing_epochs)s
    %(reject_by_annotation_epochs)s
    %(metadata_epochs)s

        .. versionadded:: 0.16
    %(event_repeated_epochs)s
    %(verbose)s

    Attributes
    ----------
    %(info_not_none)s
    %(event_id_attr)s
    ch_names : list of string
        List of channel names.
    %(selection_attr)s
    preload : bool
        Indicates whether epochs are in memory.
    drop_log : tuple of tuple
        A tuple of the same length as the event array used to initialize the
        Epochs object. If the i-th original event is still part of the
        selection, drop_log[i] will be an empty tuple; otherwise it will be
        a tuple of the reasons the event is not longer in the selection, e.g.:

        - 'IGNORED'
            If it isn't part of the current subset defined by the user
        - 'NO_DATA' or 'TOO_SHORT'
            If epoch didn't contain enough data names of channels that exceeded
            the amplitude threshold
        - 'EQUALIZED_COUNTS'
            See :meth:`~mne.Epochs.equalize_event_counts`
        - 'USER'
            For user-defined reasons (see :meth:`~mne.Epochs.drop`).

        When dropping based on flat or reject parameters the tuple of
        reasons contains a tuple of channels that satisfied the rejection
        criteria.
    filename : str
        The filename of the object.
    times :  ndarray
        Time vector in seconds. Goes from ``tmin`` to ``tmax``. Time interval
        between consecutive time samples is equal to the inverse of the
        sampling frequency.

    See Also
    --------
    mne.epochs.combine_event_ids
    mne.Epochs.equalize_event_counts

    Notes
    -----
    When accessing data, Epochs are detrended, baseline-corrected, and
    decimated, then projectors are (optionally) applied.

    For indexing and slicing using ``epochs[...]``, see
    :meth:`mne.Epochs.__getitem__`.

    All methods for iteration over objects (using :meth:`mne.Epochs.__iter__`,
    :meth:`mne.Epochs.iter_evoked` or :meth:`mne.Epochs.next`) use the same
    internal state.

    If ``event_repeated`` is set to ``'merge'``, the coinciding events
    (duplicates) will be merged into a single event_id and assigned a new
    id_number as::

        event_id['{event_id_1}/{event_id_2}/...'] = new_id_number

    For example with the event_id ``{'aud': 1, 'vis': 2}`` and the events
    ``[[0, 0, 1], [0, 0, 2]]``, the "merge" behavior will update both event_id
    and events to be: ``{'aud/vis': 3}`` and ``[[0, 0, 3]]`` respectively.

    There is limited support for :class:`~mne.Annotations` in the
    :class:`~mne.Epochs` class. Currently annotations that are present in the
    :class:`~mne.io.Raw` object will be preserved in the resulting
    :class:`~mne.Epochs` object, but:

    1. It is not yet possible to add annotations
       to the Epochs object programmatically (via code) or interactively
       (through the plot window)
    2. Concatenating :class:`~mne.Epochs` objects
       that contain annotations is not supported, and any annotations will
       be dropped when concatenating.
    3. Annotations will be lost on save.
    """

    @verbose
    def __init__(
        self,
        raw,
        events=None,
        event_id=None,
        tmin=-0.2,
        tmax=0.5,
        baseline=(None, 0),
        picks=None,
        preload=False,
        reject=None,
        flat=None,
        proj=True,
        decim=1,
        reject_tmin=None,
        reject_tmax=None,
        detrend=None,
        on_missing="raise",
        reject_by_annotation=True,
        metadata=None,
        event_repeated="error",
        verbose=None,
    ):
        from .io import BaseRaw

        if not isinstance(raw, BaseRaw):
            raise ValueError(
                "The first argument to `Epochs` must be an instance of mne.io.BaseRaw"
            )
        info = deepcopy(raw.info)
        annotations = raw.annotations.copy()

        # proj is on when applied in Raw
        proj = proj or raw.proj

        self.reject_by_annotation = reject_by_annotation

        # keep track of original sfreq (needed for annotations)
        raw_sfreq = raw.info["sfreq"]

        # get events from annotations if no events given
        if events is None:
            events, event_id, annotations = _events_from_annotations(
                raw, events, event_id, annotations, on_missing
            )

        # call BaseEpochs constructor
        super().__init__(
            info,
            None,
            events,
            event_id,
            tmin,
            tmax,
            metadata=metadata,
            baseline=baseline,
            raw=raw,
            picks=picks,
            reject=reject,
            flat=flat,
            decim=decim,
            reject_tmin=reject_tmin,
            reject_tmax=reject_tmax,
            detrend=detrend,
            proj=proj,
            on_missing=on_missing,
            preload_at_end=preload,
            event_repeated=event_repeated,
            verbose=verbose,
            raw_sfreq=raw_sfreq,
            annotations=annotations,
        )

    @verbose
    def _get_epoch_from_raw(self, idx, verbose=None):
        """Load one epoch from disk.

        Returns
        -------
        data : array | str | None
            If string, it's details on rejection reason.
            If array, it's the data in the desired range (good segment)
            If None, it means no data is available.
        """
        if self._raw is None:
            # This should never happen, as raw=None only if preload=True
            raise ValueError(
                "An error has occurred, no valid raw file found. "
                "Please report this to the mne-python "
                "developers."
            )
        sfreq = self._raw.info["sfreq"]
        event_samp = self.events[idx, 0]
        # Read a data segment from "start" to "stop" in samples
        first_samp = self._raw.first_samp
        start = int(round(event_samp + self._raw_times[0] * sfreq))
        start -= first_samp
        stop = start + len(self._raw_times)

        # reject_tmin, and reject_tmax need to be converted to samples to
        # check the reject_by_annotation boundaries: reject_start, reject_stop
        reject_tmin = self.reject_tmin
        if reject_tmin is None:
            reject_tmin = self._raw_times[0]
        reject_start = int(round(event_samp + reject_tmin * sfreq))
        reject_start -= first_samp

        reject_tmax = self.reject_tmax
        if reject_tmax is None:
            reject_tmax = self._raw_times[-1]
        diff = int(round((self._raw_times[-1] - reject_tmax) * sfreq))
        reject_stop = stop - diff

        logger.debug(f"    Getting epoch for {start}-{stop}")
        data = self._raw._check_bad_segment(
            start,
            stop,
            self.picks,
            reject_start,
            reject_stop,
            self.reject_by_annotation,
        )
        return data


@fill_doc
class EpochsArray(BaseEpochs):
    """Epochs object from numpy array.

    Parameters
    ----------
    data : array, shape (n_epochs, n_channels, n_times)
        The channels' time series for each epoch. See notes for proper units of
        measure.
    %(info_not_none)s Consider using :func:`mne.create_info` to populate this
        structure.
    %(events_epochs)s
    %(tmin_epochs)s
    %(event_id)s
    %(reject_epochs)s
    %(flat)s
    %(epochs_reject_tmin_tmax)s
    %(baseline_epochs)s
        Defaults to ``None``, i.e. no baseline correction.
    %(proj_epochs)s
    %(on_missing_epochs)s
    %(metadata_epochs)s

        .. versionadded:: 0.16
    %(selection)s
    %(drop_log)s

        .. versionadded:: 1.3
    %(raw_sfreq)s

        .. versionadded:: 1.3
    %(verbose)s

    See Also
    --------
    create_info
    EvokedArray
    io.RawArray

    Notes
    -----
    Proper units of measure:

    * V: eeg, eog, seeg, dbs, emg, ecg, bio, ecog
    * T: mag
    * T/m: grad
    * M: hbo, hbr
    * Am: dipole
    * AU: misc

    EpochsArray does not set `Annotations`. If you would like to create
    simulated data with Annotations that are then preserved in the Epochs
    object, you would use `mne.io.RawArray` first and then create an
    `mne.Epochs` object.
    """

    @verbose
    def __init__(
        self,
        data,
        info,
        events=None,
        tmin=0.0,
        event_id=None,
        reject=None,
        flat=None,
        reject_tmin=None,
        reject_tmax=None,
        baseline=None,
        proj=True,
        on_missing="raise",
        metadata=None,
        selection=None,
        *,
        drop_log=None,
        raw_sfreq=None,
        verbose=None,
    ):
        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

        super().__init__(
            info,
            data,
            events,
            event_id,
            tmin,
            tmax,
            baseline,
            reject=reject,
            flat=flat,
            reject_tmin=reject_tmin,
            reject_tmax=reject_tmax,
            decim=1,
            metadata=metadata,
            selection=selection,
            proj=proj,
            on_missing=on_missing,
            drop_log=drop_log,
            raw_sfreq=raw_sfreq,
            verbose=verbose,
        )
        if self.baseline is not None:
            self._do_baseline = True
        if (
            len(events)
            != np.isin(self.events[:, 2], list(self.event_id.values())).sum()
        ):
            raise ValueError("The events must only contain event numbers from event_id")
        detrend_picks = self._detrend_picks
        for e in self._data:
            # This is safe without assignment b/c there is no decim
            self._detrend_offset_decim(e, detrend_picks)
        self.drop_bad()


def combine_event_ids(epochs, old_event_ids, new_event_id, copy=True):
    """Collapse event_ids from an epochs instance into a new event_id.

    Parameters
    ----------
    epochs : instance of Epochs
        The epochs to operate on.
    old_event_ids : str, or list
        Conditions to collapse together.
    new_event_id : dict, or int
        A one-element dict (or a single integer) for the new
        condition. Note that for safety, this cannot be any
        existing id (in epochs.event_id.values()).
    copy : bool
        Whether to return a new instance or modify in place.

    Returns
    -------
    epochs : instance of Epochs
        The modified epochs.

    Notes
    -----
    This For example (if epochs.event_id was ``{'Left': 1, 'Right': 2}``::

        combine_event_ids(epochs, ['Left', 'Right'], {'Directional': 12})

    would create a 'Directional' entry in epochs.event_id replacing
    'Left' and 'Right' (combining their trials).
    """
    epochs = epochs.copy() if copy else epochs
    old_event_ids = np.asanyarray(old_event_ids)
    if isinstance(new_event_id, int):
        new_event_id = {str(new_event_id): new_event_id}
    else:
        if not isinstance(new_event_id, dict):
            raise ValueError("new_event_id must be a dict or int")
        if not len(list(new_event_id.keys())) == 1:
            raise ValueError("new_event_id dict must have one entry")
    new_event_num = list(new_event_id.values())[0]
    new_event_num = operator.index(new_event_num)
    if new_event_num in epochs.event_id.values():
        raise ValueError("new_event_id value must not already exist")
    # could use .pop() here, but if a latter one doesn't exist, we're
    # in trouble, so run them all here and pop() later
    old_event_nums = np.array([epochs.event_id[key] for key in old_event_ids])
    # find the ones to replace
    inds = np.any(
        epochs.events[:, 2][:, np.newaxis] == old_event_nums[np.newaxis, :], axis=1
    )
    # replace the event numbers in the events list
    epochs.events[inds, 2] = new_event_num
    # delete old entries
    for key in old_event_ids:
        epochs.event_id.pop(key)
    # add the new entry
    epochs.event_id.update(new_event_id)
    return epochs


@fill_doc
def equalize_epoch_counts(epochs_list, method="mintime", *, random_state=None):
    """Equalize the number of trials in multiple Epochs or EpochsTFR instances.

    Parameters
    ----------
    epochs_list : list of Epochs instances
        The Epochs instances to equalize trial counts for.
    %(equalize_events_method)s
    %(random_state)s Used only if ``method='random'``.

    Notes
    -----
    The method ``'mintime'`` tries to make the remaining epochs occurring as close as
    possible in time. This method is motivated by the possibility that if there happened
    to be some time-varying (like on the scale of minutes) noise characteristics during
    a recording, they could be compensated for (to some extent) in the
    equalization process. This method thus seeks to reduce any of those effects
    by minimizing the differences in the times of the events in the two sets of
    epochs. For example, if one had event times [1, 2, 3, 4, 120, 121] and the
    other one had [3.5, 4.5, 120.5, 121.5], it would remove events at times
    [1, 2] in the first epochs and not [120, 121].

    Examples
    --------
    >>> equalize_epoch_counts([epochs1, epochs2])  # doctest: +SKIP
    """
    if not all(isinstance(epoch, BaseEpochs | EpochsTFR) for epoch in epochs_list):
        raise ValueError("All inputs must be Epochs instances")
    # make sure bad epochs are dropped
    for epoch in epochs_list:
        if not epoch._bad_dropped:
            epoch.drop_bad()
    sample_nums = [epoch.events[:, 0] for epoch in epochs_list]
    indices = _get_drop_indices(sample_nums, method, random_state)
    for epoch, inds in zip(epochs_list, indices):
        epoch.drop(inds, reason="EQUALIZED_COUNT")


def _get_drop_indices(sample_nums, method, random_state):
    """Get indices to drop from multiple event timing lists."""
    small_idx = np.argmin([e.size for e in sample_nums])
    small_epoch_indices = sample_nums[small_idx]
    _check_option("method", method, ["mintime", "truncate", "random"])
    indices = list()
    for event in sample_nums:
        if method == "mintime":
            mask = _minimize_time_diff(small_epoch_indices, event)
        elif method == "truncate":
            mask = np.ones(event.size, dtype=bool)
            mask[small_epoch_indices.size :] = False
        elif method == "random":
            rng = check_random_state(random_state)
            mask = np.zeros(event.size, dtype=bool)
            idx = rng.choice(
                np.arange(event.size), size=small_epoch_indices.size, replace=False
            )
            mask[idx] = True
        indices.append(np.where(np.logical_not(mask))[0])
    return indices


def _minimize_time_diff(t_shorter, t_longer):
    """Find a boolean mask to minimize timing differences."""
    keep = np.ones((len(t_longer)), dtype=bool)
    # special case: length zero or one
    if len(t_shorter) < 2:  # interp1d won't work
        keep.fill(False)
        if len(t_shorter) == 1:
            idx = np.argmin(np.abs(t_longer - t_shorter))
            keep[idx] = True
        return keep
    scores = np.ones(len(t_longer))
    x1 = np.arange(len(t_shorter))
    # The first set of keep masks to test
    kwargs = dict(copy=False, bounds_error=False, assume_sorted=True)
    shorter_interp = interp1d(x1, t_shorter, fill_value=t_shorter[-1], **kwargs)
    for ii in range(len(t_longer) - len(t_shorter)):
        scores.fill(np.inf)
        # set up the keep masks to test, eliminating any rows that are already
        # gone
        keep_mask = ~np.eye(len(t_longer), dtype=bool)[keep]
        keep_mask[:, ~keep] = False
        # Check every possible removal to see if it minimizes
        x2 = np.arange(len(t_longer) - ii - 1)
        t_keeps = np.array([t_longer[km] for km in keep_mask])
        longer_interp = interp1d(
            x2, t_keeps, axis=1, fill_value=t_keeps[:, -1], **kwargs
        )
        d1 = longer_interp(x1) - t_shorter
        d2 = shorter_interp(x2) - t_keeps
        scores[keep] = np.abs(d1, d1).sum(axis=1) + np.abs(d2, d2).sum(axis=1)
        keep[np.argmin(scores)] = False
    return keep


@verbose
def _is_good(
    e,
    ch_names,
    channel_type_idx,
    reject,
    flat,
    full_report=False,
    ignore_chs=(),
    verbose=None,
):
    """Test if data segment e is good according to reject and flat.

    The reject and flat parameters can accept functions as values.

    If full_report=True, it will give True/False as well as a list of all
    offending channels.
    """
    bad_tuple = tuple()
    has_printed = False
    checkable = np.ones(len(ch_names), dtype=bool)
    checkable[np.array([c in ignore_chs for c in ch_names], dtype=bool)] = False

    for refl, f, t in zip([reject, flat], [np.greater, np.less], ["", "flat"]):
        if refl is not None:
            for key, refl in refl.items():
                criterion = refl
                idx = channel_type_idx[key]
                name = key.upper()
                if len(idx) > 0:
                    e_idx = e[idx]
                    checkable_idx = checkable[idx]
                    # Check if criterion is a function and apply it
                    if callable(criterion):
                        result = criterion(e_idx)
                        _validate_type(result, tuple, "reject/flat output")
                        if len(result) != 2:
                            raise TypeError(
                                "Function criterion must return a tuple of length 2"
                            )
                        cri_truth, reasons = result
                        _validate_type(cri_truth, (bool, np.bool_), cri_truth, "bool")
                        _validate_type(
                            reasons, (str, list, tuple), reasons, "str, list, or tuple"
                        )
                        idx_deltas = np.where(np.logical_and(cri_truth, checkable_idx))[
                            0
                        ]
                    else:
                        deltas = np.max(e_idx, axis=1) - np.min(e_idx, axis=1)
                        idx_deltas = np.where(
                            np.logical_and(f(deltas, criterion), checkable_idx)
                        )[0]

                    if len(idx_deltas) > 0:
                        # Check to verify that refl is a callable that returns
                        # (bool, reason). Reason must be a str/list/tuple.
                        # If using tuple
                        if callable(refl):
                            if isinstance(reasons, str):
                                reasons = (reasons,)
                            for idx, reason in enumerate(reasons):
                                _validate_type(reason, str, reason)
                            bad_tuple += tuple(reasons)
                        else:
                            bad_names = [ch_names[idx[i]] for i in idx_deltas]
                            if not has_printed:
                                logger.info(
                                    f"    Rejecting {t} epoch based on {name} : "
                                    f"{bad_names}"
                                )
                                has_printed = True
                            if not full_report:
                                return False
                            else:
                                bad_tuple += tuple(bad_names)

    if not full_report:
        return True
    else:
        if bad_tuple == ():
            return True, None
        else:
            return False, bad_tuple


def _read_one_epoch_file(f, tree, preload):
    """Read a single FIF file."""
    with f as fid:
        #   Read the measurement info
        info, meas = read_meas_info(fid, tree, clean_bads=True)

        # read in the Annotations if they exist
        annotations = _read_annotations_fif(fid, tree)
        try:
            events, mappings = _read_events_fif(fid, tree)
        except ValueError as e:
            # Allow reading empty epochs (ToDo: Maybe not anymore in the future)
            if str(e) == "Could not find any events":
                events = np.empty((0, 3), dtype=np.int32)
                mappings = dict()
            else:
                raise
        #   Metadata
        metadata = None
        metadata_tree = dir_tree_find(tree, FIFF.FIFFB_MNE_METADATA)
        if len(metadata_tree) > 0:
            for dd in metadata_tree[0]["directory"]:
                kind = dd.kind
                pos = dd.pos
                if kind == FIFF.FIFF_DESCRIPTION:
                    metadata = read_tag(fid, pos).data
                    metadata = _prepare_read_metadata(metadata)
                    break

        #   Locate the data of interest
        processed = dir_tree_find(meas, FIFF.FIFFB_PROCESSED_DATA)
        del meas
        if len(processed) == 0:
            raise ValueError("Could not find processed data")

        epochs_node = dir_tree_find(tree, FIFF.FIFFB_MNE_EPOCHS)
        if len(epochs_node) == 0:
            # before version 0.11 we errantly saved with this tag instead of
            # an MNE tag
            epochs_node = dir_tree_find(tree, FIFF.FIFFB_MNE_EPOCHS)
            if len(epochs_node) == 0:
                epochs_node = dir_tree_find(tree, 122)  # 122 used before v0.11
                if len(epochs_node) == 0:
                    raise ValueError("Could not find epochs data")

        my_epochs = epochs_node[0]

        # Now find the data in the block
        data = None
        data_tag = None
        bmin, bmax = None, None
        baseline = None
        selection = None
        drop_log = None
        raw_sfreq = None
        reject_params = {}
        for k in range(my_epochs["nent"]):
            kind = my_epochs["directory"][k].kind
            pos = my_epochs["directory"][k].pos
            if kind == FIFF.FIFF_FIRST_SAMPLE:
                tag = read_tag(fid, pos)
                first = int(tag.data.item())
            elif kind == FIFF.FIFF_LAST_SAMPLE:
                tag = read_tag(fid, pos)
                last = int(tag.data.item())
            elif kind == FIFF.FIFF_EPOCH:
                # delay reading until later
                fid.seek(pos, 0)
                data_tag = _read_tag_header(fid, pos)
                data_tag.type = data_tag.type ^ (1 << 30)
            elif kind in [FIFF.FIFF_MNE_BASELINE_MIN, 304]:
                # Constant 304 was used before v0.11
                tag = read_tag(fid, pos)
                bmin = float(tag.data.item())
            elif kind in [FIFF.FIFF_MNE_BASELINE_MAX, 305]:
                # Constant 305 was used before v0.11
                tag = read_tag(fid, pos)
                bmax = float(tag.data.item())
            elif kind == FIFF.FIFF_MNE_EPOCHS_SELECTION:
                tag = read_tag(fid, pos)
                selection = np.array(tag.data)
            elif kind == FIFF.FIFF_MNE_EPOCHS_DROP_LOG:
                tag = read_tag(fid, pos)
                drop_log = tag.data
                drop_log = json.loads(drop_log)
                drop_log = tuple(tuple(x) for x in drop_log)
            elif kind == FIFF.FIFF_MNE_EPOCHS_REJECT_FLAT:
                tag = read_tag(fid, pos)
                reject_params = json.loads(tag.data)
            elif kind == FIFF.FIFF_MNE_EPOCHS_RAW_SFREQ:
                tag = read_tag(fid, pos)
                raw_sfreq = tag.data

        if bmin is not None or bmax is not None:
            baseline = (bmin, bmax)

        n_samp = last - first + 1
        logger.info("    Found the data of interest:")
        logger.info(
            f"        t = {1000 * first / info['sfreq']:10.2f} ... "
            f"{1000 * last / info['sfreq']:10.2f} ms"
        )
        if info["comps"] is not None:
            logger.info(
                f"        {len(info['comps'])} CTF compensation matrices available"
            )

        # Inspect the data
        if data_tag is None:
            raise ValueError("Epochs data not found")
        epoch_shape = (len(info["ch_names"]), n_samp)
        size_expected = len(events) * np.prod(epoch_shape)
        # on read double-precision is always used
        if data_tag.type == FIFF.FIFFT_FLOAT:
            datatype = np.float64
            fmt = ">f4"
        elif data_tag.type == FIFF.FIFFT_DOUBLE:
            datatype = np.float64
            fmt = ">f8"
        elif data_tag.type == FIFF.FIFFT_COMPLEX_FLOAT:
            datatype = np.complex128
            fmt = ">c8"
        elif data_tag.type == FIFF.FIFFT_COMPLEX_DOUBLE:
            datatype = np.complex128
            fmt = ">c16"
        fmt_itemsize = np.dtype(fmt).itemsize
        assert fmt_itemsize in (4, 8, 16)
        size_actual = data_tag.size // fmt_itemsize - 16 // fmt_itemsize

        if not size_actual == size_expected:
            raise ValueError(
                f"Incorrect number of samples ({size_actual} instead of "
                f"{size_expected})."
            )

        # Calibration factors
        cals = np.array(
            [
                [info["chs"][k]["cal"] * info["chs"][k].get("scale", 1.0)]
                for k in range(info["nchan"])
            ],
            np.float64,
        )

        # Read the data
        if preload:
            data = read_tag(fid, data_tag.pos).data.astype(datatype)
            data *= cals

        # Put it all together
        tmin = first / info["sfreq"]
        tmax = last / info["sfreq"]
        event_id = (
            {str(e): e for e in np.unique(events[:, 2])}
            if mappings is None
            else mappings
        )
        # In case epochs didn't have a FIFF.FIFF_MNE_EPOCHS_SELECTION tag
        # (version < 0.8):
        if selection is None:
            selection = np.arange(len(events))
        if drop_log is None:
            drop_log = ((),) * len(events)

    return (
        info,
        data,
        data_tag,
        events,
        event_id,
        metadata,
        tmin,
        tmax,
        baseline,
        selection,
        drop_log,
        epoch_shape,
        cals,
        reject_params,
        fmt,
        annotations,
        raw_sfreq,
    )


@verbose
def read_epochs(fname, proj=True, preload=True, verbose=None) -> "EpochsFIF":
    """Read epochs from a fif file.

    Parameters
    ----------
    %(fname_epochs)s
    %(proj_epochs)s
    preload : bool
        If True, read all epochs from disk immediately. If ``False``, epochs
        will be read on demand.
    %(verbose)s

    Returns
    -------
    epochs : instance of Epochs
        The epochs.
    """
    return EpochsFIF(fname, proj, preload, verbose)


class _RawContainer:
    """Helper for a raw data container."""

    def __init__(self, fid, data_tag, event_samps, epoch_shape, cals, fmt):
        self.fid = fid
        self.data_tag = data_tag
        self.event_samps = event_samps
        self.epoch_shape = epoch_shape
        self.cals = cals
        self.proj = False
        self.fmt = fmt

    def __del__(self):  # noqa: D105
        self.fid.close()


@fill_doc
class EpochsFIF(BaseEpochs):
    """Epochs read from disk.

    Parameters
    ----------
    %(fname_epochs)s
    %(proj_epochs)s
    preload : bool
        If True, read all epochs from disk immediately. If False, epochs will
        be read on demand.
    %(verbose)s

    See Also
    --------
    mne.Epochs
    mne.epochs.combine_event_ids
    mne.Epochs.equalize_event_counts
    """

    @verbose
    def __init__(self, fname, proj=True, preload=True, verbose=None):
        from .io.base import _get_fname_rep

        if _path_like(fname):
            check_fname(
                fname=fname,
                filetype="epochs",
                endings=("-epo.fif", "-epo.fif.gz", "_epo.fif", "_epo.fif.gz"),
            )
            fname = _check_fname(fname=fname, must_exist=True, overwrite="read")
        elif not preload:
            raise ValueError("preload must be used with file-like objects")

        fnames = [fname]
        fname_rep = _get_fname_rep(fname)
        ep_list = list()
        raw = list()
        for fname in fnames:
            logger.info(f"Reading {fname_rep} ...")
            fid, tree, _ = fiff_open(fname, preload=preload)
            next_fname = _get_next_fname(fid, fname, tree)
            (
                info,
                data,
                data_tag,
                events,
                event_id,
                metadata,
                tmin,
                tmax,
                baseline,
                selection,
                drop_log,
                epoch_shape,
                cals,
                reject_params,
                fmt,
                annotations,
                raw_sfreq,
            ) = _read_one_epoch_file(fid, tree, preload)

            if (events[:, 0] < 0).any():
                events = events.copy()
                warn(
                    "Incorrect events detected on disk, setting event "
                    "numbers to consecutive increasing integers"
                )
                events[:, 0] = np.arange(1, len(events) + 1)
            # here we ignore missing events, since users should already be
            # aware of missing events if they have saved data that way
            # we also retain original baseline without re-applying baseline
            # correction (data is being baseline-corrected when written to
            # disk)
            epoch = BaseEpochs(
                info,
                data,
                events,
                event_id,
                tmin,
                tmax,
                baseline=None,
                metadata=metadata,
                on_missing="ignore",
                selection=selection,
                drop_log=drop_log,
                proj=False,
                verbose=False,
                raw_sfreq=raw_sfreq,
            )
            epoch.baseline = baseline
            epoch._do_baseline = False  # might be superfluous but won't hurt
            ep_list.append(epoch)

            if not preload:
                # store everything we need to index back to the original data
                raw.append(
                    _RawContainer(
                        fiff_open(fname)[0],
                        data_tag,
                        events[:, 0].copy(),
                        epoch_shape,
                        cals,
                        fmt,
                    )
                )

            if next_fname is not None:
                fnames.append(next_fname)

        unsafe_annot_add = raw_sfreq is None
        (
            info,
            data,
            raw_sfreq,
            events,
            event_id,
            tmin,
            tmax,
            metadata,
            baseline,
            selection,
            drop_log,
        ) = _concatenate_epochs(
            ep_list,
            with_data=preload,
            add_offset=False,
            on_mismatch="raise",
        )
        # we need this uniqueness for non-preloaded data to work properly
        if len(np.unique(events[:, 0])) != len(events):
            raise RuntimeError("Event time samples were not unique")

        # correct the drop log
        assert len(drop_log) % len(fnames) == 0
        step = len(drop_log) // len(fnames)
        offsets = np.arange(step, len(drop_log) + 1, step)
        drop_log = list(drop_log)
        for i1, i2 in zip(offsets[:-1], offsets[1:]):
            other_log = drop_log[i1:i2]
            for k, (a, b) in enumerate(zip(drop_log, other_log)):
                if a == ("IGNORED",) and b != ("IGNORED",):
                    drop_log[k] = b
        drop_log = tuple(drop_log[:step])

        # call BaseEpochs constructor
        # again, ensure we're retaining the baseline period originally loaded
        # from disk without trying to re-apply baseline correction
        super().__init__(
            info,
            data,
            events,
            event_id,
            tmin,
            tmax,
            baseline=None,
            raw=raw,
            proj=proj,
            preload_at_end=False,
            on_missing="ignore",
            selection=selection,
            drop_log=drop_log,
            filename=fname_rep,
            metadata=metadata,
            verbose=verbose,
            raw_sfreq=raw_sfreq,
            annotations=annotations,
            **reject_params,
        )
        self.baseline = baseline
        self._do_baseline = False
        # use the private property instead of drop_bad so that epochs
        # are not all read from disk for preload=False
        self._bad_dropped = True
        # private property to suggest that people re-save epochs if they add
        # annotations
        self._unsafe_annot_add = unsafe_annot_add

    @verbose
    def _get_epoch_from_raw(self, idx, verbose=None):
        """Load one epoch from disk."""
        # Find the right file and offset to use
        event_samp = self.events[idx, 0]
        for raw in self._raw:
            idx = np.where(raw.event_samps == event_samp)[0]
            if len(idx) == 1:
                fmt = raw.fmt
                idx = idx[0]
                size = np.prod(raw.epoch_shape) * np.dtype(fmt).itemsize
                offset = idx * size + 16  # 16 = Tag header
                break
        else:
            # read the correct subset of the data
            raise RuntimeError(
                "Correct epoch could not be found, please contact mne-python developers"
            )
        # the following is equivalent to this, but faster:
        #
        # >>> data = read_tag(raw.fid, raw.data_tag.pos).data.astype(float)
        # >>> data *= raw.cals[np.newaxis, :, :]
        # >>> data = data[idx]
        #
        # Eventually this could be refactored in io/tag.py if other functions
        # could make use of it
        raw.fid.seek(raw.data_tag.pos + offset, 0)
        if fmt == ">c8":
            read_fmt = ">f4"
        elif fmt == ">c16":
            read_fmt = ">f8"
        else:
            read_fmt = fmt
        data = np.frombuffer(raw.fid.read(size), read_fmt)
        if read_fmt != fmt:
            data = data.view(fmt)
            data = data.astype(np.complex128)
        else:
            data = data.astype(np.float64)

        data.shape = raw.epoch_shape
        data *= raw.cals
        return data


@fill_doc
def bootstrap(epochs, random_state=None):
    """Compute epochs selected by bootstrapping.

    Parameters
    ----------
    epochs : Epochs instance
        epochs data to be bootstrapped
    %(random_state)s

    Returns
    -------
    epochs : Epochs instance
        The bootstrap samples
    """
    if not epochs.preload:
        raise RuntimeError(
            "Modifying data of epochs is only supported "
            "when preloading is used. Use preload=True "
            "in the constructor."
        )

    rng = check_random_state(random_state)
    epochs_bootstrap = epochs.copy()
    n_events = len(epochs_bootstrap.events)
    idx = rng_uniform(rng)(0, n_events, n_events)
    epochs_bootstrap = epochs_bootstrap[idx]
    return epochs_bootstrap


def _concatenate_epochs(
    epochs_list, *, with_data=True, add_offset=True, on_mismatch="raise"
):
    """Auxiliary function for concatenating epochs."""
    if not isinstance(epochs_list, list | tuple):
        raise TypeError(f"epochs_list must be a list or tuple, got {type(epochs_list)}")

    # to make warning messages only occur once during concatenation
    warned = False

    for ei, epochs in enumerate(epochs_list):
        if not isinstance(epochs, BaseEpochs):
            raise TypeError(
                f"epochs_list[{ei}] must be an instance of Epochs, got {type(epochs)}"
            )

        if (
            getattr(epochs, "annotations", None) is not None
            and len(epochs.annotations) > 0
            and not warned
        ):
            warned = True
            warn(
                "Concatenation of Annotations within Epochs is not supported yet. All "
                "annotations will be dropped."
            )

            # create a copy, so that the Annotations are not modified in place
            # from the original object
            epochs = epochs.copy()
            epochs.set_annotations(None)
    out = epochs_list[0]
    offsets = [0]
    if with_data:
        out.drop_bad()
        offsets.append(len(out))
    events = [out.events]
    metadata = [out.metadata]
    baseline, tmin, tmax = out.baseline, out.tmin, out.tmax
    raw_sfreq = out._raw_sfreq
    info = deepcopy(out.info)
    drop_log = out.drop_log
    event_id = deepcopy(out.event_id)
    selection = out.selection
    # offset is the last epoch + tmax + 10 second
    shift = np.int64((10 + tmax) * out.info["sfreq"])
    # Allow reading empty epochs (ToDo: Maybe not anymore in the future)
    if out._allow_empty:
        events_offset = 0
    else:
        events_offset = int(np.max(events[0][:, 0])) + shift
    events_offset = np.int64(events_offset)
    events_overflow = False
    warned = False
    for ii, epochs in enumerate(epochs_list[1:], 1):
        _ensure_infos_match(epochs.info, info, f"epochs[{ii}]", on_mismatch=on_mismatch)
        if not np.allclose(epochs.times, epochs_list[0].times):
            raise ValueError("Epochs must have same times")

        if epochs.baseline != baseline:
            raise ValueError("Baseline must be same for all epochs")

        if epochs._raw_sfreq != raw_sfreq and not warned:
            warned = True
            warn(
                "The original raw sampling rate of the Epochs does not "
                "match for all Epochs. Please proceed cautiously."
            )

        # compare event_id
        common_keys = list(set(event_id).intersection(set(epochs.event_id)))
        for key in common_keys:
            if not event_id[key] == epochs.event_id[key]:
                msg = (
                    "event_id values must be the same for identical keys "
                    'for all concatenated epochs. Key "{}" maps to {} in '
                    "some epochs and to {} in others."
                )
                raise ValueError(msg.format(key, event_id[key], epochs.event_id[key]))

        if with_data:
            epochs.drop_bad()
            offsets.append(len(epochs))
        evs = epochs.events.copy()
        if len(epochs.events) == 0:
            warn("One of the Epochs objects to concatenate was empty.")
        elif add_offset:
            # We need to cast to a native Python int here to detect an
            # overflow of a numpy int32 (which is the default on windows)
            max_timestamp = int(np.max(evs[:, 0]))
            evs[:, 0] += events_offset
            events_offset += max_timestamp + shift
            if events_offset > INT32_MAX:
                warn(
                    f"Event number greater than {INT32_MAX} created, "
                    "events[:, 0] will be assigned consecutive increasing "
                    "integer values"
                )
                events_overflow = True
                add_offset = False  # we no longer need to add offset
        events.append(evs)
        selection = np.concatenate((selection, epochs.selection))
        drop_log = drop_log + epochs.drop_log
        event_id.update(epochs.event_id)
        metadata.append(epochs.metadata)
    events = np.concatenate(events, axis=0)
    # check to see if we exceeded our maximum event offset
    if events_overflow:
        events[:, 0] = np.arange(1, len(events) + 1)

    # Create metadata object (or make it None)
    n_have = sum(this_meta is not None for this_meta in metadata)
    if n_have == 0:
        metadata = None
    elif n_have != len(metadata):
        raise ValueError(
            f"{n_have} of {len(metadata)} epochs instances have metadata, either "
            "all or none must have metadata"
        )
    else:
        pd = _check_pandas_installed(strict=False)
        if pd is not False:
            metadata = pd.concat(metadata)
        else:  # dict of dicts
            metadata = sum(metadata, list())
    assert len(offsets) == (len(epochs_list) if with_data else 0) + 1
    data = None
    if with_data:
        offsets = np.cumsum(offsets)
        for start, stop, epochs in zip(offsets[:-1], offsets[1:], epochs_list):
            this_data = epochs.get_data(copy=False)
            if data is None:
                data = np.empty(
                    (offsets[-1], len(out.ch_names), len(out.times)),
                    dtype=this_data.dtype,
                )
            data[start:stop] = this_data
    return (
        info,
        data,
        raw_sfreq,
        events,
        event_id,
        tmin,
        tmax,
        metadata,
        baseline,
        selection,
        drop_log,
    )


@verbose
def concatenate_epochs(
    epochs_list, add_offset=True, *, on_mismatch="raise", verbose=None
):
    """Concatenate a list of `~mne.Epochs` into one `~mne.Epochs` object.

    .. note:: Unlike `~mne.concatenate_raws`, this function does **not**
              modify any of the input data.

    Parameters
    ----------
    epochs_list : list
        List of `~mne.Epochs` instances to concatenate (in that order).
    add_offset : bool
        If True, a fixed offset is added to the event times from different
        Epochs sets, such that they are easy to distinguish after the
        concatenation.
        If False, the event times are unaltered during the concatenation.
    %(on_mismatch_info)s
    %(verbose)s

        .. versionadded:: 0.24

    Returns
    -------
    epochs : instance of EpochsArray
        The result of the concatenation. All data will be loaded into memory.

    Notes
    -----
    .. versionadded:: 0.9.0
    """
    (
        info,
        data,
        raw_sfreq,
        events,
        event_id,
        tmin,
        tmax,
        metadata,
        baseline,
        selection,
        drop_log,
    ) = _concatenate_epochs(
        epochs_list,
        with_data=True,
        add_offset=add_offset,
        on_mismatch=on_mismatch,
    )
    selection = np.where([len(d) == 0 for d in drop_log])[0]
    out = EpochsArray(
        data=data,
        info=info,
        events=events,
        event_id=event_id,
        tmin=tmin,
        baseline=baseline,
        selection=selection,
        drop_log=drop_log,
        proj=False,
        on_missing="ignore",
        metadata=metadata,
        raw_sfreq=raw_sfreq,
    )
    out.drop_bad()
    return out


@verbose
def average_movements(
    epochs,
    head_pos=None,
    orig_sfreq=None,
    picks=None,
    origin="auto",
    weight_all=True,
    int_order=8,
    ext_order=3,
    destination=None,
    ignore_ref=False,
    return_mapping=False,
    mag_scale=100.0,
    verbose=None,
):
    """Average data using Maxwell filtering, transforming using head positions.

    Parameters
    ----------
    epochs : instance of Epochs
        The epochs to operate on.
    %(head_pos_maxwell)s
    orig_sfreq : float | None
        The original sample frequency of the data (that matches the
        event sample numbers in ``epochs.events``). Can be ``None``
        if data have not been decimated or resampled.
    %(picks_all_data)s
    %(origin_maxwell)s
    weight_all : bool
        If True, all channels are weighted by the SSS basis weights.
        If False, only MEG channels are weighted, other channels
        receive uniform weight per epoch.
    %(int_order_maxwell)s
    %(ext_order_maxwell)s
    %(destination_maxwell_dest)s
    %(ignore_ref_maxwell)s
    return_mapping : bool
        If True, return the mapping matrix.
    %(mag_scale_maxwell)s

        .. versionadded:: 0.13
    %(verbose)s

    Returns
    -------
    evoked : instance of Evoked
        The averaged epochs.

    See Also
    --------
    mne.preprocessing.maxwell_filter
    mne.chpi.read_head_pos

    Notes
    -----
    The Maxwell filtering version of this algorithm is described in [1]_,
    in section V.B "Virtual signals and movement correction", equations
    40-44. For additional validation, see [2]_.

    Regularization has not been added because in testing it appears to
    decrease dipole localization accuracy relative to using all components.
    Fine calibration and cross-talk cancellation, however, could be added
    to this algorithm based on user demand.

    .. versionadded:: 0.11

    References
    ----------
    .. [1] Taulu S. and Kajola M. "Presentation of electromagnetic
           multichannel data: The signal space separation method,"
           Journal of Applied Physics, vol. 97, pp. 124905 1-10, 2005.
    .. [2] Wehner DT, Hämäläinen MS, Mody M, Ahlfors SP. "Head movements
           of children in MEG: Quantification, effects on source
           estimation, and compensation. NeuroImage 40:541–550, 2008.
    """  # noqa: E501
    from .preprocessing.maxwell import (
        _check_destination,
        _check_usable,
        _col_norm_pinv,
        _get_coil_scale,
        _get_mf_picks_fix_mags,
        _get_n_moments,
        _get_sensor_operator,
        _prep_mf_coils,
        _remove_meg_projs_comps,
        _reset_meg_bads,
        _trans_sss_basis,
    )

    if head_pos is None:
        raise TypeError("head_pos must be provided and cannot be None")
    from .chpi import head_pos_to_trans_rot_t

    if not isinstance(epochs, BaseEpochs):
        raise TypeError(f"epochs must be an instance of Epochs, not {type(epochs)}")
    orig_sfreq = epochs.info["sfreq"] if orig_sfreq is None else orig_sfreq
    orig_sfreq = float(orig_sfreq)
    if isinstance(head_pos, np.ndarray):
        head_pos = head_pos_to_trans_rot_t(head_pos)
    trn, rot, t = head_pos
    del head_pos
    _check_usable(epochs, ignore_ref)
    origin = _check_origin(origin, epochs.info, "head")
    recon_trans = _check_destination(destination, epochs.info, "head")

    logger.info(f"Aligning and averaging up to {len(epochs.events)} epochs")
    if not np.array_equal(epochs.events[:, 0], np.unique(epochs.events[:, 0])):
        raise RuntimeError("Epochs must have monotonically increasing events")
    info_to = epochs.info.copy()
    meg_picks, mag_picks, grad_picks, good_mask, _ = _get_mf_picks_fix_mags(
        info_to, int_order, ext_order, ignore_ref
    )
    coil_scale, mag_scale = _get_coil_scale(
        meg_picks, mag_picks, grad_picks, mag_scale, info_to
    )
    mult = _get_sensor_operator(epochs, meg_picks)
    n_channels, n_times = len(epochs.ch_names), len(epochs.times)
    other_picks = np.setdiff1d(np.arange(n_channels), meg_picks)
    data = np.zeros((n_channels, n_times))
    count = 0
    # keep only MEG w/bad channels marked in "info_from"
    info_from = pick_info(info_to, meg_picks[good_mask], copy=True)
    all_coils_recon = _prep_mf_coils(info_to, ignore_ref=ignore_ref)
    all_coils = _prep_mf_coils(info_from, ignore_ref=ignore_ref)
    # remove MEG bads in "to" info
    _reset_meg_bads(info_to)
    # set up variables
    w_sum = 0.0
    n_in, n_out = _get_n_moments([int_order, ext_order])
    S_decomp = 0.0  # this will end up being a weighted average
    last_trans = None
    decomp_coil_scale = coil_scale[good_mask]
    exp = dict(int_order=int_order, ext_order=ext_order, head_frame=True, origin=origin)
    n_in = _get_n_moments(int_order)
    for ei, epoch in enumerate(epochs):
        event_time = epochs.events[epochs._current - 1, 0] / orig_sfreq
        use_idx = np.where(t <= event_time)[0]
        if len(use_idx) == 0:
            trans = info_to["dev_head_t"]["trans"]
        else:
            use_idx = use_idx[-1]
            trans = np.vstack(
                [np.hstack([rot[use_idx], trn[[use_idx]].T]), [[0.0, 0.0, 0.0, 1.0]]]
            )
        loc_str = ", ".join(f"{tr:0.1f}" for tr in (trans[:3, 3] * 1000))
        if last_trans is None or not np.allclose(last_trans, trans):
            logger.info(
                f"    Processing epoch {ei + 1} (device location: {loc_str} mm)"
            )
            reuse = False
            last_trans = trans
        else:
            logger.info(f"    Processing epoch {ei + 1} (device location: same)")
            reuse = True
        epoch = epoch.copy()  # because we operate inplace
        if not reuse:
            S = _trans_sss_basis(exp, all_coils, trans, coil_scale=decomp_coil_scale)
            # Get the weight from the un-regularized version (eq. 44)
            weight = np.linalg.norm(S[:, :n_in])
            # XXX Eventually we could do cross-talk and fine-cal here
            S *= weight
        S_decomp += S  # eq. 41
        epoch[slice(None) if weight_all else meg_picks] *= weight
        data += epoch  # eq. 42
        w_sum += weight
        count += 1
    del info_from
    mapping = None
    if count == 0:
        data.fill(np.nan)
    else:
        data[meg_picks] /= w_sum
        data[other_picks] /= w_sum if weight_all else count
        # Finalize weighted average decomp matrix
        S_decomp /= w_sum
        # Get recon matrix
        # (We would need to include external here for regularization to work)
        exp["ext_order"] = 0
        S_recon = _trans_sss_basis(exp, all_coils_recon, recon_trans)
        if mult is not None:
            S_decomp = mult @ S_decomp
            S_recon = mult @ S_recon
        exp["ext_order"] = ext_order
        # We could determine regularization on basis of destination basis
        # matrix, restricted to good channels, as regularizing individual
        # matrices within the loop above does not seem to work. But in
        # testing this seemed to decrease localization quality in most cases,
        # so we do not provide the option here.
        S_recon /= coil_scale
        # Invert
        pS_ave = _col_norm_pinv(S_decomp)[0][:n_in]
        pS_ave *= decomp_coil_scale.T
        # Get mapping matrix
        mapping = np.dot(S_recon, pS_ave)
        # Apply mapping
        data[meg_picks] = np.dot(mapping, data[meg_picks[good_mask]])
    info_to["dev_head_t"] = recon_trans  # set the reconstruction transform
    evoked = epochs._evoked_from_epoch_data(
        data, info_to, picks, n_events=count, kind="average", comment=epochs._name
    )
    _remove_meg_projs_comps(evoked, ignore_ref)
    logger.info(f"Created Evoked dataset from {count} epochs")
    return (evoked, mapping) if return_mapping else evoked


@verbose
def make_fixed_length_epochs(
    raw,
    duration=1.0,
    preload=False,
    reject_by_annotation=True,
    proj=True,
    overlap=0.0,
    id=1,  # noqa: A002
    verbose=None,
):
    """Divide continuous raw data into equal-sized consecutive epochs.

    Parameters
    ----------
    raw : instance of Raw
        Raw data to divide into segments.
    duration : float
        Duration of each epoch in seconds. Defaults to 1.
    %(preload)s
    %(reject_by_annotation_epochs)s

        .. versionadded:: 0.21.0
    %(proj_epochs)s

        .. versionadded:: 0.22.0
    overlap : float
        The overlap between epochs, in seconds. Must be
        ``0 <= overlap < duration``. Default is 0, i.e., no overlap.

        .. versionadded:: 0.23.0
    id : int
        The id to use (default 1).

        .. versionadded:: 0.24.0
    %(verbose)s

    Returns
    -------
    epochs : instance of Epochs
        Segmented data.

    Notes
    -----
    .. versionadded:: 0.20
    """
    events = make_fixed_length_events(raw, id=id, duration=duration, overlap=overlap)
    delta = 1.0 / raw.info["sfreq"]
    return Epochs(
        raw,
        events,
        event_id=[id],
        tmin=0,
        tmax=duration - delta,
        baseline=None,
        preload=preload,
        reject_by_annotation=reject_by_annotation,
        proj=proj,
        verbose=verbose,
    )