{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Compute MxNE with time-frequency sparse prior\n\nThe TF-MxNE solver is a distributed inverse method (like dSPM or sLORETA)\nthat promotes focal (sparse) sources (such as dipole fitting techniques)\n:footcite:`GramfortEtAl2013b,GramfortEtAl2011`.\nThe benefit of this approach is that:\n\n  - it is spatio-temporal without assuming stationarity (sources properties\n    can vary over time)\n  - activations are localized in space, time and frequency in one step.\n  - with a built-in filtering process based on a short time Fourier\n    transform (STFT), data does not need to be low passed (just high pass\n    to make the signals zero mean).\n  - the solver solves a convex optimization problem, hence cannot be\n    trapped in local minima.\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, tf_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()\nsubjects_dir = data_path / \"subjects\"\nmeg_path = data_path / \"MEG\" / \"sample\"\nfwd_fname = meg_path / \"sample_audvis-meg-eeg-oct-6-fwd.fif\"\nave_fname = meg_path / \"sample_audvis-no-filter-ave.fif\"\ncov_fname = meg_path / \"sample_audvis-shrunk-cov.fif\"\n\n# Read noise covariance matrix\ncov = mne.read_cov(cov_fname)\n\n# Handling average file\ncondition = \"Left visual\"\nevoked = mne.read_evokeds(ave_fname, condition=condition, baseline=(None, 0))\n# We make the window slightly larger than what you'll eventually be interested\n# in ([-0.05, 0.3]) to avoid edge effects.\nevoked.crop(tmin=-0.1, tmax=0.4)\n\n# Handling forward solution\nforward = mne.read_forward_solution(fwd_fname)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run solver\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# alpha parameter is between 0 and 100 (100 gives 0 active source)\nalpha = 40.0  # general regularization parameter\n# l1_ratio parameter between 0 and 1 promotes temporal smoothness\n# (0 means no temporal regularization)\nl1_ratio = 0.03  # temporal regularization parameter\n\nloose, depth = 0.2, 0.9  # loose orientation & depth weighting\n\n# Compute dSPM solution to be used as weights in MxNE\ninverse_operator = make_inverse_operator(\n    evoked.info, forward, cov, loose=loose, depth=depth\n)\nstc_dspm = apply_inverse(evoked, inverse_operator, lambda2=1.0 / 9.0, method=\"dSPM\")\n\n# Compute TF-MxNE inverse solution with dipole output\ndipoles, residual = tf_mixed_norm(\n    evoked,\n    forward,\n    cov,\n    alpha=alpha,\n    l1_ratio=l1_ratio,\n    loose=loose,\n    depth=depth,\n    maxit=200,\n    tol=1e-6,\n    weights=stc_dspm,\n    weights_min=8.0,\n    debias=True,\n    wsize=16,\n    tstep=4,\n    window=0.05,\n    return_as_dipoles=True,\n    return_residual=True,\n)\n\n# Crop to remove edges\nfor dip in dipoles:\n    dip.crop(tmin=-0.05, tmax=0.3)\nevoked.crop(tmin=-0.05, tmax=0.3)\nresidual.crop(tmin=-0.05, tmax=0.3)"
      ]
    },
    {
      "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)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot location of the strongest dipole with MRI slices\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "idx = 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:\n# for dip in dipoles:\n#     plot_dipole_locations(dip, forward['mri_head_t'], 'sample',\n#                           subjects_dir=subjects_dir, mode='orthoview',\n#                           idx='amplitude')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Show the evoked response and the residual for gradiometers\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ylim = dict(grad=[-120, 120])\nevoked.pick(picks=\"grad\", exclude=\"bads\")\nevoked.plot(\n    titles=dict(grad=\"Evoked Response: Gradiometers\"),\n    ylim=ylim,\n    proj=True,\n    time_unit=\"s\",\n)\n\nresidual.pick(picks=\"grad\", exclude=\"bads\")\nresidual.plot(\n    titles=dict(grad=\"Residuals: Gradiometers\"), ylim=ylim, proj=True, time_unit=\"s\"\n)"
      ]
    },
    {
      "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": [
        "plot_sparse_source_estimates(\n    forward[\"src\"],\n    stc,\n    bgcolor=(1, 1, 1),\n    opacity=0.1,\n    fig_name=f\"TF-MxNE (cond {condition})\",\n    modes=[\"sphere\"],\n    scale_factors=[1.0],\n)\n\ntime_label = \"TF-MxNE time=%0.2f ms\"\nclim = dict(kind=\"value\", lims=[10e-9, 15e-9, 20e-9])\nbrain = stc.plot(\n    \"sample\",\n    \"inflated\",\n    \"rh\",\n    views=\"medial\",\n    clim=clim,\n    time_label=time_label,\n    smoothing_steps=5,\n    subjects_dir=subjects_dir,\n    initial_time=150,\n    time_unit=\"ms\",\n)\nbrain.add_label(\"V1\", color=\"yellow\", scalar_thresh=0.5, borders=True)\nbrain.add_label(\"V2\", color=\"red\", scalar_thresh=0.5, borders=True)"
      ]
    },
    {
      "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
}