Hi! We'd like to implement the entire PREP pipeline, which will include:

1. Bug fixing of Noisydata class
2. Feature enhancement by adding filtering, line noise removal and robust referencing.
3. Integrating all the steps into one function

MemoryError: Unable to allocate 9.75 GiB for an array with shape (6545, 100001) and data type complex128

It just cannot work on a big dataset, which is just 500mb.

both the pyprep and the autoreject package have an implementation of RANSAC.

We should make a comparison to:

1. see whether they provide the same results
2. potentially find bugs in one of the implementation
3. check which implementation is more computationally efficient

After the comparison we could decide to drop one implementation and only use the "better" one.

Or we keep using both implementations, but make an example that explicitly compares them, pointing out pros and cons in terms of features.

In any case, this would be an important enterprise for point 1 and 2 in my list alone.

Hi @sappelhoff, @Nick3151 and @AamnaLawrence here are working on migrating the full functionality of PREP into this forked version of pyprep. I'm overseeing the project and will help wherever they are lacking some skillsets.

This will involve a variety of factors:

1. Scikit-learn "fit" style pipeling
2. Turning things into classes where it is relevant.
3. Adding additional unit/integration tests for the above
4. Adding python3 into travis.ci

I'm just making an issue here, so you can comment on any of these above things and guide dev. For PRs on a forked version, you can check out dev at anytime at: https://github.com/NeuroDataDesign/pyprep/pulls.

Lmk your thoughts!

Currently the prepPipeline object has access to the bad channels before and after interpolation but it doesn't identify which criteria marked which channel. In the original PREP this info is in the EEG.etc.referenceOut.noisyStatisticsBeforeInterpolation and afterInterpolation.

I understand that maybe for most users this info won't be needed but for some power users and contributors may find this useful for catching bugs, inconsistencies and in general test the criterias not isolated from each other and in the context of the whole pipeline.

I approached this by making the noisy detectors attributes of the reference class as follows:

        # Phase 2: Find the bad channels and interpolate
        self.raw._data = self.EEG * 1e-6
        self.noisy_detector_before_interpolation = NoisyChannels(self.raw)
        self.noisy_detector_before_interpolation.find_all_bads(ransac=self.ransac)

        # Record Noisy channels and EEG before interpolation
        self.bad_before_interpolation = self.noisy_detector_before_interpolation.get_bads(verbose=True)
        self.EEG_before_interpolation = self.EEG.copy()

.................

        # Still noisy channels after interpolation
        self.interpolated_channels = bad_channels
        self.noisy_detector_after_interpolation = NoisyChannels(self.raw)
        self.noisy_detector_after_interpolation.find_all_bads(ransac=self.ransac)
        self.still_noisy_channels = self.noisy_detector_after_interpolation.get_bads()

And then in prepPipeline object do:

        self.noisy_detector_before_interpolation = \
        { \
        'bad_by_nan' : reference.noisy_detector_before_interpolation.bad_by_nan,
        'bad_by_flat' : reference.noisy_detector_before_interpolation.bad_by_flat,
        'bad_by_deviation' : reference.noisy_detector_before_interpolation.bad_by_deviation,
        'bad_by_hf_noise' : reference.noisy_detector_before_interpolation.bad_by_hf_noise,
        'bad_by_correlation' : reference.noisy_detector_before_interpolation.bad_by_correlation,
        'bad_by_SNR' : reference.noisy_detector_before_interpolation.bad_by_SNR, 
        'bad_by_dropout' : reference.noisy_detector_before_interpolation.bad_by_dropout,
        'bad_by_ransac' : reference.noisy_detector_before_interpolation.bad_by_ransac
        }


        self.noisy_detector_after_interpolation = \
        { \
        'bad_by_nan' : reference.noisy_detector_after_interpolation.bad_by_nan,
        'bad_by_flat' : reference.noisy_detector_after_interpolation.bad_by_flat,
        'bad_by_deviation' : reference.noisy_detector_after_interpolation.bad_by_deviation,
        'bad_by_hf_noise' : reference.noisy_detector_after_interpolation.bad_by_hf_noise,
        'bad_by_correlation' : reference.noisy_detector_after_interpolation.bad_by_correlation,
        'bad_by_SNR' : reference.noisy_detector_after_interpolation.bad_by_SNR, 
        'bad_by_dropout' : reference.noisy_detector_after_interpolation.bad_by_dropout,
        'bad_by_ransac' : reference.noisy_detector_after_interpolation.bad_by_ransac
        }

I do know that this is somewhat cumbersome to the eyes but I believe this helps greatly to examine the pyprep with more detail without having to do debug to get inside the functions.

Looking forward to your feedback.

Dear pyprep-team,

I've used pyprep for preprocessing my EEG data and now want to do additional preprocessing steps including ICA. For ICA I want to use the EOG channel as a reference channel to select ICA components. Unfortunately, during preprocessing the eog channels were dropped from the info structure and it does not seem to be possible to add the back. Is this an issue with the pyprep pipeline or is there a simple workaround?

Before pyprep:

`<Info | 8 non-empty values
bads: []
ch_names: Fp1, AF7, AF3, F1, F3, F5, F7, FT7, FC5, FC3, FC1, C1, C3, C5, ...
chs: 64 EEG, 4 EOG
custom_ref_applied: False
dig: 67 items (3 Cardinal, 64 EEG)
highpass: 0.0 Hz
lowpass: 104.0 Hz
meas_date: 2018-04-06 09:34:26 UTC
nchan: 68
projs: []
sfreq: 512.0 Hz

`

After pyprep:

`<Info | 8 non-empty values
bads: []
ch_names: Fp1, AF7, AF3, F1, F3, F5, F7, FT7, FC5, FC3, FC1, C1, C3, C5, ...
chs: 64 EEG
custom_ref_applied: False
dig: 67 items (3 Cardinal, 64 EEG)
highpass: 0.0 Hz
lowpass: 50.0 Hz
meas_date: 2018-04-06 09:34:26 UTC
nchan: 64
projs: []
sfreq: 100.0 Hz

`

Hi, I'm Yorguin I've been studying the prep code for some time now. I'm testing differences between matlab's results and the current pyprep implementation. I have noticed that when I run pyprep I get an unusually high number of bad channels and in particular for the ransac criteria. For example:

If we take the following raw file: https://dataverse.tdl.org/file.xhtml?fileId=37645&version=1.0 which correspond to a public dataset we get:

Interpolated Channels Matlab : [10, 41, 44, 50, 65] (0 index corrected)

