{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Compare simulated and estimated source activity\n\nThis example illustrates how to compare the simulated and estimated\nsource time courses (STC) by computing different metrics. Simulated\nsource is a cortical region or dipole. It is meant to be a brief\nintroduction and only highlights the simplest use case.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Author: Kostiantyn Maksymenko \n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors.\n\nfrom functools import partial\n\nimport 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\nfrom mne.simulation.metrics import (\n cosine_score,\n f1_score,\n peak_position_error,\n precision_score,\n recall_score,\n region_localization_error,\n spatial_deviation_error,\n)\n\nrandom_state = 42 # set random state to make this example deterministic\n\n# Import sample data\ndata_path = sample.data_path()\nsubjects_dir = data_path / \"subjects\"\nsubject = \"sample\"\nevoked_fname = data_path / \"MEG\" / subject / \"sample_audvis-ave.fif\"\ninfo = mne.io.read_info(evoked_fname)\ntstep = 1.0 / info[\"sfreq\"]\n\n# Import forward operator and source space\nfwd_fname = data_path / \"MEG\" / subject / \"sample_audvis-meg-eeg-oct-6-fwd.fif\"\nfwd = mne.read_forward_solution(fwd_fname)\nsrc = fwd[\"src\"]\n\n# To select source, we use the caudal middle frontal to grow\n# a region of interest.\nselected_label = mne.read_labels_from_annot(\n subject, regexp=\"caudalmiddlefrontal-lh\", subjects_dir=subjects_dir\n)[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example we simulate two types of cortical sources: a region and\na dipole sources. We will test corresponding performance metrics.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define main parameters of sources\n\nFirst we define both region and dipole sources in terms of\nWhere?, What? and When?.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# WHERE?\n\n# Region\nlocation = \"center\" # Use the center of the label as a seed.\nextent = 20.0 # Extent in mm of the region.\nlabel_region = mne.label.select_sources(\n subject,\n selected_label,\n location=location,\n extent=extent,\n subjects_dir=subjects_dir,\n random_state=random_state,\n)\n\n# Dipole\nlocation = 1915 # Use the index of the vertex as a seed\nextent = 0.0 # One dipole source\nlabel_dipole = mne.label.select_sources(\n subject,\n selected_label,\n location=location,\n extent=extent,\n subjects_dir=subjects_dir,\n random_state=random_state,\n)\n\n# WHAT?\n# Define the time course of the activity\nsource_time_series = np.sin(2.0 * np.pi * 18.0 * np.arange(100) * tstep) * 10e-9\n\n# WHEN?\n# Define when the activity occurs using events.\nn_events = 50\nevents = np.zeros((n_events, 3), int)\nevents[:, 0] = 200 * np.arange(n_events) # Events sample.\nevents[:, 2] = 1 # All events have the sample id." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create simulated source activity\n\nHere, :class:`~mne.simulation.SourceSimulator` is used.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Region\nsource_simulator_region = mne.simulation.SourceSimulator(src, tstep=tstep)\nsource_simulator_region.add_data(label_region, source_time_series, events)\n\n# Dipole\nsource_simulator_dipole = mne.simulation.SourceSimulator(src, tstep=tstep)\nsource_simulator_dipole.add_data(label_dipole, source_time_series, events)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulate raw data\n\nProject the source time series to sensor space with multivariate Gaussian\nnoise obtained from the noise covariance from the sample data.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Region\nraw_region = mne.simulation.simulate_raw(info, source_simulator_region, forward=fwd)\nraw_region = raw_region.pick(picks=[\"eeg\", \"stim\"], exclude=\"bads\")\ncov = mne.make_ad_hoc_cov(raw_region.info)\nmne.simulation.add_noise(\n raw_region, cov, iir_filter=[0.2, -0.2, 0.04], random_state=random_state\n)\n\n# Dipole\nraw_dipole = mne.simulation.simulate_raw(info, source_simulator_dipole, forward=fwd)\nraw_dipole = raw_dipole.pick(picks=[\"eeg\", \"stim\"], exclude=\"bads\")\ncov = mne.make_ad_hoc_cov(raw_dipole.info)\nmne.simulation.add_noise(\n raw_dipole, cov, iir_filter=[0.2, -0.2, 0.04], random_state=random_state\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute evoked from raw data\n\nAveraging epochs corresponding to events.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Region\nevents = mne.find_events(raw_region, initial_event=True)\ntmax = (len(source_time_series) - 1) * tstep\nepochs = mne.Epochs(raw_region, events, 1, tmin=0, tmax=tmax, baseline=None)\nevoked_region = epochs.average()\n\n# Dipole\nevents = mne.find_events(raw_dipole, initial_event=True)\ntmax = (len(source_time_series) - 1) * tstep\nepochs = mne.Epochs(raw_dipole, events, 1, tmin=0, tmax=tmax, baseline=None)\nevoked_dipole = epochs.average()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create true stcs corresponding to evoked\n\nBefore we computed stcs corresponding to raw data. To be able to compare\nit with the reconstruction, based on the evoked, true stc should have the\nsame number of time samples.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Region\nstc_true_region = source_simulator_region.get_stc(\n start_sample=0, stop_sample=len(source_time_series)\n)\n\n# Dipole\nstc_true_dipole = source_simulator_dipole.get_stc(\n start_sample=0, stop_sample=len(source_time_series)\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reconstruct simulated sources\n\nCompute inverse solution using sLORETA.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Region\nsnr = 30.0\ninv_method = \"sLORETA\"\nlambda2 = 1.0 / snr**2\n\ninverse_operator = make_inverse_operator(\n evoked_region.info, fwd, cov, loose=\"auto\", depth=0.8, fixed=True\n)\n\nstc_est_region = apply_inverse(\n evoked_region, inverse_operator, lambda2, inv_method, pick_ori=None\n)\n\n# Dipole\nsnr = 3.0\ninv_method = \"sLORETA\"\nlambda2 = 1.0 / snr**2\n\ninverse_operator = make_inverse_operator(\n evoked_dipole.info, fwd, cov, loose=\"auto\", depth=0.8, fixed=True\n)\n\nstc_est_dipole = apply_inverse(\n evoked_dipole, inverse_operator, lambda2, inv_method, pick_ori=None\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute performance scores for different source amplitude thresholds\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "thresholds = [10, 30, 50, 70, 80, 90, 95, 99]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### For region\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# create a set of scorers\nscorers = {\n \"RLE\": partial(region_localization_error, src=src),\n \"Precision\": precision_score,\n \"Recall\": recall_score,\n \"F1 score\": f1_score,\n}\n\n# compute results\nregion_results = {}\nfor name, scorer in scorers.items():\n region_results[name] = [\n scorer(stc_true_region, stc_est_region, threshold=f\"{thx}%\", per_sample=False)\n for thx in thresholds\n ]\n\n# Plot the results\nf, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=\"col\", layout=\"constrained\")\nfor ax, (title, results) in zip([ax1, ax2, ax3, ax4], region_results.items()):\n ax.plot(thresholds, results, \".-\")\n ax.set(title=title, ylabel=\"score\", xlabel=\"Threshold\", xticks=thresholds)\n\nf.suptitle(\"Performance scores per threshold\") # Add Super title\nax1.ticklabel_format(axis=\"y\", style=\"sci\", scilimits=(0, 1)) # tweak RLE\n\n# Cosine score with respect to time\nf, ax1 = plt.subplots(layout=\"constrained\")\nax1.plot(stc_true_region.times, cosine_score(stc_true_region, stc_est_region))\nax1.set(title=\"Cosine score\", xlabel=\"Time\", ylabel=\"Score\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### For Dipoles\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# create a set of scorers\nscorers = {\n \"Peak Position Error\": peak_position_error,\n \"Spatial Deviation Error\": spatial_deviation_error,\n}\n\n\n# compute results\ndipole_results = {}\nfor name, scorer in scorers.items():\n dipole_results[name] = [\n scorer(\n stc_true_dipole,\n stc_est_dipole,\n src=src,\n threshold=f\"{thx}%\",\n per_sample=False,\n )\n for thx in thresholds\n ]\n\n\n# Plot the results\nfor name, results in dipole_results.items():\n f, ax1 = plt.subplots(layout=\"constrained\")\n ax1.plot(thresholds, 100 * np.array(results), \".-\")\n ax1.set(title=name, ylabel=\"Error (cm)\", xlabel=\"Threshold\", xticks=thresholds)" ] } ], "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 }