Ransac Comparison Example

[yjmantilla]

PR Description

Basically a PR for the ransac comparison between autoreject's and pyprep's implementation.
This is a first version so is a WIP. All ideas welcome for the comparison.

Currently for the real EEG comparison they differ in two channels:

AR: ['C3', 'F7', 'F8', 'FT8', 'T9', 'P7']
PYPREP: ['C3', 'C1', 'F8', 'FT8', 'P7', 'P4']

The reason why is not entirely clear to me. We should maybe investigate further how autoreject's ransac works and whether the input parameters I give him here really do make sense for the comparison.

I did notice autoreject use an specific way to generate the channels subsets (adds a bias to the seed --> https://github.com/autoreject/autoreject/blob/dd0942438440070876e8ea4796b83ffdd4981600/autoreject/ransac.py#L198 )

We have the correlation matrices of both methods available here for comparison but I have not really used them yet.

Also Autoreject's Ransac is faster at least for these examples. Im not sure yet how is the memory consumption in autoreject.

The doc's building is failing because autoreject is not installed in whatever executes the example. I don't know where that is configured.

In the final version of this PR it should close #36

Merge Checklist

• the PR has been reviewed and all comments are resolved
• all CI checks pass
• (if applicable): the PR description includes the phrase closes #<issue-number> to automatically close an issue
• (if applicable): bug fixes, new features, or API changes are documented in whats_new.rst

[sappelhoff]

looks like a good start to me :)

I recently found out that in sphinx gallery you can replace the "bar" of # with # %% signs. I like that second option better, because then editors like VS Code can pick up on that and let you run the example similar to a notebook in the editor.

see here for more info: https://sphinx-gallery.github.io/stable/syntax.html#embed-rst-in-your-example-python-files

regarding the autoreject import, I added it to requirements-dev.txt --> that's what contains all dependencies that are needed for testing, building the docs, and development in general.

The dependencies for pyprep are in setup.cfg

[sappelhoff]

link to example as built by the CI service: https://pyprep--53.org.readthedocs.build/en/53/auto_examples/ransac_compare.html#sphx-glr-auto-examples-ransac-compare-py

[yjmantilla]

@sappelhoff

So been delving a bit into the workings of both ransac's.

• pyprep calculates 1 epoch less than AR, i fixed that, was just a limit in a np.arangue
• AR selects particular subset of channels for each (epoch,sample,channel), pyprep just does for (sample,channel)
• AR will iterate through epochs finding the correlations for all channels in each iteration
• Pyprep will iterate through channels finding the correlations for all epochs in each iteration

As of now depending on the data pyprep and AR will either consume more memory (I think...,my speculation):

• If there are more channels than epochs --> AR consumes more
• If there are more epochs than channels --> Pyprep consumes more

Memory wise :

• AR will use a (n_points_1epoch,n_channels,n_samples) array
• pyprep will use (n_points_all_epochs,1 channel,n_samples) array

If we wanted to optimize memory consumption it would be something of the sort of iterating through both epochs and channels. Or intelligently select the iteration through the greatest between n_epochs and n_channels
In that sense it may be easier to refactor AR to include the channel iteration because they already made the matlab validation.

