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 <[email protected]> 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 +@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 <[email protected]> +# +# 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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
