Repository: sappelhoff/pyprep (python)

Issue:
add channels in raw.info['bads'] to unusable channels in Reference class
Currently, PREP assumes you have no prior information about which channels are bad. This isn't necessarily a bad design choice if it's philosophically important that all bad channel identification be automatic, but currently it will result in an error if you have any channels marked as bad in `raw.info['bads']` before running PREP.

The culprit re: the error is in `Reference.__init__`:
```
self.raw.pick_types(eeg=True, eog=False, meg=False)
```
`Raw.pick_types` will not include bad channels in its output by default, and since `Reference` inherits `self.prep_params["ref_chs"] = self.ch_names_eeg` from the `PrepPipeline` class that initializes it (where channels manually marked as bad _haven't_ been removed yet), then the `self.raw.get_data(picks=self.reference_channels)` in `Reference.robust_reference` will throw a missing value error.

So a quick solution would be to change the offending line to 
```
self.raw.pick_types(eeg=True, eog=False, meg=False, exclude=[])
```
so that manually marked bad channels are not unintentionally removed. But they should probably also be added to `Reference.unusable_channels`, so users have the ability to manually tell PREP to leave out certain channels (e.g. if they know the channel is broken in their hardware). 

In the meantime, users can circumvent the problem by just setting `raw.info['bads'] = []` before running PREP. Hopefully that saves someone some trouble :)

Pull request description:
Reference interpolates instead of deletes channels manually marked as bad, closes #146
<!--Thanks for contributing-->
<!--If this is your first time, please make sure to read the contributing guideline-->
<!--https://github.com/sappelhoff/pyprep/blob/main/.github/CONTRIBUTING.md-->

# PR Description

 <!--Please describe your pull request here.-->

Fixes #146, an error that results from the `raw.pick_types` call inside `pyprep.Reference` deleting channels which have been marked as bad prior to PREP. These channels should now be interpolated instead.

A new field `'bad_manual'` is added to `noisy_channels_before_interpolation` dict so there's still a record the channel was interpolated. This happens inside `Reference`; the `NoisyChannels` class is untouched. 


# 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): the changes are documented in the changelog [changelog.rst][changelog]
- [ ] (if applicable): new contributors have added themselves to the authors list in the [`CITATION.cff`][citationcff] 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
[changelog]: https://github.com/sappelhoff/pyprep/blob/main/docs/changelog.rst
[citationcff]: https://github.com/sappelhoff/pyprep/blob/main/CITATION.cff


Patch:
diff --git a/CITATION.cff b/CITATION.cff
index dce4ee7..1da3c61 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -51,6 +51,10 @@ authors:
     - given-names: Ayush
       family-names: Agarwal
       affiliation: 'Techno India University, India'
+    - given-names: John
+      family-names: Veillette
+      affiliation: 'Department of Psychology, University of Chicago, Chicago, IL, USA'
+      orcid: 'https://orcid.org/0000-0002-0332-4372'
 type: software
 repository-code: 'https://github.com/sappelhoff/pyprep'
 license: MIT
diff --git a/docs/authors.rst b/docs/authors.rst
index 13779d0..0ddb8f6 100644
--- a/docs/authors.rst
+++ b/docs/authors.rst
@@ -10,3 +10,4 @@
 .. _Stefan Appelhoff: https://stefanappelhoff.com/
 .. _Victor Xiang: https://github.com/Nick3151
 .. _Yorguin Mantilla: https://github.com/yjmantilla
+.. _John Veillette: https://github.com/john-veillette
diff --git a/docs/changelog.rst b/docs/changelog.rst
index 331724c..d389553 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -31,7 +31,7 @@ Changelog
 
 Bug
 ~~~
-- nothing yet
+- :class:`~pyprep.Reference` now keeps and interpolates channels channels manually marked as bad before PREP, by `John Veillette`_ (:gh:`146`)
 
 Code health
 ~~~~~~~~~~~
diff --git a/pyprep/find_noisy_channels.py b/pyprep/find_noisy_channels.py
index 65431d9..44b1a34 100644
--- a/pyprep/find_noisy_channels.py
+++ b/pyprep/find_noisy_channels.py
@@ -59,6 +59,7 @@ class NoisyChannels:
 
         raw.load_data()
         self.raw_mne = raw.copy()
+        self.bad_by_manual = raw.info["bads"]
         self.raw_mne.pick_types(eeg=True)
         self.sample_rate = raw.info["sfreq"]
         if do_detrend:
@@ -162,6 +163,7 @@ class NoisyChannels:
             "bad_by_SNR": self.bad_by_SNR,
             "bad_by_dropout": self.bad_by_dropout,
             "bad_by_ransac": self.bad_by_ransac,
+            "bad_by_manual": self.bad_by_manual,
         }
 
         all_bads = set()
