lcmv_maxpower.py

lcmv_maxpower.py

# 1. Replicate test
# 2. Expand to left and right hemi
# 3. Switch to oct-6

import os.path as op
import numpy as np
from numpy.testing import assert_allclose, assert_array_less
import mne
from mne.beamformer import make_lcmv, apply_lcmv
from mne.io.constants import FIFF
from mne.transforms import apply_trans, invert_transform

s_path = op.join(mne.datasets.testing.data_path(), 'MEG', 'sample')
fname_fwd = op.join(s_path, 'sample_audvis_trunc-meg-eeg-oct-6-fwd.fif')
fname_evoked = op.join(s_path, 'sample_audvis_trunc-ave.fif')
fname_cov = op.join(s_path, 'sample_audvis_trunc-cov.fif')

evoked = mne.read_evokeds(fname_evoked, condition='Left Auditory',
                          baseline=(None, 0))
evoked.pick_types(meg=True, eeg=False, exclude='bads')

noise_cov = mne.pick_channels_cov(mne.read_cov(fname_cov), evoked.ch_names)
fwd = mne.read_forward_solution(fname_fwd)
fwd = mne.pick_channels_forward(fwd, evoked.ch_names)
fwd_fixed = mne.convert_forward_solution(
    fwd, force_fixed=True, surf_ori=True)

# Pick vertices to use: ones that are localized to the correct vertex when used
# together in a covariance matrix that is the product of their subset of the
# gain matrix with itself (determined iteratively starting from all vertices)

# Good for oct-6
use = [0, 7, 8, 9, 27, 28, 38, 42, 44, 51, 61, 67, 73, 76, 86, 110, 112, 113, 121, 134, 195, 198, 239, 249, 257, 272, 294, 297, 307, 353, 383, 426, 430, 542, 564, 584, 604, 673, 919, 981, 1007, 1033, 1090, 1228, 1313, 1403, 1564, 1622, 1651, 1703, 1712, 1731, 1771, 1842, 1881, 1988, 2363, 2398, 2486, 2568, 2581, 2672, 2742, 2885, 2964, 2980, 3104, 3224, 3250, 3330, 3347, 3378, 3397, 3433, 3460, 3475, 3510, 3517, 3529, 3555, 3572, 3579, 3778, 3780, 3826, 3831, 3832, 3867, 3869, 3878, 3887, 3956, 3957, 3958, 3960, 3965, 3967, 3968, 3969, 3972, 3978, 3985, 3994, 3997, 4003, 4019, 4030, 4039, 4073, 4077, 4078, 4089, 4097, 4098, 4102, 4127, 4177, 4192, 4203, 4219, 4239, 4242, 4276, 4306, 4370, 4385, 4449, 4471, 4485, 4521, 4522, 4618, 4716, 4745, 4789, 4816, 4823, 4851, 4879, 4894, 4898, 4930, 4949, 4952, 5074, 5324, 5358, 5489, 5495, 5514, 5550, 5617, 5619, 5732, 5806, 5876, 5968, 6060, 6177, 6193, 6267, 6300, 6323, 6446, 6473, 6491, 6563, 6627, 6677, 6723, 6736, 6829, 6883, 6899, 7067, 7150, 7177, 7298, 7300, 7332, 7421, 7460, 7504, 7610, 7615, 7629, 7793, 7803, 7806, 7808, 7817, 7881, 7888, 7891, 7903]

# Good for oct-4
# use = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 19, 20, 21, 24, 25, 26, 27, 29, 31, 34, 36, 37, 38, 39, 41, 42, 44, 45, 46, 47, 49, 51, 53, 54, 55, 57, 58, 61, 63, 65, 66, 67, 68, 73, 75, 77, 78, 80, 81, 82, 84, 85, 88, 92, 98, 99, 102, 103, 105, 109, 112, 114, 116, 117, 120, 123, 124, 127, 134, 136, 137, 138, 142, 146, 147, 148, 149, 150, 151, 154, 155, 158, 161, 163, 164, 171, 172, 173, 174, 176, 177, 179, 181, 182, 184, 185, 187, 188, 189, 190, 193, 194, 195, 198, 199, 203, 210, 211, 212, 213, 215, 216, 218, 219, 220, 223, 225, 229, 230, 232, 233, 234, 236, 237, 239, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 280, 281, 282, 283, 285, 288, 289, 291, 292, 293, 295, 296, 299, 300, 301, 302, 303, 304, 307, 308, 309, 310, 311, 312, 313, 320, 321, 325, 326, 327, 329, 330, 332, 333, 335, 338, 339, 340, 346, 348, 350, 353, 361, 363, 365, 366, 369, 371, 373, 374, 375, 377, 378, 382, 387, 388, 389, 396, 399, 401, 408, 413, 415, 419, 420, 421, 422, 425, 429, 432, 434, 435, 436, 439, 440, 441, 444, 449, 450, 458, 459, 461, 462, 467, 469, 470, 476, 478, 479, 482, 483, 484, 485, 486, 488, 489, 490, 493, 494, 495]

evoked = mne.EvokedArray(fwd_fixed['sol']['data'].copy()[:, use], evoked.info)
data_cov = noise_cov.copy()
data = fwd_fixed['sol']['data'][:, use] @ fwd_fixed['sol']['data'][:, use].T
data *= 1e-14  # 100 nAm at each source, effectively (1e-18 would be 1 nAm)
# This is rank-deficient, so let's make it actually positive semidefinite
# by regularizing a tiny bit
data.flat[::data.shape[0] + 1] += mne.make_ad_hoc_cov(evoked.info)['data']
# Do our projection
proj, _, _ = mne.io.proj.make_projector(
    data_cov['projs'], data_cov['names'])
data = proj @ data @ proj.T
data_cov['data'][:] = data
assert data_cov['data'].shape[0] == len(noise_cov['names'])
want = np.arange(fwd_fixed['sol']['data'].shape[1])[use]

filters = make_lcmv(evoked.info, fwd, data_cov, reg=0.,
                    noise_cov=noise_cov, pick_ori='max-power',
                    weight_norm='unit-noise-gain-invariant',
                    depth=True)
loc = apply_lcmv(evoked, filters).data
ori = filters['max_power_ori']
loc = np.abs(loc)
# Compute the percentage of sources for which there is no loc bias:
max_idx = np.argmax(loc, axis=0)
perc = (want == max_idx).mean() * 100
assert perc == 100, perc
# Compute the dot products of our forward normals and
assert fwd['coord_frame'] == FIFF.FIFFV_COORD_HEAD
nn = np.concatenate(
    [s['nn'][v] for s, v in zip(fwd['src'], filters['vertices'])])
nn = nn[want]
ori = ori[max_idx]
nn = apply_trans(invert_transform(fwd['mri_head_t']), nn, move=False)
assert_allclose(np.linalg.norm(nn, axis=1), 1, atol=1e-6)
assert_allclose(np.linalg.norm(ori, axis=1), 1, atol=1e-12)
dots = np.abs((nn * ori).sum(-1))
assert_array_less(dots, 1)
assert_array_less(0, dots)
got = np.mean(dots)
assert 0.78 < got < 0.8, got