Interpolated Channels Python: [11, 17, 46, 36, 65, 10, 52, 9, 6, 4, 57, 8, 29, 59, 12, 26, 18, 13, 58, 42, 35, 56, 61, 15, 43, 7, 51, 44, 62, 32, 41, 14, 20, 47]

I'm currently working on it and I think I have found the mistake (it is actually outside the ransac) so I made a fork where I will put the changes that I made. With the changes I implemented we get: [10,41,44,49, 50, 52 ,54, 65]

** This was done skipping the line noise multitaper notch (for both matlab and pyprep) because that took too long to run the tests...

Apart from this I was able to find out a way to make the ransac less memory hungry using either one of two methods:

• Automatically calculate the resample rate needed for the ransac to fit into memory.
• Calculate the ransac in a more channels wise manner.

There are also some little discrepancies in the logic of the Matlab PREP and in pyprep here and there. Particularly in performReference method.

I will be uploading those changes to the fork in a separated manner so it is easier to keep track of the changes.

Let me know what you think about all of this.

I've been picking away at refactoring the NoisyChannels unit tests so that each 'bad-by' type is covered by a separate test. In the process, though, I think I've discovered a weak spot in the original PREP logic: if I use a different test file, the 'bad-by-low-deviation test' below stops working:

This holds true even if I replace the "divide by 10" with a "divide by 100,000": the "robust channel deviation" z-score seems to plateau at around -4, never reaching the +/- 5 threshold for being bad by deviation.

Looking at the actual math, this makes sense:

1. PREP calculates the variance within each channel using 0.7413 * the interquartile range for the signal.
2. PREP calculates the median and variance (again, 0.7413 * the IQR) of those values and uses them to do a non-traditional Z-score of the channel variances (i.e., (variance - median_variance) / variance_of_the_variance).
3. PREP compares the absolute values of those Z-scores to a threshold (default = 5) and flags any channels that exceed it as "bad-by-deviation"

The problem here is that in step 2, the variances calculated in step 1 have a minimum value of zero. As such, even a channel with an amplitude of zero won't be detected as bad-by-deviation if the median isn't at least 5x the variance of the variances.

A quick-and-dirty fix is to z-score the log of the channel variances, which makes the variances more normally-distributed and thus makes the detection of low-amplitude channels much easier. However, this also increases the threshold for detecting high-amplitude bad-by-deviation channels, with the required multiplication factor going from 2.4 to 3.9 for my test channel.

Any ideas on how to better handle this?

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:

gets reused here:

whereas in MATLAB PREP the former is a local variable and doesn't affect the latter.

I just today realized that to make intersphinx and numpydoc work together, on needs to set numpydoc_xref_param_type = True in conf.py.

I did that, and it works in many other projects ... but somehow not in pyprep. when building the docs, I get these warnings (see below)

@larsoner, given your expertise with both numpydoc and sphinx - do you know what's going on? why is intersphinx not detecting the "bool" type?

link to docs where links are broken: https://pyprep.readthedocs.io/en/latest/generated/find_noisy_channels.NoisyChannels.html#find_noisy_channels.NoisyChannels

/find_noisy_channels.py:docstring of find_noisy_channels.NoisyChannels.__init__:12: WARNING: undefined label: python:bltin-boolean-values
/find_noisy_channels.py:docstring of find_noisy_channels.NoisyChannels.find_all_bads:22: WARNING: undefined label: python:bltin-boolean-values
/find_noisy_channels.py:docstring of find_noisy_channels.NoisyChannels.get_bads:23: WARNING: undefined label: python:bltin-boolean-values

/find_noisy_channels.py:docstring of find_noisy_channels.NoisyChannels.get_ransac_pred:14: WARNING: term not in glossary: numpy:array_like
/find_noisy_channels.py:docstring of find_noisy_channels.NoisyChannels.ransac_correlations:17: WARNING: term not in glossary: numpy:array_like
/find_noisy_channels.py:docstring of find_noisy_channels.NoisyChannels.run_ransac:16: WARNING: term not in glossary: numpy:array_like
/prep_pipeline.py:docstring of prep_pipeline.PrepPipeline:32: WARNING: undefined label: python:bltin-boolean-values

numpydoc config:

interspinx config:

example of what's apparently broken:

see the corresponding error:

/find_noisy_channels.py:docstring of find_noisy_channels.NoisyChannels.init:12: WARNING: undefined label: python:bltin-boolean-values

When RANSAC does its random channel selection, the random channel choices diverge between MATLAB PREP and PyPREP, even when using the same seed. This isn't super-surprising since they're different languages with different randomization functions, but it would be handy if we could PyPREP to make the same picks as MatPREP given the same seed.

In my preliminary efforts looking into this, I've learned a few things:

• Most importantly, Numpy offers an identical implementation of the "twister" randomization algorithm to MATLAB such that they should be able to produce identical random results given the same seed
• MATLAB PREP implements its own custom randsample function for selecting subsets of channels without replacement. If I re-implement the same logic in Python I identical results.
• MATLAB PREP also has a getRandomSubsets function that calls the above randsample function and produces different numbers yet again. Haven't tried re-implementing it in Python yet, but hopefully it's trivial!

With some luck I'll have a PR for this shortly. It's going to break seed compatibility with earlier PyPREP versions, but I think we've already done that since last release in #43.

Whenever I try to run PyPREP I get the following output which seems promising, however it stalls here and continues to act as if it is running and no other output comes out. Any assistance would be greatly appreciated. (Also I tried running this on both the present master version v20, and also the development version and the same results occur in each one.)

