{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Source localization with a custom inverse solver\n\nThe objective of this example is to show how to plug a custom inverse solver\nin MNE in order to facilate empirical comparison with the methods MNE already\nimplements (wMNE, dSPM, sLORETA, eLORETA, LCMV, DICS, (TF-)MxNE etc.).\n\nThis script is educational and shall be used for methods\nevaluations and new developments. It is not meant to be an example\nof good practice to analyse your data.\n\nThe example makes use of 2 functions ``apply_solver`` and ``solver``\nso changes can be limited to the ``solver`` function (which only takes three\nparameters: the whitened data, the gain matrix and the number of orientations)\nin order to try out another inverse algorithm.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Authors: The MNE-Python contributors.\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\nfrom scipy import linalg\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.viz import plot_sparse_source_estimates\n\ndata_path = sample.data_path()\nmeg_path = data_path / \"MEG\" / \"sample\"\nfwd_fname = meg_path / \"sample_audvis-meg-eeg-oct-6-fwd.fif\"\nave_fname = meg_path / \"sample_audvis-ave.fif\"\ncov_fname = meg_path / \"sample_audvis-shrunk-cov.fif\"\nsubjects_dir = data_path / \"subjects\"\ncondition = \"Left Auditory\"\n\n# Read noise covariance matrix\nnoise_cov = mne.read_cov(cov_fname)\n# Handling average file\nevoked = mne.read_evokeds(ave_fname, condition=condition, baseline=(None, 0))\nevoked.crop(tmin=0.04, tmax=0.18)\n\nevoked = evoked.pick(picks=\"meg\", exclude=\"bads\")\n# Handling forward solution\nforward = mne.read_forward_solution(fwd_fname)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Auxiliary function to run the solver\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def apply_solver(solver, evoked, forward, noise_cov, loose=0.2, depth=0.8):\n \"\"\"Call a custom solver on evoked data.\n\n This function does all the necessary computation:\n\n - to select the channels in the forward given the available ones in\n the data\n - to take into account the noise covariance and do the spatial whitening\n - to apply loose orientation constraint as MNE solvers\n - to apply a weigthing of the columns of the forward operator as in the\n weighted Minimum Norm formulation in order to limit the problem\n of depth bias.\n\n Parameters\n ----------\n solver : callable\n The solver takes 3 parameters: data M, gain matrix G, number of\n dipoles orientations per location (1 or 3). A solver shall return\n 2 variables: X which contains the time series of the active dipoles\n and an active set which is a boolean mask to specify what dipoles are\n present in X.\n evoked : instance of mne.Evoked\n The evoked data\n forward : instance of Forward\n The forward solution.\n noise_cov : instance of Covariance\n The noise covariance.\n loose : float in [0, 1] | 'auto'\n Value that weights the source variances of the dipole components\n that are parallel (tangential) to the cortical surface. If loose\n is 0 then the solution is computed with fixed orientation.\n If loose is 1, it corresponds to free orientations.\n The default value ('auto') is set to 0.2 for surface-oriented source\n space and set to 1.0 for volumic or discrete source space.\n depth : None | float in [0, 1]\n Depth weighting coefficients. If None, no depth weighting is performed.\n\n Returns\n -------\n stc : instance of SourceEstimate\n The source estimates.\n \"\"\"\n # Import the necessary private functions\n from mne.inverse_sparse.mxne_inverse import (\n _make_sparse_stc,\n _prepare_gain,\n _reapply_source_weighting,\n is_fixed_orient,\n )\n\n all_ch_names = evoked.ch_names\n\n # Handle depth weighting and whitening (here is no weights)\n forward, gain, gain_info, whitener, source_weighting, mask = _prepare_gain(\n forward,\n evoked.info,\n noise_cov,\n pca=False,\n depth=depth,\n loose=loose,\n weights=None,\n weights_min=None,\n rank=None,\n )\n\n # Select channels of interest\n sel = [all_ch_names.index(name) for name in gain_info[\"ch_names\"]]\n M = evoked.data[sel]\n\n # Whiten data\n M = np.dot(whitener, M)\n\n n_orient = 1 if is_fixed_orient(forward) else 3\n X, active_set = solver(M, gain, n_orient)\n X = _reapply_source_weighting(X, source_weighting, active_set)\n\n stc = _make_sparse_stc(\n X, active_set, forward, tmin=evoked.times[0], tstep=1.0 / evoked.info[\"sfreq\"]\n )\n\n return stc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define your solver\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def solver(M, G, n_orient):\n \"\"\"Run L2 penalized regression and keep 10 strongest locations.\n\n Parameters\n ----------\n M : array, shape (n_channels, n_times)\n The whitened data.\n G : array, shape (n_channels, n_dipoles)\n The gain matrix a.k.a. the forward operator. The number of locations\n is n_dipoles / n_orient. n_orient will be 1 for a fixed orientation\n constraint or 3 when using a free orientation model.\n n_orient : int\n Can be 1 or 3 depending if one works with fixed or free orientations.\n If n_orient is 3, then ``G[:, 2::3]`` corresponds to the dipoles that\n are normal to the cortex.\n\n Returns\n -------\n X : array, (n_active_dipoles, n_times)\n The time series of the dipoles in the active set.\n active_set : array (n_dipoles)\n Array of bool. Entry j is True if dipole j is in the active set.\n We have ``X_full[active_set] == X`` where X_full is the full X matrix\n such that ``M = G X_full``.\n \"\"\"\n inner = np.dot(G, G.T)\n trace = np.trace(inner)\n K = linalg.solve(inner + 4e-6 * trace * np.eye(G.shape[0]), G).T\n K /= np.linalg.norm(K, axis=1)[:, None]\n X = np.dot(K, M)\n\n indices = np.argsort(np.sum(X**2, axis=1))[-10:]\n active_set = np.zeros(G.shape[1], dtype=bool)\n for idx in indices:\n idx -= idx % n_orient\n active_set[idx : idx + n_orient] = True\n X = X[active_set]\n return X, active_set" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Apply your custom solver\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# loose, depth = 0.2, 0.8 # corresponds to loose orientation\nloose, depth = 1.0, 0.0 # corresponds to free orientation\nstc = apply_solver(solver, evoked, forward, noise_cov, loose, depth)" ] }, { "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(forward[\"src\"], stc, bgcolor=(1, 1, 1), opacity=0.1)" ] } ], "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 }