[med-svn] [python-mne] 133/353: ENH : refactor proj/SSP/PCA computation and adding support for computation on evoked data
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:44 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to tag 0.4
in repository python-mne.
commit 93074a19ada9dfe7f5155110f2ed34ab91931f5f
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Wed Apr 11 10:42:24 2012 +0200
ENH : refactor proj/SSP/PCA computation and adding support for computation on evoked data
---
bin/mne_compute_proj_ecg.py | 120 ++++++++++++++++++++++++++++++++++++++
mne/__init__.py | 3 +-
mne/fiff/proj.py | 45 ++------------
mne/proj.py | 119 +++++++++++++++++++++++++++++++++++--
mne/{fiff => }/tests/test_proj.py | 30 +++++-----
mne/utils.py | 92 +++++++++++++++++++++++++++++
setup.py | 2 +-
7 files changed, 351 insertions(+), 60 deletions(-)
diff --git a/bin/mne_compute_proj_ecg.py b/bin/mne_compute_proj_ecg.py
new file mode 100644
index 0000000..c679630
--- /dev/null
+++ b/bin/mne_compute_proj_ecg.py
@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+"""Compute SSP/PCA projections for ECG artifacts
+"""
+
+# Authors : Alexandre Gramfort, Ph.D.
+
+import mne
+
+
+def compute_proj_ecg(in_fif_fname, tmin, tmax, n_grad, n_mag, n_eeg, l_freq,
+ h_freq, average, preload, filter_length, n_jobs):
+ """Compute SSP/PCA projections for ECG artifacts
+
+ Parameters
+ ----------
+ in_fif_fname: string
+ Raw fif File
+ XXX
+ """
+ # Reading fif File
+ raw = mne.fiff.Raw(in_fif_fname, preload=preload)
+
+ if in_fif_fname.endswith('_raw.fif') or in_fif_fname.endswith('-raw.fif'):
+ prefix = in_fif_fname[:-8]
+ else:
+ prefix = in_fif_fname[:-4]
+
+ ecg_proj_fname = prefix + '_ecg_proj.fif'
+ ecg_event_fname = prefix + '_ecg-eve.fif'
+
+ print 'Running ECG SSP computation'
+
+ ecg_events, _, _ = mne.artifacts.find_ecg_events(raw)
+ print "Writing ECG events in %s" % ecg_event_fname
+ mne.write_events(ecg_event_fname, ecg_events)
+
+ print 'Computing ECG projector'
+
+ picks = mne.fiff.pick_types(raw.info, meg=True, eeg=True)
+ if l_freq is None and h_freq is not None:
+ raw.high_pass_filter(picks, h_freq, filter_length, n_jobs)
+ if l_freq is not None and h_freq is None:
+ raw.low_pass_filter(picks, h_freq, filter_length, n_jobs)
+ if l_freq is not None and h_freq is not None:
+ raw.band_pass_filter(picks, l_freq, h_freq, filter_length, n_jobs)
+
+ epochs = mne.Epochs(raw, ecg_events, None, tmin, tmax, baseline=None,
+ picks=picks)
+
+ if average:
+ evoked = epochs.average()
+ projs = mne.compute_proj_evoked(evoked, n_grad=n_grad, n_mag=n_mag,
+ n_eeg=n_eeg)
+ else:
+ projs = mne.compute_proj_epochs(epochs, n_grad=n_grad, n_mag=n_mag,
+ n_eeg=n_eeg)
+
+ print "Writing ECG projections in %s" % ecg_proj_fname
+ mne.write_proj(ecg_proj_fname, projs)
+ print 'Done.'
+
+
+if __name__ == '__main__':
+
+ from optparse import OptionParser
+
+ parser = OptionParser()
+ parser.add_option("-i", "--in", dest="raw_in",
+ help="Input raw FIF file", metavar="FILE")
+ parser.add_option("-b", "--tmin", dest="tmin",
+ help="time before event in seconds",
+ default=-0.2)
+ parser.add_option("-c", "--tmax", dest="tmax",
+ help="time before event in seconds",
+ default=0.4)
+ parser.add_option("-g", "--n-grad", dest="n_grad",
+ help="Number of SSP vectors for gradiometers",
+ default=2)
+ parser.add_option("-m", "--n-mag", dest="n_mag",
+ help="Number of SSP vectors for magnetometers",
+ default=2)
+ parser.add_option("-e", "--n-eeg", dest="n_eeg",
+ help="Number of SSP vectors for EEG",
+ default=2)
+ parser.add_option("-l", "--l-freq", dest="l_freq",
+ help="Filter low cut-off frequency in Hz",
+ default=None) # XXX
+ parser.add_option("-t", "--h-freq", dest="h_freq",
+ help="Filter high cut-off frequency in Hz",
+ default=None) # XXX
+ parser.add_option("-p", "--preload", dest="preload",
+ help="Temporary file used during computaion",
+ default='tmp.mmap')
+ parser.add_option("-a", "--average", dest="average", action="store_true",
+ help="Compute SSP after averaging",
+ default=False)
+ parser.add_option("-f", "--filter-length", dest="filter_length",
+ help="Number of SSP vectors for EEG",
+ default=2048)
+ parser.add_option("-j", "--n-jobs", dest="n_jobs",
+ help="Number of jobs to run in parallel",
+ default=1)
+
+ options, args = parser.parse_args()
+
+ raw_in = options.raw_in
+ tmin = options.tmin
+ tmax = options.tmax
+ n_grad = options.n_grad
+ n_mag = options.n_mag
+ n_eeg = options.n_eeg
+ l_freq = options.l_freq
+ h_freq = options.h_freq
+ average = options.average
+ preload = options.preload
+ filter_length = options.filter_length
+ n_jobs = options.n_jobs
+
+ compute_proj_ecg(raw_in, tmin, tmax, n_grad, n_mag, n_eeg, l_freq, h_freq,
+ average, preload, filter_length, n_jobs)
diff --git a/mne/__init__.py b/mne/__init__.py
index 4c2cfca..7ab8ac9 100644
--- a/mne/__init__.py
+++ b/mne/__init__.py
@@ -20,7 +20,8 @@ from .label import label_time_courses, read_label, label_sign_flip, \
write_label, stc_to_label
from .misc import parse_config, read_reject_parameters
from .transforms import transform_coordinates
-from .proj import read_proj
+from .proj import read_proj, write_proj, compute_proj_epochs, \
+ compute_proj_evoked
from . import fiff
from . import artifacts
from . import stats
diff --git a/mne/fiff/proj.py b/mne/fiff/proj.py
index 81032e7..f8ee02c 100644
--- a/mne/fiff/proj.py
+++ b/mne/fiff/proj.py
@@ -10,7 +10,7 @@ from scipy import linalg
from .tree import dir_tree_find
from .constants import FIFF
from .tag import find_tag
-from .pick import pick_types
+from ..utils import deprecated
class Projection(dict):
@@ -309,6 +309,7 @@ def make_projector_info(info):
return proj, nproj
+ at deprecated("Use mne.compute_proj_epochs")
def compute_spatial_vectors(epochs, n_grad=2, n_mag=2, n_eeg=2):
"""Compute SSP (spatial space projection) vectors
@@ -328,43 +329,5 @@ def compute_spatial_vectors(epochs, n_grad=2, n_mag=2, n_eeg=2):
projs: list
List of projection vectors
"""
- data = sum(np.dot(e, e.T) for e in epochs) # compute data covariance
-
- mag_ind = pick_types(epochs.info, meg='mag')
- grad_ind = pick_types(epochs.info, meg='grad')
- eeg_ind = pick_types(epochs.info, meg=False, eeg=True)
-
- if (n_grad > 0) and len(grad_ind) == 0:
- print "No gradiometers found. Forcing n_grad to 0"
- n_grad = 0
- if (n_mag > 0) and len(mag_ind) == 0:
- print "No magnetometers found. Forcing n_mag to 0"
- n_mag = 0
- if (n_eeg > 0) and len(eeg_ind) == 0:
- print "No EEG channels found. Forcing n_eeg to 0"
- n_eeg = 0
-
- grad_names, mag_names, eeg_names = ([epochs.ch_names[k] for k in ind]
- for ind in [grad_ind, mag_ind, eeg_ind])
-
- event_id = epochs.event_id
- projs = []
- for n, ind, names, desc in zip([n_grad, n_mag, n_eeg],
- [grad_ind, mag_ind, eeg_ind],
- [grad_names, mag_names, eeg_names],
- ['planar', 'axial', 'eeg']):
- if n == 0:
- continue
- data_ind = data[ind][:, ind]
- U = linalg.svd(data_ind, full_matrices=False,
- overwrite_a=True)[0][:, :n]
- for k, u in enumerate(U.T):
- proj_data = dict(col_names=names, row_names=None,
- data=u[np.newaxis, :], nrow=1, ncol=u.size)
- this_desc = "%s-%-d-%-.3f-%-.3f-PCA-%02d" % (desc, event_id,
- epochs.tmin, epochs.tmax, k + 1)
- print "Adding projection: %s" % this_desc
- proj = dict(active=True, data=proj_data, desc=this_desc, kind=1)
- projs.append(proj)
-
- return projs
+ import mne # XXX : ugly due to circular mess in imports
+ return mne.compute_proj_epochs(epochs, n_grad, n_mag, n_eeg)
diff --git a/mne/proj.py b/mne/proj.py
index 2bd4172..6c3c5b6 100644
--- a/mne/proj.py
+++ b/mne/proj.py
@@ -2,8 +2,11 @@
#
# License: BSD (3-clause)
-from .fiff import fiff_open
-from .fiff.proj import read_proj as fiff_read_proj
+import numpy as np
+from scipy import linalg
+
+from . import fiff
+from .fiff.pick import pick_types
def read_proj(fname):
@@ -19,6 +22,114 @@ def read_proj(fname):
projs: list
The list of projection vectors.
"""
- fid, tree, _ = fiff_open(fname)
- projs = fiff_read_proj(fid, tree)
+ fid, tree, _ = fiff.fiff_open(fname)
+ projs = fiff.proj.read_proj(fid, tree)
+ return projs
+
+
+def write_proj(fname, projs):
+ """Write projections to a FIF file.
+
+ Parameters
+ ----------
+ fname: string
+ The name of file containing the projections vectors.
+
+ projs: list
+ The list of projection vectors.
+ """
+ fid = fiff.write.start_file(fname)
+ fiff.proj.write_proj(fid, projs)
+ fiff.write.end_file(fid)
+
+
+def _compute_proj(data, info, n_grad, n_mag, n_eeg, desc_prefix):
+
+ mag_ind = pick_types(info, meg='mag')
+ grad_ind = pick_types(info, meg='grad')
+ eeg_ind = pick_types(info, meg=False, eeg=True)
+
+ if (n_grad > 0) and len(grad_ind) == 0:
+ print "No gradiometers found. Forcing n_grad to 0"
+ n_grad = 0
+ if (n_mag > 0) and len(mag_ind) == 0:
+ print "No magnetometers found. Forcing n_mag to 0"
+ n_mag = 0
+ if (n_eeg > 0) and len(eeg_ind) == 0:
+ print "No EEG channels found. Forcing n_eeg to 0"
+ n_eeg = 0
+
+ ch_names = info['ch_names']
+ grad_names, mag_names, eeg_names = ([ch_names[k] for k in ind]
+ for ind in [grad_ind, mag_ind, eeg_ind])
+
+ projs = []
+ for n, ind, names, desc in zip([n_grad, n_mag, n_eeg],
+ [grad_ind, mag_ind, eeg_ind],
+ [grad_names, mag_names, eeg_names],
+ ['planar', 'axial', 'eeg']):
+ if n == 0:
+ continue
+ data_ind = data[ind][:, ind]
+ U = linalg.svd(data_ind, full_matrices=False,
+ overwrite_a=True)[0][:, :n]
+ for k, u in enumerate(U.T):
+ proj_data = dict(col_names=names, row_names=None,
+ data=u[np.newaxis, :], nrow=1, ncol=u.size)
+ this_desc = "%s-%s-PCA-%02d" % (desc, desc_prefix, k + 1)
+ print "Adding projection: %s" % this_desc
+ proj = dict(active=True, data=proj_data, desc=this_desc, kind=1)
+ projs.append(proj)
+
return projs
+
+
+def compute_proj_epochs(epochs, n_grad=2, n_mag=2, n_eeg=2):
+ """Compute SSP (spatial space projection) vectors on Epochs
+
+ Parameters
+ ----------
+ epochs: instance of Epochs
+ The epochs containing the artifact
+ n_grad: int
+ Number of vectors for gradiometers
+ n_mag: int
+ Number of vectors for gradiometers
+ n_eeg: int
+ Number of vectors for gradiometers
+
+ Returns
+ -------
+ projs: list
+ List of projection vectors
+ """
+ data = sum(np.dot(e, e.T) for e in epochs) # compute data covariance
+ event_id = epochs.event_id
+ if event_id is None:
+ event_id = 0
+ desc_prefix = "%-d-%-.3f-%-.3f" % (event_id, epochs.tmin, epochs.tmax)
+ return _compute_proj(data, epochs.info, n_grad, n_mag, n_eeg, desc_prefix)
+
+
+def compute_proj_evoked(evoked, n_grad=2, n_mag=2, n_eeg=2):
+ """Compute SSP (spatial space projection) vectors on Evoked
+
+ Parameters
+ ----------
+ evoked: instance of Evoked
+ The Evoked obtained by averaging the artifact
+ n_grad: int
+ Number of vectors for gradiometers
+ n_mag: int
+ Number of vectors for gradiometers
+ n_eeg: int
+ Number of vectors for gradiometers
+
+ Returns
+ -------
+ projs: list
+ List of projection vectors
+ """
+ data = np.dot(evoked.data, evoked.data.T) # compute data covariance
+ desc_prefix = "%-.3f-%-.3f" % (evoked.times[0], evoked.times[-1])
+ return _compute_proj(data, evoked.info, n_grad, n_mag, n_eeg, desc_prefix)
diff --git a/mne/fiff/tests/test_proj.py b/mne/tests/test_proj.py
similarity index 69%
rename from mne/fiff/tests/test_proj.py
rename to mne/tests/test_proj.py
index 032fc27..d9968cb 100644
--- a/mne/fiff/tests/test_proj.py
+++ b/mne/tests/test_proj.py
@@ -4,17 +4,19 @@ from nose.tools import assert_true
import numpy as np
from numpy.testing import assert_array_almost_equal
-from .. import Raw, pick_types, compute_spatial_vectors
-from ..proj import make_projector, read_proj
-from ..open import fiff_open
-from ... import read_events, Epochs
+from ..fiff import Raw, pick_types
+from .. import compute_proj_epochs, compute_proj_evoked
+from ..fiff.proj import make_projector
+from ..proj import read_proj
+from .. import read_events, Epochs
-raw_fname = op.join(op.dirname(__file__), 'data', 'test_raw.fif')
-event_fname = op.join(op.dirname(__file__), 'data', 'test-eve.fif')
-proj_fname = op.join(op.dirname(__file__), 'data', 'test_proj.fif')
+base_dir = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data')
+raw_fname = op.join(base_dir, 'test_raw.fif')
+event_fname = op.join(base_dir, 'test-eve.fif')
+proj_fname = op.join(base_dir, 'test_proj.fif')
-def test_compute_spatial_vectors():
+def test_compute_proj():
"""Test SSP computation"""
event_id, tmin, tmax = 1, -0.2, 0.3
@@ -27,10 +29,10 @@ def test_compute_spatial_vectors():
epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
baseline=None, proj=False)
- projs = compute_spatial_vectors(epochs, n_grad=1, n_mag=1, n_eeg=0)
+ evoked = epochs.average()
+ projs = compute_proj_epochs(epochs, n_grad=1, n_mag=1, n_eeg=0)
- fid, tree, _ = fiff_open(proj_fname)
- projs2 = read_proj(fid, tree)
+ projs2 = read_proj(proj_fname)
for k, (p1, p2) in enumerate(zip(projs, projs2)):
assert_true(p1['desc'] == p2['desc'])
@@ -57,5 +59,7 @@ def test_compute_spatial_vectors():
evoked = epochs.average()
evoked.save('foo.fif')
- fid, tree, _ = fiff_open(proj_fname)
- projs = read_proj(fid, tree)
+ projs = read_proj(proj_fname)
+
+ projs_evoked = compute_proj_evoked(evoked, n_grad=1, n_mag=1, n_eeg=0)
+ # XXX : test something
diff --git a/mne/utils.py b/mne/utils.py
new file mode 100644
index 0000000..31f826a
--- /dev/null
+++ b/mne/utils.py
@@ -0,0 +1,92 @@
+"""Some utility functions"""
+
+# Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+import warnings
+
+
+# Following deprecated class copied from scikit-learn
+
+
+class deprecated(object):
+ """Decorator to mark a function or class as deprecated.
+
+ Issue a warning when the function is called/the class is instantiated and
+ adds a warning to the docstring.
+
+ The optional extra argument will be appended to the deprecation message
+ and the docstring. Note: to use this with the default value for extra, put
+ in an empty of parentheses:
+
+ >>> from sklearn.utils import deprecated
+ >>> deprecated() # doctest: +ELLIPSIS
+ <sklearn.utils.deprecated object at ...>
+
+ >>> @deprecated()
+ ... def some_function(): pass
+ """
+
+ # Adapted from http://wiki.python.org/moin/PythonDecoratorLibrary,
+ # but with many changes.
+
+ def __init__(self, extra=''):
+ """
+ Parameters
+ ----------
+ extra: string
+ to be added to the deprecation messages
+
+ """
+ self.extra = extra
+
+ def __call__(self, obj):
+ if isinstance(obj, type):
+ return self._decorate_class(obj)
+ else:
+ return self._decorate_fun(obj)
+
+ def _decorate_class(self, cls):
+ msg = "Class %s is deprecated" % cls.__name__
+ if self.extra:
+ msg += "; %s" % self.extra
+
+ # FIXME: we should probably reset __new__ for full generality
+ init = cls.__init__
+
+ def wrapped(*args, **kwargs):
+ warnings.warn(msg, category=DeprecationWarning)
+ return init(*args, **kwargs)
+ cls.__init__ = wrapped
+
+ wrapped.__name__ = '__init__'
+ wrapped.__doc__ = self._update_doc(init.__doc__)
+ wrapped.deprecated_original = init
+
+ return cls
+
+ def _decorate_fun(self, fun):
+ """Decorate function fun"""
+
+ msg = "Function %s is deprecated" % fun.__name__
+ if self.extra:
+ msg += "; %s" % self.extra
+
+ def wrapped(*args, **kwargs):
+ warnings.warn(msg, category=DeprecationWarning)
+ return fun(*args, **kwargs)
+
+ wrapped.__name__ = fun.__name__
+ wrapped.__dict__ = fun.__dict__
+ wrapped.__doc__ = self._update_doc(fun.__doc__)
+
+ return wrapped
+
+ def _update_doc(self, olddoc):
+ newdoc = "DEPRECATED"
+ if self.extra:
+ newdoc = "%s: %s" % (newdoc, self.extra)
+ if olddoc:
+ newdoc = "%s\n\n%s" % (newdoc, olddoc)
+ return newdoc
diff --git a/setup.py b/setup.py
index e1f513b..b6bf6e4 100755
--- a/setup.py
+++ b/setup.py
@@ -60,4 +60,4 @@ if __name__ == "__main__":
'mne.layouts',
'mne.time_frequency', 'mne.time_frequency.tests'],
scripts=['bin/mne_clean_eog_ecg.py', 'bin/mne_flash_bem_model.py',
- 'bin/mne_surf2bem.py'])
+ 'bin/mne_surf2bem.py', 'bin/mne_compute_proj_ecg.py'])
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git
More information about the debian-med-commit
mailing list