Strange errors with calculating source estimates in a frequency band

✉️ Mailing List Archive (read-only)
1

Hi everyone!

I have been trying to code a program that would enable to plot user selectable time windows of averaged source estimates in a certain frequency band in real time with FieldTrip Buffer. I'm now trying to simulate that in offline mode with looping through the epochs.get_data() and then building an Epoch from a single epoch with EpochsArray (because that's the way the real time sytem works). I end up with a couple of odd errors using compute_source_psd_epochs depending on how long my time window is. If my time window is shorter than 10 seconds, my error is this:

/usr/lib/python2.7/dist-packages/IPython/utils/py3compat.pyc in execfile(fname, *where)
    173 else:
    174 filename = fname
--> 175 __builtin__.execfile(filename, *where)

/home/silfves1/mne-python/examples/omatfunktiot/jatk_off.py in <module>()
    380
    381
--> 382 stc = compute_source_psd_epochs(epochs=ep, inverse_operator=inverse_operator, lambda2=lambda2, method=method, fmin=15.0, fmax=25.0, pick_ori=None, pca=True, bandwidth=bandwidth, n_jobs=10, return_generator=False)
    383
    384

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/minimum_norm/time_frequency.pyc in compute_source_psd_epochs(epochs, inverse_operator, lambda2, method, fmin, fmax, pick_ori, label, nave, pca, inv_split, bandwidth, adaptive, low_bias, return_generator, n_jobs, verbose, pick_normal)

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/utils.pyc in verbose(function, *args, **kwargs)
    549 return ret
    550 else:
--> 551 ret = function(*args, **kwargs)
    552 return ret
    553

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/minimum_norm/time_frequency.pyc in compute_source_psd_epochs(epochs, inverse_operator, lambda2, method, fmin, fmax, pick_ori, label, nave, pca, inv_split, bandwidth, adaptive, low_bias, return_generator, n_jobs, verbose, pick_normal)
    669 # return a list

    670 stcs = list()
--> 671 for stc in stcs_gen:
    672 stcs.append(stc)
    673

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/minimum_norm/time_frequency.pyc in _compute_source_psd_epochs(epochs, inverse_operator, lambda2, method, fmin, fmax, pick_ori, label, nave, pca, inv_split, bandwidth, adaptive, low_bias, n_jobs, verbose)
    519
    520 dpss, eigvals = dpss_windows(n_times, half_nbw, n_tapers_max,
--> 521 low_bias=low_bias)
    522 n_tapers = len(dpss)
    523

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/time_frequency/multitaper.pyc in dpss_windows(N, half_nbw, Kmax, low_bias, interp_from, interp_kind)
    202 # only calculate the highest Kmax eigenvalues

    203 w = linalg.eigvals_banded(ab, select='i',
--> 204 select_range=(N - Kmax, N - 1))
    205 w = w[::-1]
    206

/usr/lib/python2.7/dist-packages/scipy/linalg/decomp.pyc in eigvals_banded(a_band, lower, overwrite_a_band, select, select_range)
    698 return eig_banded(a_band, lower=lower, eigvals_only=1,
    699 overwrite_a_band=overwrite_a_band, select=select,
--> 700 select_range=select_range)
    701
    702 _double_precision = ['i','l','d']

/usr/lib/python2.7/dist-packages/scipy/linalg/decomp.pyc in eig_banded(a_band, lower, eigvals_only, overwrite_a_band, select, select_range, max_ev)
    490 vl, vu, il, iu = 0.0, 0.0, min(select_range), max(select_range)
    491 if min(il, iu) < 0 or max(il, iu) >= a1.shape[1]:
--> 492 raise ValueError('select_range out of bounds')
    493 max_ev = iu - il + 1
    494 else: # 1, 'v', 'value'

ValueError: select_range out of bounds

However, if my time window is 10 seconds, here's the error then:

