Repository: sappelhoff/pyprep (python)

Issue:
Bad channels by correlation differ between PyPREP and MATLAB PREP
In my recent comparison efforts, I noticed that the numbers of initial bad channels by correlation sometimes differed slightly between PyPREP and MATLAB PREP. Since this code doesn't involve any random number generation, it should always be consistent between runs and also between implementations.

Anyway, I did some digging with breakpoints and think I've figured it out. It has to do with [these lines here](https://github.com/sappelhoff/pyprep/blob/8476ea5c1f80c60e564ce8e4c2fca34464515b40/pyprep/find_noisy_channels.py#L320-L332):

```python
        for k in range(0, w_correlation):
            eeg_portion = np.transpose(np.squeeze(EEG_new_win[:, :, k]))
            data_portion = np.transpose(np.squeeze(data_win[:, :, k]))
            window_correlation = np.corrcoef(np.transpose(eeg_portion))
            abs_corr = np.abs(
                np.subtract(window_correlation, np.diag(np.diag(window_correlation)))
            )
            channel_correlation[k, :] = np.quantile(abs_corr, 0.98, axis=0)
            noiselevels[k, :] = np.divide(
                robust.mad(np.subtract(data_portion, eeg_portion), c=1),
                robust.mad(eeg_portion, c=1),
            )
            channel_deviations[k, :] = 0.7413 * iqr(data_portion, axis=0)
```
Using the same input data, we see that `eeg_portion`, `data_portion`, `window_correlation`, `abs_corr`, and `noiselevels` are all identical between PyPREP and MatPREP, but `channel_correlation` and `channel_deviations` both differ.

