plot_theiler_shit.py

import numpy as np

from scipy import linalg
import matplotlib.pyplot as plt
import mne
from mne import io
from mne.connectivity import spectral_connectivity
from mne.datasets import sample
from swish.surrogates import theilerize_raw
from mne.minimum_norm import (apply_inverse, apply_inverse_epochs,
read_inverse_operator)
from mne.connectivity import seed_target_indices, spectral_connectivity


###############################################################################
# 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'
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='grad', 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='grad', 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, eog=150e-6))

stc_gen = apply_inverse_epochs(
    epochs, inverse_operator=inverse_operator, return_generator=True)

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, forward=fwd, stc=stc).data

    
# epochs.average().plot_joint()
# 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()