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

Issue:
MNE-ICALABEL changes
There were a few breaking changes in MNE-ICALabel 0.5, which came out last week. Once is that they moved pytorch from a dependency to a weak dependency, meaning that pytorch is not installed automatically.  Users now have an option to install pytorch or install a package called onnyx.

I'm not sure what we should do here for the long term, but to get the CI's happy I will install pytorch directly in the CI workflows.

As an aside, we should add a CI workflow called something like `test latest` where we install the development version of all our hard dependencies, so that next time API breaking changes occur, we can catch them early. 




Discussion / hints:

Fix patch:
diff --git a/docs/examples/README.txt b/docs/examples/README.txt
index e727d9e..762332d 100644
--- a/docs/examples/README.txt
+++ b/docs/examples/README.txt
@@ -1,2 +1,2 @@
-Tutorials Gallery
-^^^^^^^^^^^^^^^^^
+PyLossless Tutorials
+^^^^^^^^^^^^^^^^^^^^
\ No newline at end of file
diff --git a/docs/examples/plot_0_implementation.py b/docs/examples/plot_0_implementation.py
new file mode 100644
index 0000000..93a8f8c
--- /dev/null
+++ b/docs/examples/plot_0_implementation.py
@@ -0,0 +1,191 @@
+r"""
+.. _notation:
+
+PyLossless Calculations: Understanding each step
+================================================
+
+This tutorial explains the calculations that PyLossless performs at each step of the
+pipeline. We will use example EEG data to demonstrate the
+calculations.
+
+Notation 
+--------
+Before we begin, we define some notation that will be used throughout the text:
+
+- ``s``, ``e``, and ``t`` are sensor, epochs, and time dimensions, respectively.
+
+- a 3D matrix :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 time points (i.e. samples).
+
+- We use superscripts to denote operations across a dimension, and we use subscripts to
+  denote indexing a dimension.
+
+- 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}`.
+
+- 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}`.
+
+- sensor-specific thresholds for rejecting epochs as
+  :math:`\tau^e_i \in \mathbb{R}^{S_\mathcal{G}}`
+
+- epoch-specific thresholds for rejecting sensors as
+  :math:`\tau^s_j \in \mathbb{R}^{E_\mathcal{G}}`
+
+- *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 default configuration file
+config = ll.config.Config()
+config.load_default()
+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
+
+# %%
+# 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.
+
+# %%
+# 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
+
+# %%
+# 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
+
+# %%
+# 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
+
+# %%
+# ## c) 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
+
+# %%
+# d) Identify noisy sensors part 2
+# --------------------------------
+# Next, we define a fractional 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}`.
+threshold = 0.2
+percent_outliers[percent_outliers > threshold]  # noisy sensors
diff --git a/docs/examples/usage.py b/docs/examples/usage.py
deleted file mode 100644
index b5276c9..0000000
--- a/docs/examples/usage.py
+++ /dev/null
@@ -1,22 +0,0 @@
-"""
-
-Run pyLossless using an example dataset.
-========================================
-
-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/docs/source/conf.py b/docs/source/conf.py
index 5d02672..ef1f24f 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -32,8 +32,9 @@
 extensions = [
     "sphinx.ext.intersphinx",
     "sphinx.ext.autodoc",
-    "numpydoc",
+    "sphinx.ext.mathjax",
     "sphinx.ext.todo",
+    "numpydoc",
     "sphinx_gallery.gen_gallery",
     "sphinxemoji.sphinxemoji",
     "sphinx_design",
@@ -62,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),
+    "pylossless": ("https://pylossless.readthedocs.io/en/latest/", None),
 }
 
 # -- Options for HTML output -------------------------------------------------