The reason for this is that apparently [MATLAB's quantile functions work slightly differently than Numpy's and SciPy's](https://stackoverflow.com/a/58424860): instead of having quantiles be restricted to the range [0, 1], they get restricted to the range [0.5 / n, (n - 0.5) / n], with n being the number of samples passed to the quantile function. This is apparently equivalent to treating the provided values as a sample of an unknown population rather than a full population themselves. After hacking together some new `quantile` and `iqr` functions that account for this, the `channel_correlation` and `channel_deviations` data both match the MATLAB values!

Anyway, I'll put together a PR to address this, but I figure there should be a discussion of which quantile logic actually makes the most sense for this use-case. If Numpy's default logic is actually more reasonable/justifiable here, we can make the new quantile functions optional with a `match_matlab` flag  or parameter of some kind indicating strict adherence to the original PREP logic. If MATLAB's logic makes more sense, I'm assuming we should just fully replace `np.quantile` and `iqr` with their modified equivalents?

Pull request:
Add MATLAB-style implementations of the quantile and IQR functions
# PR Description

Closes #56. As mentioned in that issue, MATLAB's quantile and IQR functions calculate quantiles slightly differently than Numpy and SciPy, leading to higher sensitivity to differences in channel correlation than in the original MATLAB PREP. This PR adds new `mat_quantile` and `mat_iqr` functions that match MATLAB's quantile logic and make PyPREP's internal correlation matrix for `NoisyChannels.find_bad_by_correlation` match the values in MATLAB PREP's equivalent matrix exactly.

### TODO

- I wasn't sure whether we'd want these functions to be part of the public API or private (i.e. `_mat_irq` and `_mat_quantile`) so I held off on writing docstrings.
- I also wasn't sure what we wanted to do here for unit testing. For the sake of this PR, would individual tests for these two functions be fine, or should we add full-pipeline testing of these changes as well? Eventually we'll want to modify the CI tests to do proper comparisons between MATLAB and PyPREP output values for each stage of NoisyChannels, but that would require figuring out some good test data we could use for both pipelines that wouldn't demand too much in terms of RAM or CPU, so I think that's a bigger problem to tackle in a separate PR.
- I need to update the changelog (which sort of depends on whether we want the functions to be public or private).

# Merge Checklist

<!--To merge your PR we need to first take the following points into account.-->
<!--Please just leave this checklist untouched-->

- [x] the PR has been reviewed and all comments are resolved
- [x] all [CI][what-is-ci] checks pass
- [x] (if applicable): the PR description includes the phrase `closes #<issue-number>` to [automatically close an issue][auto-close-documentation]
- [x] (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 8943538..cec22ab 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -42,6 +42,7 @@ Bug
 ~~~
 - Fixed RANSAC to give consistent results with a fixed seed across different chunk sizes, by `Austin Hurst`_ and `Yorguin Mantilla`_ (:gh:`43`)
 - Fixed "bad channel by flat" threshold in :meth:`NoisyChannels.find_bad_by_nan_flat` to be consistent with MATLAB PREP, by `Austin Hurst`_ (:gh:`60`)
+- Changed "bad channel by deviation" and "bad channel by correlation" detection code in :class:`NoisyChannels` to compute IQR and quantiles in the same manner as MATLAB, thus producing identical results to MATLAB PREP, by `Austin Hurst`_ (:gh:`57`)
 
 API
 ~~~
diff --git a/pyprep/find_noisy_channels.py b/pyprep/find_noisy_channels.py
index e0cd4ef..408d264 100644
--- a/pyprep/find_noisy_channels.py
+++ b/pyprep/find_noisy_channels.py
@@ -3,12 +3,11 @@ import mne
 import numpy as np
 from mne.utils import check_random_state
 from scipy import signal
-from scipy.stats import iqr
 from statsmodels import robust
 
 from pyprep.ransac import find_bad_by_ransac
 from pyprep.removeTrend import removeTrend
-from pyprep.utils import filter_design
+from pyprep.utils import filter_design, _mat_quantile, _mat_iqr
 
 
 class NoisyChannels:
@@ -195,8 +194,8 @@ class NoisyChannels:
         deviation_channel_mask = [False] * (self.new_dimensions[0])
         channel_deviation = np.zeros(self.new_dimensions[0])
         for i in range(0, self.new_dimensions[0]):
-            channel_deviation[i] = 0.7413 * iqr(self.EEGData[i, :])
-        channel_deviationSD = 0.7413 * iqr(channel_deviation)
+            channel_deviation[i] = 0.7413 * _mat_iqr(self.EEGData[i, :])
+        channel_deviationSD = 0.7413 * _mat_iqr(channel_deviation)
         channel_deviationMedian = np.nanmedian(channel_deviation)
         robust_channel_deviation = np.divide(
             np.subtract(channel_deviation, channel_deviationMedian), channel_deviationSD
@@ -325,12 +324,12 @@ class NoisyChannels:
             abs_corr = np.abs(
                 np.subtract(window_correlation, np.diag(np.diag(window_correlation)))
             )
-            channel_correlation[k, :] = np.quantile(abs_corr, 0.98, axis=0)
+            channel_correlation[k, :] = _mat_quantile(abs_corr, 0.98, axis=0)
             noiselevels[k, :] = np.divide(
                 robust.mad(np.subtract(data_portion, eeg_portion), c=1),
                 robust.mad(eeg_portion, c=1),
             )
-            channel_deviations[k, :] = 0.7413 * iqr(data_portion, axis=0)
+            channel_deviations[k, :] = 0.7413 * _mat_iqr(data_portion, axis=0)
         for i in range(0, w_correlation):
             for j in range(0, self.new_dimensions[0]):
                 drop[i, j] = np.int(
diff --git a/pyprep/utils.py b/pyprep/utils.py
index fd2a689..ae05ecb 100644
--- a/pyprep/utils.py
+++ b/pyprep/utils.py
@@ -5,6 +5,7 @@ from cmath import sqrt
 import mne
 import numpy as np
 import scipy.interpolate
+from scipy.stats import iqr
 from psutil import virtual_memory
 
 
@@ -20,6 +21,71 @@ def _intersect(list1, list2):
     return list(set(list1).intersection(set(list2)))
 
 
+def _mat_quantile(arr, q, axis=None):
+    """Calculate the numeric value at quantile (`q`) for a given distribution.
+
+    Parameters
+    ----------
+    arr : np.ndarray
+        Input array containing samples from the distribution to summarize.
+    q : float
+        The quantile to calculate for the input data. Must be between 0 and 1,
+        inclusive.
+    axis : {int, tuple of int, None}, optional
+        Axis along which quantile values should be calculated. Defaults to
+        calculating the value at the given quantile for the entire array.
+
+    Returns
+    -------
+    quantile : scalar or np.ndarray
+        If no axis is specified, returns the value at quantile (q) for the full
+        input array as a single numeric value. Otherwise, returns an
+        ``np.ndarray`` containing the values at quantile (q) for each row along
+        the specified axis.
+
+    Notes
+    -----
+    MATLAB calculates quantiles using different logic than Numpy: Numpy treats
+    the provided values as a whole population, whereas MATLAB treats them as a
+    sample from a population of unknown size and adjusts quantiles accordingly.
+    This function mimics MATLAB's logic to produce identical results.
+
+    """
+    q = np.asarray(q, dtype=np.float64)
+    n = len(arr)
+    q_adj = ((q - 0.5) * n / (n - 1)) + 0.5
+    return np.quantile(arr, np.clip(q_adj, 0, 1), axis=axis)
+
+
+def _mat_iqr(arr, axis=None):
+    """Calculate the inter-quartile range (IQR) for a given distribution.
+
+    Parameters
+    ----------
+    arr : np.ndarray
+        Input array containing samples from the distribution to summarize.
+    axis : {int, tuple of int, None}, optional
+        Axis along which IQRs should be calculated. Defaults to calculating the
+        IQR for the entire array.
+
+    Returns
+    -------
+    iqr : scalar or np.ndarray
+        If no axis is specified, returns the IQR for the full input array as a
+        single numeric value. Otherwise, returns an ``np.ndarray`` containing
+        the IQRs for each row along the specified axis.
+
+    Notes
+    -----
+    See notes for :func:`utils._mat_quantile`.
+
+    """
+    iqr_q = np.asarray([25, 75], dtype=np.float64)
+    n = len(arr)
+    iqr_adj = ((iqr_q - 50) * n / (n - 1)) + 50
+    return iqr(arr, rng=np.clip(iqr_adj, 0, 100), axis=axis)
+
+
 def filter_design(N_order, amp, freq):
     """Create FIR low-pass filter for EEG data using frequency sampling method.
 


Test patch:
diff --git a/tests/test_utils.py b/tests/test_utils.py
new file mode 100644
index 0000000..1eb36b2
--- /dev/null
+++ b/tests/test_utils.py
@@ -0,0 +1,37 @@
+"""Test various helper functions."""
+import numpy as np
+
+from pyprep.utils import _mat_quantile, _mat_iqr
+
+
+def test_mat_quantile_iqr():
+    """Test MATLAB-compatible quantile and IQR functions.
+
+    MATLAB code used to generate the comparison results:
+
+    .. code-block:: matlab
+
+       % Generate test data
+       rng(435656);
+       tst = rand(100, 3);
+
+       % Calculate IQR and 0.98 quantile for test data
+       quantile(tst, 0.98);
+       iqr(tst);
+
+    """
+    # Generate test data
+    np.random.seed(435656)
+    tst = np.transpose(np.random.rand(3, 100))
+
+    # Create arrays containing MATLAB results
+    quantile_expected = np.asarray([0.9710, 0.9876, 0.9802])
+    iqr_expected = np.asarray([0.4776, 0.5144, 0.4851])
+
+    # Test quantile equivalence with MATLAB
+    quantile_actual = _mat_quantile(tst, 0.98, axis=0)
+    assert all(np.isclose(quantile_expected, quantile_actual, atol=0.001))
+
+    # Test IQR equivalence with MATLAB
+    iqr_actual = _mat_iqr(tst, axis=0)
+    assert all(np.isclose(iqr_expected, iqr_actual, atol=0.001))