mne-tools / mne-python

MNE: Magnetoencephalography (MEG) and Electroencephalography (EEG) in Python
https://mne.tools
BSD 3-Clause "New" or "Revised" License
2.7k stars 1.31k forks source link

connectivity.spectral_connectivity drops times when given Epochs object #8824

Closed apadee closed 3 years ago

apadee commented 3 years ago

Describe the bug

I'm trying to compute EEG functional connectivity in various time intervals relative to the stimulus.
I run into a problem with the tmin and tmax arguments in mne.connectivity.spectral_connectivity.

If the data object passed to mne.connectivity.spectral_connectivity is of type mne.Epochs, then spectral._prepare_connectivity should use the times from mne.Epochs.times: Otherwise, times is generated using np.linspace() from the sampling frequency sfreq and the number of data samples (n_times_in).

This is done in _prepare_connectivity with a subset of the data, which itself tries to create the times_in through _get_and_verify_data_sizes.

Here, I want to get from -0.4 to -0.1 seconds from the event (300 ms span 100 ms before event.) Here is an evoked plot with the full time axis with the marked interval in which I want to compute connectivity.

image

However, the time points from epochs.times are lost in the execution of spectral._prepare_connectivity and replaced with generated time axis. This makes the time points provided in tmin and tmax no longer consistent with the time in the data object.

A possible fix to this issue is to pass times_in to both spectral._prepare_connectivity and spectral._get_and_verify_data_sizes, with which I then get my expected results. Something like: Definition of _prepare_connectivity :

def _prepare_connectivity(epoch_block, tmin, tmax, fmin, fmax, sfreq, indices,
                          mode, fskip, n_bands,
                          cwt_freqs, faverage, times_in=None):
    """Check and precompute dimensions of results data."""
    first_epoch = epoch_block[0]

    # get the data size and time scale
    n_signals, n_times_in, times_in = _get_and_verify_data_sizes(first_epoch, times=times_in)
...

Calling _prepare_connectivity:

...
(n_cons, times, n_times, times_in, n_times_in, tmin_idx,
             tmax_idx, n_freqs, freq_mask, freqs, freqs_bands, freq_idx_bands,
             n_signals, indices_use) = _prepare_connectivity(
                epoch_block=epoch_block, tmin=tmin, tmax=tmax, fmin=fmin,
                fmax=fmax, sfreq=sfreq, indices=indices, mode=mode,
                fskip=fskip, n_bands=n_bands,
                cwt_freqs=cwt_freqs, faverage=faverage, times_in=times_in)
...

Steps to reproduce

import mne
from mne import io
from mne.connectivity import spectral_connectivity
from mne.datasets import sample
from mne.viz import plot_connectivity_circle

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

raw = io.read_raw_fif(raw_fname)
events = mne.read_events(event_fname)
raw.info['bads'] += ['MEG 2443']
picks = mne.pick_types(raw.info, meg='grad', eeg=False, stim=False, eog=True,
                       exclude='bads')

event_id, tmin, tmax = 3, -1, 1.5
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6))

fmin, fmax = 10., 12.
sfreq = raw.info['sfreq']

tmin = -0.5
tmax = -0.1
epochs.load_data().pick_types(meg='grad')
con, freqs, times, n_epochs, n_tapers = spectral_connectivity(
    epochs, method='pli', mode='multitaper', sfreq=sfreq, fmin=fmin, fmax=fmax,
    faverage=True, tmin=tmin, tmax=tmax, mt_adaptive=False, n_jobs=1)

plot_connectivity_circle(con[:,:,0], epochs.ch_names, n_lines=300)

Expected results

I would expect the times from the epochs object to be used, rather than time points generated from 0.

Actual results

Connectivity computation...
Traceback (most recent call last):
  File "/home/ap/PycharmProjects/eeg_processing/conn_example.py", line 30, in <module>
    faverage=True, tmin=tmin, tmax=tmax, mt_adaptive=False, n_jobs=1)
  File "<decorator-gen-387>", line 21, in spectral_connectivity
  File "/home/ap/miniconda3/envs/MNE/lib/python3.7/site-packages/mne/connectivity/spectral.py", line 789, in spectral_connectivity
    cwt_freqs=cwt_freqs, faverage=faverage)
  File "/home/ap/miniconda3/envs/MNE/lib/python3.7/site-packages/mne/connectivity/spectral.py", line 948, in _prepare_connectivity
    mask = _time_mask(times_in, tmin, tmax, sfreq=sfreq)
  File "/home/ap/miniconda3/envs/MNE/lib/python3.7/site-packages/mne/utils/numerics.py", line 499, in _time_mask
    % (orig_tmin, orig_tmax, extra, times[0], times[-1]))
ValueError: No samples remain when using tmin=-0.5 and tmax=-0.1 (original time bounds are [0.0, 2.4974401644798485])

With a bit of guidance, I'm happy to make PR to resolve the issue :)

Additional information

mne.sys_info()
Backend Qt5Agg is interactive backend. Turning interactive mode on.
Platform:      Linux-5.0.0-32-generic-x86_64-with-debian-buster-sid
Python:        3.7.9 (default, Aug 31 2020, 12:42:55)  [GCC 7.3.0]
Executable:    /home/ap/miniconda3/envs/MNE2/bin/python
CPU:           x86_64: 12 cores
Memory:        Unavailable (requires "psutil" package)
mne:           0.22.0
numpy:         1.19.2 {blas=mkl_rt, lapack=mkl_rt}
scipy:         1.5.2
matplotlib:    3.3.2 {backend=Qt5Agg}
sklearn:       0.23.2
numba:         Not found
nibabel:       Not found
nilearn:       Not found
dipy:          Not found
cupy:          Not found
pandas:        1.2.1
mayavi:        Not found
pyvista:       Not found
vtk:           Not found
larsoner commented 3 years ago

I have not looked at the code but what you describe sounds like a bug. I would start by coding some simple test that fails on master -- for example a two-epoch EpochsArray where -0.5 to 0 has some signal that is coherent and 0 to 0.5 does not -- and show that it fails there. Then when you write your code you already have the test in place that should pass once things are fixed. A PR for this would be great!