[MRG] More work on the beamformer code #5447

wmvanvliet · 2018-08-16T19:03:48Z

This PR gives the _create_beamformer another refactoring pass to unify the LCMV and DICS flow.

Bugs fixed:

LCMV beamformer when pick_ori='normal', weight_norm='unit_noise_gain'. The existing implementation has a bug where the weights were first normalized to have unit length before picking the normal direction. Instead, the weights should be normalized to unit length after picking the orientation, as is the case when pick_ori='max-power'.
DICS beamformer when inversion='single', pick_ori='max-power', weight_norm='unit_noise_gain'. The computation of the power (in order to pick the orientation with max power) was bugged.
DICS beamformer when inversion='matrix', pick_ori='max-power', weight_norm='unit_noise_gain'. In this case, a computational shortcut is performed, but it was performed incorrectly. See this comment.

New features:

New fancy reg_pinv function, which now lives in utils.py and will even walk your dog.
No more branching between method='lcmv' and method='dics'.
The sign of the max power orientation is now always aligned with the sign of the surface normal (tnx @larsoner)

New features we got for free in the process:

pick_ori='max-power' now also supported when weight_norm=None
weight_norm='nai' now also supported when pick_ori=None and pick_ori='normal'
weight_norm='nai' now also supported for the DICS beamformer
rank parameter for the DICS beamformer to explicitly specify rank
rank='auto' setting to do rank-based inversion as in ENH: Ditching regularization - beamformer with rank-based inversion of cov matrix #5362 (DICS already did this)

In all other cases, the new version gives exactly the same output as the version in the current master branch:

LCMV

pick_ori	weight_norm	new == old?
None	None	True
None	'unit-noise-gain'	True
None	'nai'	Not available in old
'normal'	None	True
'normal'	'unit-noise-gain'	True
'normal'	'nai'	Not available in old
'max-power'	None	Not available in old
'max-power'	'unit-noise-gain'	True
'max-power'	'nai'	True

DICS

pick_ori	weight_norm	inversion	new == old?
None	None	'single'	True
None	None	'matrix'	True
None	'unit-noise-gain'	'single'	True
None	'unit-noise-gain'	'matrix'	True
None	'nai'	'single'	Not available in old
None	'nai'	'matrix'	Not available in old
'normal'	None	'single'	True
'normal'	None	'matrix'	True
'normal'	'unit-noise-gain'	'single'	True
'normal'	'unit-noise-gain'	'matrix'	True
'normal'	'nai'	'single'	Not available in old
'normal'	'nai'	'matrix'	Not available in old
'max-power'	None	'single'	True
'max-power'	None	'matrix'	True
'max-power'	'unit-noise-gain'	'single'	False
'max-power'	'unit-noise-gain'	'matrix'	False
'max-power'	'nai'	'single'	Not available in old
'max-power'	'nai'	'matrix'	Not available in old

larsoner · 2018-08-16T19:05:51Z

@wmvanvliet can you confirm that we should indeed be emitting warnings when using DICS with multiple channel types, e.g. mag+grad, or worse mag+grad+eeg?

wmvanvliet · 2018-08-16T19:14:09Z

for now, yes, but we should just whiten the CSD like we do with the LCMV beamformer. EIther me or @britta-wstnr will probably get to it soon.

larsoner · 2018-08-16T19:20:57Z

Okay cool. If not before 0.17 (not sure when that will be) then we'll want to make a note in the example

wmvanvliet · 2018-08-16T19:21:30Z

I'm being over-enthousiastic in sharing this, unit tests currently fail. I'll look at it tomorrow and fix things :)

codecov · 2018-08-17T10:30:37Z

Codecov Report

Merging #5447 into master will increase coverage by 0.07%.
The diff coverage is 94.53%.

@@            Coverage Diff             @@
##           master    #5447      +/-   ##
==========================================
+ Coverage   88.43%   88.51%   +0.07%     
==========================================
  Files         363      369       +6     
  Lines       68196    69136     +940     
  Branches    11524    11645     +121     
==========================================
+ Hits        60309    61193     +884     
- Misses       5045     5089      +44     
- Partials     2842     2854      +12

larsoner · 2018-08-17T14:01:38Z

tf_dics failure (sometimes):

https://travis-ci.org/mne-tools/mne-python/jobs/417213324#L2690
https://travis-ci.org/mne-tools/mne-python/jobs/417213328#L2724

