This is an automated email from the git hooks/post-receive script. yoh pushed a commit to annotated tag v0.1 in repository python-mne.
commit c36a91a1cb5cd3179f36f87944af7129baaf620a Author: Alexandre Gramfort <[email protected]> Date: Wed Jul 6 11:41:56 2011 +0200 ENH : adding EEG ssp in Epochs FIX : Evoked.times now match C code TEST : better test with C code between Epochs.average and Evoked NOTE : cov estimation + SSP estimation still needs to be fixed --- mne/epochs.py | 19 ++++++- mne/fiff/__init__.py | 1 + mne/fiff/proj.py | 57 +++++++++++++++++++++ mne/fiff/raw.py | 2 +- mne/fiff/tests/data/process_raw.sh | 17 ++++-- mne/fiff/tests/data/test-cov.fif | Bin 541025 -> 547234 bytes mne/fiff/tests/data/test-nf-ave.fif | Bin 0 -> 5546473 bytes .../tests/data/{test.ave => test-no-reject.ave} | 8 +-- mne/fiff/tests/data/test.ave | 8 +-- mne/fiff/tests/test_proj.py | 39 ++++++++++++++ mne/tests/test_epochs.py | 39 ++++++++++---- 11 files changed, 165 insertions(+), 25 deletions(-) diff --git a/mne/epochs.py b/mne/epochs.py index e62a7ec..1509a83 100644 --- a/mne/epochs.py +++ b/mne/epochs.py @@ -124,6 +124,19 @@ class Epochs(object): print '%d projection items activated' % len(self.info['projs']) + # Add EEG ref reference proj + print "Adding average EEG reference projection." + eeg_sel = pick_types(self.info, meg=False, eeg=True) + eeg_names = [self.ch_names[k] for k in eeg_sel] + n_eeg = len(eeg_sel) + if n_eeg > 0: + vec = np.ones((1, n_eeg)) / n_eeg + eeg_proj_data = dict(col_names=eeg_names, row_names=None, + data=vec, nrow=1, ncol=n_eeg) + eeg_proj = dict(active=True, data=eeg_proj_data, + desc='Average EEG reference', kind=1) + self.info['projs'].append(eeg_proj) + # Create the projector proj, nproj = fiff.proj.make_projector_info(self.info) if nproj == 0: @@ -159,9 +172,11 @@ class Epochs(object): raise ValueError('No desired events found.') # Handle times + assert tmin < tmax sfreq = raw.info['sfreq'] - self.times = np.arange(int(round(tmin * sfreq)), - int(round(tmax * sfreq)), + n_times_min = int(round(tmin * float(sfreq))) + n_times_max = int(round(tmax * float(sfreq))) + self.times = np.arange(n_times_min, n_times_max + 1, dtype=np.float) / sfreq # setup epoch rejection diff --git a/mne/fiff/__init__.py b/mne/fiff/__init__.py index 67a9b03..906d1cc 100644 --- a/mne/fiff/__init__.py +++ b/mne/fiff/__init__.py @@ -12,3 +12,4 @@ from .raw import Raw, read_raw_segment, read_raw_segment_times, \ start_writing_raw, write_raw_buffer, finish_writing_raw from .pick import pick_types, pick_channels from .compensator import get_current_comp +from .proj import compute_spatial_vectors diff --git a/mne/fiff/proj.py b/mne/fiff/proj.py index def4df3..bf2b8bb 100644 --- a/mne/fiff/proj.py +++ b/mne/fiff/proj.py @@ -10,6 +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 def read_proj(fid, node): @@ -280,3 +281,59 @@ def make_projector_info(info): proj, nproj, _ = make_projector(info['projs'], info['ch_names'], info['bads']) return proj, nproj + + +def compute_spatial_vectors(epochs, n_grad=2, n_mag=2, n_eeg=2): + """Compute SSP (spatial space projection) vectors + + 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 = epochs.get_data() + data = data.swapaxes(0, 1).reshape(data.shape[1], -1) + + 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]) + + projs = [] + for n, ind, names in zip([n_grad, n_mag, n_eeg], + [grad_ind, mag_ind, eeg_ind], + [grad_names, mag_names, eeg_names]): + if n == 0: + continue + 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) + proj = dict(active=True, data=proj_data, desc='MEG %s' % k, kind=1) + projs.append(proj) + + return projs diff --git a/mne/fiff/raw.py b/mne/fiff/raw.py index e022cc2..9603c6f 100644 --- a/mne/fiff/raw.py +++ b/mne/fiff/raw.py @@ -178,7 +178,7 @@ class Raw(dict): if step is not None: raise ValueError('step needs to be 1 : %d given' % step) - if len(sel) == 0: + if sel is not None and len(sel) == 0: raise Exception("Empty channel list") return read_raw_segment(self, start=start, stop=stop, sel=sel) diff --git a/mne/fiff/tests/data/process_raw.sh b/mne/fiff/tests/data/process_raw.sh index 2db2d95..a7d5d84 100755 --- a/mne/fiff/tests/data/process_raw.sh +++ b/mne/fiff/tests/data/process_raw.sh @@ -1,13 +1,22 @@ #!/usr/bin/env bash # Generate events -mne_process_raw --raw test_raw.fif \ - --eventsout test-eve.fif +mne_process_raw --raw test_raw.fif --eventsout test-eve.fif -# Averaging +# Averaging no filter +mne_process_raw --raw test_raw.fif --projon --filteroff \ + --saveavetag -nf-ave --ave test-no-reject.ave + +# Averaging 40Hz mne_process_raw --raw test_raw.fif --lowpass 40 --projoff \ --saveavetag -ave --ave test.ave # Compute the noise covariance matrix -mne_process_raw --raw test_raw.fif --lowpass 40 --projoff \ +mne_process_raw --raw test_raw.fif --lowpass 40 --projon \ --savecovtag -cov --cov test.cov + +# Compute projection +mne_process_raw --raw test_raw.fif --events test-eve.fif --makeproj \ + --projtmin -0.2 --projtmax 0.3 --saveprojtag _proj \ + --projnmag 1 --projngrad 1 --projevent 1 \ + --projmagrej 6000 --projgradrej 5000 diff --git a/mne/fiff/tests/data/test-cov.fif b/mne/fiff/tests/data/test-cov.fif index b84054b..9484e17 100755 Binary files a/mne/fiff/tests/data/test-cov.fif and b/mne/fiff/tests/data/test-cov.fif differ diff --git a/mne/fiff/tests/data/test-nf-ave.fif b/mne/fiff/tests/data/test-nf-ave.fif new file mode 100644 index 0000000..0ea99b2 Binary files /dev/null and b/mne/fiff/tests/data/test-nf-ave.fif differ diff --git a/mne/fiff/tests/data/test.ave b/mne/fiff/tests/data/test-no-reject.ave similarity index 84% copy from mne/fiff/tests/data/test.ave copy to mne/fiff/tests/data/test-no-reject.ave index abc6533..cfabf40 100755 --- a/mne/fiff/tests/data/test.ave +++ b/mne/fiff/tests/data/test-no-reject.ave @@ -11,10 +11,10 @@ average { # # Rejection values # - gradReject 4000e-13 - magReject 4e-12 - eegReject 40e-6 - eogReject 150e-6 + # gradReject 4000e-13 + # magReject 4e-12 + # eegReject 40e-6 + # eogReject 150e-6 # # Category specifications # diff --git a/mne/fiff/tests/data/test.ave b/mne/fiff/tests/data/test.ave index abc6533..2ee6dab 100755 --- a/mne/fiff/tests/data/test.ave +++ b/mne/fiff/tests/data/test.ave @@ -11,10 +11,10 @@ average { # # Rejection values # - gradReject 4000e-13 - magReject 4e-12 - eegReject 40e-6 - eogReject 150e-6 + gradReject 4000e-13 + magReject 4e-12 + eegReject 40e-6 + eogReject 150e-6 # # Category specifications # diff --git a/mne/fiff/tests/test_proj.py b/mne/fiff/tests/test_proj.py new file mode 100644 index 0000000..b6da52a --- /dev/null +++ b/mne/fiff/tests/test_proj.py @@ -0,0 +1,39 @@ +import os.path as op + +from numpy.testing import assert_array_almost_equal + +from .. import Raw, pick_types, compute_spatial_vectors +from ..proj import make_projector +from ..open import fiff_open +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') + +def test_compute_spatial_vectors(): + """Test SSP computation + """ + event_id, tmin, tmax = 1, -0.2, 0.3 + + raw = Raw(raw_fname) + events = read_events(event_fname) + exclude = ['MEG 2443', 'EEG 053'] + picks = pick_types(raw.info, meg=True, eeg=False, stim=False, eog=False, + exclude=exclude) + epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks, + baseline=(None, 0), proj=False, + reject=dict(mag=5000e-15, grad=16e-10)) + + projs = compute_spatial_vectors(epochs, n_grad=1, n_mag=2, n_eeg=2) + + proj, nproj, U = make_projector(projs, epochs.ch_names, bads=[]) + assert nproj == 3 + assert U.shape[1] == 3 + + epochs.info['projs'] += projs + evoked = epochs.average() + evoked.save('foo.fif') + + fid, tree, _ = fiff_open(proj_fname) + projs = read_proj(fid, tree) diff --git a/mne/tests/test_epochs.py b/mne/tests/test_epochs.py index eda6db6..1503d33 100644 --- a/mne/tests/test_epochs.py +++ b/mne/tests/test_epochs.py @@ -3,29 +3,31 @@ # License: BSD (3-clause) import os.path as op -from numpy.testing import assert_array_equal +from numpy.testing import assert_array_equal, assert_array_almost_equal -import mne -from mne import fiff +from .. import fiff, Epochs, read_events raw_fname = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data', 'test_raw.fif') event_name = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data', 'test-eve.fif') +evoked_nf_name = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data', + 'test-nf-ave.fif') event_id, tmin, tmax = 1, -0.2, 0.5 raw = fiff.Raw(raw_fname) -events = mne.read_events(event_name) -picks = fiff.pick_types(raw.info, meg=True, eeg=True, stim=False, - eog=True, include=['STI 014']) +events = read_events(event_name) +picks = fiff.pick_types(raw.info, meg=True, eeg=True, stim=True, + ecg=True, eog=True, include=['STI 014']) reject = dict(grad=1000e-12, mag=4e-12, eeg=80e-6, eog=150e-6) flat = dict(grad=1e-15, mag=1e-15) + def test_read_epochs(): """Reading epochs from raw files """ - epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks, + epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0)) epochs.average() data = epochs.get_data() @@ -40,7 +42,7 @@ def test_read_epochs(): def test_reject_epochs(): """Test of epochs rejection """ - epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks, + epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0), reject=reject, flat=flat) data = epochs.get_data() @@ -56,13 +58,30 @@ def test_reject_epochs(): def test_preload_epochs(): """Test of epochs rejection """ - epochs = mne.Epochs(raw, events[:12], event_id, tmin, tmax, picks=picks, + epochs = Epochs(raw, events[:12], event_id, tmin, tmax, picks=picks, baseline=(None, 0), preload=True, reject=reject, flat=flat) data_preload = epochs.get_data() - epochs = mne.Epochs(raw, events[:12], event_id, tmin, tmax, picks=picks, + epochs = Epochs(raw, events[:12], event_id, tmin, tmax, picks=picks, baseline=(None, 0), preload=False, reject=reject, flat=flat) data_no_preload = epochs.get_data() assert_array_equal(data_preload, data_no_preload) + + +def test_comparision_with_c(): + """Test of average obtained vs C code + """ + c_evoked = fiff.Evoked(evoked_nf_name, setno=0) + epochs = Epochs(raw, events, event_id, tmin, tmax, + baseline=None, preload=True, + reject=None, flat=None) + evoked = epochs.average() + sel = fiff.pick_channels(c_evoked.ch_names, evoked.ch_names) + evoked_data = evoked.data + c_evoked_data = c_evoked.data[sel] + + assert evoked.nave == c_evoked.nave + assert_array_almost_equal(evoked_data, c_evoked_data, 10) + assert_array_almost_equal(evoked.times, c_evoked.times, 12) -- 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
