{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Source localization with MNE, dSPM, sLORETA, and eLORETA\n\nThe aim of this tutorial is to teach you how to compute and apply a linear\nminimum-norm inverse method on evoked/raw/epochs data.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: The MNE-Python contributors.\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.minimum_norm import apply_inverse, make_inverse_operator"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Process MEG data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_path = sample.data_path()\nraw_fname = data_path / \"MEG\" / \"sample\" / \"sample_audvis_filt-0-40_raw.fif\"\n\nraw = mne.io.read_raw_fif(raw_fname)  # already has an average reference\nevents = mne.find_events(raw, stim_channel=\"STI 014\")\n\nevent_id = dict(aud_l=1)  # event trigger and conditions\ntmin = -0.2  # start of each epoch (200ms before the trigger)\ntmax = 0.5  # end of each epoch (500ms after the trigger)\nraw.info[\"bads\"] = [\"MEG 2443\", \"EEG 053\"]  # mark known bad channels\nbaseline = (None, 0)  # means from the first instant to t = 0\nreject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)\n\nepochs = mne.Epochs(\n    raw,\n    events,\n    event_id,\n    tmin,\n    tmax,\n    proj=True,\n    picks=(\"meg\", \"eog\"),\n    baseline=baseline,\n    reject=reject,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute regularized noise covariance\nFor more details see `tut-compute-covariance`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "noise_cov = mne.compute_covariance(\n    epochs, tmax=0.0, method=[\"shrunk\", \"empirical\"], rank=None, verbose=True\n)\n\nfig_cov, fig_spectra = mne.viz.plot_cov(noise_cov, raw.info)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute the evoked response\nLet's just use the MEG channels for simplicity.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "evoked = epochs.average().pick(\"meg\")\nevoked.plot(time_unit=\"s\")\nevoked.plot_topomap(times=np.linspace(0.05, 0.15, 5), ch_type=\"mag\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "It's also a good idea to look at whitened data:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "evoked.plot_white(noise_cov, time_unit=\"s\")\ndel epochs, raw  # to save memory"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Inverse modeling: MNE/dSPM on evoked and raw data\nHere we first read the forward solution. You will likely need to compute\none for your own data -- see `tut-forward` for information on how\nto do it.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fname_fwd = data_path / \"MEG\" / \"sample\" / \"sample_audvis-meg-oct-6-fwd.fif\"\nfwd = mne.read_forward_solution(fname_fwd)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we make an MEG inverse operator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "inverse_operator = make_inverse_operator(\n    evoked.info, fwd, noise_cov, loose=0.2, depth=0.8\n)\ndel fwd"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>You can write the inverse operator to disk with:\n\n```\nfrom mne.minimum_norm import write_inverse_operator\nwrite_inverse_operator(\n    \"sample_audvis-meg-oct-6-inv.fif\", inverse_operator\n)</p></div>\n```\n## Compute inverse solution\nWe can use this to compute the inverse solution and obtain source time\ncourses:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "method = \"dSPM\"  # could choose MNE, sLORETA, or eLORETA instead\nsnr = 3.0\nlambda2 = 1.0 / snr**2\nstc, residual = apply_inverse(\n    evoked,\n    inverse_operator,\n    lambda2,\n    method=method,\n    pick_ori=None,\n    return_residual=True,\n    verbose=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualization\nWe can look at different dipole activations:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nax.plot(1e3 * stc.times, stc.data[::100, :].T)\nax.set(xlabel=\"time (ms)\", ylabel=f\"{method} value\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Examine the original data and the residual after fitting:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(2, 1)\nevoked.plot(axes=axes)\nfor ax in axes:\n    for text in list(ax.texts):\n        text.remove()\n    for line in ax.lines:\n        line.set_color(\"#98df81\")\nresidual.plot(axes=axes)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here we use peak getter to move visualization to the time point of the peak\nand draw a marker at the maximum peak vertex.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "vertno_max, time_max = stc.get_peak(hemi=\"rh\")\n\nsubjects_dir = data_path / \"subjects\"\nsurfer_kwargs = dict(\n    hemi=\"rh\",\n    subjects_dir=subjects_dir,\n    clim=dict(kind=\"value\", lims=[8, 12, 15]),\n    views=\"lateral\",\n    initial_time=time_max,\n    time_unit=\"s\",\n    size=(800, 800),\n    smoothing_steps=10,\n)\nbrain = stc.plot(**surfer_kwargs)\nbrain.add_foci(\n    vertno_max,\n    coords_as_verts=True,\n    hemi=\"rh\",\n    color=\"blue\",\n    scale_factor=0.6,\n    alpha=0.5,\n)\nbrain.add_text(\n    0.1, 0.9, \"dSPM (plus location of maximal activation)\", \"title\", font_size=14\n)\n\n# The documentation website's movie is generated with:\n# brain.save_movie(..., tmin=0.05, tmax=0.15, interpolation='linear',\n#                  time_dilation=20, framerate=10, time_viewer=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "There are many other ways to visualize and work with source data, see\nfor example:\n\n- `tut-viz-stcs`\n- `ex-morph-surface`\n- `ex-morph-volume`\n- `ex-vector-mne-solution`\n- `tut-dipole-orientations`\n- `tut-mne-fixed-free`\n- `examples using apply_inverse\n  <sphx_glr_backreferences_mne.minimum_norm.apply_inverse>`.\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
}