Pyprep Ransac results wont exactly match AR until we manage to work also epoch wise for the prediction, given that the subsets of channels change epoch-wise in Autoreject. (edit: not only epoch-wise but also sample-wise which the original MATLAB RANSAC also did --> #41 (comment))

[sappelhoff]

great work! Thanks, that makes a lot of sense!

pyprep calculates 1 epoch less than AR, i fixed that, was just a limit in a np.arangue

did you push that commit yet? I can't see it in the changes

As of now depending on the data pyprep and AR will either consume more memory (I think...,my speculation):

your speculation sounds reasonable 👍

If we wanted to optimize memory consumption it would be something of the sort of iterating through both epochs and channels. Or intelligently select the iteration through the greatest between n_epochs and n_channels

That might be a bit too much optimization, wdyt? the autoreject way already uses quite little memory compared to pyprep's current one.

In that sense it may be easier to refactor AR to include the channel iteration because they already made the matlab validation.

mh yes, it would be totally fine to have one working RANSAC in python, we don't need different implementations in different packages. But the autoreject RANSAC can only be fit on epochs, it would be nice to be able to fit it on continuous data (even if then under the hood, epochs are being calculated).

[codecov-io]

Codecov Report

Merging #53 (0314f94) into master (8476ea5) will increase coverage by 0.03%.
The diff coverage is n/a.

@@            Coverage Diff             @@
##           master      #53      +/-   ##
==========================================
+ Coverage   96.80%   96.84%   +0.03%     
==========================================
  Files           6        7       +1     
  Lines         532      538       +6     
==========================================
+ Hits          515      521       +6     
  Misses         17       17              

Impacted Files	Coverage Δ
pyprep/__init__.py	100.00% <0.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 8476ea5...57043ac. Read the comment docs.

[yjmantilla]

did you push that commit yet? I can't see it in the changes

Just did, had to clean up a bit the files since I had experiments in the branch.

Thinking about it now, I wrote it wrongly. It is not giving "1 less epoch compared to AR", it is giving 1 less epoch compared to mne.make_fixed_length_epochs and I think that just happens when corr_frames divides the data evenly which was the case in the simulated example.

For reference: https://github.com/mne-tools/mne-python/blob/5909e1bb0044e558536e62d441f0d694be65377a/mne/event.py#L906

That might be a bit too much optimization, wdyt? the autoreject way already uses quite little memory compared to pyprep's current one.

I still think it is worthwhile to do. I mean being either of:

• AR (n_points_1epoch,n_channels,n_samples)
• PYPREP (n_points_all_epochs,1 channel,n_samples)

Both will consume a lot of ram and for some users it will matter. Guess one could say is not an urgent feature but it is important nonetheless.

mh yes, it would be totally fine to have one working RANSAC in python, we don't need different implementations in different packages. But the autoreject RANSAC can only be fit on epochs, it would be nice to be able to fit it on continuous data (even if then under the hood, epochs are being calculated).

Yes that is fixed under the hood; was solved by #36 (comment) , mainly:

        ############AUTO REJECT####################
        if mode == 'autoreject':
            ransac = Ransac(picks=good_idx,n_resample=n_samples, min_channels=fraction_good, min_corr=corr_thresh, unbroken_time=fraction_bad, n_jobs=njobs, random_state=435656, verbose= 'tqdm')
            info = mne.create_info(self.ch_names_new.tolist(), self.sample_rate, ch_types='eeg', verbose=None)
            # Option 1: Single Epoch
            #EEGData_epoch = self.EEGData[np.newaxis, :, :]
            #epochs = mne.EpochsArray(EEGData_epoch,info)
            
            # Option 2: Epochs of corr_window_secs seconds as windows
            raw = mne.io.RawArray(self.EEGData,info)
            epochs = mne.make_fixed_length_epochs(raw, duration=corr_window_secs, preload=True,reject_by_annotation=False, verbose=None)

            # Either way...
            epochs.set_montage(self.montage)
            epochs_clean = ransac.fit(epochs)
            #print(epochs_clean.bad_chs_)
            self.bad_by_ransac = [i for i in epochs_clean.bad_chs_]
            return None
        ############AUTO REJECT####################

Disclaimer: Not sure what behaviour it has if corr_window_secs doesnt evenly divide the data, thats worth checking. My guess is that it cuts off the tail.

Im not sure of how Autoreject does the memory checking but I dont think they do one. So for example it may be nice to introduce a verify free ram for them too.

Anyway, what approach do we take? Kill pyprep ransac (so long cowboy...) or keep both and just give the option to run either by a parameter? Do direct pull request to AR if we need changes?

I restate here the advantage of AR being already matlab-validated, that is, we could even offer a parameter for it to wrongly calculate the median if the user is a matlab believer.

[sappelhoff]

Both will consume a lot of ram and for some users it will matter. Guess one could say is not an urgent feature but it is important nonetheless.

I think I'd be convinced when I see the first reasonable memory issue for AR ransac :) ... so far I haven't seen any.

Anyway, what approach do we take? Kill pyprep ransac (so long cowboy...) or keep both and just give the option to run either by a parameter? Do direct pull request to AR if we need changes?

it would probably be more sustainable to kill pyprep ransac, and stay with the wrapped AR ransac within pyprep. The wrapping then allows us to build some functionality on top of AR ransac that would be out of scope in the AR package (such as support for continuous, as opposed to epochs). -- of course, core "enhancements" could then be done as PRs to the AR package.

I restate here the advantage of AR being already matlab-validated, that is, we could even offer a parameter for it to wrongly calculate the median if the user is a matlab believer.

yes, but is there an open comparison of AR-ransac and PREP-ransac? Or is this all buried in Mainak's PhD work? :-)

[yjmantilla]

I think I'd be convinced when I see the first reasonable memory issue for AR ransac :) ... so far I haven't seen any.

Lets see how it goes 😬

yes, but is there an open comparison of AR-ransac and PREP-ransac? Or is this all buried in Mainak's PhD work? :-)

Yeah, i guess thats buried somewhere, although Mainak has left hints along the comments he has made in pyprep's issues.

So what i take from this is that this PR will mutate to:

1 - Implement AR ransac inside pyprep
2 - Example will now be between MATLAB RANSAC and AR RANSAC rather than pyprep's RANSAC.

Is this correct?

[sappelhoff]

Implement AR ransac inside pyprep

yes (=adding autoreject as a dependency), and let's keep pyprep ransac around for some more time --- once we are sure that the AR ransac works fine for our needs, we can remove pyprep ransac.

Example will now be between MATLAB RANSAC and AR RANSAC rather than pyprep's RANSAC.

yes, either that, or comparing all three --> but in the long term then rather only AR ransac 👍

[codecov-commenter]

Codecov Report

Merging #53 (1d7e03f) into master (4614ca0) will not change coverage.
The diff coverage is n/a.

❗ Current head 1d7e03f differs from pull request most recent head 3cbfcf0. Consider uploading reports for the commit 3cbfcf0 to get more accurate results

@@           Coverage Diff           @@
##           master      #53   +/-   ##
=======================================
  Coverage   98.51%   98.51%           
=======================================
  Files           7        7           
  Lines         675      675           
=======================================
  Hits          665      665           
  Misses         10       10           

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4614ca0...3cbfcf0. Read the comment docs.

[sappelhoff]

tests are failing here because I wanted to pass a np.random.RandomState object to autoreject's RANSAC. Opened an issue about that here: autoreject/autoreject#218 has been merged and fixed :)

On a different note: This PR is now very old and A LOT of stuff has happened in the meantime, so my opinion from #53 (comment) is no longer up to date.

Back then I thought it'd be good to start depending on autoreject's RANSAC, modify it to our needs, and part by part remove pyprep's RANSAC.

Now I think we should keep our ransac with all of @a-hurst's fixes.

Still I think it would be nice to have this example of comparing pyprep's and autoreject's RANSAC. So I'd say, let's review / improve this a bit and then merge it. What do you think @a-hurst @yjmantilla ?

One last important note: I merged master into this branch in debef8c and that "removed" a fix that @yjmantilla added and described in #53 (comment)

should we add that back? Could you have a look how that fits in with your recent changes @a-hurst ?

[a-hurst]

Wow, looking at the commit that syncs this up with the current version you can really see the extent of the work that's been done on pyprep in the past few months. Nice!

So I'd say, let's review / improve this a bit and then merge it. What do you think @a-hurst @yjmantilla?

Yep, this definitely seems useful! I still need to read this PR over properly so it may have this already, but it would be great to have a section in the example that compares/contrasts the usage of PyPREP vs AutoReject (i.e., types of noise detection, full file vs per-epoch) so that understanding when/where to use each is clearer to newcomers (I remember getting confused by this myself last year when first setting up my EEG pipeline).

I can contribute a bit to this since I've stared into the depths of PyPREP long enough to roughly know how it works/behaves and I don't mind writing docs, I'd just need to get a handle on how AutoReject works enough to compare/contrast. I also wonder about whether AutoReject and PyPREP might be usable together once PyPREP's final interpolation of bad channels is made optional, i.e. using PyPREP for adaptive line noise removal and robust re-referencing, but using AutoReject for final bad channel detection/interpolation post-reference to avoid unnecessary interpolation of channels that are only bad during a certain period of time.

One last important note: I merged master into this branch in debef8c and that "removed" a fix that @yjmantilla added and described in #53 (comment)

should we add that back? Could you have a look how that fits in with your recent changes @a-hurst ?

Are you saying that this PR also had an implementation of window-wise chunking for RANSAC, that got reverted with my current implementation? I looked through the commit you linked but couldn't find it: ransac.py looks functionally identical to how it did before I started my RANSAC rewriting/refactoring. Can you point me to the lines in question?

[sappelhoff]

I meant specifically this commit: 57043ac

note the correlation_frames + 1

which we don't have in master

pyprep/pyprep/ransac.py

Lines 168 to 175 in 4614ca0

	# Calculate the size (in frames) and count of correlation windows
	correlation_frames = corr_window_secs * sample_rate
	signal_frames = data.shape[1]
	correlation_offsets = np.arange(
	0, (signal_frames - correlation_frames), correlation_frames
	)
	win_size = int(correlation_frames)
	win_count = correlation_offsets.shape[0]

[sappelhoff]

I can contribute a bit to this since I've stared into the depths of PyPREP long enough to roughly know how it works/behaves and I don't mind writing docs, I'd just need to get a handle on how AutoReject works enough to compare/contrast.

pretty sure that @yjmantilla would welcome your help on this - as would I :-)

