csd_linearity_fail.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
CSD computation linearity test (which fails).
"""

from copy import deepcopy
from numpy.testing import assert_allclose
import numpy as np
import mne
from mne.time_frequency import csd_fourier, csd_morlet, csd_multitaper

data = np.random.RandomState(0).randn(10, 5, 1000)  # epochs, chs, time
data_orig = data.copy()
epochs = mne.EpochsArray(data, mne.create_info(5, 1000., 'grad'))
epochs.apply_baseline((None, 0))
proj = mne.compute_proj_epochs(epochs)
freq = 10.
for csd_method in (csd_fourier, csd_morlet, csd_multitaper):
    if csd_method is csd_morlet:
        kwargs = dict(frequencies=[freq])
    else:
        kwargs = dict(fmin=freq - 0.1, fmax=freq + 0.1)
    data_csd = csd_method(epochs, **kwargs).get_data(freq)
    epochs_proj = epochs.copy().add_proj(deepcopy(proj)).apply_proj()
    proj_op = mne.proj.make_projector(proj, epochs.ch_names, [], True)[0]
    assert_allclose(epochs_proj._data,
                    np.einsum('oc,ect->eot', proj_op, epochs_proj._data))
    assert_allclose(epochs_proj._projector, proj_op)
    assert_allclose(epochs_proj._data,
                    np.einsum('oc,ect->eot', proj_op, epochs._data))
    data_proj_csd = csd_method(epochs_proj, **kwargs).get_data(freq)
    assert not np.allclose(data_csd, data_proj_csd, atol=1e-7, rtol=1e-2)
    data_csd_proj = np.dot(np.dot(proj_op, data_csd), proj_op.T)
    assert_allclose(data_proj_csd, data_csd_proj, atol=1e-7, rtol=1e-7)