Samu validation

test_sph.py

from __future__ import print_function

import sys
from os import path as op
import numpy as np
from scipy.io import loadmat
from mne.io import read_info
from mne import pick_types, pick_info
from mne.preprocessing.maxwell import (_concatenate_sph_coils, _sh_negate,
                                       _bases_real_to_complex, _deg_order_idx,
                                       _sph_to_cart_partials, _sph_harm_norm,
                                       _sph_harm, _alegendre_deriv,
                                       _cart_to_sph, _sss_basis)
from mne.datasets import testing
from mne.forward._make_forward import _prep_meg_channels


###############################################################################
# Validate Sin

mu0 = 1.25664e-6  # Permeability of vacuum


def Sin_vsh_enm(info, origin, magscale, Lin):
    coils = _prep_meg_channels(info, True, [], True, True, True, False)[0]
    Sin = np.empty((len(info['chs']), Lin * (Lin + 2)), np.complex128)
    rs, ws, ezs, bins = _concatenate_sph_coils(coils)
    rs -= origin
    rad, az, pol = _cart_to_sph(rs).T
    for degree in range(1, Lin + 1):
        for order in range(degree + 1):
            sph_norm = _sph_harm_norm(order, degree)
            sph = _sph_harm(order, degree, az, pol, norm=False)
            sph *= sph_norm
            az_factor = 1j * order * sph / np.sin(np.maximum(pol, 1e-16))
            pol_factor = (-sph_norm * np.sin(pol) * np.exp(1j * order * az) *
                          _alegendre_deriv(order, degree, np.cos(pol)))

            in_norm = rad ** -(degree + 2)
            g_rad = in_norm * -(degree + 1) * sph
            g_az = in_norm * az_factor
            g_pol = in_norm * pol_factor
            V = _sph_to_cart_partials(az, pol, g_rad, g_az, g_pol)
            V = np.einsum('ij,ij,i->i', V, ezs, ws)
            x = (np.bincount(bins, V.real, minlength=len(coils)) +
                 1j * np.bincount(bins, V.imag, minlength=len(coils)))
            Sin[:, _deg_order_idx(degree, order)] = -mu0 * x
            if order > 0:
                Sin[:, _deg_order_idx(degree, -order)] = \
                    -mu0 * _sh_negate(x, order)
    Sin[pick_types(info, meg='mag')] *= magscale
    SNin = Sin / np.linalg.norm(Sin, axis=0)[np.newaxis, :]
    return Sin, SNin

data_path = testing.data_path()
fname_mat = op.join(data_path, 'SSS', 'sss_data.mat')
fname_raw = op.join(data_path, 'SSS', 'test_move_anon_raw.fif')
data = loadmat(fname_mat)
origin = np.array([0., 0., 0.04])
info = read_info(fname_raw)
info_meg = pick_info(info, pick_types(info, meg=True, exclude=[]))
print('Validating bases (S_in)...', end='')
sys.stdout.flush()
Sin_mat, SNin_mat = data['Sin0040'], data['SNin0040']
Sin, SNin = Sin_vsh_enm(info_meg, origin, 100., 8)
np.testing.assert_allclose(Sin, Sin_mat, rtol=1e-6, atol=1e-9)
np.testing.assert_allclose(SNin, SNin_mat, rtol=1e-6, atol=1e-9)
print(' [Done]')
sys.stdout.flush()

print('Validating mne bases (S_in)...', end='')
coils = _prep_meg_channels(info_meg, accurate=True, elekta_defs=True,
                           head_frame=True, verbose=False)[0]
Sin, SNin = Sin_vsh_enm(info_meg, origin, 100., 8)
Sin_py = _bases_real_to_complex(_sss_basis(origin, coils, 8, 0, 100.,
                                           method='alternative'), 8, 0)
Sin_py *= mu0
np.testing.assert_allclose(Sin, Sin_py, rtol=1e-6, atol=1e-9)
print(' [Done]')
sys.stdout.flush()