ENH: Streamline ICA reporting (#899)

Co-authored-by: Richard Höchenberger <[email protected]>

.circleci/config.yml

@@ -50,7 +50,7 @@ jobs:
             pip install --upgrade --progress-bar off pip
             pip install --upgrade --progress-bar off "autoreject @ https://api.github.com/repos/autoreject/autoreject/zipball/master" "mne[hdf5] @ git+https://github.com/mne-tools/mne-python@main" "mne-bids[full] @ https://api.github.com/repos/mne-tools/mne-bids/zipball/main" numba
             pip install -ve .[tests]
-            pip install "PyQt6!=6.6.1,!=6.6.2" "PyQt6-Qt6!=6.6.1,!=6.6.2"
+            pip install "PyQt6!=6.6.1" "PyQt6-Qt6!=6.6.1,!=6.6.2,!=6.6.3"
       - run:
           name: Check Qt
           command: |

.circleci/remove_examples.sh

@@ -0,0 +1,33 @@
+#!/bin/bash
+
+set -eo pipefail
+
+VER=$1
+if [ -z "$VER" ]; then
+  echo "Usage: $0 <version>"
+  exit 1
+fi
+ROOT="$PWD/$VER/examples/"
+if [ ! -d ${ROOT} ]; then
+  echo "Version directory does not exist or appears incorrect:"
+  echo
+  echo "$ROOT"
+  echo
+  echo "Are you on the gh-pages branch and is the ds000117 directory present?"
+  exit 1
+fi
+if [ ! -d ${ROOT}ds000117 ]; then
+  echo "Directory does not exist:"
+  echo
+  echo "${ROOT}ds000117"
+  echo
+  echo "Assuming already pruned and exiting."
+  exit 0
+fi
+echo "Pruning examples in ${ROOT} ..."
+
+find $ROOT -type d -name "*" | tail -n +2 | xargs rm -Rf
+find $ROOT -name "*.html" -exec sed -i /^\<h2\ id=\"generated-output\"\>Generated/,/^\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $/{d} {} \;
+find $ROOT -name "*.html" -exec sed -i '/^  <a href="#generated-output"/,/^  <\/a>$/{d}' {} \;
+
+echo "Done"

.circleci/setup_bash.sh

@@ -55,6 +55,7 @@ echo 'export MPLBACKEND=Agg' >> "$BASH_ENV"
 echo "source ~/python_env/bin/activate" >> "$BASH_ENV"
 echo "export MNE_3D_OPTION_MULTI_SAMPLES=1" >> "$BASH_ENV"
 echo "export MNE_BIDS_PIPELINE_FORCE_TERMINAL=true" >> "$BASH_ENV"
+echo "export FORCE_COLOR=1" >> "$BASH_ENV"  # for rich to use color in logs
 mkdir -p ~/.local/bin
 if [[ ! -f ~/.local/bin/python ]]; then
     ln -s ~/python_env/bin/python ~/.local/bin/python

docs/source/examples/gen_examples.py

@@ -109,6 +109,18 @@ def _gen_demonstrated_funcs(example_config_path: Path) -> dict:
         datasets_without_html.append(dataset_name)
         continue
 
+    # For ERP_CORE, cut down on what we show otherwise our website is huge
+    if "ERP_CORE" in test_name:
+        show = ["015", "average"]
+        orig_fnames = html_report_fnames
+        html_report_fnames = [
+            f
+            for f in html_report_fnames
+            if any(f"sub-{s}" in f.parts and f"sub-{s}" in f.name for s in show)
+        ]
+        assert len(html_report_fnames), orig_fnames
+        del orig_fnames
+
     fname_iter = tqdm(
         html_report_fnames,
         desc=f"  {test_name}",

docs/source/v1.9.md.inc

@@ -8,8 +8,10 @@
 
 ### :warning: Behavior changes
 
-- Changed default for `source_info_path_update` to `None`.  In `_04_make_forward.py` 
-  and `_05_make_inverse.py`, we retrieve the info from the file from which 
+- All ICA HTML reports have been consolidated in the standard subject `*_report.html`
+  file instead of producing separate files (#899 by @larsoner).
+- Changed default for `source_info_path_update` to `None`.  In `_04_make_forward.py`
+  and `_05_make_inverse.py`, we retrieve the info from the file from which
   the `noise_cov` is computed (#919 by @SophieHerbst)
 - The [`depth`][mne_bids_pipeline._config.depth] parameter doesn't accept `None`
   anymore. Please use `0` instead. (#915 by @hoechenberger)
@@ -31,6 +33,7 @@
   bad channels from the first pipeline run. Now, we ensure that the original bad channels would be used and the
   related section is removed from the report in this case. (#902 by @larsoner)
 - Fixed group-average decoding statistics were not updated in some cases, even if relevant configuration options had been changed. (#902 by @larsoner)
+- Fixed a compatibility bug with joblib 1.4.0
 
 ### :medical_symbol: Code health and infrastructure
 

mne_bids_pipeline/_config.py

@@ -1126,6 +1126,16 @@
 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.
+
+!!! warning "ICA requires manual intervention!"
+    After the automatic ICA component detection step, review each subject's
+    `*_report.html` report file check if the set of ICA components to be removed
+    is correct. Adjustments should be made to the `*_proc-ica_components.tsv`
+    file, which will then be used in the step that is applied during ICA.
+
+    ICA component order can be considered arbitrary, so any time the ICA is
+    re-fit – i.e., if you change any parameters that affect steps prior to
+    ICA fitting – this file will need to be updated!
 """
 
 min_ecg_epochs: Annotated[int, Ge(1)] = 5

mne_bids_pipeline/_run.py

@@ -277,7 +277,12 @@ def wrapper(*args, **kwargs):
 
             # https://joblib.readthedocs.io/en/latest/memory.html#joblib.memory.MemorizedFunc.call  # noqa: E501
             if force_run or unknown_inputs or bad_out_files:
-                out_files, _ = memorized_func.call(*args, **kwargs)
+                # Joblib 1.4.0 only returns the output, but 1.3.2 returns both.
+                # Fortunately we can use tuple-ness to tell the difference (we always
+                # return None or a dict)
+                out_files = memorized_func.call(*args, **kwargs)
+                if isinstance(out_files, tuple):
+                    out_files = out_files[0]
             else:
                 out_files = memorized_func(*args, **kwargs)
             if self.require_output:

mne_bids_pipeline/steps/preprocessing/_06a1_fit_ica.py

@@ -81,9 +81,6 @@ def run_ica(
     out_files["epochs"] = (
         out_files["ica"].copy().update(suffix="epo", processing="icafit")
     )
-    out_files["report"] = bids_basename.copy().update(
-        processing="icafit", suffix="report", extension=".h5"
-    )
     del bids_basename
 
     # Generate a list of raw data paths (i.e., paths of individual runs)
@@ -247,17 +244,23 @@ def run_ica(
     logger.info(**gen_log_kwargs(message=msg))
     ica.save(out_files["ica"], overwrite=True)
 
-    # Start a report
+    # Add to report
+    tags = ("ica", "epochs")
+    title = "ICA: epochs for fitting"
     with _open_report(
         cfg=cfg,
         exec_params=exec_params,
         subject=subject,
         session=session,
         task=cfg.task,
-        fname_report=out_files["report"],
-        name="ICA.fit report",
     ) as report:
-        report.title = f"ICA – {report.title}"
+        report.add_epochs(
+            epochs=epochs,
+            title=title,
+            drop_log_ignore=(),
+            replace=True,
+            tags=tags,
+        )
         if cfg.ica_reject == "autoreject_local":
             caption = (
                 f"Autoreject was run to produce cleaner epochs before fitting ICA. "
@@ -274,19 +277,14 @@ def run_ica(
             )
             report.add_figure(
                 fig=fig,
-                title="Epochs: Autoreject cleaning",
+                title="Autoreject cleaning",
+                section=title,
                 caption=caption,
-                tags=("ica", "epochs", "autoreject"),
+                tags=tags + ("autoreject",),
                 replace=True,
             )
             plt.close(fig)
             del caption
-        report.add_epochs(
-            epochs=epochs,
-            title="Epochs used for ICA fitting",
-            drop_log_ignore=(),
-            replace=True,
-        )
     return _prep_out_files(exec_params=exec_params, out_files=out_files)
 
 

mne_bids_pipeline/steps/preprocessing/_06a2_find_ica_artifacts.py

@@ -7,7 +7,6 @@
 run the apply_ica step.
 """
 
-import shutil
 from types import SimpleNamespace
 from typing import Literal
 
@@ -113,9 +112,6 @@ def get_input_fnames_find_ica_artifacts(
         )
         _update_for_splits(in_files, key, single=True)
     in_files["ica"] = bids_basename.copy().update(processing="icafit", suffix="ica")
-    in_files["report"] = bids_basename.copy().update(
-        processing="icafit", suffix="report", extension=".h5"
-    )
     return in_files
 
 
@@ -143,9 +139,6 @@ def find_ica_artifacts(
     out_files_components = bids_basename.copy().update(
         processing="ica", suffix="components", extension=".tsv"
     )
-    out_files["report"] = bids_basename.copy().update(
-        processing="ica+components", suffix="report", extension=".h5"
-    )
     del bids_basename
     msg = "Loading ICA solution"
     logger.info(**gen_log_kwargs(message=msg))
@@ -297,41 +290,31 @@ def find_ica_artifacts(
     ecg_scores = None if len(ecg_scores) == 0 else ecg_scores
     eog_scores = None if len(eog_scores) == 0 else eog_scores
 
-    shutil.copyfile(in_files.pop("report"), out_files["report"])
+    title = "ICA: components"
     with _open_report(
         cfg=cfg,
         exec_params=exec_params,
         subject=subject,
         session=session,
         task=cfg.task,
-        fname_report=out_files["report"],
-        name="ICA report",
     ) as report:
+        logger.info(**gen_log_kwargs(message=f'Adding "{title}" to report.'))
         report.add_ica(
             ica=ica,
-            title="ICA cleaning",
+            title=title,
             inst=epochs,
             ecg_evoked=ecg_evoked,
             eog_evoked=eog_evoked,
             ecg_scores=ecg_scores,
             eog_scores=eog_scores,
             replace=True,
             n_jobs=1,  # avoid automatic parallelization
+            tags=("ica",),  # the default but be explicit
         )
 
-    msg = (
-        f"ICA completed. Please carefully review the extracted ICs in the "
-        f"report {out_files['report'].basename}, and mark all components "
-        f"you wish to reject as 'bad' in "
-        f"{out_files_components.basename}"
-    )
-    logger.info(**gen_log_kwargs(message=msg))
-
-    report.save(
-        out_files["report"],
-        overwrite=True,
-        open_browser=exec_params.interactive,
-    )
+    msg = 'Carefully review the extracted ICs and mark components "bad" in:'
+    logger.info(**gen_log_kwargs(message=msg, emoji="🛑"))
+    logger.info(**gen_log_kwargs(message=str(out_files_components), emoji="🛑"))
 
     assert len(in_files) == 0, in_files.keys()
     return _prep_out_files(exec_params=exec_params, out_files=out_files)

mne_bids_pipeline/steps/preprocessing/_08a_apply_ica.py

@@ -10,7 +10,6 @@
 import mne
 import pandas as pd
 from mne.preprocessing import read_ica
-from mne.report import Report
 from mne_bids import BIDSPath
 
 from ..._config_utils import (
@@ -21,7 +20,7 @@
 from ..._import_data import _get_run_rest_noise_path, _import_data_kwargs
 from ..._logging import gen_log_kwargs, logger
 from ..._parallel import get_parallel_backend, parallel_func
-from ..._report import _add_raw, _agg_backend, _open_report
+from ..._report import _add_raw, _open_report
 from ..._run import _prep_out_files, _update_for_splits, failsafe_run, save_logs
 
 
@@ -116,12 +115,8 @@ def apply_ica_epochs(
     session: str | None,
     in_files: dict,
 ) -> dict:
-    bids_basename = in_files["ica"].copy().update(processing=None)
     out_files = dict()
     out_files["epochs"] = in_files["epochs"].copy().update(processing="ica", split=None)
-    out_files["report"] = bids_basename.copy().update(
-        processing="ica", suffix="report", extension=".html"
-    )
 
     title = f"ICA artifact removal – sub-{subject}"
     if session is not None:
@@ -156,32 +151,6 @@ def apply_ica_epochs(
         split_size=cfg._epochs_split_size,
     )
     _update_for_splits(out_files, "epochs")
-
-    # Compare ERP/ERF before and after ICA artifact rejection. The evoked
-    # response is calculated across ALL epochs, just like ICA was run on
-    # all epochs, regardless of their respective experimental condition.
-    #
-    # We apply baseline correction here to (hopefully!) make the effects of
-    # ICA easier to see. Otherwise, individual channels might just have
-    # arbitrary DC shifts, and we wouldn't be able to easily decipher what's
-    # going on!
-    report = Report(out_files["report"], title=title, verbose=False)
-    picks = ica.exclude if ica.exclude else None
-    with _agg_backend():
-        report.add_ica(
-            ica=ica,
-            title="Effects of ICA cleaning",
-            inst=epochs.copy().apply_baseline(cfg.baseline),
-            picks=picks,
-            replace=True,
-            n_jobs=1,  # avoid automatic parallelization
-        )
-    report.save(
-        out_files["report"],
-        overwrite=True,
-        open_browser=exec_params.interactive,
-    )
-
     assert len(in_files) == 0, in_files.keys()
 
     # Report
@@ -198,18 +167,18 @@ def apply_ica_epochs(
             exec_params=exec_params,
             subject=subject,
             session=session,
-            name="ICA.apply report",
         ) as report:
             report.add_ica(
                 ica=ica,
-                title="ICA",
+                title="ICA: removals",
                 inst=epochs,
                 picks=ica.exclude,
                 # TODO upstream
                 # captions=f'Evoked response (across all epochs) '
                 # f'before and after ICA '
                 # f'({len(ica.exclude)} ICs removed)'
                 replace=True,
+                n_jobs=1,  # avoid automatic parallelization
             )
 
     return _prep_out_files(exec_params=exec_params, out_files=out_files)

mne_bids_pipeline/tests/configs/config_ERP_CORE.py

@@ -1,6 +1,6 @@
 """ERP CORE.
 
-This example demonstrate how to process 5 participants from the
+This example demonstrates how to process 5 participants from the
 [ERP CORE](https://erpinfo.org/erp-core) dataset. It shows how to obtain 7 ERP
 components from a total of 6 experimental tasks:
 

mne_bids_pipeline/tests/conftest.py

@@ -66,6 +66,10 @@ def pytest_configure(config):
     ignore:Persisting input arguments took.*:UserWarning
     # matplotlib needs to update
     ignore:Conversion of an array with ndim.*:DeprecationWarning
+    # scipy
+    ignore:nperseg .* is greater.*:UserWarning
+    # NumPy 2.0
+    ignore:__array_wrap__ must accept context.*:DeprecationWarning
     """
     for warning_line in warning_lines.split("\n"):
         warning_line = warning_line.strip()