diff --git a/docs/source/implementation.rst b/docs/source/implementation.rst
deleted file mode 100644
index 8966cdd..0000000
--- a/docs/source/implementation.rst
+++ /dev/null
@@ -1,185 +0,0 @@
-Implementation
-==============
-
-The lossless pipeline be broken down into the following steps in Figure 1a.
-
-.. figure:: _images/pipeline_schema.png
-    
-    *Figure 1a)* Steps of the lossless pipeline
-
-
-Step 0: Making your data BIDS compliant
----------------------------------------
-
-The Lossless Pipeline expects the input data to comply with the BIDS standard.
-This assures that this pipeline knows how to load/save your data, and find important
-information (like channel names and locations), regardless
-of which manufacurer you used to collect your data.
-
-If you are not familiar with the BIDS standard for EEG data, see
-`the MNE-BIDS website <https://mne.tools/mne-bids/stable/index.html>`__ for
-information on how to make your data BIDS compliant using Python.
-
-Step 1: "Staging" Your data
----------------------------
-
-The first step of the pipeline is optional, but allows the user to carry out
-any operations on the data prior to running the rest of the pipeline.
-
-The most common staging operation will be identifying **break periods** that
-should *not* be considered by pipeline methods (because they are likely
-noisy, not representative of the data, and will not be used in analyses).
-
-Specifying parameters for finding breaks can be done directly in the YAML
-config file provided by pyLossless. For other staging operations, the user
-will have to provide their own script.
-
-.. todo::
-    TODO: create a tutorial for creating/using a custom staging script.
-
-Step 2: Identifying Outlier Channels
-------------------------------------
-.. figure:: _images/pipeline_step_2.png
-
-    *Figure 1b*) Step 2: Identifying Outlier Channels
-
-    .. todo::
-        TODO: explain this step
-
-Steps 3 & 4: Identifying Channels & Epochs with large standard deviations
--------------------------------------------------------------------------
-
-.. list-table::
-   :class: borderless
-
-   * - .. figure:: _images/pipeline_step_3.png
-            :align: center
-
-            *Figure 1c*) Step 3: Identify Channels with large standard deviations
-     - .. figure:: _images/pipeline_step_4.png
-            :align: center
-
-            *Figure 1d*) Step 4: Identify Epochs with large standard deviations across
-            channels
-
-.. todo::
-    TODO: explain this step
-
-Step 5: Filtering
------------------
-By default, the pipeline will filter the data with a 1Hz highpass. This is a
-requirement for running ICA, as low frequency oscillations are known to
-adversely affect ICA decomposition.
-
-Additional filtering parameters for lowpass and notch filtering can be set in
-the YAML config file provided by pyLossless.
-
-
-Step 6: Identifying nearest neighbor channels
----------------------------------------------
-
-.. figure:: _images/pipeline_step_6,9.png
-
-    *Figure 1e)*: Steps 6: Identifying nearest neighbor channels
-
-.. todo::
-    TODO: explain this step
-
-Step 7: Identifying bridged channels
-------------------------------------
-
-.. figure:: _images/pipeline_step_7.png
-
-    *Figure 1f)*: Step 7: Identifying bridged channels
-
-.. todo::
-    TODO: explain this step
-
-Step 8: Identifying the rank channel
-------------------------------------
-
-.. figure:: _images/pipeline_step_8.png
-
-    *Figure 1g)*: Step 8: Identifying the rank channel
-
-.. todo::
-    TODO: explain this step
-
-Step 9: Identifying nearest neighbor channels
----------------------------------------------
-
-.. figure:: _images/pipeline_step_6,9.png
-
-    *Figure 1h)*: Step 9: Identifying nearest neighbor channels
-
-.. todo::
-    TODO: explain this step
-
-Step 10: Identifying small time periods between flagged epochs
----------------------------------------------------------------
-
-.. figure:: _images/pipeline_step_10.png
-
-    *Figure 1i)*: Step 10: Identifying small time periods between flagged epochs
-
-.. todo::
-    TODO: explain this step
-
-Step 11: Run initial ICA
-------------------------
-
-The pyLossless pipeline will run an initial ICA that will only be used for step 11
-(identifying epochs with large standard deviations across ICs).
-
-By default, the fastICA algorithm is used for efficiency.
-
-Step 11: Identify Epochs with large standard deviations across ICs
-------------------------------------------------------------------
-
-.. figure:: _images/pipeline_step_12.png
-
-    *Figure 1j)*: Step 12: Epochs with large standard deviations across ICs
-
-.. todo::
-    TODO: explain this step
-
-Step 13: Run final ICA
-----------------------
-
-A final ICA will be run, excluding the epochs that were identified to contain
-large deviations across ICs. These time periods are excluded to improve the
-final ICA's decomposition.
-
-MNE-ICAlabel will be run on this decomposition to classify the ICs as brain,
-or a category of artifact.
-
-By default, Extended infomax is used for the final ICA.
-
-Step 14: Identifying small time periods between flagged epochs
----------------------------------------------------------------
-
-Again, the pipeline will identify small periods between any new epochs that
-have been flagged as artefactual in step 11.
-
-.. figure:: _images/pipeline_step_10.png
-
-    *Figure 1i)*: Step 14: Identifying small time periods between flagged epochs
-
-Step 15: Save data
-------------------
-The pipeline will save the data in a folder named ``derivatives/pyLossless``
-in the bids root directory.
-
-- By default, the pipeline will save EEG data in the European Data Format (``EDF``). 
-- The pipeline will save each ICA decomposition in  the ``.fif`` format.
-- Finally, the pipeline will save the IC labels in a ``.tsv`` file.
-
-Step 16: Expert Quality Control (QCR)
--------------------------------------
-
-The last step of the pipeline is for the user to review the pipelines decisions
-for each file. For this, the pyLossless pipeline will provide a dashboard to
-the user. 
-
-.. todo::
-    TODO: Finish the dashboard and provide a tutorial and demo
\ No newline at end of file
diff --git a/docs/source/index.rst b/docs/source/index.rst
index ff6a61e..89e8e5e 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -11,7 +11,6 @@ PyLossless EEG Processing Pipeline
    :hidden:
 
    install
