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