[MRG] Add spectro-spatial decomposition SSD example #7070

dengemann · 2019-11-16T23:24:15Z

This is a first example demonstrating spatial filtering with SSD (Nikulin et al 2011, NIMG), which essentially enhances the oscillatory signal of interest by removing background signal from surrounding frequencies.
The example is autonomous as is, but if we converge on ideas, I'd extend it to have a Transformer API like SPoC or CSP.

cc @agramfort @DavidSabbagh @pierreablin @wmvanvliet @britta-wstnr

codecov · 2019-11-16T23:51:28Z

Codecov Report

Merging #7070 into master will decrease coverage by 0.04%.
The diff coverage is n/a.

@@            Coverage Diff            @@
##           master   #7070      +/-   ##
=========================================
- Coverage   89.74%   89.7%   -0.05%     
=========================================
  Files         442     444       +2     
  Lines       77783   78981    +1198     
  Branches    12620   12674      +54     
=========================================
+ Hits        69807   70850    +1043     
- Misses       5167    5319     +152     
- Partials     2809    2812       +3

dengemann · 2019-11-16T23:57:21Z

Some outputs:

agramfort

thanks @dengemann for taking a stab at this. How do we integrate this in the MNE public API?

examples/decoding/plot_ssd_spatial_filters.py

agramfort · 2019-11-17T15:36:15Z

See you can make it embrace the same API as XDawn. Basically have a transform for arrays and apply for instances. YOu would get filters and patterns using LinearModel https://mne.tools/stable/generated/mne.decoding.LinearModel.html#mne.decoding.LinearModel

…

dengemann · 2019-11-18T06:46:10Z

Another option would be an ICA like api. Allowing for exploring components and handling artifact removal, whitening etc. You just use a different model to get filters. This makes me also think that I should explore the alternative whitened eigenvalue problem to compute SSD. It should help us handling magnetometers and gradiometers.

…

On 17 Nov 2019, at 16:36, Alexandre Gramfort ***@***.***> wrote: See you can make it embrace the same API as XDawn. Basically have a transform for arrays and apply for instances. YOu would get filters and patterns using LinearModel https://mne.tools/stable/generated/mne.decoding.LinearModel.html#mne.decoding.LinearModel > — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

agramfort · 2019-11-18T09:05:14Z

yes ICA API is similar. fit, transform & apply methods with some plot methods

…

dengemann · 2019-11-18T17:34:22Z

@agramfort pushed first API draft, let me know what you think.

Todo:

switch to whitenining implementation for multiple channels
handle EEG vs MEG.
consolidate API
tests

agramfort · 2019-11-18T21:04:26Z

so far so good https://16861-1301584-gh.circle-artifacts.com/0/dev/auto_examples/decoding/plot_ssd_spatial_filters.html#sphx-glr-auto-examples-decoding-plot-ssd-spatial-filters-py !

britta-wstnr

Just caught some typos :)

examples/decoding/plot_ssd_spatial_filters.py

dengemann · 2019-11-19T11:02:41Z

@britta-wstnr @wmvanvliet thank you for your review. Any thoughts on the proposed API?

wmvanvliet · 2019-11-19T12:04:46Z

examples/decoding/plot_ssd_spatial_filters.py

+    filt_params_noise  : dict
+        Filtering for the frequencies of non-interest.
+    estimator : str
+        Which covariance estimator to use. Defaults to "oas". 


needs pointer to documentation on what estimators you can use and what they do.

examples/decoding/plot_ssd_spatial_filters.py

wmvanvliet · 2019-11-19T12:07:29Z

+1 for having a .transform() which takes an array and an .apply() that takes a Raw/Epochs/Evoked.

dengemann · 2019-11-19T12:26:05Z

@wmvanvliet thoughts on apply semantics? It could be used to reconstruct sensor signal from selected components, be it the most or the least oscillatory. We also should have an apply_cov to filter cov matrices.

wmvanvliet · 2019-11-19T12:49:24Z

Exactly. It could just act as a filter: put an Epochs object in, get an Epochs object out that was filtered for maximum oscillatory activity. Put a Covariance object in, get a Covariance object out (no need for an explicit apply_cov)

dengemann · 2019-12-05T14:00:01Z

examples/decoding/plot_ssd_spatial_filters.py

+            cov_signal = mne.compute_raw_covariance(
+                inst_signal, picks=self.picks_, method=self.estimator)
+            cov_noise = mne.compute_raw_covariance(inst_noise, picks=self.picks_,
+                method=self.estimator)