I also wonder about whether AutoReject and PyPREP might be usable together once PyPREP's final interpolation of bad channels is made optional, i.e. using PyPREP for adaptive line noise removal and robust re-referencing, but using AutoReject for final bad channel detection/interpolation post-reference to avoid unnecessary interpolation of channels that are only bad during a certain period of time.

could be a nice future example / documentation (separate from this thread / PR)

[sappelhoff]

we should add one more code block below to easily show:

• chs bad that are equal between pyprep and AR
• chs bad only AR
• chs bad only pyprep

[a-hurst]

I meant specifically this commit: 57043ac

note the correlation_frames + 1

which we don't have in master

Ahh, I see. Presumably that matches up PyPREP's RANSAC with AutoReject's, but that change would diverge from MATLAB PREP so we'd need to wrap it in matlab_strict if we wanted to add that change.

My understanding was that PREP/PyPREP don't use that last correlation window because it's at the end of the file and would thus be shorter than the rest of the windows. In AR, making use of that last window makes sense since you want all the data per epoch that you can get, but for PyPREP I think it makes sense to ignore it (e.g. the last window might only be 100 ms long, but would be weighted the same as a full 5 second window in determining whether a channel is bad or not). Does that make sense?

[sappelhoff]

ah yes, that makes sense. A feature for -non-matlab-strict pyprep may be to include the final window if it is only X% shorter than the others, where X could be set to something like 20% .... and if it's more than 20% shorter, don't calculate anything based on it.