error Traceback (most recent call last)
/usr/lib/python2.7/dist-packages/IPython/utils/py3compat.pyc in execfile(fname, *where)
    173 else:
    174 filename = fname
--> 175 __builtin__.execfile(filename, *where)

/home/silfves1/mne-python/examples/omatfunktiot/jatk_off.py in <module>()
    380
    381
--> 382 stc = compute_source_psd_epochs(epochs=ep, inverse_operator=inverse_operator, lambda2=lambda2, method=method, fmin=15.0, fmax=25.0, pick_ori=None, pca=True, bandwidth=bandwidth, n_jobs=10, return_generator=False)
    383
    384

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/minimum_norm/time_frequency.pyc in compute_source_psd_epochs(epochs, inverse_operator, lambda2, method, fmin, fmax, pick_ori, label, nave, pca, inv_split, bandwidth, adaptive, low_bias, return_generator, n_jobs, verbose, pick_normal)

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/utils.pyc in verbose(function, *args, **kwargs)
    549 return ret
    550 else:
--> 551 ret = function(*args, **kwargs)
    552 return ret
    553

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/minimum_norm/time_frequency.pyc in compute_source_psd_epochs(epochs, inverse_operator, lambda2, method, fmin, fmax, pick_ori, label, nave, pca, inv_split, bandwidth, adaptive, low_bias, return_generator, n_jobs, verbose, pick_normal)
    669 # return a list

    670 stcs = list()
--> 671 for stc in stcs_gen:
    672 stcs.append(stc)
    673

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/minimum_norm/time_frequency.pyc in _compute_source_psd_epochs(epochs, inverse_operator, lambda2, method, fmin, fmax, pick_ori, label, nave, pca, inv_split, bandwidth, adaptive, low_bias, n_jobs, verbose)
    545
    546 # compute tapered spectra in sensor space

--> 547 x_mt, freqs = _mt_spectra(data, dpss, sfreq)
    548
    549 if k == 0:

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/time_frequency/multitaper.pyc in _mt_spectra(x, dpss, sfreq, n_fft)
    439 # remove mean (do not use in-place subtraction as it may modify input x)

    440 x = x - np.mean(x, axis=-1)[:, np.newaxis]
--> 441 x_mt = fftpack.fft(x[:, np.newaxis, :] * dpss, n=n_fft)
    442
    443 # only keep positive frequencies

/usr/lib/python2.7/dist-packages/scipy/fftpack/basic.pyc in fft(x, n, axis, overwrite_x)
    213
    214 if axis == -1 or axis == len(tmp.shape) - 1:
--> 215 return work_function(tmp,n,1,0,overwrite_x)
    216
    217 tmp = swapaxes(tmp, axis, -1)

