lcmv_sim.py

import os.path as op

import numpy as np
from scipy import linalg
from numpy.testing import assert_array_equal, assert_allclose
from mayavi import mlab

import mne
from mne.beamformer import make_lcmv, apply_lcmv
from mne.simulation import simulate_evoked
from mne.minimum_norm import make_inverse_operator, apply_inverse

data_path = mne.datasets.testing.data_path()
fname_raw = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc_raw.fif')
fname_fwd = op.join(data_path, 'MEG', 'sample',
                    'sample_audvis_trunc-meg-eeg-oct-4-fwd.fif')

info = mne.io.read_raw_fif(fname_raw).info
# For speed and for rank-deficiency calculation simplicity,
# just use grads:
info = mne.pick_info(info, mne.pick_types(info, meg='grad', exclude=()))
info.update(bads=[], projs=[])
forward = mne.read_forward_solution(fname_fwd)
forward = mne.pick_channels_forward(forward, info['ch_names'])
vertices = [s['vertno'][::100] for s in forward['src']]
n_vertices = sum(len(v) for v in vertices)
assert 5 < n_vertices < 20
amplitude = 100e-9
stc = mne.SourceEstimate(amplitude * np.eye(n_vertices), vertices,
                         0, 1. / info['sfreq'])
forward_sim = mne.convert_forward_solution(forward, force_fixed=True,
                                           use_cps=True, copy=True)
forward_sim = mne.forward.restrict_forward_to_stc(forward_sim, stc)
noise_cov = mne.make_ad_hoc_cov(info)
noise_cov.update(data=np.diag(noise_cov['data']), diag=False)
evoked = simulate_evoked(forward_sim, stc, info, noise_cov, nave=1)
source_nn = forward_sim['source_nn']
source_rr = forward_sim['source_rr']
# Figure out our indices
mask = np.concatenate([np.in1d(s['vertno'], v)
                       for s, v in zip(forward['src'], vertices)])
mapping = np.where(mask)[0]
assert_array_equal(source_rr, forward['source_rr'][mapping])
# Don't check NN because we didn't rotate to surf ori
del forward_sim

#
# Let's do minimum norm as a sanity check (dipole_fit is slower)
#
inv = make_inverse_operator(info, forward, noise_cov, loose=1.)
stc_vector_mne = apply_inverse(evoked, inv, pick_ori='vector')
mne_ori = stc_vector_mne.data[mapping, :, np.arange(n_vertices)]
subjects_dir = mne.datasets.sample.data_path() + '/subjects'
mne_ori /= np.linalg.norm(mne_ori, axis=-1, keepdims=True)
mne_angles = np.rad2deg(np.arccos(np.sum(mne_ori * source_nn, axis=-1)))
assert np.mean(mne_angles) < 40

ti = np.argmin(mne_angles)
kwargs = dict(smoothing_steps=2, subjects_dir=subjects_dir,
              overlay_alpha=1., hemi='both', views='frontal')
fmax = stc_vector_mne.magnitude().data[:, 5].max()
fmid = fmax / 2.
brain_mne = stc_vector_mne.plot(figure=1, initial_time=evoked.times[5],
                                clim=dict(kind='values', lims=(0, fmid, fmax)),
                                time_label='MNE %0.3f', **kwargs)
if ti < len(vertices[0]):
    hemi = 'lh'
    vert = vertices[0][[ti]]
else:
    hemi = 'rh'
    vert = vertices[1][[ti - len(vertices[0])]]
brain_mne.add_foci(vert, hemi=hemi, coords_as_verts=True)
mlab.view(90, 90, roll=180)

#
# Now let's do LCMV
#
data_cov = mne.make_ad_hoc_cov(info)  # just a stub for later
lcmv_ori = list()

this_evoked = evoked.copy().crop(evoked.times[ti], evoked.times[ti])
data_cov['data'] = (np.outer(this_evoked.data, this_evoked.data) +
                    noise_cov['data'])
vals = linalg.svdvals(data_cov['data'])
assert vals[0] / vals[-1] < 1e5  # not rank deficient
filters = make_lcmv(info, forward, data_cov, 0.05, noise_cov)
filters_vector = make_lcmv(info, forward, data_cov, 0.05, noise_cov,
                           pick_ori='vector')
stc = apply_lcmv(this_evoked, filters)
assert isinstance(stc, mne.SourceEstimate)
stc_vector = apply_lcmv(this_evoked, filters_vector)
assert isinstance(stc_vector, mne.VectorSourceEstimate)
assert_allclose(stc.data, stc_vector.magnitude().data)
fmax = stc.data.max()
fmid = fmax / 2.
brain_lcmv = stc_vector.plot(
    figure=2, scale_factor=1e2, time_label='LCMV %0.3f',
    clim=dict(kind='values', lims=(0, fmid, fmax)), **kwargs)
brain_lcmv.add_foci(vert, hemi=hemi, coords_as_verts=True)
mlab.view(90, 90, roll=180)