compute data somato

import os
import mne
from mne.datasets import somato

data_path = somato.data_path()

raw_fname = data_path + '/MEG/somato/sef_raw_sss.fif'
trans = data_path + '/MEG/somato/sef_raw_sss-trans.fif'
src = data_path + '/subjects/somato/bem/somato-oct-6-src.fif'
bem = data_path + '/subjects/somato/bem/somato-5120-bem-sol.fif'


fwd_fname = data_path + '/MEG/somato/somato-meg-oct-6-fwd.fif'
if os.path.isfile(fwd_fname):
    fwd = mne.read_forward_solution(fwd_fname)

else:
    fwd = mne.make_forward_solution(raw_fname, trans, src, bem,
                                    meg='grad', eeg=False, mindist=5., overwrite=True)
    mne.write_forward_solution(fwd_fname, fwd)

raw = mne.io.read_raw_fif(raw_fname)
events = mne.find_events(raw, stim_channel='STI 014')
reject = dict(grad=4000e-13, eog=150e-6)
picks = mne.pick_types(raw.info, meg='grad', eog=True)

event_id, tmin, tmax = 1, -1., 3.
baseline = (None, 0)
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=baseline, reject=reject,
                    preload=True)
evoked = epochs.average()
evoked = evoked.pick_types(meg=True)
cov = mne.compute_covariance(epochs, tmax=0.)

mne.write_evokeds(data_path + '/MEG/somato/sef-ave.fif', evoked)
mne.write_cov(data_path + '/MEG/somato/sef-cov.fif', cov)

fig fixdgap_l21

from os.path import join as pjoin
import numpy as np
from numpy.linalg import norm
import time
import matplotlib.pyplot as plt
import mne
from mne.datasets import sample, somato
from mne.inverse_sparse.mxne_inverse import (_to_fixed_ori, _prepare_gain)
from mne.inverse_sparse.mxne_optim import mixed_norm_solver, norm_l2inf, norm_l21
# from mne.inverse_sparse.mxne_optim import _mixed_norm_solver_cd
from pygit2 import Repository

fig, ax = plt.subplots(figsize=(7, 3.7))
times_dict = dict()

for data in ["sample", "somato"]:
    if data == "sample":
        data_path = sample.data_path()
        fwd_fname = pjoin(data_path, "MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif")
        ave_fname = data_path + '/MEG/sample/sample_audvis-ave.fif'
        cov_fname = data_path + '/MEG/sample/sample_audvis-shrunk-cov.fif'
        condition = 'Left Auditory'
    elif data == "somato":
        data_path = somato.data_path()
        fwd_fname = data_path + '/MEG/somato/somato-meg-oct-6-fwd.fif'
        ave_fname = data_path + '/MEG/somato/sef-ave.fif'
        cov_fname = data_path + '/MEG/somato/sef-cov.fif'
        condition = 'Unknown'

    fwd = mne.read_forward_solution(fwd_fname, surf_ori=True)
    fwd = _to_fixed_ori(fwd)
    noise_cov = mne.read_cov(cov_fname)
    evoked = mne.read_evokeds(ave_fname, condition=condition, baseline=(None, 0))
    if data == 'sample':
        evoked.crop(tmin=0.04, tmax=0.18)
    else:
        evoked.crop(tmin=0.008, tmax=0.25)
    evoked = evoked.pick_types(eeg=True, meg=True)
    loose, depth = 0, 0.8

    # apply solver on whitened data
    gain, gain_info, whitener, source_weighting, mask = _prepare_gain(
        fwd, evoked.info, noise_cov, pca=False, depth=depth,
        loose=loose, weights=None, weights_min=None)

    sel = [evoked.ch_names.index(name) for name in gain_info['ch_names']]
    M = evoked.data[sel]
    M = np.dot(whitener, M)

    alpha_max = norm_l2inf(np.dot(gain.T, M), n_orient=1)
    n_alphas = 5
    alpha_divs = np.logspace(np.log10(2), 1, n_alphas)
    times = np.zeros(n_alphas)
    tol = 1e-6
    for i, alpha_div in enumerate(alpha_divs):
        alpha = alpha_max / alpha_div
        t0 = time.time()
        X, active_set, E = mixed_norm_solver(M, gain, alpha,
                                             solver='cd',
                                             debias=False, tol=tol)
        times[i] = time.time() - t0

        R = M - np.dot(gain[:, active_set], X)
        p_obj = (R ** 2).sum() / 2. + alpha * norm_l21(X, n_orient=1)
        dual_scale = max(alpha,
                         norm_l2inf(np.dot(gain.T, R), n_orient=1))
        theta = R / dual_scale
        d_obj = (M ** 2).sum() / 2. - alpha ** 2 / 2. * norm(M / alpha - theta, ord='fro') ** 2

        assert p_obj - d_obj < tol

    times_dict[data] = times.copy()
    width = 0.35
    ind = np.arange(n_alphas)
    if data == "sample":
        ax.bar(ind, times, width, label="sample")
    else:
        ax.bar(ind + width, times, width, label="somato")


ax.set_ylabel('Time (s)')
branch = Repository('.').head.shorthand
ax.set_title('Time to reach a dual gap of %.0e on branch %s' % (tol, branch))
ax.set_xticks(ind + width / 2)
ax.set_xticklabels(["%.1f" % ad for ad in alpha_divs])
ax.set_xlabel("$\lambda_{max} / \lambda$")
plt.legend(loc='best')
plt.tight_layout()
plt.show()