text: """Default settings for data processing and analysis. Do not edit this file,
  but instead create a new configuration changing only the settings you need to
  alter for your specific analysis.
  """
  
  import importlib
  import pathlib
  import functools
  import os
  import pdb
  import traceback
  import sys
  import copy
  import logging
  import time
  from typing import Optional, Union, Iterable, List, Tuple, Dict, Callable
  if sys.version_info >= (3, 8):
      from typing import Literal, TypedDict
  else:
      from typing_extensions import Literal, TypedDict
  
  import coloredlogs
  import numpy as np
  import pandas as pd
  from openpyxl import load_workbook
  import json_tricks
  import mne
  import mne_bids
  from mne_bids import BIDSPath, read_raw_bids
  
  PathLike = Union[str, pathlib.Path]
  
  
  study_name: str = ''
  """
  Specify the name of your study. It will be used to populate filenames for
  saving the analysis results.
  
  ???+ example "Example"
      ```python
  
      study_name = 'my-study'
      ```
  """
  
  bids_root: Optional[PathLike] = None
  """
  Specify the BIDS root directory. Pass an empty string or ```None`` to use
  the value specified in the ``BIDS_ROOT`` environment variable instead.
  Raises an exception if the BIDS root has not been specified.
  
  ???+ example "Example"
      ``` python
      bids_root = '/path/to/your/bids_root'  # Use this to specify a path here.
      bids_root = None  # Make use of the ``BIDS_ROOT`` environment variable.
      ```
  """
  
  deriv_root: Optional[PathLike] = None
  """
  The root of the derivatives directory in which the pipeline will store
  the processing results. If ``None``, this will be
  ``derivatives/mne-bids-pipeline`` inside the BIDS root.
  
  Note: Note
      If specified and you wish to run the source analysis steps, you must
      set [`subjects_dir`][config.subjects_dir] as well.
  """
  
  subjects_dir: Optional[PathLike] = None
  """
  Path to the directory that contains the FreeSurfer reconstructions of all
  subjects. Specifically, this defines the ``SUBJECTS_DIR`` that is used by
  FreeSurfer.
  
  - When running the ``freesurfer`` processing step to create the
    reconstructions from anatomical scans in the BIDS dataset, the
    output will be stored in this directory.
  - When running the source analysis steps, we will look for the surfaces in this
    directory and also store the BEM surfaces there.
  
  If ``None``, this will default to
  [`bids_root`][config.bids_root]`/derivatives/freesurfer/subjects`.
  
  Note: Note
      This setting is required if you specify [`deriv_root`][config.deriv_root]
      and want to run the source analysis steps.
  """
  
  interactive: bool = False
  """
  If True, the scripts will provide some interactive elements, such as
  figures. If running the scripts from a notebook or Spyder,
  run `%matplotlib qt` in the command line to open the figures in a separate
  window.
  """
  
  sessions: Union[List, Literal['all']] = 'all'
  """
  The sessions to process.
  """
  
  task: str = ''
  """
  The task to process.
  """
  
  runs: Union[Iterable, Literal['all']] = 'all'
  """
  The runs to process.
  """
  
  exclude_runs: Optional[Dict[str, List[str]]] = None
  """
  Specify runs to exclude from analysis, for each participant individually.
  
  ???+ example "Example"
      ```python
      exclude_runs = None  # Include all runs.
      exclude_runs = {'01': ['02']}  # Exclude run 02 of subject 01.
      ```
  
  ???+ info "Good Practice / Advice"
      Keep track of the criteria leading you to exclude
      a run (e.g. too many movements, missing blocks, aborted experiment,
      did not understand the instructions, etc.).
  """
  
  crop_runs: Optional[Tuple[float, float]] = None
  """
  Crop the raw data of each run to the specified time interval ``[tmin, tmax]``,
  in seconds. The runs will be cropped before Maxwell or frequency filtering is
  applied. If ``None``, do not crop the data.
  """
  
  acq: Optional[str] = None
  """
  The BIDS `acquisition` entity.
  """
  
  proc: Optional[str] = None
  """
  The BIDS `processing` entity.
  """
  
  rec: Optional[str] = None
  """
  The BIDS `recording` entity.
  """
  
  space: Optional[str] = None
  """
  The BIDS `space` entity.
  """
  
  subjects: Union[Iterable[str], Literal['all']] = 'all'
  """
  Subjects to analyze. If ``'all'``, include all subjects. To only
  include a subset of subjects, pass a list of their identifiers. Even
  if you plan on analyzing only a single subject, pass their identifier
  as a list.
  
  Please note that if you intend to EXCLUDE only a few subjects, you
  should consider setting ``subjects = 'all'`` and adding the
  identifiers of the excluded subjects to ``exclude_subjects`` (see next
  section).
  
  ???+ example "Example"
      ```python
      subjects = 'all'  # Include all subjects.
      subjects = ['05']  # Only include subject 05.
      subjects = ['01', '02']  # Only include subjects 01 and 02.
      ```
  """
  
  exclude_subjects: Iterable[str] = []
  """
  Specify subjects to exclude from analysis. The MEG empty-room mock-subject
  is automatically excluded from regular analysis.
  
  ???+ info "Good Practice / Advice"
      Keep track of the criteria leading you to exclude
      a participant (e.g. too many movements, missing blocks, aborted experiment,
      did not understand the instructions, etc, ...)
      The ``emptyroom`` subject will be excluded automatically.
  """
  
  process_er: bool = False
  """
  Whether to apply the same pre-processing steps to the empty-room data as
  to the experimental data (up until including frequency filtering). This
  is required if you wish to use the empty-room recording to estimate noise
  covariance (via ``noise_cov='emptyroom'``). The empty-room recording
  corresponding to the processed experimental data will be retrieved
  automatically.
  """
  
  ch_types: Iterable[Literal['meg', 'mag', 'grad', 'eeg']] = []
  """
  The channel types to consider.
  
  !!! info
      Currently, MEG and EEG data cannot be processed together.
  
  ???+ example "Example"
      ```python
      # Use EEG channels:
      ch_types = ['eeg']
  
      # Use magnetometer and gradiometer MEG channels:
      ch_types = ['mag', 'grad']
  
      # Currently does not work and will raise an error message:
      ch_types = ['meg', 'eeg']
      ```
  """
  
  data_type: Optional[Literal['meg', 'eeg']] = None
  """
  The BIDS data type.
  
  For MEG recordings, this will usually be 'meg'; and for EEG, 'eeg'.
  However, if your dataset contains simultaneous recordings of MEG and EEG,
  stored in a single file, you will typically need to set this to 'meg'.
  If ``None``, we will assume that the data type matches the channel type.
  
  ???+ example "Example"
      The dataset contains simultaneous recordings of MEG and EEG, and we only
      wish to process the EEG data, which is stored inside the MEG files:
  
      ```python
      ch_types = ['eeg']
      data_type = 'eeg'
      ```
  
      The dataset contains simultaneous recordings of MEG and EEG, and we only
      wish to process the gradiometer data:
  
      ```python
      ch_types = ['grad']
      data_type = 'meg'  # or data_type = None
      ```
  
      The dataset contains only EEG data:
  
      ```python
      ch_types = ['eeg']
      data_type = 'eeg'  # or data_type = None
      ```
  """
  
  eog_channels: Optional[Iterable[str]] = None
  """
  Specify EOG channels to use, or create virtual EOG channels.
  
  Allows the specification of custom channel names that shall be used as
  (virtual) EOG channels. For example, say you recorded EEG **without** dedicated
  EOG electrodes, but with some EEG electrodes placed close to the eyes, e.g.
  Fp1 and Fp2. These channels can be expected to have captured large quantities
  of ocular activity, and you might want to use them as "virtual" EOG channels,
  while also including them in the EEG analysis. By default, MNE won't know that
  these channels are suitable for recovering EOG, and hence won't be able to
  perform tasks like automated blink removal, unless a "true" EOG sensor is
  present in the data as well. Specifying channel names here allows MNE to find
  the respective EOG signals based on these channels.
  
  You can specify one or multiple channel names. Each will be treated as if it
  were a dedicated EOG channel, without excluding it from any other analyses.
  
  If ``None``, only actual EOG channels will be used for EOG recovery.
  
  If there are multiple actual EOG channels in your data, and you only specify
  a subset of them here, only this subset will be used during processing.
  
  ???+ example "Example"
      Treat ``Fp1`` as virtual EOG channel:
      ```python
      eog_channels = ['Fp1']
      ```
  
      Treat ``Fp1`` and ``Fp2`` as virtual EOG channels:
      ```python
      eog_channels = ['Fp1', 'Fp2']
      ```
  """
  
  eeg_bipolar_channels: Optional[Dict[str, Tuple[str, str]]] = None
  """
  Combine two channels into a bipolar channel, whose signal is the **difference**
  between the two combined channels, and add it to the data.
  A typical use case is the combination of two EOG channels – for example, a
  left and a right horizontal EOG – into a single, bipolar EOG channel. You need
  to pass a dictionary whose **keys** are the name of the new bipolar channel you
  wish to create, and whose **values** are tuples consisting of two strings: the
  name of the channel acting as anode and the name of the channel acting as
  cathode, i.e. `{'ch_name': ('anode', 'cathode')}`. You can request
  to construct more than one bipolar channel by specifying multiple key/value
  pairs. See the examples below.
  
  Can also be `None` if you do not want to create bipolar channels.
  
  Note: Note
      The channels used to create the bipolar channels are **not** automatically
      dropped from the data. To drop channels, set `drop_channels`.
  
  ???+ example "Example"
      Combine the existing channels `HEOG_left` and `HEOG_right` into a new,
      bipolar channel, `HEOG`:
      ```python
      eeg_add_bipolar_channels = {'HEOG': ('HEOG_left', 'HEOG_right')}
      ```
  
      Create two bipolar channels, `HEOG` and `VEOG`:
      ```python
      eeg_add_bipolar_channels = {'HEOG': ('HEOG_left', 'HEOG_right'),
                                  'VEOG': ('VEOG_lower', 'VEOG_upper')}
      ```
  """
  
  eeg_reference: Union[Literal['average'], str, Iterable['str']] = 'average'
  """
  The EEG reference to use. If ``average``, will use the average reference,
  i.e. the average across all channels. If a string, must be the name of a single
  channel. To use multiple channels as reference, set to a list of channel names.
  
  ???+ example "Example"
      Use the average reference:
      ```python
      eeg_reference = 'average'
      ```
  
      Use the `P9` channel as reference:
      ```python
      eeg_reference = 'P9'
      ```
  
      Use the average of the `P9` and `P10` channels as reference:
      ```python
      eeg_reference = ['P9', 'P10']
      ```
  """
  
  eeg_template_montage: Optional[str] = None
  """
  In situations where you wish to process EEG data and no individual
  digitization points (measured channel locations) are available, you can apply
  a "template" montage. This means we will assume the EEG cap was placed
  either according to an international system like 10/20, or as suggested by
  the cap manufacturers in their respective manual.
  
  Please be aware that the actual cap placement most likely deviated somewhat
  from the template, and, therefore, source reconstruction may be impaired.
  
  If ``None``, do not apply a template montage. If a string, must be the
  name of a built-in template montage in MNE-Python.
  You can find an overview of supported template montages at
  https://mne.tools/stable/generated/mne.channels.make_standard_montage.html
  
  ???+ example "Example"
      Do not apply template montage:
      ```python
      eeg_template_montage = None
      ```
  
      Apply 64-channel Biosemi 10/20 template montage:
      ```python
      eeg_template_montage = 'biosemi64'
      ```
  """
  
  drop_channels: Iterable[str] = []
  """
  Names of channels to remove from the data. This can be useful, for example,
  if you have added a new bipolar channel via `eeg_bipolar_channels` and now wish
  to remove the anode, cathode, or both.
  
  ???+ example "Example"
      Exclude channels `Fp1` and `Cz` from processing:
      ```python
      drop_channels = ['Fp1', 'Cz]
      ```
  """
  
  analyze_channels: Union[Literal['all'], Iterable['str']] = 'all'
  """
  The names of the channels to analyze during ERP/ERF and time-frequency analysis
  steps. For certain paradigms, e.g. EEG ERP research, it is common to contrain
  sensor-space analysis to only a few specific sensors. If `'all'`, do not
  exclude any channels (except for those selected for removal via the
  `drop_channels` setting). The constraint will be applied to all sensor-level
  analyses after the preprocessing stage, but not to the preprocessing stage
  itself, nor to the source analysis stage.
  
  ???+ example "Example"
      Only use channel `Pz` for ERP, evoked contrasts, time-by-time
      decoding, and time-frequency analysis:
      ```python
      analyze_channels = ['Pz']
      ```
  """
  
  ###############################################################################
  # BREAK DETECTION
  # ---------------
  
  find_breaks: bool = False
  """
  During an experimental run, the recording might be interrupted by breaks of
  various durations, e.g. to allow the participant to stretch, blink, and swallow
  freely. During these periods, large-scale artifacts are often picked up by the
  recording system. These artifacts can impair certain stages of processing, e.g.
  the peak-detection algorithms we use to find EOG and ECG activity. In some
  cases, even the bad channel detection algorithms might not function optimally.
  It is therefore advisable to mark such break periods for exclusion at early
  processing stages.
  
  If `True`, try to mark breaks by finding segments of the data where no
  exprimental events have occurred. This will then add annotations with the
  description `BAD_break` to the continuous data, causing these segments to be
  ignored in all following processing steps.
  
  ???+ example "Example"
      Automatically find break periods, and annotate them as `BAD_break`.
      ```python
      find_breaks = True
      ```
  
      Disable break detection.
      ```python
      find_breaks = False
      ```
  """
  
  min_break_duration: float = 15.
  """
  The minimal duration (in seconds) of a data segment without any experimental
  events for it to be considered a "break". Note that the minimal duration of the
  generated `BAD_break` annotation will typically be smaller than this, as by
  default, the annotation will not extend across the entire break.
  See [`t_break_annot_start_after_previous_event`][config.t_break_annot_start_after_previous_event]
  and [`t_break_annot_stop_before_next_event`][config.t_break_annot_stop_before_next_event]
  to control this behavior.
  
  ???+ example "Example"
      Periods between two consecutive experimental events must span at least
      `15` seconds for this period to be considered a "break".
      ```python
      min_break_duration = 15.
      ```
  """
  
  t_break_annot_start_after_previous_event: float = 5.
  """
  Once a break of at least [`min_break_duration`][config.min_break_duration]
  seconds has been discovered, we generate a `BAD_break` annotation that does not
  necessarily span the entire break period. Instead, you will typically want to
  start it some time after the last event before the break period, as to not
  unnecessarily discard brain activity immediately following that event.
  
  This parameter controls how much time (in seconds) should pass after the last
  pre-break event before we start annotating the following segment of the break
  period as bad.
  
  ???+ example "Example"
      Once a break period has been detected, add a `BAD_break` annotation to it,
      starting `5` seconds after the latest pre-break event.
      ```python
      t_break_annot_start_after_previous_event = 5.
      ```
  
      Start the `BAD_break` annotation immediately after the last pre-break
      event.
      ```python
      t_break_annot_start_after_previous_event = 0.
      ```
  """
  
  t_break_annot_stop_before_next_event: float = 5.
  """
  Similarly to how
  [`t_break_annot_start_after_previous_event`][config.t_break_annot_start_after_previous_event]
  controls the "gap" between beginning of the break period and `BAD_break`
  annotation onset,  this parameter controls how far the annotation should extend
  toward the first experimental event immediately following the break period
  (in seconds). This can help not to waste a post-break trial by marking its
  pre-stimulus period as bad.
  
  ???+ example "Example"
      Once a break period has been detected, add a `BAD_break` annotation to it,
      starting `5` seconds after the latest pre-break event.
      ```python
      t_break_annot_start_after_previous_event = 5.
      ```
  
      Start the `BAD_break` annotation immediately after the last pre-break
      event.
      ```python
      t_break_annot_start_after_previous_event = 0.
      ```
  """
  
  ###############################################################################
  # MAXWELL FILTER PARAMETERS
  # -------------------------
  # done in 01-import_and_maxfilter.py
  
  find_flat_channels_meg: bool = False
  """
  Auto-detect "flat" channels (i.e. those with unusually low variability) and
  mark them as bad.
  """
  
  find_noisy_channels_meg: bool = False
  """
  Auto-detect "noisy" channels and mark them as bad.
  """
  
  use_maxwell_filter: bool = False
  """
  Whether or not to use Maxwell filtering to preprocess the data.
  
  warning:
      If the data were recorded with internal active compensation (MaxShield),
      they need to be run through Maxwell filter to avoid distortions.
      Bad channels need to be set through BIDS channels.tsv and / or via the
      ``find_flat_channels_meg`` and ``find_noisy_channels_meg`` options above
      before applying Maxwell filter.
  """
  
  mf_st_duration: Optional[float] = None
  """
  There are two kinds of maxfiltering: SSS (signal space separation) and tSSS
  (temporal signal space separation)
  (see [Taulu et al., 2004](http://cds.cern.ch/record/709081/files/0401166.pdf)).
  
  If not None, apply spatiotemporal SSS (tSSS) with specified buffer
  duration (in seconds). MaxFilter™'s default is 10.0 seconds in v2.2.
  Spatiotemporal SSS acts as implicitly as a high-pass filter where the
  cut-off frequency is 1/st_dur Hz. For this (and other) reasons, longer
  buffers are generally better as long as your system can handle the
  higher memory usage. To ensure that each window is processed
  identically, choose a buffer length that divides evenly into your data.
  Any data at the trailing edge that doesn't fit evenly into a whole
  buffer window will be lumped into the previous buffer.
  
  ???+ info "Good Practice / Advice"
      If you are interested in low frequency activity (<0.1Hz), avoid using
      tSSS and set ``mf_st_duration`` to ``None``.
  
      If you are interested in low frequency above 0.1 Hz, you can use the
      default ``mf_st_duration`` to 10 s, meaning it acts like a 0.1 Hz
      high-pass filter.
  
  ???+ example "Example"
      ```python
      mf_st_duration = None
      mf_st_duration = 10.  # to apply tSSS with 0.1Hz highpass filter.
      ```
  """
  
  mf_head_origin = 'auto'
  """
  ``mf_head_origin`` : array-like, shape (3,) | 'auto'
  Origin of internal and external multipolar moment space in meters.
  If 'auto', it will be estimated from headshape points.
  If automatic fitting fails (e.g., due to having too few digitization
  points), consider separately calling the fitting function with different
  options or specifying the origin manually.
  
  ???+ example "Example"
      ```python
      mf_head_origin = 'auto'
      ```
  """
  
  mf_reference_run: Optional[str] = None
  """
  Despite all possible care to avoid movements in the MEG, the participant
  will likely slowly drift down from the Dewar or slightly shift the head
  around in the course of the recording session. Hence, to take this into
  account, we are realigning all data to a single position. For this, you need
  to define a reference run (typically the one in the middle of
  the recording session).
  
  Which run to take as the reference for adjusting the head position of all
  runs. If ``None``, pick the first run.
  
  ???+ example "Example"
      ```python
      mf_reference_run = '01'  # Use run "01"
      ```
  """
  
  mf_cal_fname: Optional[str] = None
  """
  warning:
       This parameter should only be used for BIDS datasets that don't store
       the fine-calibration file
       [according to BIDS](https://bids-specification.readthedocs.io/en/stable/99-appendices/06-meg-file-formats.html#cross-talk-and-fine-calibration-files).
  Path to the Maxwell Filter calibration file. If None the recommended
  location is used.
  ???+ example "Example"
      ```python
      mf_cal_fname = '/path/to/your/file/calibration_cal.dat'
      ```
  """
  
  mf_ctc_fname: Optional[str] = None
  """
  Path to the Maxwell Filter cross-talk file. If None the recommended
  location is used.
  warning:
       This parameter should only be used for BIDS datasets that don't store
       the cross-talk file
       [according to BIDS](https://bids-specification.readthedocs.io/en/stable/99-appendices/06-meg-file-formats.html#cross-talk-and-fine-calibration-files).
  ???+ example "Example"
      ```python
      mf_ctc_fname = '/path/to/your/file/crosstalk_ct.fif'
      ```
  """
  
  ###############################################################################
  # STIMULATION ARTIFACT
  # --------------------
  # used in 01-import_and_maxfilter.py
  
  fix_stim_artifact: bool = False
  """
  Apply interpolation to fix stimulation artifact.
  
  ???+ example "Example"
      ```python
      fix_stim_artifact = False
      ```
  """
  
  stim_artifact_tmin: float = 0.
  """
  Start time of the interpolation window in seconds.
  
  ???+ example "Example"
      ```python
      stim_artifact_tmin = 0.  # on stim onset
      ```
  """
  
  stim_artifact_tmax: float = 0.01
  """
  End time of the interpolation window in seconds.
  
  ???+ example "Example"
      ```python
      stim_artifact_tmax = 0.01  # up to 10ms post-stimulation
      ```
  """
  
  ###############################################################################
  # FREQUENCY FILTERING & RESAMPLING
  # --------------------------------
  # done in 02-frequency_filter.py
  
  l_freq: Optional[float] = None
  """
  The low-frequency cut-off in the highpass filtering step.
  Keep it None if no highpass filtering should be applied.
  """
  
  h_freq: Optional[float] = 40.
  """
  The high-frequency cut-off in the lowpass filtering step.
  Keep it None if no lowpass filtering should be applied.
  """
  
  resample_sfreq: Optional[float] = None
  """
  Specifies at which sampling frequency the data should be resampled.
  If `None`, then no resampling will be done.
  
  ???+ example "Example"
      ```python
      resample_sfreq = None  # no resampling
      resample_sfreq = 500  # resample to 500Hz
      ```
  """
  
  ###############################################################################
  # DECIMATION
  # ----------
  
  decim: int = 1
  """
  Says how much to decimate data at the epochs level.
  It is typically an alternative to the `resample_sfreq` parameter that
  can be used for resampling raw data. ``1`` means no decimation.
  
  ???+ info "Good Practice / Advice"
      Decimation requires to lowpass filtered the data to avoid aliasing.
      Note that using decimation is much faster than resampling.
  
  ???+ example "Example"
      ```python
      decim = 1  # no decimation
      decim = 4  # decimate by 4 ie devide sampling frequency by 4
      ```
  """
  
  
  ###############################################################################
  # RENAME EXPERIMENTAL EVENTS
  # --------------------------
  
  rename_events: dict = dict()
  """
  A dictionary specifying which events in the BIDS dataset to rename upon
  loading, and before processing begins.
  
  Pass an empty dictionary to not perform any renaming.
  
  ???+ example "Example"
      Rename ``audio_left`` in the BIDS dataset to ``audio/left`` in the
      pipeline:
      ```python
      rename_events = {'audio_left': 'audio/left'}
      ```
  """
  
  on_rename_missing_events: Literal['warn', 'raise'] = 'raise'
  """
  How to handle the situation where you specified an event to be renamed via
  ``rename_events``, but this particular event is not present in the data. By
  default, we will raise an exception to avoid accidental mistakes due to typos;
  however, if you're sure what you're doing, you may change this to ``'warn'``
  to only get a warning instead.
  """
  
  ###############################################################################
  # HANDLING OF REPEATED EVENTS
  # ---------------------------
  
  event_repeated: Literal['error', 'drop', 'merge'] = 'error'
  """
  How to handle repeated events. We call events "repeated" if more than one event
  occurred at the exact same time point. Currently, MNE-Python cannot handle
  this situation gracefully when trying to create epochs, and will throw an
  error. To only keep the event of that time point ("first" here referring to
  the order that events appear in `*_events.tsv`), pass `'drop'`. You can also
  request to create a new type of event by merging repeated events by setting
  this to `'merge'`.
  
  warning:
      The `'merge'` option is entirely untested in the MNE BIDS Pipeline as of
      April 1st, 2021.
  """
  
  ###############################################################################
  # EPOCHING
  # --------
  
  epochs_metadata_tmin: Optional[float] = None
  """
  The beginning of the time window for metadata generation, in seconds,
  relative to the time-locked event of the respective epoch. This may be less
  than or larger than the epoch's first time point. If ``None``, use the first
  time point of the epoch.
  """
  
  epochs_metadata_tmax: Optional[float] = None
  """
  Same as ``epochs_metadata_tmin``, but specifying the **end** of the time
  window for metadata generation.
  """
  
  epochs_metadata_keep_first: Optional[Iterable[str]] = None
  """
  Event groupings using hierarchical event descriptors (HEDs) for which to store
  the time of the **first** occurrence of any event of this group in a new column
  with the group name, and the **type** of that event in a column named after the
  group, but with a ``first_`` prefix. If ``None`` (default), no event
  aggregation will take place and no new columns will be created.
  
  ???+ example "Example"
      Assume you have two response events types, ``response/left`` and
      ``response/right``; in some trials, both responses occur, because the
      participant pressed both buttons. Now, you want to keep the first response
      only. To achieve this, set
      ```python
      epochs_metadata_keep_first = ['response']
      ```
      This will add two new columns to the metadata: ``response``, indicating
      the **time** relative to the time-locked event; and ``first_response``,
      depicting the **type** of event (``'left'`` or ``'right'``).
  
      You may also specify a grouping for multiple event types:
      ```python
      epochs_metadata_keep_first = ['response', 'stimulus']
      ```
      This will add the columns ``response``, ``first_response``, ``stimulus``,
      and ``first_stimulus``.
  """
  
  epochs_metadata_keep_last: Optional[Iterable[str]] = None
  """
  Same as ``epochs_metadata_keep_first``, but for keeping the **last**
  occurrence of matching event types. The columns indicating the event types
  will be named with a ``last_`` instead of a ``first_`` prefix.
  """
  
  conditions: Optional[Union[Iterable[str], Dict[str, str]]] = None
  """
  The time-locked events based on which to create evoked responses.
  This can either be name of the experimental condition as specified in the
  BIDS ``*_events.tsv`` file; or the name of condition *groups*, if the condition
  names contain the (MNE-specific) group separator, ``/``. See the [Subselecting
  epochs tutorial](https://mne.tools/stable/auto_tutorials/epochs/plot_10_epochs_overview.html#subselecting-epochs)
  for more information.
  
  Passing a dictionary allows to assign a name to map a complex condition name
  (value) to a more legible one (value).
  
  This is a **required** parameter in the configuration file, unless you are
  processing resting-state data. If left as `None` and [`task`][config.task]
  is not `'rest'`, we will raise an error.
  
  ???+ example "Example"
      Specifying conditions as lists of strings:
      ```python
      conditions = ['auditory/left', 'visual/left']
      conditions = ['auditory/left', 'auditory/right']
      conditions = ['auditory']  # All "auditory" conditions (left AND right)
      conditions = ['auditory', 'visual']
      conditions = ['left', 'right']
      conditions = None  # for a resting-state analysis
      ```
      Pass a dictionary to define a mapping:
      ```python
      conditions = {'simple_name': 'complex/condition/with_subconditions'}
      conditions = {'correct': 'response/correct',
                    'incorrect': 'response/incorrect'}
  """
  
  epochs_tmin: float = -0.2
  """
  The beginning of an epoch, relative to the respective event, in seconds.
  
  ???+ example "Example"
      ```python
      epochs_tmin = -0.2  # 200 ms before event onset
      ```
  """
  
  epochs_tmax: float = 0.5
  """
  The end of an epoch, relative to the respective event, in seconds.
  ???+ example "Example"
      ```python
      epochs_tmax = 0.5  # 500 ms after event onset
      ```
  """
  
  rest_epochs_duration: Optional[float] = None
  """
  Duration of epochs in seconds.
  """
  
  rest_epochs_overlap: Optional[float] = None
  """
  Overlap between epochs in seconds. This is used if the task is ``'rest'``
  and when the annotations do not contain any stimulation or behavior events.
  """
  
  baseline: Optional[Tuple[Optional[float], Optional[float]]] = (None, 0)
  """
  Specifies which time interval to use for baseline correction of epochs;
  if ``None``, no baseline correction is applied.
  
  ???+ example "Example"
      ```python
      baseline = (None, 0)  # beginning of epoch until time point zero
      ```
  """
  
  contrasts: Iterable[Tuple[str, str]] = []
  """
  The conditions to contrast via a subtraction of ERPs / ERFs. Each tuple
  in the list corresponds to one contrast. The condition names must be
  specified in ``conditions`` above. Pass an empty list to avoid calculation
  of contrasts.
  
  ???+ example "Example"
      Contrast the "left" and the "right" conditions by calculating
      ``left - right`` at every time point of the evoked responses:
      ```python
      conditions = ['left', 'right']
      contrasts = [('left', 'right')]  # Note we pass a tuple inside the list!
      ```
  
      Contrast the "left" and the "right" conditions within the "auditory" and
      the "visual" modality, and "auditory" vs "visual" regardless of side:
      ```python
      conditions = ['auditory/left', 'auditory/right',
                    'visual/left', 'visual/right']
      contrasts = [('auditory/left', 'auditory/right'),
                   ('visual/left', 'visual/right'),
                   ('auditory', 'visual')]
      ```
  """
  
  ###############################################################################
  # ARTIFACT REMOVAL
  # ----------------
  #
  # You can choose between ICA and SSP to remove eye and heart artifacts.
  # SSP: https://mne-tools.github.io/stable/auto_tutorials/plot_artifacts_correction_ssp.html?highlight=ssp # noqa
  # ICA: https://mne-tools.github.io/stable/auto_tutorials/plot_artifacts_correction_ica.html?highlight=ica # noqa
  # if you choose ICA, run scripts 5a and 6a
  # if you choose SSP, run scripts 5b and 6b
  #
  # Currently you cannot use both.
  
  spatial_filter: Optional[Literal['ssp', 'ica']] = None
  """
  Whether to use a spatial filter to detect and remove artifacts. The BIDS
  Pipeline offers the use of signal-space projection (SSP) and independent
  component analysis (ICA).
  
  Use `'ssp'` for SSP, `'ica'` for ICA, and `None` if you do not wish to apply
  a spatial filter for artifact removal.
  
  The Pipeline will try to automatically discover EOG and ECG artifacts. For SSP,
  it will then produce projection vectors that remove ("project out") these
  artifacts from the data. For ICA, the independent components related to
  EOG and ECG activity will be omitted during the signal reconstruction step in
  order to remove the artifacts. The ICA procedure can be configured in various
  ways using the configuration options you can find below.
  """
  
  # Rejection based on SSP
  # ~~~~~~~~~~~~~~~~~~~~~~
  
  
  n_proj_eog: Dict[str, float] = dict(n_mag=1, n_grad=1, n_eeg=1)
  """
  Number of SSP vectors to create for EOG artifacts for each channel type.
  """
  
  n_proj_ecg: Dict[str, float] = dict(n_mag=1, n_grad=1, n_eeg=1)
  """
  Number of SSP vectors to create for ECG artifacts for each channel type.
  """
  
  ecg_proj_from_average: bool = True
  """
  Whether to calculate the ECG projection vectors based on the the averaged or
  on individual ECG epochs.
  """
  
  eog_proj_from_average: bool = True
  """
  Whether to calculate the EOG projection vectors based on the the averaged or
  on individual EOG epochs.
  """
  
  ssp_reject_ecg: Optional[
      Union[
          Dict[str, float],
          Literal['autoreject_global']
      ]
  ] = None
  """
  Peak-to-peak amplitude limits of the ECG epochs to exclude from SSP fitting.
  This allows you to remove strong transient artifacts, which could negatively
  affect SSP performance.
  
  The pipeline will automatically try to detect ECG artifacts in
  your data, and remove them via SSP. For this to work properly, it is
  recommended to **not** specify rejection thresholds for ECG channels here –
  otherwise, SSP won't be able to "see" these artifacts.
  ???+ example "Example"
      ```python
      ssp_reject_ecg = {'grad': 10e-10, 'mag': 20e-12, 'eeg': 400e-6}
      ssp_reject_ecg = {'grad': 15e-10}
      ssp_reject_ecg = None
      ```
  """
  
  ssp_reject_eog: Optional[
      Union[
          Dict[str, float],
          Literal['autoreject_global']
      ]
  ] = None
  """
  Peak-to-peak amplitude limits of the EOG epochs to exclude from SSP fitting.
  This allows you to remove strong transient artifacts, which could negatively
  affect SSP performance.
  
  The pipeline will automatically try to detect EOG artifacts in
  your data, and remove them via SSP. For this to work properly, it is
  recommended to **not** specify rejection thresholds for EOG channels here –
  otherwise, SSP won't be able to "see" these artifacts.
  ???+ example "Example"
      ```python
      ssp_reject_eog = {'grad': 10e-10, 'mag': 20e-12, 'eeg': 400e-6}
      ssp_reject_eog = {'grad': 15e-10}
      ssp_reject_eog = None
      ```
  """
  
  
  # Rejection based on ICA
  # ~~~~~~~~~~~~~~~~~~~~~~
  
  
  ica_reject: Optional[Dict[str, float]] = None
  """
  Peak-to-peak amplitude limits to exclude epochs from ICA fitting.
  
  This allows you to remove strong transient artifacts, which could negatively
  affect ICA performance.
  
  This will also be applied to ECG and EOG epochs created during preprocessing.
  
  The BIDS Pipeline will automatically try to detect EOG and ECG artifacts in
  your data, and remove them. For this to work properly, it is recommended
  to **not** specify rejection thresholds for EOG and ECG channels here –
  otherwise, ICA won't be able to "see" these artifacts.
  
  If `None` (default), do not apply artifact rejection. If a dictionary,
  manually specify peak-to-peak rejection thresholds (see examples).
  
  ???+ example "Example"
      ```python
      ica_reject = {'grad': 10e-10, 'mag': 20e-12, 'eeg': 400e-6}
      ica_reject = {'grad': 15e-10}
      ica_reject = None  # no rejection
      ```
  """
  
  ica_algorithm: Literal['picard', 'fastica', 'extended_infomax'] = 'picard'
  """
  The ICA algorithm to use.
  """
  
  ica_l_freq: Optional[float] = 1.
  """
  The cutoff frequency of the high-pass filter to apply before running ICA.
  Using a relatively high cutoff like 1 Hz will remove slow drifts from the
  data, yielding improved ICA results. Must be set to 1 Hz or above.
  
  Set to ``None`` to not apply an additional high-pass filter.
  
  Note: Note
        The filter will be applied to raw data which was already filtered
        according to the ``l_freq`` and ``h_freq`` settings. After filtering, the
        data will be epoched, and the epochs will be submitted to ICA.
  
  !!! info
      The Pipeline will only allow you to perform ICA on data that has been
      high-pass filtered with a 1 Hz cutoff or higher. This is a conscious,
      opinionated (but partially data-driven) decision made by the developers.
      If you have reason to challenge this behavior, please get in touch with
      us so we can discuss.
  """
  
  ica_max_iterations: int = 500
  """
  Maximum number of iterations to decompose the data into independent
  components. A low number means to finish earlier, but the consequence is
  that the algorithm may not have finished converging. To ensure
  convergence, pick a high number here (e.g. 3000); yet the algorithm will
  terminate as soon as it determines that is has successfully converged, and
  not necessarily exhaust the maximum number of iterations. Note that the
  default of 200 seems to be sufficient for Picard in many datasets, because
  it converges quicker than the other algorithms; but e.g. for FastICA, this
  limit may be too low to achieve convergence.
  """
  
  ica_n_components: Optional[Union[float, int]] = 0.8
  """
  MNE conducts ICA as a sort of a two-step procedure: First, a PCA is run
  on the data (trying to exclude zero-valued components in rank-deficient
  data); and in the second step, the principal componenets are passed
  to the actual ICA. You can select how many of the total principal
  components to pass to ICA – it can be all or just a subset. This determines
  how many independent components to fit, and can be controlled via this
  setting.
  
  If int, specifies the number of principal components that are passed to the
  ICA algorithm, which will be the number of independent components to
  fit. It must not be greater than the rank of your data (which is typically
  the number of channels, but may be less in some cases).
  
  If float between 0 and 1, all principal components with cumulative
  explained variance less than the value specified here will be passed to
  ICA.
  
  If ``None``, **all** principal components will be used.
  
  This setting may drastically alter the time required to compute ICA.
  """
  
  ica_decim: Optional[int] = None
  """
  The decimation parameter to compute ICA. If 5 it means
  that 1 every 5 sample is used by ICA solver. The higher the faster
  it is to run but the less data you have to compute a good ICA. Set to
  ``1`` or ``None`` to not perform any decimation.
  """
  
  ica_ctps_ecg_threshold: float = 0.1
  """
  The threshold parameter passed to `find_bads_ecg` method.
  """
  
  ica_eog_threshold: float = 3.0
  """
  The threshold to use during automated EOG classification. Lower values mean
  that more ICs will be identified as EOG-related. If too low, the
  false-alarm rate increases dramatically.
  """
  
  
  # Rejection based on peak-to-peak amplitude
  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  
  reject: Optional[
      Union[Dict[str, float],
            Literal['autoreject_global']]
  ] = None
  """
  Peak-to-peak amplitude limits to mark epochs as bad. This allows you to remove
  epochs with strong transient artifacts.
  
  If `None` (default), do not apply artifact rejection. If a dictionary,
  manually specify rejection thresholds (see examples).  If
  `'autoreject_global'`, use [`autoreject`](https://autoreject.github.io) to find
  suitable "global" rejection thresholds for each channel type, i.e. `autoreject`
  will generate a dictionary with (hopefully!) optimal thresholds for each
  channel type.
  
  The thresholds provided here must be at least as stringent as those in
  [`ica_reject`][config.ica_reject] if using ICA. In case of
  `'autoreject_global'`, thresholds for any channel that do not meet this
  requirement will be automatically replaced with those used in `ica_reject`.
  
  Note: Note
        The rejection is performed **after** SSP or ICA, if any of those methods
        is used. To reject epochs **before** fitting ICA, see the
        [`ica_reject`][config.ica_reject] setting.
  
  If `None` (default), do not apply automated rejection. If a dictionary,
  manually specify rejection thresholds (see examples).  If `'auto'`, use
  [`autoreject`](https://autoreject.github.io) to find suitable "global"
  rejection thresholds for each channel type, i.e. `autoreject` will generate
  a dictionary with (hopefully!) optimal thresholds for each channel type. Note
  that using `autoreject` can be a time-consuming process.
  
  Note: Note
        `autoreject` basically offers two modes of operation: "global" and
        "local". In "global" mode, it will try to estimate one rejection
        threshold **per channel type.** In "local" mode, it will generate
        thresholds **for each individual channel.** Currently, the BIDS Pipeline
        only supports the "global" mode.
  
  ???+ example "Example"
      ```python
      reject = {'grad': 4000e-13, 'mag': 4e-12, 'eog': 150e-6}
      reject = {'eeg': 100e-6, 'eog': 250e-6}
      reject = None  # no rejection based on PTP amplitude
      ```
  """
  
  reject_tmin: Optional[float] = None
  """
  Start of the time window used to reject epochs. If ``None``, the window will
  start with the first time point.
  ???+ example "Example"
      ```python
      reject_tmin = -0.1  # 100 ms before event onset.
      ```
  """
  
  reject_tmax: Optional[float] = None
  """
  End of the time window used to reject epochs. If ``None``, the window will end
  with the last time point.
  ???+ example "Example"
      ```python
      reject_tmax = 0.3  # 300 ms after event onset.
      ```
  """
  
  ###############################################################################
  # DECODING
  # --------
  
  decode: bool = True
  """
  Whether to perform decoding (MVPA) on the contrasts specified above as
  "contrasts". MVPA will be performed on the level of individual epochs.
  """
  
  decoding_metric: str = 'roc_auc'
  """
  The metric to use for estimating classification performance. It can be
  `'roc_auc'` or `'accuracy'` – or any other metric supported by `scikit-learn`.
  
  With ROC AUC, chance level is the same regardless of class balance, that is,
  you don't need to be worried about **exactly** balancing class sizes.
  """
  
  decoding_n_splits: int = 5
  """
  The number of folds (a.k.a. "splits") to use in the cross-validation scheme.
  """
  
  n_boot: int = 5000
  """
  The number of bootstrap resamples when estimating the standard error and
  confidence interval of the mean decoding score.
  """
  
  ###############################################################################
  # GROUP AVERAGE SENSORS
  # ---------------------
  
  interpolate_bads_grand_average: bool = True
  """
  Interpolate bad sensors in each dataset before calculating the grand
  average. This parameter is passed to the `mne.grand_average` function via
  the keyword argument `interpolate_bads`. It requires to have channel
  locations set.
  
  ???+ example "Example"
      ```python
      interpolate_bads_grand_average = True
      ```
  """
  
  ###############################################################################
  # TIME-FREQUENCY
  # --------------
  
  time_frequency_conditions: Iterable[str] = []
  """
  The conditions to compute time-frequency decomposition on.
  
  ???+ example "Example"
      ```python
      time_frequency_conditions = ['left', 'right']
      ```
  """
  
  time_frequency_freq_min: Optional[float] = 10
  """
  Minimum frequency for the time frequency analysis, in Hz.
  ???+ example "Example"
      ```python
      time_frequency_freq_min = 0.3  # 0.3 Hz
      ```
  """
  
  time_frequency_freq_max: Optional[float] = 40
  """
  Maximum frequency for the time frequency analysis, in Hz.
  ???+ example "Example"
      ```python
      time_frequency_freq_max = 22.3  # 22.3 Hz
      ```
  """
  
  ###############################################################################
  # SOURCE ESTIMATION PARAMETERS
  # ----------------------------
  #
  
  run_source_estimation: bool = True
  """
  Whether to run source estimation processing steps if not explicitly requested.
  """
  
  use_template_mri: bool = False
  """
  Whether to use FreeSurfer's `fsaverage` subject as MRI template. This may
  come in handy if you don't haver individual MR scans of your participants, as
  is often the case in EEG studies.
  """
  
  bem_mri_images: Literal['FLASH', 'T1', 'auto'] = 'auto'
  """
  Which types of MRI images to use when creating the BEM model.
  If ``'FLASH'``, use FLASH MRI images, and raise an exception if they cannot be
  found.
  
  ???+ info "Advice"
      It is recommended to use the FLASH images if available, as the quality
      of the extracted BEM surfaces will be higher.
  
  If ``'T1'``, create the BEM surfaces from the T1-weighted images using the
  ``watershed`` algorithm.
  
  If ``'auto'``, use FLASH images if available, and use the ``watershed``
  algorithm with the T1-weighted images otherwise.
  
  *[FLASH MRI]: Fast low angle shot magnetic resonance imaging
  """
  
  recreate_bem: bool = False
  """
  Whether to re-create the BEM surfaces, even if existing surfaces have been
  found. If ``False``, the BEM surfaces are only created if they do not exist
  already. ``True`` forces their recreation, overwriting existing BEM surfaces.
  """
  
  recreate_scalp_surface: bool = False
  """
  Whether to re-create the high-resolution scalp surface used for visualization
  of the coregistration in the report. If ``False``, the scalp surface is only
  created if it does not exist already. If ``True``, forces a re-computation.
  """
  
  mri_t1_path_generator: Optional[Callable] = None
  """
  To perform source-level analyses, the Pipeline needs to generate a
  transformation matrix that translates coordinates from MEG and EEG sensor
  space to MRI space, and vice versa. This process, called "coregistration",
  requires access to both, the electrophyisiological recordings as well as
  T1-weighted MRI images of the same participant. If both are stored within
  the same session, the Pipeline (or, more specifically, MNE-BIDS) can find the
  respective files automatically.
  
  However, in certain situations, this is not possible. Examples include:
  
  - MRI was conducted during a different session than the electrophysiological
    recording.
  - MRI was conducted in a single session, while electrophysiological recordings
    spanned across several sessions.
  - MRI and electrophysiological data are stored in separate BIDS datasets to
    allow easier storage and distribution in certain situations.
  
  To allow the Pipeline to find the correct MRI images and perform coregistration
  automatically, we provide a "hook" that allows you to provide in a custom
  function whose output tells the Pipeline where to find the T1-weighted image.
  
  The function is expected to accept a single parameter. The Pipeline will pass
  a `BIDSPath` with the following parameters set based on the currently processed
  electrophysiological data:
  
  - the subject ID, `BIDSPath.subject`
  - the experimental session, `BIDSPath.session`
  - the BIDS root, `BIDSPath.root`
  
  This `BIDSPath` can then be modified – or an entirely new `BIDSPath` can be
  generated – and returned by the function, pointing to the T1-weighted image.
  
  Note: Note
      The function accepts and returns a single `BIDSPath`.
  
  ???+ example "Example"
      The MRI session is different than the electrophysiological session:
      ```python
      def get_t1_from_meeg(bids_path):
          bids_path.session = 'MRI'
          return bids_path
  
  
      mri_t1_path_generator = get_t1_from_meeg
      ```
  
      The MRI recording is stored in a different BIDS dataset than the
      electrophysiological data:
      ```python
      def get_t1_from_meeg(bids_path):
          bids_path.root = '/data/mri'
          return bids_path
  
  
      mri_t1_path_generator = get_t1_from_meeg
      ```
  """
  
  spacing: Union[Literal['oct5', 'oct6', 'ico4', 'ico5', 'all'], int] = 'oct6'
  """
  The spacing to use. Can be ``'ico#'`` for a recursively subdivided
  icosahedron, ``'oct#'`` for a recursively subdivided octahedron,
  ``'all'`` for all points, or an integer to use approximate
  distance-based spacing (in mm). See (the respective MNE-Python documentation)
  [https://mne.tools/dev/overview/cookbook.html#setting-up-the-source-space]
  for more info.
  """
  
  mindist: float = 5
  """
  Exclude points closer than this distance (mm) to the bounding surface.
  """
  
  # loose: Union[float, Literal['auto']] = 0.2
  # # ``loose`` : float in [0, 1] | 'auto'
  # """
  # Value that weights the source variances of the dipole components
  # that are parallel (tangential) to the cortical surface. If ``0``, then the
  # inverse solution is computed with **fixed orientation.**
  # If ``1``, it corresponds to **free orientation.**
  # The default value, ``'auto'``, is set to ``0.2`` for surface-oriented source
  # spaces, and to ``1.0`` for volumetric, discrete, or mixed source spaces,
  # unless ``fixed is True`` in which case the value 0. is used.
  # """
  
  # depth: Optional[Union[float, dict]] = 0.8
  # """
  # If float (default 0.8), it acts as the depth weighting exponent (``exp``)
  # to use (must be between 0 and 1). None is equivalent to 0, meaning no
  # depth weighting is performed. Can also be a `dict` containing additional
  # keyword arguments to pass to :func:`mne.forward.compute_depth_prior`
  # (see docstring for details and defaults).
  # """
  
  inverse_method: Literal['MNE', 'dSPM', 'sLORETA', 'eLORETA'] = 'dSPM'
  """
  Use minimum norm, dSPM (default), sLORETA, or eLORETA to calculate the inverse
  solution.
  """
  
  noise_cov: Union[Tuple[Optional[float], Optional[float]],
                   Literal['emptyroom']] = (None, 0)
  """
  Specify how to estimate the noise covariance matrix, which is used in
  inverse modeling.
  
  If a tuple, it takes the form ``(tmin, tmax)`` with the time specified in
  seconds. If the first value of the tuple is ``None``, the considered
  period starts at the beginning of the epoch. If the second value of the
  tuple is ``None``, the considered period ends at the end of the epoch.
  The default, ``(None, 0)``, includes the entire period before the event,
  which is typically the pre-stimulus period.
  
  If ``emptyroom``, the noise covariance matrix will be estimated from an
  empty-room MEG recording. The empty-room recording will be automatically
  selected based on recording date and time.
  
  Please note that when processing data that contains EEG channels, the noise
  covariance can ONLY be estimated from the pre-stimulus period.
  
  ???+ example "Example"
      Use the period from start of the epoch until 100 ms before the experimental
      event:
      ```python
      noise_cov = (None, -0.1)
      ```
  
      Use the time period from the experimental event until the end of the epoch:
      ```python
      noise_cov = (0, None)
      ```
  
      Use an empty-room recording:
      ```python
      noise_cov = 'emptyroom'
      ```
  """
  
  ###############################################################################
  # ADVANCED
  # --------
  
  l_trans_bandwidth: Union[float, Literal['auto']] = 'auto'
  """
  Specifies the transition bandwidth of the
  highpass filter. By default it's `'auto'` and uses default MNE
  parameters.
  """
  
  h_trans_bandwidth: Union[float, Literal['auto']] = 'auto'
  """
  Specifies the transition bandwidth of the
  lowpass filter. By default it's `'auto'` and uses default MNE
  parameters.
  """
  
  N_JOBS: int = 1
  """
  Specifies how many subjects you want to process in parallel.
  """
  
  random_state: Optional[int] = 42
  """
  You can specify the seed of the random number generator (RNG).
  This setting is passed to the ICA algorithm and to the decoding function,
  ensuring reproducible results. Set to ``None`` to avoid setting the RNG
  to a defined state.
  """
  
  shortest_event: int = 1
  """
  Minimum number of samples an event must last. If the
  duration is less than this, an exception will be raised.
  """
  
  log_level: Literal['info', 'error'] = 'info'
  """
  Set the pipeline logging verbosity.
  """
  
  mne_log_level: Literal['info', 'error'] = 'error'
  """
  Set the MNE-Python logging verbosity.
  """
  
  OnErrorT = Literal['continue', 'abort', 'debug']
  on_error: OnErrorT = 'abort'
  """
  Whether to abort processing as soon as an error occurs, continue with all other
  processing steps for as long as possible, or drop you into a debugger in case
  of an error.
  """
  
  ###############################################################################
  #                                                                             #
  #                      CUSTOM CONFIGURATION ENDS HERE                         #
  #                                                                             #
  ###############################################################################
  
  ###############################################################################
  # Name, version, and hosting location of the pipeline
  # ---------------------------------------------------
  
  PIPELINE_NAME = 'mne-bids-pipeline'
  VERSION = '0.1.dev0'
  CODE_URL = 'https://github.com/mne-tools/mne-bids-pipeline'
  
  os.environ['MNE_BIDS_STUDY_SCRIPT_PATH'] = str(__file__)
  
  ###############################################################################
  # Logger
  # ------
  
  logger = logging.getLogger('mne-bids-pipeline')
  
  log_fmt = '[%(asctime)s] %(step)s%(subject)s%(session)s%(run)s%(message)s'
  log_date_fmt = coloredlogs.DEFAULT_DATE_FORMAT = '%H:%M:%S'
  log_level_styles = {
      'warning': {
          'color': 202,
          'bold': True
      },
      'error': {
          'background': 'red',
          'bold': True
      },
      'critical': {
          'background': 'red',
          'bold': True
      }
  }
  log_field_styles = {
      'asctime': {
          'color': 'green'
      },
      'step': {
          'color': 'white',
          'bold': True
      },
      'subject': {
          'color': 'blue',
          'bright': True,
          'bold': True
      },
      'session': {
          'bold': True
      },
      'run': {
          'bold': True
      }
  }
  
  
  class LogFilter(logging.Filter):
      def filter(self, record):
          if not hasattr(record, 'step'):
              record.step = ''
          if not hasattr(record, 'subject'):
              record.subject = ''
          if not hasattr(record, 'session'):
              record.session = ''
          if not hasattr(record, 'run'):
              record.run = ''
  
          return True
  
  
  logger.addFilter(LogFilter())
  
  coloredlogs.install(
      level=log_level, logger=logger, fmt=log_fmt,
      date_fmt=log_date_fmt, level_styles=log_level_styles,
      field_styles=log_field_styles
  )
  
  
  mne.set_log_level(verbose=mne_log_level.upper())
  
  
  class LogKwargsT(TypedDict):
      msg: str
      extra: Dict[str, str]
  
  
  def gen_log_kwargs(
      message: str,
      subject: Optional[Union[str, int]] = None,
      session: Optional[Union[str, int]] = None,
      run: Optional[Union[str, int]] = None
  ) -> LogKwargsT:
      if subject is not None:
          subject = f' sub-{subject}'
      if session is not None:
          session = f' ses-{session}'
      if run is not None:
          run = f' run-{run}'
  
      message = f' {message}'
  
      script_path = pathlib.Path(os.environ['MNE_BIDS_STUDY_SCRIPT_PATH'])
      step_name = f'{script_path.parent.name}/{script_path.stem}'
  
      extra = {
          'step': step_name
      }
      if subject:
          extra['subject'] = subject
      if session:
          extra['session'] = session
      if run:
          extra['run'] = run
  
      kwargs: LogKwargsT = {
          'msg': message,
          'extra': extra
      }
      return kwargs
  
  
  ###############################################################################
  # Retrieve custom configuration options
  # -------------------------------------
  #
  # For testing a specific dataset, create a Python file with a name of your
  # liking (e.g., ``mydataset-template-config.py``), and set an environment
  # variable ``MNE_BIDS_STUDY_CONFIG`` to that file.
  #
  # Example
  # ~~~~~~~
  # ``export MNE_BIDS_STUDY_CONFIG=/data/mystudy/mydataset-template-config.py``
  
  if "MNE_BIDS_STUDY_CONFIG" in os.environ:
      cfg_path = pathlib.Path(os.environ['MNE_BIDS_STUDY_CONFIG'])
  
      if cfg_path.exists():
          msg = f'Using custom configuration: {cfg_path}'
          logger.info(**gen_log_kwargs(message=msg))
      else:
          msg = ('The custom configuration file specified in the '
                 'MNE_BIDS_STUDY_CONFIG environment variable could not be '
                 'found: {cfg_path}'.format(cfg_path=cfg_path))
          raise ValueError(msg)
  
      # Import configuration from an arbitrary path without having to fiddle
      # with `sys.path`.
      spec = importlib.util.spec_from_file_location(name='custom_config',
                                                    location=cfg_path)
      custom_cfg = importlib.util.module_from_spec(spec)
      spec.loader.exec_module(custom_cfg)
      del spec, cfg_path
  
      new = None
      for val in dir(custom_cfg):
          if not val.startswith('__'):
              exec("new = custom_cfg.%s" % val)
              logger.debug('Overwriting: %s -> %s' % (val, new))
              exec("%s = custom_cfg.%s" % (val, val))
  
  
  ###############################################################################
  # CHECKS
  # ------
  
  if (use_maxwell_filter and
          len(set(ch_types).intersection(('meg', 'grad', 'mag'))) == 0):
      raise ValueError('Cannot use maxwell filter without MEG channels.')
  
  if (spatial_filter == 'ica' and
          ica_algorithm not in ('picard', 'fastica', 'extended_infomax')):
      msg = (f"Invalid ICA algorithm requested. Valid values for ica_algorithm "
             f"are: 'picard', 'fastica', and 'extended_infomax', but received "
             f"{ica_algorithm}.")
      raise ValueError(msg)
  
  if (spatial_filter == 'ica' and
          ica_l_freq is None and
          l_freq is not None and l_freq < 1):
      msg = (f'You requested to high-pass filter your data with l_freq={l_freq} '
             f'Hz and to perform ICA without performing any additional '
             f'filtering first by setting ica_l_freq=None. However, ICA will '
             f'not work reliably unless slow drifts have been removed from the '
             f'data. Please either increase l_freq to 1 Hz or above, or enable '
             f'additional filtering for ICA by setting ica_l_freq to 1 Hz or '
             f'higher.')
      raise ValueError(msg)
  
  if spatial_filter == 'ica' and ica_l_freq < 1:
      msg = (f'You requested to high-pass filter the data before ICA with '
             f'ica_l_freq={ica_l_freq} Hz. Please increase this setting to '
             f'1 Hz or above to ensure reliable ICA function.')
      raise ValueError(msg)
  
  if (spatial_filter == 'ica' and
          ica_l_freq is not None and
          l_freq is not None and
          ica_l_freq < l_freq):
      msg = (f'You requested a lower high-pass filter cutoff frequency for ICA '
             f'than for your raw data: ica_l_freq = {ica_l_freq} < '
             f'l_freq = {l_freq}. Adjust the cutoffs such that ica_l_freq >= '
             f'l_freq, or set ica_l_freq to None if you do not wish to apply '
             f'an additional high-pass filter before running ICA.')
      raise ValueError(msg)
  
  if not ch_types:
      msg = 'Please specify ch_types in your configuration.'
      raise ValueError(msg)
  
  if ch_types == ['eeg']:
      pass
  elif 'eeg' in ch_types and len(ch_types) > 1:  # EEG + some other channel types
      msg = ('EEG data can only be analyzed separately from other channel '
             'types. Please adjust `ch_types` in your configuration.')
      raise ValueError(msg)
  elif any([ch_type not in ('meg', 'mag', 'grad') for ch_type in ch_types]):
      msg = ('Invalid channel type passed. Please adjust `ch_types` in your '
             'configuration.')
      raise ValueError(msg)
  
  if 'eeg' in ch_types:
      if spatial_filter == 'ssp':
          msg = ("You requested SSP for EEG data via spatial_filter='ssp'. "
                 "However, this is not presently supported. Please use ICA "
                 "instead by setting spatial_filter='ica'.")
          raise ValueError(msg)
  
  
  if on_error not in ('continue', 'abort', 'debug'):
      msg = (f"on_error must be one of 'continue', 'debug' or 'abort', "
             f"but received: {on_error}.")
      logger.info(**gen_log_kwargs(message=msg))
  
  if isinstance(noise_cov, str) and noise_cov != 'emptyroom':
      msg = (f"noise_cov must be a tuple or 'emptyroom', but received "
             f"{noise_cov}")
      raise ValueError(msg)
  
  if noise_cov == 'emptyroom' and 'eeg' in ch_types:
      msg = ('You requested to process data that contains EEG channels. In this '
             'case, noise covariance can only be estimated from the '
             'experimental data, e.g., the pre-stimulus period. Please set '
             'noise_cov to (tmin, tmax)')
      raise ValueError(msg)
  
  if noise_cov == 'emptyroom' and not process_er:
      msg = ('You requested noise covariance estimation from empty-room '
             'recordings by setting noise_cov = "emptyroom", but you did not '
             'enable empty-room data processing. Please set process_er = True')
      raise ValueError(msg)
  
  if bem_mri_images not in ('FLASH', 'T1', 'auto'):
      msg = (f'Unknown bem_mri_images: {bem_mri_images}. Valid values '
             f'are: "FLASH", "T1", and "auto".')
      raise ValueError(msg)
  
  
  def check_baseline(
      *,
      baseline: Optional[Tuple[Optional[float], Optional[float]]],
      epochs_tmin: float,
      epochs_tmax: float
  ) -> None:
      """Raises error if baseline not compatible with [epochs_tmin, epochs_tmax].
  
      Parameters
      ----------
      baseline
          Tuple indicating the beginning and end of the baseline interval.
      epochs_tmin
          Beginning of Epochs.
      epochs_tmax
          End of Epochs.
  
      Raises
      ------
      ValueError
          if baseline not contained in [epochs_tmin, epochs_tmax].
          if baseline is not a correct time-interval.
      """
      if baseline is None:
          return
  
      if ((baseline[0] is not None and baseline[0] < epochs_tmin) or
              (baseline[1] is not None and baseline[1] > epochs_tmax)):
          msg = (f'baseline {baseline} outside of epochs interval '
                 f'{[epochs_tmin, epochs_tmax]}.')
  
          raise ValueError(msg)
  
      if ((baseline[0] is not None) and
              (baseline[1] is not None) and
              (baseline[0] >= baseline[1])):
          msg = (f'The end of the baseline period must occur after its start, '
                 f'but you set baseline={baseline}')
          raise ValueError(msg)
  
  
  check_baseline(baseline=baseline, epochs_tmin=epochs_tmin,
                 epochs_tmax=epochs_tmax)
  
  
  # check PTP rejection thresholds
  if (spatial_filter == 'ica' and
          ica_reject is not None and
          reject is not None and
          reject != 'autoreject_global'):
      for ch_type in reject:
          if (ch_type in ica_reject and
                  reject[ch_type] > ica_reject[ch_type]):
              raise ValueError(
                  f'Rejection threshold in '
                  f'reject["{ch_type}"] ({reject[ch_type]}) must be at least as '
                  f'stringent as that in '
                  f'ica_reject["{ch_type}"] ({ica_reject[ch_type]})'
              )
  
  
  ###############################################################################
  # Helper functions
  # ----------------
  
  
  @functools.lru_cache(maxsize=None)
  def _get_entity_vals_cached(*args, **kwargs):
      return mne_bids.get_entity_vals(*args, **kwargs)
  
  
  def get_bids_root() -> pathlib.Path:
      # BIDS_ROOT environment variable takes precedence over any configuration file
      # values.
      root = os.getenv('BIDS_ROOT')
      if root is not None:
          return (pathlib.Path(root)
                  .expanduser()
                  .resolve(strict=True))
  
      # If we don't have a bids_root until now, raise an exception as we cannot
      # proceed.
      if not bids_root:
          msg = ('You need to specify `bids_root` in your configuration, or '
                 'define an environment variable `BIDS_ROOT` pointing to the '
                 'root folder of your BIDS dataset')
          raise ValueError(msg)
  
      return (pathlib.Path(bids_root)
              .expanduser()
              .resolve(strict=True))
  
  
  def get_datatype() -> Literal['meg', 'eeg']:
      # Content of ch_types should be sanitized already, so we don't need any
      # extra sanity checks here.
      if data_type is not None:
          return data_type
      elif data_type is None and ch_types == ['eeg']:
          return 'eeg'
      elif data_type is None and any([t in ['meg', 'mag', 'grad']
                                      for t in ch_types]):
          return 'meg'
      else:
          raise RuntimeError("This probably shouldn't happen. Please contact "
                             "the MNE-BIDS-pipeline developers. Thank you.")
  
  
  _all_datatypes = mne_bids.get_datatypes(root=get_bids_root())
  _ignore_datatypes = set(_all_datatypes) - set([get_datatype()])
  
  _valid_tasks = _get_entity_vals_cached(
      root=get_bids_root(),
      entity_key='task',
      ignore_datatypes=tuple(_ignore_datatypes)
  )
  
  _valid_subjects = _get_entity_vals_cached(
      root=get_bids_root(),
      entity_key='subject',
      ignore_datatypes=tuple(_ignore_datatypes)
  )
  
  _all_sessions = _get_entity_vals_cached(
      root=get_bids_root(),
      entity_key='session',
      ignore_datatypes=tuple(_ignore_datatypes)
  )
  
  
  def get_deriv_root() -> pathlib.Path:
      if deriv_root is None:
          return get_bids_root() / 'derivatives' / PIPELINE_NAME
      else:
          return (pathlib.Path(deriv_root)
                  .expanduser()
                  .resolve())
  
  
  def get_subjects() -> List[str]:
      global subjects
  
      env = os.environ
      if env.get('MNE_BIDS_STUDY_SUBJECT'):
          env_subject = env['MNE_BIDS_STUDY_SUBJECT']
          if env_subject not in _valid_subjects:
              raise ValueError(
                  f'Invalid subject. It can be {_valid_subjects} but '
                  f'got {env_subject}')
          s = [env_subject]
      elif subjects == 'all':
          s = _valid_subjects
      else:
          s = subjects
  
      subjects = set(s) - set(exclude_subjects)
      # Drop empty-room subject.
      subjects = subjects - set(['emptyroom'])
  
      return sorted(subjects)
  
  
  def get_sessions():
      sessions_ = copy.deepcopy(sessions)  # Avoid clash with global variable.
  
      env = os.environ
      if env.get('MNE_BIDS_STUDY_SESSION'):
          sessions_ = env['MNE_BIDS_STUDY_SESSION']
      elif sessions_ == 'all':
          sessions_ = _all_sessions
  
      if not sessions_:
          return [None]
      else:
          return sessions_
  
  
  def get_runs_all_subjects() -> dict:
      """Gives the mapping between subjects and their runs.
  
      Returns
      -------
      a dict of runs present in the bids_path
      for each subject asked in the configuration file
      (and not for each subject present in the bids_path).
      """
      # We cannot use get_subjects() because if there is just one subject
      subj_runs = dict()
      for subject in get_subjects():
          # Only traverse through the current subject's directory
          valid_runs_subj = _get_entity_vals_cached(
              get_bids_root() / f'sub-{subject}', entity_key='run',
              ignore_datatypes=tuple(_ignore_datatypes)
          )
          if exclude_runs and subject in exclude_runs:
              valid_runs_subj = [r for r in valid_runs_subj
                                 if r not in exclude_runs[subject]]
          subj_runs[subject] = valid_runs_subj
  
      return subj_runs
  
  
  def get_intersect_run() -> List[str]:
      """Returns the intersection of all the runs of all subjects."""
      subj_runs = get_runs_all_subjects()
      return list(set.intersection(*map(set, subj_runs.values())))
  
  
  def get_runs(
      *,
      subject: str,
      verbose: bool = False
  ) -> Union[List[str], List[None]]:
      """Returns a list of runs in the BIDS input data.
  
      Parameters
      ----------
      subject
          Returns a list of the runs of this subject.
      verbose
          Notify if different subjects do not share the same runs.
  
      Returns
      -------
      The list of runs of the subject. If no BIDS `run` entity could be found,
      returns `[None]`.
      """
      if subject == 'average':  # Used when creating the report
          return [None]
  
      runs_ = copy.deepcopy(runs)  # Avoid clash with global variable.
  
      subj_runs = get_runs_all_subjects()
      valid_runs = subj_runs[subject]
  
      if len(get_subjects()) > 1:
          # Notify if different subjects do not share the same runs
  
          same_runs = True
          for runs_sub_i in subj_runs.values():
              if set(runs_sub_i) != set(list(subj_runs.values())[0]):
                  same_runs = False
  
          if not same_runs and verbose:
              msg = ('Extracted all the runs. '
                     'Beware, not all subjects share the same '
                     'set of runs.')
              logger.info(**gen_log_kwargs(message=msg))
  
      env_run = os.environ.get('MNE_BIDS_STUDY_RUN')
      if env_run and env_run not in valid_runs:
          raise ValueError(
              f'Invalid run. It can be {valid_runs} but '
              f'got {env_run}')
      elif env_run:
          runs_ = [env_run]
      elif runs_ == 'all':
          runs_ = valid_runs
  
      if not runs_:
          return [None]
      else:
          inclusion = set(runs_).issubset(set(valid_runs))
          if not inclusion:
              raise ValueError(
                  f'Invalid run. It can be a subset of {valid_runs} but '
                  f'got {runs_}')
          return runs_
  
  
  def get_mf_reference_run() -> str:
      # Retrieve to run identifier (number, name) of the reference run
      if mf_reference_run is None:
          # Use the first run
          if inter_runs:
              return inter_runs[0]
          else:
              ValueError("The intersection of runs by subjects is empty. "
                         "Check the list of runs: "
                         f"{get_runs_all_subjects()}")
      else:
          return mf_reference_run
  
  
  def get_task() -> Optional[str]:
      global task
  
      env = os.environ
      if env.get('MNE_BIDS_STUDY_TASK'):
          task = env['MNE_BIDS_STUDY_TASK']
          if task not in _valid_tasks and 'MKDOCS' not in os.environ:
              raise ValueError(f'Invalid task. It can be: '
                               f'{", ".join(_valid_tasks)} but got: {task}')
  
      if not task:
          if not _valid_tasks:
              return None
          else:
              return _valid_tasks[0]
      else:
          return task
  
  
  def get_n_jobs() -> int:
      env = os.environ
      if env.get('MNE_BIDS_STUDY_NJOBS'):
          n_jobs = int(env['MNE_BIDS_STUDY_NJOBS'])
      else:
          n_jobs = N_JOBS
  
      return n_jobs
  
  
  def _get_reject(
      *,
      subject: Optional[str] = None,
      session: Optional[str] = None,
      reject: Optional[Union[Dict[str, float], Literal['autoreject_global']]],
      ch_types: Iterable[Literal['meg', 'mag', 'grad', 'eeg']],
      epochs: Optional[mne.BaseEpochs] = None
  ) -> Dict[str, float]:
      if reject is None:
          return dict()
  
      if reject == 'autoreject_global':
          # Automated threshold calculation requested
          import autoreject
  
          ch_types_autoreject = list(ch_types)
          if 'meg' in ch_types_autoreject:
              ch_types_autoreject.remove('meg')
              if 'mag' in epochs:
                  ch_types_autoreject.append('mag')
              if 'grad' in epochs:
                  ch_types_autoreject.append('grad')
  
          msg = 'Generating rejection thresholds using autoreject …'
          logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                       session=session))
          reject = autoreject.get_rejection_threshold(
              epochs=epochs, ch_types=ch_types_autoreject, decim=decim,
              verbose=False
          )
          return reject
  
      # Only keep thresholds for channel types of interest
      reject = reject.copy()
      if ch_types == ['eeg']:
          ch_types_to_remove = ('mag', 'grad')
      else:
          ch_types_to_remove = ('eeg',)
  
      for ch_type in ch_types_to_remove:
          try:
              del reject[ch_type]
          except KeyError:
              pass
  
      return reject
  
  
  def get_reject(
      *,
      epochs: Optional[mne.BaseEpochs] = None
  ) -> Dict[str, float]:
      return _get_reject(reject=reject, ch_types=ch_types, epochs=epochs)
  
  
  def get_ssp_reject(
      *,
      ssp_type: Literal['ecg', 'eog'],
      epochs: mne.BaseEpochs
  ) -> Dict[str, float]:
      if ssp_type == 'ecg':
          reject = ssp_reject_ecg
      elif ssp_type == 'eog':
          reject = ssp_reject_eog
      else:
          raise ValueError("Only 'eog' and 'ecg' are supported.")
      return _get_reject(reject=reject, ch_types=ch_types, epochs=epochs)
  
  
  def get_ica_reject() -> Dict[str, float]:
      return _get_reject(reject=ica_reject, ch_types=ch_types)
  
  
  def get_fs_subjects_dir():
      if not subjects_dir and deriv_root is not None:
          # We do this check here (and not in our regular checks section) to
          # avoid an error message when a user doesn't intend to run the source
          # analysis steps anyway.
          raise ValueError(
              'When specifying a "deriv_root", you must also supply a '
              '"subjects_dir".'
          )
  
      if not subjects_dir:
          return get_bids_root() / 'derivatives' / 'freesurfer' / 'subjects'
      else:
          return (pathlib.Path(subjects_dir)
                  .expanduser()
                  .resolve())
  
  
  def failsafe_run(
      on_error: OnErrorT,
      script_path: PathLike,
  ):
      def failsafe_run_decorator(func):
          @functools.wraps(func)  # Preserve "identity" of original function
          def wrapper(*args, **kwargs):
              os.environ['MNE_BIDS_STUDY_SCRIPT_PATH'] = str(script_path)
              kwargs_copy = copy.deepcopy(kwargs)
              t0 = time.time()
              if "cfg" in kwargs_copy:
                  kwargs_copy["cfg"] = json_tricks.dumps(
                      kwargs_copy["cfg"], sort_keys=False, indent=4
                  )
              log_info = pd.Series(kwargs_copy)
              log_info['time'] = np.nan  # put time before success & err columns
              try:
                  assert len(args) == 0  # make sure params are only kwargs
                  out = func(*args, **kwargs)
                  assert out is None  # make sure the function return None
                  log_info['success'] = True
                  log_info['error_message'] = ''
              except Exception as e:
                  message = 'A critical error occurred.'
                  log_info['success'] = False
                  log_info['error_message'] = str(e)
  
                  if on_error == 'abort':
                      logger.critical(**gen_log_kwargs(message=message))
                      raise(e)
                  elif on_error == 'debug':
                      logger.critical(**gen_log_kwargs(message=message))
                      extype, value, tb = sys.exc_info()
                      traceback.print_exc()
                      pdb.post_mortem(tb)
                  else:
                      message = f'{message} The error message was:\n{str(e)}'
                      logger.critical(**gen_log_kwargs(message=message))
              log_info['time'] = round(time.time() - t0, ndigits=1)
              return log_info
          return wrapper
      return failsafe_run_decorator
  
  
  def plot_auto_scores(auto_scores):
      # Plot scores of automated bad channel detection.
      import matplotlib.pyplot as plt
      import seaborn as sns
      import pandas as pd
  
      if ch_types == ['meg']:
          ch_types_ = ['grad', 'mag']
      else:
          ch_types_ = ch_types
  
      figs = []
      for ch_type in ch_types_:
          # Only select the data for mag or grad channels.
          ch_subset = auto_scores['ch_types'] == ch_type
          ch_names = auto_scores['ch_names'][ch_subset]
          scores = auto_scores['scores_noisy'][ch_subset]
          limits = auto_scores['limits_noisy'][ch_subset]
          bins = auto_scores['bins']  # The the windows that were evaluated.
  
          # We will label each segment by its start and stop time, with up to 3
          # digits before and 3 digits after the decimal place (1 ms precision).
          bin_labels = [f'{start:3.3f} – {stop:3.3f}'
                        for start, stop in bins]
  
          # We store the data in a Pandas DataFrame. The seaborn heatmap function
          # we will call below will then be able to automatically assign the
          # correct labels to all axes.
          data_to_plot = pd.DataFrame(data=scores,
                                      columns=pd.Index(bin_labels,
                                                       name='Time (s)'),
                                      index=pd.Index(ch_names, name='Channel'))
  
          # First, plot the "raw" scores.
          fig, ax = plt.subplots(1, 2, figsize=(12, 8))
          fig.suptitle(f'Automated noisy channel detection: {ch_type}',
                       fontsize=16, fontweight='bold')
          sns.heatmap(data=data_to_plot, cmap='Reds',
                      cbar_kws=dict(label='Score'), ax=ax[0])
          [ax[0].axvline(x, ls='dashed', lw=0.25, dashes=(25, 15), color='gray')
              for x in range(1, len(bins))]
          ax[0].set_title('All Scores', fontweight='bold')
  
          # Now, adjust the color range to highlight segments that exceeded the
          # limit.
          sns.heatmap(data=data_to_plot,
                      vmin=np.nanmin(limits),  # input data may contain NaNs
                      cmap='Reds', cbar_kws=dict(label='Score'), ax=ax[1])
          [ax[1].axvline(x, ls='dashed', lw=0.25, dashes=(25, 15), color='gray')
              for x in range(1, len(bins))]
          ax[1].set_title('Scores > Limit', fontweight='bold')
  
          # The figure title should not overlap with the subplots.
          fig.tight_layout(rect=[0, 0.03, 1, 0.95])
          figs.append(fig)
  
      return figs
  
  
  def get_channels_to_analyze(info) -> List[str]:
      # Return names of the channels of the channel types we wish to analyze.
      # We also include channels marked as "bad" here.
      # `exclude=[]`: keep "bad" channels, too.
      if get_datatype() == 'meg' and ('mag' in ch_types or 'grad' in ch_types
                                      or 'meg' in ch_types):
          pick_idx = mne.pick_types(info, eog=True, ecg=True, exclude=[])
  
          if 'mag' in ch_types:
              pick_idx = np.concatenate(
                  [pick_idx, mne.pick_types(info, meg='mag', exclude=[])])
          if 'grad' in ch_types:
              pick_idx = np.concatenate(
                  [pick_idx, mne.pick_types(info, meg='grad', exclude=[])])
          if 'meg' in ch_types:
              pick_idx = mne.pick_types(info, meg=True, eog=True, ecg=True,
                                        exclude=[])
      elif ch_types == ['eeg']:
          pick_idx = mne.pick_types(info, meg=False, eeg=True, eog=True,
                                    ecg=True, exclude=[])
      else:
          raise RuntimeError('Something unexpected happened. Please contact '
                             'the mne-bids-pipeline developers. Thank you.')
  
      ch_names = [info['ch_names'][i] for i in pick_idx]
      return ch_names
  
  
  def get_fs_subject(subject) -> str:
      subjects_dir = get_fs_subjects_dir()
  
      if use_template_mri:
          return 'fsaverage'
  
      if (pathlib.Path(subjects_dir) / subject).exists():
          return subject
      else:
          return f'sub-{subject}'
  
  
  def sanitize_cond_name(cond: str) -> str:
      cond = (cond
              .replace(os.path.sep, '')
              .replace('_', '')
              .replace('-', ''))
      return cond
  
  
  def get_mf_cal_fname(
      subject: str,
      session: str
  ) -> pathlib.Path:
      if mf_cal_fname is None:
          mf_cal_fpath = (BIDSPath(subject=subject,
                                   session=session,
                                   suffix='meg',
                                   datatype='meg',
                                   root=get_bids_root())
                          .meg_calibration_fpath)
          if mf_cal_fpath is None:
              raise ValueError('Could not find Maxwell Filter Calibration '
                               'file.')
      else:
          mf_cal_fpath = pathlib.Path(mf_cal_fname).expanduser().absolute()
          if not mf_cal_fpath.exists():
              raise ValueError(f'Could not find Maxwell Filter Calibration '
                               f'file at {str(mf_cal_fpath)}.')
  
      return mf_cal_fpath
  
  
  def get_mf_ctc_fname(
      subject: str,
      session: str
  ) -> pathlib.Path:
      if mf_ctc_fname is None:
          mf_ctc_fpath = (BIDSPath(subject=subject,
                                   session=session,
                                   suffix='meg',
                                   datatype='meg',
                                   root=get_bids_root())
                          .meg_crosstalk_fpath)
          if mf_ctc_fpath is None:
              raise ValueError('Could not find Maxwell Filter cross-talk '
                               'file.')
      else:
          mf_ctc_fpath = pathlib.Path(mf_ctc_fname).expanduser().absolute()
          if not mf_ctc_fpath.exists():
              raise ValueError(f'Could not find Maxwell Filter cross-talk '
                               f'file at {str(mf_ctc_fpath)}.')
  
      return mf_ctc_fpath
  
  
  def make_epochs(
      *,
      raw: mne.io.BaseRaw,
      event_id: Optional[Dict[str, int]] = None,
      tmin: float,
      tmax: float,
      metadata_tmin: Optional[float] = None,
      metadata_tmax: Optional[float] = None,
      metadata_keep_first: Optional[Iterable[str]] = None,
      metadata_keep_last: Optional[Iterable[str]] = None,
      event_repeated: Literal['error', 'drop', 'merge'],
      decim: int
  ) -> mne.Epochs:
      """Generate Epochs from raw data.
  
      No EEG reference will be set and no projectors will be applied. No
      rejection thresholds will be applied. No baseline-correction will be
      performed.
      """
      if get_task().lower() == 'rest':
          stop = raw.times[-1] - rest_epochs_duration
          assert epochs_tmin == 0., "epochs_tmin must be 0 for rest"
          assert rest_epochs_overlap is not None, \
              "epochs_overlap cannot be None for rest"
          events = mne.make_fixed_length_events(
              raw, id=3000, start=0,
              duration=rest_epochs_duration,
              overlap=rest_epochs_overlap,
              stop=stop)
          event_id = dict(rest=3000)
      else:  # Events for task runs
          events, event_id_from_annotations = mne.events_from_annotations(raw)
  
      if event_id is None:
          event_id = event_id_from_annotations
  
      # Construct metadata from the epochs
      if metadata_tmin is None:
          metadata_tmin = tmin
  
      if metadata_tmax is None:
          metadata_tmax = tmax
  
      metadata, _, _ = mne.epochs.make_metadata(
          events=events, event_id=event_id,
          tmin=metadata_tmin, tmax=metadata_tmax,
          keep_first=metadata_keep_first,
          keep_last=metadata_keep_last,
          sfreq=raw.info['sfreq'])
  
      # Epoch the data
      # Do not reject based on peak-to-peak or flatness thresholds at this stage
      epochs = mne.Epochs(raw, events=events, event_id=event_id,
                          tmin=tmin, tmax=tmax,
                          proj=False, baseline=None,
                          preload=False, decim=decim,
                          metadata=metadata,
                          event_repeated=event_repeated,
                          reject=None)
  
      return epochs
  
  
  def annotations_to_events(
      *,
      raw_paths: List[PathLike]
  ) -> Dict[str, int]:
      """Generate a unique event name -> event code mapping.
  
      The mapping can that can be used across all passed raws.
      """
      event_names: List[str] = []
      for raw_fname in raw_paths:
          raw = mne.io.read_raw_fif(raw_fname)
          _, event_id = mne.events_from_annotations(raw=raw)
          for event_name in event_id.keys():
              if event_name not in event_names:
                  event_names.append(event_name)
  
      event_names = sorted(event_names)
      event_name_to_code_map = {
          name: code
          for code, name in enumerate(event_names, start=1)
      }
  
      return event_name_to_code_map
  
  
  def _rename_events_func(cfg, raw, subject, session) -> None:
      """Rename events (actually, annotations descriptions) in ``raw``.
  
      Modifies ``raw`` in-place.
      """
      if not cfg.rename_events:
          return
  
      # Check if the user requested to rename events that don't exist.
      # We don't want this to go unnoticed.
      event_names_set = set(raw.annotations.description)
      rename_events_set = set(cfg.rename_events.keys())
      events_not_in_raw = rename_events_set - event_names_set
      if events_not_in_raw:
          msg = (f'You requested to rename the following events, but '
                 f'they are not present in the BIDS input data:\n'
                 f'{", ".join(sorted(list(events_not_in_raw)))}')
          if on_rename_missing_events == 'warn':
              logger.warning(msg)
          else:
              raise ValueError(msg)
  
      # Do the actual event renaming.
      msg = 'Renaming events …'
      logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                   session=session))
      descriptions = list(raw.annotations.description)
      for old_event_name, new_event_name in cfg.rename_events.items():
          msg = f'… {old_event_name} -> {new_event_name}'
          logger.info(**gen_log_kwargs(message=msg,
                                       subject=subject, session=session))
          for idx, description in enumerate(descriptions.copy()):
              if description == old_event_name:
                  descriptions[idx] = new_event_name
  
      descriptions = np.asarray(descriptions, dtype=str)
      raw.annotations.description = descriptions
  
  
  def _find_bad_channels(cfg, raw, subject, session, task, run) -> None:
      """Find and mark bad MEG channels.
  
      Modifies ``raw`` in-place.
      """
      if not (cfg.find_flat_channels_meg or cfg.find_noisy_channels_meg):
          return
  
      if (cfg.find_flat_channels_meg and
              not cfg.find_noisy_channels_meg):
          msg = 'Finding flat channels.'
      elif (cfg.find_noisy_channels_meg and
              not cfg.find_flat_channels_meg):
          msg = 'Finding noisy channels using Maxwell filtering.'
      else:
          msg = ('Finding flat channels, and noisy channels using '
                 'Maxwell filtering.')
  
      logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                   session=session))
  
      bids_path = BIDSPath(subject=subject,
                           session=session,
                           task=task,
                           run=run,
                           acquisition=acq,
                           processing=proc,  # XXX : what is proc?
                           recording=cfg.rec,
                           space=cfg.space,
                           suffix=cfg.datatype,
                           datatype=cfg.datatype,
                           root=cfg.deriv_root)
  
      auto_noisy_chs, auto_flat_chs, auto_scores = \
          mne.preprocessing.find_bad_channels_maxwell(
              raw=raw,
              calibration=cfg.mf_cal_fname,
              cross_talk=cfg.mf_ctc_fname,
              origin=mf_head_origin,
              coord_frame='head',
              return_scores=True
          )
  
      preexisting_bads = raw.info['bads'].copy()
      bads = preexisting_bads.copy()
  
      if find_flat_channels_meg:
          msg = f'Found {len(auto_flat_chs)} flat channels.'
          logger.info(**gen_log_kwargs(message=msg,
                                       subject=subject, session=session))
          bads.extend(auto_flat_chs)
      if find_noisy_channels_meg:
          msg = f'Found {len(auto_noisy_chs)} noisy channels.'
          logger.info(**gen_log_kwargs(message=msg,
                                       subject=subject, session=session))
          bads.extend(auto_noisy_chs)
  
      bads = sorted(set(bads))
      raw.info['bads'] = bads
      msg = f'Marked {len(raw.info["bads"])} channels as bad.'
      logger.info(**gen_log_kwargs(message=msg,
                                   subject=subject, session=session))
  
      if find_noisy_channels_meg:
          auto_scores_fname = bids_path.copy().update(
              suffix='scores', extension='.json', check=False)
          with open(auto_scores_fname, 'w') as f:
              json_tricks.dump(auto_scores, fp=f, allow_nan=True,
                               sort_keys=False)
  
          if interactive:
              import matplotlib.pyplot as plt
              plot_auto_scores(auto_scores)
              plt.show()
  
      # Write the bad channels to disk.
      bads_tsv_fname = bids_path.copy().update(suffix='bads',
                                               extension='.tsv',
                                               check=False)
      bads_for_tsv = []
      reasons = []
  
      if find_flat_channels_meg:
          bads_for_tsv.extend(auto_flat_chs)
          reasons.extend(['auto-flat'] * len(auto_flat_chs))
          preexisting_bads = set(preexisting_bads) - set(auto_flat_chs)
  
      if find_noisy_channels_meg:
          bads_for_tsv.extend(auto_noisy_chs)
          reasons.extend(['auto-noisy'] * len(auto_noisy_chs))
          preexisting_bads = set(preexisting_bads) - set(auto_noisy_chs)
  
      preexisting_bads = list(preexisting_bads)
      if preexisting_bads:
          bads_for_tsv.extend(preexisting_bads)
          reasons.extend(['pre-existing (before MNE-BIDS-pipeline was run)'] *
                         len(preexisting_bads))
  
      tsv_data = pd.DataFrame(dict(name=bads_for_tsv, reason=reasons))
      tsv_data = tsv_data.sort_values(by='name')
      tsv_data.to_csv(bads_tsv_fname, sep='\t', index=False)
  
  
  def _load_data(cfg, bids_path):
      # read_raw_bids automatically
      # - populates bad channels using the BIDS channels.tsv
      # - sets channels types according to BIDS channels.tsv `type` column
      # - sets raw.annotations using the BIDS events.tsv
  
      subject = bids_path.subject
      raw = read_raw_bids(bids_path=bids_path)
  
      # Save only the channel types we wish to analyze (including the
      # channels marked as "bad").
      if not use_maxwell_filter:
          picks = get_channels_to_analyze(raw.info)
          raw.pick(picks)
  
      _crop_data(cfg, raw=raw, subject=subject)
  
      raw.load_data()
      if hasattr(raw, 'fix_mag_coil_types'):
          raw.fix_mag_coil_types()
  
      return raw
  
  
  def _crop_data(cfg, raw, subject):
      """Crop the data to the desired duration.
  
      Modifies ``raw`` in-place.
      """
      if subject != 'emptyroom' and cfg.crop_runs is not None:
          raw.crop(*crop_runs)
  
  
  def _drop_channels_func(cfg, raw, subject, session) -> None:
      """Drop channels from the data.
  
      Modifies ``raw`` in-place.
      """
      if cfg.drop_channels:
          msg = f'Dropping channels: {", ".join(cfg.drop_channels)}'
          logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                       session=session))
          raw.drop_channels(cfg.drop_channels)
  
  
  def _create_bipolar_channels(cfg, raw, subject, session) -> None:
      """Create a channel from a bipolar referencing scheme..
  
      Modifies ``raw`` in-place.
      """
      if ch_types == ['eeg'] and cfg.eeg_bipolar_channels:
          msg = 'Creating bipolar channels …'
          logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                       session=session))
          raw.load_data()
          for ch_name, (anode, cathode) in cfg.eeg_bipolar_channels.items():
              msg = f'    {anode} – {cathode} -> {ch_name}'
              logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                           session=session))
              mne.set_bipolar_reference(raw, anode=anode, cathode=cathode,
                                        ch_name=ch_name, drop_refs=False,
                                        copy=False)
          # If we created a new bipolar channel that the user wishes to
          # # use as an EOG channel, it is probably a good idea to set its
          # channel type to 'eog'. Bipolar channels, by default, don't have a
          # location, so one might get unexpected results otherwise, as the
          # channel would influence e.g. in GFP calculations, but not appear on
          # topographic maps.
          if (eog_channels and
                  any([eog_ch_name in cfg.eeg_bipolar_channels
                       for eog_ch_name in eog_channels])):
              msg = 'Setting channel type of new bipolar EOG channel(s) …'
              logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                           session=session))
          for eog_ch_name in eog_channels:
              if eog_ch_name in cfg.eeg_bipolar_channels:
                  msg = f'    {eog_ch_name} -> eog'
                  logger.info(**gen_log_kwargs(message=msg,
                                               subject=subject,
                                               session=session))
                  raw.set_channel_types({eog_ch_name: 'eog'})
  
  
  def _set_eeg_montage(cfg, raw, subject, session) -> None:
      """Set an EEG template montage if requested.
  
      Modifies ``raw`` in-place.
      """
      montage = cfg.eeg_template_montage
      is_mne_montage = isinstance(montage,
                                  mne.channels.montage.DigMontage)
      montage_name = 'custom_montage' if is_mne_montage else montage
      if cfg.datatype == 'eeg' and montage:
          msg = (f'Setting EEG channel locations to template montage: '
                 f'{montage}.')
          logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                       session=session))
          if not is_mne_montage:
              montage = mne.channels.make_standard_montage(montage_name)
          raw.set_montage(montage, match_case=False, on_missing='warn')
  
  
  def _fix_stim_artifact_func(cfg: dict, raw: mne.io.BaseRaw) -> None:
      """Fix stimulation artifact in the data."""
      if not cfg.fix_stim_artifact:
          return
  
      events, _ = mne.events_from_annotations(raw)
      mne.preprocessing.fix_stim_artifact(
          raw, events=events, event_id=None,
          tmin=cfg.stim_artifact_tmin,
          tmax=cfg.stim_artifact_tmax,
          mode='linear'
      )
  
  
  def import_experimental_data(
      *,
      cfg: dict,
      subject: str,
      session: Optional[str] = None,
      run: Optional[str] = None,
      save: bool = False
  ) -> mne.io.BaseRaw:
      """Run the data import.
  
      Parameters
      ----------
      cfg : Bunch
          The local configuration.
      subject : str
          The subject to import.
      session : str | None
          The session to import.
      run : str | None
          The run to import.
      save : bool
          Whether to save the data to disk or not.
  
      Returns
      -------
      raw
          The imported data.
      """
      bids_path_in = BIDSPath(subject=subject,
                              session=session,
                              run=run,
                              task=cfg.task,
                              acquisition=cfg.acq,
                              processing=cfg.proc,
                              recording=cfg.rec,
                              space=cfg.space,
                              suffix=cfg.datatype,
                              datatype=cfg.datatype,
                              root=cfg.bids_root)
      bids_path_out = bids_path_in.copy().update(suffix='raw',
                                                 extension='.fif',
                                                 root=cfg.deriv_root,
                                                 check=False)
  
      raw = _load_data(cfg=cfg, bids_path=bids_path_in)
      _set_eeg_montage(cfg=cfg, raw=raw, subject=subject, session=session)
      _create_bipolar_channels(cfg=cfg, raw=raw, subject=subject,
                               session=session)
      _drop_channels_func(cfg=cfg, raw=raw, subject=subject, session=session)
      _rename_events_func(cfg=cfg, raw=raw, subject=subject, session=session)
      _find_breaks_func(cfg=cfg, raw=raw, subject=subject, session=session,
                        run=run)
      _fix_stim_artifact_func(cfg=cfg, raw=raw)
      _find_bad_channels(cfg=cfg, raw=raw, subject=subject, session=session,
                         task=get_task(), run=run)
  
      # Save the data.
      if save:
          raw.save(fname=bids_path_out, split_naming='bids', overwrite=True)
  
      return raw
  
  
  def import_er_data(
      *,
      cfg: dict,
      subject: str,
      session: Optional[str] = None,
      bads: List[str],
      save: bool = False
  ) -> mne.io.BaseRaw:
      """Import empty-room data.
  
      Parameters
      ----------
      cfg : Bunch
          The local configuration.
      subject
          The subject for whom to import the empty-room data.
      session
          The session for which to import the empty-room data.
      bads
          The selection of bad channels from the corresponding experimental
          recording.
      save
          Whether to save the data to disk or not.
  
      Returns
      -------
      raw_er
          The imported data.
      """
      bids_path_er_in = BIDSPath(
          subject=subject,
          session=session,
          run=cfg.runs[0],
          task=cfg.task,
          acquisition=cfg.acq,
          processing=cfg.proc,
          recording=cfg.rec,
          space=cfg.space,
          suffix=cfg.datatype,
          datatype=cfg.datatype,
          root=cfg.bids_root
      ).find_empty_room()
  
      bids_path_er_out = (bids_path_er_in.copy()
                          .update(task='noise',
                                  run=None,
                                  suffix='raw',
                                  extension='.fif',
                                  root=cfg.deriv_root,
                                  check=False))
  
      if bads is None:
          bads = []
  
      raw_er = _load_data(cfg, bids_path_er_in)
      _drop_channels_func(cfg, raw=raw_er, subject='emptyroom', session=session)
  
      # Set same set of bads as in the experimental run, but only for MEG
      # channels (because we won't have any others in empty-room recordings)
      raw_er.info['bads'] = [ch for ch in bads if ch.startswith('MEG')]
  
      # Save the data.
      if save:
          raw_er.save(fname=bids_path_er_out, split_naming='bids',
                      overwrite=True)
  
      return raw_er
  
  
  def get_reference_run_info(
      *,
      subject: str,
      session: Optional[str] = None,
      run: str
  ) -> mne.Info:
  
      msg = f'Loading info for run: {run}.'
      logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                   session=session, run=run))
  
      bids_path = BIDSPath(
          subject=subject,
          session=session,
          run=run,
          task=get_task(),
          acquisition=acq,
          recording=rec,
          space=space,
          suffix='meg',
          extension='.fif',
          datatype=get_datatype(),
          root=get_bids_root(),
      )
  
      info = mne.io.read_info(bids_path)
      return info
  
  
  def _find_breaks_func(
      *,
      cfg,
      raw: mne.io.BaseRaw,
      subject: str,
      session: Optional[str],
      run: Optional[str],
  ) -> None:
      if not cfg.find_breaks:
          msg = 'Finding breaks has been disabled by the user.'
          logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                       session=session, run=run))
          return
  
      msg = (f'Finding breaks with a mininum duration of '
             f'{cfg.min_break_duration} seconds.')
      logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                   session=session, run=run))
  
      break_annots = mne.preprocessing.annotate_break(
          raw=raw,
          min_break_duration=cfg.min_break_duration,
          t_start_after_previous=cfg.t_break_annot_start_after_previous_event,
          t_stop_before_next=cfg.t_break_annot_stop_before_next_event
      )
  
      msg = (f'Found and annotated '
             f'{len(break_annots) if break_annots else "no"} break periods.')
      logger.info(**gen_log_kwargs(message=msg, subject=subject,
                                   session=session, run=run))
  
      raw.set_annotations(raw.annotations + break_annots)  # add to existing
  
  
  def get_eeg_reference() -> Union[Literal['average'], Iterable[str]]:
      if eeg_reference == 'average':
          return eeg_reference
      elif isinstance(eeg_reference, str):
          return [eeg_reference]
      else:
          return eeg_reference
  
  
  def save_logs(logs):
      fname = get_deriv_root() / f'task-{get_task()}_log.xlsx'
  
      # Get the script from which the function is called for logging
      script_path = pathlib.Path(os.environ['MNE_BIDS_STUDY_SCRIPT_PATH'])
      sheet_name = f'{script_path.parent.name}-{script_path.stem}'
      sheet_name = sheet_name[-30:]  # shorten due to limit of excel format
  
      df = pd.DataFrame(logs)
  
      columns = df.columns
      if "cfg" in columns:
          columns = list(columns)
          idx = columns.index("cfg")
          del columns[idx]
          columns.insert(-3, "cfg")  # put it before time, success & err cols
  
      df = df[columns]
  
      if fname.exists():
          book = load_workbook(fname)
          if sheet_name in book:
              book.remove(book[sheet_name])
          writer = pd.ExcelWriter(fname, engine='openpyxl')
          writer.book = book
      else:
          writer = pd.ExcelWriter(fname, engine='openpyxl')
  
      df.to_excel(writer, sheet_name=sheet_name, index=False)
      writer.save()
      writer.close()
  
  
  # XXX This check should actually go into the CHECKS section, but it depends
  # XXX on get_runs(), which is defined after that section.
  if 'MKDOCS' not in os.environ:
      inter_runs = get_intersect_run()
      mf_ref_error = (
          (mf_reference_run is not None) and
          (mf_reference_run not in inter_runs)
      )
      if mf_ref_error:
          msg = (f'You set mf_reference_run={mf_reference_run}, but your '
                 f'dataset only contains the following runs: {inter_runs}')
          raise ValueError(msg)
  
  
  # Another check that depends on some of the functions defined above
  if (get_task() is not None and
          get_task().lower() != 'rest' and
          conditions is None and
          'MKDOCS' not in os.environ):
      msg = ('Please indicate the name of your conditions in your '
             'configuration. Currently the `conditions` parameter is empty. '
             'This is only allowed for resting-state analysis.')
      raise ValueError(msg)
  