fail_reweighted_right_audi.py

"""
================================================================
Compute sparse inverse solution with mixed norm: MxNE and irMxNE
================================================================

Runs an (ir)MxNE (L1/L2 [1]_ or L0.5/L2 [2]_ mixed norm) inverse solver.
L0.5/L2 is done with irMxNE which allows for sparser
source estimates with less amplitude bias due to the non-convexity
of the L0.5/L2 mixed norm penalty.

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
"""
#
# License: BSD (3-clause)

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'

# Read noise covariance matrix
cov = mne.read_cov(cov_fname)
# Handling average file
condition = 'Right Auditory'
# 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)

###############################################################################
# Run solver
if condition =='Left Auditory':
    alpha = 55  # regularization parameter between 0 and 100 (100 is high)
else:
    alpha = 29.5  # regularization parameter between 0 and 100 (100 is high)
loose, depth = 0.2, 0.9  # loose orientation & depth weighting
n_mxne_iter = 10  # if > 1 use L0.5/L2 reweighted mixed norm solver
# if n_mxne_iter > 1 dSPM weighting can be avoided.

# Compute dSPM solution to be used as weights in MxNE
inverse_operator = make_inverse_operator(evoked.info, forward, cov,
                                         depth=depth, fixed=True,
                                         use_cps=True)
stc_dspm = apply_inverse(evoked, inverse_operator, lambda2=1. / 9.,
                         method='dSPM')

# Compute (ir)MxNE inverse solution with dipole output
dipoles, residual = mixed_norm(
    evoked, forward, cov, alpha, loose=loose, depth=depth, maxit=3000,
    tol=1e-4, active_set_size=10, debias=True, weights=stc_dspm,
    weights_min=8., n_mxne_iter=n_mxne_iter, return_residual=True,
    return_as_dipoles=True)

###############################################################################
# Plot dipole location of the strongest dipole with MRI slices
idx = np.argmax([np.max(np.abs(dip.amplitude)) for dip in dipoles])

# Plot residual
ylim = dict(eeg=[-10, 10], grad=[-400, 400], mag=[-600, 600])

###############################################################################
# Generate stc from dipoles
stc = make_stc_from_dipoles(dipoles, forward['src'])

###############################################################################
# View in 2D and 3D ("glass" brain like 3D plot)
solver = "MxNE" if n_mxne_iter == 1 else "irMxNE"
plot_sparse_source_estimates(forward['src'], stc, bgcolor=(1, 1, 1),
                             fig_name="%s (cond %s)" % (solver, condition),
                             opacity=0.1)

@QB3

forwards_dip = np.zeros([forward["sol"]["data"].shape[0], 9])
l_fwd = []
acts_dip = np.zeros([181, 9])
for i, dip in enumerate(dipoles):
    lh_dist = np.linalg.norm(forward["src"][0]["rr"] - dip.pos[0], axis=1)
    rh_dist = np.linalg.norm(forward["src"][1]["rr"] - dip.pos[0], axis=1)
    lh = np.argmin(lh_dist)
    rh = np.argmin(rh_dist)

    if lh_dist[lh] == 0:
        hemi_idx = 0
        vert_no = lh
    elif rh_dist[rh] == 0:
        hemi_idx = 1
        vert_no = rh
    else:
        raise ValueError("did not exactly find position in LH or RH sources")

    vert_idx = np.sum(forward["src"][hemi_idx]["inuse"][:vert_no])
    if hemi_idx == 1:
        vert_idx += np.sum(forward["src"][0]["inuse"])
    print(vert_idx)

    g = 3 * vert_idx
    forwards_dip[:, 3 * i:3 * i + 3] = forward["sol"]["data"][:, g:g + 3]
    l_fwd.append(forward["sol"]["data"][:, g:g + 3])
    acts_dip[:, 3 * i:3 * i + 3] = dip.ori * dip.amplitude[:, None]
    # wanted to check residuals but 2 channels are discarded

print((l_fwd[0] * l_fwd[1]).sum() / norm(l_fwd[0]) / norm(l_fwd[1]))
print(np.corrcoef(forwards_dip.reshape(-1, 3, order='F').T))  # more elegant

and you need to import norm. Does the code to obtain the forward fields for the active dipoles seem correct?