MRG: Add joint projection for SSS

[larsoner]

It's probably not great that N projectors reduce the rank by N for non-SSS'ed data, but generally by < N (often 0) for SSS'ed data. In other words, once SSS makes mag and grad correlated, computing SSP projectors for them separately might not make sense.

This PR adds support for sss='separate' (default, current master behavior) and 'joint' for SSS'ed data. I didn't do anything special about the scaling between mag+grad, but I think it shouldn't matter because of the correlation structure in the data.

@agramfort @dengemann WDYT? You can look at the test I added at the bottom of the diff to see how this should work, or recreated here in fuller form in case you want to try it:

import mne
data_path = mne.datasets.sample.data_path()
raw = mne.io.read_raw_fif(data_path + '/MEG/sample/sample_audvis_raw.fif')
raw.crop(0, 1.0).load_data().pick_types(exclude=())
raw_sss = mne.preprocessing.maxwell_filter(raw, regularize=None)
sss_rank = 80
assert len(raw_sss.info['projs']) == 0
for sss, n_proj, want_rank in (('separate', 6, sss_rank),
                               ('joint', 3, sss_rank - 3)):
    proj = mne.compute_proj_raw(raw_sss, n_grad=3, n_mag=3, sss=sss,
                                verbose='error')
    this_raw = raw_sss.copy().add_proj(proj).apply_proj()
    assert len(this_raw.info['projs']) == n_proj
    sss_proj_rank = this_raw.estimate_rank()
    cov = mne.compute_raw_covariance(this_raw, verbose='error')
    W, ch_names, rank = mne.cov.compute_whitener(cov, this_raw.info,
                                                 return_rank=True)
    assert ch_names == this_raw.ch_names
    assert want_rank == sss_proj_rank == rank  # proper reduction
    if sss == 'joint':
        assert this_raw.info['projs'][0]['data']['col_names'] == ch_names
    else:
        mag_names = ch_names[2::3]
        assert this_raw.info['projs'][3]['data']['col_names'] == mag_names

[larsoner]

Incidentally, I'm not sure the best way to show that this makes much difference in practice. Maybe there will be some difference in the reconstructed signals. I'll think on it.

[codecov]

Codecov Report

Merging #5146 into master will decrease coverage by <.01%.
The diff coverage is 100%.

@@            Coverage Diff            @@
##           master   #5146      +/-   ##
=========================================
- Coverage    88.8%   88.8%   -0.01%     
=========================================
  Files         401     401              
  Lines       72687   72748      +61     
  Branches    12153   12160       +7     
=========================================
+ Hits        64552   64604      +52     
- Misses       5211    5216       +5     
- Partials     2924    2928       +4

[agramfort]

i am sorry I don't understand what this does.

[larsoner]

For data processed with SSS, gives the option to estimate projectors for mag+grad jointly (i.e. like a single channel type) rather than separately.

[larsoner]

... implemented for now for only a single projector calculation function as an example for discussion (if people agree it's useful I'll add it to other projector calculation functions)

[agramfort]

do you think it makes any difference?

[larsoner]

I have the following theoretical arguments:

1. We treat mag+grad jointly in whitening/inversion, so we should do the same thing in projection.
2. We should have rank_after_proj = rank_before_proj - n_proj. Currently, rank_after_proj = rank_before_proj - n_proj currently depends on how closely the projectors of mag and grad end up aligning to the maxwell_filter bases.

After playing with sample with EOG and ECG projectors I haven't yet found a convincing problematic case. With synthetic ones I can show that rank of compute_whitener depends on how closely the mag and grad projectors align to the same SSS basis component, but all this does is slightly change the GFP plots (and for including or excluding this component, does not appear to change the whitening). So maybe it won't matter...?

[dengemann]

@larsoner this brings us closer to the covariance refactory doesn't it :) I think it's a terrific idea!

[dengemann]

@larsoner does this give any advantage in practice apart from speed?

[agramfort]

@dengemann please review then if you have time. my plate is too full to allocate time for this now. My MNE time now is really restricted to blockers of the next release.

[larsoner]

I haven't managed to produce a case where 'separate' (master) behavior is noticeably better or worse, aside from the fact that the rank values that we get are not consistent. This at least does not seem to break anything in terms of the whitening.

[larsoner]

@dengemann do you have time to look / think about this?

[larsoner]

Okay, rebased, added to whats_new.rst and changed the .. versionaddeds from 0.16 to 0.18.

Ready for review/merge from my end.

[larsoner]

@agramfort this is ready for review/merge (see how it fixes tests). The remaining question is do we deprecate meg='separate' default to get to meg='auto' meaning "joint for SSS'ed data, separate for non-SSS'ed data".

[agramfort]

it's hard for me to judge if meg='joint' should become the default with SSS data without more empirical evaluation of the impact.

[larsoner]

At least in terms of the rank numbers, joint does what you'd expect and separate does not. I've shown this now in the tests.

But I'm okay not changing the default for a bit. In that case I'll just update what's new

[agramfort]

let's just open the what's new for now. Sorry for being many too conservative :)

…

[larsoner]

I also updated the naming from meg='joint' to meg='combined' since that's how we have it in our codebase elsewhere (in cov stuff mostly) and it's less confusing with plot_evoked_joint and similar functions

[larsoner]

@agramfort okay to merge once CIs come back happy?

[agramfort]

thx @larsoner