Error when compensation channel missing if compensation grade is applied

[romquentin]

Hi,

I guess this is related to #5133 and #5218

With an epoch from CTF data with compensation grade 3 applied, it is not possible to use the tfr_morlet function:

freqs = np.logspace(*np.log10([4, 50]), num=20)
n_cycles = freqs / 2.  # different number of cycle per frequency
power, itc = tfr_morlet(epochs, freqs=freqs, n_cycles=n_cycles, use_fft=True,
                        return_itc=True, n_jobs=1)

gives:

RuntimeError: Compensation grade 3 has been applied, but compensation channels are missing: [u'BG1-1609', u'BG2-1609', u'BG3-1609', u'G11-1609', u'G12-1609', u'G13-1609', u'G22-1609', u'G23-1609', u'P11-1609', u'Q11-1609', u'Q13-1609', u'R11-1609', u'R12-1609', u'R13-1609', u'R22-1609', u'R23-1609']
Either remove compensation or pick compensation channels

Obviously, the same error occurs if trying manually to remove the compensation channels: epochs.pick_types(ref_meg=False) Is it possible to make this error simply a warning? Thanks a lot!

Romain

[larsoner]

For now a workaround is probably to do:

picks = pick_types(epcohs.info, meg=True, ref_meg=True)
tfr_morlet(..., picks=picks)

Can you see if it works?

If so, we can adjust _prepare_picks to include ref_meg here:

https://github.com/mne-tools/mne-python/blob/master/mne/time_frequency/tfr.py#L2092

It should really use _pick_data_channels

[romquentin]

Thank you! Using picks = pick_types(epochs.info, meg=True, ref_meg=True) works.

Nevertheless, I found very convenient in previous MNE version to be able to manually remove the ref_meg with epochs.pick_types (ref_meg=False) and then X = epochs.get_data() to play with the data for example. Why just a warning is not enough?

Thanks a lot!

[larsoner]

Why just a warning is not enough?

If you go to do forward/inverse computation it will not work.

[romquentin]

Ok. I think it is too much to prevent doing that only for the case of a later forward/inverse model. A warning saying that forward/inverse will not work seems more appropriate for me. In any case, thanks a lot for your help!

[larsoner]

@bloyl WDYT?

[agramfort]

as the output of tfr_morlet should not depend on the comps (I think) I would favor seeing how we can avoid the error. Same should apply for all sensor space data analysis. Unless it causes silent issues on inverse models?

[larsoner]

@agramfort your original suggestion was to avoid storing data without compensation channels:

https://github.com/mne-tools/mne-python/issues/4241#issuecomment-351407141

It sounds like you think we should allow it now? I don't see a big difference between allowing people to pick during tfr_morlet(..., picks=picks) and just using epochs.pick_types(...). Both allow saving an object with info eventually.

The alternative it sounds like people advocate for now would be to allow saving these objects, but discard the comps that are no longer un-applyable, leaving the user's data in a "can't source localize" state?

[bloyl]

My 2cnts is that we should avoid storing data without compensation channels particularly if compensation matrices have been applied.

I agree that default arguments shouldn't throw exceptions so In this particular case I would change the default picks that are used to include ref_meg if comp matrices are present.

so I would change https://github.com/mne-tools/mne-python/blob/55da4e49709ca6d5189cdceffba84c107b0ac517/mne/time_frequency/tfr.py#L2122-L2124

to something like

    if picks is None:
        if info['comps']:
            picks = pick_types(info, meg=True, eeg=True, ref_meg=True,
                               exclude='bads')
        else:
            picks = pick_types(info, meg=True, eeg=True, ref_meg=False,
                               exclude='bads')

[romquentin]

Ok. I understand. Thank you!

[larsoner]

I don't even think we need to triage really, just use _pick_data_channels(..., ref_meg=True) in all cases.

[larsoner]

My 2cnts is that we should avoid storing data without compensation channels particularly if compensation matrices have been applied.

@bloyl I chatted with @agramfort about this a bit, and one idea is to proceed as:

1. Add a tutorial about working with CTF data (or maybe a general one about system-specific considerations?) that mentions that if you want to do source loc, you need to ensure that you keep your reference channels.
2. Change picks defaults if necessary for objects that can be source localized to include reference channels.
3. In forward/inverse code, if compensation is applied (which we can tell from the MEG coil_type) but appropriate ref_meg are not present, raise an error.
4. Because these make it "safe enough" and facilitates source-space analysis, allow picking with ref_meg=False even if comps are applied (this means dropping info['comps']) with just a logger.info rather than raise RuntimeError.
5. Remove in read_epochs / read_info where the coil_type / info['comps'] consistency check is done (this more or less just needs to be moved to the fwd/inv code for step (3)).

I know this is different from what we decided on earlier, but based on the quick rate of post-release bug reports, it seems like this change would be more helpful for facilitating users' data analyses. Thinking about e.g. decoding on CTF data with compensators applied, you might reasonably watn to do it on compensated data, but not include ref_meg in the Epochs object that. WDYT?

