[MRG+2] New DICS implementation (mne-tools#5066)

* New DICS implementation

Taking @britta-wstnr 's implementation in mne-tools#4430 as basis, we implemented
a new make_dics function. Notable differences with the original is that
it no longer requires a noise-CSD matrix (unit noise is assumed instead,
we can in a later PR tackle noise normalization again if need be), and
pick_ori='max-power' is now supported.

Also note that we made an enhancement to _lcmv._reg_pinv, namely the
rcond='auto' mode, to match the original MATLAB implementation. This
mode is necessary to make DICS behave properly and may be of some use to
LCMV beamformers as well.

In tutorials/plot_dics.py the first part of a DICS connectivity analysis
tutorial is provided. This part covers using DICS for finding two
simulated sources. In the upcoming PR's for connectivity analysis this
tutorial will be extended with one-to-all connectivity and finally
all-to-all connectivity analysis.

Applying a DICS beamformer to raw time-courses does have some practical
purpose, but it will not give you spectral source power over time (only
dics_source_power and tf_dics give you that). I'm afraid the current
plot_dics_beamformer.py example gives users the wrong idea on how to use
DICS. We can later add a proper example after deliberation with Jan
Kujala.

* Choose and document DICS defaults.

* PEP8

* Fix plot_tf_dics example

* Make tf_dics less verbose. Tweak plot_tf_dics example a bit.

* Rename "mode" parameter of make_dics to "inversion"

* Address comments

* Make tf_dics parameter list more similar to the original

* DOC: Comments

* Tweak docstrings of make_dics and apply_dics_csd

* Small nitpicks

* Remove unused import

* Update what's new

doc/documentation.rst

@@ -385,6 +385,7 @@ There are also **examples**, which contain a short use-case to highlight MNE-fun
     auto_tutorials/plot_dipole_fit.rst
     auto_tutorials/plot_point_spread.rst
     auto_tutorials/plot_dipole_orientations.rst
+    auto_tutorials/plot_dics.rst
 
 
 .. raw:: html

doc/python_reference.rst

@@ -712,9 +712,10 @@ Inverse Solutions
    apply_lcmv
    apply_lcmv_epochs
    apply_lcmv_raw
-   dics
-   dics_epochs
-   dics_source_power
+   make_dics
+   apply_dics
+   apply_dics_csd
+   apply_dics_epochs
    rap_music
    tf_dics
    tf_lcmv

doc/whats_new.rst

@@ -105,6 +105,8 @@ Changelog
 
 - Add multidictionary time-frequency support to :func:`mne.inverse_sparse.tf_mixed_norm` by `Mathurin Massias`_ and `Daniel Strohmeier`_
 
+- Add new DICS implementation as :func:`mne.beamformer.make_dics`, :func:`mne.beamformer.apply_dics`, :func:`mne.beamformer.apply_dics_csd` and :func:`mne.beamformer.apply_dics_epochs`, by `Marijn van Vliet`_ and `Susanna Aro`_
+
 Bug
 ~~~
 
@@ -201,6 +203,8 @@ API
 
 - The functions ``lcmv``, ``lcmv_epochs``, and ``lcmv_raw`` are now deprecated in favor of :func:`mne.beamformer.make_lcmv` and :func:`mne.beamformer.apply_lcmv`, :func:`mne.beamformer.apply_lcmv_epochs`, and :func:`mne.beamformer.apply_lcmv_raw`, by `Britta Westner`_  
 
+- The functions ``mne.beamformer.dics``, ``mne.beamformer.dics_epochs`` and ``mne.beamformer.dics_source_power`` are now deprecated in favor of :func:`mne.beamformer.make_dics`, :func:`mne.beamformer.apply_dics`, and :func:`mne.beamformer.apply_dics_csd`, by `Marijn van Vliet`_
+
 .. _changes_0_15:
 
 Version 0.15
@@ -2652,3 +2656,5 @@ of commits):
 .. _Pierre Ablin: https://pierreablin.com
 
 .. _Oleh Kozynets: https://github.com/OlehKSS
+
+.. _Susanna Aro: https://www.linkedin.com/in/susanna-aro

examples/inverse/plot_dics_source_power.py

@@ -4,23 +4,23 @@
 =========================================
 
 Compute a Dynamic Imaging of Coherent Sources (DICS) [1]_ filter from