-   implementation
    API/API_index
    generated/index
    contributing
diff --git a/pylossless/__init__.py b/pylossless/__init__.py
index c12cdbe..16a49f6 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, flagging, utils
+from . import pipeline, bids, config, datasets, flagging, utils
 from .pipeline import LosslessPipeline
diff --git a/pylossless/datasets/__init__.py b/pylossless/datasets/__init__.py
new file mode 100644
index 0000000..8cb532a
--- /dev/null
+++ b/pylossless/datasets/__init__.py
@@ -0,0 +1,1 @@
+from .simulated import load_simulated_raw
diff --git a/pylossless/datasets/simulated.py b/pylossless/datasets/simulated.py
new file mode 100644
index 0000000..e120058
--- /dev/null
+++ b/pylossless/datasets/simulated.py
@@ -0,0 +1,127 @@
+import numpy as np
+
+import mne
+
+from mne import make_ad_hoc_cov
+from mne.datasets import sample
+
+from mne.simulation import (
+    simulate_sparse_stc,
+    simulate_raw,
+    add_noise,
+    add_ecg,
+    add_eog,
+)
+
+
+def load_simulated_raw():
+    """Load simulated raw data.
+
+    This function returns a modified version of the MNE sample data.
+
+    Returns
+    -------
+    raw_ : instance of Raw
+        Simulated raw data based on the MNE sample data.
+
+    Notes
+    -----
+    This function is based on the code in
+    https://mne.tools/stable/auto_examples/simulation/simulate_raw_data.html
+    """
+    # LOAD DATA
+    data_path = sample.data_path()
+    meg_path = data_path / "MEG" / "sample"
+    raw_fname = meg_path / "sample_audvis_raw.fif"
+    fwd_fname = meg_path / "sample_audvis-meg-eeg-oct-6-fwd.fif"
+
+    # Load real data as the template
+    raw = mne.io.read_raw_fif(raw_fname)
+    raw.set_eeg_reference(projection=True)
+
+    # GENERATE DIPOLE TIME SERIES
+    n_dipoles = 4  # number of dipoles to create
+    epoch_duration = 2.0  # duration of each epoch/event
+    rng = np.random.RandomState(0)  # random state (make reproducible)
+    np.random.seed(5)
+
+    def data_fun(times, n=0):
+        """Generate time-staggered sinusoids at harmonics of 10Hz."""
+        n_samp = len(times)
+        window = np.zeros(n_samp)
+        start, stop = [
+            int(ii * float(n_samp) / (2 * n_dipoles)) for ii in (2 * n, 2 * n + 1)
+        ]
+        window[start:stop] = 1.0
+        n += 1
+        data = 25e-9 * np.sin(2.0 * np.pi * 10.0 * n * times)
+        data *= window
+        return data
+
+    times = raw.times[: int(raw.info["sfreq"] * epoch_duration)]
+    fwd = mne.read_forward_solution(fwd_fname)
+    src = fwd["src"]
+    stc = simulate_sparse_stc(
+        src, n_dipoles=n_dipoles, times=times, data_fun=data_fun, random_state=rng
+    )
+
+    # SIMULATE RAW DATA
+    raw_sim = simulate_raw(raw.info, [stc] * 10, forward=fwd, verbose=True)
+    cov = make_ad_hoc_cov(raw_sim.info)
+    add_noise(raw_sim, cov, iir_filter=[0.2, -0.2, 0.04], random_state=rng)
+    add_ecg(raw_sim, random_state=rng)
+    add_eog(raw_sim, random_state=rng)
+    raw_sim.pick("eeg")
+
+    # Save Info and Montage for later reuse
+    montage = raw_sim.get_montage()
+    info = mne.create_info(
+        ch_names=raw_sim.ch_names,
+        sfreq=raw_sim.info["sfreq"],
+        ch_types=raw_sim.get_channel_types(),
+    )
+
+    # MAKE A VERY NOISY TIME PERIOD
+    raw_selection1 = raw_sim.copy().crop(tmin=0, tmax=2, include_tmax=False)
+    raw_selection2 = raw_sim.copy().crop(tmin=2, tmax=3, include_tmax=False)
+    raw_selection3 = raw_sim.copy().crop(tmin=3, tmax=19.994505956825666)
+
+    cov_noisy_period = make_ad_hoc_cov(raw_selection2.info, std=dict(eeg=0.000002))
+    add_noise(
+        raw_selection2, cov_noisy_period, iir_filter=[0.2, -0.2, 0.04], random_state=rng
+    )
+
+    raw_selection1.append([raw_selection2, raw_selection3])
+    raw_selection1.set_annotations(None)
+    raw_sim = raw_selection1
+
+    # MAKE SOME VERY NOISY CHANNELS
+
+    make_these_noisy = ["EEG 001", "EEG 002"]
+    raw_noisy_chs = raw_sim.copy().pick(make_these_noisy)
+    cov_noisy = make_ad_hoc_cov(raw_noisy_chs.info, std=dict(eeg=0.000002))
+    add_noise(raw_noisy_chs, cov_noisy, iir_filter=[0.2, -0.2, 0.04], random_state=rng)
+
+    raw_sim.drop_channels(make_these_noisy)
+    raw_noisy_chs.add_channels([raw_sim], force_update_info=True)
+    raw_sim = raw_noisy_chs
+
+    # MAKE BRIDGED CHANNELS AND 1 FLAT CHANNEL
+    data = raw_sim.get_data()  # ch x times
+    data[52, :] = data[53, :]  # duplicate ch 53 and 54
+
+    # MAKE AN UNCORRELATED CH
+    min_val = data[28, :].min()
+    max_val = data[28, :].min() + 0.0000065
+    low_correlated_ch = np.random.uniform(
+        low=min_val, high=max_val, size=len(data[28, :])
+    )
+    data[28] = low_correlated_ch
+    # Shuffle it Again. in-place
+    # np.random.shuffle(data[23])
+
+    # Make new raw out of data
+    raw_sim = mne.io.RawArray(data, info)
+    # Re-set the montage
+    raw_sim.set_montage(montage)
+    return raw_sim


