[med-svn] [python-mne] 175/376: ENH: refactoring covariance estimation. TODO: reject when cov estimation
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:31 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 25a16a735324064cfad8ca858d0c334cf3851f1c
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Wed Mar 30 18:09:07 2011 -0400
ENH: refactoring covariance estimation. TODO: reject when cov estimation
---
examples/plot_estimate_covariance_matrix.py | 27 +---
examples/plot_from_raw_to_epochs_to_evoked.py | 4 +-
examples/plot_minimum_norm_estimate.py | 3 +-
examples/plot_read_epochs.py | 2 +-
examples/plot_read_noise_covariance_matrix.py | 4 +-
examples/plot_whitened_evoked_data.py | 7 +-
.../plot_cluster_1samp_test_time_frequency.py | 11 +-
examples/stats/plot_cluster_stats_evoked.py | 6 +-
.../stats/plot_cluster_stats_time_frequency.py | 11 +-
examples/stats/plot_sensor_permutation_test.py | 4 +-
examples/time_frequency/plot_time_frequency.py | 9 +-
mne/__init__.py | 5 +-
mne/cov.py | 163 ++++++++++++++++-----
mne/epochs.py | 10 +-
mne/fiff/tests/data/test-cov.fif | Bin 541025 -> 541025 bytes
mne/fiff/tests/data/test.cov | 12 +-
mne/fiff/tests/data/test_empty_room.cov | 44 ++++++
mne/fiff/tests/data/test_erm-cov.fif | Bin 0 -> 541025 bytes
mne/misc.py | 39 +++--
mne/tests/test_cov.py | 48 +++---
mne/tests/test_epochs.py | 7 +
21 files changed, 282 insertions(+), 134 deletions(-)
diff --git a/examples/plot_estimate_covariance_matrix.py b/examples/plot_estimate_covariance_matrix.py
index 65c62bc..5b41bb5 100644
--- a/examples/plot_estimate_covariance_matrix.py
+++ b/examples/plot_estimate_covariance_matrix.py
@@ -20,32 +20,13 @@ fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
raw = fiff.Raw(fname)
# Set up pick list: MEG + STI 014 - bad channels
-want_meg = True
-want_eeg = False
-want_stim = False
-
-picks = fiff.pick_types(raw.info, meg=want_meg, eeg=want_eeg,
- stim=want_stim, exclude=raw.info['bads'])
-
-print "Number of picked channels : %d" % len(picks)
-
-full_cov = mne.Covariance(kind='full')
-full_cov.estimate_from_raw(raw, picks=picks)
-print full_cov
-
-diagonal_cov = mne.Covariance(kind='diagonal')
-diagonal_cov.estimate_from_raw(raw, picks=picks)
-print diagonal_cov
+cov = mne.noise_covariance_segment(raw)
+print cov
###############################################################################
# Show covariance
import pylab as pl
-pl.figure(figsize=(8, 4))
-pl.subplot(1, 2, 1)
-pl.imshow(full_cov.data, interpolation="nearest")
+pl.figure()
+pl.imshow(cov.data, interpolation="nearest")
pl.title('Full covariance matrix')
-pl.subplot(1, 2, 2)
-pl.imshow(diagonal_cov.data, interpolation="nearest")
-pl.title('Diagonal covariance matrix')
-pl.subplots_adjust(0.06, 0.02, 0.98, 0.94, 0.16, 0.26)
pl.show()
diff --git a/examples/plot_from_raw_to_epochs_to_evoked.py b/examples/plot_from_raw_to_epochs_to_evoked.py
index 629f70b..ac77f86 100644
--- a/examples/plot_from_raw_to_epochs_to_evoked.py
+++ b/examples/plot_from_raw_to_epochs_to_evoked.py
@@ -36,12 +36,12 @@ include = [] # or stim channels ['STI 014']
exclude = raw.info['bads'] + ['EEG 053'] # bads + 1 more
# pick EEG channels
-picks = fiff.pick_types(raw.info, meg=False, eeg=True, stim=False,
+picks = fiff.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=True,
include=include, exclude=exclude)
# Read epochs
epochs = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
-epochs.reject(grad=4000e-13, mag=4e-12, eeg=40e-6, eog=150e-6)
+epochs.reject(eeg=40e-6, eog=150e-6)
evoked = epochs.average() # average epochs and get an Evoked dataset.
evoked.save('sample_audvis_eeg-ave.fif') # save evoked data to disk
diff --git a/examples/plot_minimum_norm_estimate.py b/examples/plot_minimum_norm_estimate.py
index 9ae2d4e..b6c88de 100644
--- a/examples/plot_minimum_norm_estimate.py
+++ b/examples/plot_minimum_norm_estimate.py
@@ -33,8 +33,7 @@ dSPM = True
# Load data
evoked = Evoked(fname_evoked, setno=setno, baseline=(None, 0))
forward = mne.read_forward_solution(fname_fwd)
-noise_cov = mne.Covariance()
-noise_cov.load(fname_cov)
+noise_cov = mne.Covariance(fname_cov)
# Compute inverse solution
stc, K, W = mne.minimum_norm(evoked, forward, noise_cov, orientation='loose',
diff --git a/examples/plot_read_epochs.py b/examples/plot_read_epochs.py
index d10cca6..6f30c51 100644
--- a/examples/plot_read_epochs.py
+++ b/examples/plot_read_epochs.py
@@ -33,7 +33,7 @@ events = mne.read_events(event_fname)
# Set up pick list: EEG + MEG - bad channels (modify to your needs)
exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more
-picks = fiff.pick_types(raw.info, meg=True, eeg=False, stim=True,
+picks = fiff.pick_types(raw.info, meg=True, eeg=False, stim=True, eog=True,
exclude=exclude)
# Read epochs
diff --git a/examples/plot_read_noise_covariance_matrix.py b/examples/plot_read_noise_covariance_matrix.py
index 8e9b0ab..9efd639 100644
--- a/examples/plot_read_noise_covariance_matrix.py
+++ b/examples/plot_read_noise_covariance_matrix.py
@@ -15,9 +15,7 @@ from mne.datasets import sample
data_path = sample.data_path('.')
fname = data_path + '/MEG/sample/sample_audvis-cov.fif'
-cov = mne.Covariance(kind='full')
-cov.load(fname)
-
+cov = mne.Covariance(fname)
print cov
###############################################################################
diff --git a/examples/plot_whitened_evoked_data.py b/examples/plot_whitened_evoked_data.py
index 01985ac..890df5c 100644
--- a/examples/plot_whitened_evoked_data.py
+++ b/examples/plot_whitened_evoked_data.py
@@ -35,15 +35,14 @@ events = mne.find_events(raw)
# pick EEG channels - bad channels (modify to your needs)
exclude = raw.info['bads'] + ['EEG 053'] # bads + 1 more
-picks = fiff.pick_types(raw.info, meg=False, eeg=True, stim=False,
+picks = fiff.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=True,
exclude=exclude)
epochs = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
-epochs.reject(grad=4000e-13, mag=4e-12, eeg=40e-6, eog=150e-6)
+epochs.reject(eeg=40e-6, eog=150e-6)
evoked = epochs.average() # average epochs and get an Evoked dataset.
-cov = mne.Covariance()
-cov.load(cov_fname)
+cov = mne.Covariance(cov_fname)
# Whiten data
W, ch_names = cov.whitener(evoked.info, pca=False) # get whitening matrix
diff --git a/examples/stats/plot_cluster_1samp_test_time_frequency.py b/examples/stats/plot_cluster_1samp_test_time_frequency.py
index a7fad3a..29312fd 100644
--- a/examples/stats/plot_cluster_1samp_test_time_frequency.py
+++ b/examples/stats/plot_cluster_1samp_test_time_frequency.py
@@ -46,22 +46,23 @@ include = []
exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more
# picks MEG gradiometers
-picks = fiff.pick_types(raw.info, meg='grad', eeg=False,
+picks = fiff.pick_types(raw.info, meg='grad', eeg=False, eog=True,
stim=False, include=include, exclude=exclude)
-picks = [picks[97]]
-ch_name = raw.info['ch_names'][picks[0]]
-
# Load condition 1
event_id = 1
epochs = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
-epochs.reject(grad=4000e-13, mag=4e-12, eeg=40e-6, eog=150e-6)
+epochs.reject(grad=4000e-13, eog=150e-6)
data = epochs.get_data() # as 3D matrix
data *= 1e13 # change unit to fT / cm
# Time vector
times = 1e3 * epochs.times # change unit to ms
+# Take only one channel
+ch_name = raw.info['ch_names'][97]
+data = data[:,97:98,:]
+
evoked_data = np.mean(data, 0)
# data -= evoked_data[None,:,:] # remove evoked component
diff --git a/examples/stats/plot_cluster_stats_evoked.py b/examples/stats/plot_cluster_stats_evoked.py
index 08fbe55..05bf10d 100644
--- a/examples/stats/plot_cluster_stats_evoked.py
+++ b/examples/stats/plot_cluster_stats_evoked.py
@@ -38,17 +38,17 @@ include = [channel]
###############################################################################
# Read epochs for the channel of interest
-picks = fiff.pick_types(raw.info, meg=False, include=include)
+picks = fiff.pick_types(raw.info, meg=False, eog=True, include=include)
event_id = 1
epochs1 = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
-epochs1.reject(grad=4000e-13, mag=4e-12, eeg=40e-6, eog=150e-6)
+epochs1.reject(grad=4000e-13, eog=150e-6)
condition1 = epochs1.get_data() # as 3D matrix
event_id = 2
epochs2 = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
-epochs2.reject(grad=4000e-13, mag=4e-12, eeg=40e-6, eog=150e-6)
+epochs2.reject(grad=4000e-13, eog=150e-6)
condition2 = epochs2.get_data() # as 3D matrix
condition1 = condition1[:,0,:] # take only one channel to get a 2D array
diff --git a/examples/stats/plot_cluster_stats_time_frequency.py b/examples/stats/plot_cluster_stats_time_frequency.py
index 55c0b06..ac3eb05 100644
--- a/examples/stats/plot_cluster_stats_time_frequency.py
+++ b/examples/stats/plot_cluster_stats_time_frequency.py
@@ -47,17 +47,16 @@ include = []
exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more
# picks MEG gradiometers
-picks = fiff.pick_types(raw.info, meg='grad', eeg=False,
+picks = fiff.pick_types(raw.info, meg='grad', eeg=False, eog=True,
stim=False, include=include, exclude=exclude)
-picks = [picks[97]]
ch_name = raw.info['ch_names'][picks[0]]
# Load condition 1
event_id = 1
epochs_condition_1 = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
-epochs_condition_1.reject(grad=4000e-13, mag=4e-12, eeg=40e-6, eog=150e-6)
+epochs_condition_1.reject(grad=4000e-13, eog=150e-6)
data_condition_1 = epochs_condition_1.get_data() # as 3D matrix
data_condition_1 *= 1e13 # change unit to fT / cm
@@ -65,10 +64,14 @@ data_condition_1 *= 1e13 # change unit to fT / cm
event_id = 2
epochs_condition_2 = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
-epochs_condition_2.reject(grad=4000e-13, mag=4e-12, eeg=40e-6, eog=150e-6)
+epochs_condition_2.reject(grad=4000e-13, eog=150e-6)
data_condition_2 = epochs_condition_2.get_data() # as 3D matrix
data_condition_2 *= 1e13 # change unit to fT / cm
+# Take only one channel
+data_condition_1 = data_condition_1[:,97:98,:]
+data_condition_2 = data_condition_2[:,97:98,:]
+
# Time vector
times = 1e3 * epochs_condition_1.times # change unit to ms
diff --git a/examples/stats/plot_sensor_permutation_test.py b/examples/stats/plot_sensor_permutation_test.py
index 17d3d1a..7dd8e9c 100644
--- a/examples/stats/plot_sensor_permutation_test.py
+++ b/examples/stats/plot_sensor_permutation_test.py
@@ -40,11 +40,11 @@ include = [] # or stim channel ['STI 014']
exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more
# pick MEG Magnetometers
-picks = fiff.pick_types(raw.info, meg='grad', eeg=False, stim=False,
+picks = fiff.pick_types(raw.info, meg='grad', eeg=False, stim=False, eog=True,
include=include, exclude=exclude)
epochs = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
-epochs.reject(grad=4000e-13, mag=4e-12, eog=150e-6)
+epochs.reject(grad=4000e-13, eog=150e-6)
data = epochs.get_data()
times = epochs.times
diff --git a/examples/time_frequency/plot_time_frequency.py b/examples/time_frequency/plot_time_frequency.py
index 7933d40..7133d52 100644
--- a/examples/time_frequency/plot_time_frequency.py
+++ b/examples/time_frequency/plot_time_frequency.py
@@ -38,19 +38,22 @@ include = []
exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more
# picks MEG gradiometers
-picks = fiff.pick_types(raw.info, meg='grad', eeg=False,
+picks = fiff.pick_types(raw.info, meg='grad', eeg=False, eog=True,
stim=False, include=include, exclude=exclude)
-picks = [picks[97]]
epochs = mne.Epochs(raw, events, event_id,
tmin, tmax, picks=picks, baseline=(None, 0))
-epochs.reject(grad=4000e-13, mag=4e-12, eeg=40e-6, eog=150e-6)
+epochs.reject(grad=4000e-13, eog=150e-6)
data = epochs.get_data() # as 3D matrix
evoked = epochs.average() # compute evoked fields
times = 1e3 * epochs.times # change unit to ms
evoked_data = evoked.data * 1e13 # change unit to fT / cm
+# Take only one channel
+data = data[:,97:98,:]
+evoked_data = evoked_data[97:98,:]
+
frequencies = np.arange(7, 30, 3) # define frequencies of interest
Fs = raw.info['sfreq'] # sampling in Hz
power, phase_lock = induced_power(data, Fs=Fs, frequencies=frequencies,
diff --git a/mne/__init__.py b/mne/__init__.py
index d78e642..c471729 100644
--- a/mne/__init__.py
+++ b/mne/__init__.py
@@ -1,6 +1,7 @@
__version__ = '0.1.git'
-from .cov import read_cov, write_cov, write_cov_file, Covariance
+from .cov import read_cov, write_cov, write_cov_file, Covariance, \
+ noise_covariance_segment, noise_covariance
from .event import read_events, write_events, find_events
from .forward import read_forward_solution
from .stc import read_stc, write_stc
@@ -9,5 +10,5 @@ from .source_space import read_source_spaces
from .inverse import read_inverse_operator, apply_inverse, minimum_norm
from .epochs import Epochs
from .label import label_time_courses, read_label
-from .misc import parse_config
+from .misc import parse_config, read_reject_parameters
import fiff
diff --git a/mne/cov.py b/mne/cov.py
index f0d24af..2873a5c 100644
--- a/mne/cov.py
+++ b/mne/cov.py
@@ -4,6 +4,7 @@
# License: BSD (3-clause)
import os
+from math import floor, ceil
import numpy as np
from scipy import linalg
@@ -18,7 +19,7 @@ from .fiff.write import start_block, end_block, write_int, write_name_list, \
from .fiff.proj import write_proj, make_projector
from .fiff import fiff_open
from .fiff.pick import pick_types
-
+from .epochs import Epochs
def rank(A, tol=1e-8):
s = linalg.svd(A, compute_uv=0)
@@ -51,11 +52,11 @@ class Covariance(object):
_kind_to_id = dict(full=1, sparse=2, diagonal=3) # XXX : check
_id_to_kind = {1: 'full', 2: 'sparse', 3: 'diagonal'} # XXX : check
- def __init__(self, kind='full'):
+ def __init__(self, fname, kind='full'):
self.kind = kind
- def load(self, fname):
- """load covariance matrix from FIF file"""
+ if fname is None:
+ return
if self.kind in Covariance._kind_to_id:
cov_kind = Covariance._kind_to_id[self.kind]
@@ -194,43 +195,14 @@ class Covariance(object):
return W, names_meg + names_eeg
- def estimate_from_raw(self, raw, picks=None, quantum_sec=10):
- """Estimate noise covariance matrix from a raw FIF file
- """
- # Set up the reading parameters
- start = raw.first_samp
- stop = raw.last_samp + 1
- quantum = int(quantum_sec * raw.info['sfreq'])
-
- cov = 0
- n_samples = 0
-
- # Read data
- for first in range(start, stop, quantum):
- last = first + quantum
- if last >= stop:
- last = stop
-
- data, times = raw[picks, first:last]
-
- if self.kind is 'full':
- cov += np.dot(data, data.T)
- elif self.kind is 'diagonal':
- cov += np.diag(np.sum(data ** 2, axis=1))
- else:
- raise ValueError("Unsupported covariance kind")
-
- n_samples += data.shape[1]
-
- self.data = cov / n_samples # XXX : check
- print '[done]'
-
def __repr__(self):
s = "kind : %s" % self.kind
s += ", size : %s x %s" % self.data.shape
s += ", data : %s" % self.data
return "Covariance (%s)" % s
+###############################################################################
+# IO
def read_cov(fid, node, cov_kind):
"""Read a noise covariance matrix
@@ -259,7 +231,8 @@ def read_cov(fid, node, cov_kind):
# Is any of the covariance matrices a noise covariance
for p in range(len(covs)):
tag = find_tag(fid, covs[p], FIFF.FIFF_MNE_COV_KIND)
- if tag is not None and tag.data == cov_kind:
+
+ if tag is not None and int(tag.data) == cov_kind:
this = covs[p]
# Find all the necessary data
@@ -340,6 +313,124 @@ def read_cov(fid, node, cov_kind):
return None
###############################################################################
+# Estimate from data
+
+def _estimate_noise_covariance_from_epochs(epochs, bmin, bmax,
+ keep_sample_mean):
+ """Estimate noise covariance matrix from epochs
+
+ XXX : doc missing
+ """
+ picks_no_eog = pick_types(epochs.info, meg=True, eeg=True, eog=False)
+ n_channels = len(picks_no_eog)
+ ch_names = [epochs.ch_names[k] for k in picks_no_eog]
+ data = np.zeros((n_channels, n_channels))
+ n_samples = 0
+ if bmin is None:
+ bmin = epochs.times[0]
+ if bmax is None:
+ bmax = epochs.times[-1]
+ bmask = (epochs.times >= bmin) & (epochs.times <= bmax)
+ for e in epochs:
+ e = e[picks_no_eog]
+ mu = e[:,bmask].mean(axis=1)
+ e -= mu[:,None]
+ if not keep_sample_mean:
+ e -= np.mean(e, axis=0)
+ data += np.dot(e, e.T)
+ n_samples += e.shape[1]
+ print "Number of samples used : %d" % n_samples
+ print '[done]'
+ return data, n_samples, ch_names
+
+
+def noise_covariance_segment(raw, reject_params=None, keep_sample_mean=True,
+ tmin=None, tmax=None, tstep=0.2):
+ """Estimate noise covariance matrix from a continuous segment of raw data
+
+ XXX : doc missing
+
+ Parameters
+ ----------
+ Returns
+ -------
+ """
+ sfreq = raw.info['sfreq']
+
+ # Convert to samples
+ start = raw.first_samp if tmin is None else int(floor(tmin * sfreq))
+ stop = raw.last_samp if tmax is None else int(ceil(tmax * sfreq))
+ step = int(ceil(tstep * raw.info['sfreq']))
+
+ picks = pick_types(raw.info, meg=True, eeg=True, eog=True)
+ picks_data = pick_types(raw.info, meg=True, eeg=True, eog=False)
+ idx = [list(picks).index(k) for k in picks_data]
+
+ data = 0
+ n_samples = 0
+ mu = 0
+
+ # Read data in chuncks
+ for first in range(start, stop, step):
+ last = first + step
+ if last >= stop:
+ last = stop
+ raw_segment, times = raw[picks, first:last]
+ # XXX : check for artefacts
+ mu += raw_segment[idx].sum(axis=1)
+ data += np.dot(raw_segment[idx], raw_segment[idx].T)
+ n_samples += raw_segment.shape[1]
+
+ mu /= n_samples
+ data -= n_samples * mu[:,None] * mu[None,:]
+ data /= (n_samples - 1.0)
+ print "Number of samples used : %d" % n_samples
+ print '[done]'
+
+ cov = Covariance(None)
+ cov.data = data
+ cov.ch_names = [raw.info['ch_names'][k] for k in picks_data]
+ return cov
+
+
+def noise_covariance(raw, events, event_ids, tmin, tmax,
+ bmin=None, bmax=None, reject_params=None,
+ keep_sample_mean=True):
+ """Estimate noise covariance matrix from raw file
+
+ XXX : doc missing
+
+ Parameters
+ ----------
+ Returns
+ -------
+ """
+ # Pick all channels
+ picks = pick_types(raw.info, meg=True, eeg=True, eog=True)
+ if isinstance(event_ids, int):
+ event_ids = list(event_ids)
+ data = 0.0
+ n_samples = 0
+
+ for event_id in event_ids:
+ print "Processing event : %d" % event_id
+ epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+ baseline=(None, 0))
+ d, n, ch_names = _estimate_noise_covariance_from_epochs(epochs,
+ bmin=bmin, bmax=bmax, keep_sample_mean=keep_sample_mean)
+
+ # XXX : check artefacts
+ data += d
+ n_samples += n
+
+ data /= n_samples - 1
+ cov = Covariance(None)
+ cov.data = data
+ cov.ch_names = ch_names
+ return cov
+
+
+###############################################################################
# Writing
def write_cov(fid, cov):
diff --git a/mne/epochs.py b/mne/epochs.py
index ef167ee..ea65192 100644
--- a/mne/epochs.py
+++ b/mne/epochs.py
@@ -81,9 +81,6 @@ class Epochs(object):
self.baseline = baseline
self.preload = preload
- if len(picks) == 0:
- raise ValueError, "Picks cannot be empty."
-
# Handle measurement info
self.info = copy.copy(raw.info)
if picks is not None:
@@ -97,6 +94,9 @@ class Epochs(object):
else:
self.ch_names = [raw.info['ch_names'][k] for k in picks]
+ if len(picks) == 0:
+ raise ValueError, "Picks cannot be empty."
+
# Set up projection
if raw.info['projs'] is None:
print 'No projector specified for these data'
@@ -144,7 +144,7 @@ class Epochs(object):
# Handle times
sfreq = raw.info['sfreq']
- self.times = np.arange(int(tmin*sfreq), int(tmax*sfreq),
+ self.times = np.arange(int(round(tmin*sfreq)), int(round(tmax*sfreq)),
dtype=np.float) / sfreq
if self.preload:
@@ -172,7 +172,7 @@ class Epochs(object):
event_samp = self.events[idx, 0]
# Read a data segment
- start = int(event_samp + self.tmin*sfreq)
+ start = int(round(event_samp + self.tmin*sfreq))
stop = start + len(self.times)
epoch, _ = self.raw[self.picks, start:stop]
diff --git a/mne/fiff/tests/data/test-cov.fif b/mne/fiff/tests/data/test-cov.fif
index 1a2da11..574ac2a 100644
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.cov b/mne/fiff/tests/data/test.cov
index 1f9fae2..a9d7782 100644
--- a/mne/fiff/tests/data/test.cov
+++ b/mne/fiff/tests/data/test.cov
@@ -5,15 +5,15 @@ cov {
#
# Output files
#
- outfile sample_audvis-cov.fif
- logfile sample_audvis-cov.log
+ outfile test-cov.fif
+ logfile test-cov.log
#
# Rejection values
#
- gradReject 4000e-13
- magReject 4e-12
- eegReject 40e-6
- eogReject 150e-6
+ gradReject 10000e-13
+ magReject 4e-12
+ eegReject 80e-6
+ eogReject 150e-6
#
# What to include in the covariance matrix?
#
diff --git a/mne/fiff/tests/data/test_empty_room.cov b/mne/fiff/tests/data/test_empty_room.cov
new file mode 100644
index 0000000..c038c76
--- /dev/null
+++ b/mne/fiff/tests/data/test_empty_room.cov
@@ -0,0 +1,44 @@
+cov {
+# name "Empty Room"
+#
+# Output files
+# The log file is useful for debugging and
+# selection of interesting events using 'eventfile'
+#
+ outfile test_erm-cov.fif
+ logfile test_erm-cov.log
+#
+# Rejection limits
+#
+# stimIgnore is optional to omit a stimulus artefact from
+# the rejection
+#
+# fixSkew
+# logfile erm-ave.log
+ # gradReject 10000e-13
+ # magReject 3e-12
+ # magFlat 1e-14
+ # gradflat 1000e-15
+
+# Additional rejection parameters
+#
+# eegReject 20e-6
+# ecgReject 10e-3
+#
+# The first definition follows
+#
+ def {
+#
+# The name of the category (condition) is irrelevant
+# but useful as a comment
+#
+# 'event' can be left out to compute covariance matrix
+# from continuous data
+#
+# 'ignore' is a mask to apply to the trigger line
+# before searching for 'event' (default = 0)
+#
+ tmin 0
+ tmax 99999
+ }
+}
\ No newline at end of file
diff --git a/mne/fiff/tests/data/test_erm-cov.fif b/mne/fiff/tests/data/test_erm-cov.fif
new file mode 100644
index 0000000..cd637f3
Binary files /dev/null and b/mne/fiff/tests/data/test_erm-cov.fif differ
diff --git a/mne/misc.py b/mne/misc.py
index 9c25a5b..271776c 100644
--- a/mne/misc.py
+++ b/mne/misc.py
@@ -19,26 +19,19 @@ def parse_config(fname):
tmin, tmax, name, grad_reject, mag_reject,
eeg_reject, eog_reject
"""
+ reject_params = read_reject_parameters(fname)
+
try:
with open(fname, 'r') as f:
- ave_lines = f.readlines()
+ lines = f.readlines()
except:
print("Error while reading %s" % fname)
- reject_names = ['gradReject', 'magReject', 'eegReject', 'eogReject']
- reject_pynames = ['grad_reject', 'mag_reject', 'eeg_reject', 'eog_reject']
- reject_params = dict()
- for line in ave_lines:
- words = line.split()
- if words[0] in reject_names:
- reject_params[reject_pynames[reject_names.index(words[0])]] = \
- float(words[1])
-
- cat_ind = [i for i, x in enumerate(ave_lines) if "category {" in x]
+ cat_ind = [i for i, x in enumerate(lines) if "category {" in x]
event_dict = dict()
for ind in cat_ind:
for k in range(ind+1, ind+7):
- words = ave_lines[k].split()
+ words = lines[k].split()
if len(words) >= 2:
key = words[0]
if key == 'event':
@@ -48,7 +41,7 @@ def parse_config(fname):
raise ValueError('Could not find event id.')
event_dict[event] = dict(**reject_params)
for k in range(ind+1, ind+7):
- words = ave_lines[k].split()
+ words = lines[k].split()
if len(words) >= 2:
key = words[0]
if key == 'name':
@@ -62,3 +55,23 @@ def parse_config(fname):
event_dict[event][key] = float(words[1])
return event_dict
+def read_reject_parameters(fname):
+ """Read rejection parameters from .cov or .ave config file"""
+
+ try:
+ with open(fname, 'r') as f:
+ lines = f.readlines()
+ except:
+ print("Error while reading %s" % fname)
+
+ reject_names = ['gradReject', 'magReject', 'eegReject', 'eogReject']
+ reject_pynames = ['grad', 'mag', 'eeg', 'eog']
+ reject_params = dict()
+ for line in lines:
+ words = line.split()
+ print words
+ if words[0] in reject_names:
+ reject_params[reject_pynames[reject_names.index(words[0])]] = \
+ float(words[1])
+
+ return reject_params
\ No newline at end of file
diff --git a/mne/tests/test_cov.py b/mne/tests/test_cov.py
index 08cf415..6ba2e31 100644
--- a/mne/tests/test_cov.py
+++ b/mne/tests/test_cov.py
@@ -1,21 +1,24 @@
import os.path as op
from numpy.testing import assert_array_almost_equal
+from scipy import linalg
import mne
-from ..fiff import fiff_open, read_evoked, pick_types
+from ..fiff import fiff_open, read_evoked, Raw
from ..datasets import sample
-fname = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data',
+cov_fname = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data',
'test-cov.fif')
raw_fname = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data',
'test_raw.fif')
+erm_cov_fname = op.join('mne', 'fiff', 'tests', 'data',
+ 'test_erm-cov.fif')
def test_io_cov():
"""Test IO for noise covariance matrices
"""
- fid, tree, _ = fiff_open(fname)
+ fid, tree, _ = fiff_open(cov_fname)
cov_type = 1
cov = mne.read_cov(fid, tree, cov_type)
fid.close()
@@ -29,24 +32,30 @@ def test_io_cov():
assert_array_almost_equal(cov['data'], cov2['data'])
-def test_cov_estimation():
- """Test estimation of noise covariance from raw data
+def test_cov_estimation_on_raw_segment():
+ """Estimate raw on continuous recordings (typically empty room)
"""
- raw = mne.fiff.Raw(raw_fname)
- # Set up pick list: MEG + STI 014 - bad channels
- want_meg = True
- want_eeg = False
- want_stim = False
+ raw = Raw(raw_fname)
+ cov = mne.noise_covariance_segment(raw)
+ cov_mne = mne.Covariance(erm_cov_fname)
+ assert cov_mne.ch_names == cov.ch_names
+ assert (linalg.norm(cov.data - cov_mne.data, ord='fro')
+ / linalg.norm(cov.data, ord='fro')) < 1e-6
- picks = pick_types(raw.info, meg=want_meg, eeg=want_eeg,
- stim=want_stim, exclude=raw.info['bads'])
- full_cov = mne.Covariance(kind='full')
- full_cov.estimate_from_raw(raw, picks=picks)
-
- diagonal_cov = mne.Covariance(kind='diagonal')
- diagonal_cov.estimate_from_raw(raw, picks=picks)
- # XXX : test something
+def test_cov_estimation_with_triggers():
+ """Estimate raw with triggers
+ """
+ raw = Raw(raw_fname)
+ events = mne.find_events(raw)
+ event_ids = [1, 2, 3, 4]
+ cov = mne.noise_covariance(raw, events, event_ids,
+ tmin=-0.2, tmax=0, keep_sample_mean=True)
+ cov_mne = mne.Covariance(cov_fname)
+ assert cov_mne.ch_names == cov.ch_names
+ # XXX : check something
+ # assert (linalg.norm(cov.data - cov_mne.data, ord='fro')
+ # / linalg.norm(cov.data, ord='fro')) < 0.1 # XXX : fix
def test_whitening_cov():
@@ -61,8 +70,7 @@ def test_whitening_cov():
# Reading
evoked = read_evoked(ave_fname, setno=0, baseline=(None, 0))
- cov = mne.Covariance()
- cov.load(cov_fname)
+ cov = mne.Covariance(cov_fname)
cov.whitener(evoked.info)
# XXX : test something
diff --git a/mne/tests/test_epochs.py b/mne/tests/test_epochs.py
index ee0c851..df1170f 100644
--- a/mne/tests/test_epochs.py
+++ b/mne/tests/test_epochs.py
@@ -14,6 +14,8 @@ event_name = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data',
def test_read_epochs():
+ """Reading epochs from raw files
+ """
event_id = 1
tmin = -0.2
tmax = 0.5
@@ -29,6 +31,8 @@ def test_read_epochs():
def test_reject_epochs():
+ """Test of epochs rejection
+ """
event_id = 1
tmin = -0.2
tmax = 0.5
@@ -44,5 +48,8 @@ def test_reject_epochs():
n_epochs = len(epochs)
epochs.reject(grad=1000e-12, mag=4e-12, eeg=80e-6, eog=150e-6)
n_clean_epochs = len(epochs)
+ # Should match
+ # mne_process_raw --raw test_raw.fif --projoff \
+ # --saveavetag -ave --ave test.ave --filteroff
assert n_epochs > n_clean_epochs
assert n_clean_epochs == 3
--
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