wmvanvliet · 2018-08-17T17:06:52Z

goody, that looks like fun to debug...

wmvanvliet · 2018-08-22T13:35:44Z

CircleCI is failing because sphynx-gallery refuses to build or something.

wmvanvliet · 2018-08-22T14:32:52Z

@britta-wstnr are you happy with the weight_norm='nai' as it's implemented in this PR?

wmvanvliet · 2018-08-27T12:17:13Z

Ok, tests all pass and this PR is ready for comments. @britta-wstnr WDYT?

larsoner

Other than versionadded, LGTM. I'm trusting @britta-wstnr to check the math

mne/beamformer/_dics.py

britta-wstnr · 2018-08-29T13:32:41Z

@wmvanvliet I was hoping to find some time while at Biomag - but realized that I need some more time to go over the code! So excuse my slowness here, I will give it a look when back!

wmvanvliet · 2018-08-29T20:46:54Z

Enjoy Biomag!

britta-wstnr · 2018-09-11T10:04:27Z

Hi @wmvanvliet - finally had time looking over the code. However, it is basically impossible to compare old and new code snippets for the things that change the output of LCMV or DICS due to the refactoring. Could you please split the PR into 2, probably best code changes first and then refactoring?

agramfort · 2018-09-11T18:49:05Z

@wmvanvliet diff is indeed quite big and it is harder to review if part is refactoring and part is bug fix / change of behavior. Can you allocate some time to help @britta-wstnr review this work?

If you have doubts let me share a quote. "Along you go fast, together we go far" !

wmvanvliet · 2018-09-12T07:26:19Z

@britta-wstnr here are the changes to the LCMV beamformer, in a separate PR: #5511

wmvanvliet · 2018-09-28T14:52:27Z

Rebased after merging #5511

larsoner · 2018-09-28T14:57:31Z

@wmvanvliet in the meantime can you see if you agree wmvanvliet#6 is correct?

wmvanvliet · 2018-09-28T15:00:07Z

I will take a look!

larsoner · 2018-10-04T13:55:34Z

@britta-wstnr this would be great to get into 0.17. Can you take a look in the next week or so?

larsoner · 2018-10-04T13:58:16Z

@wmvanvliet does this close #5362?

larsoner · 2018-10-04T14:03:50Z

... and also @wmvanvliet we have to figure out how to deal with whats_new.rst. Perhaps in the entry you should add a URL link to this PR since you have a full description of the changes. Something like:

API
---
- Modified and unified DICS and LCMV computations, ... . A table listing the backward-incompatible changes from 0.16 to 0.17 can be viewed `here <https://github.com/mne-tools/mne-python/pull/5447>`__. By `Marijn van Vliet`_

larsoner · 2018-11-03T19:50:23Z

mne/beamformer/_lcmv.py

               if k in info['ch_names']])

+    # XXX this could maybe use pca=True to avoid needing to use
+    # _reg_pinv(..., rank=rank) later