it would be good to add that explanation within the codebase in some way, so that future devs can easily know the reasoning.

[yjmantilla]

@sappelhoff @a-hurst

Sorry for the late reply. You guys have been really working hard on pyprep and currently I cant keep up with your fast work hehe.

Now I think we should keep our ransac with all of @a-hurst's fixes.

I agree, after all that was the initial idea.

Still I think it would be nice to have this example of comparing pyprep's and autoreject's RANSAC. So I'd say, let's review / improve this a bit and then merge it. What do you think @a-hurst @yjmantilla ?

I agree, although it wouldnt be high priority imo.

One last important note: I merged master into this branch in debef8c and that "removed" a fix that @yjmantilla added and described in #53 (comment)
Ahh, I see. Presumably that matches up PyPREP's RANSAC with AutoReject's, but that change would diverge from MATLAB PREP so we'd need to wrap it in matlab_strict if we wanted to add that change.

Yes,as @a-hurst said, this was a "matching" I did so that the matrix of correlations had the same shape (and presumably the windows in both implementations of pyprep and ar ransac were the same). That way I could numerically compare the matrix of correlations (or whatever matrix got outside of ransac, I dont quite remember).

My understanding was that PREP/PyPREP don't use that last correlation window because it's at the end of the file and would thus be shorter than the rest of the windows. In AR, making use of that last window makes sense since you want all the data per epoch that you can get, but for PyPREP I think it makes sense to ignore it (e.g. the last window might only be 100 ms long, but would be weighted the same as a full 5 second window in determining whether a channel is bad or not). Does that make sense?

Nice logic. Maybe for the comparison of the correlation matrix between Pyprep's and AR's RANSAC one can just use every window but the last on the AR's Ransac side. I wouldnt know at this moment if the windows align exactly u.u , but I think at the time of this PR they did. After all is a matter of checking the window chunking logic and maybe numerically checking they align if they both have that info in an array like for example the 'correlation offsets' one.

it would be good to add that explanation within the codebase in some way, so that future devs can easily know the reasoning.

I agree.