{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Compute evoked ERS source power using DICS, LCMV beamformer, and dSPM\n\nHere we examine 3 ways of localizing event-related synchronization (ERS) of\nbeta band activity in this dataset: `somato-dataset` using\n:term:`DICS`, :term:`LCMV beamformer`, and :term:`dSPM` applied to active and\nbaseline covariance matrices.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: Luke Bloy <luke.bloy@gmail.com>\n#          Eric Larson <larson.eric.d@gmail.com>\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.beamformer import apply_dics_csd, apply_lcmv_cov, make_dics, make_lcmv\nfrom mne.cov import compute_covariance\nfrom mne.datasets import somato\nfrom mne.minimum_norm import apply_inverse_cov, make_inverse_operator\nfrom mne.time_frequency import csd_morlet\n\nprint(__doc__)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Reading the raw data and creating epochs:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_path = somato.data_path()\nsubject = \"01\"\ntask = \"somato\"\nraw_fname = data_path / f\"sub-{subject}\" / \"meg\" / f\"sub-{subject}_task-{task}_meg.fif\"\n\n# crop to 5 minutes to save memory\nraw = mne.io.read_raw_fif(raw_fname).crop(0, 300)\n\n# We are interested in the beta band (12-30 Hz)\nraw.load_data().filter(12, 30)\n\n# The DICS beamformer currently only supports a single sensor type.\n# We'll use the gradiometers in this example.\npicks = mne.pick_types(raw.info, meg=\"grad\", exclude=\"bads\")\n\n# Read epochs\nevents = mne.find_events(raw)\nepochs = mne.Epochs(\n    raw, events, event_id=1, tmin=-1.5, tmax=2, picks=picks, preload=True, decim=3\n)\n\n# Read forward operator and point to freesurfer subject directory\nfname_fwd = (\n    data_path / \"derivatives\" / f\"sub-{subject}\" / f\"sub-{subject}_task-{task}-fwd.fif\"\n)\nsubjects_dir = data_path / \"derivatives\" / \"freesurfer\" / \"subjects\"\n\nfwd = mne.read_forward_solution(fname_fwd)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute covariances\nERS activity starts at 0.5 seconds after stimulus onset. Because these\ndata have been processed by MaxFilter directly (rather than MNE-Python's\nversion), we have to be careful to compute the rank with a more conservative\nthreshold in order to get the correct data rank (64). Once this is used in\ncombination with an advanced covariance estimator like \"shrunk\", the rank\nwill be correctly preserved.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rank = mne.compute_rank(epochs, tol=1e-6, tol_kind=\"relative\")\nactive_win = (0.5, 1.5)\nbaseline_win = (-1, 0)\nbaseline_cov = compute_covariance(\n    epochs,\n    tmin=baseline_win[0],\n    tmax=baseline_win[1],\n    method=\"shrunk\",\n    rank=rank,\n    verbose=True,\n)\nactive_cov = compute_covariance(\n    epochs,\n    tmin=active_win[0],\n    tmax=active_win[1],\n    method=\"shrunk\",\n    rank=rank,\n    verbose=True,\n)\n\n# Weighted averaging is already in the addition of covariance objects.\ncommon_cov = baseline_cov + active_cov\nbaseline_cov.plot(epochs.info)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute some source estimates\nHere we will use DICS, LCMV beamformer, and dSPM.\n\nSee `ex-inverse-source-power` for more information about DICS.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def _gen_dics(active_win, baseline_win, epochs):\n    freqs = np.logspace(np.log10(12), np.log10(30), 9)\n    csd = csd_morlet(epochs, freqs, tmin=-1, tmax=1.5, decim=20)\n    csd_baseline = csd_morlet(\n        epochs, freqs, tmin=baseline_win[0], tmax=baseline_win[1], decim=20\n    )\n    csd_ers = csd_morlet(\n        epochs, freqs, tmin=active_win[0], tmax=active_win[1], decim=20\n    )\n    filters = make_dics(\n        epochs.info,\n        fwd,\n        csd.mean(),\n        pick_ori=\"max-power\",\n        reduce_rank=True,\n        real_filter=True,\n        rank=rank,\n    )\n    stc_base, freqs = apply_dics_csd(csd_baseline.mean(), filters)\n    stc_act, freqs = apply_dics_csd(csd_ers.mean(), filters)\n    stc_act /= stc_base\n    return stc_act\n\n\n# generate lcmv source estimate\ndef _gen_lcmv(active_cov, baseline_cov, common_cov):\n    filters = make_lcmv(\n        epochs.info, fwd, common_cov, reg=0.05, noise_cov=None, pick_ori=\"max-power\"\n    )\n    stc_base = apply_lcmv_cov(baseline_cov, filters)\n    stc_act = apply_lcmv_cov(active_cov, filters)\n    stc_act /= stc_base\n    return stc_act\n\n\n# generate mne/dSPM source estimate\ndef _gen_mne(active_cov, baseline_cov, common_cov, fwd, info, method=\"dSPM\"):\n    inverse_operator = make_inverse_operator(info, fwd, common_cov)\n    stc_act = apply_inverse_cov(\n        active_cov, info, inverse_operator, method=method, verbose=True\n    )\n    stc_base = apply_inverse_cov(\n        baseline_cov, info, inverse_operator, method=method, verbose=True\n    )\n    stc_act /= stc_base\n    return stc_act\n\n\n# Compute source estimates\nstc_dics = _gen_dics(active_win, baseline_win, epochs)\nstc_lcmv = _gen_lcmv(active_cov, baseline_cov, common_cov)\nstc_dspm = _gen_mne(active_cov, baseline_cov, common_cov, fwd, epochs.info)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot source estimates\nDICS:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "brain_dics = stc_dics.plot(\n    hemi=\"rh\",\n    subjects_dir=subjects_dir,\n    subject=subject,\n    time_label=\"DICS source power in the 12-30 Hz frequency band\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "LCMV:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "brain_lcmv = stc_lcmv.plot(\n    hemi=\"rh\",\n    subjects_dir=subjects_dir,\n    subject=subject,\n    time_label=\"LCMV source power in the 12-30 Hz frequency band\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "dSPM:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "brain_dspm = stc_dspm.plot(\n    hemi=\"rh\",\n    subjects_dir=subjects_dir,\n    subject=subject,\n    time_label=\"dSPM source power in the 12-30 Hz frequency band\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For more advanced usage, see\n`mne-gui-addons:sphx_glr_auto_examples_evoked_ers_source_power.py`.\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
}