apply_inverse_raw() calls _check_reference() even though no EEG channels are used

[ckiefer0]

This bug report is based on this discussion in the MNE group.

I am trying to compute the inverse solution for continuous raw data. I am only interested in MEG channels so I set eeg to False during the computation of the inverse operator. However, apply_inverse_raw() still calls _check_reference().

Unfortunately, the flag custom_ref_applied was set to True for the HCP data (see discussion for details), which caused _check_reference() to raise an error despite me not even using EEG channels.

apply_inverse_raw() should either only call _check_reference() if EEG channels are being used or _check_reference() should check whether EEG channels are being used before raising an error. Personally, I think the second option is better since this seems relevant for other methods calling _check_reference() as well. A solution could look like this:

def _check_reference(inst, ch_names=None):
    """Check for EEG ref."""
    info = inst.info
    if ch_names is not None:
        picks = [ci for ci, ch_name in enumerate(info['ch_names'])
                 if ch_name in ch_names]
        info = pick_info(info, sel=picks)
    if 'eeg' not in info.get_channel_types():
        return
    if _needs_eeg_average_ref_proj(info):
        raise ValueError(
            'EEG average reference (using a projector) is mandatory for '
            'modeling, use the method set_eeg_reference(projection=True)')
    if info.get('custom_ref_applied', False):
        raise ValueError('Custom EEG reference is not allowed for inverse '
                         'modeling.')

Also as a side question, why does apply_inverse_epochs() not call _check_reference()?

[larsoner]

@ckiefer0 I could not reproduce the _check_reference issue for apply_inverse_raw using the testing dataset but I did fix the issue for apply_inverse_epochs where it was not checked at all. Can you make a minimal example using the sample data perhaps? See #11182

[ckiefer0]

@larsoner sorry for the late reply. I was on vacation.

My question is the following: Why would the apply_inverse functions call _check_reference() in the first place if none of the EEG channels are used for the inverse solution?

Here is a short example, which should fail during the computation of the inverse solution for evoked2.

import mne
from mne.datasets import sample
from mne.minimum_norm import make_inverse_operator, apply_inverse

# process MEG data

data_path = sample.data_path()
raw_fname = data_path / 'MEG' / 'sample' / 'sample_audvis_filt-0-40_raw.fif'

raw = mne.io.read_raw_fif(raw_fname, preload=True)  # already has an average reference
events = mne.find_events(raw, stim_channel='STI 014')

event_id = dict(aud_l=1)  # event trigger and conditions
tmin = -0.2  # start of each epoch (200ms before the trigger)
tmax = 0.5  # end of each epoch (500ms after the trigger)
raw.info['bads'] = ['MEG 2443', 'EEG 053']
baseline = (None, 0)  # means from the first instant to t = 0
reject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)

epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                    picks=('meg', 'eog'), baseline=baseline, reject=reject)

# compute noise cov
noise_cov = mne.compute_covariance(epochs, tmax=0., method=['shrunk', 'empirical'], rank=None, verbose=True)

epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                    picks=('meg', 'eog', 'eeg'), preload=True,
                    baseline=baseline, reject=reject)

evoked = epochs.average()
evoked2 = epochs.set_eeg_reference().average()

evoked = evoked.pick('meg')
evoked2 = evoked2.pick('meg')

fname_fwd = data_path / 'MEG' / 'sample' / 'sample_audvis-meg-oct-6-fwd.fif'
fwd = mne.read_forward_solution(fname_fwd)

inv = make_inverse_operator(evoked.info, fwd, noise_cov, loose=0.2, depth=0.8)
del fwd

method = "dSPM"
snr = 3.
lambda2 = 1. / snr ** 2
stc = apply_inverse(evoked, inv, lambda2, method=method, pick_ori=None, verbose=True)
stc2 = apply_inverse(evoked2, inv, lambda2, method=method, pick_ori=None, verbose=True)

[larsoner]

Okay I can replicate with that and am pushing a fix to #11212: https://github.com/mne-tools/mne-python/pull/11212/commits/bcebd3c7f422d8650dcf68dea157093524b7e9c9 (ignore the bem.py change, it's unrelated)

My question is the following: Why would the apply_inverse functions call _check_reference() in the first place if none of the EEG channels are used for the inverse solution?

I think I see what you mean -- that it seems confusing to call _check_reference if no EEG are in the data. At the same time, the logic for "only check for ref if EEG actually exists in these data" has to be in the code somewhere. In principle this logic can be inside the _check_reference function, or outside it. It's a bit simpler to have it inside actually because it ends up being more DRY that way. Otherwise we'll end up duplicating the "if EEG actually exists in these data" check everywhere we want to possibly call _check_reference. You can see I've gone with the "inside the function" approach, which now has this updated line:

    if _electrode_types(info) and info.get('custom_ref_applied', False):

[ckiefer0]

Seems like a good solution!