# Authors: Alexandre Gramfort # Matti Hamalainen # Denis Engemann # Andrew Dykstra # # License: BSD (3-clause) from copy import deepcopy import numpy as np import warnings from .baseline import rescale from .channels import ContainsMixin, PickDropChannelsMixin from .filter import resample, detrend from .fixes import in1d from .utils import (_check_pandas_installed, check_fname, logger, verbose, deprecated, object_hash) from .viz import plot_evoked, plot_evoked_topomap, _mutable_defaults from .viz import plot_evoked_field from .viz import plot_evoked_image from .externals.six import string_types from .io.constants import FIFF from .io.open import fiff_open from .io.tag import read_tag from .io.tree import dir_tree_find from .io.pick import channel_type, pick_types from .io.meas_info import read_meas_info, write_meas_info from .io.proj import ProjMixin from .io.write import (start_file, start_block, end_file, end_block, write_int, write_string, write_float_matrix, write_id) aspect_dict = {'average': FIFF.FIFFV_ASPECT_AVERAGE, 'standard_error': FIFF.FIFFV_ASPECT_STD_ERR} aspect_rev = {str(FIFF.FIFFV_ASPECT_AVERAGE): 'average', str(FIFF.FIFFV_ASPECT_STD_ERR): 'standard_error'} class Evoked(ProjMixin, ContainsMixin, PickDropChannelsMixin): """Evoked data Parameters ---------- fname : string Name of evoked/average FIF file to load. If None no data is loaded. condition : int, or str Dataset ID number (int) or comment/name (str). Optional if there is only one data set in file. baseline : tuple or list of length 2, or None The time interval to apply rescaling / baseline correction. If None do not apply it. If baseline is (a, b) the interval is between "a (s)" and "b (s)". If a is None the beginning of the data is used and if b is None then b is set to the end of the interval. If baseline is equal ot (None, None) all the time interval is used. If None, no correction is applied. proj : bool, optional Apply SSP projection vectors kind : str Either 'average' or 'standard_error'. The type of data to read. Only used if 'condition' is a str. verbose : bool, str, int, or None If not None, override default verbose level (see mne.verbose). Attributes ---------- info : dict Measurement info. ch_names : list of string List of channels' names. nave : int Number of averaged epochs. kind : str Type of data, either average or standard_error. first : int First time sample. last : int Last time sample. comment : string Comment on dataset. Can be the condition. times : array Array of time instants in seconds. data : array of shape (n_channels, n_times) Evoked response. verbose : bool, str, int, or None. See above. """ @verbose def __init__(self, fname, condition=None, baseline=None, proj=True, kind='average', verbose=None, setno=None): if fname is None: raise ValueError('No evoked filename specified') if condition is None and setno is not None: condition = setno msg = ("'setno' will be deprecated in 0.9. Use 'condition' " "instead.") warnings.warn(msg, DeprecationWarning) self.verbose = verbose logger.info('Reading %s ...' % fname) fid, tree, _ = fiff_open(fname) if not isinstance(proj, bool): raise ValueError(r"'proj' must be 'True' or 'False'") # Read the measurement info info, meas = read_meas_info(fid, tree) info['filename'] = fname # Locate the data of interest processed = dir_tree_find(meas, FIFF.FIFFB_PROCESSED_DATA) if len(processed) == 0: fid.close() raise ValueError('Could not find processed data') evoked_node = dir_tree_find(meas, FIFF.FIFFB_EVOKED) if len(evoked_node) == 0: fid.close() raise ValueError('Could not find evoked data') # find string-based entry if isinstance(condition, string_types): if not kind in aspect_dict.keys(): fid.close() raise ValueError('kind must be "average" or ' '"standard_error"') comments, aspect_kinds, t = _get_entries(fid, evoked_node) goods = np.logical_and(in1d(comments, [condition]), in1d(aspect_kinds, [aspect_dict[kind]])) found_cond = np.where(goods)[0] if len(found_cond) != 1: fid.close() raise ValueError('condition "%s" (%s) not found, out of found ' 'datasets:\n %s' % (condition, kind, t)) condition = found_cond[0] if condition >= len(evoked_node) or condition < 0: fid.close() raise ValueError('Data set selector out of range') my_evoked = evoked_node[condition] # Identify the aspects aspects = dir_tree_find(my_evoked, FIFF.FIFFB_ASPECT) if len(aspects) > 1: logger.info('Multiple aspects found. Taking first one.') my_aspect = aspects[0] # Now find the data in the evoked block nchan = 0 sfreq = -1 chs = [] comment = None for k in range(my_evoked['nent']): my_kind = my_evoked['directory'][k].kind pos = my_evoked['directory'][k].pos if my_kind == FIFF.FIFF_COMMENT: tag = read_tag(fid, pos) comment = tag.data elif my_kind == FIFF.FIFF_FIRST_SAMPLE: tag = read_tag(fid, pos) first = int(tag.data) elif my_kind == FIFF.FIFF_LAST_SAMPLE: tag = read_tag(fid, pos) last = int(tag.data) elif my_kind == FIFF.FIFF_NCHAN: tag = read_tag(fid, pos) nchan = int(tag.data) elif my_kind == FIFF.FIFF_SFREQ: tag = read_tag(fid, pos) sfreq = float(tag.data) elif my_kind == FIFF.FIFF_CH_INFO: tag = read_tag(fid, pos) chs.append(tag.data) if comment is None: comment = 'No comment' # Local channel information? if nchan > 0: if chs is None: fid.close() raise ValueError('Local channel information was not found ' 'when it was expected.') if len(chs) != nchan: fid.close() raise ValueError('Number of channels and number of ' 'channel definitions are different') info['chs'] = chs info['nchan'] = nchan logger.info(' Found channel information in evoked data. ' 'nchan = %d' % nchan) if sfreq > 0: info['sfreq'] = sfreq nsamp = last - first + 1 logger.info(' Found the data of interest:') logger.info(' t = %10.2f ... %10.2f ms (%s)' % (1000 * first / info['sfreq'], 1000 * last / info['sfreq'], comment)) if info['comps'] is not None: logger.info(' %d CTF compensation matrices available' % len(info['comps'])) # Read the data in the aspect block nave = 1 epoch = [] for k in range(my_aspect['nent']): kind = my_aspect['directory'][k].kind pos = my_aspect['directory'][k].pos if kind == FIFF.FIFF_COMMENT: tag = read_tag(fid, pos) comment = tag.data elif kind == FIFF.FIFF_ASPECT_KIND: tag = read_tag(fid, pos) aspect_kind = int(tag.data) elif kind == FIFF.FIFF_NAVE: tag = read_tag(fid, pos) nave = int(tag.data) elif kind == FIFF.FIFF_EPOCH: tag = read_tag(fid, pos) epoch.append(tag) logger.info(' nave = %d - aspect type = %d' % (nave, aspect_kind)) nepoch = len(epoch) if nepoch != 1 and nepoch != info['nchan']: fid.close() raise ValueError('Number of epoch tags is unreasonable ' '(nepoch = %d nchan = %d)' % (nepoch, info['nchan'])) if nepoch == 1: # Only one epoch all_data = epoch[0].data.astype(np.float) # May need a transpose if the number of channels is one if all_data.shape[1] == 1 and info['nchan'] == 1: all_data = all_data.T.astype(np.float) else: # Put the old style epochs together all_data = np.concatenate([e.data[None, :] for e in epoch], axis=0).astype(np.float) if all_data.shape[1] != nsamp: fid.close() raise ValueError('Incorrect number of samples (%d instead of %d)' % (all_data.shape[1], nsamp)) # Calibrate cals = np.array([info['chs'][k]['cal'] * info['chs'][k].get('scale', 1.0) for k in range(info['nchan'])]) all_data *= cals[:, np.newaxis] times = np.arange(first, last + 1, dtype=np.float) / info['sfreq'] self.info = info # Put the rest together all together self.nave = nave self._aspect_kind = aspect_kind self.kind = aspect_rev.get(str(self._aspect_kind), 'Unknown') self.first = first self.last = last self.comment = comment self.times = times self.data = all_data # bind info, proj, data to self so apply_proj can be used self.data = all_data self.proj = False if proj: self.apply_proj() # Run baseline correction self.data = rescale(self.data, times, baseline, 'mean', copy=False) fid.close() def save(self, fname): """Save dataset to file. Parameters ---------- fname : string Name of the file where to save the data. """ write_evokeds(fname, self) def __repr__(self): s = "comment : '%s'" % self.comment s += ", time : [%f, %f]" % (self.times[0], self.times[-1]) s += ", n_epochs : %d" % self.nave s += ", n_channels x n_times : %s x %s" % self.data.shape return "" % s @property def ch_names(self): return self.info['ch_names'] def crop(self, tmin=None, tmax=None): """Crop data to a given time interval """ times = self.times mask = np.ones(len(times), dtype=np.bool) if tmin is not None: mask = mask & (times >= tmin) if tmax is not None: mask = mask & (times <= tmax) self.times = times[mask] self.first = int(self.times[0] * self.info['sfreq']) self.last = len(self.times) + self.first - 1 self.data = self.data[:, mask] def shift_time(self, tshift, relative=True): """Shift time scale in evoked data Parameters ---------- tshift : float The amount of time shift to be applied if relative is True else the first time point. When relative is True, positive value of tshift moves the data forward while negative tshift moves it backward. relative : bool If true, move the time backwards or forwards by specified amount. Else, set the starting time point to the value of tshift. Notes ----- Maximum accuracy of time shift is 1 / evoked.info['sfreq'] """ times = self.times sfreq = self.info['sfreq'] offset = self.first if relative else 0 self.first = int(tshift * sfreq) + offset self.last = self.first + len(times) - 1 self.times = np.arange(self.first, self.last + 1, dtype=np.float) / sfreq def plot(self, picks=None, exclude='bads', unit=True, show=True, ylim=None, proj=False, xlim='tight', hline=None, units=None, scalings=None, titles=None, axes=None): """Plot evoked data as butterfly plots Note: If bad channels are not excluded they are shown in red. Parameters ---------- picks : array-like of int | None The indices of channels to plot. If None show all. exclude : list of str | 'bads' Channels names to exclude from being shown. If 'bads', the bad channels are excluded. unit : bool Scale plot with channel (SI) unit. show : bool Call pyplot.show() at the end or not. ylim : dict ylim for plots. e.g. ylim = dict(eeg=[-200e-6, 200e6]) Valid keys are eeg, mag, grad xlim : 'tight' | tuple | None xlim for plots. proj : bool | 'interactive' If true SSP projections are applied before display. If 'interactive', a check box for reversible selection of SSP projection vectors will be shown. hline : list of floats | None The values at which show an horizontal line. units : dict | None The units of the channel types used for axes lables. If None, defaults to `dict(eeg='uV', grad='fT/cm', mag='fT')`. scalings : dict | None The scalings of the channel types to be applied for plotting. If None, defaults to `dict(eeg=1e6, grad=1e13, mag=1e15)`. titles : dict | None The titles associated with the channels. If None, defaults to `dict(eeg='EEG', grad='Gradiometers', mag='Magnetometers')`. axes : instance of Axes | list | None The axes to plot to. If list, the list must be a list of Axes of the same length as the number of channel types. If instance of Axes, there must be only one channel type plotted. """ return plot_evoked(self, picks=picks, exclude=exclude, unit=unit, show=show, ylim=ylim, proj=proj, xlim=xlim, hline=hline, units=units, scalings=scalings, titles=titles, axes=axes) def plot_image(self, picks=None, exclude='bads', unit=True, show=True, clim=None, proj=False, xlim='tight', units=None, scalings=None, titles=None, axes=None, cmap='RdBu_r'): """Plot evoked data as images Parameters ---------- picks : array-like of int | None The indices of channels to plot. If None show all. exclude : list of str | 'bads' Channels names to exclude from being shown. If 'bads', the bad channels are excluded. unit : bool Scale plot with channel (SI) unit. show : bool Call pyplot.show() at the end or not. clim : dict clim for images. e.g. clim = dict(eeg=[-200e-6, 200e6]) Valid keys are eeg, mag, grad xlim : 'tight' | tuple | None xlim for plots. proj : bool | 'interactive' If true SSP projections are applied before display. If 'interactive', a check box for reversible selection of SSP projection vectors will be shown. units : dict | None The units of the channel types used for axes lables. If None, defaults to `dict(eeg='uV', grad='fT/cm', mag='fT')`. scalings : dict | None The scalings of the channel types to be applied for plotting. If None, defaults to `dict(eeg=1e6, grad=1e13, mag=1e15)`. titles : dict | None The titles associated with the channels. If None, defaults to `dict(eeg='EEG', grad='Gradiometers', mag='Magnetometers')`. axes : instance of Axes | list | None The axes to plot to. If list, the list must be a list of Axes of the same length as the number of channel types. If instance of Axes, there must be only one channel type plotted. cmap : matplotlib colormap Colormap. """ return plot_evoked_image(self, picks=picks, exclude=exclude, unit=unit, show=show, clim=clim, proj=proj, xlim=xlim, units=units, scalings=scalings, titles=titles, axes=axes, cmap=cmap) def plot_topomap(self, times=None, ch_type='mag', layout=None, vmin=None, vmax=None, cmap='RdBu_r', sensors='k,', colorbar=True, scale=None, scale_time=1e3, unit=None, res=64, size=1, format="%3.1f", time_format='%01d ms', proj=False, show=True, show_names=False, title=None, mask=None, mask_params=None, outlines='head', contours=6, image_interp='bilinear'): """Plot topographic maps of specific time points Parameters ---------- times : float | array of floats | None. The time point(s) to plot. If None, 10 topographies will be shown will a regular time spacing between the first and last time instant. ch_type : 'mag' | 'grad' | 'planar1' | 'planar2' | 'eeg' The channel type to plot. For 'grad', the gradiometers are collec- ted in pairs and the RMS for each pair is plotted. layout : None | Layout Layout instance specifying sensor positions (does not need to be specified for Neuromag data). If possible, the correct layout file is inferred from the data; if no appropriate layout file was found, the layout is automatically generated from the sensor locations. vmin : float | callable The value specfying the lower bound of the color range. If None, and vmax is None, -vmax is used. Else np.min(data). If callable, the output equals vmin(data). vmax : float | callable The value specfying the upper bound of the color range. If None, the maximum absolute value is used. If vmin is None, but vmax is not, defaults to np.min(data). If callable, the output equals vmax(data). cmap : matplotlib colormap Colormap. Defaults to 'RdBu_r' sensors : bool | str Add markers for sensor locations to the plot. Accepts matplotlib plot format string (e.g., 'r+' for red plusses). colorbar : bool Plot a colorbar. scale : float | None Scale the data for plotting. If None, defaults to 1e6 for eeg, 1e13 for grad and 1e15 for mag. scale_time : float | None Scale the time labels. Defaults to 1e3 (ms). units : str | None The units of the channel types used for colorbar lables. If scale == None the unit is automatically determined. res : int The resolution of the topomap image (n pixels along each side). size : scalar Side length of the topomaps in inches (only applies when plotting multiple topomaps at a time). format : str String format for colorbar values. time_format : str String format for topomap values. Defaults to "%01d ms" proj : bool | 'interactive' If true SSP projections are applied before display. If 'interactive', a check box for reversible selection of SSP projection vectors will be shown. show : bool Call pyplot.show() at the end. show_names : bool | callable If True, show channel names on top of the map. If a callable is passed, channel names will be formatted using the callable; e.g., to delete the prefix 'MEG ' from all channel names, pass the function lambda x: x.replace('MEG ', ''). If `mask` is not None, only significant sensors will be shown. title : str | None Title. If None (default), no title is displayed. mask : ndarray of bool, shape (n_channels, n_times) | None The channels to be marked as significant at a given time point. Indicies set to `True` will be considered. Defaults to None. mask_params : dict | None Additional plotting parameters for plotting significant sensors. Default (None) equals: dict(marker='o', markerfacecolor='w', markeredgecolor='k', linewidth=0, markersize=4) outlines : 'head' | dict | None The outlines to be drawn. If 'head', a head scheme will be drawn. If dict, each key refers to a tuple of x and y positions. The values in 'mask_pos' will serve as image mask. If None, nothing will be drawn. Defaults to 'head'. image_interp : str The image interpolation to be used. All matplotlib options are accepted. """ return plot_evoked_topomap(self, times=times, ch_type=ch_type, layout=layout, vmin=vmin, vmax=vmax, cmap=cmap, sensors=sensors, colorbar=colorbar, scale=scale, scale_time=scale_time, unit=unit, res=res, proj=proj, size=size, format=format, time_format=time_format, show=show, show_names=show_names, title=title, mask=mask, mask_params=mask_params, outlines=outlines, contours=contours, image_interp=image_interp) def plot_field(self, surf_maps, time=None, time_label='t = %0.0f ms', n_jobs=1): """Plot MEG/EEG fields on head surface and helmet in 3D Parameters ---------- surf_maps : list The surface mapping information obtained with make_field_map. time : float | None The time point at which the field map shall be displayed. If None, the average peak latency (across sensor types) is used. time_label : str How to print info about the time instant visualized. n_jobs : int Number of jobs to run in parallel. Returns ------- fig : instance of mlab.Figure The mayavi figure. """ return plot_evoked_field(self, surf_maps, time=time, time_label=time_label, n_jobs=n_jobs) def to_nitime(self, picks=None): """Export Evoked object to NiTime Parameters ---------- picks : array-like of int | None Indices of channels to apply. If None, all channels will be exported. Returns ------- evoked_ts : instance of nitime.TimeSeries The TimeSeries instance """ try: from nitime import TimeSeries # to avoid strong dependency except ImportError: raise Exception('the nitime package is missing') evoked_ts = TimeSeries(self.data if picks is None else self.data[picks], sampling_rate=self.info['sfreq']) return evoked_ts def as_data_frame(self, picks=None, scale_time=1e3, scalings=None, use_time_index=True, copy=True): """Get the Evoked object as a Pandas DataFrame Export data in tabular structure: each row corresponds to a time point, and each column to a channel. Parameters ---------- picks : array-like of int | None If None all channels are kept, otherwise the channels indices in picks are kept. scale_time : float Scaling to be applied to time units. scalings : dict | None Scaling to be applied to the channels picked. If None, defaults to ``scalings=dict(eeg=1e6, grad=1e13, mag=1e15, misc=1.0)`. use_time_index : bool If False, times will be included as in the data table, else it will be used as index object. copy : bool If true, evoked will be copied. Else data may be modified in place. Returns ------- df : instance of pandas.core.DataFrame Evoked data exported into tabular data structure. """ pd = _check_pandas_installed() if picks is None: picks = list(range(self.info['nchan'])) else: if not in1d(picks, np.arange(len(self.ch_names))).all(): raise ValueError('At least one picked channel is not present ' 'in this Evoked instance.') data, times = self.data, self.times if copy is True: data = data.copy() types = [channel_type(self.info, idx) for idx in picks] n_channel_types = 0 ch_types_used = [] scalings = _mutable_defaults(('scalings', scalings))[0] for t in scalings.keys(): if t in types: n_channel_types += 1 ch_types_used.append(t) for t in ch_types_used: scaling = scalings[t] idx = [picks[i] for i in range(len(picks)) if types[i] == t] if len(idx) > 0: data[idx] *= scaling assert times.shape[0] == data.shape[1] col_names = [self.ch_names[k] for k in picks] df = pd.DataFrame(data.T, columns=col_names) df.insert(0, 'time', times * scale_time) if use_time_index is True: if 'time' in df: df['time'] = df['time'].astype(np.int64) with warnings.catch_warnings(record=True): df.set_index('time', inplace=True) return df def resample(self, sfreq, npad=100, window='boxcar'): """Resample data This function operates in-place. Parameters ---------- sfreq : float New sample rate to use npad : int Amount to pad the start and end of the data. window : string or tuple Window to use in resampling. See scipy.signal.resample. """ o_sfreq = self.info['sfreq'] self.data = resample(self.data, sfreq, o_sfreq, npad, -1, window) # adjust indirectly affected variables self.info['sfreq'] = sfreq self.times = (np.arange(self.data.shape[1], dtype=np.float) / sfreq + self.times[0]) self.first = int(self.times[0] * self.info['sfreq']) self.last = len(self.times) + self.first - 1 def detrend(self, order=1, picks=None): """Detrend data This function operates in-place. Parameters ---------- order : int Either 0 or 1, the order of the detrending. 0 is a constant (DC) detrend, 1 is a linear detrend. picks : array-like of int | None If None only MEG and EEG channels are detrended. """ if picks is None: picks = pick_types(self.info, meg=True, eeg=True, ref_meg=False, stim=False, eog=False, ecg=False, emg=False, exclude='bads') self.data[picks] = detrend(self.data[picks], order, axis=-1) def copy(self): """Copy the instance of evoked Returns ------- evoked : instance of Evoked """ evoked = deepcopy(self) return evoked def __add__(self, evoked): """Add evoked taking into account number of epochs""" out = merge_evoked([self, evoked]) out.comment = self.comment + " + " + evoked.comment return out def __sub__(self, evoked): """Add evoked taking into account number of epochs""" this_evoked = deepcopy(evoked) this_evoked.data *= -1. out = merge_evoked([self, this_evoked]) out.comment = self.comment + " - " + this_evoked.comment return out def __hash__(self): return object_hash(dict(info=self.info, data=self.data)) def get_peak(self, ch_type=None, tmin=None, tmax=None, mode='abs', time_as_index=False): """Get location and latency of peak amplitude Parameters ---------- ch_type : {'mag', 'grad', 'eeg', 'misc', None} The channel type to use. Defaults to None. If more than one sensor Type is present in the data the channel type has to be explicitly set. tmin : float | None The minimum point in time to be considered for peak getting. tmax : float | None The maximum point in time to be considered for peak getting. mode : {'pos', 'neg', 'abs'} How to deal with the sign of the data. If 'pos' only positive values will be considered. If 'neg' only negative values will be considered. If 'abs' absolute values will be considered. Defaults to 'abs'. time_as_index : bool Whether to return the time index instead of the latency in seconds. Returns ------- ch_name : str The channel exhibiting the maximum response. latency : float | int The time point of the maximum response, either latency in seconds or index. """ supported = ('mag', 'grad', 'eeg', 'misc', 'None') data_picks = pick_types(self.info, meg=True, eeg=True, ref_meg=False) types_used = set([channel_type(self.info, idx) for idx in data_picks]) if str(ch_type) not in supported: raise ValueError('Channel type must be `{supported}`. You gave me ' '`{ch_type}` instead.' .format(ch_type=ch_type, supported='` or `'.join(supported))) elif ch_type is not None and ch_type not in types_used: raise ValueError('Channel type `{ch_type}` not found in this ' 'evoked object.' .format(ch_type=ch_type)) elif len(types_used) > 1 and ch_type is None: raise RuntimeError('More than one sensor type found. `ch_type` ' 'must not be `None`, pass a sensor type ' 'value instead') meg, eeg, misc, picks = False, False, False, None if ch_type == 'mag': meg = ch_type elif ch_type == 'grad': meg = ch_type elif ch_type == 'eeg': eeg = True elif ch_type == 'misc': misc = True if ch_type is not None: picks = pick_types(self.info, meg=meg, eeg=eeg, misc=misc, ref_meg=False) data = self.data if picks is None else self.data[picks] ch_idx, time_idx = _get_peak(data, self.times, tmin, tmax, mode) return (self.ch_names[ch_idx], time_idx if time_as_index else self.times[time_idx]) class EvokedArray(Evoked): """Evoked object from numpy array Parameters ---------- data : array of shape (n_channels, n_times) The channels' evoked response. info : instance of Info Info dictionary. Consider using ``create_info`` to populate this structure. tmin : float Start time before event. comment : string Comment on dataset. Can be the condition. Defaults to ''. nave : int Number of averaged epochs. Defaults to 1. kind : str Type of data, either average or standard_error. Defaults to 'average'. verbose : bool, str, int, or None If not None, override default verbose level (see mne.verbose). Defaults to raw.verbose. """ @verbose def __init__(self, data, info, tmin, comment='', nave=1, kind='average', verbose=None): dtype = np.complex128 if np.any(np.iscomplex(data)) else np.float64 data = np.asanyarray(data, dtype=dtype) if data.ndim != 2: raise ValueError('Data must be a 2D array of shape (n_channels, ' 'n_samples)') if len(info['ch_names']) != np.shape(data)[0]: raise ValueError('Info and data must have same number of ' 'channels.') self.data = data # XXX: this should use round and be tested self.first = int(tmin * info['sfreq']) self.last = self.first + np.shape(data)[-1] - 1 self.times = np.arange(self.first, self.last + 1, dtype=np.float) self.times /= info['sfreq'] self.info = info self.nave = nave self.kind = kind self.comment = comment self.proj = None self.picks = None self.verbose = verbose self._projector = None if self.kind == 'average': self._aspect_kind = aspect_dict['average'] else: self._aspect_kind = aspect_dict['standard_error'] def _get_entries(fid, evoked_node): """Helper to get all evoked entries""" comments = list() aspect_kinds = list() for ev in evoked_node: for k in range(ev['nent']): my_kind = ev['directory'][k].kind pos = ev['directory'][k].pos if my_kind == FIFF.FIFF_COMMENT: tag = read_tag(fid, pos) comments.append(tag.data) my_aspect = dir_tree_find(ev, FIFF.FIFFB_ASPECT)[0] for k in range(my_aspect['nent']): my_kind = my_aspect['directory'][k].kind pos = my_aspect['directory'][k].pos if my_kind == FIFF.FIFF_ASPECT_KIND: tag = read_tag(fid, pos) aspect_kinds.append(int(tag.data)) comments = np.atleast_1d(comments) aspect_kinds = np.atleast_1d(aspect_kinds) if len(comments) != len(aspect_kinds) or len(comments) == 0: fid.close() raise ValueError('Dataset names in FIF file ' 'could not be found.') t = [aspect_rev.get(str(a), 'Unknown') for a in aspect_kinds] t = ['"' + c + '" (' + tt + ')' for tt, c in zip(t, comments)] t = ' ' + '\n '.join(t) return comments, aspect_kinds, t def _get_evoked_node(fname): """Helper to get info in evoked file""" f, tree, _ = fiff_open(fname) with f as fid: _, meas = read_meas_info(fid, tree) evoked_node = dir_tree_find(meas, FIFF.FIFFB_EVOKED) return evoked_node def merge_evoked(all_evoked): """Merge/concat evoked data Data should have the same channels and the same time instants. Parameters ---------- all_evoked : list of Evoked The evoked datasets Returns ------- evoked : Evoked The merged evoked data """ evoked = deepcopy(all_evoked[0]) ch_names = evoked.ch_names for e in all_evoked[1:]: assert e.ch_names == ch_names, ValueError("%s and %s do not contain " "the same channels" % (evoked, e)) assert np.max(np.abs(e.times - evoked.times)) < 1e-7, \ ValueError("%s and %s do not contain the same time " "instants" % (evoked, e)) # use union of bad channels bads = list(set(evoked.info['bads']).union(*(ev.info['bads'] for ev in all_evoked[1:]))) evoked.info['bads'] = bads all_nave = sum(e.nave for e in all_evoked) evoked.data = sum(e.nave * e.data for e in all_evoked) / all_nave evoked.nave = all_nave return evoked @deprecated("'read_evoked' will be removed in v0.9. Use 'read_evokeds.'") def read_evoked(fname, setno=None, baseline=None, kind='average', proj=True): """Read an evoked dataset Parameters ---------- fname : string The file name. setno : int or str | list of int or str | None The index or list of indices of the evoked dataset to read. FIF files can contain multiple datasets. If None and there is only one dataset in the file, this dataset is loaded. baseline : None (default) or tuple of length 2 The time interval to apply baseline correction. If None do not apply it. If baseline is (a, b) the interval is between "a (s)" and "b (s)". If a is None the beginning of the data is used and if b is None then b is set to the end of the interval. If baseline is equal ot (None, None) all the time interval is used. kind : str Either 'average' or 'standard_error', the type of data to read. proj : bool If False, available projectors won't be applied to the data. Returns ------- evoked : instance of Evoked or list of Evoked The evoked datasets. """ evoked_node = _get_evoked_node(fname) if setno is None and len(evoked_node) > 1: fid, _, _ = fiff_open(fname) try: _, _, t = _get_entries(fid, evoked_node) except: t = 'None found, must use integer' else: fid.close() raise ValueError('%d datasets present, setno parameter must be set.' 'Candidate setno names:\n%s' % (len(evoked_node), t)) elif isinstance(setno, list): return [Evoked(fname, s, baseline=baseline, kind=kind, proj=proj) for s in setno] else: if setno is None: setno = 0 return Evoked(fname, setno, baseline=baseline, kind=kind, proj=proj) @verbose def read_evokeds(fname, condition=None, baseline=None, kind='average', proj=True, verbose=None): """Read evoked dataset(s) Parameters ---------- fname : string The file name, which should end with -ave.fif or -ave.fif.gz. condition : int or str | list of int or str | None The index or list of indices of the evoked dataset to read. FIF files can contain multiple datasets. If None, all datasets are returned as a list. baseline : None (default) or tuple of length 2 The time interval to apply baseline correction. If None do not apply it. If baseline is (a, b) the interval is between "a (s)" and "b (s)". If a is None the beginning of the data is used and if b is None then b is set to the end of the interval. If baseline is equal ot (None, None) all the time interval is used. kind : str Either 'average' or 'standard_error', the type of data to read. proj : bool If False, available projectors won't be applied to the data. verbose : bool, str, int, or None If not None, override default verbose level (see mne.verbose). Returns ------- evoked : Evoked (if condition is int or str) or list of Evoked (if condition is None or list) The evoked dataset(s). """ check_fname(fname, 'evoked', ('-ave.fif', '-ave.fif.gz')) return_list = True if condition is None: evoked_node = _get_evoked_node(fname) condition = range(len(evoked_node)) elif not isinstance(condition, list): condition = [condition] return_list = False out = [Evoked(fname, c, baseline=baseline, kind=kind, proj=proj, verbose=verbose) for c in condition] return out if return_list else out[0] @deprecated("'write_evoked' will be removed in v0.9. Use 'write_evokeds.'") def write_evoked(fname, evoked): """Write an evoked dataset to a file Parameters ---------- fname : string The file name. evoked : instance of Evoked, or list of Evoked The evoked dataset to save, or a list of evoked datasets to save in one file. Note that the measurement info from the first evoked instance is used, so be sure that information matches. """ if not isinstance(evoked, list): evoked = [evoked] # Create the file and save the essentials fid = start_file(fname) start_block(fid, FIFF.FIFFB_MEAS) write_id(fid, FIFF.FIFF_BLOCK_ID) if evoked[0].info['meas_id'] is not None: write_id(fid, FIFF.FIFF_PARENT_BLOCK_ID, evoked[0].info['meas_id']) # Write measurement info write_meas_info(fid, evoked[0].info) # One or more evoked data sets start_block(fid, FIFF.FIFFB_PROCESSED_DATA) for e in evoked: start_block(fid, FIFF.FIFFB_EVOKED) # Comment is optional if e.comment is not None and len(e.comment) > 0: write_string(fid, FIFF.FIFF_COMMENT, e.comment) # First and last sample write_int(fid, FIFF.FIFF_FIRST_SAMPLE, e.first) write_int(fid, FIFF.FIFF_LAST_SAMPLE, e.last) # The epoch itself start_block(fid, FIFF.FIFFB_ASPECT) write_int(fid, FIFF.FIFF_ASPECT_KIND, e._aspect_kind) write_int(fid, FIFF.FIFF_NAVE, e.nave) decal = np.zeros((e.info['nchan'], 1)) for k in range(e.info['nchan']): decal[k] = 1.0 / (e.info['chs'][k]['cal'] * e.info['chs'][k].get('scale', 1.0)) write_float_matrix(fid, FIFF.FIFF_EPOCH, decal * e.data) end_block(fid, FIFF.FIFFB_ASPECT) end_block(fid, FIFF.FIFFB_EVOKED) end_block(fid, FIFF.FIFFB_PROCESSED_DATA) end_block(fid, FIFF.FIFFB_MEAS) end_file(fid) def write_evokeds(fname, evoked): """Write an evoked dataset to a file Parameters ---------- fname : string The file name, which should end with -ave.fif or -ave.fif.gz. evoked : Evoked instance, or list of Evoked instances The evoked dataset, or list of evoked datasets, to save in one file. Note that the measurement info from the first evoked instance is used, so be sure that information matches. """ check_fname(fname, 'evoked', ('-ave.fif', '-ave.fif.gz')) if not isinstance(evoked, list): evoked = [evoked] # Create the file and save the essentials with start_file(fname) as fid: start_block(fid, FIFF.FIFFB_MEAS) write_id(fid, FIFF.FIFF_BLOCK_ID) if evoked[0].info['meas_id'] is not None: write_id(fid, FIFF.FIFF_PARENT_BLOCK_ID, evoked[0].info['meas_id']) # Write measurement info write_meas_info(fid, evoked[0].info) # One or more evoked data sets start_block(fid, FIFF.FIFFB_PROCESSED_DATA) for e in evoked: start_block(fid, FIFF.FIFFB_EVOKED) # Comment is optional if e.comment is not None and len(e.comment) > 0: write_string(fid, FIFF.FIFF_COMMENT, e.comment) # First and last sample write_int(fid, FIFF.FIFF_FIRST_SAMPLE, e.first) write_int(fid, FIFF.FIFF_LAST_SAMPLE, e.last) # The epoch itself start_block(fid, FIFF.FIFFB_ASPECT) write_int(fid, FIFF.FIFF_ASPECT_KIND, e._aspect_kind) write_int(fid, FIFF.FIFF_NAVE, e.nave) decal = np.zeros((e.info['nchan'], 1)) for k in range(e.info['nchan']): decal[k] = 1.0 / (e.info['chs'][k]['cal'] * e.info['chs'][k].get('scale', 1.0)) write_float_matrix(fid, FIFF.FIFF_EPOCH, decal * e.data) end_block(fid, FIFF.FIFFB_ASPECT) end_block(fid, FIFF.FIFFB_EVOKED) end_block(fid, FIFF.FIFFB_PROCESSED_DATA) end_block(fid, FIFF.FIFFB_MEAS) end_file(fid) def _get_peak(data, times, tmin=None, tmax=None, mode='abs'): """Get feature-index and time of maximum signal from 2D array Note. This is a 'getter', not a 'finder'. For non-evoked type data and continuous signals, please use proper peak detection algorithms. Parameters ---------- data : instance of numpy.ndarray (n_locations, n_times) The data, either evoked in sensor or source space. times : instance of numpy.ndarray (n_times) The times in seconds. tmin : float | None The minimum point in time to be considered for peak getting. tmax : float | None The maximum point in time to be considered for peak getting. mode : {'pos', 'neg', 'abs'} How to deal with the sign of the data. If 'pos' only positive values will be considered. If 'neg' only negative values will be considered. If 'abs' absolute values will be considered. Defaults to 'abs'. Returns ------- max_loc : int The index of the feature with the maximum value. max_time : int The time point of the maximum response, index. """ modes = ('abs', 'neg', 'pos') if mode not in modes: raise ValueError('The `mode` parameter must be `{modes}`. You gave ' 'me `{mode}`'.format(modes='` or `'.join(modes), mode=mode)) if tmin is None: tmin = times[0] if tmax is None: tmax = times[-1] if tmin < times.min(): raise ValueError('The tmin value is out of bounds. It must be ' 'within {0} and {1}'.format(times.min(), times.max())) if tmax > times.max(): raise ValueError('The tmin value is out of bounds. It must be ' 'within {0} and {1}'.format(times.min(), times.max())) if tmin >= tmax: raise ValueError('The tmin must be smaller than tmax') time_win = (times >= tmin) & (times <= tmax) mask = np.ones_like(data).astype(np.bool) mask[:, time_win] = False maxfun = np.argmax if mode == 'pos': if not np.any(data > 0): raise ValueError('No positive values encountered. Cannot ' 'operate in pos mode.') elif mode == 'neg': if not np.any(data < 0): raise ValueError('No negative values encountered. Cannot ' 'operate in neg mode.') maxfun = np.argmin masked_index = np.ma.array(np.abs(data) if mode == 'abs' else data, mask=mask) max_loc, max_time = np.unravel_index(maxfun(masked_index), data.shape) return max_loc, max_time