plot_compare_label_extraction.py

""" check single trial morphing + time series extraction

The Problem 
-----------
establish equivalence across morphing + label extraction paths

mode : single trial, single trial averaged, evoked

morphing : sample,  fsaverage

Batch-mode
----------

import os
cmd = 'ipython --pylab osx plot_compare_label_morphing_extraction.py {} {}'
methods = ['dSPM', 'sLORETA', 'MNE']
oris = ['normal', None]
[os.system(cmd.format(i, j)) for i in methods for j in oris]
"""
# Author: Denis Engemann <denis.engemann@gmail.com>
#
# License: BSD (3-clause)

import sys
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import zscore
import mne
from mne.datasets import sample
from mne.fiff import Evoked
from mne import fiff
from mne.minimum_norm import (make_inverse_operator, apply_inverse,
                              apply_inverse_epochs)

method = sys.argv[1] if len(sys.argv) > 1 else 'dSPM'
pick_ori = sys.argv[2] if len(sys.argv) > 2 else None
print(method, pick_ori)
if pick_ori == 'None':
    pick_ori = None

data_path = sample.data_path()
fname_fwd_meeg = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif'
fname_fwd_eeg = data_path + '/MEG/sample/sample_audvis-eeg-oct-6-fwd.fif'
fname_cov = data_path + '/MEG/sample/sample_audvis-cov.fif'
subject = 'sample'
subjects_dir = data_path + '/subjects'
spatial_smooth = 7
snr = 3.0
lambda2 = 1.0 / snr ** 2

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'
event_ids = {'AudL': 1, 'AudR': 2, 'VisL': 3, 'VisR': 4}
tmin = -0.2
tmax = 0.5
# xfun = lambda x: zscore(np.abs(x).mean(0))  # xtraction function

def xfun(x, times):
    x = np.abs(x).mean(0)
    baseline = times < 0
    x -= x[baseline].mean(0)[None]
    x /= x[baseline].std(0)[None]
    return x

label_fname = 'meg-social-gaze/custom_tpj-rh.label'
raw = fiff.Raw(raw_fname)
events = mne.read_events(event_fname)

fsaverage_label = [l for l in
                   mne.labels_from_parc('fsaverage', 'PALS_B12_Brodmann',
                                        subjects_dir=subjects_dir)[0]
                   if 'lh' in l.name and '22' in l.name][0]

fsaverage_srcs = mne.read_source_spaces(subjects_dir +
                                        '/fsaverage/bem/fsaverage-ico-5-src.fif')
fsaverage_label.subject = 'fsaverage'
fsaverage_label.values.fill(1.)
subject_label = fsaverage_label.morph(subject_to=subject, 
                                      subjects_dir=subjects_dir)
subject_srcs = mne.read_source_spaces(subjects_dir + 
                                      '/{0}/bem/{0}-oct-6-src.fif'
                                          .format(subject))

include = []  # or stim channels ['STI 014']
picks = fiff.pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True,
                        include=include, exclude='bads')
epochs = mne.Epochs(raw, events, event_ids, tmin, tmax, picks=picks,
                    baseline=(None, 0), reject=dict(eog=150e-6))

epochs.equalize_event_counts(['AudL', 'AudR', 'VisL', 'VisR'], copy=False)

this_epochs = epochs['AudL'][:20]  # grab a few trials
evoked = this_epochs.average()

fwd = mne.read_forward_solution(fname_fwd_meeg, surf_ori=True)
noise_cov = mne.read_cov(fname_cov)

# regularize noise covariance
noise_cov = mne.cov.regularize(noise_cov, evoked.info,
                               mag=0.05, grad=0.05, eeg=0.1, proj=True)

# Restrict forward solution as necessary for MEG
fwd = mne.fiff.pick_types_forward(fwd, meg=True, eeg=False)

# make an M/EEG, MEG-only, and EEG-only inverse operators
info = evoked.info
inv_op = make_inverse_operator(info, fwd, noise_cov,
                               loose=0.2, depth=0.8)

times = epochs.times * 1e3

c1 = 'single trials, morphed before averaging'
stcs_ = apply_inverse_epochs(this_epochs, inv_op, lambda2, method,
                             pick_ori=pick_ori, return_generator=True)

morph_mat = None
vertices_to = [np.arange(10242), np.arange(10242)]
single_trial_fsaverage_before_average = []
for ii, stc in enumerate(stcs_):
   if ii == 0:
       morph_mat = mne.compute_morph_matrix(subject, 
                                           'fsaverage',
                                            stc.vertno,
                                            vertices_to,
                                            smooth=spatial_smooth,
                                            subjects_dir=subjects_dir)
   stc = stc.morph_precomputed('fsaverage', morph_mat=morph_mat,
                               vertices_to=vertices_to)
   single_trial_fsaverage_before_average.append(stc)

mean_stc1 = sum(single_trial_fsaverage_before_average)
mean_stc1._data /= len(single_trial_fsaverage_before_average)
time_course1 = xfun(mean_stc1.in_label(fsaverage_label).data, times)

c2 = 'single trials, un-morphed'
stcs2 = apply_inverse_epochs(this_epochs, inv_op, lambda2, method,
                             pick_ori=pick_ori)
mean_stc2 = sum(stcs2)
mean_stc2._data /= len(stcs2)
time_course2 = xfun(mean_stc2.in_label(subject_label).data, times)

c3 = 'single trials, morphed after averaging'
mean_stc3 = mean_stc2.morph_precomputed('fsaverage', morph_mat=morph_mat,
                                        vertices_to=vertices_to)
time_course3 = xfun(mean_stc3.in_label(fsaverage_label).data, times)

c4 = 'evoked, un-morphed'
evoked = this_epochs.average()
mean_stc4 = apply_inverse(evoked, inv_op, lambda2, method,
                          pick_ori=pick_ori)
time_course4 = xfun(mean_stc4.in_label(subject_label).data, times)

c5 = 'evoked, un-morphed'
mean_stc5 = mean_stc4.morph('fsaverage', n_jobs=4,
                            smooth=spatial_smooth,
                            subjects_dir=subjects_dir)
time_course5 = xfun(mean_stc5.in_label(fsaverage_label).data, times)

fig = plt.figure()
for i in range(1, 6, 1):
    plt.plot(times, eval('time_course%i' % i).T,
             label=eval('c%i' % i))
plt.legend(loc='best', fontsize=10)
plt.xlabel('Time (ms)')
plt.ylabel('Activation ({})'.format(method))
plt.title('pick_ori: {}, method: {}'.format(pick_ori, method))
mne.viz.tight_layout(fig=fig)
plt.show()
plt.savefig('fig-{}-{}.png'.format(pick_ori, method))

Could you tell me a bit more about this code?
mean_stc2 = sum(stcs2) #this seems to be a bunch of stc objects, so I'm surprised sum doesnt error out
mean_stc2._data /= len(stcs2) #how is _data different from .data or rh_data or lh_data?

mean_stc2 ends up being quite different than if I get one label timecourse from each trial and then average that. Any ideas on why that might be?