reminder to self: expose flat/reject values.

mmagnuski · 2019-12-11T15:16:50Z

Hi, if I understand it correctly this is the same GED that is used internally in common spatial pattern and that also has multiple other applications for EEG/MEG (see for example Mike X Cohen's publications on GED).
In that case I think it would be good to first have the basic object that does GED on two covariance matrices (and handles patterns, .apply() etc.) and only then use it for more specialized cases like the one here or CSP.

vpeterson

Hi I have been playing around and cheking this fuction implementation. I have some suggestions

examples/decoding/plot_ssd_spatial_filters.py

vpeterson · 2020-08-21T15:50:40Z

examples/decoding/plot_ssd_spatial_filters.py

+        else:
+            raise NotImplementedError()
+
+        self.eigvals_, self.filters_ = eigh(cov_signal.data, cov_noise.data)


Here the filters are sorted by eigvals ascending order, could it be more helpful to haven them already sorted in descending order?

Suggested change

self.eigvals_, self.filters_ = eigh(cov_signal.data, cov_noise.data)

self.eigvals_, self.filters_ = eigh(cov_signal.data, cov_noise.data)

#sort in descencing order

ix = np.argsort(eigvals_)[::-1]

self.eigvals_=eigvals_[ix]

self.filters_ = eigvects_[:, ix]

vpeterson · 2020-08-21T15:50:40Z

examples/decoding/plot_ssd_spatial_filters.py

+        sptial filter.
+    sort_by_spectral_ratio : bool
+        Whether to sort by the spectral ratio. Defaults to True.
+


this input is not used then, since I guess what was the goal of such input variable, I suggest how to use it while transform is called.

vpeterson · 2020-08-21T15:50:40Z

examples/decoding/plot_ssd_spatial_filters.py

+
+        self.eigvals_, self.filters_ = eigh(cov_signal.data, cov_noise.data)
+        self.patterns_ = np.linalg.pinv(self.filters_)
+


here I used the function spectral_ratio_ssd provided in the example within SSD to be called during ssd.transform

Suggested change

def spectral_ratio_ssd(self,ssd_sources):

"""Spectral ratio measure for best n_components selection

See Nikulin 2011, Eq. (24)

"""

psd, freqs = mne.time_frequency.psd_array_welch(

ssd_sources, sfreq=self.sfreq,n_fft=4096)

sig_idx = freq_mask(freqs, *self.freqs_signal)

noise_idx = freq_mask(freqs, *self.freqs_noise)

spec_ratio = psd[:, sig_idx].mean(axis=1) / psd[:, noise_idx].mean(axis=1)

sorter_spec = spec_ratio.argsort()[::-1]

return spec_ratio, sorter_spec

vpeterson · 2020-08-21T15:50:40Z

examples/decoding/plot_ssd_spatial_filters.py

+
+spec_ratio = spectral_ratio_ssd(psd, freqs, freqs_sig, freqs_noise)
+
+sorter = spec_ratio.argsort()


also this guy

Suggested change

sorter = spec_ratio.argsort()

sorter = ssd.sorter_spec

examples/decoding/plot_ssd_spatial_filters.py

dengemann · 2020-08-24T07:45:51Z

Thank you @vpeterson for your review, I shall get back to you during this week. I have marked this PR as a potential candidate for our next release, hoping that we can merge a first version within the next 3 weeks.

larsoner · 2020-09-09T14:25:32Z

Have any time for this @dengemann ? If not let's push to 0.22

dengemann · 2020-09-09T14:40:15Z

Yes I will get to it this week.

…

On 9 Sep 2020, at 16:25, Eric Larson ***@***.***> wrote: Have any time for this @dengemann ? If not let's push to 0.22 — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

dengemann · 2020-09-14T08:47:29Z

@vpeterson I rebased the pull request to allow us to continue working on this. Shall we give it another try to let you do a PR into this one? We can work through the necessary steps. It would be nice to give credit to your work through commits. I'd also have few discussion points regarding the code you have shared with me.

larsoner · 2020-09-14T13:40:17Z

Other than this PR things are looking on track for release tomorrow

dengemann · 2020-11-12T15:54:34Z

should be good now, carried forward your change @vpeterson

larsoner · 2020-11-12T17:00:32Z

Sorry @dengemann @vpeterson , one more rebase necessary because of latest.inc :(

dengemann · 2020-11-12T17:22:49Z

ok looking into that shortly

…

On 12 Nov 2020, at 18:00, Eric Larson ***@***.***> wrote: Sorry @dengemann @vpeterson , one more rebase necessary because of latest.inc :( — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

@vpeterson

cosmits cosmits draft API added option to handle epoched data, sort component with respect to spectral information and implement an apply function to go back to the signal space Big changes: 1. handle both raw and epoched data. 2. added a boolean input which allows for the components to be sorted with respect to the spectral information. 3. added the apply function, to back project into the signal space Minnor changes: 1. The optinal subtract_signal input was deleted, since this should always be done for estimating the noise. 2. Added a band-fill to the picks plot to highlight the freq-band of interest update version, still missing: bib Update plot_ssd_spatial_filters.py style-checker passed, added check_input and filter_data added :footcite: added link to funcs and handle Epochs Update plot_ssd_spatial_filters.py backticks for references integrate work by @vpeterson + bugfixes + enhance API Update examples/decoding/plot_ssd_spatial_filters.py Co-authored-by: Alexandre Gramfort <alexandre.gramfort@m4x.org> typos make info keyword, handle channel types

clean up sorting revert default move SSD to decoding module add SSD work in progress test on raw and simulated data changes at transform and apply in .transform solved the boolean bug when sort_by_spectral_power was set to False. .apply change to .inverse_transform clean pick variable (no usage anymore) picks was being used! back to it update test_ssd.py added true components detection, error between true and estimated mixing matrix, skelearn pipeline and raw vs epoch analysis update version acording to comments solving code style, run pipeline, add channels type update requirements files added white space added changes accroding to pytest fails Update version according to comments version with pytest passed small changes code style added picks documentation add reference changes according to tests typos checked by codespell-error Apply suggestions from code review Co-authored-by: Eric Larson <larson.eric.d@gmail.com> check documentation building avoid usage of private functions change cross-ref and add __hash__ added decoding.SSD to conf.py file

Co-authored-by: Eric Larson <larson.eric.d@gmail.com>

dengemann · 2020-11-12T19:29:37Z

Rebased.

larsoner

One remaining issue with the existing code, otherwise LGTM!

We should probably add it to tutorials/machine-learning/plot_sensors_decoding.py, though, since that is more or less our directory of machine learning techniques. Can you add a section in there? If it's easy to quickly run on sample data and it makes sense to do so feel free, otherwise just add a new section like we have for Xdawn linking to the new example.

doc/changes/latest.inc

dengemann · 2020-11-12T20:38:56Z

We should probably add it to tutorials/machine-learning/plot_sensors_decoding.py, though, since that is more or less our directory of machine learning techniques. Can you add a section in there?

@larsoner the idea is to merge the basics here that provide a solid foundation for doing SSD. There is more work planned and we want for sure a good machine learning example. I would love to do this in several PRs though. Ok if we postpone additional examples?

larsoner · 2020-11-12T20:41:22Z

I can push it quickly if you are both burned out on making changes -- I think for each new ML method it should have a corresponding entry in tutorials/machine-learning/plot_sensors_decoding.py just like we do nowadays for datasets with doc/overview/datasets_index.rst

dengemann · 2020-11-12T20:46:26Z

Not sure if we need to rush with this. But if you think we need this before merge I don't object if you try the code. One could add it as preprocessor to SPoC or CSP. But I'd really like not to rush with this. We still have some time in this release cycle.

larsoner · 2020-11-12T20:49:34Z

One could add it as preprocessor to SPoC or CSP.

Ahh okay I was thinking about it more as a method meant to be used directly like Xdawn. We can add it later when the bigger picture is clearer then if necessary

dengemann · 2020-11-12T20:52:41Z

Ahh okay I was thinking about it more as a method meant to be used directly like Xdawn.

It's actually unsupervised, if you want. It's more like SSP/ICA actually. I see two important use cases, boosting SNR in noisy data and adding specificity to the interpretation of results as you can kick out non-oscillatory components or signal from surrounding bands. Exploring good examples might take some time. If you plug it in as is you may or may not see obvious benefits for prediction performance.

vpeterson · 2020-11-12T21:08:16Z

Ahh okay I was thinking about it more as a method meant to be used directly like Xdawn.

It's actually unsupervised, if you want. It's more like SSP/ICA actually. I see two important use cases, boosting SNR in noisy data and adding specificity to the interpretation of results as you can kick out non-oscillatory components or signal from surrounding bands. Exploring good examples might take some time. If you plug it in as is you may or may not see obvious benefits for prediction performance.

I agree with this. I just add a comment, in the test_ssd.py we also test how SSD can be used before CSP within the makepipeline sklearn structure. So, yeah, future improvement will be to show the usage of SSD in a good example, but it may take some time.

dengemann · 2020-11-12T22:11:47Z

One Azure build is timing out, unrelated failure.

vpeterson · 2020-11-13T11:46:42Z

now all green @dengemann !!

larsoner · 2020-11-13T13:04:17Z

Thanks @vpeterson @dengemann !

agramfort

sorry I was slow on this one

@dengemann can you have a look?

happy to discuss if you want

thanks for this cool feature

agramfort · 2020-11-13T16:11:19Z

examples/decoding/plot_ssd_spatial_filters.py

@@ -0,0 +1,138 @@
+"""
+===========================================================
+Compute Sepctro-Spatial Decomposition (SDD) spatial filters


agramfort · 2020-11-13T16:11:19Z

examples/decoding/plot_ssd_spatial_filters.py

+
+ssd = SSD(info=raw.info,
+          reg='oas',
+          sort_by_spectral_ratio=False,  # True is recommended here.


why do you use False and in the recommended way you say it should be True?

We need it below for the purpose of the example.

agramfort · 2020-11-13T16:11:19Z

mne/decoding/ssd.py

+    n_components : int | None (default None)
+        The number of components to extract from the signal.
+        If n_components is None, no dimensionality reduction is applied, and
+        the transformed data is projected in the whole source space.


source space?

misnomer, indeed. should be fixed.

agramfort · 2020-11-13T16:11:19Z

mne/decoding/ssd.py

+    picks : array of int | None (default None)
+        The indices of good channels.


agramfort · 2020-11-13T16:11:19Z

mne/decoding/ssd.py

+
+    Attributes
+    ----------
+    filters_ : array, shape(n_channels, n_components)


space missing before (

agramfort · 2020-11-13T16:11:19Z

mne/decoding/ssd.py

+        self.freqs_signal = (filt_params_signal['l_freq'],
+                             filt_params_signal['h_freq'])
+        self.freqs_noise = (filt_params_noise['l_freq'],
+                            filt_params_noise['h_freq'])


why storing freqs_signal and freqs_noise as it's already in the filt_params_signal and filt_params_noise attributes?

agramfort · 2020-11-13T16:11:19Z

mne/decoding/ssd.py

+
+        return self
+
+    def transform(self, X, y=None):


please no y in transform

agramfort · 2020-11-13T16:11:19Z

mne/decoding/ssd.py

+            obtained from continuous data or 3D array obtained from epoched
+            data.
+        y : None | array, shape (n_samples,)
+            Used for scikit-learn compatibility.


there is no y in transform in sklearn

agramfort · 2020-11-13T16:11:19Z

mne/decoding/ssd.py

+        raise NotImplementedError('See ssd.inverse_transform()')
+
+    def inverse_transform(self, X):
+        """Remove selected components from the signal.


this is not an inverse transform

if you follow the API of ICA apply will denoise the data and transform will bring you into "source space"

this deviates from our CSP API

agramfort · 2020-11-13T16:11:19Z

mne/decoding/tests/test_ssd.py

+    Data are simulated in the statistical source space, where n=n_components
+    sources contain the peak of interest.
+    """
+    rs = np.random.RandomState(random_state)


we call it rng usually

dengemann changed the title ~~Add spectro-spatial decomposision SSD example~~ Add spectro-spatial decomposition SSD example Nov 16, 2019

agramfort reviewed Nov 17, 2019

View changes

britta-wstnr reviewed Nov 19, 2019

View changes

wmvanvliet reviewed Nov 19, 2019

View changes

examples/decoding/plot_ssd_spatial_filters.py Show resolved Hide resolved

wmvanvliet reviewed Nov 19, 2019

View changes

examples/decoding/plot_ssd_spatial_filters.py Outdated Show resolved Hide resolved

dengemann commented Dec 5, 2019

View changes

vpeterson suggested changes Aug 21, 2020

View changes

dengemann added this to the 0.21 milestone Aug 24, 2020

dengemann force-pushed the ssd branch from 5e6a5b5 to a397ee6 Compare September 14, 2020 08:41

vpeterson mentioned this pull request Sep 15, 2020

added option to handle epoched data dengemann/mne-python#22

Merged

larsoner modified the milestones: 0.21, 0.22 Sep 24, 2020

dengemann and others added 10 commits November 12, 2020 20:26

remove unnecessary lines

267547f

Update names.inc

a142f0a

add class mne.decoding.ssd

10471d5

fix reg param default + naming

abdbbfe

Apply suggestions from code review by @larsoner

31f15af

Co-authored-by: Eric Larson <larson.eric.d@gmail.com>

add work by @vpeterson

d803081

fix import and remove unnecessary private function in test

82e824b

new style for announcing latest changes cc@vpeterson

4c97e95

dengemann force-pushed the ssd branch from b3f6883 to 4c97e95 Compare November 12, 2020 19:28

larsoner reviewed Nov 12, 2020

View changes

doc/changes/latest.inc Outdated Show resolved Hide resolved

push @vpeterson to the top

ccbc68a

larsoner merged commit d495b13 into mne-tools:master Nov 13, 2020
5 checks passed

agramfort reviewed Nov 13, 2020

View changes

dengemann mentioned this pull request Nov 15, 2020

[FIX] some post-merge fixes for SSD #8527

Merged

vpeterson mentioned this pull request Dec 11, 2020

[FIX]: sort by spectral ratio part of fit + returned_filter not used #8648

Merged

-        self.eigvals_, self.filters_ = eigh(cov_signal.data, cov_noise.data)
+        self.eigvals_, self.filters_ = eigh(cov_signal.data, cov_noise.data)
+        #sort in descencing order
+        ix = np.argsort(eigvals_)[::-1]
+        self.eigvals_=eigvals_[ix]
+        self.filters_ = eigvects_[:, ix]


		self.eigvals_, self.filters_ = eigh(cov_signal.data, cov_noise.data)
		self.patterns_ = np.linalg.pinv(self.filters_)

+def spectral_ratio_ssd(self,ssd_sources):
+       """Spectral ratio measure for best n_components selection
+        See Nikulin 2011, Eq. (24)
+        """
+        psd, freqs = mne.time_frequency.psd_array_welch(
+        ssd_sources, sfreq=self.sfreq,n_fft=4096)
+        sig_idx = freq_mask(freqs, *self.freqs_signal)
+        noise_idx = freq_mask(freqs, *self.freqs_noise)
+        spec_ratio = psd[:, sig_idx].mean(axis=1) / psd[:, noise_idx].mean(axis=1)
+        sorter_spec = spec_ratio.argsort()[::-1]
+        return spec_ratio, sorter_spec


		spec_ratio = spectral_ratio_ssd(psd, freqs, freqs_sig, freqs_noise)

		sorter = spec_ratio.argsort()

		picks : array of int \| None (default None)
		The indices of good channels.

[MRG] Add spectro-spatial decomposition SSD example #7070

[MRG] Add spectro-spatial decomposition SSD example #7070

Conversation

dengemann commented Nov 16, 2019

codecov bot commented Nov 16, 2019 • edited

Codecov Report

dengemann commented Nov 16, 2019

agramfort left a comment

Choose a reason for hiding this comment

agramfort commented Nov 17, 2019 via email

dengemann commented Nov 18, 2019 via email

agramfort commented Nov 18, 2019 via email

dengemann commented Nov 18, 2019

agramfort commented Nov 18, 2019

britta-wstnr left a comment

Choose a reason for hiding this comment

dengemann commented Nov 19, 2019

Choose a reason for hiding this comment

wmvanvliet commented Nov 19, 2019

dengemann commented Nov 19, 2019

wmvanvliet commented Nov 19, 2019

Choose a reason for hiding this comment

mmagnuski commented Dec 11, 2019

vpeterson left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dengemann commented Aug 24, 2020

larsoner commented Sep 9, 2020

dengemann commented Sep 9, 2020 via email

dengemann commented Sep 14, 2020

larsoner commented Sep 14, 2020

dengemann commented Nov 12, 2020

larsoner commented Nov 12, 2020

dengemann commented Nov 12, 2020 via email

dengemann commented Nov 12, 2020

larsoner left a comment

Choose a reason for hiding this comment

dengemann commented Nov 12, 2020

larsoner commented Nov 12, 2020

dengemann commented Nov 12, 2020

larsoner commented Nov 12, 2020

dengemann commented Nov 12, 2020

vpeterson commented Nov 12, 2020

dengemann commented Nov 12, 2020

vpeterson commented Nov 13, 2020

larsoner commented Nov 13, 2020

agramfort left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov bot commented Nov 16, 2019 •

edited