[med-svn] [python-mne] 309/376: 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
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:23:12 UTC 2015
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 <alexandre.gramfort at inria.fr>
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
More information about the debian-med-commit
mailing list