mixed_norm_visual_left.py

mixed_norm_visual_left.py

import numpy as np

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

data_path = sample.data_path()
meg_path = data_path / "MEG" / "sample"
fwd_fname = meg_path / "sample_audvis-meg-eeg-oct-6-fwd.fif"
ave_fname = meg_path / "sample_audvis-ave.fif"
cov_fname = meg_path / "sample_audvis-shrunk-cov.fif"
subjects_dir = data_path / "subjects"
cov = mne.read_cov(cov_fname)
condition = "Left visual"
evoked = mne.read_evokeds(ave_fname, condition=condition, baseline=(None, 0))
evoked.crop(tmin=0, tmax=0.3)
evoked.plot
forward = mne.read_forward_solution(fwd_fname)
alpha = 30.  # 40 fits 2 dipoles, 30 gives 3 dipoles
loose, depth = 0.9, 0.9  # loose orientation & depth weighting
n_mxne_iter = 10  # if > 1 use L0.5/L2 reweighted mixed norm solver
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.0 / 9.0, method="dSPM")
dipoles, residual = mixed_norm(
    evoked,
    forward,
    cov,
    alpha=alpha,
    loose=loose,
    depth=depth,
    maxit=3000,
    tol=1e-4,
    active_set_size=10,
    debias=False,
    weights=stc_dspm,
    weights_min=8.0,
    n_mxne_iter=n_mxne_iter,
    return_residual=True,
    return_as_dipoles=True,
    verbose=True,
    random_state=0,
)

for di, dip in enumerate(dipoles, 1):
    tidx = dip.gof.argmax()
    t = dip.times[tidx]
    print(f"Dipole #{di} GOF at {1000 * t:0.1f} ms: {float(dip.gof[tidx]):0.1f}%")

plot_dipole_amplitudes(dipoles)
idx = np.argmax([np.max(np.abs(dip.amplitude)) for dip in dipoles])
plot_dipole_locations(
    dipoles[idx],
    forward["mri_head_t"],
    "sample",
    subjects_dir=subjects_dir,
    mode="orthoview",
    idx="amplitude",
)
for dip in dipoles:
    plot_dipole_locations(
        dip,
        forward["mri_head_t"],
        "sample",
        subjects_dir=subjects_dir,
        mode="orthoview",
        idx="amplitude",
    )

ylim = dict(eeg=[-10, 10], grad=[-400, 400], mag=[-600, 600])
evoked.pick(picks=["meg", "eeg"], exclude="bads")
fig = evoked.plot(ylim=ylim, proj=True, time_unit="s")
fig.suptitle("Evoked data", fontsize=20)
residual.pick(picks=["meg", "eeg"], exclude="bads")
fig = residual.plot(ylim=ylim, proj=True, time_unit="s")
fig.suptitle("Residual", fontsize=20)
stc = make_stc_from_dipoles(dipoles, forward["src"])
solver = "MxNE" if n_mxne_iter == 1 else "irMxNE"
plot_sparse_source_estimates(
    forward["src"],
    stc,
    bgcolor=(1, 1, 1),
    fig_name=f"{solver} (cond {condition})",
    opacity=0.1,
)