error: failed in converting 1st argument `x' of _fftpack.zrfft to C/Fortran array

I don't know, why these error come up. Do some of you have any ideas to help? Thank you so much!

Best regards

Sebastian
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.nmr.mgh.harvard.edu/pipermail/mne_analysis/attachments/20140731/55d6286b/attachment.html

2

Hi Sebastian,

It seems to me that you are somehow iterating over the epochs and stc
objects incorrectly. You do not need to use EpochsArray because you can
iterate over the Epochs object directly (same with RtEpochs I think). For
stc, if you have a single epoch, can you check if your iterations are
correct?

A small code snippet / gist that reproduces this issue would be great.

Mainak

3

This next snippet is from the offline version I've been testing:

for ii, data in enumerate(epochs.get_data()):

    # make the 2D epoch to 3D

    data = data[newaxis, ...]

    ep = EpochsArray(data=data, info=raw.info, events=events, tmin=0,event_id=event_id, reject=reject)

    stc = compute_source_psd_epochs(epochs=ep, inverse_operator=inverse_operator, lambda2=lambda2,method=method, fmin=15.0, fmax=25.0, pick_ori=None, pca=True, bandwidth=bandwidth, n_jobs=10, return_generator=False)

and the error comes in the last command.

What is the best way to iterate over rt_epochs? rt.epochs.next()? I'm now doing just

for ii, data in enumerate(rt_epochs):

Sebastian

4

What are you using for `events`? It should have the same number of entries
as the number of epochs (1), I suspect it might not. (We should probably
add a check in mne-python that this is the case...)

Eric

5

Also, can you try something like this?

for ii, ep in enumerate(epochs):
    stc = compute_source_psd_epochs(epochs=ep ...)

which would not require EpochsArray.

Mainak

6

Eric:

I have tried to use the events with which I have determined the epochs which has 99 events. That didn't work so I tried to make another events with just one event with this:

          # make the 2D epoch to 3D (because the EpochsArray has to be given 3D data

    ep = ep[newaxis, ...]

    new_event_list = np.empty([1, 3])
    new_event_list[0] = [event_list[ii][0], 0, 10] # put here the measurement recording start time.
    # Make your own filepath to which to make the event
    fname_new_event = 'FILEPATH'

    # Write the event
    write_events(fname_new_event, new_event_list)

    # Finally make the event
    new_events = read_events(fname_new_event)

and from here on I use the EpochsArray and compute_source_psd_epochs.

Mainak:

I tried your suggestion and that doesn't work either:

AttributeError Traceback (most recent call last)
/usr/lib/python2.7/dist-packages/IPython/utils/py3compat.pyc in execfile(fname, *where)
    173 else:
    174 filename = fname
--> 175 __builtin__.execfile(filename, *where)

/home/silfves1/mne-python/examples/omatfunktiot/jatk_off.py in <module>()
    345 # for ii, ep in enumerate(epochs.get_data()):

    346 for ii, ep in enumerate(epochs):
--> 347 stc = compute_source_psd_epochs(epochs=ep, inverse_operator=inverse_operator, lambda2=lambda2, method=method, fmin=15.0, fmax=25.0, pick_ori=None, pca=True, bandwidth=bandwidth, n_jobs=10, return_generator=False)
    348
    349 # make the 2D epoch to 3D

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/minimum_norm/time_frequency.pyc in compute_source_psd_epochs(epochs, inverse_operator, lambda2, method, fmin, fmax, pick_ori, label, nave, pca, inv_split, bandwidth, adaptive, low_bias, return_generator, n_jobs, verbose, pick_normal)

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/utils.pyc in verbose(function, *args, **kwargs)
    549 return ret
    550 else:
--> 551 ret = function(*args, **kwargs)
    552 return ret
    553

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/minimum_norm/time_frequency.pyc in compute_source_psd_epochs(epochs, inverse_operator, lambda2, method, fmin, fmax, pick_ori, label, nave, pca, inv_split, bandwidth, adaptive, low_bias, return_generator, n_jobs, verbose, pick_normal)
    669 # return a list

    670 stcs = list()
--> 671 for stc in stcs_gen:
    672 stcs.append(stc)
    673

/home/silfves1/.local/lib/python2.7/site-packages/mne-0.8.git-py2.7.egg/mne/minimum_norm/time_frequency.pyc in _compute_source_psd_epochs(epochs, inverse_operator, lambda2, method, fmin, fmax, pick_ori, label, nave, pca, inv_split, bandwidth, adaptive, low_bias, n_jobs, verbose)
    483 # Pick the correct channels from the data

    484 #

--> 485 sel = _pick_channels_inverse_operator(epochs.ch_names, inv)
    486 logger.info('Picked %d channels from the data' % len(sel))
    487 logger.info('Computing inverse...')

AttributeError: 'numpy.ndarray' object has no attribute 'ch_names'

because then the compute_source_psd_epochs uses a numpy array as an epoch.

Sebastian

7

I had a look at the script together with Sebastian. This works:

for ii, ep in enumerate(epochs):
    stc = compute_source_psd_epochs(epochs=ep ..., return_generator=False)

but setting return_generator=True does not work. I suspect it returns an
empty generator but maybe others have more insight.

Mainak