EDIT by @sappelhoff : 2020-06-06 --> I added backticks ````` to format the following as a codeblock. Please see the guidelines here @barlehmann

Extracting EDF parameters from /content/drive/My Drive/Colab Notebooks/Jacob5MeO2min.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 200499  =      0.000 ...   801.996 secs...
<Info | 7 non-empty values
 bads: []
 ch_names: fp1, fp2, f3, f4, c3, c4, p3, p4, o1, o2, f7, f8, t3, t4, t5, ...
 chs: 19 EEG
 custom_ref_applied: False
 highpass: 0.0 Hz
 lowpass: 128.0 Hz
 meas_date: unspecified
 nchan: 19
 projs: []
 sfreq: 256.0 Hz
>
[<DigPoint |        LPA : (-82.5, -0.0, 0.0) mm     : head frame>, <DigPoint |     Nasion : (0.0, 102.7, 0.0) mm      : head frame>, <DigPoint |        RPA : (82.2, 0.0, 0.0) mm       : head frame>, <DigPoint |     EEG #1 : (-28.2, 102.3, 31.7) mm   : head frame>, <DigPoint |     EEG #3 : (28.6, 103.2, 31.6) mm    : head frame>, <DigPoint |    EEG #16 : (-67.4, 62.3, 30.5) mm    : head frame>, <DigPoint |    EEG #18 : (-48.2, 76.3, 80.9) mm    : head frame>, <DigPoint |    EEG #20 : (0.3, 83.3, 103.8) mm     : head frame>, <DigPoint |    EEG #22 : (49.7, 77.4, 79.6) mm     : head frame>, <DigPoint |    EEG #24 : (70.0, 64.2, 29.8) mm     : head frame>, <DigPoint |    EEG #40 : (-62.7, 16.1, 106.8) mm   : head frame>, <DigPoint |    EEG #42 : (0.4, 21.0, 140.9) mm     : head frame>, <DigPoint |    EEG #44 : (64.3, 16.7, 106.0) mm    : head frame>, <DigPoint |    EEG #62 : (-50.8, -48.7, 103.6) mm  : head frame>, <DigPoint |    EEG #64 : (0.3, -49.0, 129.2) mm    : head frame>, <DigPoint |    EEG #66 : (53.3, -48.5, 104.2) mm   : head frame>, <DigPoint |    EEG #81 : (-28.2, -84.3, 61.0) mm   : head frame>, <DigPoint |    EEG #83 : (28.6, -84.0, 60.9) mm    : head frame>, <DigPoint |    EEG #87 : (-80.7, 6.6, 36.6) mm     : head frame>, <DigPoint |    EEG #88 : (-69.4, -47.8, 47.3) mm   : head frame>, <DigPoint |    EEG #89 : (81.5, 7.5, 36.5) mm      : head frame>, <DigPoint |    EEG #90 : (70.0, -47.4, 47.3) mm    : head frame>]
[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18]
['fp1', 'fp2', 'f3', 'f4', 'c3', 'c4', 'p3', 'p4', 'o1', 'o2', 'f7', 'f8', 't3', 't4', 't5', 't6', 'fz', 'cz', 'pz']
Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 845 samples (3.301 sec)

<ipython-input-15-246072ead167>:34: DeprecationWarning: Using ``raise_if_subset`` to ``set_montage``  is deprecated and ``set_dig`` will be  removed in 0.21
  raw.set_montage(montage, match_case=False, raise_if_subset=False)
Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 845 samples (3.301 sec)

After fixing the correlation, deviation, and channel picks code to match MATLAB PREP, I finally started trying to dig into differences in RANSAC interpolation between implementations. In the process, however, I noticed that the channel picks still weren't lining up between PyPREP and MATLAB PREP despite my earlier efforts.

Some deeper digging revealed that this is because PyPREP and MATLAB PREP are starting with different pools of "good" channels: PyPREP with 58 and MATLAB PREP with 60! Apparently, MATLAB PREP doesn't exclude channels bad by HF noise prior to RANSAC: only channels bad by NaN, flat, deviation, correlation, or dropout (note: NaN and flat channels already removed earlier in code):

%% Bad so far by amplitude and correlation (take these out before doing ransac)
noisyChannels = union(noisyOut.noisyChannels.badChannelsFromDeviation, ...
    union(noisyOut.noisyChannels.badChannelsFromCorrelation, ...
          noisyOut.noisyChannels.badChannelsFromDropOuts));

Presumably this is because RANSAC applies a 50 Hz lowpass filter before doing correlations anyway, so HF noise wouldn't necessarily be a priori grounds for exclusion. Does this make sense, or does it look like a bug in the MatPREP code? Is this a behaviour we'd like to mimic in PyPREP?

This bit has been harder to compare than earlier parts, but I think I can say with a high degree of confidence now that RANSAC correlations differ between PyPREP and MATLAB PREP.

Before I get into the details, some other relevant stuff I've learned:

• MatPREP's internal spherical_interpolation function and MNE's internal _make_interpolation_matrix functions produce identical outputs, given the same inputs! This is a huge relief, and one less thing we have to worry about or figure out.
• Conversely, there's something weird about how MNE's mne.channels.read_custom_montage imports .elp files, like the one bundled with EEGLAB we've been using for the MatPREP data. Specifically, if I just load my pre-cleanlined and montage'd .set/.fdt files into MNE, the Cartesian X/Y/Z coordinates for the electrodes more-or-less line up with MATLAB (the X and Y axes are swapped) and the interpolation matrices match exactly. However, if I re-read the standard-10-5-cap385.elp montage into MNE and apply it to the file, the coordinates and interpolation matrices both differ quite a bit. Not sure what this is all about, but regardless I don't think it's a PyPREP bug.

Anyway, the gist of the difference (and why it's been so hard to dig into): MatPREP and PyPREP have two very different approaches of getting the same RANSAC correlation results. Here's a quick overview:

MATLAB PREP's RANSAC:

1. Choose random channels for each RANSAC sample.
2. Generate interpolation matrix for all channels for each RANSAC sample.
3. For each 5 second RANSAC window in the data, calculate the predicted signal for each channel across each of the 50 RANSAC samples, and then correlate the predicted values with the actual values to get a correlation for each channel in that window.
4. After calculating a matrix of correlations of size [number_of_ransac_windows, number_of_channels], apply thresholding and identify channels that fall below target.

PyPREP's RANSAC:

1. Choose random channels for each RANSAC sample.
2. Determine the amount of available RAM and split the channels-to-predict into chunks accordingly.
3. For each chunk of channels to predict:
   • For each of the 50 RANSAC samples, calculate the interpolation matrix and predict the full signal for those channels, saving the output into one large eeg_predictions matrix.
   • After interpolating the signal for all channels in the chunk for all 50 samples, calculate the median EEG signal for each channel from those 50 individual predictions.
   • For each 5 second RANSAC window in the data, calculate the correlation between the median predicted signal and the actual signal for the channels in that chunk.
4. After calculating a matrix of correlations of size [number_of_ransac_windows, number_of_channels], apply thresholding and identify channels that fall below target.

Anyway, the output correlations of each implementation differ, even with the same input. I know there's a deliberate difference somewhere regarding a median calculation, but I couldn't remember what it was and given how different the code layout is it hasn't popped out to me. @sappelhoff @yjmantilla @jasmainak do any of you remember what this difference was, so I can try reversing it and see if it makes things match up?

Also regardless of other factors, the MatPREP approach seems to have much lower RAM demands since at no point are they estimating signals across the whole duration of the file. Should we consider changing PyPREP to work in a similar manner?

While chugging along on the PyPREP/MatPREP comparison unit tests, I had a thought about the PyPREP API: since there are more-or-less three distinct stages to the PREP pipeline (adaptive line noise removal, robust re-referencing, interpolation of remaining bad channels), what if the PrepPipeline object had separate methods for each of these? The fit() method would still exist, but just run all of these in sequence.

Separate methods would be handy for situations like my own, where I basically copy/pasted the guts of fit() into my own script in order to do progress messages and plotting between CleanLine and robust referencing. Also for my script the project lead wanted the option to ignore bad channels instead of interpolating them, so I had to write extra code to revert interpolation and the names of bad channels instead of having the option to not interpolate in the first place.

Here's my basic proposal for PrepPipeline methods:

• fit: Performs the whole PREP pipeline, like ususal
• remove_line_noise: Performs adaptive line noise removal.
• robust_reference: Performs robust re-referencing without interpolation. Will default to printing a warning if remove_line_noise isn't done first.
• interpolate_bads: Interpolates any remaining bad channels and adjusts the average reference accordingly.

Let me know what you think!

EDIT: With the addition of #96, a PR addressing this issue should also add MATLAB comparison tests for final bad channel interpolation.

flake8-docstrings is not catching issues like

def tmp(foo):
    """Calculate an example.

    Parameters
    ----------
    foo: int
        This is an example.
    """
    return foo + 1

where a space is missing between foo and :

clicking on classes, methods, or functions documented in the changelog, should directly lead users to the API docs.

But it's not working ... the properly documented classes, methods, etc. do show up in bold, but their underlying html link is missing.

https://pyprep.readthedocs.io/en/latest/whats_new.html

Similar to #17 the MNE notch filter consumes quite a lot of RAM.

For example in my dataset I have 58 channels @ 1000Hz and 5 minutes. This will need somewhere around 9gb of free ram. For personal computers it is quite a lot. Now this would require either looking a way to optimize the notch of MNE and fixing it overthere or make a custom one for pyprep.

Nevertheless in my case my data has already a notch filter for the line noise and maybe some user will already have their data which such filters so I propose for that case to add the feature of optional line noise filtering. This is easily done in prep_pipeline.py by assessing whether line_freqs is empty and in that case to skip directly to the referencing process.

This would look like the following:

    if len(self.prep_params["line_freqs"]) != 0:

            # Step 1: 1Hz high pass filtering
            self.EEG_new = removeTrend(self.EEG_raw, sample_rate=self.sfreq)

            # Step 2: Removing line noise
            linenoise = self.prep_params["line_freqs"]

            self.EEG_clean = mne.filter.notch_filter(
                self.EEG_new,
                Fs=self.sfreq,
                freqs=linenoise,
                method="spectrum_fit",
                mt_bandwidth=2,
                p_value=0.01,
            )

            # Add Trend back
            self.EEG = self.EEG_raw - self.EEG_new + self.EEG_clean
            self.raw._data = self.EEG * 1e-6

And then the referencing code.

It is important to note that this can be done since the first detrend is temporal and only affects the notch filter since the trend is added back at the end. If we don't do the notch then the first detrend doesn't do anything so thats why it is inside the if.

In MATLAB PREP, robust re-referencing starts by detecting noisy channels, calculating the median signal excluding all channels bad by NaN, bad by flat, or bad by low SNR, and then re-referencing the data using that signal. PyPREP does almost the same thing, except it leaves in channels bad by low SNR which results in numeric differences when those sorts of channels are present:

Based on the above code & comment, it looks like this difference is deliberate? Is there an intentional reason to include these channels when calculating a median reference? If so, I'll make it a matlab_strict thing and document it in the "deliberate differences" page. If not, I'll just make a PR to add those in.

Looking at the MatPREP Git history, it looks like the SNR channels were a deliberate addition to the list of exclusions from the initial reference, added when the default method of calculating the initial reference signal was changed from "huber mean" (?) to median. Maybe it's worth asking Dr. Robbins about?

Anyway, if I add "bad-by-SNR" channels to the list of initial unusables in PyPREP, the initial average-referenced signal is identical to the one in MATLAB PREP.

So I guess we could test that by replacing the current interpolation with the interpolation method of raw.

Basically ,

• use the self.raw_mne copy we have of the object, copy it
• drop bads channels indicated by prior methods
• replace the data with the one that the ransac uses
• Mark as bad the current channels that are being reconstructed
• raw.interpolate_bads()
• Get the data out and adequate it for ransac to compute the correlations

We could also imitate the procedure to get the origin parameter instead of implementing raw.interpolate_bads() directly but the advantage of the former is that it would in essence be exactly what we do outside in the robust_reference algorithm.

I may give it a try when I get out of an exam I have on Thursday, or at worst after the 18 of Dec when Im finally on vacation.

Originally posted by @yjmantilla in #41 (comment)

Before the contributions of @AamnaLawrence @adam2392 and @Nick3151 in #6 there was only the single module noisy.py implementing the "find_noisy_channels" functionality.

Now the full PREP is implemented, including a new "find_noisy_channels" functionality, see: https://github.com/sappelhoff/pyprep/blob/master/pyprep/find_noisy_channels.py

noisy.py should be removed and the find_noisy_channels.py module should be overhauled so that it can be used standalone, without having to run the full PREP.

Haven't really dug into it much yet, but I've verified that adaptive line noise removal (a.k.a. CleanLine) differs somehow between PyPREP and Matlab PREP.

This is sort of a weird one, where the MNE notch filter that PyPREP uses seems to be very aggressive and distorts the waveband, whereas the MatPREP version doesn't seem to do much of anything (there's a very subtle difference around the 50Hz and 100Hz bands if you a/b the images). I should probably try different test files though, since average-referencing seems to handle 99% of the line noise distortions in the PSD for this data (is this a feature of CMS/DRL I didn't know about?)

Pre-cleanline data:

Post-CleanLine data (MatPREP):

Post-CleanLine data (PyPREP):

Anyway, I think this the final component I hadn't formally tested yet for #58!

Dear pyprep group
I try to compare the design difference between prep pipeline between pyprep and original on in eeglab. I am a little confused about whether pyprep work like original version, which contain a postprocess.m file including postprocess(interpolate the bad channels if you want to and high-pass filter if you want to).

This is not something to do right now, but as far as I heard it's now possible to test actual MATLAB in some CI services, see for reference: bids-standard/bids-matlab#68 (comment)

That could mean we can in our CI

1. download some lightweight data
2. run PREP on it (matlab), save outputs as build artifacts
3. run pyprep on it (python), save outputs as build artifacts
4. run another python script that compares the two build artifacts

In reference.py we have the method perform_reference() which states the following:

        # Phase 1: Estimate the true signal mean with robust referencing
        self.robust_reference()
        if self.noisy_channels["bad_all"]:
            self.raw.info["bads"] = self.noisy_channels["bad_all"]
            self.raw.interpolate_bads() <----------------------------------NOTICE THIS
        self.reference_signal = (
            np.nanmean(self.raw.get_data(picks=self.reference_channels), axis=0) * 1e6
        )
        rereferenced_index = [
            self.ch_names_eeg.index(ch) for ch in self.rereferenced_channels
        ]
        self.EEG = self.remove_reference(
            self.EEG, self.reference_signal, rereferenced_index
        )

Now in the original Matlab code this is equivalent to doRobustPost() inside performReference.m:

    function [] = doRobustPost()
        % Robust reference with interpolation afterwards
        referenceOut = robustReference(signal, referenceOut);
        noisy = referenceOut.badChannels.all;
        if isempty(noisy)   %No noisy channels -- ordinary ref
            referenceOut.referenceSignal = ...
                nanmean(signal.data(referenceOut.referenceChannels, :), 1);
        else
            sourceChannels = setdiff(referenceOut.evaluationChannels, noisy);
            signalNew = interpolateChannels(signal, noisy, sourceChannels); <----- NOTICE THIS
            referenceOut.referenceSignal = ...
                nanmean(signalNew.data(referenceOut.referenceChannels, :), 1);
            clear signalNew;
        end
        signal = removeReference(signal, referenceOut.referenceSignal, ...
            referenceOut.rereferencedChannels);
        referenceOut.noisyStatistics  = ...
            findNoisyChannels(removeTrend(signal, referenceOut), referenceOut);
        

So the basic difference is that robustReference produce a list of bad channels and this are interpolated in a dummy variable then averaged to produce a robust reference. In pyprep the interpolation is done on the raw object directly, not on a dummy copy of it.

I'm not so sure if this affects greatly or not but I would like your opinion on this.

If we wanted to make the code more alike the original one we would:

        # Phase 1: Estimate the true signal mean with robust referencing
        self.robust_reference()
        if self.noisy_channels["bad_all"]:
            dummy = self.raw.copy()
            dummy.info["bads"] = self.noisy_channels["bad_all"]
            dummy.interpolate_bads()
            self.reference_signal = (
                np.nanmean(dummy.get_data(picks=self.reference_channels), axis=0) * 1e6
            )
            del dummy
        else:
            self.reference_signal = (
                np.nanmean(self.raw.get_data(picks=self.reference_channels), axis=0)
                * 1e6
            )
        rereferenced_index = [
            self.ch_names_eeg.index(ch) for ch in self.rereferenced_channels
        ]
        self.EEG = self.remove_reference(
            self.EEG, self.reference_signal, rereferenced_index
        )

I have been testing pyprep against our datasets (GRUNECO, University of Antioquia, Colombia). For example I have 58 channels @ 1000Hz for about 5 minutes. This will need like 10 gb of RAM awarded from the OS to python. This is quite a lot for personal computers.

I have two possible proposals to solve this:

• Run RANSAC in a channel wise manner.
• Automatically identify the resample needed for the data to fit memory

I like more the first proposal and @sappelhoff agrees with this approach. This was discussed in #9 .

Once #9 detrend_correction is merged into master I will upload this code so it may be reviewed and tested. It considerably changes find_bad_by_ransac() and a tiny change in run_ransac()

Working on #79, I modified matlab_artifacts to simulate a "dropout" channel (I set the second 25% and last 25% of the channel to 0). This channel is successfully flagged as a bad-by-dropout by MatPREP, but in PyPREP the addition of this channel breaks NoisyChannels completely, seriously messing up the correlation windows and leaving too few usable channels to use for RANSAC.

I've been looking into it for a bit and it seems to be due to differences in how NaNs are generated during correlation between Numpy and MATLAB, but haven't figured out the exact source yet.

Splitting this off from #49 to keep things organized. Basically, the removeTrend function (which essentially just runs a highpass filter to remove frequencies below 1 Hz) produces different values in PyPREP's implementation (by default, a wrapper around mne.filter.filter_data) than in MATLAB PREP's implementation (by default, a wrapper around EEGLAB's pop_eegfiltnew).

For comparison I loaded the same .set file into both EEGLAB and MNE, applied removeTrend to both, then re-saved the .set from EEGLAB so I could load it into MNE and compare. The PSD plots for both look identical to me, but the numeric values themselves via raw._data look quite different:

signal.data(:, 1)

ans =

   1.0e+04 *

   -0.3335
    1.0037
   -0.3581
   -0.1795
    0.9913
   -0.7414
   -0.7198
   -0.5665
    0.5557
   -1.2237

signal2.data(:, 1)

ans =

   -3.7068
   -4.8012
   -0.8008
   -2.5809
   -1.5205
    0.1093
    1.9928
   -4.6906
  -10.2859
    1.6814

raw_copy._data[:10, 0]
array([-0.00333541,  0.01003745, -0.00358109, -0.00179516,  0.00991251,
       -0.00741443, -0.00719762, -0.00566537,  0.00555661, -0.01223686])

raw_filtered._data[:10, 0]
array([-4.11658012e-19, -1.79973325e-18,  6.98802181e-19,  2.89261751e-19,
        2.50128829e-18, -9.28348110e-19,  7.57247455e-19,  3.55753838e-19,
        3.82858892e-19, -7.01343280e-19])

Since the plots look the same the differences might not practically matter all that much, but if we want to be able to do formal (or informal) 1:1 comparisons of PyPREP results to MatPREP, having identical results from this function would be helpful, especially since this always happens before any of the NoisyChannels operations.

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:

        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: 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?

In MATLAB PREP, the random seed for RANSAC during findNoisyChannels is always set to the same value (435656) immediately before RANSAC is run. This means that, given the same set of 'good' channels, two consecutive runs of RANSAC during re-referencing will have identical channel picks.

In PyPREP, the random seed is set at the onset of PrepPipeline or Reference and gets made into a Numpy RandomState object. This means that the random state is going to be different for every RANSAC run during re-referencing, resulting in different re-referencing results than MATLAB PREP and potentially higher variability between iterations (resulting in more channels being flagged as bad). Actually, given how random subsets are calculated in PyPREP/MatPREP (generating random floats from 0 to 1, multiplying them by the number of channels to choose from, and then popping the channel with the nearest index repeatedly in a loop), using the same seed across runs would mean that similar channels would get chosen across RANSAC runs, assuming that adjacent channels in the data are often also going to be adjacent spatially.

This can be easily fixed with a PR, but before I do I just want to double-check whether this is a behaviour we want to match? I mean, MATLAB's approach would lead to less variability in NoisyChannels results during re-referencing, but it would also mean that RANSAC's choices aren't fully random and would probably cause problems in an instance where a random seed makes really bad random picks the first time around.

What do you think, should this be default behaviour or another one for matlab_strict?

In Matlab PREP, the criteria for a "flat" channel are that it either has a median absolute deviation below 10e-10 or a standard deviation below 10e-10.

In PyPREP this is mis-copied over as 1e-10, but there's an additional difference in that MNE and EEGLAB scale their raw data values differently, such that the first 2 samples in MATLAB are

[-3.7068, -4.8012, ...]

whereas in MNE they're

[-3.70680070e-06, -4.80121040e-06, ...].

Unless I'm doing the math wrong in my head, the correct equivalent threshold to use for PyPREP would therefore be 10e-16.

In the full PREP example there is a validation against the Matlab PREP, for which we need the Matlab Results stored here.

That's reasonable, I think.

HOWEVER, the test data stored here could easily be downloaded and not stored in the repo. E.g., via the MNE datasets module.

Once this is done, Manifest.in will have to be adjusted.

I seem to be having an issue when running this pipeline. Below is the traceback:

RuntimeWarning: invalid value encountered in greater
  mask = (size > 0)
Traceback (most recent call last):
  File "/Users/opt/anaconda3/envs/mytfenv/lib/python3.7/site-packages/IPython/core/interactiveshell.py", line 3319, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-4-997c912b5fa5>", line 23, in <module>
    prep_pipeline.fit()
  File "/Users/Documents/preprocessing/prep_pipeline.py", line 80, in fit
    reference.perform_reference()
  File "/Users/Documents/preprocessing/reference.py", line 66, in perform_reference
    self.robust_reference()
  File "/Users/Documents/preprocessing/reference.py", line 144, in robust_reference
    noisy_detector.find_all_bad_channels(ransac=self.ransac)
  File "/Users/Documents/preprocessing/find_noisy_channels.py", line 140, in find_all_bad_channels
    self.find_bad_by_ransac()
  File "/Users/Documents/preprocessing/find_noisy_channels.py", line 476, in find_bad_by_ransac
    n_samples=n_samples,
  File "/Users/Documents/preprocessing/find_noisy_channels.py", line 581, in run_ransac
    chn_pos, chn_pos_good, good_chn_labs, n_pred_predicted, data
  File "/Users/Documents/preprocessing/find_noisy_channels.py", line 632, in get_ransac_pred
    interpol_mat = _make_interpolation_matrix(reconstructed_positions, chn_pos)
  File "/Users/opt/anaconda3/envs/mytfenv/lib/python3.7/site-packages/mne/channels/interpolation.py", line 85, in _make_interpolation_matrix
    C_inv = linalg.pinv(C)
  File "/Users/opt/anaconda3/envs/mytfenv/lib/python3.7/site-packages/scipy/linalg/basic.py", line 1304, in pinv
    a = _asarray_validated(a, check_finite=check_finite)
  File "/Users/opt/anaconda3/envs/mytfenv/lib/python3.7/site-packages/scipy/_lib/_util.py", line 246, in _asarray_validated
    a = toarray(a)
  File "/Users/opt/anaconda3/envs/mytfenv/lib/python3.7/site-packages/numpy/lib/function_base.py", line 499, in asarray_chkfinite
    "array must not contain infs or NaNs")
ValueError: array must not contain infs or NaNs

I understand that ValueError: array must not contain infs or NaNs means that there are some non numerical values in my array. However, when I check whether there are in fact NaNs or infs values in my array using this answer I get False as the result (i.e there are no infs or NaNs). Is there anyway to fix this issue?

In PyPREP, we do a first pass of NoisyChannels prior to initial average referencing to

a) keep track of what the original channels flagged as "bad" were, and

b) figure out which channels are "unusable" for median-referencing before we start the actual NoisyChannels -> interpolate bads -> re-reference loop.

In PyPREP, we use the bad channels from this initial pass of NoisyChannels as the starting state for the subsequent loop (and once a channel is considered 'bad', it's permanently excluded from being a reference channel during re-referencing). See the last line here, where self.noisy_channels gets initialized with a copy of self.noisy_channels_original:

MatPREP, however, only adds the 'unusable' channels (i.e. those bad by NaN, flat, or SNR) to the starting state of permanently bad channels, leaving the deviation, HF noise, correlation, dropout, and RANSAC initial bads empty.

Personally I think MatPREP's behaviour makes more sense here, given that the signal hasn't necessarily been average-referenced coming into PyPREP (e.g. could be referenced to Cz) and the initial average referencing could change the "bad" status for a lot of channels (e.g. if Cz was particularly noisy). In fact, I'd go one further and say that "bad-by-flat" channels from the first NoisyChannels pass shouldn't be permanently excluded either for this same reason: if it's flat because it's a reference channel, it won't be flat anymore after the initial average referencing.

Finally coming back to PyPREP after an extended break (got super-busy with other stuff), and just noticed a weirdness I think might be a bug: during robust re-referencing, after the initial NoisyChannels call, the noisy channel detection loop that's supposed to run a maximum of 4 iterations will actually run a maximum of 5, and requires at least 3 loops to complete before breaking.

The reason for this is that iterations is set to 0 at the start of the loop, but is only incremented at the very end (after the "check if able to break the loop or if too many iterations" part as already happened). I noticed this because for one file, the 'bad_all' for my first two iterations were identical but PyPREP still went on to do a 3rd loop, which seemed weird to me.

Looking at the corresponding MATLAB PREP code, their logic has the same issue: they still allow 5 iterations when the maximum is 4, and require a minimum of 3 RANSAC loops before no change in bad channels can break the loop. However, looking at the pseudocode in the original PREP publication, it doesn't say anything about a minimum number of iterations, just "break from loop if badChannels didn't change or iteration criteria has been met" which suggests that it might be an honest mistake.

My question is: should PyPREP mimic this behaviour? For the purpose of PyPREP vs MATLAB PREP comparison it should probably be left in as a configurable setting (e.g. a matlab_compat flag or something that tells PREP to strictly follow the original logic), but beyond that is there a rationale for this extra iteration I'm not thinking of?

If it's confirmed this is unwanted, I'll put together a PR to fix it!

Dear Pyprep team,

I'm currently using the pyprep pipeline and it has been quite helpful so far. However, I am confused about which of the the attributes in the prep output is the final preprocessed data? Is it EEG_clean or is it EEG_new? In addition, how can the cleaned data be integrated back into the mne data structure?

Thanks for your help,

Sebastian

The Pull Request Template

is written completely in Markdown ... but it contains parts that contributors are supposed to delete, such as

Thanks for contributing. If this is your first time, please make sure to read the contributing guideline.

These parts should be written as an html comment, so that if these parts are accidentally not deleted, they do at least not show up in the opened PR, like here:

Html comment:

 <!--This is a comment. Comments are not displayed in the browser-->

Furthermore, there should be an additional point in the template, asking contributors to please pay attention to the "development install instructions" --> https://github.com/sappelhoff/pyprep#installation

Regarding the validation example. I saw the plot containing the differences across channels and times. It is helpful to see the overall differences in all the channels.

Now I think it would be helpful to see the signals directly to get a more intuitive grasp of how much difference is there. Why this?

Back when I tried to do my own implementation of prep in python I developed some functions that would compare two signals. Following the example they made I see they did a kind of relative difference against the max value of the EEG. When I did mine I divided by the value of the matlab EEG in that point, this caused that when the signal of matlab was near zero I would get big errors, to solve this problems I searched for various methods and found this thread: https://stats.stackexchange.com/questions/86708/how-to-calculate-relative-error-when-the-true-value-is-zero . So I added the relative percent difference and some others.

Point is, I was still not getting a real sense of what the difference was so I made a function that would get the difference vector and automatically plot the signals in the neighborhood of the worst difference. This would help me assess better whether the differences were significant or not.

The original plot implemented in the example that I get when running it:

What follows are the plots I get using my error plotting functions (the data itself is the same as the one used by the previous plot). I use the relative percentage difference (RPD) method from the link I mentioned. I provide for each stage the error matrix image (channels x time) and an example of the signals plotted near the greatest difference points.

NOTICE: RPD value is always between 2 and -2 thus should be interpreted in that range

Raw Signal:

Here the signal is exactly alike so no signals plot here. (Notice the 1e-16 scaling).

New Signal: (Detrend)

Clean Signal (Line Noise Filter)

Ref Signal (Referencing): (maybe related to #19 ?--> Apparently NO, This pipeline was run with a code that already fixed that one)

Interpolated Final Signal :

In this last one the difference is quite high across the channels x time points. Some high differences are expected because prep and pyprep didn't interpolate exactly the same channels, so those channels not agreed by both will have a high difference.

Disclaimer

Depending on how one measures the difference the error matrix image (channels x time) and the choice of the color map can look quite different and one may freak out or not depending on that.

For example in the final stage (interpolation) we have the following plots that doesn't look at scary as the relative percentage difference one, these were obtained by other methods.

Unsigned RPD (0 to 2)

Unsigned Relative to the mean value (this is %):

Notice axis goes up to 700%

Unsigned Relative to the max value (this is %)

Notice axis goes up to 80%

Squared Error divided by max (Similar to the method used in the example)

Nevertheless all refer to the same signals so because of that I think is useful to take a look at the signals themselves to evaluate the significance of the differences.

Currently the full PREP example does both, showing how it can be used, and validating the results against the Matlab results.

https://pyprep.readthedocs.io/en/latest/auto_examples/index.html

I think we should rather have two examples:

• one example where PREP usage is shown more extensively and with more explanations than currently
• and another example with a focus on validation against the Matlab results (here, the plotting code should be improved, because currently it's very messy to read and probably lots of unneeded code duplication)

I've been using PyPREP to set up an EEG preprocessing pipeline for a project I'm working on (~40 minute sessions, 1000Hz, 32 channels), and for most of it I've been using a fixed random seed. However, I got curious about the variability in RANSAC between runs with random seeds, so I went through my full dataset a few times with a script, writing out the names and numbers of initial bad, interpolated, and remaining bad channels per session to a .csv. On some files, I saw differences as big as 4 channels between two consecutive runs. Overall, there was at least one difference in interpolated channels on a majority of the 85 files.

Given that RANSAC is based on thresholds and random sample consensus, I'm not necessarily surprised by this, but I was wondering if there were any parameters I could adjust to make RANSAC more stable between runs (presumably at the cost of speed/RAM)?

we need to refactor much of the code to be in smaller chunks (functions, preferably) that can be properly tested with pytest.

that'd make everyone more comfortable when we change pyprep behavior, like suggested in #19

For whatever reason, MATLAB PREP uses a different method of interpolation for channel prediction during RANSAC than it does for interpolating channels during and after re-referencing.

For RANSAC, they use their own internal interpolation code (a copy of the private interpolation code from EEGLAB's clean_rawdata plugin), which happens to match MNE's _make_interpolation_matrix exactly. For all other interpolation, they use a wrapper around EEGLAB's public interpolation code, which seems to use a different method of spherical spline interpolation (different numeric results, no reference to Perrin '89 in the docstring).

Anyway, I can't think of any a priori justification for using different interpolation logic in one place over the other and my gut is to trust MNE's code more than EEGLAB's, but regardless it's something that needs to be fixed for Reference to produce numerically-identical results to MATLAB PREP for #58. According to the original manuscript it seems like they comparison-tested the methods and found they gave similar correlations with actual signals?

The PREP pipeline uses the spherical option of EEGLAB eeg_interp function for channel interpolation. This function uses Legendre polynomials up through degree 7. To test this choice, we compared this interpolation method with the v4 option and two other spherical interpolation functions, including the one used in RANSAC. We applied a rank sum test on the correlation of the interpolated channel with the actual channel and found no significant difference between the three spherical interpolation methods. The v4 option performed significantly worse. The results were not sensitive to whether a block of channels was removed in the immediate vicinity of the bad channel or whether the interpolation was done in ICA space and projected back.

Sort of confused by "there were no significant differences, but v4 was significantly worse", though. May be worth doing our own testing.

(came first up in #93 (comment) and following comment in that thread.)

This isn't a PyPREP bug so much as a weirdness in the original PREP that PyPREP now matches. Essentially, during robust referencing, PREP performs an initial NoisyChannels pass and flags any bad-by-NaN, bad-by-flat, and bad-by-SNR channels as "unusable", permanently excluding them from being included in the calculated average reference signals:

However, PREP doesn't include initial bad-by-SNR channels in the full set of bad channels to interpolate (see #91), so a channel that's initially bad-by-SNR prior to initial average referencing but is no longer bad afterwards will still get used for interpolating bads during the reference loop, despite being fully excluded from the average reference signal itself. This seems wrong, and like something worth fixing.

My gut instinct here would be to default to only marking initial bad-by-flat and bad-by-NaN channels as fully "unusable". I'd imagine that the original intent was to removing channels that contribute more noise than signal from contaminating the average reference, but given how sensitive the 'bad-by-correlation' default settings are (if over 1% of all correlation windows are below threshold, the whole channel is considered bad) and how that can change after average referencing, I don't think initial bad-by-SNR channels should be treated as permanently bad.

Another difference I caught trying to unit-test my code for #77: during re-referencing, after the initial pass that detects and removes bad-by-NaN/flat/low-SNR channels from the initial average reference, MATLAB PREP's rereferencing loop doesn't interpolate channels that are bad by SNR or dropout (including any bad-by-SNR channels that were flagged on the first pass).

Basically, in PyPREP we always update the list of total bad channels after each re-reference iteration using NoisyChannels.get_bads(), which returns all bad channel types. In MatPREP, it's updated by the following function:

function ref = updateBadChannels(ref, noisy)
% Update the bad channel lists from ref based on bad channels in noisy
    ref.badChannelsFromNaNs = union(ref.badChannelsFromNaNs, ...
        noisy.badChannelsFromNaNs);
    ref.badChannelsFromNoData = union(ref.badChannelsFromNoData, ...
        noisy.badChannelsFromNoData);
    ref.badChannelsFromHFNoise = union(ref.badChannelsFromHFNoise, ...
        noisy.badChannelsFromHFNoise);
    ref.badChannelsFromCorrelation = union(ref.badChannelsFromCorrelation, ...
        noisy.badChannelsFromCorrelation);
    ref.badChannelsFromDeviation = union(ref.badChannelsFromDeviation, ...
        noisy.badChannelsFromDeviation);
    ref.badChannelsFromRansac = union(ref.badChannelsFromRansac, ...
        noisy.badChannelsFromRansac);
    ref.badChannelsFromDropOuts =  union(ref.badChannelsFromDropOuts, ...
        noisy.badChannelsFromDropOuts);
    ref.all = union(...
              union(ref.badChannelsFromNaNs, ref.badChannelsFromNoData), ...
              union( ...
                union(ref.badChannelsFromHFNoise, ref.badChannelsFromCorrelation), ...
              union( ...  
              union(ref.badChannelsFromDeviation,ref.badChannelsFromRansac), ...
                    ref.badChannelsFromRansac)));

As you can see,

1. Bad-by-SNR channels aren't updated with each iteration of re-referencing (make sense, since they're just the union of bad-by-HF and bad-by-correlation channels).
2. The set of all bads is the combination of all updated bad-by-NaN/flat/deviation/HF-noise/correlation/RANSAC channels instead of a union between ref.all and noisy.all. This means that any bad-by-SNR channels flagged during the first pass don't get carried over into the ref.all for future iterations if they aren't flagged bad by some other metric.
3. Dropout channels are updated with each iteration but aren't included in ref.all, possibly by mistake (see: the duplicate bad-by-RANSAC at the end).

Looking back at the MatPREP commit where that function was introduced, it looks like this function replaces simple updating of the full set of noisy channels so at least some of it seems to be deliberate?

Anyway, I'm not entirely sure how to best fix this because I can't easily tell whether PyPREP's logic or MatPREP's makes the most sense, and thus whether this should be fully or partly wrapped in matlab_strict. 3 certainly seems to be undesirable, since you probably wouldn't want a channel with intermittent flat regions being included in your estimated clean reference signal. What do you think?

once #9 and some other fixes that @yjmantilla wanted to do are merged, we can release a new version of pyprep.

I added a .zenodo.json file to the repo in preparation for the next release, which will be archived on Zenodo.

I filled in what I already knew, but for some of you (the authors so far), the affiliation and ORCID is missing.

Please post to this thread so that we can complete the document. Please follow this example when you post here:

I suppose most people would want to use PyPrep as an early preprocessing pipeline to get ride of line noise and get a robust average reference (at least, that is my understanding from the documentation). To do that, it seems that PyPrep should make available a raw object that has been cleaned and re-referenced. My understanding is that it does, as the .raw attribute of the PrepPipeline class. However, in this Raw object, some channels (like EOG channels) are dropped by PyPrep. I don't think that they should be, since they are likely to be necessary for the following artifact rejection steps that will be run after PyPrep.

This is a general tracking issue for the numeric equivalence of PyPREP vs Matlab PREP at each stage of the pipeline. I figure this will be useful for keeping a bird's eye view of what's been tested and what hasn't, as well as what parts have been confirmed to differ.

Guide: items with a checkmark have been tested and confirmed as producing identical results to MATLAB. Items with no checkmark and a corresponding issue # have been tested and confirmed as diverging. Items with no checkmark and no issue number have yet to be formally tested.

NoisyChannels Methods

• Bad by NaN
• Bad by flat (fixed by #60)
• Bad by deviation (fixed by #57)
• Bad by HF noise
• Bad by correlation (fixed by #57)
• Bad by low SNR (fixed by #57)
• Bad by dropout (fixed by #81)
• Bad by RANSAC
  • 'Bad' channel exclusion (fixed by #64)
  • Channel picks (with fixed seed) (fixed by #62)
  • Spherical interpolation algorithm (with identical channel locations, produces identical matrix 🎉)
  • Channel interpolation (fixed by #70)
  • Channel correlation (fixed by #67 and #70)
  • Channel thresholding

Other Functions

• removeTrend (fixed by #71)
• CleanLine (mne.filter.notch_filter) (#82: PyPREP's line noise removal is different, seems to be more aggressive)
• Robust re-referencing
  • Initial average referencing (fixed by #78)
  • Reference removal
  • Bad channel identification (fixed by #78, #89, and #93)
  • Bad channel interpolation (fixed by #96)
  • Final re-referencing (fixed by #92)
• Final bad channel interpolation (fixed by #92 and #96)

Through the last few PRs, I've noticed there's quite a bit of configurable options through PyPREP that aren't actually exposed at the level of the PrepPipeline object (e.g. all of the NoisyChannels threshold parameters). In addition, the arguments/options that are currently exposed are divided between the params dict and actual keyword arguments in a mix of MatPREP's struct-based config style and a more Pythonic approach.

As an improved way of handling things, I was wondering if a PrepSettings object might be a good idea: instead of a difficult-to-document dict or a bunch of duplicated kwargs passed down multiple levels, PyPREP would have a properly-documented class that would initialize defaults, be able to validate parameter values, and expose methods and getters/setters for things to configure.

For example, if you wanted to run the Pipeline with a maximum of 3 robust reference iterations and a relaxed bad-by-deviation threshold of 3.29, on a file with a powerline frequency 50 Hz:

config = PrepSettings(line_freq=50)
config.reference_settings(max_iterations=3)
config.bad_by_deviation_settings(threshold=3.29)

pipeline = PrepPipeline(raw, config, montage)
pipeline.fit()

What do you think, would something like this make for a better API going forward?