Repository: lina-usc/pylossless (language: python)

Issue:
Implement rejection policy
Conceptually the pipeline has the following parts:
- Processing a file (generate flags)
     - Automated processing
     - QC
-  Loading a pylossless output (reject epochs, channels, and ICs based on the flags)

The rejection should be based on an established "policy" which describes how the flags should be mapped to rejections. We should have an existing default policy that can be changed using a config file. We should probably create a Policy class to manage this. 

E.g. of default policy could be :
- reject automated components not brain or other with probability larger than 0.3 (from the automated process)
- reject manual labels (from the QC process).


Discussion / hints:
@christian-oreilly  Based off todays dev meeting, here is a summary of my understanding, and a proposal. 

We should implement a class, something like `RejectionPolicy`.

Proposal: 

A pipeline instance can have an attribute `rejection_policy`, which is an instance of the `RejectionPolicy` class.

`RejectionPolicy` can have attributes for `channels` and `ICs `. 

Each attribute can br a list of the flag types to exclude. When constructing the `RejectionPolicy` instance, if channels or ic are `None` (default), 
Then `pipeline.rejection_policy`can follow the historical matlab lossless use-case:

`['ch_sd', 'low_r', 'bridged'] `for the channels attribute,

and ['muscle', 'heart', 'ch_noise','eye'] - basically all the iclabels except for brain and other - for the ics attribute. 


(we can make this configurable in the config or something). 


