{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# DICS for power mapping\n\nIn this tutorial, we'll simulate two signals originating from two\nlocations on the cortex. These signals will be sinusoids, so we'll be looking\nat oscillatory activity (as opposed to evoked activity).\n\nWe'll use dynamic imaging of coherent sources (DICS) :footcite:`GrossEtAl2001`\nto map out spectral power along the cortex. Let's see if we can find our two\nsimulated sources.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Author: Marijn van Vliet <w.m.vanvliet@gmail.com>\n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Setup\nWe first import the required packages to run this tutorial and define a list\nof filenames for various things we'll be using.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom matplotlib import pyplot as plt\nfrom scipy.signal import coherence, unit_impulse, welch\n\nimport mne\nfrom mne.beamformer import apply_dics_csd, make_dics\nfrom mne.datasets import sample\nfrom mne.minimum_norm import apply_inverse, make_inverse_operator\nfrom mne.simulation import add_noise, simulate_raw\nfrom mne.time_frequency import csd_morlet\n\n# We use the MEG and MRI setup from the MNE-sample dataset\ndata_path = sample.data_path(download=False)\nsubjects_dir = data_path / \"subjects\"\n\n# Filenames for various files we'll be using\nmeg_path = data_path / \"MEG\" / \"sample\"\nraw_fname = meg_path / \"sample_audvis_raw.fif\"\nfwd_fname = meg_path / \"sample_audvis-meg-eeg-oct-6-fwd.fif\"\ncov_fname = meg_path / \"sample_audvis-cov.fif\"\nfwd = mne.read_forward_solution(fwd_fname)\n\n# Seed for the random number generator\nrand = np.random.RandomState(42)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Data simulation\n\nThe following function generates a timeseries that contains an oscillator,\nwhose frequency fluctuates a little over time, but stays close to 10 Hz.\nWe'll use this function to generate our two signals.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sfreq = 50.0  # Sampling frequency of the generated signal\nn_samp = int(round(10.0 * sfreq))\ntimes = np.arange(n_samp) / sfreq  # 10 seconds of signal\nn_times = len(times)\n\n\ndef coh_signal_gen():\n    \"\"\"Generate an oscillating signal.\n\n    Returns\n    -------\n    signal : ndarray\n        The generated signal.\n    \"\"\"\n    t_rand = 0.001  # Variation in the instantaneous frequency of the signal\n    std = 0.1  # Std-dev of the random fluctuations added to the signal\n    base_freq = 10.0  # Base frequency of the oscillators in Hertz\n    n_times = len(times)\n\n    # Generate an oscillator with varying frequency and phase lag.\n    signal = np.sin(\n        2.0\n        * np.pi\n        * (\n            base_freq * np.arange(n_times) / sfreq\n            + np.cumsum(t_rand * rand.randn(n_times))\n        )\n    )\n\n    # Add some random fluctuations to the signal.\n    signal += std * rand.randn(n_times)\n\n    # Scale the signal to be in the right order of magnitude (~100 nAm)\n    # for MEG data.\n    signal *= 100e-9\n\n    return signal"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's simulate two timeseries and plot some basic information about them.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "signal1 = coh_signal_gen()\nsignal2 = coh_signal_gen()\n\nfig, axes = plt.subplots(2, 2, figsize=(8, 4), layout=\"constrained\")\n\n# Plot the timeseries\nax = axes[0][0]\nax.plot(times, 1e9 * signal1, lw=0.5)\nax.set(\n    xlabel=\"Time (s)\", xlim=times[[0, -1]], ylabel=\"Amplitude (Am)\", title=\"Signal 1\"\n)\nax = axes[0][1]\nax.plot(times, 1e9 * signal2, lw=0.5)\nax.set(xlabel=\"Time (s)\", xlim=times[[0, -1]], title=\"Signal 2\")\n\n# Power spectrum of the first timeseries\nf, p = welch(signal1, fs=sfreq, nperseg=128, nfft=256)\nax = axes[1][0]\n# Only plot the first 100 frequencies\nax.plot(f[:100], 20 * np.log10(p[:100]), lw=1.0)\nax.set(\n    xlabel=\"Frequency (Hz)\",\n    xlim=f[[0, 99]],\n    ylabel=\"Power (dB)\",\n    title=\"Power spectrum of signal 1\",\n)\n\n# Compute the coherence between the two timeseries\nf, coh = coherence(signal1, signal2, fs=sfreq, nperseg=100, noverlap=64)\nax = axes[1][1]\nax.plot(f[:50], coh[:50], lw=1.0)\nax.set(\n    xlabel=\"Frequency (Hz)\",\n    xlim=f[[0, 49]],\n    ylabel=\"Coherence\",\n    title=\"Coherence between the timeseries\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we put the signals at two locations on the cortex. We construct a\n:class:`mne.SourceEstimate` object to store them in.\n\nThe timeseries will have a part where the signal is active and a part where\nit is not. The techniques we'll be using in this tutorial depend on being\nable to contrast data that contains the signal of interest versus data that\ndoes not (i.e. it contains only noise).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# The locations on the cortex where the signal will originate from. These\n# locations are indicated as vertex numbers.\nvertices = [[146374], [33830]]\n\n# Construct SourceEstimates that describe the signals at the cortical level.\ndata = np.vstack((signal1, signal2))\nstc_signal = mne.SourceEstimate(\n    data, vertices, tmin=0, tstep=1.0 / sfreq, subject=\"sample\"\n)\nstc_noise = stc_signal * 0.0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Before we simulate the sensor-level data, let's define a signal-to-noise\nratio. You are encouraged to play with this parameter and see the effect of\nnoise on our results.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "snr = 1.0  # Signal-to-noise ratio. Decrease to add more noise."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we run the signal through the forward model to obtain simulated sensor\ndata. To save computation time, we'll only simulate gradiometer data. You can\ntry simulating other types of sensors as well.\n\nSome noise is added based on the baseline noise covariance matrix from the\nsample dataset, scaled to implement the desired SNR.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Read the info from the sample dataset. This defines the location of the\n# sensors and such.\ninfo = mne.io.read_raw(raw_fname).crop(0, 1).resample(50).info\n\n# Only use gradiometers\npicks = mne.pick_types(info, meg=\"grad\", stim=True, exclude=())\nmne.pick_info(info, picks, copy=False)  # modifies info in-place\n\n# Define a covariance matrix for the simulated noise. In this tutorial, we use\n# a simple diagonal matrix.\ncov = mne.cov.make_ad_hoc_cov(info)\ncov[\"data\"] *= (20.0 / snr) ** 2  # Scale the noise to achieve the desired SNR\n\n# Simulate the raw data, with a lowpass filter on the noise\nstcs = [\n    (stc_signal, unit_impulse(n_samp, dtype=int) * 1),\n    (stc_noise, unit_impulse(n_samp, dtype=int) * 2),\n]  # stacked in time\nduration = (len(stc_signal.times) * 2) / sfreq\nraw = simulate_raw(info, stcs, forward=fwd)\nadd_noise(raw, cov, iir_filter=[4, -4, 0.8], random_state=rand)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create an :class:`mne.Epochs` object containing two trials: one with\nboth noise and signal and one with just noise\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "events = mne.find_events(raw, initial_event=True)\ntmax = (len(stc_signal.times) - 1) / sfreq\nepochs = mne.Epochs(\n    raw,\n    events,\n    event_id=dict(signal=1, noise=2),\n    tmin=0,\n    tmax=tmax,\n    baseline=None,\n    preload=True,\n)\nassert len(epochs) == 2  # ensure that we got the two expected events\n\n# Plot some of the channels of the simulated data that are situated above one\n# of our simulated sources.\npicks = mne.read_vectorview_selection(\"Left-frontal\")  # contains both mag and grad\npicks = [p for p in picks if p in epochs.ch_names]  # now only grads\nepochs.plot(picks=picks, events=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Power mapping\nWith our simulated dataset ready, we can now pretend to be researchers that\nhave just recorded this from a real subject and are going to study what parts\nof the brain communicate with each other.\n\nFirst, we'll create a source estimate of the MEG data. We'll use both a\nstraightforward MNE-dSPM inverse solution for this, and the DICS beamformer\nwhich is specifically designed to work with oscillatory data.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Computing the inverse using MNE-dSPM:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Compute the inverse operator\nfwd = mne.read_forward_solution(fwd_fname)\ninv = make_inverse_operator(epochs.info, fwd, cov)\n\n# Apply the inverse model to the trial that also contains the signal.\ns = apply_inverse(epochs[\"signal\"].average(), inv)\n\n# Take the root-mean square along the time dimension and plot the result.\ns_rms = np.sqrt((s**2).mean())\ntitle = \"MNE-dSPM inverse (RMS)\"\nbrain = s_rms.plot(\n    \"sample\",\n    subjects_dir=subjects_dir,\n    hemi=\"both\",\n    figure=1,\n    size=600,\n    time_label=title,\n    title=title,\n)\n\n# Indicate the true locations of the source activity on the plot.\nbrain.add_foci(vertices[0][0], coords_as_verts=True, hemi=\"lh\")\nbrain.add_foci(vertices[1][0], coords_as_verts=True, hemi=\"rh\")\n\n# Rotate the view and add a title.\nbrain.show_view(azimuth=0, elevation=0, distance=550, focalpoint=(0, 0, 0))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will now compute the cortical power map at 10 Hz. using a DICS beamformer.\nA beamformer will construct for each vertex a spatial filter that aims to\npass activity originating from the vertex, while dampening activity from\nother sources as much as possible.\n\nThe :func:`mne.beamformer.make_dics` function has many switches that offer\nprecise control\nover the way the filter weights are computed. Currently, there is no clear\nconsensus regarding the best approach. This is why we will demonstrate two\napproaches here:\n\n 1. The approach as described in :footcite:`vanVlietEtAl2018`, which first\n    normalizes the forward solution and computes a vector beamformer.\n 2. The scalar beamforming approach based on\n    :footcite:`SekiharaNagarajan2008`, which uses weight normalization\n    instead of normalizing the forward solution.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Estimate the cross-spectral density (CSD) matrix on the trial containing the\n# signal.\ncsd_signal = csd_morlet(epochs[\"signal\"], frequencies=[10])\n\n# Compute the spatial filters for each vertex, using two approaches.\nfilters_approach1 = make_dics(\n    info,\n    fwd,\n    csd_signal,\n    reg=0.05,\n    pick_ori=\"max-power\",\n    depth=1.0,\n    inversion=\"single\",\n    weight_norm=None,\n    real_filter=True,\n)\nprint(filters_approach1)\n\nfilters_approach2 = make_dics(\n    info,\n    fwd,\n    csd_signal,\n    reg=0.05,\n    pick_ori=\"max-power\",\n    depth=None,\n    inversion=\"matrix\",\n    weight_norm=\"unit-noise-gain\",\n    real_filter=True,\n)\nprint(filters_approach2)\n\n# You can save these to disk with:\n# filters_approach1.save('filters_1-dics.h5')\n\n# Compute the DICS power map by applying the spatial filters to the CSD matrix.\npower_approach1, f = apply_dics_csd(csd_signal, filters_approach1)\npower_approach2, f = apply_dics_csd(csd_signal, filters_approach2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the DICS power maps for both approaches, starting with the first:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def plot_approach(power, n):\n    \"\"\"Plot the results on a brain.\"\"\"\n    title = f\"DICS power map, approach {n}\"\n    brain = power_approach1.plot(\n        \"sample\",\n        subjects_dir=subjects_dir,\n        hemi=\"both\",\n        size=600,\n        time_label=title,\n        title=title,\n    )\n    # Indicate the true locations of the source activity on the plot.\n    brain.add_foci(vertices[0][0], coords_as_verts=True, hemi=\"lh\", color=\"b\")\n    brain.add_foci(vertices[1][0], coords_as_verts=True, hemi=\"rh\", color=\"b\")\n    # Rotate the view and add a title.\n    brain.show_view(azimuth=0, elevation=0, distance=550, focalpoint=(0, 0, 0))\n    return brain\n\n\nbrain1 = plot_approach(power_approach1, 1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now the second:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "brain2 = plot_approach(power_approach2, 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Excellent! All methods found our two simulated sources. Of course, with a\nsignal-to-noise ratio (SNR) of 1, is isn't very hard to find them. You can\ntry playing with the SNR and see how the MNE-dSPM and DICS approaches hold up\nin the presence of increasing noise. In the presence of more noise, you may\nneed to increase the regularization parameter of the DICS beamformer.\n\n## 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
}