plot_theiler_shit.py

import numpy as np
import matplotlib.pyplot as plt
import mne
from mne import io
from mne.connectivity import spectral_connectivity
from mne.datasets import sample

# Set parameters
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'
print(__doc__)

data_path = sample.data_path()
subjects_dir = data_path + '/subjects'
fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
fname_raw = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
fname_event = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
fname_cov = data_path + '/MEG/sample/sample_audvis-cov.fif'
label_name_lh = 'Aud-lh'
fname_label_lh = data_path + '/MEG/sample/labels/%s.label' % label_name_lh

fwd_fname = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif'
fwd = mne.read_forward_solution(fwd_fname, surf_ori=True)

event_id, tmin, tmax = 1, -0.2, 0.5
method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)

# Load data
inverse_operator = read_inverse_operator(fname_inv)


# Setup for reading the raw data
raw = io.Raw(raw_fname, preload=True)

# for pick in picks:
# Add a bad channel
raw.info['bads'] += ['MEG 2443']


picks = mne.pick_types(raw.info, meg=True, eeg=False, stim=False, eog=False,
                       exclude='bads')

# rng = np.random.RandomState(42)
# raw._data[picks] = raw._data[rng.permutation(picks)]
# randomize_phase_raw(raw)
# theilerize_raw(raw, random_seed=1, pi_factor=200)


events = mne.read_events(event_fname)

# Pick MEG gradiometers
picks = mne.pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True,
                       exclude='bads')

# Create epochs for the visual condition
event_id, tmin, tmax = 3, -0.2, 0.5
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=(None, 0), reject=dict(grad=4000e-13, mag=5e-12, eog=150e-6), preload=True)

fwd = mne.pick_channels_forward(fwd, include=epochs.ch_names)
fwd = mne.convert_forward_solution(fwd, force_fixed=True)

noise_cov = mne.read_cov(fname_cov)
inverse_operator = make_inverse_operator(info=epochs.info, forward=fwd,
                                         noise_cov=noise_cov, fixed=True,
                                         loose=None, depth=None)

epochs.drop_channels(['EOG 061'])
stc_gen = apply_inverse_epochs(
    epochs,
    inverse_operator=inverse_operator,
    return_generator=True,
    lambda2=1. / 3 ** 2,
    method='MNE',
    pick_ori=None)

rng = np.random.RandomState(42)
new_epochs = epochs.copy()  # XXX hack
for ii, stc in enumerate(stc_gen):
    # new_epochs.append(
    # XXX insert theilerization here
    new_epochs._data[ii] = mne.apply_forward(
            info=epochs.info, fwd=fwd, stc=stc).data


# epochs.average().plot_joint()
new_epochs.average().plot_joint()
stop
# stop
# Compute connectivity for band containing the evoked response.
# We exclude the baseline period
fmin, fmax = 10., 20.
sfreq = raw.info['sfreq']  # the sampling frequency
tmin = 0.0  # exclude the baseline period
con, freqs, times, n_epochs, n_tapers = spectral_connectivity(
    epochs, method='pli', mode='multitaper', sfreq=sfreq, fmin=fmin, fmax=fmax,
    faverage=True, tmin=tmin, mt_adaptive=False, n_jobs=1)

# the epochs contain an EOG channel, which we remove now
ch_names = epochs.ch_names
idx = [ch_names.index(name) for name in ch_names if name.startswith('MEG')]
con = con[idx][:, idx]

# con is a 3D array where the last dimension is size one since we averaged
# over frequencies in a single band. Here we make it 2D
con = con[:, :, 0]

show_con = con + con.T
plt.matshow(show_con, cmap='viridis')
plt.colorbar()
plt.show()