[ENH] Add a bad channel detection method using LOF algorithm (#66)

Co-authored-by: Nicolas Barascud <10333715+nbara@users.noreply.github.com>

Author: vpKumaravel

README.md

@@ -3,7 +3,7 @@
 [![codecov](https://codecov.io/gh/nbara/python-meegkit/branch/master/graph/badge.svg)](https://codecov.io/gh/nbara/python-meegkit)
 [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/nbara/python-meegkit/master)
 [![DOI](https://zenodo.org/badge/117451752.svg)](https://zenodo.org/badge/latestdoi/117451752)
-[![twitter](https://img.shields.io/twitter/follow/lebababa?label=Twitter&style=flat&logo=Twitter)](https://twitter.com/intent/follow?screen_name=lebababa)
+[![twitter](https://img.shields.io/twitter/follow/lebababa?style=flat&logo=Twitter)](https://twitter.com/intent/follow?screen_name=lebababa)
 
 # MEEGkit
 
@@ -54,7 +54,7 @@ Other available options are `[docs]` (which installs dependencies required to bu
 
 ## References
 
-### 1. CCA, STAR, SNS, DSS, ZapLine, and robust detrending
+### 1. CCA, STAR, SNS, DSS, ZapLine, and Robust Detrending
 
 This is mostly a translation of Matlab code from the [NoiseTools toolbox](http://audition.ens.fr/adc/NoiseTools/) by Alain de Cheveigné. It builds on an initial python implementation by [Pedro Alcocer](https://github.com/pealco).
 
@@ -83,10 +83,9 @@ If you use this code, you should cite the relevant methods from the original art
  Journal of Neuroscience Methods, 168(1), 195–202. https://doi.org/10.1016/j.jneumeth.2007.09.012
 [10] de Cheveigné, A., & Simon, J. Z. (2007). Denoising based on time-shift PCA.
  Journal of Neuroscience Methods, 165(2), 297–305. https://doi.org/10.1016/j.jneumeth.2007.06.003
-
 ```
 
-### 2. Artifact subspace reconstruction (ASR)
+### 2. Artifact Subspace Reconstruction (ASR)
 
 The base code is inspired from the original [EEGLAB inplementation](https://github.com/sccn/clean_rawdata) [1], while the riemannian variant [2] was adapted from the [rASR toolbox](https://github.com/s4rify/rASRMatlab) by Sarah Blum.
 
@@ -101,7 +100,7 @@ If you use this code, you should cite the relevant methods from the original art
  13, 141.
 ```
 
-### 3. Rhythmic entrainment source separation (RESS)
+### 3. Rhythmic Entrainment Source Separation (RESS)
 
 The code is based on [Matlab code from Mike X. Cohen](https://mikexcohen.com/data/) [1]
 
@@ -130,3 +129,15 @@ If you use this, you should cite the following articles:
  "High-speed spelling with a noninvasive brain-computer interface",
  Proc. Int. Natl. Acad. Sci. U. S. A, 112(44): E6058-6067, 2015.
 ```
+
+### 5. Local Outlier Factor (LOF)
+
+If you use this, you should cite the following article:
+
+```sql
+[1] Breunig M, Kriegel HP, Ng RT, Sander J. 2000. LOF: identifying density-based 
+ local outliers. SIGMOD Rec. 29, 2, 93-104. https://doi.org/10.1145/335191.335388
+[2] Kumaravel VP, Buiatti M, Parise E, Farella E. 2022. Adaptable and Robust 
+ EEG Bad Channel Detection Using Local Outlier Factor (LOF). Sensors (Basel). 
+ 2022 Sep 27;22(19):7314. https://doi.org/10.3390/s22197314.
+```

doc/index.rst

@@ -37,6 +37,7 @@ Here is a list of the methods and techniques available in ``meegkit``:
  ~meegkit.cca
  ~meegkit.dss
  ~meegkit.detrend
+ ~meegkit.lof
  ~meegkit.ress
  ~meegkit.sns
  ~meegkit.star

doc/modules/meegkit.lof.rst

@@ -0,0 +1,23 @@
+meegkit.lof
+===========
+
+.. automodule:: meegkit.lof
+
+ .. rubric:: Functions
+
+ .. autosummary::
+
+ LOF
+
+
+
+
+
+
+
+
+
+
+
+
+

meegkit/__init__.py

@@ -1,7 +1,7 @@
 """M/EEG denoising utilities in python."""
 __version__ = '0.1.3'
 
-from . import asr, cca, detrend, dss, sns, star, ress, trca, tspca, utils
+from . import asr, cca, detrend, dss, lof, sns, star, ress, trca, tspca, utils
 
-__all__ = ['asr', 'cca', 'detrend', 'dss', 'ress', 'sns', 'star', 'trca',
+__all__ = ['asr', 'cca', 'detrend', 'dss', 'lof', 'ress', 'sns', 'star', 'trca',
  'tspca', 'utils']

meegkit/lof.py

@@ -0,0 +1,85 @@
+"""Local Outlier Factor (LOF)."""
+# Authors: Velu Prabhakar Kumaravel <vkumaravel@fbk.eu>
+# License: BSD-3-Clause
+
+import logging
+from sklearn.neighbors import LocalOutlierFactor
+
+
+class LOF():
+ """Local Outlier Factor.
+
+ Local Outlier Factor (LOF) is an automatic, density-based outlier detection
+ algorithm based on [1]_ and [2]_.
+
+ Parameters
+ ----------
+ n_neighbours : int
+ Number of neighbours defining the local neighbourhood.
+ metric: str in {'euclidean', 'nan_euclidean', 'cosine',
+ 'cityblock', 'manhattan'}
+ Metric to use for distance computation. Default is “euclidean”
+ threshold : float
+ Threshold to define outliers. Theoretical threshold ranges anywhere
+ between 1.0 and any integer. Default: 1.5
+
+ Notes
+ -----
+ It is recommended to perform a CV (e.g., 10-fold) on training set to
+ calibrate this parameter for the given M/EEG dataset.
+
+ See [2]_ for details.
+
+ References
+ ----------
+ .. [1] Breunig M, Kriegel HP, Ng RT, Sander J.
+ 2000. LOF: identifying density-based local outliers.
+ SIGMOD Rec. 29, 2, 93-104. https://doi.org/10.1145/335191.335388
+ .. [2] Kumaravel VP, Buiatti M, Parise E, Farella E.
+ 2022. Adaptable and Robust EEG Bad Channel Detection Using
+ Local Outlier Factor (LOF). Sensors (Basel). 2022 Sep 27;22(19):7314.
+ doi: 10.3390/s22197314. PMID: 36236413; PMCID: PMC9571252.
+
+ """
+
+ def __init__(self, n_neighbors=20, metric='euclidean',
+ threshold=1.5, **kwargs):
+
+ self.n_neighbors = n_neighbors
+ self.metric = metric
+ self.threshold = threshold
+
+ def predict(self, X):
+ """Detect bad channels using Local Outlier Factor algorithm.
+
+ Parameters
+ ----------
+ X : array, shape=(n_channels, n_samples)
+ The data X should have been high-pass filtered.
+
+ Returns
+ -------
+ bad_channel_indices : Detected bad channel indices.
+
+ """
+ if X.ndim == 3: # in case the input data is epoched
+ logging.warning('Expected input data with shape '
+ '(n_channels, n_samples)')
+ return []
+
+ if self.n_neighbors >= X.shape[0]:
+ logging.warning('Number of neighbours cannot be greater than the '
+ 'number of channels')
+ return []
+
+ if self.threshold < 1.0:
+ logging.warning('Invalid threshold. Try a positive integer >= 1.0')
+ return []
+
+ clf = LocalOutlierFactor(self.n_neighbors)
+ logging.debug('[LOF] Predicting bad channels')
+ clf.fit_predict(X)
+ lof_scores = clf.negative_outlier_factor_
+ bad_channel_indices = -lof_scores >= self.threshold
+
+ return bad_channel_indices

tests/data/lofdata.mat

@@ -9,8 +9,6 @@
 from meegkit.utils.matrix import sliding_window
 from scipy import signal
 
-np.random.seed(9)
-
 # Data files
 THIS_FOLDER = os.path.dirname(os.path.abspath(__file__))
 # file = os.path.join(THIS_FOLDER, 'data', 'eeg_raw.fif')
@@ -175,6 +173,7 @@ def test_asr_functions(show=False, method='riemann'):
 @pytest.mark.parametrize(argnames='reref', argvalues=(False, True))
 def test_asr_class(method, reref, show=False):
  """Test ASR class (simulate online use)."""
+ np.random.default_rng(9)
  raw = np.load(os.path.join(THIS_FOLDER, 'data', 'eeg_raw.npy'))
  sfreq = 250
  # Train on a clean portion of data

tests/test_asr.py

@@ -0,0 +1,37 @@
+"""LOF test."""
+import os
+
+import numpy as np
+import pytest
+import scipy.io as sio
+
+from meegkit.lof import LOF
+
+np.random.seed(9)
+
+# Data files
+THIS_FOLDER = os.path.dirname(os.path.abspath(__file__)) # data folder of MEEGKIT
+
+
+@pytest.mark.parametrize(argnames='n_neighbors', argvalues=(8, 20, 40, 2048))
+def test_lof(n_neighbors, show=False):
+ mat = sio.loadmat(os.path.join(THIS_FOLDER, 'data', 'lofdata.mat'))
+ X = mat['X']
+ lof = LOF(n_neighbors)
+ bad_channel_indices = lof.predict(X)
+ print(bad_channel_indices)
+
+@pytest.mark.parametrize(argnames='metric',
+ argvalues=('euclidean', 'nan_euclidean',
+ 'cosine', 'cityblock', 'manhattan'))
+def test_lof2(metric, show=False):
+ mat = sio.loadmat(os.path.join(THIS_FOLDER, 'data', 'lofdata.mat'))
+ X = mat['X']
+ lof = LOF(20, metric)
+ bad_channel_indices = lof.predict(X)
+ print(bad_channel_indices)
+
+if __name__ == "__main__":
+ pytest.main([__file__])
+ #test_lof(20, True)
+ #test_lof(metric='euclidean')

tests/test_lof.py

@@ -51,7 +51,8 @@ def test_tspca_sns_dss(): # TODO
  y_tspca_sns_dss = fold(
  np.dot(unfold(y_tspca_sns), todss), y_tspca_sns.shape[0])
 
- return y_tspca, y_tspca_sns, y_tspca_sns_dss
+ # TODO do something with it
+ assert y_tspca_sns_dss.shape == noisy_data.shape
 
 
 def test_tsr(show=True):
@@ -107,6 +108,7 @@ def test_tsr(show=True):
  plt.show()
 
 if __name__ == '__main__':
- # import pytest
- # pytest.main([__file__])
- test_tsr()
+ import pytest
+ pytest.main([__file__])
+ # test_tspca_sns_dss()
+ # test_tsr()