-single-trial activity to estimate source power for two frequencies of
-interest.
+single-trial activity to estimate source power across a frequency band.
 
 References
 ----------
 .. [1] Gross et al. Dynamic imaging of coherent sources: Studying neural
        interactions in the human brain. PNAS (2001) vol. 98 (2) pp. 694-699
 """
-# Author: Roman Goj <[email protected]>
+# Author: Marijn van Vliet <[email protected]>
+#         Roman Goj <[email protected]>
 #         Denis Engemann <[email protected]>
-#         Marijn van Vliet <[email protected]>
+#
 # License: BSD (3-clause)
-
+import numpy as np
 import mne
 from mne.datasets import sample
 from mne.time_frequency import csd_morlet
-from mne.beamformer import dics_source_power
+from mne.beamformer import make_dics, apply_dics_csd
 
 print(__doc__)
 
@@ -31,7 +31,7 @@
 subjects_dir = data_path + '/subjects'
 
 ###############################################################################
-# Read raw data
+# Reading the raw data:
 raw = mne.io.read_raw_fif(raw_fname)
 raw.info['bads'] = ['MEG 2443']  # 1 bad MEG channel
 
@@ -50,22 +50,17 @@
 # Read forward operator
 forward = mne.read_forward_solution(fname_fwd)
 
-# Computing the data and noise cross-spectral density matrices
-# The time-frequency window was chosen on the basis of spectrograms from
-# example time_frequency/plot_time_frequency.py
-# We use Morlet wavelets to estimate the CSD for two specific frequencies.
-data_csds = csd_morlet(epochs, tmin=0.04, tmax=0.15, frequencies=[18, 27])
-noise_csds = csd_morlet(epochs, tmin=-0.11, tmax=-0.001, frequencies=[18, 27])
+###############################################################################
+# Computing the cross-spectral density matrix at 4 evenly spaced frequencies
+# from 6 to 10 Hz. We use a decim value of 20 to speed up the computation in
+# this example at the loss of accuracy.
+csd = csd_morlet(epochs, tmin=0, tmax=0.5, decim=20,
+                 frequencies=np.linspace(6, 10, 4))
 
-# Compute DICS spatial filter and estimate source power
-stc = dics_source_power(epochs.info, forward, noise_csds, data_csds)
+# Compute DICS spatial filter and estimate source power.
+filters = make_dics(epochs.info, forward, csd, reg=0.5)
+stc, freqs = apply_dics_csd(csd, filters)
 
-# Plot the power maps at both frequencies
-for i, freq in enumerate(data_csds.frequencies):
-    message = 'DICS source power at %0.1f Hz' % freq
-    brain = stc.plot(surface='inflated', hemi='rh', subjects_dir=subjects_dir,
-                     time_label=message, figure=i)
-    brain.set_data_time_index(i)
-    brain.show_view('lateral')
-    # Uncomment line below to save images
-    # brain.save_image('DICS_source_power_freq_%d.png' % csd.freqs[0])
+message = 'DICS source power in the 8-12 Hz frequency band'
+brain = stc.plot(surface='inflated', hemi='rh', subjects_dir=subjects_dir,
+                 time_label=message)

examples/inverse/plot_tf_dics.py

@@ -13,7 +13,7 @@
        NeuroImage (2008) vol. 40 (4) pp. 1686-1700
 """
 # Author: Roman Goj <[email protected]>
-#         Marijn van Vliet <[email protected]>
+#
 # License: BSD (3-clause)
 
 import mne
@@ -108,15 +108,13 @@
 for freq_bin, win_length, n_fft in zip(freq_bins, win_lengths, n_ffts):
     noise_csd = csd_fourier(epochs_noise, fmin=freq_bin[0], fmax=freq_bin[1],
                             tmin=-win_length, tmax=0, n_fft=n_fft)
-    # Sum the noise CSD across the frequencies in the bin
-    noise_csd = noise_csd.sum()
-    noise_csds.append(noise_csd)
+    noise_csds.append(noise_csd.sum())
 
 # Computing DICS solutions for time-frequency windows in a label in source
 # space for faster computation, use label=None for full solution
 stcs = tf_dics(epochs, forward, noise_csds, tmin, tmax, tstep, win_lengths,
                freq_bins=freq_bins, subtract_evoked=subtract_evoked,
-               n_ffts=n_ffts, reg=0.001, label=label)
+               n_ffts=n_ffts, reg=0.05, label=label, inversion='matrix')
 
 # Plotting source spectrogram for source with maximum activity
 # Note that tmin and tmax are set to display a time range that is smaller than

mne/beamformer/__init__.py

@@ -2,5 +2,6 @@
 
 from ._lcmv import (make_lcmv, apply_lcmv, apply_lcmv_epochs, apply_lcmv_raw,
                     lcmv, lcmv_epochs, lcmv_raw, tf_lcmv)
-from ._dics import dics, dics_epochs, dics_source_power, tf_dics
+from ._dics import (make_dics, apply_dics, apply_dics_epochs, apply_dics_csd,
+                    dics, dics_epochs, dics_source_power, tf_dics)
 from ._rap_music import rap_music