Day 27: Depth weighting and optimization

Today I'm trying to improve upon yesterday's results. I've added in the concept of depth-weighting, which weights source
contributions according to their depth. My strategy is to apply different depth weighting values for each method and try
to find an optimal set of parameters. Right now I'm using a simple dot product to score sources. I then add the dot
products for all the cortical sources and subtract them from the dot product of the actual sources.

Oddly enough, the smallest depth-weighting parameters yield the best results. Tomorrow I'll try a wider range of
parameters.

In [1]:

%load spm_compare_models.py

In []:

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import mne
from mne.datasets import spm_face
from surfer import Brain, TimeViewer

mne.set_log_level(False)

# set data files
data_path = spm_face.data_path()
labels_path = data_path + '/MEG/spm/labels/%s.label'
fstring = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces1_3D'
fwd_fname = fstring + '_oct-6-amyg-vol-fwd.fif'
evo_fname = fstring + '-ave.fif'
cov_fname = fstring + '-cov.fif'

# read evoked data
evo = mne.read_evokeds(evo_fname, condition=0)
sfreq = evo.info['sfreq']
times = evo.times

# read the cortical parcellation for left and right VLPFC
label_names = ['parstriangularis-lh', 'parstriangularis-rh']
labels = [mne.read_label(labels_path % ln) for ln in label_names]

# read the forward model
fwd = mne.read_forward_solution(fwd_fname, force_fixed=False)
fwd = mne.convert_forward_solution(fwd, surf_ori=True)

# extract the source space
src = fwd['src']

# read noise covariance
cov = mne.read_cov(cov_fname)

# generate two waveforms
stc_data = np.zeros((2, len(times)))
ws = mne.time_frequency.morlet(sfreq, [3, 10], n_cycles=[1, 1.5])
stc_data[0][:len(ws[0])] = np.real(ws[0])  # right amygdala
stc_data[1][:len(ws[1])] = np.real(ws[1])  # left amygdala
stc_data[0] = np.roll(stc_data[0], 100)
stc_data[1] = np.roll(stc_data[1], 200)
stc_data *= 100 * 1e-9

# generate simulated sparse source estimates
v = [np.array([]), np.array([]), np.array([0., 3000.])]
sim_stc = mne.source_estimate.SourceEstimate(stc_data, vertices=v,
                                             tmin=times.min(), tstep=1./sfreq)

# set parameters for inverse operator
snr = 10
lambda2 = 1. / snr**2

# generate simulated evoked data
gain = fwd['sol']['data'][:, (0, 3000)]
sim_data = np.dot(gain, sim_stc.data)
sim = evo.copy()
sim.nave = 1
sim.data = sim_data
sim.info = mne.forward._fill_measurement_info(sim.info, fwd, sfreq)

high_score = -100

for method in ['MNE', 'dSPM', 'sLORETA']:

    for depth in np.arange(0.01, 1, 0.01):

        # make the inverse operator
        inv = mne.minimum_norm.make_inverse_operator(evo.info, fwd, cov, loose=None, 
                                                     depth=depth, fixed=False)
        
        # apply the inverse operator
        stc = mne.minimum_norm.apply_inverse(sim, inv, lambda2, method=method)

        # compute the dot product of source data with input signal
        d = np.dot(stc.data, stc_data.T)
    
        score = d[4098, 0] - d[4098, 1] + d[11098, 1] - d[11098, 0] - d[:4098].sum()

        if score > high_score:
            high_score = score
            print method, depth, score   

In [2]:

'MNE 0.01 -3.44334180329e-14'

Out[2]:

'MNE 0.01 -3.44334180329e-14'