mne-tools / mne-connectivity

Connectivity algorithms that leverage the MNE-Python API.
https://mne.tools/mne-connectivity/dev/index.html
BSD 3-Clause "New" or "Revised" License
66 stars 34 forks source link

Error when using mne_connectivity.spectral_connectivity_time #95

Closed SezanMert closed 2 years ago

SezanMert commented 2 years ago

I have a one-minute-long 12-channel intracerebral EEG data. It is a single epoch, single-trial data and has no event markers. It is simply a measurement taken from a sleeping patient for a minute. I wanted to compute the time-resolved coherence of channels pairwise. For the time-frequency decomposition, I used Continuous Morlet Wavelet Transform. I wanted to obtain the coherence for specific frequency bands, so I used foi parameter of this function. I encountered the following error:

100%|██████████|  : 78/78 [00:05<00:00,   13.02it/s]
Traceback (most recent call last):
  File "C:\Users\user1\Documents\ccn\intracerebral\coherence_analysis\codes\main.py", line 22, in <module>
    cohs = mne_connectivity.spectral_connectivity_time(epochs_onemin, names=None, method='coh', indices=None,
  File "C:\Users\user1\Documents\ccn\intracerebral\coherence_analysis\codes\venv\lib\site-packages\mne_connectivity\spectral\time.py", line 223, in spectral_connectivity_time
    conn = EpochSpectroTemporalConnectivity(
  File "C:\Users\user1\Documents\ccn\intracerebral\coherence_analysis\codes\venv\lib\site-packages\mne_connectivity\base.py", line 1036, in __init__
    super(EpochSpectroTemporalConnectivity, self).__init__(
  File "C:\Users\user1\Documents\ccn\intracerebral\coherence_analysis\codes\venv\lib\site-packages\mne_connectivity\base.py", line 943, in __init__
    super(SpectroTemporalConnectivity, self).__init__(
  File "C:\Users\user1\Documents\ccn\intracerebral\coherence_analysis\codes\venv\lib\site-packages\mne_connectivity\base.py", line 414, in __init__
    self._prepare_xarray(data, names=names, indices=indices,
  File "C:\Users\user1\Documents\ccn\intracerebral\coherence_analysis\codes\venv\lib\site-packages\mne_connectivity\base.py", line 495, in _prepare_xarray
    xarray_obj = xr.DataArray(
  File "C:\Users\user1\Documents\ccn\intracerebral\coherence_analysis\codes\venv\lib\site-packages\xarray\core\dataarray.py", line 402, in __init__
    coords, dims = _infer_coords_and_dims(data.shape, coords, dims)
  File "C:\Users\user1\Documents\ccn\intracerebral\coherence_analysis\codes\venv\lib\site-packages\xarray\core\dataarray.py", line 152, in _infer_coords_and_dims
    raise ValueError(
ValueError: conflicting sizes for dimension 'freqs': length 7 on the data but length 63 on coordinate 'freqs'

When I trace back in the error log, and tried to understand how time.py works, I saw that the function computes the spectral connectivity (inside for loop in line 214). However, once it creates an EpochSpectroTemporalConnectivity object from the result (in line 223), it gives the dimension error that is seen at the end of the error log.

I believe the reason for this error is this: In time.py, line 202, variable “conn” is initialized as a zeros numpy array of size n_epochs x n_pairs x n_freqs x len(times). Here, n_freqs is defined as the frequencies of interest. In my case, its length is equal to 7. But, later on, at line 223, to create “EpochSpectroTemporalConnectivity” class object, “freq” parameter given to the object is equal to the frequencies of the time-frequency decomposition. In my case, the length of it is equal to 63. I think, once the connectivity is computed for the frequencies of interest, the frequencies given to the “EpochSpectroTemporalConnectivity” class object should have been the freqeuncies of interest, rather than the decomposition frequencies.

I am using,

I have created a sample edf file to put here a working example. Inside the created edf file, I have put a single frequency sine wave of same number of samples of the original data and I made sure that the same error occurs for both the original and created sample data. And here is my working minimal code:

import mne
import numpy as np
import mne_connectivity

# read the edf file as a raw type
raw_edf = mne.io.read_raw_edf(input_fname='created_sample.edf')

# To obtain Epochs class object without actually dividing my time-series data, I will create an Epochs
# object of duration 60 seconds
epochs_onemin = mne.make_fixed_length_epochs(raw_edf, duration=60, preload=True)

# sampling rate of my data
fs = 200

# frequency bands of interest
fois = np.array([[0.5, 1], [1, 2], [2, 4], [4, 8], [8, 12], [12, 16], [16, 32]])

# frequencies of Continuous Morlet Wavelet Transform
freqs = np.arange(0.5, 32., 0.5)

# compute coherence
cohs = mne_connectivity.spectral_connectivity_time(epochs_onemin, names=None, method='coh', indices=None,
                                                   sfreq=fs, foi=fois, sm_times=0.5, sm_freqs=1, sm_kernel='hanning',
                                                   mode='cwt_morlet', mt_bandwidth=None, freqs=freqs, n_cycles=5)

created_sample.zip

drammock commented 2 years ago

Thank you for the very detailed issue description, sample data, and code to reproduce! Using your code, I see the same thing on my system.

To me it seems like line 202 is OK though. I ran your example interactively in ipython --pdb and (after hitting u several times to reach the right stack level) I see this:

ipdb> p n_freqs
7
ipdb> p conn.shape
(1, 78, 7, 12000)

So it looks like conn has the correct shape. I think rather the problem is line 223, where the xarray container is initialized with conn data, where instead of freqs=freqs it should probably say freqs=f_vec. For me, changing that in the source code makes your example run without error; can you try it locally (and run on your real data, not the created_sample) and see if the results look sensible?

If that works, would you be up for opening a Pull Request that implements this proposed fix?

SezanMert commented 2 years ago

Thank you for your kind reply! I also thought that line 223, EpochSpectroTemporalConnectivity object creation was faulty. And yes, changing freqs=f_vec indeed solved the problem. Thank you for your time. I am going to open a Pull Request now.

adam2392 commented 2 years ago

Closed by #98