From b14e1747017da80d0ded0fdd3f3f3fd489460515 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Thu, 22 Jul 2021 14:43:10 -0400 Subject: [PATCH 1/3] ENH: Add copy-and-skip for EEG-only data --- mnefun/_sss.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/mnefun/_sss.py b/mnefun/_sss.py index b43abb3..3ebab6e 100644 --- a/mnefun/_sss.py +++ b/mnefun/_sss.py @@ -49,6 +49,24 @@ def mark_flat(raw, bad_percent=5., min_duration=0.005, picks=None, return raw +def _check_skip(p, subj, run_idx): + raw_files = get_raw_fnames(p, subj, 'raw', erm=False, + run_indices=run_idx) + if all('meg' in read_raw_fif(fname, allow_maxshield='yes') + for fname in raw_files): + return False # don't skip SSS + print(f' Skipping SSS for {subj} as no MEG data found') + # copy files + raw_files_out = get_raw_fnames(p, subj, 'sss', erm=False, + run_indices=run_idx) + raw_files += get_raw_fnames(p, subj, 'raw', 'only') + raw_files_out += get_raw_fnames(p, subj, 'sss', 'only') + os.makedirs(op.dirname(raw_files_out[0]), exist_ok=True) + for a, b in zip(raw_files, raw_files_out): + shutil.copyfile(a, b) + return True + + def run_sss(p, subjects, run_indices): """Run SSS preprocessing remotely (only designed for *nix platforms) or locally using Maxwell filtering in mne-python""" @@ -58,6 +76,8 @@ def run_sss(p, subjects, run_indices): run_sss_locally(p, subjects, run_indices) else: for si, subj in enumerate(subjects): + if _check_skip(p, subj, run_indices[si]): + continue files = get_raw_fnames(p, subj, 'raw', False, True, run_indices[si]) n_files = len(files) @@ -414,6 +434,8 @@ def run_sss_locally(p, subjects, run_indices): reg = p.sss_regularize for si, subj in enumerate(subjects): + if _check_skip(subj, subj, run_indices[si]): + continue if p.disp_files: print(' Maxwell filtering subject %g/%g (%s).' % (si + 1, len(subjects), subj)) From 3ebd1b18d98e04ef6fe9fe89c31d5265d1fe9147 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Mon, 26 Jul 2021 12:34:45 -0400 Subject: [PATCH 2/3] ENH: Add last parts for EEG-only --- doc/overview.rst | 27 ++++++++++++++-- mnefun/_cov.py | 37 +++++++++++----------- mnefun/_fix.py | 21 +++++++------ mnefun/_forward.py | 8 ++--- mnefun/_inverse.py | 13 +++++--- mnefun/_mnefun.py | 5 +++ mnefun/_paths.py | 23 +++++++++++--- mnefun/_report.py | 66 +++++++++++++++++++++++---------------- mnefun/_ssp.py | 22 ++++++++----- mnefun/_status.py | 3 +- mnefun/data/canonical.yml | 5 +++ 11 files changed, 152 insertions(+), 78 deletions(-) diff --git a/doc/overview.rst b/doc/overview.rst index 6cd543a..ff5eacc 100644 --- a/doc/overview.rst +++ b/doc/overview.rst @@ -120,6 +120,10 @@ subject_run_indices : list of array-like | dict | None (must be same length as ``params.subjects``) or a dict (keys are subject strings, values are the run indices) where missing subjects get all runs. None is an alias for "all runs". +fix_eeg_order : bool | dict + Whether or not to fix the EEG order for a given subject (dict) or all + subjects (bool). Default is True, which is appropriate for use with the + mis-ordered Neuromag caps. 2. do_score ----------- @@ -318,6 +322,13 @@ fir_window : str See :func:`mne.filter.create_filter`. phase : str See :func:`mne.filter.create_filter`. +notch_filter : None | dict + If None (default), perform no FIR-based notch filtering (MEG data will + by default still have a notch filter applied during the Maxwell filtering + step). If ``dict``, arguments are passed directly to + :meth:`mne.io.Raw.notch_filter`, e.g.:: + + params.notch_filter = dict(freqs=[60], notch_widths=1, trans_bandwidth=0.5) .. _preprocessing_auto_bads: @@ -383,6 +394,10 @@ ssp_eog_reject : dict | None ssp_ecg_reject : dict | None Amplitude rejection criteria for ECG SSP computation. None will use the mne-python default. +ssp_eog_baseline : None | tuple + The baseline to use when constructing EOG epochs (default None). +ssp_ecg_baseline : None | tuple + The baseline to use when constructing ECG epochs (default None). eog_channel : str | dict | None The channel to use to detect blink events. None will use EOG* channels. In lieu of an EOG recording, MEG1411 may work. @@ -564,7 +579,10 @@ force_erm_cov_rank_full : bool If True, force the ERM cov to be full rank. Usually not needed, but might help when the empty-room data is short and/or there are a lot of head movements. - +erm_cov_from_task : bool + If True (default False), use the task runs to compute the "empty-room" + covariance rather than any empty-room run(s). This should only be used + with EEG data, and should be interpreted with caution. 9. gen_fwd ---------- @@ -625,8 +643,11 @@ good_hpi_count : bool Number of good HPI coils (default True). head_movement : bool Head movement (default True). -raw_segments : bool - 10 evenly spaced raw data segments (default True). +raw_segments : bool | dict + 10 evenly spaced raw data segments (default True). Can also be a dict + of keyword arguments to pass to :meth:`mne.Epochs.plot`, e.g.:: + + raw_segments=dict(scalings=dict(eeg=100e-6)) psd : bool Raw PSDs, often slow (default True). ssp_topomaps : bool diff --git a/mnefun/_cov.py b/mnefun/_cov.py index 27bf454..8ec4ef8 100644 --- a/mnefun/_cov.py +++ b/mnefun/_cov.py @@ -14,7 +14,8 @@ from mne.viz import plot_cov from ._epoching import _concat_resamp_raws -from ._paths import get_epochs_evokeds_fnames, get_raw_fnames, safe_inserter +from ._paths import (get_epochs_evokeds_fnames, get_raw_fnames, safe_inserter, + get_erm_cov_fname) from ._scoring import _read_events from ._utils import (get_args, _get_baseline, _restrict_reject_flat, _handle_dict, _handle_decim) @@ -72,11 +73,9 @@ def gen_covariances(p, subjects, run_indices, decim): for si, subj in enumerate(subjects): print(' Subject %2d/%2d...' % (si + 1, len(subjects)), end='') cov_dir = op.join(p.work_dir, subj, p.cov_dir) - if not op.isdir(cov_dir): - os.mkdir(cov_dir) + os.makedirs(cov_dir, exist_ok=True) has_rank_arg = 'rank' in get_args(compute_covariance) - kwargs = dict() - kwargs_erm = dict() + kwargs = dict(method=p.cov_method) if p.cov_rank == 'full': # backward compat if has_rank_arg: kwargs['rank'] = 'full' @@ -90,7 +89,7 @@ def gen_covariances(p, subjects, run_indices, decim): kwargs['rank'] = _compute_rank(p, subj, run_indices[si]) else: kwargs['rank'] = p.cov_rank - kwargs_erm['rank'] = kwargs['rank'] + kwargs_erm = kwargs.copy() if p.force_erm_cov_rank_full and has_rank_arg: kwargs_erm['rank'] = 'full' # Use the same thresholds we used for primary Epochs @@ -103,22 +102,25 @@ def gen_covariances(p, subjects, run_indices, decim): flat = _handle_dict(p.flat, subj) # Make empty room cov - if p.runs_empty: - if len(p.runs_empty) > 1: - raise ValueError('Too many empty rooms; undefined output!') - new_run = safe_inserter(p.runs_empty[0], subj) - empty_cov_name = op.join(cov_dir, new_run + p.pca_extra + - p.inv_tag + '-cov.fif') - empty_fif = get_raw_fnames(p, subj, 'pca', 'only', False)[0] - raw = read_raw_fif(empty_fif, preload=True) - raw.pick_types(meg=True, eog=True, exclude='bads') + empty_cov_name = get_erm_cov_fname(p, subj) + if empty_cov_name is not None: + if p.erm_cov_from_task: + # read in raw files + raw_fnames = get_raw_fnames( + p, subj, 'pca', False, False, run_indices[si]) + raw, ratios = _concat_resamp_raws(p, subj, raw_fnames) + raw.pick_types(meg=True, eeg=True, eog=True, exclude='bads') + else: + empty_fif = get_raw_fnames(p, subj, 'pca', 'only', False)[0] + raw = read_raw_fif(empty_fif, preload=True) + raw.pick_types(meg=True, eog=True, exclude='bads') use_reject, use_flat = _restrict_reject_flat(reject, flat, raw) if 'eeg' in use_reject: del use_reject['eeg'] if 'eeg' in use_flat: del use_flat['eeg'] cov = compute_raw_covariance(raw, reject=use_reject, flat=use_flat, - method=p.cov_method, **kwargs_erm) + **kwargs_erm) write_cov(empty_cov_name, cov) # Make evoked covariances @@ -153,8 +155,7 @@ def gen_covariances(p, subjects, run_indices, decim): on_missing=p.on_missing, reject_by_annotation=p.reject_epochs_by_annot) epochs.pick_types(meg=True, eeg=True, exclude=[]) - cov = compute_covariance(epochs, method=p.cov_method, - **kwargs) + cov = compute_covariance(epochs, **kwargs) if kwargs.get('rank', None) not in (None, 'full'): want_rank = sum(kwargs['rank'].values()) out_rank = compute_whitener( diff --git a/mnefun/_fix.py b/mnefun/_fix.py index ba25adf..5bf3fe7 100644 --- a/mnefun/_fix.py +++ b/mnefun/_fix.py @@ -13,6 +13,7 @@ from mne import pick_types from ._paths import get_raw_fnames +from ._utils import _handle_dict _1_ORDER = (1, 2, 3, 5, 6, 7, 9, 10, @@ -46,6 +47,7 @@ def fix_eeg_files(p, subjects, structurals=None, dates=None, run_indices=None): if run_indices is None: run_indices = [None] * len(subjects) for si, subj in enumerate(subjects): + fix_eeg = _handle_dict(p.fix_eeg_order, subj) if p.disp_files: print(' Fixing subject %g/%g.' % (si + 1, len(subjects))) raw_names = get_raw_fnames(p, subj, 'sss', True, False, @@ -64,10 +66,10 @@ def fix_eeg_files(p, subjects, structurals=None, dates=None, run_indices=None): birthday=dates[si]) else: anon = None - fix_eeg_channels(names, anon) + fix_eeg_channels(names, anon, fix_eeg=fix_eeg) -def fix_eeg_channels(raw_files, anon=None, verbose=True): +def fix_eeg_channels(raw_files, anon=None, fix_eeg=True, verbose=True): """Reorder EEG channels based on UW cap setup Parameters @@ -88,7 +90,7 @@ def fix_eeg_channels(raw_files, anon=None, verbose=True): # actually do the reordering for ri, raw_file in enumerate(raw_files): need_reorder, need_anon, write_key, anon_key, picks, order = \ - _is_file_unfixed(raw_file, anon) + _is_file_unfixed(raw_file, anon, fix_eeg) if need_anon or need_reorder: to_do = [] if need_reorder: @@ -123,13 +125,14 @@ def fix_eeg_channels(raw_files, anon=None, verbose=True): print(' File %i already corrected' % (ri + 1)) -def _all_files_fixed(p, subj, type_='pca'): - """Determine if all files have been fixed for a subject""" - return all(op.isfile(fname) and not any(_is_file_unfixed(fname)[:2]) +def _all_files_fixed(p, subj, type_='pca', fix_eeg=True): + """Determine if all files have been fixed for a subject.""" + kw = dict(fix_eeg=_handle_dict(p.fix_eeg_order, subj)) + return all(op.isfile(fname) and not any(_is_file_unfixed(fname, **kw)[:2]) for fname in get_raw_fnames(p, subj, type_)) -def _is_file_unfixed(fname, anon=None): +def _is_file_unfixed(fname, anon=None, fix_eeg=True): """Determine if a file needs reordering or anonymization.""" order = np.array(_1_ORDER) - 1 assert len(order) == 60 @@ -142,7 +145,7 @@ def _is_file_unfixed(fname, anon=None): raw = read_raw_fif(fname, preload=False, allow_maxshield='yes') picks = pick_types(raw.info, meg=False, eeg=True, exclude=[]) need_reorder = False - if len(picks) > 0: + if len(picks) > 0 and fix_eeg: if len(picks) == len(order): need_reorder = (write_key not in raw.info['description']) else: @@ -152,7 +155,7 @@ def _is_file_unfixed(fname, anon=None): warnings.warn(msg) else: raise RuntimeError(msg) - need_anon = (anon_key not in raw.info['description']) + need_anon = (anon_key not in (raw.info['description'] or '')) return need_reorder, need_anon, write_key, anon_key, picks, order diff --git a/mnefun/_forward.py b/mnefun/_forward.py index 7300bc5..408b0d1 100644 --- a/mnefun/_forward.py +++ b/mnefun/_forward.py @@ -46,8 +46,8 @@ def gen_forwards(p, subjects, structurals, run_indices): if not getattr(p, 'translate_positions', True): raise RuntimeError('Not translating positions is no longer ' 'supported') - print(' Creating forward solution(s) using a %s for %s...' - % (bem_type, subj), end='') + print(f' Creating forward solution(s) using a {bem_type} for ' + f'{subj} ({struc})...', end='') # XXX Don't actually need to generate a different fwd for each inv # anymore, since all runs are included, but changing the filename # would break a lot of existing pipelines :( @@ -110,8 +110,8 @@ def _get_bem_src_trans(p, info, subj, struc): if op.isfile(src_space_file): break else: # if neither exists, use last filename - print(' Creating %s%s source space for %s...' - % (kind, num, subj)) + print(f' Creating {kind}{num} source space for {subj} ' + f'({struc})...') if kind == 'oct': src = setup_source_space( struc, spacing='%s%s' % (kind, num), diff --git a/mnefun/_inverse.py b/mnefun/_inverse.py index 7fa8c61..ad4c4dc 100644 --- a/mnefun/_inverse.py +++ b/mnefun/_inverse.py @@ -15,7 +15,7 @@ from ._cov import _compute_rank from ._paths import (get_epochs_evokeds_fnames, safe_inserter, - get_cov_fwd_inv_fnames) + get_cov_fwd_inv_fnames, get_erm_cov_fname) try: from mne import spatial_src_adjacency @@ -46,7 +46,6 @@ def gen_inverses(p, subjects, run_indices): cov_dir = op.join(p.work_dir, subj, p.cov_dir) if not op.isdir(inv_dir): os.mkdir(inv_dir) - make_erm_inv = len(p.runs_empty) > 0 epochs_fnames, _ = get_epochs_evokeds_fnames(p, subj, p.analyses) _, fif_file = epochs_fnames @@ -71,11 +70,12 @@ def gen_inverses(p, subjects, run_indices): rank = _compute_rank(p, subj, run_indices[si]) else: rank = None # should be safe from our gen_covariances step + erm_name = get_erm_cov_fname(p, subj) + empty_cov = None + make_erm_inv = erm_name is not None if make_erm_inv: # We now process the empty room with "movement # compensation" so it should get the same rank! - erm_name = op.join(cov_dir, safe_inserter(p.runs_empty[0], subj) + - p.pca_extra + p.inv_tag + '-cov.fif') empty_cov = read_cov(erm_name) if p.force_erm_cov_rank_full and p.cov_method == 'empirical': empty_cov = regularize( @@ -100,6 +100,7 @@ def gen_inverses(p, subjects, run_indices): temp_name = subj cov = empty_cov tag = p.inv_erm_tag + assert cov is not None else: s_name = safe_inserter(name, subj) temp_name = s_name + ('-%d' % p.lp_cut) + p.inv_tag @@ -117,7 +118,9 @@ def gen_inverses(p, subjects, run_indices): inv_dir, temp_name + f + tag + s + '-inv.fif') kwargs = dict(loose=l_, depth=d, fixed=x, use_cps=True, verbose='error') - if name is not True or not e: + if (name is not True or + (not e) or + p.erm_cov_from_task): inv = make_inverse_operator( epochs.info, fwd_restricted, cov, rank=rank, **kwargs) diff --git a/mnefun/_mnefun.py b/mnefun/_mnefun.py index 07317a6..cb64944 100644 --- a/mnefun/_mnefun.py +++ b/mnefun/_mnefun.py @@ -274,6 +274,11 @@ def __init__(self, tmin=None, tmax=None, t_adjust=0, bmin=-0.2, bmax=0.0, self.cont_as_esss = False self.cont_reject = None self.allow_resample = False + self.fix_eeg_order = True + self.erm_cov_from_task = False + self.notch_filter = None + self.ssp_ecg_baseline = None + self.ssp_eog_baseline = None self.freeze() # Read static-able paraws from config file _set_static(self) diff --git a/mnefun/_paths.py b/mnefun/_paths.py index 73272b9..cf2f186 100644 --- a/mnefun/_paths.py +++ b/mnefun/_paths.py @@ -135,7 +135,8 @@ def get_cov_fwd_inv_fnames(p, subj, run_indices): inv_dir = op.join(p.work_dir, subj, p.inverse_dir) fwd_dir = op.join(p.work_dir, subj, p.forward_dir) cov_dir = op.join(p.work_dir, subj, p.cov_dir) - make_erm_inv = len(p.runs_empty) > 0 + erm_cov_fname = get_erm_cov_fname(p, subj) + make_erm_inv = erm_cov_fname is not None # Shouldn't matter which raw file we use raw_fname = get_raw_fnames(p, subj, 'pca', True, False, run_indices)[0] @@ -169,11 +170,10 @@ def get_cov_fwd_inv_fnames(p, subj, run_indices): inv_fnames += [op.join(inv_dir, temp_name + f + s + '-inv.fif')] if make_erm_inv: - cov_fnames += [op.join(cov_dir, safe_inserter(p.runs_empty[0], subj) + - p.pca_extra + p.inv_tag + '-cov.fif')] + cov_fnames += [erm_cov_fname] for f, m, e in zip(out_flags, meg_bools, eeg_bools): for s in [p.inv_fixed_tag, '']: - if (not e): + if (p.erm_cov_from_task or not e): inv_fnames += [op.join(inv_dir, subj + f + p.inv_erm_tag + s + '-inv.fif')] return cov_fnames, fwd_fnames, inv_fnames @@ -242,3 +242,18 @@ def _is_dir(d): def _get_config_file(): return op.expanduser(op.join('~', '.mnefun', 'mnefun.json')) + + +def get_erm_cov_fname(p, subj): + """Get the continuous (empty-room) covariance filename.""" + cov_dir = op.join(p.work_dir, subj, p.cov_dir) + fname = None + if p.erm_cov_from_task: + fname = f'{subj}-cont' + elif len(p.runs_empty): + if len(p.runs_empty) > 1: + raise ValueError('Too many empty rooms; undefined output!') + fname = safe_inserter(p.runs_empty[0], subj) + if fname is not None: + fname = op.join(cov_dir, f'{fname}{p.pca_extra}{p.inv_tag}-cov.fif') + return fname diff --git a/mnefun/_report.py b/mnefun/_report.py index ee19abe..6e594a6 100644 --- a/mnefun/_report.py +++ b/mnefun/_report.py @@ -13,6 +13,7 @@ import mne from mne import read_proj, read_epochs, find_events, pick_types +from mne.channels.channels import _get_ch_info from mne.io import BaseRaw from mne.viz import (plot_projs_topomap, plot_cov, plot_snr_estimate, plot_events) @@ -25,7 +26,8 @@ from ._epoching import _concat_resamp_raws from ._forward import _get_bem_src_trans from ._paths import (get_raw_fnames, get_proj_fnames, get_report_fnames, - get_bad_fname, get_epochs_evokeds_fnames, safe_inserter) + get_bad_fname, get_epochs_evokeds_fnames, safe_inserter, + get_cov_fwd_inv_fnames) from ._ssp import _proj_nums from ._sss import (_load_trans_to, _read_raw_prebad, _get_t_window, _get_fit_data) @@ -179,7 +181,7 @@ def _report_events(report, fnames, p=None, subj=None): print('%5.1f sec' % ((time.time() - t0),)) -def _report_raw_segments(report, raw, lowpass=None): +def _report_raw_segments(report, raw, lowpass=None, **kwargs): t0 = time.time() section = 'Raw segments' print((' %s ... ' % section).ljust(LJUST), end='') @@ -202,9 +204,14 @@ def _report_raw_segments(report, raw, lowpass=None): new_events = np.array([new_events, np.zeros_like(new_events), np.ones_like(new_events)]).T + if _get_ch_info(raw_plot.info)[1]: # has_vv_grad + group_by = 'selection' + else: + group_by = 'position' with mne.utils.use_log_level('error'): - fig = raw_plot.plot(group_by='selection', butterfly=True, - events=new_events, lowpass=lowpass) + fig = raw_plot.plot(group_by=group_by, butterfly=True, + events=new_events, lowpass=lowpass, + **kwargs) fig.axes[0].lines[-1].set_zorder(10) # events fig.axes[0].set(xticks=np.arange(0, len(times)) + 0.5) xticklabels = ['%0.1f' % t for t in times] @@ -344,22 +351,26 @@ def gen_html_report(p, subjects, structurals, run_indices=None): for si, subj in enumerate(subjects): struc = structurals[si] if structurals is not None else None report = Report(verbose=False) + run_idx = run_indices[si] print(' Processing subject %s/%s (%s)' % (si + 1, len(subjects), subj)) # raw fnames = get_raw_fnames(p, subj, 'raw', erm=False, add_splits=False, - run_indices=run_indices[si]) + run_indices=run_idx) for fname in fnames: if not op.isfile(fname): raise RuntimeError('Cannot create reports until raw data ' 'exist, missing:\n%s' % fname) raw, _ = _concat_resamp_raws( p, subj, fnames, fix='all', prebad=True, preload=preload) + raw.info['projs'] = [] + if 'eeg' in raw: + raw.set_eeg_reference(projection=True, verbose='error') + raw.apply_proj() # sss - sss_fnames = get_raw_fnames(p, subj, 'sss', False, False, - run_indices[si]) + sss_fnames = get_raw_fnames(p, subj, 'sss', False, False, run_idx) has_sss = all(op.isfile(fname) for fname in sss_fnames) sss_info = mne.io.read_raw_fif(sss_fnames[0]) if has_sss else None bad_file = get_bad_fname(p, subj) @@ -375,7 +386,7 @@ def gen_html_report(p, subjects, structurals, run_indices=None): ('raw_erm', 'raw', 'only'), ('raw_erm_pca', 'pca', 'only')]: these_fnames = get_raw_fnames( - p, subj, which, erm, False, run_indices[si]) + p, subj, which, erm, False, run_idx) if len(these_fnames) and all(op.isfile(f) for f in these_fnames): extra_raws[key], _ = _concat_resamp_raws( p, subj, these_fnames, 'all', preload=True) @@ -427,16 +438,18 @@ def gen_html_report(p, subjects, structurals, run_indices=None): # Head movement # if p.report_params.get('head_movement', True) and p.movecomp: - _report_head_movement(report, fnames, p, subj, run_indices[si]) + _report_head_movement(report, fnames, p, subj, run_idx) else: print(' Head movement skipped') # # Raw segments # - if p.report_params.get('raw_segments', True) and \ - raw_pca is not None: - _report_raw_segments(report, raw_pca) + raw_segments = p.report_params.get('raw_segments', True) + if raw_segments and raw_pca is not None: + kw = {} if raw_segments is True else raw_segments + assert isinstance(kw, dict), type(raw_segments) + _report_raw_segments(report, raw_pca, lowpass=None, **kw) else: print(' Raw segments skipped') @@ -651,7 +664,7 @@ def gen_html_report(p, subjects, structurals, run_indices=None): t0 = time.time() print((' %s ... ' % section).ljust(LJUST), end='') cov_name = p.report_params.get('covariance', None) - cov_name = _get_cov_name(p, subj, cov_name) + cov_name = _get_cov_name(p, subj, cov_name, run_idx) if cov_name is None: print(' Missing covariance: %s' % op.basename(cov_name), end='') @@ -683,7 +696,8 @@ def gen_html_report(p, subjects, structurals, run_indices=None): assert isinstance(whitening, dict) analysis = whitening['analysis'] name = whitening['name'] - cov_name = _get_cov_name(p, subj, whitening.get('cov')) + cov_name = _get_cov_name( + p, subj, whitening.get('cov'), run_idx) # Load the inverse fname_evoked = op.join(inv_dir, '%s_%d%s_%s_%s-ave.fif' % (analysis, p.lp_cut, p.inv_tag, @@ -874,7 +888,8 @@ def gen_html_report(p, subjects, structurals, run_indices=None): # Define the time slices to include times = sensor.get('times', [0.1, 0.2]) if isinstance(times, str) and times == 'peaks': - cov_name = _get_cov_name(p, subj, sensor.get('cov')) + cov_name = _get_cov_name( + p, subj, sensor.get('cov'), run_idx) times = _get_memo_times( this_evoked, cov_name, (fname_evoked, name), times_memo) @@ -1028,7 +1043,7 @@ def gen_html_report(p, subjects, structurals, run_indices=None): times = source.get('times', [0.1, 0.2]) if isinstance(times, str) and times == 'peaks': cov_name = _get_cov_name( - p, subj, source.get('cov')) + p, subj, source.get('cov'), run_idx) times = _get_memo_times( this_evoked, cov_name, (fname_evoked, name), times_memo) @@ -1203,7 +1218,8 @@ def _proj_fig(fname, info, proj_nums, proj_meg, kind, use_ch, duration): t.remove() scale = 0.8 * np.abs(ax.get_ylim()).max() hs, labels = list(), list() - traces /= np.abs(traces).max() # uniformly scaled + with np.errstate(invalid='raise'): + traces /= np.abs(traces).max() # uniformly scaled for ti, trace in enumerate(traces): hs.append(ax.plot( this_evoked.times, trace * scale, @@ -1254,18 +1270,14 @@ def _proj_fig(fname, info, proj_nums, proj_meg, kind, use_ch, duration): return fig -def _get_cov_name(p, subj, cov_name=None): +def _get_cov_name(p, subj, cov_name, run_idx): # just the first for now if cov_name is None or cov_name is True: - if p.inv_names: - cov_name = (safe_inserter(p.inv_names[0], subj) + - ('-%d' % p.lp_cut) + p.inv_tag + '-cov.fif') - elif p.runs_empty: # erm cov - new_run = safe_inserter(p.runs_empty[0], subj) - cov_name = new_run + p.pca_extra + p.inv_tag + '-cov.fif' - if cov_name is not None: + cov_name = get_cov_fwd_inv_fnames(p, subj, run_idx)[0][0] + else: # user-provided string (relative path) cov_dir = op.join(p.work_dir, subj, p.cov_dir) cov_name = op.join(cov_dir, safe_inserter(cov_name, subj)) - if not op.isfile(cov_name): - cov_name = None + if not op.isfile(cov_name): + warnings.warn(f'Could not find covariance: {cov_name}') + cov_name = None return cov_name diff --git a/mnefun/_ssp.py b/mnefun/_ssp.py index e4bbb02..e657475 100644 --- a/mnefun/_ssp.py +++ b/mnefun/_ssp.py @@ -41,7 +41,7 @@ def _raw_LRFCP(raw_names, sfreq, l_freq, h_freq, n_jobs, n_jobs_resample, force_bads=False, l_trans=0.5, h_trans=0.5, allow_maxshield=False, phase='zero-double', fir_window='hann', fir_design='firwin2', pick=True, - skip_by_annotation=('bad', 'skip')): + skip_by_annotation=('bad', 'skip'), notch=None): """Helper to load, filter, concatenate, then project raw files""" from mne.io.proj import _needs_eeg_average_ref_proj from ._sss import _read_raw_prebad @@ -71,6 +71,9 @@ def _raw_LRFCP(raw_names, sfreq, l_freq, h_freq, n_jobs, n_jobs_resample, filter_length=filter_length, phase=phase, l_trans_bandwidth=l_trans, h_trans_bandwidth=h_trans, fir_window=fir_window, **fir_kwargs) + if notch is not None: + assert isinstance(notch, dict), type(notch) + r.notch_filter(**notch) raw.append(r) _fix_raw_eog_cals(raw) raws_del = raw[1:] @@ -130,7 +133,7 @@ def _compute_erm_proj(p, subj, projs, kind, bad_file, remove_existing=False, bad_file=bad_file, disp_files=disp_files, method='fir', filter_length=p.filter_length, force_bads=True, l_trans=p.cont_hp_trans, h_trans=p.cont_lp_trans, - phase=p.phase, fir_window=p.fir_window, + phase=p.phase, fir_window=p.fir_window, notch=p.notch_filter, skip_by_annotation='edge', **fir_kwargs) if remove_existing: raw.del_proj() @@ -203,6 +206,7 @@ def do_preprocessing_combined(p, subjects, run_indices): raise NameError('File not found (' + r + ')') fir_kwargs, old_kwargs = _get_fir_kwargs(p.fir_design) + notch = _handle_dict(p.notch_filter, subj) if isinstance(p.auto_bad, float): print(' Creating post SSS bad channel file:\n' ' %s' % bad_file) @@ -214,6 +218,7 @@ def do_preprocessing_combined(p, subjects, run_indices): l_trans=p.hp_trans, h_trans=p.lp_trans, phase=p.phase, fir_window=p.fir_window, pick=True, skip_by_annotation='edge', + notch=notch, **fir_kwargs) events = fixed_len_events(p, raw) rtmin = p.reject_tmin \ @@ -325,7 +330,7 @@ def do_preprocessing_combined(p, subjects, run_indices): method='fir', filter_length=p.filter_length, force_bads=False, l_trans=p.hp_trans, h_trans=p.lp_trans, phase=p.phase, fir_window=p.fir_window, pick=True, skip_by_annotation='edge', - **fir_kwargs) + notch=notch, **fir_kwargs) # Apply any user-supplied extra projectors if p.proj_extra is not None: @@ -382,7 +387,8 @@ def do_preprocessing_combined(p, subjects, run_indices): _handle_dict(p.ssp_ecg_reject, subj), flat, raw) ecg_epochs = Epochs( raw, ecg_events, 999, ecg_t_lims[0], ecg_t_lims[1], - baseline=None, reject=use_reject, flat=use_flat, preload=True) + baseline=p.ssp_ecg_baseline, reject=use_reject, flat=use_flat, + preload=True) print(' obtained %d epochs from %d events.' % (len(ecg_epochs), len(ecg_events))) if len(ecg_epochs) >= 20: @@ -487,7 +493,8 @@ def _compute_add_eog(p, subj, raw_orig, projs, eog_nums, kind, pca_dir, _handle_dict(p.ssp_eog_reject, subj), flat, raw) eog_epochs = Epochs( raw, eog_events, 998, eog_t_lims[0], eog_t_lims[1], - baseline=None, reject=use_reject, flat=use_flat, preload=True) + baseline=p.ssp_eog_baseline, reject=use_reject, flat=use_flat, + preload=True) print(' obtained %d epochs from %d events.' % (len(eog_epochs), len(eog_events))) del eog_events @@ -542,6 +549,7 @@ def apply_preprocessing_combined(p, subjects, run_indices): all_proj = op.join(pca_dir, 'preproc_all-proj.fif') projs = read_proj(all_proj) fir_kwargs = _get_fir_kwargs(p.fir_design)[0] + notch = _handle_dict(p.notch_filter, subj) if len(erm_in) > 0: for ii, (r, o) in enumerate(zip(erm_in, erm_out)): if p.disp_files: @@ -554,7 +562,7 @@ def apply_preprocessing_combined(p, subjects, run_indices): apply_proj=False, filter_length=p.filter_length, force_bads=True, l_trans=p.hp_trans, h_trans=p.lp_trans, phase=p.phase, fir_window=p.fir_window, pick=False, - **fir_kwargs) + notch=notch, **fir_kwargs) raw.save(o, overwrite=True, buffer_size_sec=None) for ii, (r, o) in enumerate(zip(names_in, names_out)): if p.disp_files: @@ -567,7 +575,7 @@ def apply_preprocessing_combined(p, subjects, run_indices): apply_proj=False, filter_length=p.filter_length, force_bads=False, l_trans=p.hp_trans, h_trans=p.lp_trans, phase=p.phase, fir_window=p.fir_window, pick=False, - **fir_kwargs) + notch=notch, **fir_kwargs) raw.save(o, overwrite=True, buffer_size_sec=None) # look at raw_clean for ExG events if p.plot_raw: diff --git a/mnefun/_status.py b/mnefun/_status.py index 9ff3433..5e6a6ef 100644 --- a/mnefun/_status.py +++ b/mnefun/_status.py @@ -49,7 +49,8 @@ def print_proc_status(p, subjects, structurals, analyses, run_indices): # check if coreg has been done if struc is None or \ - op.isfile(op.join(subj, 'trans', subj + '-trans.fif')): + op.isfile(op.join(p.work_dir, subj, 'trans', + f'{subj}-trans.fif')): coreg = 'complete' # check if sss has been fetched (+1 is for erm) diff --git a/mnefun/data/canonical.yml b/mnefun/data/canonical.yml index f1be082..f4ea324 100644 --- a/mnefun/data/canonical.yml +++ b/mnefun/data/canonical.yml @@ -15,6 +15,7 @@ fetch_raw: run_names: runs_empty: subject_run_indices: + fix_eeg_order: scoring: score: @@ -68,6 +69,7 @@ preprocessing: fir_design: fir_window: phase: + notch_filter: ssp: auto_bad: auto_bad_reject: @@ -81,6 +83,8 @@ preprocessing: plot_raw: ssp_eog_reject: ssp_ecg_reject: + ssp_eog_baseline: + ssp_ecg_baseline: eog_channel: heog_channel: veog_channel: @@ -143,6 +147,7 @@ covariance: cov_rank_method: cov_rank_tol: force_erm_cov_rank_full: + erm_cov_from_task: forward: bem_type: From 1408e11f049d48059288bde8443456fba833bc9b Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Mon, 26 Jul 2021 15:24:29 -0400 Subject: [PATCH 3/3] FIX: More --- mnefun/_report.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/mnefun/_report.py b/mnefun/_report.py index 6e594a6..1f76da9 100644 --- a/mnefun/_report.py +++ b/mnefun/_report.py @@ -730,8 +730,10 @@ def gen_html_report(p, subjects, structurals, run_indices=None): evo.data *= SQRT_2 sl = slice(ei, n_e * n_s, n_e) these_axes = list(axes[sl]) + [axes[-1]] - evo.plot_white( - noise_cov, verbose='error', axes=these_axes) + if evo.nave > 0: + evo.plot_white( + noise_cov, axes=these_axes, + verbose='error') for ax in these_axes[:-1]: n_text = 'N=%d' % (evo.nave,) if ei != 0: @@ -748,9 +750,13 @@ def gen_html_report(p, subjects, structurals, run_indices=None): ax = axes[-1:] else: ax = axes[si * n_e:(si + 1) * n_e] - this_max = max(np.nanmax(np.abs(line.get_ydata())) - for a in ax - for line in a.lines) + lines = sum((list(a.lines) for a in ax), []) + if len(lines): + this_max = max( + np.nanmax(np.abs(line.get_ydata())) + for line in lines) + else: + this_max = np.nan this_max = 1 if np.isnan(this_max) else this_max if si == n_s: ax[0].set(ylim=[0, this_max], xlim=xlim) @@ -769,7 +775,8 @@ def gen_html_report(p, subjects, structurals, run_indices=None): else: line.set(alpha=0.5, linewidth=1) n_real += 1 - assert n_real == n_e * n_s + # not true when some conditions are empty + # assert n_real == n_e * n_s axes[-1].legend(hs, labels) if n_e > 1: axes[-1].set_title(