Just a minor comment here that eventually we could make compute_beamformer pass on the pca argument (we currently don't publicly expose it), meaning that it would return a non-rank-deficient multplier of shape (rank, 306) instead of (306, 306). It might simplify things at some point. I've also wanted this parameter before when analyzing data.

larsoner · 2018-11-03T20:12:28Z

Okay @wmvanvliet done, feel free to look to make sure it looks like I did not break anything. I touched the beamformer examples so we could ensure that they are not broken:

@wmvanvliet it looks like the DICS example changed, maybe not for the better. It looks like it's also that way on your commit when I try a8f4f1f locally. Do you think it's okay, or is there something to fix?

EDIT: Updated to reflect CircleCI build for 05813cc

wmvanvliet · 2018-11-03T20:30:33Z

I'll have to take a look. Don't merge this yet. Currently busy with kids but I'll try to get to it soon.

mne/beamformer/_compute_beamformer.py

mne/beamformer/_dics.py

wmvanvliet · 2018-11-04T09:32:40Z

@larsoner how about this? rank='full' makes more sense to me than rank=False.

larsoner · 2018-11-04T13:13:56Z

Yeah that makes sense. Travis is red, though, can you look into it? If not I can do it in a few hours

larsoner · 2018-11-04T18:00:28Z

@wmvanvliet so for lcmv it looks like rank=None was being passed to compute_whitener before by default. I have made it so that rank='full' still does this to maintain backward-compatible behavior, even though it's a bit confusing. compute_whitener still gives an array of full size (n_ch, n_ch), though, so in the end the reg will still make this full rank, so I think it's okay.

larsoner · 2018-11-04T18:32:39Z

doc/whats_new.rst

 API
 ~~~

+- The parameter ``rank`` was added to :func:`mne.beamformer.make_dics` and :func:`mne.beamformer.make_lcmv`` with a default of ``'full'`` that will change to ``None`` in 0.18 to auto-compute the rank of the covariance matrix before regularization by `Marijn van Vliet`_


@wmvanvliet I realized that we should probably do this. We need a deprecation cycle for make_lcmv to do this. I assume we also need it for make_dics but if not let me know and I can set the defaults there back to None instead of '' (the internal signal I used for checking deprecation).

The default for make_dics is already None. For the LCMV, I don't like changing the default without some good simulation study backing it up. Both the LCMV and DICS defaults are carefully selected to yield generally stable beamformers.

@britta-wstnr care to weigh in on this?

FYI this commit is outdated, now DICS has no deprecation but yes LCMV still has a deprecation cycle scheduled to go from 'full' to None. I agree it might make sense to wait a release cycle to do this deprecation, though, to be safe. I'll revert it.

Ok, if you do that, I'm happy with everything.

larsoner · 2018-11-04T19:16:03Z

Updated example links above to reflect their current status. @wmvanvliet if you are happy with the deprecation updates I made and the examples, then I think this is finally ready for merge!

mne/beamformer/tests/test_lcmv.py

larsoner · 2018-11-04T20:31:31Z

Okay reverted the deprecation. It seems like it should be okay to have it in there but I agree it's safer not to. And we want to release ASAP, so let's not do the deprecation "full"->None in 0.17.

@wmvanvliet do you want to try to fix the "Singular matrix" vs "complex spectrum" error depending on build type? To me it seems not worth the effort, unless you see that it's obvious and easy, e.g. some default tolerance is bit too high and fixing it makes it work.

wmvanvliet · 2018-11-04T21:06:20Z

Hang on... something changed in the output of the plot_lcmv_beamformer.py example. This shouldn't happen, as the max-power LCMV beamformer should not have changed at all (see table at the top). I'll need to see what exactly caused the change, it may be a bug. Sorry for yet another delay :)

larsoner · 2018-11-04T21:49:08Z

I did change the example to:

plot the signed values
use the max(abs(...)) instead of max(...) to find the STC max value to plot
pass rank=None, but this only changes the estimate changes a tiny bit

It looks okay to me with these differences in mind. (It actually looks better, since (2) makes it so the same vertex is now found in free orientation as in fixed.)

larsoner · 2018-11-05T00:29:33Z

I can confirm that this (reverting) diff:

-    max_idx = np.argmax(np.abs(stc.data[:, time_idx]))
+    max_idx = np.argmax(stc.data[:, time_idx])
... 
-brain = stc.plot(hemi='lh', views='lat', subjects_dir=subjects_dir,
-                 initial_time=0.1, time_unit='s', smoothing_steps=5)
+brain = abs(stc).plot(hemi='lh', views='lat', subjects_dir=subjects_dir,
+                      initial_time=0.1, time_unit='s')

Puts the plot_lcmv_beamformer.py brain image back as:

But I think the one in this PR is more correct (using max(abs(...)); all three max(abs(...)) vertices are the same!) and more informative (plotting signed data so you can see the signed-ness).

wmvanvliet · 2018-11-05T05:18:56Z

ah ok. Thanks for looking into it.

larsoner · 2018-11-05T11:57:43Z

Thanks @wmvanvliet !

mmagnuski · 2018-11-12T14:39:47Z

🎈 🍾 🎉

wmvanvliet force-pushed the cleanup_beamformer branch from 391ff25 to d56849e Compare August 17, 2018 07:57

wmvanvliet force-pushed the cleanup_beamformer branch from 0a7116a to 8d22f79 Compare August 23, 2018 13:39

wmvanvliet changed the title ~~More work on the beamformer code~~ [MRG] More work on the beamformer code Aug 27, 2018

larsoner reviewed Aug 27, 2018

View reviewed changes

mne/beamformer/_dics.py Outdated Show resolved Hide resolved

mne/beamformer/_dics.py Show resolved Hide resolved

wmvanvliet mentioned this pull request Sep 12, 2018

[MRG+1] Fix weight normalization for LCMV beamformer #5511

Merged

larsoner mentioned this pull request Sep 25, 2018

ENH: agenda for beamformer module #5365

Open

15 tasks

wmvanvliet force-pushed the cleanup_beamformer branch 2 times, most recently from 833e64b to efa3c40 Compare September 28, 2018 14:51

larsoner added this to the 0.17 milestone Oct 4, 2018

larsoner reviewed Nov 3, 2018

View reviewed changes

wmvanvliet commented Nov 4, 2018

View reviewed changes

mne/beamformer/_compute_beamformer.py Outdated Show resolved Hide resolved

wmvanvliet commented Nov 4, 2018

View reviewed changes

mne/beamformer/_dics.py Outdated Show resolved Hide resolved

wmvanvliet added 2 commits November 4, 2018 10:31

Change rank=False to rank='full' in _reg_pinv

f1784cc

Fix docstring of default rank for dics

4346639

Actually implement the change from rank=False => rank='full'

85f3ea2

larsoner added 2 commits November 4, 2018 13:01

Really actually implement the change from rank=False => rank='full'

ebb3de5

FIX: Deprecation

2abe2ab

larsoner reviewed Nov 4, 2018

View reviewed changes

larsoner added 3 commits November 4, 2018 13:37

FIX: Undo deprecation

05813cc

FIX: Spelling

308fc0c

FIX: Raise

9b8cad1

wmvanvliet commented Nov 4, 2018

View reviewed changes

mne/beamformer/tests/test_lcmv.py Show resolved Hide resolved

FIX: No deprecation for LCMV yet

c54e407

wmvanvliet merged commit 552b433 into mne-tools:master Nov 5, 2018

massich mentioned this pull request Feb 15, 2019

FIX: check weight_norm parameter #5946

Merged

wmvanvliet deleted the cleanup_beamformer branch October 18, 2021 06:06

Uh oh!

[MRG] More work on the beamformer code #5447

[MRG] More work on the beamformer code #5447

Uh oh!

Conversation

wmvanvliet commented Aug 16, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

LCMV

DICS

Uh oh!

larsoner commented Aug 16, 2018

Uh oh!

wmvanvliet commented Aug 16, 2018

Uh oh!

larsoner commented Aug 16, 2018

Uh oh!

wmvanvliet commented Aug 16, 2018

Uh oh!

codecov bot commented Aug 17, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

larsoner commented Aug 17, 2018

Uh oh!

wmvanvliet commented Aug 17, 2018

Uh oh!

wmvanvliet commented Aug 22, 2018

Uh oh!

wmvanvliet commented Aug 22, 2018

Uh oh!

wmvanvliet commented Aug 27, 2018

Uh oh!

larsoner left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

britta-wstnr commented Aug 29, 2018

Uh oh!

wmvanvliet commented Aug 29, 2018

Uh oh!

britta-wstnr commented Sep 11, 2018

Uh oh!

agramfort commented Sep 11, 2018

Uh oh!

wmvanvliet commented Sep 12, 2018

Uh oh!

wmvanvliet commented Sep 28, 2018

Uh oh!

larsoner commented Sep 28, 2018

Uh oh!

wmvanvliet commented Sep 28, 2018

Uh oh!

larsoner commented Oct 4, 2018

Uh oh!

larsoner commented Oct 4, 2018

Uh oh!

larsoner commented Oct 4, 2018

Uh oh!

larsoner Nov 3, 2018

Choose a reason for hiding this comment

Uh oh!

larsoner commented Nov 3, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

wmvanvliet commented Nov 3, 2018

Uh oh!

Uh oh!

Uh oh!

wmvanvliet commented Nov 4, 2018

Uh oh!

larsoner commented Nov 4, 2018

Uh oh!

larsoner commented Nov 4, 2018

Uh oh!

larsoner Nov 4, 2018

Choose a reason for hiding this comment

Uh oh!

wmvanvliet Nov 4, 2018

Choose a reason for hiding this comment

wmvanvliet commented Aug 16, 2018 •

edited

Loading

codecov bot commented Aug 17, 2018 •

edited

Loading

larsoner commented Nov 3, 2018 •

edited

Loading

wmvanvliet Nov 4, 2018 •

edited

Loading

larsoner commented Nov 4, 2018 •

edited

Loading

larsoner commented Nov 4, 2018 •

edited

Loading