Note
Click here to download the full example code
Compute a sparse inverse solution using the Gamma-MAP empirical Bayesian method¶
See 1 for details.
# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>
# Daniel Strohmeier <daniel.strohmeier@tu-ilmenau.de>
#
# License: BSD (3-clause)
import numpy as np
import mne
from mne.datasets import sample
from mne.inverse_sparse import gamma_map, make_stc_from_dipoles
from mne.viz import (plot_sparse_source_estimates,
plot_dipole_locations, plot_dipole_amplitudes)
print(__doc__)
data_path = sample.data_path()
subjects_dir = data_path + '/subjects'
fwd_fname = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif'
evoked_fname = data_path + '/MEG/sample/sample_audvis-ave.fif'
cov_fname = data_path + '/MEG/sample/sample_audvis-cov.fif'
# Read the evoked response and crop it
condition = 'Left visual'
evoked = mne.read_evokeds(evoked_fname, condition=condition,
baseline=(None, 0))
evoked.crop(tmin=-50e-3, tmax=300e-3)
# Read the forward solution
forward = mne.read_forward_solution(fwd_fname)
# Read noise noise covariance matrix and regularize it
cov = mne.read_cov(cov_fname)
cov = mne.cov.regularize(cov, evoked.info, rank=None)
# Run the Gamma-MAP method with dipole output
alpha = 0.5
dipoles, residual = gamma_map(
evoked, forward, cov, alpha, xyz_same_gamma=True, return_residual=True,
return_as_dipoles=True)
Out:
Reading /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis-ave.fif ...
Read a total of 4 projection items:
PCA-v1 (1 x 102) active
PCA-v2 (1 x 102) active
PCA-v3 (1 x 102) active
Average EEG reference (1 x 60) active
Found the data of interest:
t = -199.80 ... 499.49 ms (Left visual)
0 CTF compensation matrices available
nave = 67 - aspect type = 100
Projections have already been applied. Setting proj attribute to True.
Applying baseline correction (mode: mean)
Reading forward solution from /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif...
Reading a source space...
Computing patch statistics...
Patch information added...
Distance information added...
[done]
Reading a source space...
Computing patch statistics...
Patch information added...
Distance information added...
[done]
2 source spaces read
Desired named matrix (kind = 3523) not available
Read MEG forward solution (7498 sources, 306 channels, free orientations)
Desired named matrix (kind = 3523) not available
Read EEG forward solution (7498 sources, 60 channels, free orientations)
MEG and EEG forward solutions combined
Source spaces transformed to the forward solution coordinate frame
366 x 366 full covariance (kind = 1) found.
Read a total of 4 projection items:
PCA-v1 (1 x 102) active
PCA-v2 (1 x 102) active
PCA-v3 (1 x 102) active
Average EEG reference (1 x 60) active
Computing rank from covariance with rank=None
Using tolerance 3.3e-13 (2.2e-16 eps * 305 dim * 4.8 max singular value)
Estimated rank (mag + grad): 302
MEG: rank 302 computed from 305 data channels with 3 projectors
Using tolerance 4.7e-14 (2.2e-16 eps * 59 dim * 3.6 max singular value)
Estimated rank (eeg): 58
EEG: rank 58 computed from 59 data channels with 1 projector
8 projection items activated
MAG regularization : 0.1
Created an SSP operator (subspace dimension = 3)
Computing rank from covariance with rank={'meg': 302, 'eeg': 58}
Using tolerance 2.5e-14 (2.2e-16 eps * 102 dim * 1.1 max singular value)
Estimated rank (mag): 99
MAG: rank 99 computed from 102 data channels with 3 projectors
Setting small MAG eigenvalues to zero (without PCA)
GRAD regularization : 0.1
Computing rank from covariance with rank={'meg': 302, 'eeg': 58, 'mag': 99}
Using tolerance 1.8e-13 (2.2e-16 eps * 203 dim * 3.9 max singular value)
Estimated rank (grad): 203
GRAD: rank 203 computed from 203 data channels with 0 projectors
Setting small GRAD eigenvalues to zero (without PCA)
EEG regularization : 0.1
Created an SSP operator (subspace dimension = 1)
Computing rank from covariance with rank={'meg': 302, 'eeg': 58, 'mag': 99, 'grad': 203}
Setting small EEG eigenvalues to zero (without PCA)
Converting forward solution to surface orientation
Average patch normals will be employed in the rotation to the local surface coordinates....
Converting to surface-based source orientations...
[done]
Computing inverse operator with 364 channels.
364 out of 366 channels remain after picking
Selected 364 channels
Creating the depth weighting matrix...
Applying loose dipole orientations to surface source spaces: 0.2
Whitening the forward solution.
Created an SSP operator (subspace dimension = 4)
Computing rank from covariance with rank=None
Using tolerance 3.3e-13 (2.2e-16 eps * 305 dim * 4.8 max singular value)
Estimated rank (mag + grad): 302
MEG: rank 302 computed from 305 data channels with 3 projectors
Using tolerance 4.7e-14 (2.2e-16 eps * 59 dim * 3.6 max singular value)
Estimated rank (eeg): 58
EEG: rank 58 computed from 59 data channels with 1 projector
Setting small MEG eigenvalues to zero (without PCA)
Setting small EEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Whitening data matrix.
Iteration: 0 active set size: 22494 convergence: 9.054e-01
Iteration: 8 active set size: 22236 convergence: 1.864e-01
Iteration: 9 active set size: 20715 convergence: 1.377e-01
Iteration: 10 active set size: 17181 convergence: 1.063e-01
Iteration: 11 active set size: 13002 convergence: 8.293e-02
Iteration: 12 active set size: 10110 convergence: 6.591e-02
Iteration: 13 active set size: 8313 convergence: 5.335e-02
Iteration: 14 active set size: 7209 convergence: 4.383e-02
Iteration: 15 active set size: 6354 convergence: 3.644e-02
Iteration: 16 active set size: 5700 convergence: 3.060e-02
Iteration: 17 active set size: 5145 convergence: 2.591e-02
Iteration: 18 active set size: 4701 convergence: 2.211e-02
Iteration: 19 active set size: 4296 convergence: 1.900e-02
Iteration: 20 active set size: 4005 convergence: 1.644e-02
Iteration: 21 active set size: 3693 convergence: 1.432e-02
Iteration: 22 active set size: 3432 convergence: 1.255e-02
Iteration: 23 active set size: 3153 convergence: 1.107e-02
Iteration: 24 active set size: 2958 convergence: 9.830e-03
Iteration: 25 active set size: 2745 convergence: 8.781e-03
Iteration: 26 active set size: 2562 convergence: 7.891e-03
Iteration: 27 active set size: 2397 convergence: 7.133e-03
Iteration: 28 active set size: 2226 convergence: 6.484e-03
Iteration: 29 active set size: 2070 convergence: 6.061e-03
Iteration: 30 active set size: 1935 convergence: 5.884e-03
Iteration: 31 active set size: 1839 convergence: 5.719e-03
Iteration: 32 active set size: 1725 convergence: 5.563e-03
Iteration: 33 active set size: 1614 convergence: 5.431e-03
Iteration: 34 active set size: 1512 convergence: 5.326e-03
Iteration: 35 active set size: 1437 convergence: 5.220e-03
Iteration: 36 active set size: 1350 convergence: 5.113e-03
Iteration: 37 active set size: 1245 convergence: 5.004e-03
Iteration: 38 active set size: 1200 convergence: 4.891e-03
Iteration: 39 active set size: 1128 convergence: 4.775e-03
Iteration: 40 active set size: 1077 convergence: 4.654e-03
Iteration: 41 active set size: 1011 convergence: 4.530e-03
Iteration: 42 active set size: 975 convergence: 4.402e-03
Iteration: 43 active set size: 927 convergence: 4.271e-03
Iteration: 44 active set size: 876 convergence: 4.138e-03
Iteration: 45 active set size: 840 convergence: 4.003e-03
Iteration: 46 active set size: 798 convergence: 3.867e-03
Iteration: 47 active set size: 744 convergence: 3.731e-03
Iteration: 48 active set size: 705 convergence: 3.595e-03
Iteration: 49 active set size: 687 convergence: 3.459e-03
Iteration: 50 active set size: 657 convergence: 3.326e-03
Iteration: 51 active set size: 627 convergence: 3.194e-03
Iteration: 52 active set size: 600 convergence: 3.065e-03
Iteration: 53 active set size: 588 convergence: 2.940e-03
Iteration: 54 active set size: 570 convergence: 2.817e-03
Iteration: 55 active set size: 552 convergence: 2.698e-03
Iteration: 56 active set size: 528 convergence: 2.583e-03
Iteration: 57 active set size: 501 convergence: 2.472e-03
Iteration: 58 active set size: 486 convergence: 2.365e-03
Iteration: 59 active set size: 462 convergence: 2.263e-03
Iteration: 60 active set size: 444 convergence: 2.164e-03
Iteration: 61 active set size: 417 convergence: 2.070e-03
Iteration: 62 active set size: 408 convergence: 1.980e-03
Iteration: 63 active set size: 402 convergence: 1.894e-03
Iteration: 64 active set size: 387 convergence: 1.813e-03
Iteration: 65 active set size: 375 convergence: 1.735e-03
Iteration: 66 active set size: 345 convergence: 1.660e-03
Iteration: 67 active set size: 339 convergence: 1.590e-03
Iteration: 68 active set size: 333 convergence: 1.523e-03
Iteration: 69 active set size: 330 convergence: 1.459e-03
Iteration: 70 active set size: 327 convergence: 1.398e-03
Iteration: 71 active set size: 312 convergence: 1.341e-03
Iteration: 72 active set size: 306 convergence: 1.286e-03
Iteration: 73 active set size: 297 convergence: 1.235e-03
Iteration: 74 active set size: 288 convergence: 1.185e-03
Iteration: 75 active set size: 273 convergence: 1.139e-03
Iteration: 76 active set size: 261 convergence: 1.094e-03
Iteration: 77 active set size: 252 convergence: 1.052e-03
Iteration: 78 active set size: 240 convergence: 1.012e-03
Iteration: 79 active set size: 234 convergence: 9.744e-04
Iteration: 80 active set size: 222 convergence: 9.383e-04
Iteration: 81 active set size: 219 convergence: 9.040e-04
Iteration: 82 active set size: 213 convergence: 8.714e-04
Iteration: 83 active set size: 210 convergence: 8.404e-04
Iteration: 84 active set size: 207 convergence: 8.108e-04
Iteration: 85 active set size: 198 convergence: 7.827e-04
Iteration: 86 active set size: 192 convergence: 7.559e-04
Iteration: 87 active set size: 186 convergence: 7.303e-04
Iteration: 88 active set size: 180 convergence: 7.060e-04
Iteration: 89 active set size: 174 convergence: 6.827e-04
Iteration: 91 active set size: 168 convergence: 6.393e-04
Iteration: 92 active set size: 159 convergence: 6.191e-04
Iteration: 93 active set size: 153 convergence: 5.997e-04
Iteration: 94 active set size: 150 convergence: 5.812e-04
Iteration: 97 active set size: 147 convergence: 5.302e-04
Iteration: 98 active set size: 141 convergence: 5.146e-04
Iteration: 99 active set size: 132 convergence: 4.996e-04
Iteration: 100 active set size: 129 convergence: 4.853e-04
Iteration: 102 active set size: 126 convergence: 4.582e-04
Iteration: 103 active set size: 123 convergence: 4.455e-04
Iteration: 104 active set size: 120 convergence: 4.332e-04
Iteration: 105 active set size: 117 convergence: 4.215e-04
Iteration: 106 active set size: 111 convergence: 4.101e-04
Iteration: 107 active set size: 108 convergence: 3.992e-04
Iteration: 108 active set size: 105 convergence: 3.887e-04
Iteration: 110 active set size: 102 convergence: 3.688e-04
Iteration: 113 active set size: 99 convergence: 3.414e-04
Iteration: 114 active set size: 93 convergence: 3.329e-04
Iteration: 115 active set size: 90 convergence: 3.247e-04
Iteration: 119 active set size: 87 convergence: 2.945e-04
Iteration: 120 active set size: 84 convergence: 2.875e-04
Iteration: 123 active set size: 75 convergence: 2.679e-04
Iteration: 125 active set size: 69 convergence: 2.558e-04
Iteration: 127 active set size: 66 convergence: 2.444e-04
Iteration: 133 active set size: 63 convergence: 2.139e-04
Iteration: 136 active set size: 60 convergence: 2.005e-04
Iteration: 154 active set size: 57 convergence: 1.384e-04
Iteration: 155 active set size: 54 convergence: 1.357e-04
Iteration: 168 active set size: 51 convergence: 1.054e-04
Iteration: 177 active set size: 48 convergence: 8.890e-05
Iteration: 178 active set size: 45 convergence: 8.725e-05
Iteration: 181 active set size: 42 convergence: 8.248e-05
Iteration: 211 active set size: 39 convergence: 4.750e-05
Iteration: 225 active set size: 36 convergence: 3.685e-05
Iteration: 252 active set size: 33 convergence: 2.265e-05
Iteration: 292 active set size: 30 convergence: 1.105e-05
Iteration: 303 active set size: 27 convergence: 9.077e-06
Iteration: 310 active set size: 24 convergence: 8.008e-06
Iteration: 313 active set size: 21 convergence: 7.589e-06
Iteration: 340 active set size: 18 convergence: 4.680e-06
Iteration: 342 active set size: 15 convergence: 4.515e-06
Iteration: 381 active set size: 12 convergence: 2.247e-06
Iteration: 427 active set size: 12 convergence: 9.865e-07
Convergence reached !
Explained 22.3% variance
[done]
Plot dipole activations
plot_dipole_amplitudes(dipoles)
# Plot dipole location of the strongest dipole with MRI slices
idx = np.argmax([np.max(np.abs(dip.amplitude)) for dip in dipoles])
plot_dipole_locations(dipoles[idx], forward['mri_head_t'], 'sample',
subjects_dir=subjects_dir, mode='orthoview',
idx='amplitude')
# # Plot dipole locations of all dipoles with MRI slices
# for dip in dipoles:
# plot_dipole_locations(dip, forward['mri_head_t'], 'sample',
# subjects_dir=subjects_dir, mode='orthoview',
# idx='amplitude')
Show the evoked response and the residual for gradiometers
ylim = dict(grad=[-120, 120])
evoked.pick_types(meg='grad', exclude='bads')
evoked.plot(titles=dict(grad='Evoked Response Gradiometers'), ylim=ylim,
proj=True, time_unit='s')
residual.pick_types(meg='grad', exclude='bads')
residual.plot(titles=dict(grad='Residuals Gradiometers'), ylim=ylim,
proj=True, time_unit='s')
Out:
Removing projector <Projection | PCA-v1, active : True, n_channels : 102>
Removing projector <Projection | PCA-v2, active : True, n_channels : 102>
Removing projector <Projection | PCA-v3, active : True, n_channels : 102>
Removing projector <Projection | Average EEG reference, active : True, n_channels : 60>
Removing projector <Projection | PCA-v1, active : True, n_channels : 102>
Removing projector <Projection | PCA-v2, active : True, n_channels : 102>
Removing projector <Projection | PCA-v3, active : True, n_channels : 102>
Removing projector <Projection | Average EEG reference, active : True, n_channels : 60>
Generate stc from dipoles
stc = make_stc_from_dipoles(dipoles, forward['src'])
Out:
Converting dipoles into a SourceEstimate.
[done]
View in 2D and 3D (“glass” brain like 3D plot) Show the sources as spheres scaled by their strength
scale_factors = np.max(np.abs(stc.data), axis=1)
scale_factors = 0.5 * (1 + scale_factors / np.max(scale_factors))
plot_sparse_source_estimates(
forward['src'], stc, bgcolor=(1, 1, 1),
modes=['sphere'], opacity=0.1, scale_factors=(scale_factors, None),
fig_name="Gamma-MAP")
Out:
Total number of active sources: 4
References¶
- 1
David Wipf and Srikantan Nagarajan. A unified Bayesian framework for MEG/EEG source imaging. NeuroImage, 44(3):947–966, 2009. doi:10.1016/j.neuroimage.2008.02.059.
Total running time of the script: ( 1 minutes 2.750 seconds)
Estimated memory usage: 449 MB