bench_mixed_norm.py

"""
================================================================
Benchmark pure python vs blas Block Coordinate Descent
================================================================

References
----------
.. [1] Gramfort A., Kowalski M. and Hamalainen, M.
   "Mixed-norm estimates for the M/EEG inverse problem using accelerated
   gradient methods", Physics in Medicine and Biology, 2012.
   https://doi.org/10.1088/0031-9155/57/7/1937.

.. [2] Strohmeier D., Haueisen J., and Gramfort A.
   "Improved MEG/EEG source localization with reweighted mixed-norms",
   4th International Workshop on Pattern Recognition in Neuroimaging,
   Tuebingen, 2014. 10.1109/PRNI.2014.6858545
"""
# Author: Quentin Bertrand <quentin.bertrand@inria.fr>
#
# License: BSD (3-clause)

import time
import numpy as np

import mne
from mne.datasets import sample
from mne.inverse_sparse import mixed_norm, make_stc_from_dipoles
from mne.minimum_norm import make_inverse_operator, apply_inverse
from mne.viz import (plot_sparse_source_estimates,
                     plot_dipole_locations, plot_dipole_amplitudes)


print(__doc__)

data_path = sample.data_path()
fwd_fname = 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'
subjects_dir = data_path + '/subjects'

# parameters of the data
force_fixed = False
# force_fixed = True

# parameters of the solvers
tol = 1e-4
# solver = 'auto'
solver = 'bcd'
maxit = 3000
# Read noise covariance matrix
cov = mne.read_cov(cov_fname)
# Handling average file
condition = 'Left Auditory'
evoked = mne.read_evokeds(ave_fname, condition=condition, baseline=(None, 0))
evoked.crop(tmin=0, tmax=0.3)
# Handling forward solution
forward = mne.read_forward_solution(fwd_fname)
if force_fixed:
    forward = mne.convert_forward_solution(forward, force_fixed=forced_fixed)

###############################################################################
# Run solver
alpha = 15  # regularization parameter between 0 and 100 (100 is high)
# loose, depth = 0, None  # loose orientation & depth weighting
loose, depth = 0.2, 0.9  # loose orientation & depth weighting
n_mxne_iter = 10  # if > 1 use L0.5/L2 reweighted mixed norm solver

t0 = time.time()
dipoles, residual = mixed_norm(
    evoked, forward, cov, alpha, loose=loose, depth=depth, maxit=maxit,
    tol=tol, active_set_size=10, debias=False, weights=None,
    weights_min=8., n_mxne_iter=n_mxne_iter, return_residual=True,
    return_as_dipoles=True, solver=solver)
tF = time.time() - t0

print("Time ellapsed", tF)