Test patch:
diff --git a/.github/workflows/test_pipeline.yml b/.github/workflows/test_pipeline.yml
index 6178ac1..43770ec 100644
--- a/.github/workflows/test_pipeline.yml
+++ b/.github/workflows/test_pipeline.yml
@@ -25,6 +25,8 @@ jobs:
         run: pip install -r requirements_qc.txt
       - name: install openneuro
         run: pip install openneuro-py
+      - name: install pytorch
+        run: pip install torch
       - name: Test Pipeline
         run: |
           coverage run -m pytest
diff --git a/pylossless/tests/test_simulated.py b/pylossless/tests/test_simulated.py
index efa54f5..d5ed959 100644
--- a/pylossless/tests/test_simulated.py
+++ b/pylossless/tests/test_simulated.py
@@ -2,123 +2,12 @@
 
 from pathlib import Path
 
-import numpy as np
 from numpy.testing import assert_array_equal
 
-import mne
-from mne import make_ad_hoc_cov
-from mne.datasets import sample
-
-from mne.simulation import (
-    simulate_sparse_stc,
-    simulate_raw,
-    add_noise,
-    add_ecg,
-    add_eog,
-)
 
 import pylossless as ll
 
-# LOAD DATA
-data_path = sample.data_path()
-meg_path = data_path / "MEG" / "sample"
-raw_fname = meg_path / "sample_audvis_raw.fif"
-fwd_fname = meg_path / "sample_audvis-meg-eeg-oct-6-fwd.fif"
-
-
-# Load real data as the template
-raw = mne.io.read_raw_fif(raw_fname)
-raw.set_eeg_reference(projection=True)
-
-
-# GENERATE DIPOLE TIME SERIES
-n_dipoles = 4  # number of dipoles to create
-epoch_duration = 2.0  # duration of each epoch/event
-n = 0  # harmonic number
-rng = np.random.RandomState(0)  # random state (make reproducible)
-np.random.seed(5)
-
-
-def data_fun(times):
-    """Generate time-staggered sinusoids at harmonics of 10Hz."""
-    global n
-    n_samp = len(times)
-    window = np.zeros(n_samp)
-    start, stop = [
-        int(ii * float(n_samp) / (2 * n_dipoles)) for ii in (2 * n, 2 * n + 1)
-    ]
-    window[start:stop] = 1.0
-    n += 1
-    data = 25e-9 * np.sin(2.0 * np.pi * 10.0 * n * times)
-    data *= window
-    return data
-
-
-times = raw.times[: int(raw.info["sfreq"] * epoch_duration)]
-fwd = mne.read_forward_solution(fwd_fname)
-src = fwd["src"]
-stc = simulate_sparse_stc(
-    src, n_dipoles=n_dipoles, times=times, data_fun=data_fun, random_state=rng
-)
-
-# SIMULATE RAW DATA
-raw_sim = simulate_raw(raw.info, [stc] * 10, forward=fwd, verbose=True)
-cov = make_ad_hoc_cov(raw_sim.info)
-add_noise(raw_sim, cov, iir_filter=[0.2, -0.2, 0.04], random_state=rng)
-add_ecg(raw_sim, random_state=rng)
-add_eog(raw_sim, random_state=rng)
-raw_sim.pick("eeg")
-
-# Save Info and Montage for later reuse
-montage = raw_sim.get_montage()
-info = mne.create_info(
-    ch_names=raw_sim.ch_names,
-    sfreq=raw_sim.info["sfreq"],
-    ch_types=raw_sim.get_channel_types(),
-)
-
-
-# MAKE A VERY NOISY TIME PERIOD
-raw_selection1 = raw_sim.copy().crop(tmin=0, tmax=2, include_tmax=False)
-raw_selection2 = raw_sim.copy().crop(tmin=2, tmax=3, include_tmax=False)
-raw_selection3 = raw_sim.copy().crop(tmin=3, tmax=19.994505956825666)
-
-cov_noisy_period = make_ad_hoc_cov(raw_selection2.info, std=dict(eeg=0.000002))
-add_noise(
-    raw_selection2, cov_noisy_period, iir_filter=[0.2, -0.2, 0.04], random_state=rng
-)
-
-raw_selection1.append([raw_selection2, raw_selection3])
-raw_selection1.set_annotations(None)
-raw_sim = raw_selection1
-
-# MAKE SOME VERY NOISY CHANNELS
-
-make_these_noisy = ["EEG 001", "EEG 002"]
-raw_noisy_chs = raw_sim.copy().pick(make_these_noisy)
-cov_noisy = make_ad_hoc_cov(raw_noisy_chs.info, std=dict(eeg=0.000002))
-add_noise(raw_noisy_chs, cov_noisy, iir_filter=[0.2, -0.2, 0.04], random_state=rng)
-
-raw_sim.drop_channels(make_these_noisy)
-raw_noisy_chs.add_channels([raw_sim], force_update_info=True)
-raw_sim = raw_noisy_chs
-
-# MAKE BRIDGED CHANNELS AND 1 FLAT CHANNEL
-data = raw_sim.get_data()  # ch x times
-data[52, :] = data[53, :]  # duplicate ch 53 and 54
-
-# MAKE AN UNCORRELATED CH
-min_val = data[28, :].min()
-max_val = data[28, :].min() + 0.0000065
-low_correlated_ch = np.random.uniform(low=min_val, high=max_val, size=len(data[28, :]))
-data[28] = low_correlated_ch
-# Shuffle it Again. in-place
-# np.random.shuffle(data[23])
-
-# Make new raw out of data
-raw_sim = mne.io.RawArray(data, info)
-# Re-set the montage
-raw_sim.set_montage(montage)
+raw_sim = ll.datasets.load_simulated_raw()
 
 # LOAD DEFAULT CONFIG
 config = ll.config.Config()