In principle we could add a configurable threshold iclabel confidence rate, (so users t don't have to accept the 30% confidence threshold), which must be met for the ic to be rejected, but thinking that about 95% of use cases probably won't change it maybe we wait until this is needed?


finally, when the user runs functions like `launch_QC` or `apply_flags`, those functions can take the rejection policy instance and use it to color code flagged ICs and flagged channels (QC dashboard) or add pipeline flags to MNE bads (`apply_flags`).

Is that what you had in mind?


Fix patch:
diff --git a/.github/workflows/build_doc.yml b/.github/workflows/build_doc.yml
index 4846900..f9076f5 100644
--- a/.github/workflows/build_doc.yml
+++ b/.github/workflows/build_doc.yml
@@ -23,6 +23,9 @@ jobs:
       - name: Install Pylossless
         run: pip install -e .
 
+      - name: install openneuro-py
+        run: pip install openneuro-py
+
       - name: Install doc dependencies
         run: pip install -r docs/requirements_doc.txt
 
diff --git a/.gitignore b/.gitignore
index 4d18513..f657447 100644
--- a/.gitignore
+++ b/.gitignore
@@ -135,4 +135,7 @@ my_project_ll_config.yaml
 project_ll_config_face13_egi.yaml
 
 # VSCode ignore
-.vscode
\ No newline at end of file
+.vscode
+
+# sphinx gallery
+./docs/source/auto_examples/
\ No newline at end of file
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 3777a26..1470db4 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -7,9 +7,9 @@ version: 2
 
 # Set the version of Python and other tools you might need
 build:
-  os: ubuntu-22.04
-  tools:
-    python: "3.10"
+   os: ubuntu-22.04
+   tools:
+      python: "3.10"
 
 # Build documentation in the docs/ directory with Sphinx
 sphinx:
@@ -18,5 +18,5 @@ sphinx:
 # Optionally declare the Python requirements required to build your docs
 python:
    install:
-   - requirements: requirements.txt
-   - requirements: docs/requirements_doc.txt
\ No newline at end of file
+      - requirements: requirements.txt
+      - requirements: docs/requirements_doc.txt
diff --git a/docs/Makefile b/docs/Makefile
index d0c3cbf..e34a181 100644
--- a/docs/Makefile
+++ b/docs/Makefile
@@ -8,6 +8,10 @@ SPHINXBUILD   ?= sphinx-build
 SOURCEDIR     = source
 BUILDDIR      = build
 
+clean:
+	-rm -rf $(BUILDDIR)/*
+	-rm -rf $(SOURCEDIR)/auto_examples/
+
 # Put it first so that "make" without argument is like "make help".
 help:
 	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
diff --git a/docs/examples/plot_0_implementation.py b/docs/examples/plot_0_implementation.py
deleted file mode 100644
index f03fc2c..0000000
--- a/docs/examples/plot_0_implementation.py
+++ /dev/null
@@ -1,750 +0,0 @@
-r"""
-
-PyLossless Algorithms
-=====================
-
-This tutorial explains the calculations that PyLossless performs at each step of the
-pipeline. We will use example EEG data to demonstrate the
-calculations.
-
-.. note::
-    You can open this notebook in
-    `Google Colab <https://colab.research.google.com/drive/1ecyNo10oFgpbVNuD7Ztgr2fs8XkpfOYo?usp=sharing>`_!
-
-.. _notation:
-
-Notation 
---------
-Before we begin, we define some notation that will be used throughout the text:
-
-- We start with a 3D matrix of EEG data,
-  :math:`X \in \mathbb{R}^{S_\mathcal{G} \times E_\mathcal{G} \times T}`,
-  where :math:`S_\mathcal{G}` and :math:`E_\mathcal{G}` are the sets of good sensors and
-  epochs, respectively, and :math:`T`, is the number of samples(i.e. time-points).
-
-- ``s``, ``e``, and ``t`` are sensor, epochs, and samples, respectively.
-
-- We use superscripts to denote operations across a dimension, and we use subscripts to
-  denote indexing a dimension.
-
-- We refer to a single sensor :math:`i` as
-  :math:`X\big|_{s=i} \in \mathbb{R}^{E_\mathcal{G} \times T}`,
-  with :math:`i \in S_\mathcal{G}`.
-
-- We refer to a single epoch :math:`j` as
-  :math:`X\big|_{e=j} \in \mathbb{R}^{S_\mathcal{G} \times T}`,
-  with :math:`j \in E_\mathcal{G}`.
-
-- We denote sensor-specific thresholds for rejecting epochs as
-  :math:`\tau^e_i \in \mathbb{R}^{S_\mathcal{G}}`
-
-- We denote epoch-specific thresholds for rejecting sensors as
-  :math:`\tau^s_j \in \mathbb{R}^{E_\mathcal{G}}`
-
-- We denote *quantiles* as :math:`Q\#^{dim}`: i.e. :math:`Q75^s` is the 75th *quantile*
-  along the sensor dimension. The function :math:`Q75^s(X)` computes the 75th quantile
-  along the :math:`s` dimension of matrix :math:`X`, resulting in a matrix noted
-  :math:`X^{Q75^s} \in \mathbb{R}^{E \times T}`.
-
-Throughout the text, we use capital letters for matrices and lowercase letters for
-scalars. For example, the data point for sensor :math:`i`, epoch :math:`j`, and
-time :math:`k` is denoted as :math:`X\big|_{s=i; e=j; t=k} = x_{ijk} \in \mathbb{R}`,
-and :math:`X=\{x_{ijk}\}`.
-"""
-
-# %%
-# Imports and data loading
-# ------------------------
-from pathlib import Path
-
-import numpy as np
-
-import mne
-from mne.datasets import sample
-
-import pylossless as ll
-
-# Load example mne data
-raw = ll.datasets.load_simulated_raw()
-
-# %%
-# Load a PyLossless configuration file
-# ------------------------------------
-# Let's load a PyLossless configuration file. This file contains the parameters that
-# will be used for each step of the pipeline. For example, the ``ch_ch_sd`` section
-# contains the parameters for the :ref:`noisy_sensors` step. We can modify these
-# parameters to change the behavior of the pipeline. For example, we can change the
-# percent of epochs that a sensor must be noisy for it to be flagged via the
-# ``flag_crit`` parameter.
-config = ll.config.Config()
-config.load_default()
-config["ch_ch_sd"]["outliers_kwargs"]["lower"] = 0.25  # lower quantile
-config["ch_ch_sd"]["outliers_kwargs"]["upper"] = 0.75  # upper quantile
-config["ch_ch_sd"]["flag_crit"] = 0.30  # percent of epochs that a sensor must be noisy
-config.save("test_config.yaml")
-
-# %%
-# Create a pipeline instance
-# --------------------------
-pipeline = ll.LosslessPipeline("test_config.yaml")
-pipeline.raw = raw
-raw.plot()
-
-# %%
-# Input Data
-# ----------
-#
-# First, we epoch the data to be used for subsequent steps.
-# Let our 3D matrix below be defined as :math:`X \in \mathbb{R}^{S \times E \times T}`
-# where :math:`X` is a matrix of real numbers and of dimension :math:`S` sensors
-# :math:`\times$ E` epochs `\times T` times.
-#
-epochs = pipeline.get_epochs()
-
-# %%
-#
-# Let's convert our epochs object into a named :class:`xarray.DataArray` object.
-from pylossless.pipeline import epochs_to_xr
-
-#
-epochs_xr = epochs_to_xr(epochs, kind="ch")
-epochs_xr  # 277 epochs, 50 sensors, 602 samples per epoch
-
-
-# %%
-# .. _robust_reference:
-#
-# Robust Average Reference
-# ------------------------
-#
-# .. figure:: https://raw.githubusercontent.com/scott-huberty/wip_pipeline-figures/main/robust_rereference.png
-#    :align: center
-#    :alt: Robust Average Reference graphic.
-#
-#    Robust Average Reference. The figure shows the steps for robust average referencing.
-#    See the text below for descriptions of mathematical notation.
-#
-# Before the pipeline can begin, we must average reference the data. This is because
-# the pipeline uses data distributions to identify noisy sensors, and For EEG data that
-# uses an online reference to a single electrode, sensors that are further from the
-# reference will have a higher voltage variance, and the pipeline will be biased to
-# flag these sensors as noisy. The average reference, which subtracts the average
-# signal across sensors from each individual sensor, will ensure an even playing field.
-# Howeer, we dont want to include noisy sensors in the average reference signal. So we
-# will identify noisy sensors and and leave them out of the average reference signal.
-
-# %%
-sample_std = epochs_xr.std("time")
-q25_ch = sample_std.quantile(0.25, dim="ch")
-q50_ch = sample_std.median(dim="ch")
-q75_ch = sample_std.quantile(0.75, dim="ch")
-ch_dist = sample_std - q50_ch  # center the data
-ch_dist /= q75_ch - q25_ch  # shape (chans, epoch)
-
-mean_ch_dist = ch_dist.mean(dim="epoch")  # shape (chans)
-
-# find the median and 25 and 75 percentiles
-# of the mean of the channel distributions
-mdn = np.median(mean_ch_dist)
-deviation = np.diff(np.quantile(mean_ch_dist, [0.25, 0.75]))
-
-leave_out = mean_ch_dist.ch[mean_ch_dist > mdn + 6 * deviation].values.tolist()
-leave_out
-
-# %%
-ref_chans = [ch for ch in epochs.pick("eeg").ch_names if ch not in leave_out]
-pipeline.raw.set_eeg_reference(ref_channels=ref_chans)
-
-# %%
-#
-# .. _noisy_sensors:
-#
-# Flag Noisy Sensors
-# ------------------
-# .. figure:: https://raw.githubusercontent.com/scott-huberty/wip_pipeline-figures/main/Flag_noisy_sensors.png
-#    :align: center
-#    :alt: Flag Noisy Sensors graphic.
-#
-#    Flag Noisy Sensors. The figure shows the steps for flagging noisy sensors. See the text below
-#    for descriptions of mathematical notation.
-#
-
-# %%
-# Since we applied a robust average reference to the raw data, we will need to re-epoch
-# the data:
-epochs = pipeline.get_epochs()
-epochs_xr = epochs_to_xr(epochs, kind="ch")
-
-# First we take standard deviation of
-# :math:`X \in \mathbb{R}^{S \times E \times T}` across the samples dimension :math:`t`
-# resulting in a 2D matrix :math:`X^{\sigma_{t}} \in \mathbb{R}^{S \times E}`
-trim_ch_sd = epochs_xr.std("time")
-trim_ch_sd
-
-# %%
-# a) Take the 50th and 75th quantile across dimension sensor of :math:`X^{\sigma_{t}}`
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# This operation results in two 1D vectors of size :math:`E`:
-#
-# .. math::
-#    X^{{\sigma}_t{Q50^s}} = Q50^s(X^{\sigma_{t}}) \in \mathbb{R}^{E}
-# .. math::
-#    X^{{\sigma}_t{Q75^s}} = Q75^s(X^{\sigma_{t}}) \in \mathbb{R}^{E}
-
-# %%
-q50, q75 = trim_ch_sd.quantile([0.5, 0.75], dim="ch")
-q50  # a 1D array of median standard deviation values across channels for each epoch
-
-# %%
-# b) Define an Upper Quantile Range as :math:`Q75 - Q50`
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# .. math::
-#    UQR^s = X^{{\sigma}_T{Q75}^s} - X^{{\sigma}_T{Q50}^s}
-#
-# This operation results in a 1D vector of size :math:`E`.
-uqr = q75 - q50
-uqr
-
-# %%
-# c) Identify outlier Indices :math:`(i, j)`
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# We multiply a constant :math:`k` by the :math:`UQR` to define a measure for the
-# spread of the right tail of the distribution of :math:`X^{\sigma_{t}}` values and
-# add it to the median of :math:`X^{\sigma_{t}}` to obtain epoch-specific standard
-# deviation threshold for outliers:
-#
-# .. math::
-#    \tau^s_j = X^{{\sigma}_T{Q50}^S} + UQR^s\times k
-#
-# That is, :math:`\tau^s_j` is the epoch-specific threshold for the epoch :math:`j`
-k = 3
-upper_threshold = q50 + q75 * k
-upper_threshold  # epoch specific thresholds
-
-# %%
-# Now, we compare our 2D standard deviation matrix to the threshold vector of
-# :math:`\tau^e_j`:
-#
-# .. math::
-#    X^{\sigma_{t}} \big|_{e=j}  > \tau^s_j
-#
-# resulting in the indicator matrix :math:`C \in \{0, 1\}^{S \times E}=\{c_{ij}\}`:
-#
-# .. math::
-#    c_{ij} =
-#    \begin{cases}
-#    0 & \text{if } x^{\sigma_{t}}_{ij} < \tau^s_j \\
-#    1 & \text{if } x^{\sigma_{t}}_{ij} \geq \tau^s_j
-#    \end{cases}
-#
-# Each element of this matrix indicates whether sensor :math:`i` is an outlier at an epoch
-# :math:`j`.
-outlier_mask = trim_ch_sd > upper_threshold
-outlier_mask
-
-# %%
-# d) Identify noisy sensors part 1
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# To identify outlier sensors, we average across the epoch dimension of our indicator
-# matrix :math:`C` and obtain :math:`C^{\mu_e} \in \mathbb{R}^{S_\mathcal{G}}`, which
-# is a vector of fractional numbers :math:`c^{\mu_e}_i` representing the percentage of
-# epochs for which that sensor is an outlier.
-percent_outliers = outlier_mask.astype(float).mean("epoch")
-percent_outliers  # percent of epochs that sensor is an outlier
-
-# %%
-# e) Identify noisy sensors part 2
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# Next, we define a threshold :math:`\tau^{p}` (:math:`p` for percentile;
-# default, ``.20``) as a cutoff point for determining if a sensor should be marked
-# artifactual. The sensor :math:`i` is flagged as noisy if
-# :math:`c^{\mu_e}_i > \tau^{p}`. That is, if the sensor is an outlier for more than
-# :math:`\tau^{p}` percent of the epochs, it is flagged as noisy.
-p_threshold = config["ch_ch_sd"]["flag_crit"]  # 0.3, or 30%
-noisy_chs = percent_outliers[percent_outliers > p_threshold].coords.to_index().values
-noisy_chs
-
-# %%
-# f) Add the noisy sensors to the pipeline flags
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# Let's add the noisy sensors to the pipeline flags.
-pipeline.flags["ch"].add_flag_cat(kind="ch_sd", bad_ch_names=noisy_chs)
-pipeline.raw.info["bads"].extend(pipeline.flags["ch"]["ch_sd"].tolist())
-pipeline.flags["ch"]
-
-# %%
-#
-# .. _noisy_epochs:
-#
-# Flag Noisy Epochs
-# -----------------
-#
-# This step closely resembles the :ref:`noisy_sensors` step. For the sake of brevity
-# we will be more concise in the documentation.
-
-# %%
-#
-# .. figure:: https://raw.githubusercontent.com/scott-huberty/wip_pipeline-figures/main/Flag_noisy_epochs.png
-#    :align: center
-#
-#    Flag Noisy Epochs. The figure shows the steps for flagging noisy epochs. See the text below
-#    for descriptions of mathematical notation.
-#
-
-# %%
-# a) Take standard deviation across the samples dimension :math:`t`
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# Take a moment below to notice that the sensors flagged in the prior setp are not in
-# ``epochs_xr`` below:
-epochs = pipeline.get_epochs()
-# Let's make our epochs array into a named Array
-epochs_xr = epochs_to_xr(epochs, kind="ch")
-trim_ch_sd = epochs_xr.std("time")
-trim_ch_sd.coords["ch"]
-
-# %%
-# b) Compute 50th and 75th quantile across epochs and the UQR
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# Like before, We Take the median and 70th quantile, but now we operate across epochs,
-# resulting in two 1D vector's of size ``n_good_sensors`` :math:`S_\mathcal{G}`
-#
-# .. math::
-#    X^{{\sigma}_t{Q50^e}} = Q50^e(X^{\sigma_{t}}) \in \mathbb{R}^{S_\mathcal{G}}
-# .. math::
-#    X^{{\sigma}_t{Q75^e}} = Q75^e(X^{\sigma_{t}}) \in \mathbb{R}^{S_\mathcal{G}}
-# .. math::
-#    UQR^e = (X^{{\sigma}_T{Q75}^e} - X^{{\sigma}_T{Q50}^e})
-
-
-# %%
-# c) Define sensor-specific thresholds for rejecting epochs :math:`\tau^e_i`
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# Our sensor-specifc threshold for rejecting epochs is defined by:
-#
-# .. math::
-#    \tau^e_i = X^{{\sigma}_T{Q50}^e} + UQR^e\times k
-q50, q75 = trim_ch_sd.quantile([0.5, 0.75], dim="epoch")
-uqr_epoch = q75 - q50
-uqr_epoch
-
-# %%
-k = 8
-upper_threshold = q50 + uqr_epoch * k
-upper_threshold
-
-# %%
-# d) Identify Outlier indices
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# The indicator matrix is defined by:
-#
-# .. math::
-#    c_{ij} =
-#    \begin{cases}
-#    0 & \text{if } x^{\sigma_{t}}_{ij} < \tau^e_i \\
-#    1 & \text{if } x^{\sigma_{t}}_{ij} \geq \tau^e_i
-#    \end{cases}
-#
-#
-# To identify outlier **epochs**, we average across the **sensor** dimension of our
-# indicator matrix :math:`C` and obtain
-# :math:`C^{\mu_s} \in \mathbb{R}^{E_\mathcal{G}}`, which is a vector of numbers
-# :math:`c^{\mu_s}_j` representing the percentage of **sensors** for which that epoch
-# is an outlier.
-outlier_mask = trim_ch_sd > upper_threshold
-outlier_mask
-
-# %%
-percent_outliers = outlier_mask.astype(float).mean("ch")
-percent_outliers
-
-# %%
-# e) Identify noisy epochs
-# ^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# Next, we define a fractional threshold :math:`\tau^{p}` as a cutoff point for
-# determining if a epoch should be marked artifactual. The epoch :math:`j` is flagged
-# as noisy if :math:`c^{\mu_s}_j > \tau^{p}`.
-bad_epochs = percent_outliers[percent_outliers > p_threshold].coords.to_index().values
-bad_epochs
-
-# %%
-# f) Add the noisy epochs to the pipeline flags
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# Let's add the outlier epochs to our flags
-# These will be added directly as :class:`mne.Annotations` to the raw data.
-pipeline.flags["epoch"].add_flag_cat(
-    kind="ch_sd", bad_epoch_inds=bad_epochs, epochs=epochs
-)
-pipeline.raw.annotations.description
-
-# %%
-pipeline.raw.plot()
-
-# %%
-# Filtering
-# ---------
-#
-# After flagging noisy sensors and epochs, we filter the data. By default,
-# The pipeline uses a 1-100Hz bandpass filter. This is because 1), ICA decompositions
-# are more stable when low frequency drifts are removed, and 2) the ICLabel classifier
-# is trained on data that has been filtered between 1-100Hz. A notch filter can also be
-# optionally specified.
-pipeline.config["filtering"]["notch_filter_args"]["freqs"] = [50]
-pipeline.filter()
-
-# %%
-# Find Nearest Neighbours & return Maximum Correlation
-# ----------------------------------------------------
-#
-# .. figure:: https://raw.githubusercontent.com/scott-huberty/wip_pipeline-figures/main/Nearest_neighbors.png
-#    :align: center
-#    :alt: Nearest Neighbors graphic.
-#
-#    Nearest Neighbors. The figure shows the steps for finding nearest neighbors. See the text below
-#    for descriptions of mathematical notation.
-#
-
-# %%
-# Whereas :ref:`noisy_sensors` and :ref:`noisy_epochs` operated on a 2D matrix of
-# standard deviation values, The next few steps will operate on correlation
-# coefficients. Here we describe the procedure for defining the 2D matrix of correlation
-# coefficients.
-
-# %%
-from pylossless.pipeline import chan_neighbour_r
-
-# %%
-#
-# Notice that our flagged epochs are dropped.
-epochs = pipeline.get_epochs()
-
-# %%
-# a) Calculate Correlation Coefficients between each Sensor and its neighboring eighbors
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# - For each good sensor $i$ in :math:`S_{\mathcal{G}}`, we select its :math:`N` nearest
-#   neighbors. I.e. the :math:`N` sensors that are closest to it.
-#
-# - We call the sensor :math:`i` the *origin*, and its nearest neighbors :math`\hat{s_l}`
-#   with :math:`l \in \{1, 2, \ldots, N\}`
-#
-# - Then, for each epoch :math:`j`, we calculate the correlation coefficient
-#   :math:`\rho^t_{(i,\hat{s_l}),j}` between origin sensor :math:`i` and each neighbor
-#   :math:`\hat{s_l}` across dimension :math:`t` (samples), returning a 3D matrice of
-#   correlation coefficients:
-#
-# .. math::
-#    \mathrm{P}^t = \{\rho^t_{(i, \hat{s_l}),j}\} \in \mathbb{R}^{S_G \times E_G \times n}
-#
-# Finally, we select the maximum correlation coefficient across the neighbor dimension
-# :math:`n`:
-#
-# .. math::
-#    \mathrm{P}^{t,{\text{max}}^n}= \max\limits_{\hat{s_l}}  \rho^t_{(i, \hat{s_l}),j}
-#
-# Returning a 2D matrix where each value at :math:`(i, j)` is the maximum correlation
-# coefficient between sensor :math:`i` and its :math:`N` nearest neighbors, at each epoch
-# :math:`j`
-
-# %%
-data_r_ch = chan_neighbour_r(epochs, nneigbr=3, method="max")
-# maximum correlation out of correlations between ch and its 3 neighbors
-data_r_ch
-
-# %%
-# This matrix :math:`\mathrm{P}^{t,{\text{max}}^n}`  will be used in the steps below.
-
-# %%
-# Flag Bridged Sensors
-# --------------------
-#
-# .. figure:: https://raw.githubusercontent.com/scott-huberty/wip_pipeline-figures/main/Flag_bridged_sensors.png
-#    :align: center
-#    :alt: Flag Bridged Sensors graphic.
-#
-#    Flag Bridged Sensors. The figure shows the steps for flagging bridged sensors.
-#    See the text below for descriptions of mathematical notation.
-#
-
-# %%
-# a) Calculate the 50th, 75th quantile and IQR across epochs
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# .. math::
-#    IQR^e = \mathrm{P}^{t,{\text{max}}^nQ75^e} - \mathrm{P}^{t,{\text{max}}^nQ25^e}
-#
-# For each sensor, divide the median across epochs by the IQR across epochs. Bridged
-# channels should have a high median correlation but a low IQR of the correlation.
-# We call this measure the bridge-indicator.
-#
-# .. math::
-#    \mathcal{B}_s = \frac{\mathrm{P}^{t,{\text{max}}^nQ50^e}}{IQR^e}
-#
-# b) Define a bridging threshold
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-# Now, take the 25th, 50th, and 75th quantile of :math:`\mathcal{B}_s` across sensors,
-# And calculate the :math:`IQR^s`. A channel :math:`i` is bridged if
-#
-# .. math::
-#    \mathcal{B}_i > B^{Q50^s} +k \times IQR^s
-#
-
-# %%
-import scipy
-from functools import partial
-
-# %%
-msr = data_r_ch.median("epoch") / data_r_ch.reduce(scipy.stats.iqr, dim="epoch")
-# msr is a 1D vector of size n_sensors
-config_trim = 40
-config_bridge_z = 6
-#
-trim = config_trim
-if trim >= 1:
-    trim /= 100
-trim /= 2  # .20 and will be used as (.20, .20)
-#
-trim_mean = partial(scipy.stats.mstats.trimmed_mean, limits=(trim, trim))
-trim_std = partial(scipy.stats.mstats.trimmed_std, limits=(trim, trim))
-#
-z_val = config_bridge_z  # 6
-mask = msr > msr.reduce(trim_mean, dim="ch") + z_val * msr.reduce(
-    trim_std, dim="ch"
-)  # bridged chans
-#
-bridged_ch_names = data_r_ch.ch.values[mask]
-bridged_ch_names
-
-# %%
-# Let's add the outlier channels to our flags
-bad_chs = bridged_ch_names
-pipeline.flags["ch"].add_flag_cat(kind="bridge", bad_ch_names=bad_chs)
-pipeline.flags["ch"]
-
-# %%
-# Identify the Rank Channel
-# -------------------------
-#
-# Because the pipeline uses an average reference before the ICA decomposition, it is
-# necessary to account for rank deficiency (i.e., every sensor in the montage is
-# linearly dependent on the other channels due to the common average reference). To
-# account for this, the pipeline flags the sensor (out of the remaining good sensors)
-# with the highest median of the max correlation coefficient with its neighbors
-# (across epochs):
-#
-# .. math::
-#    \begin{equation}
-#    i = \text{arg}\max\limits_i \rho_{i}^{t,{\text{max}}^n,median^j}
-#    \end{equation}
-#
-# This sensor has the least unique time-series out of the remaining set of good sensors
-# :math:`S_\mathcal{G}` and is flagged by the pipeline as ``”rank”``. Note that this
-# sensor is not flagged because it contains artifact, but only because one of the
-# remaining sensors needs to be removed to address rank deficiency before ICA
-# decomposition is performed. By choosing this sensor, we are likely to lose little
-# information because of its high correlation with its neighbors. This sensor can be
-# reintroduced after the ICA has been applied for artifact corrections.
-
-# %%
-good_chs = [
-    ch for ch in data_r_ch.ch.values if ch not in pipeline.flags["ch"].get_flagged()
-]
-data_r_ch_good = data_r_ch.sel(ch=good_chs)
-
-flag_ch = [str(data_r_ch_good.median("epoch").idxmax(dim="ch").to_numpy())]
-pipeline.flags["ch"].add_flag_cat(kind="rank", bad_ch_names=flag_ch)
-pipeline.flags["ch"]
-
-# %%
-# Flag low correlation Epochs
-# ---------------------------
-#
-# This step is designed to identify time periods in which many sensors are
-# uncorrelated with neighboring sensors. It is similar to the :ref:`noisy_sensors` step,
-#
-# Again we calculate the 25th and 50th quantile
-# of :math:`\mathrm{P}^{t,{\text{max}}^n}`, across the epochs dimension, and calculate
-# the lower quantile range :math:`LQR^s`. This results in vectors
-# :math:`\mathrm{P}^{t,{\text{max}}^nQ25^e}` and
-# :math:`\mathrm{P}^{t,{\text{max}}^nQ50^e}` of size :math:`S_\mathcal{G}`. As for previous
-# steps, we define  sensor-specific thresholds for flagging epochs:
-#
-# .. math::
-#    \begin{equation}
-#    \tau^e = \mathrm{P}^{t,{\text{max}}^nQ50^e} - LQR^e\times k
-#    \end{equation}
-#
-# And the corresponding indicator matrix:
-#
-# .. math::
-#    \begin{equation}
-#    c_{ij} =
-#    \begin{cases}
-#    1 & \text{if } \rho^{t,{\text{max}}^n}_{ij} < \tau^e_i \\
-#    0 & \text{if } \rho^{t,{\text{max}}^n}_{ij} \geq \tau^e_i
-#    \end{cases}
-#    \end{equation}
-#
-# We average the indicator matrix across sensors and obtain a vector :math:`C^{\mu_s}`
-# that we use to flag uncorrelated epochs using the following criterion:
-#
-# .. math::
-#    c^{\mu_e}_i > \tau^{p}.
-#
-
-# %%
-# Step a
-q25, q50 = data_r_ch.quantile([0.25, 0.5], dim="epoch")
-#
-# Define the LQR
-lqr = q50 - q25
-#
-# define a threshold
-k = 3
-lower_threshold = q50 - lqr * k
-#
-outlier_mask = data_r_ch < lower_threshold
-#
-percent_outliers = outlier_mask.astype(float).mean("ch")
-#
-p_threshold = 0.2
-bad_epochs = percent_outliers[percent_outliers > p_threshold].coords.to_index().values
-#
-# Add the outlier epochs to our flags
-pipeline.flags["epoch"].add_flag_cat(
-    kind="uncorrelated", bad_epoch_inds=bad_epochs, epochs=epochs
-)
-pipeline.raw.annotations.description
-
-# %%
-# in this case, no epochs were flagged as uncorrelated.
-#
-
-# %%
-# Flag low correlation Sensors
-# -----------------------------
-#
-# .. figure:: https://raw.githubusercontent.com/scott-huberty/wip_pipeline-figures/main/Flag_uncorrelated_sensors.png
-#    :align: center
-#    :alt: Flag Uncorrelated Sensors graphic.
-#
-#    Flag Uncorrelated Sensors. The figure shows the steps for flagging uncorrelated
-#    sensors. See the text below for descriptions of mathematical notation.
-#
-# This step is designed to identify sensors that have an unusually low correlation with
-# neighboring sensors. The operations involved by this step are similar to those of the
-# :ref:`noisy_sensors` step, except we use maximal nearest neighbor correlations instead
-# of dispersion and the left instead of the right tail of the distribution to set
-# the threshold for outliers.
-
-# %%
-# a) Take lower quantile range and defined sensor-specific thresholds
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# We get the indicator matrix as described previously, using
-#
-# .. math::
-#    \tau^e_i = \mathrm{P}^{t,{\text{max}}^nQ50^e} - LQR^e\times k
-#
-# and
-#
-# .. math::
-#    c_{ij} =
-#    \begin{cases}
-#    1 & \text{if } \rho^{t,{\text{max}}^n}_{ij} < \tau^e_i \\
-#    0 & \text{if } \rho^{t,{\text{max}}^n}_{ij} \geq \tau^e_i
-#    \end{cases}
-
-# %%
-# b) Identify uncorrelated sensors
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-#
-# We define a threshold as we did in the previous step and flag uncorrelated epochs
-# :math:`j` if :math:`c^{\mu_s}_j > \tau^{p}`.
-
-# %%
-q25, q50 = data_r_ch.quantile([0.25, 0.5], dim="ch")
-#
-# Define LQR
-lqr = q50 - q25
-#
-# define a threshold
-k = 3
-lower_threshold = q50 - lqr * k
-#
-# Identify correlations less than the threshold
-outlier_mask = data_r_ch < lower_threshold
-percent_outliers = outlier_mask.astype(float).mean("epoch")
-#
-p_threshold = 0.2
-bad_chs = percent_outliers[percent_outliers > p_threshold].coords.to_index().values
-#
-# Add the outlier channels to our flags
-pipeline.flags["ch"].add_flag_cat(kind="uncorrelated", bad_ch_names=bad_chs)
-pipeline.flags["ch"]
-
-# %%
-# In this case, no sensors were flagged as uncorrelated.
-
-# %%
-# Run Initial ICA
-# ---------------
-#
-# The pipeline by runs ICA two times. The first ICA is only used to identify
-# noisy periods in its IC activation time-series. For this reason, the pipeline
-# uses the FastICA algorithm for speed.
-
-# %%
-pipeline.run_ica("run1")
-
-# %%
-# Flag Noisy IC Activation time-periods
-# -------------------------------------
-#
-# This step follows the same procedure as the :ref:`noisy_sensors` step, except that
-# the data is now the IC activation time-series. thus we start with a 3D matrix
-# :math:`X_{ica} \in \mathbb{R}^{I_\mathcal{G} \times E_\mathcal{G} \times T}` of
-# IC time-courses rather than scalp EEG data and where :math:`I` is the set of
-# independent components.
-#
-
-# %%
-pipeline.flag_epoch_ic_sd1()
-
-# %%
-pipeline.raw.annotations.description
-
-# %%
-# Run Final ICA
-# -------------
-#
-# Now The pipeline runs the final ICA decomposition, this time using the extended
-# Infomax algorithm. Note that any sensors or time-periods that have been flagged
-# up to this point will not be passed into the ICA decomposition. For the sake of
-# time, we will not run the second ICA here, as there are no more pipeline calculations.
-
-# %%
-# Run ICLabel Classifier
-# ----------------------
-#
-# The pipeline will run the ICLabel classifier on the final ICA, which will produce a
-# label for each IC, one of ``"brain"``, ``"muscle"``, ``"eog"`` (eye), ``"ecg"``
-# (heart), ``line_noise``, or ``"channel_noise"``.
-#
-#
-# Conclusion
-# ----------
-# And that's all! See the other pylossless tutorials for brief examples on running the
-# pipeline on your own data, and rejecting the flagged data.
-#
diff --git a/docs/requirements_doc.txt b/docs/requirements_doc.txt
index 884bce8..02ec517 100644
--- a/docs/requirements_doc.txt
+++ b/docs/requirements_doc.txt
@@ -1,3 +1,4 @@
+openneuro-py
 sphinx>=6
 shibuya
 sphinx-design
diff --git a/docs/source/API/API_index.rst b/docs/source/API/API_index.rst
index da7eeda..6f8100d 100644
--- a/docs/source/API/API_index.rst
+++ b/docs/source/API/API_index.rst
@@ -14,12 +14,14 @@ the ``lossless`` namespace.
 
    config_class
    bids_funcs
+   datasets
    lossless_class
    pipeline_funcs
    flag_classes
 
 - :doc:`bids_funcs`
 - :doc:`config_class`
+- :doc:`datasets`
 - :doc:`lossless_class`
 - :doc:`pipeline_funcs`
 - :doc:`flag_classes`
diff --git a/docs/source/API/config_class.rst b/docs/source/API/config_class.rst
index e54830b..6bef3c6 100644
--- a/docs/source/API/config_class.rst
+++ b/docs/source/API/config_class.rst
@@ -1,11 +1,21 @@
-PyLossless Config Class
-=======================
+PyLossless Configuration Classes
+================================
 
 This is the docstring for the Lossless Config class. The methods and
 attributes detailed below are used to parse the YAML config file that is read
 by the Lossless Pipeline, and which contains the parameters to use for the
 pipeline.
 
-.. todo::
-    Format config docstring following Numpy Format.
+PyLossless Config Class
+=======================
+
+.. autoclass:: pylossless.config.Config
+    :members:
+
+Rejection Config Class
+======================
+
+.. autoclass:: pylossless.config.RejectionPolicy
+    :members:
+
 
diff --git a/docs/source/API/datasets.rst b/docs/source/API/datasets.rst
new file mode 100644
index 0000000..f4664c3
--- /dev/null
+++ b/docs/source/API/datasets.rst
@@ -0,0 +1,5 @@
+Datasets
+========
+
+.. automodule:: pylossless.datasets.datasets
+    :members:
\ No newline at end of file
diff --git a/docs/source/conf.py b/docs/source/conf.py
index ef1f24f..20ab533 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -47,8 +47,8 @@
 # Source directory of python file tutorials and the target
 # directory for the converted rST files
 sphinx_gallery_conf = {
-    "examples_dirs": "../examples",  # path to your example scripts
-    "gallery_dirs": "./generated/auto_tutorials",  # path to save tutorials
+    "examples_dirs": "../../examples",  # path to your example scripts
+    "gallery_dirs": "auto_examples",  # path to save tutorials
 }
 
 templates_path = ["_templates"]
@@ -63,6 +63,7 @@
     "xarray": ("https://docs.xarray.dev/en/stable/", None),
     "mne": ("https://mne.tools/dev", None),
     "mne_icalabel": ("https://mne.tools/mne-icalabel/dev", None),
+    "mne_bids": ("https://mne.tools/mne-bids/dev", None),
     "pylossless": ("https://pylossless.readthedocs.io/en/latest/", None),
 }
 
@@ -88,16 +89,13 @@
 # NumPyDoc configuration -----------------------------------------------------
 
 numpydoc_class_members_toctree = False
-numpydoc_show_inherited_class_members = {
-    "FlaggedChs": False,
-    "FlaggedEpochs": False,
-    "FlaggedICs": False,
-}
+numpydoc_show_inherited_class_members = False
 numpydoc_attributes_as_param_list = True
 numpydoc_xref_param_type = True
 numpydoc_validate = True
 # Only generate documentation for public members
 autodoc_default_flags = ["members", "undoc-members", "inherited-members"]
+
 numpydoc_class_members_toctree = False
 
 numpydoc_xref_aliases = {
diff --git a/docs/source/generated/index.rst b/docs/source/generated/index.rst
deleted file mode 100644
index 53b2c07..0000000
--- a/docs/source/generated/index.rst
+++ /dev/null
@@ -1,11 +0,0 @@
-PyLossless Tutorials
-====================
-
-View the gallery below to see example tutorials using the pyLossless Pipeline.
-
-.. toctree::
-   :maxdepth: 1
-   :hidden:
-
-.. include:: auto_tutorials/index.rst
-
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 89e8e5e..f97c335 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -11,8 +11,8 @@ PyLossless EEG Processing Pipeline
    :hidden:
 
    install
+   auto_examples/index.rst
    API/API_index
-   generated/index
    contributing
 
 
diff --git a/docs/examples/README.txt b/examples/README.txt
similarity index 100%
rename from docs/examples/README.txt
rename to examples/README.txt
diff --git a/examples/plot_10_run_pipeline.py b/examples/plot_10_run_pipeline.py
new file mode 100644
index 0000000..9d50a28
--- /dev/null
+++ b/examples/plot_10_run_pipeline.py
@@ -0,0 +1,107 @@
+"""
+
+Run pyLossless on a BIDS dataset.
+=================================
+
+In this notebook, we will run the pyLossless pipeline on a publicly available dataset
+that is already in BIDS format.
+"""
+
+# %%
+# Imports
+# -------
+from pathlib import Path
+import shutil
+import pylossless as ll
+
+# %%
+# Get the data
+# ------------
+raw, config, bids_path = ll.datasets.load_openneuro_bids()
+
+# %%
+# Prep the Raw object
+# -------------------
+#
+# This data has EOG channels that are not labeled as such. We will manually set the
+# channel types to be "eog" for these channels (i.e. "EXG1"). We will also crop the
+# data to 60 seconds for speed, and load the data in memory, which is required for
+# running the pipeline.
+raw.set_channel_types({ch: "eog" for ch in raw.ch_names if ch.startswith("EX")})
+raw.load_data().crop(0, 60)
+
+# %%
+# Initialize the pipeline
+# -----------------------
+#
+# The :class:`~pylossless.LosslessPipeline` instance is the main object that will
+# run the pipeline. It takes a file path to a :class:`~pylossless.config.Config` object
+# as input. :func:`~pylossless.datasets.load_openneuro_bids` returned a
+# :class:`~pylossless.config.Config` object, so we will save it to disk and pass the
+# file path to the :class:`~pylossless.LosslessPipeline` constructor.
+config_path = Path("lossless_config.yaml")
+config["filtering"]["notch_filter_args"]["freqs"] = [60]
+config.save(config_path)
+pipeline = ll.LosslessPipeline(config_path)
+
+# %%
+# Run the pipeline
+# ----------------
+#
+# The :class:`~pylossless.LosslessPipeline` object has a
+# :meth:`~pylossless.LosslessPipeline.run_with_raw` method that takes a
+# :class:`~mne.io.Raw` object as input.
+# We will use the :class:`~mne.io.Raw` object that was returned by
+# :func:`~pylossless.datasets.load_openneuro_bids` with the pipeline.
+pipeline.run_with_raw(raw)
+
+# %%
+# View the results
+# ----------------
+#
+# The :class:`~pylossless.LosslessPipeline` object stores flagged channels and ICs in
+# the :attr:`~pylossless.LosslessPipeline.flags` attribute:
+print(f"flagged channels: {pipeline.flags['ch']}")
+print(f"flagged ICs: {pipeline.flags['ic']}")
+
+# %%
+# Get the cleaned data
+# --------------------
+#
+# The :class:`~pylossless.LosslessPipeline` by default does not modify the
+# :class:`~mne.io.Raw` object that is passed to it, so none of the flagged channels
+# or ICs are removed from the :class:`~mne.io.Raw` object yet. To get the cleaned
+# :class:`~mne.io.Raw` object, we need to call the
+# :meth:`~pylossless.config.RejectionPolicy.apply` method. This method takes a
+# :class:`~pylossless.LosslessPipeline` as input, which specifies how to apply the flags
+# to generate a new :class:`~mne.io.Raw` object.
+rejection_policy = ll.RejectionPolicy()
+rejection_policy
+
+# %%
+# Note that we are using the default  channel cleaning mode, which is ``None``, meaning
+# that that the flagged channels will simply be added to ``raw.info["bads"]``. We also
+# could have specified ``"interpolate"``, which means that the flagged
+# channels would be interpolated using :func:`~mne.io.Raw.interpolate_bads`.
+cleaned_raw = rejection_policy.apply(pipeline)
+cleaned_raw.plot()
+
+# %%
+# Save the PyLossless Derivative
+# ------------------------------
+#
+# Let's save our pipeline output to disk. We need to use our
+# :class:`~mne_bids.BIDSPath` object to set up a derivatives path to save the
+# pipeline output to:
+derivatives_path = pipeline.get_derivative_path(bids_path)
+derivatives_path
+
+# %%
+pipeline.save(derivatives_path)
+
+# %%
+# Clean up
+# --------
+#
+shutil.rmtree(bids_path.root)
+config_path.unlink()
diff --git a/examples/usage.py b/examples/usage.py
new file mode 100644
index 0000000..75b9250
--- /dev/null
+++ b/examples/usage.py
@@ -0,0 +1,22 @@
+"""
+
+Open The QCR dashboard in a jupyter notebook.
+=============================================
+
+This notebook demonstrate how the PyLossless Quality Control interface can be
+integrated in a Jupyter Notebook. This can even be done on a virtualized
+platform, like Google Colab! To try it, you can run this notebook from
+`here!
+<https://colab.research.google.com/github/lina-usc/pylossless/blob/main/
+notebooks/qc_example.ipynb>`__.
+"""
+
+# %%
+# This tutorial assumes you have pylossless installed.
+# ----------------------------------------------------
+# Then, we import the function we need.
+# For use in jupyter notebooks, We just need to import a single function.
+from pylossless.dash.app import get_app
+
+app = get_app(kind="jupyter")
+app.run_server(mode="inline")
diff --git a/pylossless/__init__.py b/pylossless/__init__.py
index 16a49f6..873b977 100644
--- a/pylossless/__init__.py
+++ b/pylossless/__init__.py
@@ -1,4 +1,4 @@
-"""Python port of EEG-IP-L pipeline for preprocessing EEG."""
-
 from . import pipeline, bids, config, datasets, flagging, utils
+
 from .pipeline import LosslessPipeline
+from .config import Config, RejectionPolicy
diff --git a/pylossless/config/__init__.py b/pylossless/config/__init__.py
new file mode 100644
index 0000000..38c7aea
--- /dev/null
+++ b/pylossless/config/__init__.py
@@ -0,0 +1,2 @@
+from .config import Config
+from .rejection import RejectionPolicy
diff --git a/pylossless/config.py b/pylossless/config/config.py
similarity index 63%
rename from pylossless/config.py
rename to pylossless/config/config.py
index 3436bc5..d9f9e7f 100644
--- a/pylossless/config.py
+++ b/pylossless/config/config.py
@@ -4,8 +4,6 @@
 #
 # License: MIT
 
-"""Module file contains Config class which defines procedures."""
-
 import sys
 from pathlib import Path
 import yaml
@@ -13,13 +11,15 @@
 
 # TODO: Add parameters and references for below functions.
 # TODO: Add magic methods and constructors for this class.
-class Config(dict):
-    """Representation of configuration files for a pipeline procedure."""
+class ConfigMixin(dict):
+    """Base configuration file class for pipeline procedures."""
 
-    DEFAULT_CONFIG_PATH = Path(__file__).parent / "assets" / "ll_default_config.yaml"
+    DEFAULT_CONFIG_PATH = (
+        Path(__file__).parent.parent / "assets" / "ll_default_config.yaml"
+    )
 
     def read(self, file_name):
-        """Read a saved pylossless config file."""
+        """Read a saved pylossless config YAML file."""
         file_name = Path(file_name)
         if not file_name.exists():
             raise FileExistsError(
@@ -31,13 +31,14 @@ def read(self, file_name):
 
         return self
 
-    def load_default(self):
-        """Get the default pylossless config file."""
-        self.read(Config.DEFAULT_CONFIG_PATH)
-        return self
-
     def save(self, file_name):
-        """Save the current Config object to disk."""
+        """Save the current config object to disk as a YAML file.
+
+        Parameters
+        ----------
+        file_name : str | pathlib.Path
+            The file name to save the :class:`~pylossless.Config` object to.
+        """
         file_name = Path(file_name)
         with file_name.open("w") as init_variables_file:
             yaml.dump(dict(self), init_variables_file, indent=4, sort_keys=True)
@@ -45,3 +46,16 @@ def save(self, file_name):
     def print(self):
         """Print the Config contents."""
         yaml.dump(dict(self), sys.stdout, indent=4, sort_keys=True)
+
+
+class Config(ConfigMixin):
+    """Representation of configuration file for running the pipeline."""
+
+    DEFAULT_CONFIG_PATH = (
+        Path(__file__).parent.parent / "assets" / "ll_default_config.yaml"
+    )
+
+    def load_default(self):
+        """Get the default pylossless config file."""
+        self.read(Config.DEFAULT_CONFIG_PATH)
+        return self
diff --git a/pylossless/config/rejection.py b/pylossless/config/rejection.py
new file mode 100644
index 0000000..835bd9d
--- /dev/null
+++ b/pylossless/config/rejection.py
@@ -0,0 +1,156 @@
+# Authors: Christian O'Reilly <christian.oreilly@sc.edu>
+#          Scott Huberty <seh33@uw.edu>
+#
+# License: MIT
+
+import numpy as np
+
+from .config import ConfigMixin
+
+
+class RejectionPolicy(ConfigMixin):
+    """Class used to implement a rejection policy for a pipeline output.
+
+    Parameters
+    ----------
+    config_fname : pathlib.Path
+        path to config file specifying the parameters to be used
+        in for the rejection policy.
+    ch_flags_to_reject : list of str
+        List of channel flags to apply. "all" is the same as
+        ``["ch_sd", "low_r", "bridge"]``, meaning that any channels that
+        have been flagged ``"ch_sd"``, ``"low_r"``, or ``"bridge"`` will
+        be applied. Defaults to "all".
+    ic_flags_to_reject : list of str
+        List of IC flags to apply. "all" is the same as
+        ``["muscle", "heart", "eye", "channel noise", "line noise"]``, meaning that
+        any ICs that have been flagged ``"muscle"``, ``"ecg"``, ``"eog"``
+        ``"channel_noise"``, or ``"line_noise"`` will be applied, if their label
+        confidence is greater than ``ic_rejection_threshold``.
+        Note this list does NOT include ``"brain"`` and ``"other"``. Defaults to
+        "all".
+    ic_rejection_threshold : float
+        threshold (between 0 and 1), representing the minimum confidence
+        percentage in a given label (i.e., .30 is 30% confidence). For any labels
+        passed to ``ic_flags_to_reject`` (i.e. ``"ecg"``, ``"eog"``, etc.), ICs with
+        that label and a confidence percentage greater than this threshold will
+        be rejected. Defaults to ``0.3``.
+    ch_cleaning_mode : str
+        Must be one of ``None``, ``"drop"``, or ``"interpolate"``. ``None`` adds the
+        channels to ``raw.info["bads"]``. ``"drop"`` drops the channels from the
+        :class:`~mne.io.Raw` object. ``"interpolate"`` interpolates the bad channels.
+        Defaults to ``None``.
+    interpolate_bads_kwargs : None | dict
+        If ``ch_cleaning_mode`` is ``"interpolate"``, these keyword arguments
+        will be passed to ``raw.interpolate_bads(**interpolate_bads_kwargs)``.
+        Must be a dictionary of valid keyword arguments for
+        :meth:`~mne.io.Raw.interpolate_bads`. Defaults to ``None``, which means no
+        keyword arguments will be passed.
+    remove_flagged_ics : bool
+        If ``True``, subtracts the signal accounted for by the flagged ICs
+        from the ``raw`` object, via :meth:`~mne.preprocessing.ICA.apply`.
+        If ``False``, does nothing. Defaults to ``True``.
+
+    """
+
+    def __init__(
+        self,
+        *,
+        config_fname=None,
+        ch_flags_to_reject="all",
+        ic_flags_to_reject="all",
+        ic_rejection_threshold=0.3,
+        ch_cleaning_mode=None,
+        interpolate_bads_kwargs=None,
+        remove_flagged_ics=True,
+    ):
+        if ch_flags_to_reject == "all":
+            ch_flags_to_reject = ["ch_sd", "low_r", "bridge"]
+        if ic_flags_to_reject == "all":
+            ic_flags = ["muscle", "ecg", "eog", "channel_noise", "line_noise"]
+            ic_flags_to_reject = ic_flags
+
+        if interpolate_bads_kwargs is None:
+            interpolate_bads_kwargs = {}
+
+        if config_fname is not None:
+            config = ConfigMixin().read(config_fname)
+            for key, value in config.items():
+                if hasattr(self, key):
+                    setattr(self, key, value)
+
+        super().__init__(
+            config_fname=config_fname,
+            ch_flags_to_reject=ch_flags_to_reject,
+            ic_flags_to_reject=ic_flags_to_reject,
+            ic_rejection_threshold=ic_rejection_threshold,
+            ch_cleaning_mode=ch_cleaning_mode,
+            interpolate_bads_kwargs=interpolate_bads_kwargs,
+            remove_flagged_ics=remove_flagged_ics,
+        )
+
+    def __repr__(self):
+        """Return a summary of the Calibration object."""
+        return (
+            f"RejectionPolicy: |\n"
+            f"  config_fname: {self['config_fname']}\n"
+            f"  ch_flags_to_reject: {self['ch_flags_to_reject']}\n"
+            f"  ic_flags_to_reject: {self['ic_flags_to_reject']}\n"
+            f"  ic_rejection_threshold: {self['ic_rejection_threshold']}\n"
+            f"  ch_cleaning_mode: {self['ch_cleaning_mode']}\n"
+            f"  remove_flagged_ics: {self['remove_flagged_ics']}\n"
+        )
+
+    def apply(self, pipeline, return_ica=False):
+        """Return a cleaned new raw object based on the rejection policy.
+
+        Parameters
+        ----------
+        pipeline : LosslessPipeline
+            An instance of LosslessPipeline with after the pipeline has
+            been ran.
+        return_ica : bool
+            If ``True``, returns the :class:`~mne.preprocessing.ica` object used to
+            clean the :class:`~mne.io.Raw` object. Defaults to ``False``.
+
+        Returns
+        -------
+        mne.io.Raw
+            An :class:`~mne.io.Raw` instance with the appropriate channels and ICs
+            added to mne bads, interpolated, or dropped.
+        """
+        # Get the raw object
+        raw = pipeline.raw.copy()
+
+        # Add channels to be rejected as bads
+        for key in self["ch_flags_to_reject"]:
+            if key in pipeline.flags["ch"]:
+                raw.info["bads"] += pipeline.flags["ch"][key].tolist()
+
+        # Clean the channels
+        if self["ch_cleaning_mode"] == "drop":
+            raw.drop_channels(raw.info["bads"])
+        elif self["ch_cleaning_mode"] == "interpolate":
+            raw.interpolate_bads(**self["interpolate_bads_kwargs"])
+
+        # Clean the epochs
+        # TODO: Not sure where we landed on having these prefixed as bad_
+        #       or not by the pipeline. If not prefixed, this would be the
+        #       step that add select types of flags as bad_ annotations.
+
+        # Clean the ics
+        ic_labels = pipeline.flags["ic"]
+        mask = np.array([False] * len(ic_labels["confidence"]))
+        for label in self["ic_flags_to_reject"]:
+            mask |= ic_labels["ic_type"] == label
+        mask &= ic_labels["confidence"] > self["ic_rejection_threshold"]
+
+        flagged_ics = ic_labels.loc[mask]
+        if not flagged_ics.empty:
+            flagged_ics = flagged_ics.index.tolist()
+            pipeline.ica2.exclude.extend(flagged_ics)
+            pipeline.ica2.apply(raw)
+
+        if return_ica:
+            return raw, pipeline.ica2
+        return raw
diff --git a/pylossless/dash/__init__.py b/pylossless/dash/__init__.py
index f983d13..8e3537f 100644
--- a/pylossless/dash/__init__.py
+++ b/pylossless/dash/__init__.py
@@ -1,7 +1,11 @@
-"""Dash based helper functions for Lossless QC procedure."""
-
-from mne_icalabel.config import ICA_LABELS_TO_MNE
+from mne.utils import check_version
+import mne_icalabel.config as iclabel_config
 
+if not check_version("mne_icalabel", "0.5.0"):
+    raise ImportError(
+        "mne_icalabel version 0.5.0 or higher is required. "
+        "Please upgrade mne_icalabel."
+    )
 IC_COLORS = [
     "#2c2c2c",
     "#003e83",
@@ -16,4 +20,4 @@
     "plum",
 ]
 
-ic_label_cmap = dict(zip(ICA_LABELS_TO_MNE.values(), IC_COLORS))
+ic_label_cmap = dict(zip(iclabel_config.ICA_LABELS_TO_MNE.values(), IC_COLORS))
diff --git a/pylossless/datasets/__init__.py b/pylossless/datasets/__init__.py
index 8cb532a..f4efce1 100644
--- a/pylossless/datasets/__init__.py
+++ b/pylossless/datasets/__init__.py
@@ -1,1 +1,2 @@
+from .datasets import load_openneuro_bids
 from .simulated import load_simulated_raw
diff --git a/pylossless/datasets/datasets.py b/pylossless/datasets/datasets.py
new file mode 100644
index 0000000..9cfb519
--- /dev/null
+++ b/pylossless/datasets/datasets.py
@@ -0,0 +1,75 @@
+from time import sleep
+from pathlib import Path
+
+import mne_bids
+import openneuro
+
+import pylossless as ll
+
+
+def load_openneuro_bids(subject="pd6"):
+    """Download and Load BIDS dataset ds002778 from OpenNeuro.
+
+    Parameters
+    ----------
+    subject : str
+        The subject to load. Defaults to ``"pd6"``.
+
+    Returns
+    -------
+    raw : mne.io.Raw
+        The :class:`~mne.io.Raw` data.
+    config : pylossless.config.Config
+        The :class:`~pylossless.config.Config` object that can be used to run the
+        pipeline.
+    bids_path : mne_bids.BIDSPath
+        The :class:``~mne_bids.BIDSPath`` instance for the dataset.
+
+    Notes
+    -----
+    This function will download the dataset to the current working directory. The
+    dataset is from the "UC San Diego Resting State EEG Data from Patients with
+    Parkinson's Disease" study, curated by Alexander Rockhill.
+    """
+    config = ll.config.Config()
+    config.load_default()
+    config["project"]["bids_montage"] = ""
+    config["project"]["analysis_montage"] = "standard_1020"
+    config["project"]["set_montage_kwargs"]["on_missing"] = "warn"
+
+    # Shamelessly copied from
+    # https://mne.tools/mne-bids/stable/auto_examples/read_bids_datasets.html
+    # pip install openneuro-py
+
+    dataset = "ds002778"
+    subject = subject
+
+    # Download one subject's data from each dataset
+    bids_root = Path(".") / dataset
+    # TODO: Delete this directory after test otherwise MNE will think the
+    # sample directory is outdated, and will re-download it the next time
+    # data_path() is called, which is annoying for users.
+    bids_root.mkdir(exist_ok=True)
+
+    openneuro.download(
+        dataset=dataset, target_dir=bids_root, include=[f"sub-{subject}"]
+    )
+
+    datatype = "eeg"
+    session = "off"
+    task = "rest"
+    suffix = "eeg"
+    bids_path = mne_bids.BIDSPath(
+        subject=subject,
+        session=session,
+        task=task,
+        suffix=suffix,
+        datatype=datatype,
+        root=bids_root,
+    )
+
+    while not bids_path.fpath.with_suffix(".bdf").exists():
+        print(list(bids_path.fpath.glob("*")))
+        sleep(1)
+    raw = mne_bids.read_raw_bids(bids_path, verbose="ERROR")
+    return raw, config, bids_path


Test patch:
diff --git a/.github/workflows/test_pipeline.yml b/.github/workflows/test_pipeline.yml
index 43770ec..7a2c226 100644
--- a/.github/workflows/test_pipeline.yml
+++ b/.github/workflows/test_pipeline.yml
@@ -13,10 +13,10 @@ jobs:
     runs-on: macos-latest
     steps:
       - uses: actions/checkout@v3
-      - name: Set up Python 3.x
+      - name: Set up Python 3.9
         uses: actions/setup-python@v4
         with:
-          python-version: "3.x"
+          python-version: "3.9"
       - name: Install Pylossless & Deps
         run: pip install -e .
       - name: Install testing dependencies
diff --git a/docs/examples/test_config.yaml b/examples/test_config.yaml
similarity index 100%
rename from docs/examples/test_config.yaml
rename to examples/test_config.yaml
diff --git a/pylossless/conftest.py b/pylossless/conftest.py
index 37e3479..a6ae964 100644
--- a/pylossless/conftest.py
+++ b/pylossless/conftest.py
@@ -6,8 +6,10 @@
 from pathlib import Path
 import shutil
 
+from mne import Annotations
+
 import pylossless as ll
-from pylossless.tests.test_pipeline import load_openneuro_bids
+from pylossless.datasets import load_openneuro_bids
 
 import pytest
 
@@ -15,8 +17,13 @@
 @pytest.fixture(scope="session")
 def pipeline_fixture():
     """Return a namedTuple containing MNE eyetracking raw data and events."""
-    raw, config, bids_root = load_openneuro_bids()
+    raw, config, bids_path = load_openneuro_bids()
     raw.crop(tmin=0, tmax=60)  # take 60 seconds for speed
+    annots = Annotations(
+        onset=[1, 15], duration=[1, 1], description=["test_annot", "test_annot"]
+    )
+    raw.set_annotations(annots)
+
     config["find_breaks"] = {}
     config["find_breaks"]["min_break_duration"] = 9
     config["find_breaks"]["t_start_after_previous"] = 1
@@ -28,5 +35,5 @@ def pipeline_fixture():
     pipeline.run_with_raw(pipeline.raw)
 
     Path("test_config.yaml").unlink()  # delete config file
-    shutil.rmtree(bids_root)
+    shutil.rmtree(bids_path.root)
     return pipeline
diff --git a/pylossless/dash/tests/test_topo_viz.py b/pylossless/dash/tests/test_topo_viz.py
index a7f924e..3ada77c 100644
--- a/pylossless/dash/tests/test_topo_viz.py
+++ b/pylossless/dash/tests/test_topo_viz.py
@@ -68,9 +68,7 @@ def test_GridTopoPlot():
 
 # chromedriver: https://chromedriver.storage.googleapis.com/
 #               index.html?path=114.0.5735.90/
-@pytest.mark.skip(
-    reason="an issue with chromedriver causes failure, and test hangs for ~10 minutes"
-)
+@pytest.mark.xfail(reason="an issue with chromedriver causes failure. Need to debug.")
 def test_TopoViz(dash_duo):
     """Test TopoViz."""
     raw, ica = get_raw_ica()
diff --git a/pylossless/tests/test_pipeline.py b/pylossless/tests/test_pipeline.py
index af5b092..5d51cea 100644
--- a/pylossless/tests/test_pipeline.py
+++ b/pylossless/tests/test_pipeline.py
@@ -1,63 +1,10 @@
 from pathlib import Path
-from time import sleep
+
 import pytest
 
 import pylossless as ll
 
 import mne
-import mne_bids
-
-import openneuro
-
-
-def load_openneuro_bids():
-    """Load a BIDS dataset from OpenNeuro."""
-    config = ll.config.Config()
-    config.load_default()
-    config["project"]["bids_montage"] = ""
-    config["project"]["analysis_montage"] = "standard_1020"
-    config["project"]["set_montage_kwargs"]["on_missing"] = "warn"
-
-    # Shamelessly copied from
-    # https://mne.tools/mne-bids/stable/auto_examples/read_bids_datasets.html
-    # pip install openneuro-py
-
-    dataset = "ds002778"
-    subject = "pd6"
-
-    # Download one subject's data from each dataset
-    bids_root = Path(".") / dataset
-    # TODO: Delete this directory after test otherwise MNE will think the
-    # sample directory is outdated, and will re-download it the next time
-    # data_path() is called, which is annoying for users.
-    bids_root.mkdir(exist_ok=True)
-
-    openneuro.download(
-        dataset=dataset, target_dir=bids_root, include=[f"sub-{subject}"]
-    )
-
-    datatype = "eeg"
-    session = "off"
-    task = "rest"
-    suffix = "eeg"
-    bids_path = mne_bids.BIDSPath(
-        subject=subject,
-        session=session,
-        task=task,
-        suffix=suffix,
-        datatype=datatype,
-        root=bids_root,
-    )
-
-    while not bids_path.fpath.with_suffix(".bdf").exists():
-        print(list(bids_path.fpath.glob("*")))
-        sleep(1)
-    raw = mne_bids.read_raw_bids(bids_path, verbose="ERROR")
-    annots = mne.Annotations(
-        onset=[1, 15], duration=[1, 1], description=["test_annot", "test_annot"]
-    )
-    raw.set_annotations(annots)
-    return raw, config, bids_root
 
 
 def test_empty_repr(tmp_path):
diff --git a/pylossless/tests/test_rejection.py b/pylossless/tests/test_rejection.py
new file mode 100644
index 0000000..994340c
--- /dev/null
+++ b/pylossless/tests/test_rejection.py
@@ -0,0 +1,40 @@
+from pathlib import Path
+
+from numpy.testing import assert_array_equal
+
+import pytest
+
+import pylossless as ll
+
+
+@pytest.mark.parametrize("clean_ch_mode", [None, "drop", "interpolate"])
+def test_rejection_policy(clean_ch_mode, pipeline_fixture):
+    """Test the rejection policy."""
+    rejection_config_path = Path("test_rejection_config.yaml")
+    rejection_config = ll.RejectionPolicy(ch_cleaning_mode=clean_ch_mode)
+    rejection_config.save(rejection_config_path)
+    assert rejection_config["ch_flags_to_reject"] == ["ch_sd", "low_r", "bridge"]
+
+    raw, ica = rejection_config.apply(pipeline_fixture, return_ica=True)
+    flagged_chs = []
+    for key in rejection_config["ch_flags_to_reject"]:
+        flagged_chs.extend(pipeline_fixture.flags["ch"][key].tolist())
+    assert flagged_chs == ["PO3", "Oz", "O2"]
+    if clean_ch_mode is None:
+        assert_array_equal(flagged_chs, raw.info["bads"])
+    elif clean_ch_mode == "drop":
+        assert len(list(set(flagged_chs) - set(raw.ch_names))) == 3
+    elif clean_ch_mode == "interpolate":
+        # interpolate_bads will drop ch_names from raw.info["bads"]
+        # so to make sure that the channel was interpolated, lets check it
+        assert len(list(set(flagged_chs) - set(raw.info["bads"]))) == 3
+
+    df = pipeline_fixture.flags["ic"]
+    assert not df.loc[ica.exclude]["ic_type"].str.contains("brain").any()
+    assert not df.loc[ica.exclude]["ic_type"].str.contains("other").any()
+    assert df.loc[ica.exclude]["ic_type"].str.contains("muscle").any()
+    assert df.loc[ica.exclude]["ic_type"].str.contains("eog").any()
+    assert df.loc[ica.exclude]["ic_type"].str.contains("line_noise").any()
+    threshold = rejection_config["ic_rejection_threshold"]
+    assert (df.loc[ica.exclude]["confidence"] > threshold).all()
+    rejection_config_path.unlink()