Repository: sappelhoff/pyprep (python)

Issue:
PyPREP excludes interpolated NaN/flat/low-SNR channels from average referencing
While pulling together a fix for #77, I ended up catching another PyPREP/MatPREP difference I didn't notice before. Essentially, PyPREP re-uses a class attribute (`self.reference_channels`) to do the same job as two separate variables in MatPREP's `performReference` and `robustReference`, which ends up causing any channels flagged as bad-by-SNR, bad-by-NaN, or bad-by-flat on the initial NoisyChannels pass to be excluded from the average reference signal even after they've been interpolated.

Specifically, the `self.reference_channels` here:

https://github.com/sappelhoff/pyprep/blob/6abe0bcf1d1d40aefd2fe4cc41b200fbb0924ee4/pyprep/reference.py#L231-L238

gets reused here:

https://github.com/sappelhoff/pyprep/blob/6abe0bcf1d1d40aefd2fe4cc41b200fbb0924ee4/pyprep/reference.py#L114-L119

whereas in MATLAB PREP the former is a local variable and doesn't affect the latter.

Pull request:
Ensure all reference channels are used during final average referencing
<!--Thanks for contributing-->
<!--If this is your first time, please make sure to read the contributing guideline-->
<!--https://github.com/sappelhoff/pyprep/blob/master/.github/CONTRIBUTING.md-->

# PR Description

Closes #90. This PR fixes the issue by making `reference_channels` a local variable within `Reference.robust_reference()` so that it doesn't affect the `self.reference_channels` variable that's used to select reference channels during `Reference.perform_reference()`.

As an aside, is there any situation where it would make sense to only use a subset of all channels for final average reference calculation, post-interpolation? Wondering if the whole `reference_channels` parameter is a carry-over from MatPREP that only exists because EEGLAB doesn't have as nice a method of marking non-EEG channels as MNE. In MatPREP, channels are excluded from being sources for interpolation if they're not specified in the "reference channels" parameter, whereas PyPREP just uses all non-bad EEG channels, but I'm not sure if that would ever actually be a problem.

# Merge Checklist

<!--To merge your PR we need to first take the following points into account.-->
<!--Please just leave this checklist untouched-->

- [ ] the PR has been reviewed and all comments are resolved
- [ ] all [CI][what-is-ci] checks pass
- [ ] (if applicable): the PR description includes the phrase `closes #<issue-number>` to [automatically close an issue][auto-close-documentation]
- [ ] (if applicable): bug fixes, new features, or [API][what-is-api] changes are documented in [whats_new.rst][whats-new-file]


[what-is-ci]: https://help.github.com/en/actions/building-and-testing-code-with-continuous-integration/about-continuous-integration
[auto-close-documentation]: https://help.github.com/en/github/managing-your-work-on-github/linking-a-pull-request-to-an-issue#linking-a-pull-request-to-an-issue-using-a-keyword
[what-is-api]: https://en.wikipedia.org/wiki/Application_programming_interface
[whats-new-file]: https://github.com/sappelhoff/pyprep/blob/master/docs/whats_new.rst


Patch:
diff --git a/docs/whats_new.rst b/docs/whats_new.rst
index 82a4cda..ed0d7e2 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -59,6 +59,7 @@ Bug
 - Fixed RANSAC to avoid making unnecessary signal predictions for known-bad channels, matching MATLAB behaviour and reducing RAM requirements, by `Austin Hurst`_ (:gh:`72`)
 - Fixed a bug in :meth:`NoisyChannels.find_bad_by_correlation` that prevented it from being able to handle channels with dropouts (intermittent flat regions), by `Austin Hurst`_ (:gh:`81`).
 - Fixed :class:`~pyprep.NoisyChannels` so that it always runs "bad channel by NaN" and "bad channel by flat" detection, preventing these channels from causing problems with other :class:`~pyprep.NoisyChannels` methods, by `Austin Hurst`_ (:gh:`79`)
+- Fixed :class:`~pyprep.Reference` so that channels are no longer excluded from final average reference calculation if they were originally bad by NaN, flat, or low SNR, by `Austin Hurst`_ (:gh:`92`)
 
 API
 ~~~
diff --git a/pyprep/reference.py b/pyprep/reference.py
index 6b3619d..8aa4b04 100644
--- a/pyprep/reference.py
+++ b/pyprep/reference.py
@@ -4,7 +4,6 @@ import logging
 import numpy as np
 from mne.utils import check_random_state
 
-# from pyprep.noisy import Noisydata
 from pyprep.find_noisy_channels import NoisyChannels
 from pyprep.removeTrend import removeTrend
 from pyprep.utils import _set_diff, _union
@@ -233,9 +232,7 @@ class Reference:
             noisy_detector.bad_by_nan + noisy_detector.bad_by_flat,
             noisy_detector.bad_by_SNR,
         )
-        self.reference_channels = _set_diff(
-            self.reference_channels, self.unusable_channels
-        )
+        reference_channels = _set_diff(self.reference_channels, self.unusable_channels)
 
         # Initialize channels to permanently flag as bad during referencing
         self.noisy_channels = {
@@ -253,11 +250,9 @@ class Reference:
         # Get initial estimate of the reference by the specified method
         signal = raw.get_data() * 1e6
         self.reference_signal = (
-            np.nanmedian(raw.get_data(picks=self.reference_channels), axis=0) * 1e6
+            np.nanmedian(raw.get_data(picks=reference_channels), axis=0) * 1e6
         )
-        reference_index = [
-            self.ch_names_eeg.index(ch) for ch in self.reference_channels
-        ]
+        reference_index = [self.ch_names_eeg.index(ch) for ch in reference_channels]
         signal_tmp = self.remove_reference(
             signal, self.reference_signal, reference_index
         )
@@ -327,8 +322,7 @@ class Reference:
             else:
                 signal_tmp = signal
             self.reference_signal = (
-                np.nanmean(raw_tmp.get_data(picks=self.reference_channels), axis=0)
-                * 1e6
+                np.nanmean(raw_tmp.get_data(picks=reference_channels), axis=0) * 1e6
             )
 
             signal_tmp = self.remove_reference(


Test patch:
diff --git a/tests/test_reference.py b/tests/test_reference.py
index 5c87df1..18f0c71 100644
--- a/tests/test_reference.py
+++ b/tests/test_reference.py
@@ -26,6 +26,9 @@ def test_basic_input(raw, montage):
     assert type(reference.still_noisy_channels) == list
     assert type(reference.raw) == mne.io.edf.edf.RawEDF
 
+    # Make sure the set of reference channels weren't modified by re-referencing
+    assert params["ref_chs"] == reference.reference_channels
+
 
 @pytest.mark.usefixtures("raw", "montage")
 def test_all_bad_input(raw, montage):