# -*- coding: utf-8 -*- # Author: Olaf Hauk # Alexandre Gramfort # Eric Larson # # License: BSD-3-Clause import os.path as op import numpy as np from numpy.testing import (assert_equal, assert_array_almost_equal, assert_array_equal, assert_allclose) import mne from mne.datasets import testing from mne.minimum_norm.resolution_matrix import (make_inverse_resolution_matrix, get_cross_talk, get_point_spread, _vertices_for_get_psf_ctf) data_path = testing.data_path(download=False) subjects_dir = op.join(data_path, 'subjects') fname_inv = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-meg-eeg-oct-6-meg-inv.fif') fname_evoked = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-ave.fif') fname_raw = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc_raw.fif') fname_t1 = op.join(data_path, 'subjects', 'sample', 'mri', 'T1.mgz') fname_src = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-meg-eeg-oct-6-fwd.fif') fname_src_fs = op.join(data_path, 'subjects', 'fsaverage', 'bem', 'fsaverage-ico-5-src.fif') fname_src_3 = op.join(data_path, 'subjects', 'sample', 'bem', 'sample-oct-4-src.fif') fname_stc = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-meg') fname_vol = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-grad-vol-7-fwd-sensmap-vol.w') fname_vsrc = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-meg-vol-7-fwd.fif') fname_inv_vol = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-meg-vol-7-meg-inv.fif') rng = np.random.RandomState(0) fname_fwd = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-meg-eeg-oct-4-fwd.fif') fname_cov = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-cov.fif') fname_label = op.join(data_path, 'subjects', 'sample', 'label', 'lh.V1.label') @testing.requires_testing_data def test_resolution_matrix(): """Test make_inverse_resolution_matrix() function.""" # read forward solution forward = mne.read_forward_solution(fname_fwd) # forward operator with fixed source orientations forward_fxd = mne.convert_forward_solution(forward, surf_ori=True, force_fixed=True) # noise covariance matrix noise_cov = mne.read_cov(fname_cov) # evoked data for info evoked = mne.read_evokeds(fname_evoked, 0) # make inverse operator from forward solution # free source orientation inverse_operator = mne.minimum_norm.make_inverse_operator( info=evoked.info, forward=forward, noise_cov=noise_cov, loose=1., depth=None) # fixed source orientation inverse_operator_fxd = mne.minimum_norm.make_inverse_operator( info=evoked.info, forward=forward, noise_cov=noise_cov, loose=0., depth=None, fixed=True) # regularisation parameter based on SNR snr = 3.0 lambda2 = 1.0 / snr ** 2 # resolution matrices for free source orientation # compute resolution matrix for MNE with free source orientations rm_mne_free = make_inverse_resolution_matrix(forward, inverse_operator, method='MNE', lambda2=lambda2) # compute resolution matrix for MNE, fwd fixed and inv free rm_mne_fxdfree = make_inverse_resolution_matrix(forward_fxd, inverse_operator, method='MNE', lambda2=lambda2) # resolution matrices for fixed source orientation # compute resolution matrix for MNE rm_mne = make_inverse_resolution_matrix(forward_fxd, inverse_operator_fxd, method='MNE', lambda2=lambda2) # compute resolution matrix for sLORETA rm_lor = make_inverse_resolution_matrix(forward_fxd, inverse_operator_fxd, method='sLORETA', lambda2=lambda2) # rectify resolution matrix for sLORETA before determining maxima rm_lor_abs = np.abs(rm_lor) # get maxima per column maxidxs = rm_lor_abs.argmax(axis=0) # create array with the expected stepwise increase in maximum indices goodidxs = np.arange(0, len(maxidxs), 1) # Tests # Does sLORETA have zero dipole localization error for columns/PSFs? assert_array_equal(maxidxs, goodidxs) # MNE resolution matrices symmetric? assert_array_almost_equal(rm_mne, rm_mne.T) assert_array_almost_equal(rm_mne_free, rm_mne_free.T) # Some arbitrary vertex numbers idx = [1, 100, 400] # check various summary and normalisation options for mode in [None, 'sum', 'mean', 'maxval', 'maxnorm', 'pca']: n_comps = [1, 3] if mode in [None, 'sum', 'mean']: n_comps = [1] for n_comp in n_comps: for norm in [None, 'max', 'norm', True]: stc_psf = get_point_spread( rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp, norm=norm, return_pca_vars=False) stc_ctf = get_cross_talk( rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp, norm=norm, return_pca_vars=False) # for MNE, PSF/CTFs for same vertices should be the same assert_array_almost_equal(stc_psf.data, stc_ctf.data) # check SVD variances stc_psf, s_vars_psf = get_point_spread( rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp, norm='norm', return_pca_vars=True) stc_ctf, s_vars_ctf = get_cross_talk( rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp, norm='norm', return_pca_vars=True) assert_array_almost_equal(s_vars_psf, s_vars_ctf) # variances for SVD components should be ordered assert s_vars_psf[0] > s_vars_psf[1] > s_vars_psf[2] # all variances should sum up to 100 assert_allclose(s_vars_psf.sum(), 100.) # Test application of free inv to fixed fwd assert_equal(rm_mne_fxdfree.shape, (3 * rm_mne.shape[0], rm_mne.shape[0])) # Test PSF/CTF for labels label = mne.read_label(fname_label) # must be list of Label label = [label] label2 = 2 * label # get relevant vertices in source space verts = _vertices_for_get_psf_ctf(label, forward_fxd['src'])[0] stc_psf_label = get_point_spread(rm_mne, forward_fxd['src'], label, norm='max') # for list of indices stc_psf_idx = get_point_spread(rm_mne, forward_fxd['src'], verts, norm='max') stc_ctf_label = get_cross_talk(rm_mne, forward_fxd['src'], label, norm='max') # For MNE, PSF and CTF for same vertices should be the same assert_array_almost_equal(stc_psf_label.data, stc_ctf_label.data) # test multiple labels stc_psf_label2 = get_point_spread(rm_mne, forward_fxd['src'], label2, norm='max') m, n = stc_psf_label.data.shape assert_array_equal( stc_psf_label.data, stc_psf_label2[0].data) assert_array_equal( stc_psf_label.data, stc_psf_label2[1].data) assert_array_equal( stc_psf_label.data, stc_psf_idx.data)