CSD with working example and tests

added some tests WIP CSD MVP example working with example weird cherry-picking error again Update examples/preprocessing/plot_csd.py Co-Authored-By: Alexandre Gramfort <[email protected]> Update examples/preprocessing/plot_csd.py Co-Authored-By: Alexandre Gramfort <[email protected]> Update examples/preprocessing/plot_csd.py Co-Authored-By: Alexandre Gramfort <[email protected]> Update mne/preprocessing/_csd.py Co-Authored-By: Alexandre Gramfort <[email protected]> Update examples/preprocessing/plot_csd.py Co-Authored-By: Alexandre Gramfort <[email protected]> typos/small fixes incorporated Eric comments more Eric comments refactoring core comp per Eric realized I did that wrong Eric comments 3 plot fix Fix julian date conversions (mne-tools#6911) * Fix julian date conversions * remove tabs * use UTC and private functions. closes mne-tools#6901 [skip travis][skip azp] (mne-tools#6914) ENH: Add source space SNR estimate (mne-tools#6908) * source-snr-edit * FIX move estimate_snr to be a method * ENH: name change * Make it work * SNR works, now clean up code * fixed unnecessary stc copy, returns stc object * added plot source space SNR tutorial file * updated author names * fixed period for pydocsty and suggested documenation edit * edits suggested from @drammock * made edits suggested by reviewers * made edits suggested by reviewers * N -> n_channels * Fixed formatting of Reference in comment * fixed formatting in doc and in code * fixed snr_stc creation * further minor edits * fixed formatting issues and colormap choice in eg * removed colon after Reference * Corrected SNR calculation to take unwhitened sensor signal b instead of whitened, updated tutorial file * Minor changes to formatting * FIX: transparent * FIX: Rename * FIX? Added dSPM (mne-tools#6913) * Added dSPM * removed term signal covariance * Minor modification of dSPM definition * Corrected according to Mattis suggestions [MRG]: read *.dat electrode position files (mne-tools#6909) * ADD: channels.read_dig_dat() * DOC * DOC * TEST * DOC, __all__ * DOC: intro * DOC * pep & py3.5 compatibility * py3.5 MRG: allow users to annotate epochs with colors (mne-tools#6855) * ENH: copy over diff from autoreject * Make it work * More fixes * More fixes * Make horizontal patch * Cleanup * TST: add tests * Address Eric's comments * Address alex comments * DOC: update whats new * FIX travis + circle * FIX For ICA plot Fix infomax (mne-tools#6921) * fix extended infomax * update whats new * cleaner way, using setdefault DOC: Document authors better [ci skip] FIX: Comments and order umlauts MRG, MAINT: Simplify code for vol ravel (mne-tools#6920) * MAINT: Simplify code for vol ravel * STY: Flake ENH: Convert NIRS raw data to optical density (mne-tools#6827) * Initial framework for fnirs raw to optical density conversion * Add optical density to docs * Update mne/preprocessing/tests/test_optical_density.py Co-Authored-By: Eric Larson <[email protected]> * Update mne/preprocessing/tests/test_optical_density.py Co-Authored-By: Eric Larson <[email protected]> * Import testing * Put optical density docs in preprocessing section * Implement optical density conversion * Fix optical density code style issues * Add plotting of optical density data * Add defaults for optical density plotting * Update mne/preprocessing/optical_density.py Co-Authored-By: Eric Larson <[email protected]> * Handle negative intensities in optical density conversion * Add in -1* that was removed in previous commit * Add tests and fix imports for optical density * Add tests for optical density conversion * Use assert_allclose for testing optical density * Update test_optical_density.py * Update mne/preprocessing/_optical_density.py Co-Authored-By: Alexandre Gramfort <[email protected]> * Do not operate in place for optical density * Force data to load and modify test to include preload=False condition
alexrockhill · Oct 28, 2019 · ea59670 · ea59670
1 parent 799c2ec
commit ea59670
Show file tree

Hide file tree

Showing 15 changed files with 532 additions and 13 deletions.
diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc
@@ -42,15 +42,19 @@ Changelog
 
 - Add reader for ``*.dat`` electrode position files :func:`mne.channels.read_dig_dat` by `Christian Brodbeck`_
 
+<<<<<<< HEAD
 - Improved :ref:`limo-dataset` usage and :ref:`example <ex-limo-data>` for usage of :func:`mne.stats.linear_regression` by `Jose Alanis`_
 
+=======
+>>>>>>> CSD with working example and tests
 - For KIT systems without built-in layout, :func:`mne.channels.find_layout` now falls back on an automatically generated layout, by `Christian Brodbeck`_
 
 - :meth:`mne.Epochs.plot` now takes a ``epochs_colors`` parameter to color specific epoch segments by `Mainak Jas`_
 
 Bug
 ~~~
 
+<<<<<<< HEAD
 - Fix :meth:`mne.Epochs.shift_time` and :meth:`mne.Evoked.shift_time` to return the modified :class:`~mne.Epochs` or :class:`~mne.Evoked` instance (instead of ``None``) by `Daniel McCloy`_.
 
 - Fix bug in :class:`~mne.preprocessing.ICA` where requesting extended infomax via ``fit_params={'extended': True}`` was overridden, by `Daniel McCloy`_.
@@ -59,6 +63,10 @@ Bug
 
 - Fix bug in :func:`mne.preprocessing.mark_flat` where acquisition skips were not handled proeprly, by `Eric Larson`_
 
+=======
+- Fix bug in :class:`~mne.preprocessing.ICA` where requesting extended infomax via ``fit_params={'extended': True}`` was overridden, by `Daniel McCloy`_.
+
+>>>>>>> CSD with working example and tests
 - Fix TAL channel parsing (annotations) for EDF-D files by `Clemens Brunner`_
 
 - Fix bug with :func:`mne.viz.plot_dipole_locations` when plotting in head coordinates by `Eric Larson`_

diff --git a/doc/python_reference.rst b/doc/python_reference.rst
@@ -343,6 +343,7 @@ Projections:
 
    ICA
    Xdawn
+   compute_current_source_density
    compute_proj_ecg
    compute_proj_eog
    create_ecg_epochs

diff --git a/examples/inverse/plot_mixed_norm_inverse.py b/examples/inverse/plot_mixed_norm_inverse.py
@@ -7,6 +7,21 @@
 L0.5/L2 is done with irMxNE which allows for sparser
 source estimates with less amplitude bias due to the non-convexity
 of the L0.5/L2 mixed norm penalty.
+<<<<<<< HEAD
+=======
+
+References
+----------
+.. [1] Gramfort A., Kowalski M. and Hämäläinen, M.
+   "Mixed-norm estimates for the M/EEG inverse problem using accelerated
+   gradient methods", Physics in Medicine and Biology, 2012.
+   https://doi.org/10.1088/0031-9155/57/7/1937.
+
+.. [2] Strohmeier D., Haueisen J., and Gramfort A.
+   "Improved MEG/EEG source localization with reweighted mixed-norms",
+   4th International Workshop on Pattern Recognition in Neuroimaging,
+   Tuebingen, 2014. 10.1109/PRNI.2014.6858545
+>>>>>>> CSD with working example and tests
 """
 # Author: Alexandre Gramfort <[email protected]>
 #         Daniel Strohmeier <[email protected]>

diff --git a/examples/preprocessing/plot_csd.py b/examples/preprocessing/plot_csd.py
@@ -0,0 +1,99 @@
+# -*- coding: utf-8 -*-
+"""
+.. _ex-: plot_csd
+
+=================================================================
+Analyze Data Using Current Source Density (CSD)
+=================================================================
+
+This script shows an example of how to use CSD. CSD
+takes the spatial Laplacian of the sensor signal (derivative in both
+x and y). It does what a planar gradiometer does in MEG. Spatial derivative
+reduces point spread. CSD transformed data have a sharper or more distinct topography,
+reducing the negative impact of volume conduction.
+"""
+# Authors: Alex Rockhill <[email protected]>
+#
+# https://gist.github.com/sherdim/28b069c8ebca17073f24322e8e721e1f
+
+# License: BSD (3-clause)
+import matplotlib.pyplot as plt
+
+import numpy as np
+import mne
+from mne.datasets import sample
+
+from mne.preprocessing import compute_current_source_density
+
+print(__doc__)
+
+###############################################################################
+# Load sample subject data
+data_path = sample.data_path()
+
+raw = mne.io.read_raw_fif(data_path + '/MEG/sample/sample_audvis_raw.fif',
+                          preload=True)
+raw = raw.pick_types(meg=False, eeg=True, eog=True, ecg=True, stim=True,
+                     exclude=raw.info['bads'])
+event_id = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3,
+            'visual/right': 4, 'smiley': 5, 'button': 32}
+events = mne.find_events(raw)
+epochs = mne.Epochs(raw, events, event_id=event_id, tmin=-0.2, tmax=.5,
+                    preload=True)
+evo = epochs.average()
+
+###############################################################################
+# Let's look at the topography of CSD compared to average
+
+fig, axes = plt.subplots(1, 6)
+times = np.array([-0.1, 0., 0.05, 0.1, 0.15])
+evo.plot_topomap(times=times, cmap='Spectral_r', axes=axes[:5],
+                 outlines='skirt', contours=4, time_unit='s',
+                 colorbar=True, show=False, title='Average Reference')
+fig, axes = plt.subplots(1, 6)
+csd_evo = compute_current_source_density(evo, copy=True)
+csd_evo.plot_topomap(times=times, axes=axes[:5], cmap='Spectral_r',
+                     outlines='skirt', contours=4, time_unit='s',
+                     colorbar=True, title='Current Source Density')
+
+###############################################################################
+# Let's add the evoked plot
+
+evo.plot_joint(title='Average Reference', show=False)
+csd_evo.plot_joint(title='Current Source Density')
+
+###############################################################################
+# Let's look at the effect of smoothing and spline flexibility
+
+fig, ax = plt.subplots(4, 4)
+fig.set_size_inches(10, 10)
+for i, lambda2 in enumerate([1e-7, 1e-5, 1e-3, 0]):
+    for j, m in enumerate([5, 4, 3, 2]):
+        this_csd_evo = compute_current_source_density(evo, stiffness=m,
+                                                      lambda2=lambda2,
+                                                      copy=True)
+        this_csd_evo.plot_topomap(0.1, axes=ax[i, j],
+                                  outlines='skirt', contours=4, time_unit='s',
+                                  colorbar=False, show=False)
+        ax[i, j].set_title('m=%i, lambda=%s' % (m, lambda2))
+
+plt.show()
+
+# References
+# ----------
+#
+# Perrin et al. (1989, 1990)
+#
+# \MATLAB\eeglab14_1_1b\functions\popfunc\eeg_laplac.m
+#    Juan Sebastian Gonzalez, DFI / Universidad de Chile
+#
+# MATLAB\eeglab14_1_1b\plugins\PrepPipeline0.55.3\utilities\private\spherical_interpolate.m
+#    Copyright 2009-     by Jason D.R. Farquhar ([email protected])
+#
+# Wang K, Begleiter, 1999
+#
+# MATLAB\eeglab14_1_1b\plugins\erplab6.1.3\functions\csd\current_source_density.m
+#    Copyright (C) 2003 by Jürgen Kayser (Email: [email protected])
+#    (published in appendix of Kayser J, Tenke CE,
+#     Clin Neurophysiol 2006;117(2):348-368)
+#
diff --git a/mne/channels/interpolation.py b/mne/channels/interpolation.py
@@ -14,6 +14,25 @@
 from ..forward import _map_meg_channels
 
 
+def _calc_h(cosang, stiffness=4, num_lterms=50):
+    """Calculate spherical spline h function between points on a sphere.
+
+    Parameters
+    ----------
+    cosang : array-like | float
+        cosine of angles between pairs of points on a spherical surface. This
+        is equivalent to the dot product of unit vectors.
+    stiffness : float
+        stiffnes of the spline. Also referred to as `m`.
+    num_lterms : int
+        number of Legendre terms to evaluate.
+    """
+    factors = [(2 * n + 1) /
+               (n ** (stiffness - 1) * (n + 1) ** (stiffness - 1) * 4 * np.pi)
+               for n in range(1, num_lterms + 1)]
+    return legval(cosang, [0] + factors)
+
+
 def _calc_g(cosang, stiffness=4, num_lterms=50):
     """Calculate spherical spline g function between points on a sphere.
 

diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py
@@ -238,7 +238,11 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
     path = _get_path(path, key, name)
     # To update the testing or misc dataset, push commits, then make a new
     # release on GitHub. Then update the "releases" variable:
+<<<<<<< HEAD
     releases = dict(testing='0.75', misc='0.5')
+=======
+    releases = dict(testing='0.74', misc='0.3')
+>>>>>>> CSD with working example and tests
     # And also update the "md5_hashes['testing']" variable below.
 
     # To update any other dataset, update the data archive itself (upload
@@ -320,7 +324,11 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
         sample='fc2d5b9eb0a144b1d6ba84dc3b983602',
         somato='f08f17924e23c57a751b3bed4a05fe02',
         spm='9f43f67150e3b694b523a21eb929ea75',
+<<<<<<< HEAD
         testing='51d2ca16c449bbb27a54f53590283f16',
+=======
+        testing='5275c0989dff8007811103b951dfb487',
+>>>>>>> CSD with working example and tests
         multimodal='26ec847ae9ab80f58f204d09e2c08367',
         opm='370ad1dcfd5c47e029e692c85358a374',
         visual_92_categories=['74f50bbeb65740903eadc229c9fa759f',

diff --git a/mne/io/meas_info.py b/mne/io/meas_info.py
@@ -2028,6 +2028,9 @@ def _bad_chans_comp(info, ch_names):
         # should this be thought of as a bug?
         return False, []
 
+    if 'csd' in info['comps']:  # added for csd 0.20
+        return True, []
+
     # only include compensation channels that would affect selected channels
     ch_names_s = set(ch_names)
     comp_names = []

diff --git a/mne/io/pick.py b/mne/io/pick.py
@@ -477,17 +477,18 @@ def pick_info(info, sel=(), copy=True, verbose=None):
     info['bads'] = [ch for ch in info['bads'] if ch in info['ch_names']]
 
     if 'comps' in info:
-        comps = deepcopy(info['comps'])
-        for c in comps:
-            row_idx = [k for k, n in enumerate(c['data']['row_names'])
-                       if n in info['ch_names']]
-            row_names = [c['data']['row_names'][i] for i in row_idx]
-            rowcals = c['rowcals'][row_idx]
-            c['rowcals'] = rowcals
-            c['data']['nrow'] = len(row_names)
-            c['data']['row_names'] = row_names
-            c['data']['data'] = c['data']['data'][row_idx]
-        info['comps'] = comps
+        if 'csd' not in info['comps']:  # added for csd 0.20
+            comps = deepcopy(info['comps'])
+            for c in comps:
+                row_idx = [k for k, n in enumerate(c['data']['row_names'])
+                           if n in info['ch_names']]
+                row_names = [c['data']['row_names'][i] for i in row_idx]
+                rowcals = c['rowcals'][row_idx]
+                c['rowcals'] = rowcals
+                c['data']['nrow'] = len(row_names)
+                c['data']['row_names'] = row_names
+                c['data']['data'] = c['data']['data'][row_idx]
+            info['comps'] = comps
     info._check_consistency()
     return info
 

diff --git a/mne/preprocessing/__init__.py b/mne/preprocessing/__init__.py
@@ -21,4 +21,8 @@
 from .stim import fix_stim_artifact
 from .maxwell import maxwell_filter
 from .xdawn import Xdawn
+<<<<<<< HEAD
+=======
+from ._csd import compute_current_source_density
+>>>>>>> CSD with working example and tests
 from ._optical_density import optical_density