One point of disagreement between @agramfort and I is whether ref_meg should be included by default in functions that return non-source-localizable data, such as tfr_morlet. I think they probably should be included for consistency with source-localizable objects (where we will want to have ref_meg=True by default), but @agramfort thinks that they aren't really "data channels" (e.g., for decoding of brain activity) so his vote is no.

@romquentin @bloyl thoughts?

[romquentin]

👍

Thinking about e.g. decoding on CTF data with compensators applied, you might reasonably watn to do it on compensated data, but not include ref_meg in the Epochs object that. WDYT?

Yes, this is the main reason why I thought raising the error even without going to the source space was not convenient.

One point of disagreement between @agramfort and I is whether ref_meg should be included by default in functions that return non-source-localizable data, such as tfr_morlet. I think they probably should be included for consistency with source-localizable objects (where we will want to have ref_meg=True by default), but @agramfort thinks that they aren't really "data channels" (e.g., for decoding of brain activity) so his vote is no.

No strong feeling but I would vote to not add the ref channels by default to non-source-localizable data (like @agramfort).

[bloyl]

@larsoner I'll defer to you and @agramfort on this.

I think the above steps seem fine particularly if they are coupled with an increase in code coverage with CTF data (which is pretty sparse I think). I think the benefit of the hard error is to raise the bug reports that have been raised. So that we find out what parts of the code base potentially cause these problems. So far its been a plotting function and tfr_morlets. That being said if its more trouble then its worth a warning is better than nothing.

I also don't think its always so easy to know what is "non-source-localizable data" for instance frequency domain data potentially like what is returned tfr_morlet is source localizable (http://www.fieldtriptoolbox.org/tutorial/beamformer)

;------------------------------------------------------------------------------- With a somewhat longer view,: Is there no room in the fif info structure to decouple the channels used in the data array from what is stored in info['chs']? If we could maintain the locations of all channels (including bad and dropped channels) I think many of these problems would go away.

[larsoner]

Is there no room in the fif info structure to decouple the channels used in the data array from what is stored in info['chs']? If we could maintain the locations of all channels (including bad and dropped channels) I think many of these problems would go away.

The FIF format does not support this unfortunately, and it would be quite a big / invasive change to do it.

However, I think we can essentially work around this limitation for source imaging. This pipeline should work in principle:

1. apply compensation
2. compute forward on full set of channels
3. pick channels in data (maybe dropping compensation channels)
4. compute inverse

In other words, because fwd is computed on all channels, we can apply the comp operator to the lead fields, then subselect channels.

This makes me think the right approach might be to:

1. always keep info['comps']
2. only check for correspondence between info['comps'] and data channels in the forward object

I think we can make this work. @bloyl do you agree? If so we'll want to roll back some of the changes we've made having to do with dropping comps. I can give it a shot.

[bloyl]

That seems like it should work.

I think the only additional check should be to enforce that the compensation hasn't changed between fwd computation and application/generation of the inverse operator.

Are you also planning to provide methods to project out artifact components from the full forward, both projs and maybe ICA components? Seems like this approach could avoid situations such as when you compute an EOG proj but only want to source localize right hemisphere channels. Currently this generates warnings and potentially a less ideal correction of the fwd model. (Sorry I couldn't find where this happens in the code)

On Fri, May 25, 2018 at 9:29 AM Eric Larson notifications@github.com wrote:

Is there no room in the fif info structure to decouple the channels used in the data array from what is stored in info['chs']? If we could maintain the locations of all channels (including bad and dropped channels) I think many of these problems would go away.

The FIF format does not support this unfortunately, and it would be quite a big / invasive change to do it.

However, I think we can essentially work around this limitation for source imaging. This pipeline should work in principle:

1. apply compensation
2. compute forward on full set of channels
3. pick channels in data (maybe dropping compensation channels)
4. compute inverse

In other words, because fwd is computed on all channels, we can apply the comp operator to the lead fields, then subselect channels.

This makes me think the right approach might be to:

1. always keep info['comps']
2. only check for correspondence between info['comps'] and data channels in the forward object

I think we can make this work. @bloyl https://github.com/bloyl do you agree? If so we'll want to roll back some of the changes we've made having to do with dropping comps. I can give it a shot.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/mne-tools/mne-python/issues/5221#issuecomment-392057909, or mute the thread https://github.com/notifications/unsubscribe-auth/AFU7nS0yOxXn7DqmMeNDiNSmMjPHdCLTks5t2AccgaJpZM4UHrJM .

[larsoner]

Are you also planning to provide methods to project out artifact components from the full forward, both projs and maybe ICA components? Seems like this approach could avoid situations such as when you compute an EOG proj but only want to source localize right hemisphere channels.

I hadn't thought about it, but I think that this would make sense. It would take care of problems like #2310, where we errantly re-apply projectors using a subset of the original channels used to compute the projector.

This would probably take a lot of careful coding and testing to get right, but it should be doable. It seems like the cleanest option.

[larsoner]

In the meantime we could merge #5225 since it has some fixes that will work either way. Then we roll back the info['comps'] removals and move the checks to the forward code.

[bloyl]

That seems like a good way forward. I think we would want to keep some of the additional testing done in #5133 .