{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Compute sparse inverse solution with mixed norm: MxNE and irMxNE\n\nRuns an (ir)MxNE (L1/L2 :footcite:`GramfortEtAl2012` or L0.5/L2\n:footcite:`StrohmeierEtAl2014` mixed norm) inverse solver.\nL0.5/L2 is done with irMxNE which allows for sparser source estimates with less\namplitude bias due to the non-convexity of the L0.5/L2 mixed norm penalty.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Author: Alexandre Gramfort <alexandre.gramfort@inria.fr>\n#         Daniel Strohmeier <daniel.strohmeier@tu-ilmenau.de>\n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.inverse_sparse import make_stc_from_dipoles, mixed_norm\nfrom mne.minimum_norm import apply_inverse, make_inverse_operator\nfrom mne.viz import (\n    plot_dipole_amplitudes,\n    plot_dipole_locations,\n    plot_sparse_source_estimates,\n)\n\nprint(__doc__)\n\ndata_path = sample.data_path()\nmeg_path = data_path / \"MEG\" / \"sample\"\nfwd_fname = meg_path / \"sample_audvis-meg-eeg-oct-6-fwd.fif\"\nave_fname = meg_path / \"sample_audvis-ave.fif\"\ncov_fname = meg_path / \"sample_audvis-shrunk-cov.fif\"\nsubjects_dir = data_path / \"subjects\"\n\n# Read noise covariance matrix\ncov = mne.read_cov(cov_fname)\n# Handling average file\ncondition = \"Left Auditory\"\nevoked = mne.read_evokeds(ave_fname, condition=condition, baseline=(None, 0))\nevoked.crop(tmin=0, tmax=0.3)\n# Handling forward solution\nforward = mne.read_forward_solution(fwd_fname)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run solver with SURE criterion :footcite:`DeledalleEtAl2014`\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "alpha = \"sure\"  # regularization parameter between 0 and 100 or SURE criterion\nloose, depth = 0.9, 0.9  # loose orientation & depth weighting\nn_mxne_iter = 10  # if > 1 use L0.5/L2 reweighted mixed norm solver\n# if n_mxne_iter > 1 dSPM weighting can be avoided.\n\n# Compute dSPM solution to be used as weights in MxNE\ninverse_operator = make_inverse_operator(\n    evoked.info, forward, cov, depth=depth, fixed=True, use_cps=True\n)\nstc_dspm = apply_inverse(evoked, inverse_operator, lambda2=1.0 / 9.0, method=\"dSPM\")\n\n# Compute (ir)MxNE inverse solution with dipole output\ndipoles, residual = mixed_norm(\n    evoked,\n    forward,\n    cov,\n    alpha,\n    loose=loose,\n    depth=depth,\n    maxit=3000,\n    tol=1e-4,\n    active_set_size=10,\n    debias=False,\n    weights=stc_dspm,\n    weights_min=8.0,\n    n_mxne_iter=n_mxne_iter,\n    return_residual=True,\n    return_as_dipoles=True,\n    verbose=True,\n    random_state=0,\n    # for this dataset we know we should use a high alpha, so avoid some\n    # of the slower (lower) alpha values\n    sure_alpha_grid=np.linspace(100, 40, 10),\n)\n\nt = 0.083\ntidx = evoked.time_as_index(t).item()\nfor di, dip in enumerate(dipoles, 1):\n    print(f\"Dipole #{di} GOF at {1000 * t:0.1f} ms: {float(dip.gof[tidx]):0.1f}%\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot dipole activations\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_dipole_amplitudes(dipoles)\n\n# Plot dipole location of the strongest dipole with MRI slices\nidx = np.argmax([np.max(np.abs(dip.amplitude)) for dip in dipoles])\nplot_dipole_locations(\n    dipoles[idx],\n    forward[\"mri_head_t\"],\n    \"sample\",\n    subjects_dir=subjects_dir,\n    mode=\"orthoview\",\n    idx=\"amplitude\",\n)\n\n# Plot dipole locations of all dipoles with MRI slices\nfor dip in dipoles:\n    plot_dipole_locations(\n        dip,\n        forward[\"mri_head_t\"],\n        \"sample\",\n        subjects_dir=subjects_dir,\n        mode=\"orthoview\",\n        idx=\"amplitude\",\n    )"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot residual\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ylim = dict(eeg=[-10, 10], grad=[-400, 400], mag=[-600, 600])\nevoked.pick(picks=[\"meg\", \"eeg\"], exclude=\"bads\")\nevoked.plot(ylim=ylim, proj=True, time_unit=\"s\")\nresidual.pick(picks=[\"meg\", \"eeg\"], exclude=\"bads\")\nresidual.plot(ylim=ylim, proj=True, time_unit=\"s\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Generate stc from dipoles\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "stc = make_stc_from_dipoles(dipoles, forward[\"src\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "View in 2D and 3D (\"glass\" brain like 3D plot)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "solver = \"MxNE\" if n_mxne_iter == 1 else \"irMxNE\"\nplot_sparse_source_estimates(\n    forward[\"src\"],\n    stc,\n    bgcolor=(1, 1, 1),\n    fig_name=f\"{solver} (cond {condition})\",\n    opacity=0.1,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Morph onto fsaverage brain and view\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "morph = mne.compute_source_morph(\n    stc,\n    subject_from=\"sample\",\n    subject_to=\"fsaverage\",\n    spacing=None,\n    sparse=True,\n    subjects_dir=subjects_dir,\n)\nstc_fsaverage = morph.apply(stc)\nsrc_fsaverage_fname = subjects_dir / \"fsaverage\" / \"bem\" / \"fsaverage-ico-5-src.fif\"\nsrc_fsaverage = mne.read_source_spaces(src_fsaverage_fname)\n\nplot_sparse_source_estimates(\n    src_fsaverage,\n    stc_fsaverage,\n    bgcolor=(1, 1, 1),\n    fig_name=f\"Morphed {solver} (cond {condition})\",\n    opacity=0.1,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n.. footbibliography::\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}