diff --git a/pyprep/reference.py b/pyprep/reference.py
index 8fea7cc..c754673 100644
--- a/pyprep/reference.py
+++ b/pyprep/reference.py
@@ -83,7 +83,7 @@ class Reference:
         raw.load_data()
         self.raw = raw.copy()
         self.ch_names = self.raw.ch_names
-        self.raw.pick_types(eeg=True, eog=False, meg=False)
+        self.raw.pick_types(eeg=True, eog=False, meg=False, exclude=[])
         self.ch_names_eeg = self.raw.ch_names
         self.EEG = self.raw.get_data()
         self.reference_channels = params["ref_chs"]
@@ -97,6 +97,7 @@ class Reference:
         self.random_state = check_random_state(random_state)
         self._extra_info = {}
         self.matlab_strict = matlab_strict
+        self.bads_manual = raw.info["bads"]
 
     def perform_reference(self, max_iterations=4):
         """Estimate the true signal mean and interpolate bad channels.
@@ -149,6 +150,7 @@ class Reference:
         self.bad_before_interpolation = noisy_detector.get_bads(verbose=True)
         self.EEG_before_interpolation = self.EEG.copy()
         self.noisy_channels_before_interpolation = noisy_detector.get_bads(as_dict=True)
+        self.noisy_channels_before_interpolation["bad_by_manual"] = self.bads_manual
         self._extra_info["interpolated"] = noisy_detector._extra_info
 
         bad_channels = _union(self.bad_before_interpolation, self.unusable_channels)
@@ -223,7 +225,7 @@ class Reference:
         # Determine channels to use/exclude from initial reference estimation
         self.unusable_channels = _union(
             noisy_detector.bad_by_nan + noisy_detector.bad_by_flat,
-            noisy_detector.bad_by_SNR,
+            noisy_detector.bad_by_SNR + self.bads_manual,
         )
         reference_channels = _set_diff(self.reference_channels, self.unusable_channels)
 
@@ -237,6 +239,7 @@ class Reference:
             "bad_by_SNR": [],
             "bad_by_dropout": [],
             "bad_by_ransac": [],
+            "bad_by_manual": self.bads_manual,
             "bad_all": [],
         }
 
@@ -282,7 +285,8 @@ class Reference:
                 noisy[bad_type] = _union(noisy[bad_type], noisy_new[bad_type])
                 if bad_type not in ignore:
                     bad_chans.update(noisy[bad_type])
-            noisy["bad_all"] = list(bad_chans)
+            noisy["bad_by_manual"] = self.bads_manual
+            noisy["bad_all"] = list(bad_chans) + self.bads_manual
             logger.info(f"Bad channels: {noisy}")
 
             if (


Test patch:
diff --git a/tests/test_find_noisy_channels.py b/tests/test_find_noisy_channels.py
index be2e6c1..594fbc5 100644
--- a/tests/test_find_noisy_channels.py
+++ b/tests/test_find_noisy_channels.py
@@ -74,6 +74,20 @@ def test_bad_by_nan(raw_tmp):
     assert nd.bad_by_nan == [raw_tmp.ch_names[nan_idx]]
 
 
+def test_bad_by_manual(raw_tmp):
+    """Test the detection of channels marked bad a priori."""
+    n_chans = raw_tmp.get_data().shape[0]
+    nan_idx = int(rng.integers(0, n_chans, 1)[0])
+    raw_tmp._data[nan_idx, 3] = np.nan
+    raw_tmp.info["bads"] = [raw_tmp.ch_names[0]]
+
+    # Test record of a priori bad channels on NoisyChannels init
+    nd = NoisyChannels(raw_tmp, do_detrend=False)
+    assert nd.bad_by_manual == [raw_tmp.ch_names[0]]
+    assert raw_tmp.ch_names[0] in nd.get_bads(as_dict=False)
+    raw_tmp.info["bads"] = []
+
+
 def test_bad_by_flat(raw_tmp):
     """Test the detection of channels with flat or very weak signals."""
     # Make the signal for a random channel extremely weak
diff --git a/tests/test_reference.py b/tests/test_reference.py
index 5f01ed1..44ab3d7 100644
--- a/tests/test_reference.py
+++ b/tests/test_reference.py
@@ -23,6 +23,7 @@ def test_basic_input(raw, montage):
     reference = Reference(raw_tmp, params, ransac=False)
     reference.perform_reference()
     assert isinstance(reference.noisy_channels, dict)
+    assert isinstance(reference.noisy_channels["bad_by_manual"], list)
     assert isinstance(reference.noisy_channels_original, dict)
     assert isinstance(reference.bad_before_interpolation, list)
     assert isinstance(reference.reference_signal, np.ndarray)