{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Compute MNE-dSPM inverse solution on single epochs\n\nCompute dSPM inverse solution on single trial epochs restricted\nto a brain label.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Author: Alexandre Gramfort <alexandre.gramfort@inria.fr>\n#\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, apply_inverse_epochs, read_inverse_operator\n\nprint(__doc__)\n\ndata_path = sample.data_path()\nmeg_path = data_path / \"MEG\" / \"sample\"\nfname_inv = meg_path / \"sample_audvis-meg-oct-6-meg-inv.fif\"\nfname_raw = meg_path / \"sample_audvis_filt-0-40_raw.fif\"\nfname_event = meg_path / \"sample_audvis_filt-0-40_raw-eve.fif\"\nlabel_name = \"Aud-lh\"\nfname_label = meg_path / \"labels\" / f\"{label_name}.label\"\n\nevent_id, tmin, tmax = 1, -0.2, 0.5\n\n# Using the same inverse operator when inspecting single trials Vs. evoked\nsnr = 3.0  # Standard assumption for average data but using it for single trial\nlambda2 = 1.0 / snr**2\n\nmethod = \"dSPM\"  # use dSPM method (could also be MNE or sLORETA)\n\n# Load data\ninverse_operator = read_inverse_operator(fname_inv)\nlabel = mne.read_label(fname_label)\nraw = mne.io.read_raw_fif(fname_raw)\nevents = mne.read_events(fname_event)\n\n# Set up pick list\ninclude = []\n\n# Add a bad channel\nraw.info[\"bads\"] += [\"EEG 053\"]  # bads + 1 more\n\n# pick MEG channels\npicks = mne.pick_types(\n    raw.info, meg=True, eeg=False, stim=False, eog=True, include=include, exclude=\"bads\"\n)\n# Read epochs\nepochs = mne.Epochs(\n    raw,\n    events,\n    event_id,\n    tmin,\n    tmax,\n    picks=picks,\n    baseline=(None, 0),\n    reject=dict(mag=4e-12, grad=4000e-13, eog=150e-6),\n)\n\n# Get evoked data (averaging across trials in sensor space)\nevoked = epochs.average()\n\n# Compute inverse solution and stcs for each epoch\n# Use the same inverse operator as with evoked data (i.e., set nave)\n# If you use a different nave, dSPM just scales by a factor sqrt(nave)\nstcs = apply_inverse_epochs(\n    epochs,\n    inverse_operator,\n    lambda2,\n    method,\n    label,\n    pick_ori=\"normal\",\n    nave=evoked.nave,\n)\n\n# Mean across trials but not across vertices in label\nmean_stc = sum(stcs) / len(stcs)\n\n# compute sign flip to avoid signal cancellation when averaging signed values\nflip = mne.label_sign_flip(label, inverse_operator[\"src\"])\n\nlabel_mean = np.mean(mean_stc.data, axis=0)\nlabel_mean_flip = np.mean(flip[:, np.newaxis] * mean_stc.data, axis=0)\n\n# Get inverse solution by inverting evoked data\nstc_evoked = apply_inverse(evoked, inverse_operator, lambda2, method, pick_ori=\"normal\")\n\n# apply_inverse() does whole brain, so sub-select label of interest\nstc_evoked_label = stc_evoked.in_label(label)\n\n# Average over label (not caring to align polarities here)\nlabel_mean_evoked = np.mean(stc_evoked_label.data, axis=0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "View activation time-series to illustrate the benefit of aligning/flipping\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "times = 1e3 * stcs[0].times  # times in ms\n\nplt.figure()\nh0 = plt.plot(times, mean_stc.data.T, \"k\")\n(h1,) = plt.plot(times, label_mean, \"r\", linewidth=3)\n(h2,) = plt.plot(times, label_mean_flip, \"g\", linewidth=3)\nplt.legend((h0[0], h1, h2), (\"all dipoles in label\", \"mean\", \"mean with sign flip\"))\nplt.xlabel(\"time (ms)\")\nplt.ylabel(\"dSPM value\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Viewing single trial dSPM and average dSPM for unflipped pooling over label\nCompare to (1) Inverse (dSPM) then average, (2) Evoked then dSPM\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Single trial\nplt.figure()\nfor k, stc_trial in enumerate(stcs):\n    plt.plot(\n        times,\n        np.mean(stc_trial.data, axis=0).T,\n        \"k--\",\n        label=\"Single Trials\" if k == 0 else \"_nolegend_\",\n        alpha=0.5,\n    )\n\n# Single trial inverse then average.. making linewidth large to not be masked\nplt.plot(times, label_mean, \"b\", linewidth=6, label=\"dSPM first, then average\")\n\n# Evoked and then inverse\nplt.plot(times, label_mean_evoked, \"r\", linewidth=2, label=\"Average first, then dSPM\")\n\nplt.xlabel(\"time (ms)\")\nplt.ylabel(\"dSPM value\")\nplt.legend()\nplt.show()"
      ]
    }
  ],
  "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
}