From 41df1f4abbb7a294fa3e1571897043c0d03fa064 Mon Sep 17 00:00:00 2001 From: CesarCaballeroGaudes Date: Mon, 11 May 2020 13:06:38 +0200 Subject: [PATCH 01/52] implementing retroicor --- phys2denoise/metrics/retroicor.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 phys2denoise/metrics/retroicor.py diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py new file mode 100644 index 0000000..df9ee3c --- /dev/null +++ b/phys2denoise/metrics/retroicor.py @@ -0,0 +1,12 @@ +"""implementation of retroicor""" + +import argparse as ap +import numpy as np +import scipy as sp +import json +import string +import random +import matplotlib as mpl; mpl.use('TkAgg') +import matplotlib.pyplot as plt +from scipy.signal import argrelmax + From 4ad18ba408801976ec0373ff813252499106b72d Mon Sep 17 00:00:00 2001 From: CesarCaballeroGaudes Date: Wed, 17 Jun 2020 13:31:14 +0200 Subject: [PATCH 02/52] draft implementation of retroicor --- phys2denoise/metrics/retroicor.py | 75 +++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py index df9ee3c..e63e929 100644 --- a/phys2denoise/metrics/retroicor.py +++ b/phys2denoise/metrics/retroicor.py @@ -10,3 +10,78 @@ import matplotlib.pyplot as plt from scipy.signal import argrelmax +""" +This function computes RETROICOR regressors (Glover et al. 2000) +""" + + +def compute_card_phase(card_peaks_timings,slice_timings,nscans,TR): + """ + + This function creates cardiac phase from cardiac peaks. + Assumes that timing of cardiac events are given in same units + as slice timings, for example seconds. + + """ + + nscans = np.shape(slice_timings) + phase_card = np.zeros(nscans) + for ii in range(nscans): + # find previous cardiac peaks + previous_card_peaks = np.asarray(np.nonzero(card_peaks_timings < slice_timings[ii])) + if np.size(previous_card_peaks) == 0: + t1 = 0 + else: + last_peak = previous_card_peaks[0][-1] + t1 = card_peaks_timings[last_peak] + + # find posterior cardiac peaks + next_card_peaks = np.asarray(np.nonzero(card_peaks_timings > slice_timings[ii])) + if np.size(next_card_peaks) == 0: + t2 = nscans * TR + else: + next_peak = next_card_peaks[0][0] + t2 = card_peaks_timings[next_peak] + + # compute cardiac phase + phase_card[ii] = (2*np.math.pi*(slice_timings[ii] - t1))/(t2-t1) + + return phase_card + +def compute_resp_phase(resp,sampling_time): + + """ + This function creates respiration phase from resp trace. + """ + + +def compute_retroicor_regressors(physio,TR,nscans,slice_timings,n_harmonics,card=FALSE,resp=FALSE): + nslices = np.shape(slice_timings) # number of slices + + # if respiration, compute histogram and temporal derivative of respiration signal + if resp: + resp_hist, resp_hist_bins = plt.hist(physio,bins=100) + resp_diff = np.diff(physio,n=1) + + #initialize output variables + retroicor_regressors = [] + phase = np.empty((nscans,nslices)) + + for jj in range(nslices): + # Initialize slice timings for current slice + crslice_timings = TR * np.arange(nscans)+slice_timings[jj] + # Compute physiological phases using the timings of physio events (e.g. peaks) slice sampling times + if card: + phase[,jj] = compute_phase_card(physio,crslice_timings) + if resp: + phase[,jj] = compute_phase_resp(resp_diff,resp_hist,resp_hist_bins,crslice_timings) + # Compute retroicor regressors + for nn in range(n_harmonics): + retroicor_regressors[jj][:,2*nn] = np.cos((nn+1)*phase[jj]) + retricor_regressor[jj][:,2*nn+1] = np.sin((nn+1)*phase[jj]) + + return retroicor_regressors,phase + + + + From ab66ba82985e05da0e3ab4cbdd41d3586a772c14 Mon Sep 17 00:00:00 2001 From: smoia Date: Wed, 11 Nov 2020 15:33:28 +0100 Subject: [PATCH 03/52] Add first skeleton draft --- phys2denoise.py | 404 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 404 insertions(+) create mode 100644 phys2denoise.py diff --git a/phys2denoise.py b/phys2denoise.py new file mode 100644 index 0000000..6bf7a47 --- /dev/null +++ b/phys2denoise.py @@ -0,0 +1,404 @@ +#!/usr/bin/env python3 + +""" +Phys2denoise is a python3 library meant to prepare physiological regressors for fMRI denoising. + +The project is under development. + +Copyright 2020, The Phys2BIDS community. +Please scroll to bottom to read full license. + +""" + +import datetime +import logging +import os +import sys +from copy import deepcopy +from shutil import copy as cp + +import numpy as np + +from phys2denoise import utils, viz, _version +from phys2denoise.cli.run import _get_parser +# from phys2denoise.metrics import cardiac, chest_belt, retroicor + +from . import __version__ +# from .due import due, Doi + +LGR = logging.getLogger(__name__) + + +def print_json(outfile, samp_freq, time_offset, ch_name): + """ + Print the json required by BIDS format. + + Parameters + ---------- + outfile: str or path + Fullpath to output file. + samp_freq: float + Frequency of sampling for the output file. + time_offset: float + Difference between beginning of file and first TR. + ch_name: list of str + List of channel names, as specified by BIDS format. + + Notes + ----- + Outcome: + outfile: .json file + File containing information for BIDS. + """ + start_time = -time_offset + summary = dict(SamplingFrequency=samp_freq, + StartTime=round(start_time, 4), + Columns=ch_name) + utils.writejson(outfile, summary, indent=4, sort_keys=False) + + +@due.dcite( + Doi('10.5281/zenodo.3470091'), + path='phys2bids', + description='Conversion of physiological trace data to BIDS format', + version=__version__, + cite_module=True) +@due.dcite( + Doi('10.1038/sdata.2016.44'), + path='phys2bids', + description='The BIDS specification', + cite_module=True) +def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None, + sub=None, ses=None, chtrig=1, chsel=None, num_timepoints_expected=None, + tr=None, thr=None, pad=9, ch_name=[], yml='', debug=False, quiet=False): + """ + Run main workflow of phys2bids. + + Runs the parser, does some checks on input, then imports + the right interface file to read the input. If only info is required, + it returns a summary onscreen. + Otherwise, it operates on the input to return a .tsv.gz file, possibly + in BIDS format. + + Raises + ------ + NotImplementedError + If the file extension is not supported yet. + """ + # Check options to make them internally coherent pt. I + # #!# This can probably be done while parsing? + outdir = utils.check_input_dir(outdir) + utils.path_exists_or_make_it(outdir) + utils.path_exists_or_make_it(os.path.join(outdir, 'code')) + conversion_path = os.path.join(outdir, 'code', 'conversion') + utils.path_exists_or_make_it(conversion_path) + + # Create logfile name + basename = 'phys2bids_' + extension = 'tsv' + isotime = datetime.datetime.now().strftime('%Y-%m-%dT%H%M%S') + logname = os.path.join(conversion_path, (basename + isotime + '.' + extension)) + + # Set logging format + log_formatter = logging.Formatter( + '%(asctime)s\t%(name)-12s\t%(levelname)-8s\t%(message)s', + datefmt='%Y-%m-%dT%H:%M:%S') + + # Set up logging file and open it for writing + log_handler = logging.FileHandler(logname) + log_handler.setFormatter(log_formatter) + sh = logging.StreamHandler() + + if quiet: + logging.basicConfig(level=logging.WARNING, + handlers=[log_handler, sh], format='%(levelname)-10s %(message)s') + elif debug: + logging.basicConfig(level=logging.DEBUG, + handlers=[log_handler, sh], format='%(levelname)-10s %(message)s') + else: + logging.basicConfig(level=logging.INFO, + handlers=[log_handler, sh], format='%(levelname)-10s %(message)s') + + version_number = _version.get_versions()['version'] + LGR.info(f'Currently running phys2bids version {version_number}') + LGR.info(f'Input file is {filename}') + + # Save call.sh + arg_str = ' '.join(sys.argv[1:]) + call_str = f'phys2bids {arg_str}' + f = open(os.path.join(conversion_path, 'call.sh'), "a") + f.write(f'#!bin/bash \n{call_str}') + f.close() + + # Check options to make them internally coherent pt. II + # #!# This can probably be done while parsing? + indir = utils.check_input_dir(indir) + if chtrig < 1: + raise Exception('Wrong trigger channel. Channel indexing starts with 1!') + + filename, ftype = utils.check_input_type(filename, + indir) + + if heur_file: + heur_file = utils.check_input_ext(heur_file, '.py') + utils.check_file_exists(heur_file) + + infile = os.path.join(indir, filename) + utils.check_file_exists(infile) + + if isinstance(num_timepoints_expected, int): + num_timepoints_expected = [num_timepoints_expected] + if isinstance(tr, (int, float)): + tr = [tr] + + if tr is not None and num_timepoints_expected is not None: + # If tr and ntp were specified, check that tr is either length one or ntp. + if len(num_timepoints_expected) != len(tr) and len(tr) != 1: + raise Exception('Number of sequence types listed with TR ' + 'doesn\'t match expected number of runs in ' + 'the session') + + # Read file! + if ftype == 'acq': + from phys2bids.interfaces.acq import populate_phys_input + elif ftype == 'txt': + from phys2bids.interfaces.txt import populate_phys_input + + LGR.info(f'Reading the file {infile}') + phys_in = populate_phys_input(infile, chtrig) + + LGR.info('Checking that units of measure are BIDS compatible') + for index, unit in enumerate(phys_in.units): + phys_in.units[index] = bids.bidsify_units(unit) + + LGR.info('Reading infos') + phys_in.print_info(filename) + # #!# Here the function viz.plot_channel should be called + viz.plot_all(phys_in.ch_name, phys_in.timeseries, phys_in.units, + phys_in.freq, infile, conversion_path) + # If only info were asked, end here. + if info: + return + + # The next few lines remove the undesired channels from phys_in. + if chsel: + LGR.info('Dropping unselected channels') + for i in reversed(range(0, phys_in.ch_amount)): + if i not in chsel: + phys_in.delete_at_index(i) + + # If requested, change channel names. + if ch_name: + LGR.info('Renaming channels with given names') + phys_in.rename_channels(ch_name) + + # Checking acquisition type via user's input + if tr is not None and num_timepoints_expected is not None: + + # Multi-run acquisition type section + # Check list length, more than 1 means multi-run + if len(num_timepoints_expected) > 1: + # if multi-run of same sequence type, pad list with ones + # and multiply array with user's input + if len(tr) == 1: + tr = np.ones(len(num_timepoints_expected)) * tr[0] + + # Sum of values in ntp_list should be equivalent to num_timepoints_found + phys_in.check_trigger_amount(thr=thr, + num_timepoints_expected=sum(num_timepoints_expected), + tr=1) + + # Check that sum of tp expected is equivalent to num_timepoints_found, + # if it passes call slice4phys + if phys_in.num_timepoints_found != sum(num_timepoints_expected): + raise Exception('The number of triggers found is different ' + 'than expected. Better stop now than break ' + 'something.') + + # slice the recording based on user's entries + # !!! ATTENTION: PHYS_IN GETS OVERWRITTEN AS DICTIONARY + phys_in = slice4phys(phys_in, num_timepoints_expected, tr, + phys_in.thr, pad) + # returns a dictionary in the form {run_idx: phys_in[startpoint, endpoint]} + + # save a figure for each run | give the right acquisition parameters for runs + fileprefix = os.path.join(conversion_path, + os.path.splitext(os.path.basename(filename))[0]) + for i, run in enumerate(phys_in.keys()): + plot_fileprefix = f'{fileprefix}_{run}' + viz.export_trigger_plot(phys_in[run], chtrig, plot_fileprefix, tr[i], + num_timepoints_expected[i], filename, + sub, ses) + + # Single run acquisition type, or : nothing to split workflow + else: + # Run analysis on trigger channel to get first timepoint + # and the time offset. + phys_in.check_trigger_amount(thr, num_timepoints_expected[0], tr[0]) + # save a figure of the trigger + fileprefix = os.path.join(conversion_path, + os.path.splitext(os.path.basename(filename))[0]) + viz.export_trigger_plot(phys_in, chtrig, fileprefix, tr[0], + num_timepoints_expected[0], filename, + sub, ses) + + # Reassign phys_in as dictionary + # !!! ATTENTION: PHYS_IN GETS OVERWRITTEN AS DICTIONARY + phys_in = {1: phys_in} + + else: + LGR.warning('Skipping trigger pulse count. If you want to run it, ' + 'call phys2bids using both "-ntp" and "-tr" arguments') + # !!! ATTENTION: PHYS_IN GETS OVERWRITTEN AS DICTIONARY + phys_in = {1: phys_in} + + # The next few lines create a dictionary of different BlueprintInput + # objects, one for each unique frequency for each run in phys_in + # they also save the amount of runs and unique frequencies + run_amount = len(phys_in) + uniq_freq_list = set(phys_in[1].freq) + freq_amount = len(uniq_freq_list) + if freq_amount > 1: + LGR.info(f'Found {freq_amount} different frequencies in input!') + + if run_amount > 1: + LGR.info(f'Found {run_amount} different scans in input!') + + LGR.info(f'Preparing {freq_amount*run_amount} output files.') + # Create phys_out dict that will have a blueprint object for each different frequency + phys_out = {} + + if heur_file is not None and sub is not None: + LGR.info(f'Preparing BIDS output using {heur_file}') + # If heuristics are used, init a dict of arguments to pass to use_heuristic + heur_args = {'heur_file': heur_file, 'sub': sub, 'ses': ses, + 'filename': filename, 'outdir': outdir, 'run': '', + 'record_label': ''} + # Generate participants.tsv file if it doesn't exist already. + # Update the file if the subject is not in the file. + # Do not update if the subject is already in the file. + bids.participants_file(outdir, yml, sub) + # Generate dataset_description.json file if it doesn't exist already. + bids.dataset_description_file(outdir) + # Generate README file if it doesn't exist already. + bids.readme_file(outdir) + cp(heur_file, os.path.join(conversion_path, + os.path.splitext(os.path.basename(heur_file))[0] + '.py')) + elif heur_file is not None and sub is None: + LGR.warning('While "-heur" was specified, option "-sub" was not.\n' + 'Skipping BIDS formatting.') + + # Export a (set of) phys_out for each element in phys_in + # run keys start from 1 (human friendly) + for run in phys_in.keys(): + for uniq_freq in uniq_freq_list: + # Initialise the key for the (possibly huge amount of) dictionary entries + key = f'{run}_{uniq_freq}' + # copy the phys_in object to the new dict entry + phys_out[key] = deepcopy(phys_in[run]) + # this counter will take into account how many channels are eliminated + count = 0 + # for each channel in the original phys_in object + # take the frequency + for idx, i in enumerate(phys_in[run].freq): + # if that frequency is different than the frequency of the phys_obj entry + if i != uniq_freq: + # eliminate that channel from the dict since we only want channels + # with the same frequency + phys_out[key].delete_at_index(idx - count) + # take into acount the elimination so in the next eliminated channel we + # eliminate correctly + count += 1 + # Also create a BlueprintOutput object for each unique frequency found. + # Populate it with the corresponding blueprint input and replace it + # in the dictionary. + # Add time channel in the proper frequency. + if uniq_freq != phys_in[run].freq[0]: + phys_out[key].ch_name.insert(0, phys_in[run].ch_name[0]) + phys_out[key].units.insert(0, phys_in[run].units[0]) + phys_out[key].timeseries.insert(0, np.linspace(phys_in[run].timeseries[0][0], + phys_in[run].timeseries[0][-1], + num=phys_out[key].timeseries[0].shape[0])) + # Add trigger channel in the proper frequency. + if uniq_freq != phys_in[run].freq[chtrig]: + phys_out[key].ch_name.insert(1, phys_in[run].ch_name[chtrig]) + phys_out[key].units.insert(1, phys_in[run].units[chtrig]) + phys_out[key].timeseries.insert(1, np.interp(phys_out[key].timeseries[0], + phys_in[run].timeseries[0], + phys_in[run].timeseries[chtrig])) + phys_out[key] = BlueprintOutput.init_from_blueprint(phys_out[key]) + + # Preparing output parameters: name and folder. + for uniq_freq in uniq_freq_list: + key = f'{run}_{uniq_freq}' + # If possible, prepare bids renaming. + if heur_file is not None and sub is not None: + # Add run info to heur_args if more than one run is present + if run_amount > 1: + heur_args['run'] = f'{run:02d}' + + # Append "recording-freq" to filename if more than one freq + if freq_amount > 1: + heur_args['record_label'] = f'{uniq_freq:.0f}Hz' + + phys_out[key].filename = bids.use_heuristic(**heur_args) + + # If any filename exists already because of multirun, append labels + # But warn about the non-validity of this BIDS-like name. + if run_amount > 1: + if any([phys.filename == phys_out[key].filename + for phys in phys_out.values()]): + phys_out[key].filename = (f'{phys_out[key].filename}' + f'_take-{run}') + LGR.warning('Identified multiple outputs with the same name.\n' + 'Appending fake label to avoid overwriting.\n' + '!!! ATTENTION !!! the output is not BIDS compliant.\n' + 'Please check heuristics to solve the problem.') + + else: + phys_out[key].filename = os.path.join(outdir, + os.path.splitext(os.path.basename(filename) + )[0]) + # Append "run" to filename if more than one run + if run_amount > 1: + phys_out[key].filename = f'{phys_out[key].filename}_{run:02d}' + # Append "freq" to filename if more than one freq + if freq_amount > 1: + phys_out[key].filename = f'{phys_out[key].filename}_{uniq_freq:.0f}Hz' + + LGR.info(f'Exporting files for run {run} freq {uniq_freq}') + np.savetxt(phys_out[key].filename + '.tsv.gz', + phys_out[key].timeseries, fmt='%.8e', delimiter='\t') + print_json(phys_out[key].filename, phys_out[key].freq, + phys_out[key].start_time, phys_out[key].ch_name) + print_summary(filename, num_timepoints_expected, + phys_in[run].num_timepoints_found, uniq_freq, + phys_out[key].start_time, + os.path.join(conversion_path, + os.path.splitext(os.path.basename(phys_out[key].filename) + )[0])) + + +def _main(argv=None): + options = _get_parser().parse_args(argv) + phys2bids(**vars(options)) + + +if __name__ == '__main__': + _main(sys.argv[1:]) + +""" +Copyright 2019, The Phys2BIDS community. + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + +http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +""" From 5f685eb7f73a7b62c976fd99c5e4c89738c5e5fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?In=C3=A9s=20Chavarr=C3=ADa?= <72545702+ineschh@users.noreply.github.com> Date: Thu, 12 Nov 2020 15:53:34 +0100 Subject: [PATCH 04/52] add parser --- phys2denoise/cli/run.py | 196 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 196 insertions(+) create mode 100644 phys2denoise/cli/run.py diff --git a/phys2denoise/cli/run.py b/phys2denoise/cli/run.py new file mode 100644 index 0000000..7f79b40 --- /dev/null +++ b/phys2denoise/cli/run.py @@ -0,0 +1,196 @@ +# -*- coding: utf-8 -*- +"""Parser for phys2denoise.""" + + +import argparse + +from phys2denoise import __version__ + + +def _get_parser(): + """ + Parse command line inputs for this function. + + Returns + ------- + parser.parse_args() : argparse dict + + Notes + ----- + # Argument parser follow template provided by RalphyZ. + # https://stackoverflow.com/a/43456577 + """ + parser = argparse.ArgumentParser() + optional = parser._action_groups.pop() + metric = parser._action_groups.pop() + required = parser.add_argument_group('Required Argument:') + required.add_argument('-in', '--input-file', + dest='filename', + type=str, + help='Path/name of the file containing physiological ' + 'data, with or without extension.', + required=True) + metric.add_argument('-crf', '--crf', + dest='metric_list', + action='append_const', + const='crf', + help='Cardiac response function. Needs the following ' + 'inputs:sr, os, tl, onset and tr.', + default=False) + metric.add_argument('-rpv', '--rpv', + dest='metric_list', + action='append_const', + const='rpv', + help='REspiratory pattern variability. Needs the following ' + 'inputs: bts and win.', + default=False) + metric.add_argument('-env', '--env', + dest='metric_list', + action='append_const', + const='env', + help='Respiratory pattern variability calculated across a sliding ' + 'window. Needs the following inputs: bts, sr, osr, win and lags.', + default=False) + metric.add_argument('-rv', '--rv', + dest='metric_list', + action='append_const', + const='rv', + help='Respiratory variance. Needs the following inputs: ' + 'bts, sr, osr, win and lags.', + default=False) + metric.add_argument('-rvt', '--rvt', + dest='metric_list', + action='append_const', + const='rvt', + help='Respiratory volume-per-time. Needs the following inputs: ' + 'bts, sr, osr, win and lags.', + default=False) + metric.add_argument('-rrf', '--rrf', + dest='metric_list', + action='append_const', + const='rrf', + help='Respiratory response function. Needs the following inputs: ' + 'sr, os, tl, onset and tr.', + default=False) + metric.add_argument('-rcard', '--retroicor-card', + dest='metric_list', + action='append_const', + const='r_card', + help='Computes regressors for cardiac signal. Needs the following ' + 'inputs: tr, nscans, slt and n_harm.', + default=False) + metric.add_argument('-rresp', '--retroicor-resp', + dest='metric_list', + action='append_const', + const='r_resp', + help='Computes regressors for respiratory signal. Needs the following ' + 'inputs: tr, nscans, slt and n_harm.', + default=False) + optional.add_argument('-outdir', '--output-dir', + dest='outdir', + type=str, + help='Folder where output should be placed. ' + 'Default is current folder.', + default='.') + optional.add_argument('-sr', '--sample-rate', + dest='sample_rate', + type=float, + help='Sampling rate of the physiological data in Hz.', + default=None) + optional.add_argument('-pk', '--peaks', + dest='peaks', + type=str, + help='Filename of the list with the indexed peaks\' positions' + ' of the physiological data.', + default=None) + optional.add_argument('-thr', '--throughts', + dest='throughts', + type=str, + help='Filename of the list with the indexed peaks\' positions' + ' of the physiological data.', + default=None) + optional.add_argument('-os', '--oversampling', + dest='oversampling', + type=int, + help='Temporal oversampling factor in seconds. ' + 'Default is 50.', + default=50) + optional.add_argument('-tl', '--time-length', + dest='time_length', + type=int, + help='RRF Kernel length in seconds.', + default=None) + optional.add_argument('-onset', '--onset', + dest='onset', + type=float, + help='Onset of the response in seconds. ' + 'Default is 0.', + default=0) + optional.add_argument('-tr', '--tr', + dest='tr', + type=float, + help='TR of sequence in seconds.', + default=None) + optional.add_argument('-bts', '--belt-ts', + dest='belt_ts', + type=str, + help='Filename of the 1D array containing the .' + 'respiratory belt time series.', + default=None) + optional.add_argument('-win', '--window', + dest='window', + type=int, + help='Size of the sliding window in seconds. ' + 'Default is 6 seconds.', + default=6) + optional.add_argument('-osr', '--out-samplerate', + dest='out_samplerate', + type=float, + help='Sampling rate for the output time series ' + 'in seconds. Corresponds to TR in fMRI data.', + default=None) + optional.add_argument('-lags', '--lags', + dest='lags', + nargs='*', + type=int, + action='append', + help='List of lags to apply to the rv estimate ' + 'in seconds.', + default=None) + optional.add_argument('-nscans', '--nscans', + dest='nscans', + type=int, + help='Number of scans. Default is 1.', + default=1) + optional.add_argument('-slt', '--slice-timings', + dest='slice_timings', + type=str, + help='Filename with the slice timings.', + default=None) + optional.add_argument('-nharm', '--number-harmonics', + dest='n_harm', + type=int, + help='Number of harmonics. ', + default=None) + optional.add_argument('-debug', '--debug', + dest='debug', + action='store_true', + help='Only print debugging info to log file. Default is False.', + default=False) + optional.add_argument('-quiet', '--quiet', + dest='quiet', + action='store_true', + help='Only print warnings to log file. Default is False.', + default=False) + optional.add_argument('-v', '--version', action='version', + version=('%(prog)s ' + __version__)) + + parser._action_groups.append(optional) + + return parser + + +if __name__ == '__main__': + raise RuntimeError('phys2denoise/cli/run.py should not be run directly;\n' + 'Please `pip install` phys2denoise and use the ' + '`phys2denoise` command') From 2fa576faf3afee16d97be66ed1ad987028a15238 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?In=C3=A9s=20Chavarr=C3=ADa?= <72545702+ineschh@users.noreply.github.com> Date: Thu, 12 Nov 2020 19:09:48 +0100 Subject: [PATCH 05/52] Update phys2denoise/cli/run.py Co-authored-by: Stefano Moia --- phys2denoise/cli/run.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phys2denoise/cli/run.py b/phys2denoise/cli/run.py index 7f79b40..d3f3384 100644 --- a/phys2denoise/cli/run.py +++ b/phys2denoise/cli/run.py @@ -41,7 +41,7 @@ def _get_parser(): dest='metric_list', action='append_const', const='rpv', - help='REspiratory pattern variability. Needs the following ' + help='Respiratory pattern variability. Needs the following ' 'inputs: bts and win.', default=False) metric.add_argument('-env', '--env', From 6142d30787c63fc5d9a090ef6beb34cc6d93ab49 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?In=C3=A9s=20Chavarr=C3=ADa?= <72545702+ineschh@users.noreply.github.com> Date: Thu, 12 Nov 2020 20:31:26 +0100 Subject: [PATCH 06/52] Update --- phys2denoise/cli/run.py | 245 +++++++++++++++++++--------------------- 1 file changed, 115 insertions(+), 130 deletions(-) diff --git a/phys2denoise/cli/run.py b/phys2denoise/cli/run.py index 7f79b40..bfee4ac 100644 --- a/phys2denoise/cli/run.py +++ b/phys2denoise/cli/run.py @@ -23,174 +23,159 @@ def _get_parser(): parser = argparse.ArgumentParser() optional = parser._action_groups.pop() metric = parser._action_groups.pop() - required = parser.add_argument_group('Required Argument:') - required.add_argument('-in', '--input-file', - dest='filename', + required = parser.add_argument_group("Required Argument:") + required.add_argument("-in", "--input-file", + dest="filename", type=str, - help='Path/name of the file containing physiological ' - 'data, with or without extension.', + help="Full path and name of the file containing " + "physiological data, with or without extension.", required=True) - metric.add_argument('-crf', '--crf', - dest='metric_list', - action='append_const', - const='crf', - help='Cardiac response function. Needs the following ' - 'inputs:sr, os, tl, onset and tr.', + metric.add_argument("-crf", "--cardiac-response-function", + dest="metrics", + action="append_const", + const="crf", + help="Cardiac response function. Needs the following " + "inputs:sample-rate, oversampling, time-length, " + "onset and tr.", default=False) - metric.add_argument('-rpv', '--rpv', - dest='metric_list', - action='append_const', - const='rpv', - help='REspiratory pattern variability. Needs the following ' - 'inputs: bts and win.', + metric.add_argument("-rpv", "--respiratory-pattern-variability", + dest="metrics", + action="append_const", + const="rpv", + help="Respiratory pattern variability. Needs the following " + "input: window.", default=False) - metric.add_argument('-env', '--env', - dest='metric_list', - action='append_const', - const='env', - help='Respiratory pattern variability calculated across a sliding ' - 'window. Needs the following inputs: bts, sr, osr, win and lags.', + metric.add_argument("-env", "--envelope", + dest="metrics", + action="append_const", + const="env", + help="Respiratory pattern variability calculated across a sliding " + "window. Needs the following inputs: sample-rate, window and lags.", default=False) - metric.add_argument('-rv', '--rv', - dest='metric_list', - action='append_const', - const='rv', - help='Respiratory variance. Needs the following inputs: ' - 'bts, sr, osr, win and lags.', + metric.add_argument("-rv", "--respiratory-variance", + dest="metrics", + action="append_const", + const="rv", + help="Respiratory variance. Needs the following inputs: " + "sample-rate, window and lags.", default=False) - metric.add_argument('-rvt', '--rvt', - dest='metric_list', - action='append_const', - const='rvt', - help='Respiratory volume-per-time. Needs the following inputs: ' - 'bts, sr, osr, win and lags.', + metric.add_argument("-rvt", "--respiratory-volume-per-time", + dest="metrics", + action="append_const", + const="rvt", + help="Respiratory volume-per-time. Needs the following inputs: " + "sample-rate, window and lags.", default=False) - metric.add_argument('-rrf', '--rrf', - dest='metric_list', - action='append_const', - const='rrf', - help='Respiratory response function. Needs the following inputs: ' - 'sr, os, tl, onset and tr.', + metric.add_argument("-rrf", "--respiratory-response-function", + dest="metrics", + action="append_const", + const="rrf", + help="Respiratory response function. Needs the following inputs: " + "sample-rate, oversampling, time-length, onset and tr.", default=False) - metric.add_argument('-rcard', '--retroicor-card', - dest='metric_list', - action='append_const', - const='r_card', - help='Computes regressors for cardiac signal. Needs the following ' - 'inputs: tr, nscans, slt and n_harm.', + metric.add_argument("-rcard", "--retroicor-card", + dest="metrics", + action="append_const", + const="r_card", + help="Computes regressors for cardiac signal. Needs the following " + "inputs: tr, nscans and n_harm.", default=False) - metric.add_argument('-rresp', '--retroicor-resp', - dest='metric_list', - action='append_const', - const='r_resp', - help='Computes regressors for respiratory signal. Needs the following ' - 'inputs: tr, nscans, slt and n_harm.', + metric.add_argument("-rresp", "--retroicor-resp", + dest="metrics", + action="append_const", + const="r_resp", + help="Computes regressors for respiratory signal. Needs the following " + "inputs: tr, nscans and n_harm.", default=False) - optional.add_argument('-outdir', '--output-dir', - dest='outdir', + optional.add_argument("-outdir", "--output-dir", + dest="outdir", type=str, - help='Folder where output should be placed. ' - 'Default is current folder.', - default='.') - optional.add_argument('-sr', '--sample-rate', - dest='sample_rate', + help="Folder where output should be placed. " + "Default is current folder.", + default=".") + optional.add_argument("-sr", "--sample-rate", + dest="sample_rate", type=float, - help='Sampling rate of the physiological data in Hz.', + help="Sampling rate of the physiological data in Hz.", default=None) - optional.add_argument('-pk', '--peaks', - dest='peaks', + optional.add_argument("-pk", "--peaks", + dest="peaks", type=str, - help='Filename of the list with the indexed peaks\' positions' - ' of the physiological data.', + help="Full path and filename of the list with the indexed peaks' " + "positions of the physiological data.", default=None) - optional.add_argument('-thr', '--throughts', - dest='throughts', + optional.add_argument("-tg", "--throughts", + dest="throughts", type=str, - help='Filename of the list with the indexed peaks\' positions' - ' of the physiological data.', + help="Full path and filename of the list with the indexed peaks' " + "positions of the physiological data.", default=None) - optional.add_argument('-os', '--oversampling', - dest='oversampling', + optional.add_argument("-os", "--oversampling", + dest="oversampling", type=int, - help='Temporal oversampling factor in seconds. ' - 'Default is 50.', + help="Temporal oversampling factor in seconds. " + "Default is 50.", default=50) - optional.add_argument('-tl', '--time-length', - dest='time_length', + optional.add_argument("-tl", "--time-length", + dest="time_length", type=int, - help='RRF Kernel length in seconds.', + help="RRF Kernel length in seconds.", default=None) - optional.add_argument('-onset', '--onset', - dest='onset', + optional.add_argument("-onset", "--onset", + dest="onset", type=float, - help='Onset of the response in seconds. ' - 'Default is 0.', + help="Onset of the response in seconds. " + "Default is 0.", default=0) - optional.add_argument('-tr', '--tr', - dest='tr', + optional.add_argument("-tr", "--tr", + dest="tr", type=float, - help='TR of sequence in seconds.', + help="TR of sequence in seconds.", default=None) - optional.add_argument('-bts', '--belt-ts', - dest='belt_ts', - type=str, - help='Filename of the 1D array containing the .' - 'respiratory belt time series.', - default=None) - optional.add_argument('-win', '--window', - dest='window', + optional.add_argument("-win", "--window", + dest="window", type=int, - help='Size of the sliding window in seconds. ' - 'Default is 6 seconds.', + help="Size of the sliding window in seconds. " + "Default is 6 seconds.", default=6) - optional.add_argument('-osr', '--out-samplerate', - dest='out_samplerate', - type=float, - help='Sampling rate for the output time series ' - 'in seconds. Corresponds to TR in fMRI data.', - default=None) - optional.add_argument('-lags', '--lags', - dest='lags', - nargs='*', + optional.add_argument("-lags", "--lags", + dest="lags", + nargs="*", type=int, - action='append', - help='List of lags to apply to the rv estimate ' - 'in seconds.', + action="append", + help="List of lags to apply to the RV estimate " + "in seconds.", default=None) - optional.add_argument('-nscans', '--nscans', - dest='nscans', + optional.add_argument("-nscans", "--number-scans", + dest="nscans", type=int, - help='Number of scans. Default is 1.', + help="Number of scans. Default is 1.", default=1) - optional.add_argument('-slt', '--slice-timings', - dest='slice_timings', - type=str, - help='Filename with the slice timings.', - default=None) - optional.add_argument('-nharm', '--number-harmonics', - dest='n_harm', + optional.add_argument("-nharm", "--number-harmonics", + dest="n_harm", type=int, - help='Number of harmonics. ', + help="Number of harmonics.", default=None) - optional.add_argument('-debug', '--debug', - dest='debug', - action='store_true', - help='Only print debugging info to log file. Default is False.', + optional.add_argument("-debug", "--debug", + dest="debug", + action="store_true", + help="Only print debugging info to log file. Default is False.", default=False) - optional.add_argument('-quiet', '--quiet', - dest='quiet', - action='store_true', - help='Only print warnings to log file. Default is False.', + optional.add_argument("-quiet", "--quiet", + dest="quiet", + action="store_true", + help="Only print warnings to log file. Default is False.", default=False) - optional.add_argument('-v', '--version', action='version', - version=('%(prog)s ' + __version__)) + optional.add_argument("-v", "--version", action="version", + version=("%(prog)s " + __version__)) parser._action_groups.append(optional) + parser._action_groups.append(metric) return parser -if __name__ == '__main__': - raise RuntimeError('phys2denoise/cli/run.py should not be run directly;\n' - 'Please `pip install` phys2denoise and use the ' - '`phys2denoise` command') +if __name__ == "__main__": + raise RuntimeError("phys2denoise/cli/run.py should not be run directly;\n" + "Please `pip install` phys2denoise and use the " + "`phys2denoise` command") From 8e521f3208ba8eea3b73a478b1c19e52e0f952d1 Mon Sep 17 00:00:00 2001 From: smoia Date: Fri, 13 Nov 2020 00:51:33 +0100 Subject: [PATCH 07/52] More skeleton --- phys2denoise.py | 314 ++++++------------------------------------------ 1 file changed, 36 insertions(+), 278 deletions(-) diff --git a/phys2denoise.py b/phys2denoise.py index 6bf7a47..17a6ed5 100644 --- a/phys2denoise.py +++ b/phys2denoise.py @@ -5,7 +5,7 @@ The project is under development. -Copyright 2020, The Phys2BIDS community. +Copyright 2020, The physiopy community. Please scroll to bottom to read full license. """ @@ -24,7 +24,7 @@ # from phys2denoise.metrics import cardiac, chest_belt, retroicor from . import __version__ -# from .due import due, Doi +from .due import due, Doi LGR = logging.getLogger(__name__) @@ -58,43 +58,32 @@ def print_json(outfile, samp_freq, time_offset, ch_name): @due.dcite( - Doi('10.5281/zenodo.3470091'), - path='phys2bids', - description='Conversion of physiological trace data to BIDS format', + Doi(''), + path='phys2denoise', + description='Creation of regressors for physiological denoising', version=__version__, cite_module=True) -@due.dcite( - Doi('10.1038/sdata.2016.44'), - path='phys2bids', - description='The BIDS specification', - cite_module=True) -def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None, - sub=None, ses=None, chtrig=1, chsel=None, num_timepoints_expected=None, - tr=None, thr=None, pad=9, ch_name=[], yml='', debug=False, quiet=False): +def phys2denoise(filename, outdir='.', debug=False, quiet=False): """ - Run main workflow of phys2bids. - - Runs the parser, does some checks on input, then imports - the right interface file to read the input. If only info is required, - it returns a summary onscreen. - Otherwise, it operates on the input to return a .tsv.gz file, possibly - in BIDS format. - - Raises - ------ - NotImplementedError - If the file extension is not supported yet. + Run main workflow of phys2denoise. + + Runs the parser, does some checks on input, then computes the required metrics. + + Notes + ----- + The code was greatly copied from phys2bids (copyright the physiopy community) + """ # Check options to make them internally coherent pt. I # #!# This can probably be done while parsing? - outdir = utils.check_input_dir(outdir) - utils.path_exists_or_make_it(outdir) - utils.path_exists_or_make_it(os.path.join(outdir, 'code')) + outdir = os.path.abspath(outdir) + os.makedirs(outdir) + os.makedirs(os.path.join(outdir, 'code')) conversion_path = os.path.join(outdir, 'code', 'conversion') - utils.path_exists_or_make_it(conversion_path) + os.makedirs(conversion_path) # Create logfile name - basename = 'phys2bids_' + basename = 'phys2denoise_' extension = 'tsv' isotime = datetime.datetime.now().strftime('%Y-%m-%dT%H%M%S') logname = os.path.join(conversion_path, (basename + isotime + '.' + extension)) @@ -120,275 +109,44 @@ def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None, handlers=[log_handler, sh], format='%(levelname)-10s %(message)s') version_number = _version.get_versions()['version'] - LGR.info(f'Currently running phys2bids version {version_number}') + LGR.info(f'Currently running phys2denoise version {version_number}') LGR.info(f'Input file is {filename}') # Save call.sh arg_str = ' '.join(sys.argv[1:]) - call_str = f'phys2bids {arg_str}' + call_str = f'phys2denoise {arg_str}' f = open(os.path.join(conversion_path, 'call.sh'), "a") f.write(f'#!bin/bash \n{call_str}') f.close() # Check options to make them internally coherent pt. II # #!# This can probably be done while parsing? - indir = utils.check_input_dir(indir) - if chtrig < 1: - raise Exception('Wrong trigger channel. Channel indexing starts with 1!') - - filename, ftype = utils.check_input_type(filename, - indir) - - if heur_file: - heur_file = utils.check_input_ext(heur_file, '.py') - utils.check_file_exists(heur_file) - - infile = os.path.join(indir, filename) - utils.check_file_exists(infile) - - if isinstance(num_timepoints_expected, int): - num_timepoints_expected = [num_timepoints_expected] - if isinstance(tr, (int, float)): - tr = [tr] - - if tr is not None and num_timepoints_expected is not None: - # If tr and ntp were specified, check that tr is either length one or ntp. - if len(num_timepoints_expected) != len(tr) and len(tr) != 1: - raise Exception('Number of sequence types listed with TR ' - 'doesn\'t match expected number of runs in ' - 'the session') - - # Read file! - if ftype == 'acq': - from phys2bids.interfaces.acq import populate_phys_input - elif ftype == 'txt': - from phys2bids.interfaces.txt import populate_phys_input - - LGR.info(f'Reading the file {infile}') - phys_in = populate_phys_input(infile, chtrig) - - LGR.info('Checking that units of measure are BIDS compatible') - for index, unit in enumerate(phys_in.units): - phys_in.units[index] = bids.bidsify_units(unit) - - LGR.info('Reading infos') - phys_in.print_info(filename) - # #!# Here the function viz.plot_channel should be called - viz.plot_all(phys_in.ch_name, phys_in.timeseries, phys_in.units, - phys_in.freq, infile, conversion_path) - # If only info were asked, end here. - if info: - return - - # The next few lines remove the undesired channels from phys_in. - if chsel: - LGR.info('Dropping unselected channels') - for i in reversed(range(0, phys_in.ch_amount)): - if i not in chsel: - phys_in.delete_at_index(i) - - # If requested, change channel names. - if ch_name: - LGR.info('Renaming channels with given names') - phys_in.rename_channels(ch_name) - - # Checking acquisition type via user's input - if tr is not None and num_timepoints_expected is not None: - - # Multi-run acquisition type section - # Check list length, more than 1 means multi-run - if len(num_timepoints_expected) > 1: - # if multi-run of same sequence type, pad list with ones - # and multiply array with user's input - if len(tr) == 1: - tr = np.ones(len(num_timepoints_expected)) * tr[0] - - # Sum of values in ntp_list should be equivalent to num_timepoints_found - phys_in.check_trigger_amount(thr=thr, - num_timepoints_expected=sum(num_timepoints_expected), - tr=1) - - # Check that sum of tp expected is equivalent to num_timepoints_found, - # if it passes call slice4phys - if phys_in.num_timepoints_found != sum(num_timepoints_expected): - raise Exception('The number of triggers found is different ' - 'than expected. Better stop now than break ' - 'something.') - - # slice the recording based on user's entries - # !!! ATTENTION: PHYS_IN GETS OVERWRITTEN AS DICTIONARY - phys_in = slice4phys(phys_in, num_timepoints_expected, tr, - phys_in.thr, pad) - # returns a dictionary in the form {run_idx: phys_in[startpoint, endpoint]} - - # save a figure for each run | give the right acquisition parameters for runs - fileprefix = os.path.join(conversion_path, - os.path.splitext(os.path.basename(filename))[0]) - for i, run in enumerate(phys_in.keys()): - plot_fileprefix = f'{fileprefix}_{run}' - viz.export_trigger_plot(phys_in[run], chtrig, plot_fileprefix, tr[i], - num_timepoints_expected[i], filename, - sub, ses) - - # Single run acquisition type, or : nothing to split workflow - else: - # Run analysis on trigger channel to get first timepoint - # and the time offset. - phys_in.check_trigger_amount(thr, num_timepoints_expected[0], tr[0]) - # save a figure of the trigger - fileprefix = os.path.join(conversion_path, - os.path.splitext(os.path.basename(filename))[0]) - viz.export_trigger_plot(phys_in, chtrig, fileprefix, tr[0], - num_timepoints_expected[0], filename, - sub, ses) - - # Reassign phys_in as dictionary - # !!! ATTENTION: PHYS_IN GETS OVERWRITTEN AS DICTIONARY - phys_in = {1: phys_in} + # filename, ftype = utils.check_input_type(filename) + + if not os.path.isfile(filename) and filename is not None: + raise FileNotFoundError(f'The file {filename} does not exist!') + + # Read input file + phys_in = np.genfromtxt(filename) + + + + + + - else: - LGR.warning('Skipping trigger pulse count. If you want to run it, ' - 'call phys2bids using both "-ntp" and "-tr" arguments') - # !!! ATTENTION: PHYS_IN GETS OVERWRITTEN AS DICTIONARY - phys_in = {1: phys_in} - - # The next few lines create a dictionary of different BlueprintInput - # objects, one for each unique frequency for each run in phys_in - # they also save the amount of runs and unique frequencies - run_amount = len(phys_in) - uniq_freq_list = set(phys_in[1].freq) - freq_amount = len(uniq_freq_list) - if freq_amount > 1: - LGR.info(f'Found {freq_amount} different frequencies in input!') - - if run_amount > 1: - LGR.info(f'Found {run_amount} different scans in input!') - - LGR.info(f'Preparing {freq_amount*run_amount} output files.') - # Create phys_out dict that will have a blueprint object for each different frequency - phys_out = {} - - if heur_file is not None and sub is not None: - LGR.info(f'Preparing BIDS output using {heur_file}') - # If heuristics are used, init a dict of arguments to pass to use_heuristic - heur_args = {'heur_file': heur_file, 'sub': sub, 'ses': ses, - 'filename': filename, 'outdir': outdir, 'run': '', - 'record_label': ''} - # Generate participants.tsv file if it doesn't exist already. - # Update the file if the subject is not in the file. - # Do not update if the subject is already in the file. - bids.participants_file(outdir, yml, sub) - # Generate dataset_description.json file if it doesn't exist already. - bids.dataset_description_file(outdir) - # Generate README file if it doesn't exist already. - bids.readme_file(outdir) - cp(heur_file, os.path.join(conversion_path, - os.path.splitext(os.path.basename(heur_file))[0] + '.py')) - elif heur_file is not None and sub is None: - LGR.warning('While "-heur" was specified, option "-sub" was not.\n' - 'Skipping BIDS formatting.') - - # Export a (set of) phys_out for each element in phys_in - # run keys start from 1 (human friendly) - for run in phys_in.keys(): - for uniq_freq in uniq_freq_list: - # Initialise the key for the (possibly huge amount of) dictionary entries - key = f'{run}_{uniq_freq}' - # copy the phys_in object to the new dict entry - phys_out[key] = deepcopy(phys_in[run]) - # this counter will take into account how many channels are eliminated - count = 0 - # for each channel in the original phys_in object - # take the frequency - for idx, i in enumerate(phys_in[run].freq): - # if that frequency is different than the frequency of the phys_obj entry - if i != uniq_freq: - # eliminate that channel from the dict since we only want channels - # with the same frequency - phys_out[key].delete_at_index(idx - count) - # take into acount the elimination so in the next eliminated channel we - # eliminate correctly - count += 1 - # Also create a BlueprintOutput object for each unique frequency found. - # Populate it with the corresponding blueprint input and replace it - # in the dictionary. - # Add time channel in the proper frequency. - if uniq_freq != phys_in[run].freq[0]: - phys_out[key].ch_name.insert(0, phys_in[run].ch_name[0]) - phys_out[key].units.insert(0, phys_in[run].units[0]) - phys_out[key].timeseries.insert(0, np.linspace(phys_in[run].timeseries[0][0], - phys_in[run].timeseries[0][-1], - num=phys_out[key].timeseries[0].shape[0])) - # Add trigger channel in the proper frequency. - if uniq_freq != phys_in[run].freq[chtrig]: - phys_out[key].ch_name.insert(1, phys_in[run].ch_name[chtrig]) - phys_out[key].units.insert(1, phys_in[run].units[chtrig]) - phys_out[key].timeseries.insert(1, np.interp(phys_out[key].timeseries[0], - phys_in[run].timeseries[0], - phys_in[run].timeseries[chtrig])) - phys_out[key] = BlueprintOutput.init_from_blueprint(phys_out[key]) - - # Preparing output parameters: name and folder. - for uniq_freq in uniq_freq_list: - key = f'{run}_{uniq_freq}' - # If possible, prepare bids renaming. - if heur_file is not None and sub is not None: - # Add run info to heur_args if more than one run is present - if run_amount > 1: - heur_args['run'] = f'{run:02d}' - - # Append "recording-freq" to filename if more than one freq - if freq_amount > 1: - heur_args['record_label'] = f'{uniq_freq:.0f}Hz' - - phys_out[key].filename = bids.use_heuristic(**heur_args) - - # If any filename exists already because of multirun, append labels - # But warn about the non-validity of this BIDS-like name. - if run_amount > 1: - if any([phys.filename == phys_out[key].filename - for phys in phys_out.values()]): - phys_out[key].filename = (f'{phys_out[key].filename}' - f'_take-{run}') - LGR.warning('Identified multiple outputs with the same name.\n' - 'Appending fake label to avoid overwriting.\n' - '!!! ATTENTION !!! the output is not BIDS compliant.\n' - 'Please check heuristics to solve the problem.') - - else: - phys_out[key].filename = os.path.join(outdir, - os.path.splitext(os.path.basename(filename) - )[0]) - # Append "run" to filename if more than one run - if run_amount > 1: - phys_out[key].filename = f'{phys_out[key].filename}_{run:02d}' - # Append "freq" to filename if more than one freq - if freq_amount > 1: - phys_out[key].filename = f'{phys_out[key].filename}_{uniq_freq:.0f}Hz' - - LGR.info(f'Exporting files for run {run} freq {uniq_freq}') - np.savetxt(phys_out[key].filename + '.tsv.gz', - phys_out[key].timeseries, fmt='%.8e', delimiter='\t') - print_json(phys_out[key].filename, phys_out[key].freq, - phys_out[key].start_time, phys_out[key].ch_name) - print_summary(filename, num_timepoints_expected, - phys_in[run].num_timepoints_found, uniq_freq, - phys_out[key].start_time, - os.path.join(conversion_path, - os.path.splitext(os.path.basename(phys_out[key].filename) - )[0])) def _main(argv=None): options = _get_parser().parse_args(argv) - phys2bids(**vars(options)) + phys2denoise(**vars(options)) if __name__ == '__main__': _main(sys.argv[1:]) """ -Copyright 2019, The Phys2BIDS community. +Copyright 2019, The phys2denoise community. Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. From 4c265e3cca715104b9f155a0e58c6b9831d059fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?In=C3=A9s=20Chavarr=C3=ADa?= <72545702+ineschh@users.noreply.github.com> Date: Fri, 13 Nov 2020 10:10:53 +0100 Subject: [PATCH 08/52] More updates --- phys2denoise/cli/run.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/phys2denoise/cli/run.py b/phys2denoise/cli/run.py index bfee4ac..4a75ac7 100644 --- a/phys2denoise/cli/run.py +++ b/phys2denoise/cli/run.py @@ -113,13 +113,13 @@ def _get_parser(): optional.add_argument("-os", "--oversampling", dest="oversampling", type=int, - help="Temporal oversampling factor in seconds. " + help="Temporal oversampling factor. " "Default is 50.", default=50) optional.add_argument("-tl", "--time-length", dest="time_length", type=int, - help="RRF Kernel length in seconds.", + help="RRF or CRF Kernel length in seconds.", default=None) optional.add_argument("-onset", "--onset", dest="onset", From 8bd4b8f0fa1351cab026bf180b3bc6f8ed14a6d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?In=C3=A9s=20Chavarr=C3=ADa?= <72545702+ineschh@users.noreply.github.com> Date: Fri, 13 Nov 2020 13:34:56 +0100 Subject: [PATCH 09/52] Suggested changes --- phys2denoise/cli/run.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/phys2denoise/cli/run.py b/phys2denoise/cli/run.py index 4a75ac7..09e0ae5 100644 --- a/phys2denoise/cli/run.py +++ b/phys2denoise/cli/run.py @@ -59,13 +59,15 @@ def _get_parser(): help="Respiratory variance. Needs the following inputs: " "sample-rate, window and lags.", default=False) + """ metric.add_argument("-rvt", "--respiratory-volume-per-time", dest="metrics", action="append_const", const="rvt", help="Respiratory volume-per-time. Needs the following inputs: " - "sample-rate, window and lags.", + "sample-rate, window, lags, peaks and troughs.", default=False) + """ metric.add_argument("-rrf", "--respiratory-response-function", dest="metrics", action="append_const", @@ -104,10 +106,10 @@ def _get_parser(): help="Full path and filename of the list with the indexed peaks' " "positions of the physiological data.", default=None) - optional.add_argument("-tg", "--throughts", - dest="throughts", + optional.add_argument("-tg", "--troughs", + dest="troughs", type=str, - help="Full path and filename of the list with the indexed peaks' " + help="Full path and filename of the list with the indexed troughs' " "positions of the physiological data.", default=None) optional.add_argument("-os", "--oversampling", From d877f83877f5fd53de3c53c4d75db32369534b3b Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 13 Nov 2020 09:41:08 -0500 Subject: [PATCH 10/52] Add duecredit and some references. --- phys2denoise/due.py | 74 ++++++++++++++++++++++++++++++++++++++ phys2denoise/references.py | 10 ++++++ 2 files changed, 84 insertions(+) create mode 100644 phys2denoise/due.py create mode 100644 phys2denoise/references.py diff --git a/phys2denoise/due.py b/phys2denoise/due.py new file mode 100644 index 0000000..9a1c4dd --- /dev/null +++ b/phys2denoise/due.py @@ -0,0 +1,74 @@ +# emacs: at the end of the file +# ex: set sts=4 ts=4 sw=4 et: +# ## ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### # +""" + +Stub file for a guaranteed safe import of duecredit constructs: if duecredit +is not available. + +To use it, place it into your project codebase to be imported, e.g. copy as + + cp stub.py /path/tomodule/module/due.py + +Note that it might be better to avoid naming it duecredit.py to avoid shadowing +installed duecredit. + +Then use in your code as + + from .due import due, Doi, BibTeX, Text + +See https://github.com/duecredit/duecredit/blob/master/README.md for examples. + +Origin: Originally a part of the duecredit +Copyright: 2015-2019 DueCredit developers +License: BSD-2 +""" + +__version__ = '0.0.8' + + +class InactiveDueCreditCollector(object): + """Just a stub at the Collector which would not do anything""" + def _donothing(self, *args, **kwargs): + """Perform no good and no bad""" + pass + + def dcite(self, *args, **kwargs): + """If I could cite I would""" + def nondecorating_decorator(func): + return func + return nondecorating_decorator + + active = False + activate = add = cite = dump = load = _donothing + + def __repr__(self): + return self.__class__.__name__ + '()' + + +def _donothing_func(*args, **kwargs): + """Perform no good and no bad""" + pass + + +try: + from duecredit import due, BibTeX, Doi, Url, Text + if 'due' in locals() and not hasattr(due, 'cite'): + raise RuntimeError( + "Imported due lacks .cite. DueCredit is now disabled") +except Exception as e: + if not isinstance(e, ImportError): + import logging + logging.getLogger("duecredit").error( + "Failed to import duecredit due to %s" % str(e)) + # Initiate due stub + due = InactiveDueCreditCollector() + BibTeX = Doi = Url = Text = _donothing_func + +# Emacs mode definitions +# Local Variables: +# mode: python +# py-indent-offset: 4 +# tab-width: 4 +# indent-tabs-mode: nil +# End: diff --git a/phys2denoise/references.py b/phys2denoise/references.py new file mode 100644 index 0000000..305d231 --- /dev/null +++ b/phys2denoise/references.py @@ -0,0 +1,10 @@ +""" +References to be imported and injected throughout the package +""" +from .due import Doi + +CHANG_GLOVER_2009 = Doi('10.1016/j.neuroimage.2009.04.048') + +POWER_2018 = Doi('10.1073/pnas.1720985115') + +POWER_2020 = Doi('10.1016/j.neuroimage.2019.116234') From 869ea49a8680fc78ebff76ff59c75ebe7d1d3809 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 13 Nov 2020 09:41:15 -0500 Subject: [PATCH 11/52] Draft some metrics. --- phys2denoise/metrics/cardiac.py | 62 ++++++++ phys2denoise/metrics/chest_belt.py | 222 +++++++++++++++++++++++++++++ phys2denoise/metrics/gas.py | 2 + phys2denoise/metrics/utils.py | 93 ++++++++++++ 4 files changed, 379 insertions(+) create mode 100644 phys2denoise/metrics/cardiac.py create mode 100644 phys2denoise/metrics/chest_belt.py create mode 100644 phys2denoise/metrics/gas.py create mode 100644 phys2denoise/metrics/utils.py diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py new file mode 100644 index 0000000..a0a0b2b --- /dev/null +++ b/phys2denoise/metrics/cardiac.py @@ -0,0 +1,62 @@ +"""Denoising metrics for cardio recordings +""" +import numpy as np + +from ..due import due +from .. import references + + +def iht(): + """Instantaneous heart rate + """ + pass + + +@due.dcite(references.CHANG_GLOVER_2009) +def crf(samplerate, oversampling=50, time_length=32, onset=0., tr=2.): + """ + Calculate the cardiac response function using the definition + supplied in Chang and Glover (2009). + + Inputs + ------ + samplerate : :obj:`float` + Sampling rate of data, in seconds. + oversampling : :obj:`int`, optional + Temporal oversampling factor, in seconds. Default is 50. + time_length : :obj:`int`, optional + RRF kernel length, in seconds. Default is 32. + onset : :obj:`float`, optional + Onset of the response, in seconds. Default is 0. + + Outputs + ------- + crf: array-like + Cardiac or "heart" response function + + Notes + ----- + This cardiac response function was defined in [1]_, Appendix A. + + The core code for this function comes from metco2, while several of the + parameters, including oversampling, time_length, and onset, are modeled on + nistats' HRF functions. + + References + ---------- + .. [1] C. Chang & G. H. Glover, "Relationship between respiration, + end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage, + issue 47, vol. 4, pp. 1381-1393, 2009. + """ + def _crf(t): + rf = (0.6 * t ** 2.7 * np.exp(-t / 1.6) - 16 * (1 / np.sqrt(2 * np.pi * 9)) + * np.exp(-0.5 * (((t - 12) ** 2) / 9))) + return rf + + dt = tr / oversampling + time_stamps = np.linspace(0, time_length, + np.rint(float(time_length) / dt).astype(np.int)) + time_stamps -= onset + crf_arr = _crf(time_stamps) + crf_arr = crf_arr / max(abs(crf_arr)) + return crf_arr diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py new file mode 100644 index 0000000..f488a3d --- /dev/null +++ b/phys2denoise/metrics/chest_belt.py @@ -0,0 +1,222 @@ +"""Denoising metrics for chest belt recordings +""" +import numpy as np +import pandas as pd +from scipy.ndimage.filters import convolve1d +from scipy.signal import resample, detrend +from scipy.stats import zscore + +from . import utils +from ..due import due +from .. import references + + +@due.dcite(references.POWER_2018) +def rpv(belt_ts, window): + """Respiratory pattern variability + + Parameters + ---------- + belt_ts + window + + Returns + ------- + rpv_arr + + Notes + ----- + This metric was first introduced in [1]_. + + 1. Z-score respiratory belt signal + 2. Calculate upper envelope + 3. Calculate standard deviation of envelope + + References + ---------- + .. [1] J. D. Power et al., "Ridding fMRI data of motion-related influences: + Removal of signals with distinct spatial and physical bases in multiecho + data," Proceedings of the National Academy of Sciences, issue 9, vol. + 115, pp. 2105-2114, 2018. + """ + # First, z-score respiratory traces + resp_z = zscore(belt_ts) + + # Collect upper envelope + rpv_upper_env = utils.rms_envelope_1d(resp_z, window) + + # Calculate standard deviation + rpv_arr = np.std(rpv_upper_env) + return rpv_arr + + +@due.dcite(references.POWER_2020) +def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)): + """Respiratory pattern variability calculated across a sliding window + + Parameters + ---------- + belt_ts : (X,) :obj:`numpy.ndarray` + A 1D array with the respiratory belt time series. + samplerate : :obj:`float` + Sampling rate for belt_ts, in Hertz. + out_samplerate : :obj:`float` + Sampling rate for the output time series in seconds. + Corresponds to TR in fMRI data. + window : :obj:`int`, optional + Size of the sliding window, in the same units as out_samplerate. + Default is 6. + lags : (Y,) :obj:`tuple` of :obj:`int`, optional + List of lags to apply to the rv estimate. In the same units as + out_samplerate. + + Returns + ------- + env_arr + + Notes + ----- + This metric was first introduced in [1]_. + + Across a sliding window, do the following: + 1. Z-score respiratory belt signal + 2. Calculate upper envelope + 3. Calculate standard deviation of envelope + + References + ---------- + .. [1] J. D. Power et al., "Characteristics of respiratory measures in + young adults scanned at rest, including systematic changes and 'missed' + deep breaths," Neuroimage, vol. 204, 2020. + """ + window = window * samplerate / out_samplerate + # Calculate RPV across a rolling window + env_arr = pd.Series(belt_ts).rolling(window=window, center=True).apply( + rpv, window=window) + env_arr[np.isnan(env_arr)] = 0. + return env_arr + + +@due.dcite(references.CHANG_GLOVER_2009) +def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)): + """Respiratory variance + + Parameters + ---------- + belt_ts : (X,) :obj:`numpy.ndarray` + A 1D array with the respiratory belt time series. + samplerate : :obj:`float` + Sampling rate for belt_ts, in Hertz. + out_samplerate : :obj:`float` + Sampling rate for the output time series in seconds. + Corresponds to TR in fMRI data. + window : :obj:`int`, optional + Size of the sliding window, in the same units as out_samplerate. + Default is 6. + lags : (Y,) :obj:`tuple` of :obj:`int`, optional + List of lags to apply to the rv estimate. In the same units as + out_samplerate. + + Returns + ------- + rv_out : (Z, 2Y) :obj:`numpy.ndarray` + Respiratory variance values, with lags applied, downsampled to + out_samplerate, convolved with an RRF, and detrended/normalized. + The first Y columns are not convolved with the RRF, while the second Y + columns are. + + Notes + ----- + Respiratory variance (RV) was introduced in [1]_, and consists of the + standard deviation of the respiratory trace within a 6-second window. + + This metric is often lagged back and/or forward in time and convolved with + a respiratory response function before being included in a GLM. + Regressors also often have mean and linear trends removed and are + standardized prior to regressions. + + References + ---------- + .. [1] C. Chang & G. H. Glover, "Relationship between respiration, + end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage, + issue 4, vol. 47, pp. 1381-1393, 2009. + """ + # Raw respiratory variance + rv_arr = pd.Series(belt_ts).rolling(window=window, center=True).std() + rv_arr[np.isnan(rv_arr)] = 0. + + # Apply lags + n_out_samples = int((belt_ts.shape[0] / samplerate) / out_samplerate) + # convert lags from out_samplerate to samplerate + delays = [abs(int(lag * samplerate)) for lag in lags] + rv_with_lags = utils.apply_lags(rv_arr, lags=delays) + + # Downsample to out_samplerate + rv_with_lags = resample(rv_with_lags, num=n_out_samples, axis=0) + + # Convolve with rrf + rrf_arr = rrf(out_samplerate, oversampling=1) + rv_convolved = convolve1d(rv_with_lags, rrf_arr, axis=0) + + # Concatenate the raw and convolved versions + rv_combined = np.hstack((rv_with_lags, rv_convolved)) + + # Detrend and normalize + rv_combined = rv_combined - np.mean(rv_combined, axis=0) + rv_combined = detrend(rv_combined, axis=0) + rv_out = zscore(rv_combined, axis=0) + return rv_out + + +def rvt(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)): + """Respiratory volume-per-time + """ + pass + + +@due.dcite(references.CHANG_GLOVER_2009) +def rrf(samplerate, oversampling=50, time_length=50, onset=0., tr=2.): + """ + Calculate the respiratory response function using the definition + supplied in Chang and Glover (2009). + + Inputs + ------ + samplerate : :obj:`float` + Sampling rate of data, in seconds. + oversampling : :obj:`int`, optional + Temporal oversampling factor, in seconds. Default is 50. + time_length : :obj:`int`, optional + RRF kernel length, in seconds. Default is 50. + onset : :obj:`float`, optional + Onset of the response, in seconds. Default is 0. + + Outputs + ------- + rrf: array-like + respiratory response function + + Notes + ----- + This respiratory response function was defined in [1]_, Appendix A. + + The core code for this function comes from metco2, while several of the + parameters, including oversampling, time_length, and onset, are modeled on + nistats' HRF functions. + + References + ---------- + .. [1] C. Chang & G. H. Glover, "Relationship between respiration, + end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage, + issue 47, vol. 4, pp. 1381-1393, 2009. + """ + def _rrf(t): + rf = (0.6 * t ** 2.1 * np.exp(-t / 1.6) - 0.0023 * t ** 3.54 * np.exp(-t / 4.25)) + return rf + dt = tr / oversampling + time_stamps = np.linspace(0, time_length, + np.rint(float(time_length) / dt).astype(np.int)) + time_stamps -= onset + rrf_arr = _rrf(time_stamps) + rrf_arr = rrf_arr / max(abs(rrf_arr)) + return rrf_arr diff --git a/phys2denoise/metrics/gas.py b/phys2denoise/metrics/gas.py new file mode 100644 index 0000000..12ea419 --- /dev/null +++ b/phys2denoise/metrics/gas.py @@ -0,0 +1,2 @@ +"""Denoising metrics for gas recordings +""" diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py new file mode 100644 index 0000000..a55976b --- /dev/null +++ b/phys2denoise/metrics/utils.py @@ -0,0 +1,93 @@ +import numpy as np + + +def mirrorpad_1d(arr, buffer=250): + """ + Pad both sides of array with flipped values from array of length 'buffer'. + + Parameters + ---------- + arr + buffer + + Returns + ------- + arr_out + """ + mirror = np.flip(arr, axis=0) + idx = range(arr.shape[0] - buffer, arr.shape[0]) + pre_mirror = np.take(mirror, idx, axis=0) + idx = range(0, buffer) + post_mirror = np.take(mirror, idx, axis=0) + arr_out = np.concatenate((pre_mirror, arr, post_mirror), axis=0) + return arr_out + + +def rms_envelope_1d(arr, window=500): + """ + Conceptual translation of MATLAB 2017b's envelope(X, x, 'rms') function. + https://www.mathworks.com/help/signal/ref/envelope.html + + Parameters + ---------- + arr + window + + Returns + ------- + rms_env : numpy.ndarray + The upper envelope. + """ + assert arr.ndim == 1, 'Input data must be 1D' + assert window % 2 == 0, 'Window must be even' + n_t = arr.shape[0] + buf = int(window / 2) + + # Pad array at both ends + arr = np.copy(arr).astype(float) + mean = np.mean(arr) + arr -= mean + arr = mirrorpad_1d(arr, buffer=buf) + rms_env = np.empty(n_t) + for i in range(n_t): + # to match matlab + start_idx = i + buf + stop_idx = i + buf + window + + # but this is probably more appropriate + # start_idx = i + buf - 1 + # stop_idx = i + buf + window + window_arr = arr[start_idx:stop_idx] + rms = np.sqrt(np.mean(window_arr ** 2)) + rms_env[i] = rms + rms_env += mean + return rms_env + + +def apply_lags(arr1d, lags): + """ + Apply delays (lags) to an array. + + Parameters + ---------- + arr1d : (X,) :obj:`numpy.ndarray` + One-dimensional array to apply delays to. + lags : (Y,) :obj:`tuple` or :obj:`int` + Delays, in the same units as arr1d, to apply to arr1d. Can be negative, + zero, or positive integers. + + Returns + ------- + arr_with_lags : (X, Y) :obj:`numpy.ndarray` + arr1d shifted according to lags. Each column corresponds to a lag. + """ + arr_with_lags = np.zeros((arr1d.shape[0], len(lags))) + for i_lag, lag in enumerate(lags): + if lag < 0: + arr_delayed = np.hstack((arr1d[lag:], np.zeros(lag))) + elif lag > 0: + arr_delayed = np.hstack((np.zeros(lag), arr1d[lag:])) + else: + arr_delayed = arr1d.copy() + arr_with_lags[:, i_lag] = arr_delayed + return arr_with_lags From 7b950b62754a3d5105d3d382d2814527a9e7972c Mon Sep 17 00:00:00 2001 From: smoia Date: Fri, 13 Nov 2020 15:49:13 +0100 Subject: [PATCH 12/52] Import metrics --- phys2denoise.py | 40 +++++++++------------------------------- 1 file changed, 9 insertions(+), 31 deletions(-) diff --git a/phys2denoise.py b/phys2denoise.py index 17a6ed5..760eb8f 100644 --- a/phys2denoise.py +++ b/phys2denoise.py @@ -19,9 +19,11 @@ import numpy as np -from phys2denoise import utils, viz, _version +from phys2denoise import utils, _version from phys2denoise.cli.run import _get_parser -# from phys2denoise.metrics import cardiac, chest_belt, retroicor +from phys2denoise.metrics.cardiac import crf +from phys2denoise.metrics.chest_belt import rpv, rv, rvt, rrf +from phys2denoise.metrics.retroicor import compute_retroicor_regressors from . import __version__ from .due import due, Doi @@ -29,41 +31,13 @@ LGR = logging.getLogger(__name__) -def print_json(outfile, samp_freq, time_offset, ch_name): - """ - Print the json required by BIDS format. - - Parameters - ---------- - outfile: str or path - Fullpath to output file. - samp_freq: float - Frequency of sampling for the output file. - time_offset: float - Difference between beginning of file and first TR. - ch_name: list of str - List of channel names, as specified by BIDS format. - - Notes - ----- - Outcome: - outfile: .json file - File containing information for BIDS. - """ - start_time = -time_offset - summary = dict(SamplingFrequency=samp_freq, - StartTime=round(start_time, 4), - Columns=ch_name) - utils.writejson(outfile, summary, indent=4, sort_keys=False) - - @due.dcite( Doi(''), path='phys2denoise', description='Creation of regressors for physiological denoising', version=__version__, cite_module=True) -def phys2denoise(filename, outdir='.', debug=False, quiet=False): +def phys2denoise(filename, outdir='.', metrics=[], debug=False, quiet=False): """ Run main workflow of phys2denoise. @@ -129,7 +103,11 @@ def phys2denoise(filename, outdir='.', debug=False, quiet=False): # Read input file phys_in = np.genfromtxt(filename) + # Goes through the list of metrics and calls them + if not metrics: + metrics = ['crf', 'rpv', 'rv', 'rvt', 'rrf', 'rcard', 'r'] + for From 6090b6f5bfd746ccf6a50cdf08df48b1258ab9df Mon Sep 17 00:00:00 2001 From: smoia Date: Sat, 14 Nov 2020 11:38:06 +0100 Subject: [PATCH 13/52] First skeleton version --- phys2denoise.py | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/phys2denoise.py b/phys2denoise.py index 760eb8f..f89d005 100644 --- a/phys2denoise.py +++ b/phys2denoise.py @@ -18,6 +18,7 @@ from shutil import copy as cp import numpy as np +import pandas as pd from phys2denoise import utils, _version from phys2denoise.cli.run import _get_parser @@ -101,18 +102,34 @@ def phys2denoise(filename, outdir='.', metrics=[], debug=False, quiet=False): raise FileNotFoundError(f'The file {filename} does not exist!') # Read input file - phys_in = np.genfromtxt(filename) - - # Goes through the list of metrics and calls them - if not metrics: - metrics = ['crf', 'rpv', 'rv', 'rvt', 'rrf', 'rcard', 'r'] - - for - - + physio = np.genfromtxt(filename) + # Prepare pandas dataset + regr = pd.DataFrame() + # If no metrics was specified, calls all of them. + if not metrics: + metrics = ['crf', 'rpv', 'rv', 'rvt', 'rrf', 'retroicor_card', 'retroicor_resp'] + # Goes through the list of metrics and calls them + for metric in metrics: + if metrics == 'retroicor_card': + regr['retroicor_card'] = compute_retroicor_regressors(physio, + vars(metric_args), + card=True) + elif metrics == 'retroicor_resp': + regr['retroicor_resp'] = compute_retroicor_regressors(physio, + vars(metric_args), + resp=True) + else: + regr[f'{metric}'] = metric(physio, vars(metric_args)) + + #!# Add regressors visualisation + + # Export regressors and sidecar + out_filename = os.join(outdir, 'derivatives', filename) + regr.to_csv(out_filename, sep='\t', index=False, float_format='%.6e') + #!# Add sidecar export def _main(argv=None): From 195c369bbe0761a9624b561fba16ff005f3aabc2 Mon Sep 17 00:00:00 2001 From: smoia Date: Sat, 14 Nov 2020 13:35:30 +0100 Subject: [PATCH 14/52] Remove unused libraries and improve function argument --- phys2denoise.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/phys2denoise.py b/phys2denoise.py index f89d005..74db30a 100644 --- a/phys2denoise.py +++ b/phys2denoise.py @@ -14,13 +14,11 @@ import logging import os import sys -from copy import deepcopy -from shutil import copy as cp import numpy as np import pandas as pd -from phys2denoise import utils, _version +from phys2denoise import _version from phys2denoise.cli.run import _get_parser from phys2denoise.metrics.cardiac import crf from phys2denoise.metrics.chest_belt import rpv, rv, rvt, rrf @@ -38,7 +36,9 @@ description='Creation of regressors for physiological denoising', version=__version__, cite_module=True) -def phys2denoise(filename, outdir='.', metrics=[], debug=False, quiet=False): +def phys2denoise(filename, outdir='.', + metrics=[crf, rpv, rv, rvt, rrf, 'retroicor_card', 'retroicor_resp'], + debug=False, quiet=False): """ Run main workflow of phys2denoise. @@ -107,10 +107,6 @@ def phys2denoise(filename, outdir='.', metrics=[], debug=False, quiet=False): # Prepare pandas dataset regr = pd.DataFrame() - # If no metrics was specified, calls all of them. - if not metrics: - metrics = ['crf', 'rpv', 'rv', 'rvt', 'rrf', 'retroicor_card', 'retroicor_resp'] - # Goes through the list of metrics and calls them for metric in metrics: if metrics == 'retroicor_card': From 298acaf4162996db41f5b2747c2984351accc06c Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Wed, 25 Nov 2020 12:21:24 +0100 Subject: [PATCH 15/52] Input args to metrics using inspect.signature() --- phys2denoise.py | 44 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 41 insertions(+), 3 deletions(-) diff --git a/phys2denoise.py b/phys2denoise.py index 74db30a..514ba07 100644 --- a/phys2denoise.py +++ b/phys2denoise.py @@ -14,6 +14,7 @@ import logging import os import sys +from inspect import signature import numpy as np import pandas as pd @@ -30,6 +31,40 @@ LGR = logging.getLogger(__name__) +def select_input_args(metric, metric_args): + """ + Retrieve required args for metric from a dictionary of possible arguments. + + Parameters + ---------- + metric : function + Metric function to retrieve arguments for + metric_args : dict + Dictionary containing all arguments for all functions requested by the + user + + Returns + ------- + args : dict + Arguments to provide as input to metric + + Raises + ------ + ValueError + If a required argument is missing + + """ + req_args = [str(arg) for arg in signature(metric).parameters.values()] + + for arg in req_args: + if arg not in metric_args: + raise ValueError(f'Missing parameter {arg} required to run {metric}') + + args = {arg: metric_args[arg] for arg in req_args} + + return args + + @due.dcite( Doi(''), path='phys2denoise', @@ -110,15 +145,18 @@ def phys2denoise(filename, outdir='.', # Goes through the list of metrics and calls them for metric in metrics: if metrics == 'retroicor_card': + args = select_input_args(compute_retroicor_regressors, metric_args) regr['retroicor_card'] = compute_retroicor_regressors(physio, - vars(metric_args), + **args, card=True) elif metrics == 'retroicor_resp': + args = select_input_args(compute_retroicor_regressors, metric_args) regr['retroicor_resp'] = compute_retroicor_regressors(physio, - vars(metric_args), + **args, resp=True) else: - regr[f'{metric}'] = metric(physio, vars(metric_args)) + args = select_input_args(metric, metric_args) + regr[f'{metric}'] = metric(physio, **args) #!# Add regressors visualisation From bd340be1324203496a4a16e42141b536f8e67c22 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Wed, 25 Nov 2020 12:27:33 +0100 Subject: [PATCH 16/52] Change log path --- phys2denoise.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/phys2denoise.py b/phys2denoise.py index 514ba07..262b395 100644 --- a/phys2denoise.py +++ b/phys2denoise.py @@ -71,7 +71,7 @@ def select_input_args(metric, metric_args): description='Creation of regressors for physiological denoising', version=__version__, cite_module=True) -def phys2denoise(filename, outdir='.', +def phys2denoise(filename, metric_args, outdir='.', metrics=[crf, rpv, rv, rvt, rrf, 'retroicor_card', 'retroicor_resp'], debug=False, quiet=False): """ @@ -89,14 +89,14 @@ def phys2denoise(filename, outdir='.', outdir = os.path.abspath(outdir) os.makedirs(outdir) os.makedirs(os.path.join(outdir, 'code')) - conversion_path = os.path.join(outdir, 'code', 'conversion') - os.makedirs(conversion_path) + log_path = os.path.join(outdir, 'code', 'logs') + os.makedirs(log_path) # Create logfile name basename = 'phys2denoise_' extension = 'tsv' isotime = datetime.datetime.now().strftime('%Y-%m-%dT%H%M%S') - logname = os.path.join(conversion_path, (basename + isotime + '.' + extension)) + logname = os.path.join(log_path, (basename + isotime + '.' + extension)) # Set logging format log_formatter = logging.Formatter( @@ -110,13 +110,16 @@ def phys2denoise(filename, outdir='.', if quiet: logging.basicConfig(level=logging.WARNING, - handlers=[log_handler, sh], format='%(levelname)-10s %(message)s') + handlers=[log_handler, sh], + format='%(levelname)-10s %(message)s') elif debug: logging.basicConfig(level=logging.DEBUG, - handlers=[log_handler, sh], format='%(levelname)-10s %(message)s') + handlers=[log_handler, sh], + format='%(levelname)-10s %(message)s') else: logging.basicConfig(level=logging.INFO, - handlers=[log_handler, sh], format='%(levelname)-10s %(message)s') + handlers=[log_handler, sh], + format='%(levelname)-10s %(message)s') version_number = _version.get_versions()['version'] LGR.info(f'Currently running phys2denoise version {version_number}') @@ -125,7 +128,7 @@ def phys2denoise(filename, outdir='.', # Save call.sh arg_str = ' '.join(sys.argv[1:]) call_str = f'phys2denoise {arg_str}' - f = open(os.path.join(conversion_path, 'call.sh'), "a") + f = open(os.path.join(log_path, 'call.sh'), "a") f.write(f'#!bin/bash \n{call_str}') f.close() From cb8b42ef99870fe5d46931073bcbcaeb747267d4 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Wed, 25 Nov 2020 23:13:00 +0100 Subject: [PATCH 17/52] Change inspection of metrics and add a logger function --- phys2denoise.py | 62 ++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 51 insertions(+), 11 deletions(-) diff --git a/phys2denoise.py b/phys2denoise.py index 262b395..3e2354d 100644 --- a/phys2denoise.py +++ b/phys2denoise.py @@ -14,7 +14,7 @@ import logging import os import sys -from inspect import signature +from inspect import signature, _empty import numpy as np import pandas as pd @@ -35,6 +35,12 @@ def select_input_args(metric, metric_args): """ Retrieve required args for metric from a dictionary of possible arguments. + This function checks what parameters are accepted by a metric. + Then, for each parameter, check if the user provided it or not. + If they did not, but the parameter is required, throw an error - + unless it's "physio", reserved name for the timeseries input to a metric. + Otherwise, use the default. + Parameters ---------- metric : function @@ -54,15 +60,46 @@ def select_input_args(metric, metric_args): If a required argument is missing """ - req_args = [str(arg) for arg in signature(metric).parameters.values()] + args = {} + + # Check the parameters required by the metric and given by the user (see docstring) + for param in signature(metric).parameters.values(): + if param.name not in metric_args: + if param.default == _empty and param.name != 'physio': + raise ValueError(f'Missing parameter {param} required ' + f'to run {metric}') + else: + args[param.name] = param.default + else: + args[param.name] = metric_args[param.name] - for arg in req_args: - if arg not in metric_args: - raise ValueError(f'Missing parameter {arg} required to run {metric}') + return args - args = {arg: metric_args[arg] for arg in req_args} - return args +def print_metric_call(metric, args): + """ + Log a message to describe how a metric is being called. + + Parameters + ---------- + metric : function + Metric function that is being called + args : dict + Dictionary containing all arguments that are used to parametrise metric + + Notes + ----- + Outcome + An info-level message for the logger. + """ + msg = f'The {metric} regressor will be computed using the following parameters:' + + for arg in args: + msg = f'{msg}\n {arg} = {args[arg]}' + + msg = f'{msg}\n' + + LGR.info(msg) @due.dcite( @@ -149,16 +186,19 @@ def phys2denoise(filename, metric_args, outdir='.', for metric in metrics: if metrics == 'retroicor_card': args = select_input_args(compute_retroicor_regressors, metric_args) + args['card'] = True + print_metric_call(metric, args) regr['retroicor_card'] = compute_retroicor_regressors(physio, - **args, - card=True) + **args) elif metrics == 'retroicor_resp': args = select_input_args(compute_retroicor_regressors, metric_args) + args['resp'] = True + print_metric_call(metric, args) regr['retroicor_resp'] = compute_retroicor_regressors(physio, - **args, - resp=True) + **args) else: args = select_input_args(metric, metric_args) + print_metric_call(metric, args) regr[f'{metric}'] = metric(physio, **args) #!# Add regressors visualisation From 9c1bd601fe6abaac97c6e8ceb6f001449fd8a0a8 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Wed, 25 Nov 2020 23:17:21 +0100 Subject: [PATCH 18/52] Recognising planning contribution to @62442katieb Co-authored-by: @62442katieb From 1541149d95b7b7268af88b85604d274ceb738579 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?In=C3=A9s=20Chavarr=C3=ADa?= <72545702+ineschh@users.noreply.github.com> Date: Thu, 26 Nov 2020 15:03:48 +0100 Subject: [PATCH 19/52] metric_args into dict --- phys2denoise/cli/run.py | 188 +++++++++++++++++++++++----------------- 1 file changed, 109 insertions(+), 79 deletions(-) diff --git a/phys2denoise/cli/run.py b/phys2denoise/cli/run.py index 09e0ae5..2776756 100644 --- a/phys2denoise/cli/run.py +++ b/phys2denoise/cli/run.py @@ -5,6 +5,24 @@ import argparse from phys2denoise import __version__ +from phys2denoise.metrics.cardiac import crf +from phys2denoise.metrics.chest_belt import rpv, rv, rvt, rrf, env + + +class MetricsArgDict(argparse.Action): + """ + Custom Argparse Action to create a dictionary with the metrics' arguments in parser's output. + + """ + def __call__(self, parser, namespace, values, option_strings): + if not hasattr(namespace, "metrics_arg"): + setattr(namespace, "metrics_arg", dict()) + Keys = ["sample_rate", "peaks", "throughs", "oversampling", "time_length", "onset", + "tr", "window", "lags", "nscans", "nharm"] + Vals = ["None", "None", "None", "50", "None", "0", "None", "6", "None", "1", "None"] + for k, v in zip(Keys, Vals): + getattr(namespace, "metrics_arg")[k] = v + getattr(namespace, "metrics_arg")[self.dest] = values def _get_parser(): @@ -17,13 +35,16 @@ def _get_parser(): Notes ----- + Default values must be updated in __call__ method from MetricsArgDict class. # Argument parser follow template provided by RalphyZ. # https://stackoverflow.com/a/43456577 """ + parser = argparse.ArgumentParser() optional = parser._action_groups.pop() - metric = parser._action_groups.pop() - required = parser.add_argument_group("Required Argument:") + required = parser.add_argument_group("Required Argument") + metric = parser.add_argument_group("Metrics") + metric_arg = parser.add_argument_group("Metrics Arguments") required.add_argument("-in", "--input-file", dest="filename", type=str, @@ -33,32 +54,32 @@ def _get_parser(): metric.add_argument("-crf", "--cardiac-response-function", dest="metrics", action="append_const", - const="crf", + const=crf, help="Cardiac response function. Needs the following " "inputs:sample-rate, oversampling, time-length, " "onset and tr.", - default=False) + default=[]) metric.add_argument("-rpv", "--respiratory-pattern-variability", dest="metrics", action="append_const", - const="rpv", + const=rpv, help="Respiratory pattern variability. Needs the following " "input: window.", - default=False) + default=[]) metric.add_argument("-env", "--envelope", dest="metrics", action="append_const", - const="env", + const=env, help="Respiratory pattern variability calculated across a sliding " "window. Needs the following inputs: sample-rate, window and lags.", - default=False) + default=[]) metric.add_argument("-rv", "--respiratory-variance", dest="metrics", action="append_const", - const="rv", + const=rv, help="Respiratory variance. Needs the following inputs: " "sample-rate, window and lags.", - default=False) + default=[]) """ metric.add_argument("-rvt", "--respiratory-volume-per-time", dest="metrics", @@ -66,98 +87,108 @@ def _get_parser(): const="rvt", help="Respiratory volume-per-time. Needs the following inputs: " "sample-rate, window, lags, peaks and troughs.", - default=False) + default=[]) """ metric.add_argument("-rrf", "--respiratory-response-function", dest="metrics", action="append_const", - const="rrf", + const=rrf, help="Respiratory response function. Needs the following inputs: " "sample-rate, oversampling, time-length, onset and tr.", - default=False) + default=[]) metric.add_argument("-rcard", "--retroicor-card", dest="metrics", action="append_const", const="r_card", help="Computes regressors for cardiac signal. Needs the following " "inputs: tr, nscans and n_harm.", - default=False) + default=[]) metric.add_argument("-rresp", "--retroicor-resp", dest="metrics", action="append_const", const="r_resp", help="Computes regressors for respiratory signal. Needs the following " "inputs: tr, nscans and n_harm.", - default=False) + default=[]) optional.add_argument("-outdir", "--output-dir", dest="outdir", type=str, help="Folder where output should be placed. " "Default is current folder.", default=".") - optional.add_argument("-sr", "--sample-rate", - dest="sample_rate", - type=float, - help="Sampling rate of the physiological data in Hz.", - default=None) - optional.add_argument("-pk", "--peaks", - dest="peaks", - type=str, - help="Full path and filename of the list with the indexed peaks' " - "positions of the physiological data.", - default=None) - optional.add_argument("-tg", "--troughs", - dest="troughs", - type=str, - help="Full path and filename of the list with the indexed troughs' " - "positions of the physiological data.", - default=None) - optional.add_argument("-os", "--oversampling", - dest="oversampling", - type=int, - help="Temporal oversampling factor. " - "Default is 50.", - default=50) - optional.add_argument("-tl", "--time-length", - dest="time_length", - type=int, - help="RRF or CRF Kernel length in seconds.", - default=None) - optional.add_argument("-onset", "--onset", - dest="onset", - type=float, - help="Onset of the response in seconds. " - "Default is 0.", - default=0) - optional.add_argument("-tr", "--tr", - dest="tr", - type=float, - help="TR of sequence in seconds.", - default=None) - optional.add_argument("-win", "--window", - dest="window", - type=int, - help="Size of the sliding window in seconds. " - "Default is 6 seconds.", - default=6) - optional.add_argument("-lags", "--lags", - dest="lags", - nargs="*", - type=int, - action="append", - help="List of lags to apply to the RV estimate " - "in seconds.", - default=None) - optional.add_argument("-nscans", "--number-scans", - dest="nscans", - type=int, - help="Number of scans. Default is 1.", - default=1) - optional.add_argument("-nharm", "--number-harmonics", - dest="n_harm", - type=int, - help="Number of harmonics.", - default=None) + metric_arg.add_argument("-sr", "--sample-rate", + dest="sample_rate", + type=float, + action=MetricsArgDict, + help="Sampling rate of the physiological data in Hz.", + default=argparse.SUPPRESS) + metric_arg.add_argument("-pk", "--peaks", + dest="peaks", + type=str, + action=MetricsArgDict, + help="Full path and filename of the list with the indexed peaks' " + "positions of the physiological data.", + default=argparse.SUPPRESS) + metric_arg.add_argument("-tg", "--troughs", + dest="troughs", + type=str, + action=MetricsArgDict, + help="Full path and filename of the list with the indexed troughs' " + "positions of the physiological data.", + default=argparse.SUPPRESS) + metric_arg.add_argument("-os", "--oversampling", + dest="oversampling", + type=int, + action=MetricsArgDict, + help="Temporal oversampling factor. " + "Default is 50.", + default=argparse.SUPPRESS) + metric_arg.add_argument("-tl", "--time-length", + dest="time_length", + type=int, + action=MetricsArgDict, + help="RRF or CRF Kernel length in seconds.", + default=argparse.SUPPRESS) + metric_arg.add_argument("-onset", "--onset", + dest="onset", + type=float, + action=MetricsArgDict, + help="Onset of the response in seconds. " + "Default is 0.", + default=argparse.SUPPRESS) + metric_arg.add_argument("-tr", "--tr", + dest="tr", + type=float, + action=MetricsArgDict, + help="TR of sequence in seconds.", + default=argparse.SUPPRESS) + metric_arg.add_argument("-win", "--window", + dest="window", + type=int, + action=MetricsArgDict, + help="Size of the sliding window in seconds. " + "Default is 6 seconds.", + default=argparse.SUPPRESS) + metric_arg.add_argument("-lags", "--lags", + dest="lags", + nargs="*", + type=int, + action=MetricsArgDict, + help="List of lags to apply to the RV estimate " + "in seconds.", + default=argparse.SUPPRESS) + metric_arg.add_argument("-nscans", "--number-scans", + dest="nscans", + type=int, + action=MetricsArgDict, + help="Number of scans. Default is 1.", + default=argparse.SUPPRESS) + metric_arg.add_argument("-nharm", "--number-harmonics", + dest="n_harm", + type=int, + action=MetricsArgDict, + help="Number of harmonics.", + default=argparse.SUPPRESS) optional.add_argument("-debug", "--debug", dest="debug", action="store_true", @@ -172,7 +203,6 @@ def _get_parser(): version=("%(prog)s " + __version__)) parser._action_groups.append(optional) - parser._action_groups.append(metric) return parser From 23f851fb1f842b6a88bd5ac987c6a14a5d30002d Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 30 Dec 2020 12:16:10 -0500 Subject: [PATCH 20/52] Run black and isort. --- phys2denoise/due.py | 19 +++++++++++------- phys2denoise/metrics/cardiac.py | 18 +++++++++-------- phys2denoise/metrics/chest_belt.py | 32 ++++++++++++++---------------- phys2denoise/metrics/utils.py | 4 ++-- phys2denoise/references.py | 6 +++--- 5 files changed, 42 insertions(+), 37 deletions(-) diff --git a/phys2denoise/due.py b/phys2denoise/due.py index 9a1c4dd..6013ed7 100644 --- a/phys2denoise/due.py +++ b/phys2denoise/due.py @@ -24,26 +24,29 @@ License: BSD-2 """ -__version__ = '0.0.8' +__version__ = "0.0.8" class InactiveDueCreditCollector(object): """Just a stub at the Collector which would not do anything""" + def _donothing(self, *args, **kwargs): """Perform no good and no bad""" pass def dcite(self, *args, **kwargs): """If I could cite I would""" + def nondecorating_decorator(func): return func + return nondecorating_decorator active = False activate = add = cite = dump = load = _donothing def __repr__(self): - return self.__class__.__name__ + '()' + return self.__class__.__name__ + "()" def _donothing_func(*args, **kwargs): @@ -52,15 +55,17 @@ def _donothing_func(*args, **kwargs): try: - from duecredit import due, BibTeX, Doi, Url, Text - if 'due' in locals() and not hasattr(due, 'cite'): - raise RuntimeError( - "Imported due lacks .cite. DueCredit is now disabled") + from duecredit import BibTeX, Doi, Text, Url, due + + if "due" in locals() and not hasattr(due, "cite"): + raise RuntimeError("Imported due lacks .cite. DueCredit is now disabled") except Exception as e: if not isinstance(e, ImportError): import logging + logging.getLogger("duecredit").error( - "Failed to import duecredit due to %s" % str(e)) + "Failed to import duecredit due to %s" % str(e) + ) # Initiate due stub due = InactiveDueCreditCollector() BibTeX = Doi = Url = Text = _donothing_func diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py index a0a0b2b..19ab81a 100644 --- a/phys2denoise/metrics/cardiac.py +++ b/phys2denoise/metrics/cardiac.py @@ -2,18 +2,17 @@ """ import numpy as np -from ..due import due from .. import references +from ..due import due def iht(): - """Instantaneous heart rate - """ + """Instantaneous heart rate""" pass @due.dcite(references.CHANG_GLOVER_2009) -def crf(samplerate, oversampling=50, time_length=32, onset=0., tr=2.): +def crf(samplerate, oversampling=50, time_length=32, onset=0.0, tr=2.0): """ Calculate the cardiac response function using the definition supplied in Chang and Glover (2009). @@ -48,14 +47,17 @@ def crf(samplerate, oversampling=50, time_length=32, onset=0., tr=2.): end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage, issue 47, vol. 4, pp. 1381-1393, 2009. """ + def _crf(t): - rf = (0.6 * t ** 2.7 * np.exp(-t / 1.6) - 16 * (1 / np.sqrt(2 * np.pi * 9)) - * np.exp(-0.5 * (((t - 12) ** 2) / 9))) + rf = 0.6 * t ** 2.7 * np.exp(-t / 1.6) - 16 * ( + 1 / np.sqrt(2 * np.pi * 9) + ) * np.exp(-0.5 * (((t - 12) ** 2) / 9)) return rf dt = tr / oversampling - time_stamps = np.linspace(0, time_length, - np.rint(float(time_length) / dt).astype(np.int)) + time_stamps = np.linspace( + 0, time_length, np.rint(float(time_length) / dt).astype(np.int) + ) time_stamps -= onset crf_arr = _crf(time_stamps) crf_arr = crf_arr / max(abs(crf_arr)) diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py index f488a3d..3d0e905 100644 --- a/phys2denoise/metrics/chest_belt.py +++ b/phys2denoise/metrics/chest_belt.py @@ -3,12 +3,12 @@ import numpy as np import pandas as pd from scipy.ndimage.filters import convolve1d -from scipy.signal import resample, detrend +from scipy.signal import detrend, resample from scipy.stats import zscore -from . import utils -from ..due import due from .. import references +from ..due import due +from . import utils @due.dcite(references.POWER_2018) @@ -91,9 +91,10 @@ def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)): """ window = window * samplerate / out_samplerate # Calculate RPV across a rolling window - env_arr = pd.Series(belt_ts).rolling(window=window, center=True).apply( - rpv, window=window) - env_arr[np.isnan(env_arr)] = 0. + env_arr = ( + pd.Series(belt_ts).rolling(window=window, center=True).apply(rpv, window=window) + ) + env_arr[np.isnan(env_arr)] = 0.0 return env_arr @@ -143,7 +144,7 @@ def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)): """ # Raw respiratory variance rv_arr = pd.Series(belt_ts).rolling(window=window, center=True).std() - rv_arr[np.isnan(rv_arr)] = 0. + rv_arr[np.isnan(rv_arr)] = 0.0 # Apply lags n_out_samples = int((belt_ts.shape[0] / samplerate) / out_samplerate) @@ -168,14 +169,8 @@ def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)): return rv_out -def rvt(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)): - """Respiratory volume-per-time - """ - pass - - @due.dcite(references.CHANG_GLOVER_2009) -def rrf(samplerate, oversampling=50, time_length=50, onset=0., tr=2.): +def rrf(samplerate, oversampling=50, time_length=50, onset=0.0, tr=2.0): """ Calculate the respiratory response function using the definition supplied in Chang and Glover (2009). @@ -210,12 +205,15 @@ def rrf(samplerate, oversampling=50, time_length=50, onset=0., tr=2.): end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage, issue 47, vol. 4, pp. 1381-1393, 2009. """ + def _rrf(t): - rf = (0.6 * t ** 2.1 * np.exp(-t / 1.6) - 0.0023 * t ** 3.54 * np.exp(-t / 4.25)) + rf = 0.6 * t ** 2.1 * np.exp(-t / 1.6) - 0.0023 * t ** 3.54 * np.exp(-t / 4.25) return rf + dt = tr / oversampling - time_stamps = np.linspace(0, time_length, - np.rint(float(time_length) / dt).astype(np.int)) + time_stamps = np.linspace( + 0, time_length, np.rint(float(time_length) / dt).astype(np.int) + ) time_stamps -= onset rrf_arr = _rrf(time_stamps) rrf_arr = rrf_arr / max(abs(rrf_arr)) diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py index a55976b..e44c2a5 100644 --- a/phys2denoise/metrics/utils.py +++ b/phys2denoise/metrics/utils.py @@ -38,8 +38,8 @@ def rms_envelope_1d(arr, window=500): rms_env : numpy.ndarray The upper envelope. """ - assert arr.ndim == 1, 'Input data must be 1D' - assert window % 2 == 0, 'Window must be even' + assert arr.ndim == 1, "Input data must be 1D" + assert window % 2 == 0, "Window must be even" n_t = arr.shape[0] buf = int(window / 2) diff --git a/phys2denoise/references.py b/phys2denoise/references.py index 305d231..c4b557f 100644 --- a/phys2denoise/references.py +++ b/phys2denoise/references.py @@ -3,8 +3,8 @@ """ from .due import Doi -CHANG_GLOVER_2009 = Doi('10.1016/j.neuroimage.2009.04.048') +CHANG_GLOVER_2009 = Doi("10.1016/j.neuroimage.2009.04.048") -POWER_2018 = Doi('10.1073/pnas.1720985115') +POWER_2018 = Doi("10.1073/pnas.1720985115") -POWER_2020 = Doi('10.1016/j.neuroimage.2019.116234') +POWER_2020 = Doi("10.1016/j.neuroimage.2019.116234") From a4454e296535bb1093e5328279239c5ade3c617a Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 30 Dec 2020 12:17:47 -0500 Subject: [PATCH 21/52] Minor changes. --- phys2denoise/due.py | 13 +++++-------- phys2denoise/references.py | 4 +--- 2 files changed, 6 insertions(+), 11 deletions(-) diff --git a/phys2denoise/due.py b/phys2denoise/due.py index 6013ed7..17071d4 100644 --- a/phys2denoise/due.py +++ b/phys2denoise/due.py @@ -1,10 +1,7 @@ # emacs: at the end of the file # ex: set sts=4 ts=4 sw=4 et: # ## ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### # -""" - -Stub file for a guaranteed safe import of duecredit constructs: if duecredit -is not available. +"""Stub file for a guaranteed safe import of duecredit constructs if duecredit is not available. To use it, place it into your project codebase to be imported, e.g. copy as @@ -28,14 +25,14 @@ class InactiveDueCreditCollector(object): - """Just a stub at the Collector which would not do anything""" + """Just a stub at the Collector which would not do anything.""" def _donothing(self, *args, **kwargs): - """Perform no good and no bad""" + """Perform no good and no bad.""" pass def dcite(self, *args, **kwargs): - """If I could cite I would""" + """If I could cite I would.""" def nondecorating_decorator(func): return func @@ -50,7 +47,7 @@ def __repr__(self): def _donothing_func(*args, **kwargs): - """Perform no good and no bad""" + """Perform no good and no bad.""" pass diff --git a/phys2denoise/references.py b/phys2denoise/references.py index c4b557f..297ba3b 100644 --- a/phys2denoise/references.py +++ b/phys2denoise/references.py @@ -1,6 +1,4 @@ -""" -References to be imported and injected throughout the package -""" +"""References to be imported and injected throughout the package.""" from .due import Doi CHANG_GLOVER_2009 = Doi("10.1016/j.neuroimage.2009.04.048") From 4612d4fe4a557dc88f256053baa400b076c8519c Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 30 Dec 2020 12:25:47 -0500 Subject: [PATCH 22/52] More fixing. --- phys2denoise/due.py | 1 + phys2denoise/metrics/cardiac.py | 9 +++------ phys2denoise/metrics/chest_belt.py | 13 +++++-------- phys2denoise/metrics/gas.py | 3 +-- phys2denoise/metrics/utils.py | 15 ++++++++------- 5 files changed, 18 insertions(+), 23 deletions(-) diff --git a/phys2denoise/due.py b/phys2denoise/due.py index 17071d4..4ff2457 100644 --- a/phys2denoise/due.py +++ b/phys2denoise/due.py @@ -43,6 +43,7 @@ def nondecorating_decorator(func): activate = add = cite = dump = load = _donothing def __repr__(self): + """Perform magic.""" return self.__class__.__name__ + "()" diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py index 19ab81a..98a5f23 100644 --- a/phys2denoise/metrics/cardiac.py +++ b/phys2denoise/metrics/cardiac.py @@ -1,5 +1,4 @@ -"""Denoising metrics for cardio recordings -""" +"""Denoising metrics for cardio recordings.""" import numpy as np from .. import references @@ -7,15 +6,13 @@ def iht(): - """Instantaneous heart rate""" + """Calculate instantaneous heart rate.""" pass @due.dcite(references.CHANG_GLOVER_2009) def crf(samplerate, oversampling=50, time_length=32, onset=0.0, tr=2.0): - """ - Calculate the cardiac response function using the definition - supplied in Chang and Glover (2009). + """Calculate the cardiac response function using Chang and Glover's definition. Inputs ------ diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py index 3d0e905..dc24794 100644 --- a/phys2denoise/metrics/chest_belt.py +++ b/phys2denoise/metrics/chest_belt.py @@ -1,5 +1,4 @@ -"""Denoising metrics for chest belt recordings -""" +"""Denoising metrics for chest belt recordings.""" import numpy as np import pandas as pd from scipy.ndimage.filters import convolve1d @@ -13,7 +12,7 @@ @due.dcite(references.POWER_2018) def rpv(belt_ts, window): - """Respiratory pattern variability + """Calculate respiratory pattern variability. Parameters ---------- @@ -52,7 +51,7 @@ def rpv(belt_ts, window): @due.dcite(references.POWER_2020) def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)): - """Respiratory pattern variability calculated across a sliding window + """Calculate respiratory pattern variability across a sliding window. Parameters ---------- @@ -100,7 +99,7 @@ def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)): @due.dcite(references.CHANG_GLOVER_2009) def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)): - """Respiratory variance + """Calculate respiratory variance. Parameters ---------- @@ -171,9 +170,7 @@ def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)): @due.dcite(references.CHANG_GLOVER_2009) def rrf(samplerate, oversampling=50, time_length=50, onset=0.0, tr=2.0): - """ - Calculate the respiratory response function using the definition - supplied in Chang and Glover (2009). + """Calculate the respiratory response function using Chang and Glover's definition. Inputs ------ diff --git a/phys2denoise/metrics/gas.py b/phys2denoise/metrics/gas.py index 12ea419..a4e3a59 100644 --- a/phys2denoise/metrics/gas.py +++ b/phys2denoise/metrics/gas.py @@ -1,2 +1 @@ -"""Denoising metrics for gas recordings -""" +"""Denoising metrics for gas recordings.""" diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py index e44c2a5..5369081 100644 --- a/phys2denoise/metrics/utils.py +++ b/phys2denoise/metrics/utils.py @@ -1,9 +1,9 @@ +"""Miscellaneous utility functions for metric calculation.""" import numpy as np def mirrorpad_1d(arr, buffer=250): - """ - Pad both sides of array with flipped values from array of length 'buffer'. + """Pad both sides of array with flipped values from array of length 'buffer'. Parameters ---------- @@ -24,9 +24,7 @@ def mirrorpad_1d(arr, buffer=250): def rms_envelope_1d(arr, window=500): - """ - Conceptual translation of MATLAB 2017b's envelope(X, x, 'rms') function. - https://www.mathworks.com/help/signal/ref/envelope.html + """Conceptual translation of MATLAB 2017b's envelope(X, x, 'rms') function. Parameters ---------- @@ -37,6 +35,10 @@ def rms_envelope_1d(arr, window=500): ------- rms_env : numpy.ndarray The upper envelope. + + Notes + ----- + https://www.mathworks.com/help/signal/ref/envelope.html """ assert arr.ndim == 1, "Input data must be 1D" assert window % 2 == 0, "Window must be even" @@ -65,8 +67,7 @@ def rms_envelope_1d(arr, window=500): def apply_lags(arr1d, lags): - """ - Apply delays (lags) to an array. + """Apply delays (lags) to an array. Parameters ---------- From 4202c2ec425d1c4eca34a35d1f956ef9af431d03 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 30 Dec 2020 14:26:02 -0500 Subject: [PATCH 23/52] Fix crf/rrf docstrings. --- phys2denoise/metrics/cardiac.py | 6 +++--- phys2denoise/metrics/chest_belt.py | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py index 98a5f23..df13d4d 100644 --- a/phys2denoise/metrics/cardiac.py +++ b/phys2denoise/metrics/cardiac.py @@ -14,8 +14,8 @@ def iht(): def crf(samplerate, oversampling=50, time_length=32, onset=0.0, tr=2.0): """Calculate the cardiac response function using Chang and Glover's definition. - Inputs - ------ + Parameters + ---------- samplerate : :obj:`float` Sampling rate of data, in seconds. oversampling : :obj:`int`, optional @@ -25,7 +25,7 @@ def crf(samplerate, oversampling=50, time_length=32, onset=0.0, tr=2.0): onset : :obj:`float`, optional Onset of the response, in seconds. Default is 0. - Outputs + Returns ------- crf: array-like Cardiac or "heart" response function diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py index dc24794..26c3c3f 100644 --- a/phys2denoise/metrics/chest_belt.py +++ b/phys2denoise/metrics/chest_belt.py @@ -172,18 +172,18 @@ def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)): def rrf(samplerate, oversampling=50, time_length=50, onset=0.0, tr=2.0): """Calculate the respiratory response function using Chang and Glover's definition. - Inputs - ------ + Parameters + ---------- samplerate : :obj:`float` Sampling rate of data, in seconds. oversampling : :obj:`int`, optional - Temporal oversampling factor, in seconds. Default is 50. + Temporal oversampling factor. Default is 50. time_length : :obj:`int`, optional RRF kernel length, in seconds. Default is 50. onset : :obj:`float`, optional Onset of the response, in seconds. Default is 0. - Outputs + Returns ------- rrf: array-like respiratory response function From e2116dc0951a01b2ddad7f8e5f66d80d87dca545 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 30 Dec 2020 14:26:58 -0500 Subject: [PATCH 24/52] Fix again. --- phys2denoise/metrics/cardiac.py | 2 +- phys2denoise/metrics/chest_belt.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py index df13d4d..f295848 100644 --- a/phys2denoise/metrics/cardiac.py +++ b/phys2denoise/metrics/cardiac.py @@ -27,7 +27,7 @@ def crf(samplerate, oversampling=50, time_length=32, onset=0.0, tr=2.0): Returns ------- - crf: array-like + crf : array-like Cardiac or "heart" response function Notes diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py index 26c3c3f..9ff8225 100644 --- a/phys2denoise/metrics/chest_belt.py +++ b/phys2denoise/metrics/chest_belt.py @@ -185,7 +185,7 @@ def rrf(samplerate, oversampling=50, time_length=50, onset=0.0, tr=2.0): Returns ------- - rrf: array-like + rrf : array-like respiratory response function Notes From 636f585e8c4cfb1b5f65971b36eae46886515342 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 15:46:43 +0100 Subject: [PATCH 25/52] Remove extra "version" imports --- phys2denoise.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/phys2denoise.py b/phys2denoise.py index 3e2354d..b1cc125 100644 --- a/phys2denoise.py +++ b/phys2denoise.py @@ -19,7 +19,6 @@ import numpy as np import pandas as pd -from phys2denoise import _version from phys2denoise.cli.run import _get_parser from phys2denoise.metrics.cardiac import crf from phys2denoise.metrics.chest_belt import rpv, rv, rvt, rrf @@ -158,7 +157,7 @@ def phys2denoise(filename, metric_args, outdir='.', handlers=[log_handler, sh], format='%(levelname)-10s %(message)s') - version_number = _version.get_versions()['version'] + version_number = __version__ LGR.info(f'Currently running phys2denoise version {version_number}') LGR.info(f'Input file is {filename}') From 36abfb80fac870bf569be8c90ed823d53fb80207 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 15:46:57 +0100 Subject: [PATCH 26/52] Change metric_args into **kwargs --- phys2denoise.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/phys2denoise.py b/phys2denoise.py index b1cc125..2fe6bd1 100644 --- a/phys2denoise.py +++ b/phys2denoise.py @@ -107,9 +107,9 @@ def print_metric_call(metric, args): description='Creation of regressors for physiological denoising', version=__version__, cite_module=True) -def phys2denoise(filename, metric_args, outdir='.', +def phys2denoise(filename, outdir='.', metrics=[crf, rpv, rv, rvt, rrf, 'retroicor_card', 'retroicor_resp'], - debug=False, quiet=False): + debug=False, quiet=False, **kwargs): """ Run main workflow of phys2denoise. @@ -117,14 +117,13 @@ def phys2denoise(filename, metric_args, outdir='.', Notes ----- + Any metric argument should go into kwargs! The code was greatly copied from phys2bids (copyright the physiopy community) """ # Check options to make them internally coherent pt. I # #!# This can probably be done while parsing? outdir = os.path.abspath(outdir) - os.makedirs(outdir) - os.makedirs(os.path.join(outdir, 'code')) log_path = os.path.join(outdir, 'code', 'logs') os.makedirs(log_path) @@ -184,19 +183,19 @@ def phys2denoise(filename, metric_args, outdir='.', # Goes through the list of metrics and calls them for metric in metrics: if metrics == 'retroicor_card': - args = select_input_args(compute_retroicor_regressors, metric_args) + args = select_input_args(compute_retroicor_regressors, kwargs) args['card'] = True print_metric_call(metric, args) regr['retroicor_card'] = compute_retroicor_regressors(physio, **args) elif metrics == 'retroicor_resp': - args = select_input_args(compute_retroicor_regressors, metric_args) + args = select_input_args(compute_retroicor_regressors, kwargs) args['resp'] = True print_metric_call(metric, args) regr['retroicor_resp'] = compute_retroicor_regressors(physio, **args) else: - args = select_input_args(metric, metric_args) + args = select_input_args(metric, kwargs) print_metric_call(metric, args) regr[f'{metric}'] = metric(physio, **args) From 889d3e57b9f17eedf127703a57138741f0c04de6 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 15:52:02 +0100 Subject: [PATCH 27/52] Move metric call function to metrics.utils --- phys2denoise/metrics/utils.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py index 5369081..7426a5c 100644 --- a/phys2denoise/metrics/utils.py +++ b/phys2denoise/metrics/utils.py @@ -1,6 +1,37 @@ """Miscellaneous utility functions for metric calculation.""" +import logging + import numpy as np +LGR = logging.getLogger(__name__) +LGR.setLevel(logging.INFO) + + +def print_metric_call(metric, args): + """ + Log a message to describe how a metric is being called. + + Parameters + ---------- + metric : function + Metric function that is being called + args : dict + Dictionary containing all arguments that are used to parametrise metric + + Notes + ----- + Outcome + An info-level message for the logger. + """ + msg = f'The {metric} regressor will be computed using the following parameters:' + + for arg in args: + msg = f'{msg}\n {arg} = {args[arg]}' + + msg = f'{msg}\n' + + LGR.info(msg) + def mirrorpad_1d(arr, buffer=250): """Pad both sides of array with flipped values from array of length 'buffer'. From 46c1e0c9c4c2ebe9e80f3bc7d35cb1ed23de47c2 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 15:57:29 +0100 Subject: [PATCH 28/52] Move metric call log to metrics.utils instead of in the main workflow. --- phys2denoise.py | 34 +++------------------------------- 1 file changed, 3 insertions(+), 31 deletions(-) diff --git a/phys2denoise.py b/phys2denoise.py index 2fe6bd1..37325ac 100644 --- a/phys2denoise.py +++ b/phys2denoise.py @@ -28,6 +28,7 @@ from .due import due, Doi LGR = logging.getLogger(__name__) +LGR.setLevel(logging.INFO) def select_input_args(metric, metric_args): @@ -75,32 +76,6 @@ def select_input_args(metric, metric_args): return args -def print_metric_call(metric, args): - """ - Log a message to describe how a metric is being called. - - Parameters - ---------- - metric : function - Metric function that is being called - args : dict - Dictionary containing all arguments that are used to parametrise metric - - Notes - ----- - Outcome - An info-level message for the logger. - """ - msg = f'The {metric} regressor will be computed using the following parameters:' - - for arg in args: - msg = f'{msg}\n {arg} = {args[arg]}' - - msg = f'{msg}\n' - - LGR.info(msg) - - @due.dcite( Doi(''), path='phys2denoise', @@ -182,21 +157,18 @@ def phys2denoise(filename, outdir='.', # Goes through the list of metrics and calls them for metric in metrics: - if metrics == 'retroicor_card': + if metric == 'retroicor_card': args = select_input_args(compute_retroicor_regressors, kwargs) args['card'] = True - print_metric_call(metric, args) regr['retroicor_card'] = compute_retroicor_regressors(physio, **args) - elif metrics == 'retroicor_resp': + elif metric == 'retroicor_resp': args = select_input_args(compute_retroicor_regressors, kwargs) args['resp'] = True - print_metric_call(metric, args) regr['retroicor_resp'] = compute_retroicor_regressors(physio, **args) else: args = select_input_args(metric, kwargs) - print_metric_call(metric, args) regr[f'{metric}'] = metric(physio, **args) #!# Add regressors visualisation From 44fde8c54f745214487c74f3d7e7632ed301ff4f Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 16:45:29 +0100 Subject: [PATCH 29/52] Move the bash call saving to a function on its own out of the main workflow --- phys2denoise.py | 32 +++++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/phys2denoise.py b/phys2denoise.py index 37325ac..b9cd65e 100644 --- a/phys2denoise.py +++ b/phys2denoise.py @@ -31,6 +31,28 @@ LGR.setLevel(logging.INFO) +def save_bash_call(outdir): + """ + Save the bash call into file `p2d_call.sh`. + + Parameters + ---------- + metric : function + Metric function to retrieve arguments for + metric_args : dict + Dictionary containing all arguments for all functions requested by the + user + """ + arg_str = ' '.join(sys.argv[1:]) + call_str = f'phys2denoise {arg_str}' + outdir = os.path.abspath(outdir) + log_path = os.path.join(outdir, 'code', 'logs') + os.makedirs(log_path) + f = open(os.path.join(log_path, 'p2d_call.sh'), "a") + f.write(f'#!bin/bash \n{call_str}') + f.close() + + def select_input_args(metric, metric_args): """ Retrieve required args for metric from a dictionary of possible arguments. @@ -135,13 +157,6 @@ def phys2denoise(filename, outdir='.', LGR.info(f'Currently running phys2denoise version {version_number}') LGR.info(f'Input file is {filename}') - # Save call.sh - arg_str = ' '.join(sys.argv[1:]) - call_str = f'phys2denoise {arg_str}' - f = open(os.path.join(log_path, 'call.sh'), "a") - f.write(f'#!bin/bash \n{call_str}') - f.close() - # Check options to make them internally coherent pt. II # #!# This can probably be done while parsing? # filename, ftype = utils.check_input_type(filename) @@ -181,6 +196,9 @@ def phys2denoise(filename, outdir='.', def _main(argv=None): options = _get_parser().parse_args(argv) + + save_bash_call(options.outdir) + phys2denoise(**vars(options)) From 570619dc0119b3958296094219198ba10cb074e6 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 16:49:04 +0100 Subject: [PATCH 30/52] Move main workflow in the right place --- phys2denoise.py => phys2denoise/phys2denoise.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename phys2denoise.py => phys2denoise/phys2denoise.py (100%) diff --git a/phys2denoise.py b/phys2denoise/phys2denoise.py similarity index 100% rename from phys2denoise.py rename to phys2denoise/phys2denoise.py From 3a66afacb181e3ba267a675c5c7098225619b42c Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 17:07:59 +0100 Subject: [PATCH 31/52] Add parser --- phys2denoise/cli/run.py | 164 ++++++++++++++++++++-------------------- 1 file changed, 81 insertions(+), 83 deletions(-) diff --git a/phys2denoise/cli/run.py b/phys2denoise/cli/run.py index 2776756..b6a7248 100644 --- a/phys2denoise/cli/run.py +++ b/phys2denoise/cli/run.py @@ -43,152 +43,150 @@ def _get_parser(): parser = argparse.ArgumentParser() optional = parser._action_groups.pop() required = parser.add_argument_group("Required Argument") - metric = parser.add_argument_group("Metrics") + metrics = parser.add_argument_group("Metrics") metric_arg = parser.add_argument_group("Metrics Arguments") + # Required arguments required.add_argument("-in", "--input-file", dest="filename", type=str, help="Full path and name of the file containing " "physiological data, with or without extension.", required=True) - metric.add_argument("-crf", "--cardiac-response-function", - dest="metrics", - action="append_const", - const=crf, - help="Cardiac response function. Needs the following " - "inputs:sample-rate, oversampling, time-length, " - "onset and tr.", - default=[]) - metric.add_argument("-rpv", "--respiratory-pattern-variability", - dest="metrics", - action="append_const", - const=rpv, - help="Respiratory pattern variability. Needs the following " - "input: window.", - default=[]) - metric.add_argument("-env", "--envelope", - dest="metrics", - action="append_const", - const=env, - help="Respiratory pattern variability calculated across a sliding " - "window. Needs the following inputs: sample-rate, window and lags.", - default=[]) - metric.add_argument("-rv", "--respiratory-variance", - dest="metrics", - action="append_const", - const=rv, - help="Respiratory variance. Needs the following inputs: " - "sample-rate, window and lags.", - default=[]) - """ - metric.add_argument("-rvt", "--respiratory-volume-per-time", - dest="metrics", - action="append_const", - const="rvt", - help="Respiratory volume-per-time. Needs the following inputs: " - "sample-rate, window, lags, peaks and troughs.", - default=[]) - """ - metric.add_argument("-rrf", "--respiratory-response-function", - dest="metrics", - action="append_const", - const=rrf, - help="Respiratory response function. Needs the following inputs: " - "sample-rate, oversampling, time-length, onset and tr.", - default=[]) - metric.add_argument("-rcard", "--retroicor-card", - dest="metrics", - action="append_const", - const="r_card", - help="Computes regressors for cardiac signal. Needs the following " - "inputs: tr, nscans and n_harm.", - default=[]) - metric.add_argument("-rresp", "--retroicor-resp", - dest="metrics", - action="append_const", - const="r_resp", - help="Computes regressors for respiratory signal. Needs the following " - "inputs: tr, nscans and n_harm.", - default=[]) + # Important optional arguments optional.add_argument("-outdir", "--output-dir", dest="outdir", type=str, help="Folder where output should be placed. " "Default is current folder.", default=".") + # Metric selection + metrics.add_argument("-crf", "--cardiac-response-function", + dest="metrics", + action="append_const", + const=crf, + help="Cardiac response function. Requires the following " + "inputs:sample-rate, oversampling, time-length, " + "onset and tr.", + default=[]) + metrics.add_argument("-rpv", "--respiratory-pattern-variability", + dest="metrics", + action="append_const", + const=rpv, + help="Respiratory pattern variability. Requires the following " + "input: window.", + default=[]) + metrics.add_argument("-env", "--envelope", + dest="metrics", + action="append_const", + const=env, + help="Respiratory pattern variability calculated across a sliding " + "window. Requires the following inputs: sample-rate, window and lags.", + default=[]) + metrics.add_argument("-rv", "--respiratory-variance", + dest="metrics", + action="append_const", + const=rv, + help="Respiratory variance. Requires the following inputs: " + "sample-rate, window and lags. If the input file " + "not a .phys file, it also requires peaks and troughs", + default=[]) + """ + metrics.add_argument("-rvt", "--respiratory-volume-per-time", + dest="metrics", + action="append_const", + const="rvt", + help="Respiratory volume-per-time. Requires the following inputs: " + "sample-rate, window, lags, peaks and troughs.", + default=[]) + """ + metrics.add_argument("-rrf", "--respiratory-response-function", + dest="metrics", + action="append_const", + const=rrf, + help="Respiratory response function. Requires the following inputs: " + "sample-rate, oversampling, time-length, onset and tr.", + default=[]) + metrics.add_argument("-rcard", "--retroicor-card", + dest="metrics", + action="append_const", + const="r_card", + help="Computes regressors for cardiac signal. Requires the following " + "inputs: tr, nscans and n_harm.", + default=[]) + metrics.add_argument("-rresp", "--retroicor-resp", + dest="metrics", + action="append_const", + const="r_resp", + help="Computes regressors for respiratory signal. Requires the following " + "inputs: tr, nscans and n_harm.", + default=[]) + # Metric arguments metric_arg.add_argument("-sr", "--sample-rate", dest="sample_rate", type=float, - action=MetricsArgDict, help="Sampling rate of the physiological data in Hz.", - default=argparse.SUPPRESS) + default=None) metric_arg.add_argument("-pk", "--peaks", dest="peaks", type=str, - action=MetricsArgDict, help="Full path and filename of the list with the indexed peaks' " "positions of the physiological data.", - default=argparse.SUPPRESS) + default=None) metric_arg.add_argument("-tg", "--troughs", dest="troughs", type=str, - action=MetricsArgDict, help="Full path and filename of the list with the indexed troughs' " "positions of the physiological data.", - default=argparse.SUPPRESS) + default=None) metric_arg.add_argument("-os", "--oversampling", dest="oversampling", type=int, - action=MetricsArgDict, help="Temporal oversampling factor. " "Default is 50.", - default=argparse.SUPPRESS) + default=50) metric_arg.add_argument("-tl", "--time-length", dest="time_length", type=int, - action=MetricsArgDict, help="RRF or CRF Kernel length in seconds.", - default=argparse.SUPPRESS) + default=None) metric_arg.add_argument("-onset", "--onset", dest="onset", type=float, - action=MetricsArgDict, help="Onset of the response in seconds. " "Default is 0.", - default=argparse.SUPPRESS) + default=0) metric_arg.add_argument("-tr", "--tr", dest="tr", type=float, - action=MetricsArgDict, help="TR of sequence in seconds.", - default=argparse.SUPPRESS) + default=None) metric_arg.add_argument("-win", "--window", dest="window", type=int, - action=MetricsArgDict, help="Size of the sliding window in seconds. " "Default is 6 seconds.", - default=argparse.SUPPRESS) + default=6) metric_arg.add_argument("-lags", "--lags", dest="lags", nargs="*", type=int, - action=MetricsArgDict, help="List of lags to apply to the RV estimate " "in seconds.", - default=argparse.SUPPRESS) + default=None) metric_arg.add_argument("-nscans", "--number-scans", dest="nscans", type=int, - action=MetricsArgDict, - help="Number of scans. Default is 1.", - default=argparse.SUPPRESS) + help="Number of timepoints in the imaging data. " + "Also called sub-bricks, TRs, scans, volumes." + "Default is 1.", + default=1) metric_arg.add_argument("-nharm", "--number-harmonics", dest="n_harm", type=int, - action=MetricsArgDict, help="Number of harmonics.", - default=argparse.SUPPRESS) + default=None) + + # Other optional arguments optional.add_argument("-debug", "--debug", dest="debug", action="store_true", From 6d8d0e90c1965c8b34f9d4eec50f4e42dbda9d5e Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 25 Feb 2021 11:08:59 -0500 Subject: [PATCH 32/52] Refactor and document RETROICOR code a bit. --- phys2denoise/metrics/retroicor.py | 147 ++++++++++++++++++++---------- phys2denoise/references.py | 2 + 2 files changed, 101 insertions(+), 48 deletions(-) diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py index e63e929..8b9fa6e 100644 --- a/phys2denoise/metrics/retroicor.py +++ b/phys2denoise/metrics/retroicor.py @@ -1,87 +1,138 @@ -"""implementation of retroicor""" +"""These functions compute RETROICOR regressors (Glover et al. 2000).""" -import argparse as ap import numpy as np -import scipy as sp -import json -import string -import random -import matplotlib as mpl; mpl.use('TkAgg') +import matplotlib as mpl + +mpl.use("TkAgg") import matplotlib.pyplot as plt -from scipy.signal import argrelmax -""" -This function computes RETROICOR regressors (Glover et al. 2000) -""" +from .. import references +from ..due import due -def compute_card_phase(card_peaks_timings,slice_timings,nscans,TR): - """ +def compute_phase_card(card_peaks_timings, slice_timings, n_scans, t_r): + """Calculate cardiac phase from cardiac peaks. - This function creates cardiac phase from cardiac peaks. Assumes that timing of cardiac events are given in same units as slice timings, for example seconds. - - """ - nscans = np.shape(slice_timings) - phase_card = np.zeros(nscans) - for ii in range(nscans): + Parameters + ---------- + peaks : 1D array_like + slice_timings : 1D array_like + n_scans : int + t_r : float + + Returns + ------- + phase_card : array_like + Cardiac phase signal. + """ + n_scans = np.shape(slice_timings) + phase_card = np.zeros(n_scans) + for ii in range(n_scans): # find previous cardiac peaks - previous_card_peaks = np.asarray(np.nonzero(card_peaks_timings < slice_timings[ii])) + previous_card_peaks = np.asarray( + np.nonzero(card_peaks_timings < slice_timings[ii]) + ) if np.size(previous_card_peaks) == 0: t1 = 0 else: last_peak = previous_card_peaks[0][-1] t1 = card_peaks_timings[last_peak] - + # find posterior cardiac peaks next_card_peaks = np.asarray(np.nonzero(card_peaks_timings > slice_timings[ii])) if np.size(next_card_peaks) == 0: - t2 = nscans * TR + t2 = n_scans * t_r else: next_peak = next_card_peaks[0][0] t2 = card_peaks_timings[next_peak] - + # compute cardiac phase - phase_card[ii] = (2*np.math.pi*(slice_timings[ii] - t1))/(t2-t1) + phase_card[ii] = (2 * np.math.pi * (slice_timings[ii] - t1)) / (t2 - t1) return phase_card -def compute_resp_phase(resp,sampling_time): - """ - This function creates respiration phase from resp trace. - """ +def compute_phase_resp(resp, sampling_time): + """Calculate respiratory phase from respiratory signal. + Parameters + ---------- + resp + sampling_time -def compute_retroicor_regressors(physio,TR,nscans,slice_timings,n_harmonics,card=FALSE,resp=FALSE): - nslices = np.shape(slice_timings) # number of slices + Returns + ------- + phase_resp : array_like + Respiratory phase signal. + """ + pass + + +@due.dcite(references.GLOVER_2000) +def compute_retroicor_regressors( + physio, t_r, n_scans, slice_timings, n_harmonics, card=False, resp=False +): + """Compute RETROICOR regressors. + + Parameters + ---------- + physio : array_like + 1D array, whether cardiac or respiratory signal. + t_r : float + n_scans : int + slice_timings + n_harmonics : int + card : bool, optional + resp : bool, optional + + Returns + ------- + retroicor_regressors : list + phase : array_like + 2D array of shape (n_scans, n_slices) + + References + ---------- + * Glover, G. H., Li, T. Q., & Ress, D. (2000). + Image‐based method for retrospective correction of physiological + motion effects in fMRI: RETROICOR. + Magnetic Resonance in Medicine: + An Official Journal of the International Society for Magnetic Resonance in Medicine, + 44(1), 162-167. + """ + n_slices = np.shape(slice_timings) # number of slices # if respiration, compute histogram and temporal derivative of respiration signal if resp: - resp_hist, resp_hist_bins = plt.hist(physio,bins=100) - resp_diff = np.diff(physio,n=1) - - #initialize output variables + # TODO: Replace with numpy.histogram + resp_hist, resp_hist_bins = plt.hist(physio, bins=100) + resp_diff = np.diff(physio, n=1) + + # initialize output variables retroicor_regressors = [] - phase = np.empty((nscans,nslices)) + phase = np.empty((n_scans, n_slices)) - for jj in range(nslices): + for i_slice in range(n_slices): # Initialize slice timings for current slice - crslice_timings = TR * np.arange(nscans)+slice_timings[jj] + crslice_timings = t_r * np.arange(n_scans) + slice_timings[i_slice] # Compute physiological phases using the timings of physio events (e.g. peaks) slice sampling times if card: - phase[,jj] = compute_phase_card(physio,crslice_timings) + phase[:, i_slice] = compute_phase_card(physio, crslice_timings) if resp: - phase[,jj] = compute_phase_resp(resp_diff,resp_hist,resp_hist_bins,crslice_timings) - # Compute retroicor regressors - for nn in range(n_harmonics): - retroicor_regressors[jj][:,2*nn] = np.cos((nn+1)*phase[jj]) - retricor_regressor[jj][:,2*nn+1] = np.sin((nn+1)*phase[jj]) - - return retroicor_regressors,phase - - - + phase[:, i_slice] = compute_phase_resp( + resp_diff, resp_hist, resp_hist_bins, crslice_timings + ) + # Compute retroicor regressors + for j_harm in range(n_harmonics): + retroicor_regressors[i_slice][:, 2 * j_harm] = np.cos( + (j_harm + 1) * phase[i_slice] + ) + retroicor_regressors[i_slice][:, 2 * j_harm + 1] = np.sin( + (j_harm + 1) * phase[i_slice] + ) + + return retroicor_regressors, phase diff --git a/phys2denoise/references.py b/phys2denoise/references.py index 297ba3b..eeb914a 100644 --- a/phys2denoise/references.py +++ b/phys2denoise/references.py @@ -6,3 +6,5 @@ POWER_2018 = Doi("10.1073/pnas.1720985115") POWER_2020 = Doi("10.1016/j.neuroimage.2019.116234") + +GLOVER_2000 = Doi("10.1002/1522-2594(200007)44:1<162::AID-MRM23>3.0.CO;2-E") From 5adc4f50a9d8ae479363bb7d9b86bfd93fee1bf9 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 17:33:08 +0100 Subject: [PATCH 33/52] Change call to retroicor into argument selection of type --- phys2denoise/cli/run.py | 40 +++++++++++++----------------------- phys2denoise/phys2denoise.py | 15 ++------------ 2 files changed, 16 insertions(+), 39 deletions(-) diff --git a/phys2denoise/cli/run.py b/phys2denoise/cli/run.py index b6a7248..eaea2d1 100644 --- a/phys2denoise/cli/run.py +++ b/phys2denoise/cli/run.py @@ -9,22 +9,6 @@ from phys2denoise.metrics.chest_belt import rpv, rv, rvt, rrf, env -class MetricsArgDict(argparse.Action): - """ - Custom Argparse Action to create a dictionary with the metrics' arguments in parser's output. - - """ - def __call__(self, parser, namespace, values, option_strings): - if not hasattr(namespace, "metrics_arg"): - setattr(namespace, "metrics_arg", dict()) - Keys = ["sample_rate", "peaks", "throughs", "oversampling", "time_length", "onset", - "tr", "window", "lags", "nscans", "nharm"] - Vals = ["None", "None", "None", "50", "None", "0", "None", "6", "None", "1", "None"] - for k, v in zip(Keys, Vals): - getattr(namespace, "metrics_arg")[k] = v - getattr(namespace, "metrics_arg")[self.dest] = values - - def _get_parser(): """ Parse command line inputs for this function. @@ -39,7 +23,6 @@ def _get_parser(): # Argument parser follow template provided by RalphyZ. # https://stackoverflow.com/a/43456577 """ - parser = argparse.ArgumentParser() optional = parser._action_groups.pop() required = parser.add_argument_group("Required Argument") @@ -106,19 +89,12 @@ def _get_parser(): help="Respiratory response function. Requires the following inputs: " "sample-rate, oversampling, time-length, onset and tr.", default=[]) - metrics.add_argument("-rcard", "--retroicor-card", + metrics.add_argument("-rcor", "--retroicor", dest="metrics", action="append_const", const="r_card", help="Computes regressors for cardiac signal. Requires the following " - "inputs: tr, nscans and n_harm.", - default=[]) - metrics.add_argument("-rresp", "--retroicor-resp", - dest="metrics", - action="append_const", - const="r_resp", - help="Computes regressors for respiratory signal. Requires the following " - "inputs: tr, nscans and n_harm.", + "inputs: either card or resp, tr, nscans and n_harm.", default=[]) # Metric arguments metric_arg.add_argument("-sr", "--sample-rate", @@ -185,6 +161,18 @@ def _get_parser(): type=int, help="Number of harmonics.", default=None) + metric_arg.add_argument("-card", "--cardiac", + dest="card", + type=bool, + action="store_true", + help="Compute *cardiac* RETROICOR.", + default=False) + metric_arg.add_argument("-resp", "--resp", + dest="resp", + type=bool, + action="store_true", + help="Compute *respiratory* RETROICOR.", + default=False) # Other optional arguments optional.add_argument("-debug", "--debug", diff --git a/phys2denoise/phys2denoise.py b/phys2denoise/phys2denoise.py index b9cd65e..9a6a9b3 100644 --- a/phys2denoise/phys2denoise.py +++ b/phys2denoise/phys2denoise.py @@ -172,19 +172,8 @@ def phys2denoise(filename, outdir='.', # Goes through the list of metrics and calls them for metric in metrics: - if metric == 'retroicor_card': - args = select_input_args(compute_retroicor_regressors, kwargs) - args['card'] = True - regr['retroicor_card'] = compute_retroicor_regressors(physio, - **args) - elif metric == 'retroicor_resp': - args = select_input_args(compute_retroicor_regressors, kwargs) - args['resp'] = True - regr['retroicor_resp'] = compute_retroicor_regressors(physio, - **args) - else: - args = select_input_args(metric, kwargs) - regr[f'{metric}'] = metric(physio, **args) + args = select_input_args(metric, kwargs) + regr[f'{metric}'] = metric(physio, **args) #!# Add regressors visualisation From b98e19c245c2ea10f20140a43ef973a9df32f816 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 25 Feb 2021 11:37:10 -0500 Subject: [PATCH 34/52] Resolve missing arguments and some bugs. --- phys2denoise/metrics/retroicor.py | 105 +++++++++++++++++++++++------- 1 file changed, 81 insertions(+), 24 deletions(-) diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py index 8b9fa6e..deb57d9 100644 --- a/phys2denoise/metrics/retroicor.py +++ b/phys2denoise/metrics/retroicor.py @@ -10,7 +10,7 @@ from ..due import due -def compute_phase_card(card_peaks_timings, slice_timings, n_scans, t_r): +def compute_phase_card(peaks, slice_timings, n_scans, t_r): """Calculate cardiac phase from cardiac peaks. Assumes that timing of cardiac events are given in same units @@ -19,61 +19,108 @@ def compute_phase_card(card_peaks_timings, slice_timings, n_scans, t_r): Parameters ---------- peaks : 1D array_like + Cardiac peak times, in seconds. slice_timings : 1D array_like + Slice times, in seconds. n_scans : int + Number of volumes in the imaging run. t_r : float + Sampling rate of the imaging run, in seconds. Returns ------- phase_card : array_like Cardiac phase signal. """ - n_scans = np.shape(slice_timings) + n_slices = np.shape(slice_timings) phase_card = np.zeros(n_scans) - for ii in range(n_scans): + + for i_slice in range(n_slices): # find previous cardiac peaks - previous_card_peaks = np.asarray( - np.nonzero(card_peaks_timings < slice_timings[ii]) - ) + previous_card_peaks = np.asarray(np.nonzero(peaks < slice_timings[i_slice])) if np.size(previous_card_peaks) == 0: t1 = 0 else: last_peak = previous_card_peaks[0][-1] - t1 = card_peaks_timings[last_peak] + t1 = peaks[last_peak] # find posterior cardiac peaks - next_card_peaks = np.asarray(np.nonzero(card_peaks_timings > slice_timings[ii])) + next_card_peaks = np.asarray(np.nonzero(peaks > slice_timings[i_slice])) if np.size(next_card_peaks) == 0: t2 = n_scans * t_r else: next_peak = next_card_peaks[0][0] - t2 = card_peaks_timings[next_peak] + t2 = peaks[next_peak] # compute cardiac phase - phase_card[ii] = (2 * np.math.pi * (slice_timings[ii] - t1)) / (t2 - t1) + phase_card[i_slice] = (2 * np.math.pi * (slice_timings[i_slice] - t1)) / ( + t2 - t1 + ) return phase_card -def compute_phase_resp(resp, sampling_time): +def compute_phase_resp(resp, sample_rate, n_scans, slice_timings, t_r): """Calculate respiratory phase from respiratory signal. Parameters ---------- - resp - sampling_time + resp : 1D array_like + Respiratory signal. + sample_rate : float + Sample rate of physio, in Hertz. + n_scans + Number of volumes in the imaging run. + slice_timings + Slice times, in seconds. + t_r + Sample rate of the imaging run, in seconds. Returns ------- phase_resp : array_like Respiratory phase signal. """ - pass + n_slices = np.shape(slice_timings) + phase_resp = np.zeros((n_scans, n_slices)) + + # generate histogram from respiratory signal + # TODO: Replace with numpy.histogram + resp_hist, resp_hist_bins = plt.hist(resp, bins=100) + + # first compute derivative of respiration signal + resp_diff = np.diff(resp, n=1) + + for i_slice in range(n_slices): + # generate slice acquisition timings across all scans + times_crSlice = t_r * np.arange(n_scans) + slice_timings[i_slice] + phase_resp_crSlice = np.zeros(n_scans) + for j_scan in range(n_scans): + iphys = int( + max([1, round(times_crSlice[j_scan] * sample_rate)]) + ) # closest idx in resp waveform + iphys = min([iphys, len(resp_diff)]) # cannot be longer than resp_diff + thisBin = np.argmin(abs(resp[iphys] - resp_hist_bins)) + numerator = np.sum(resp_hist[0:thisBin]) + phase_resp_crSlice[j_scan] = ( + np.math.pi * np.sign(resp_diff[iphys]) * (numerator / len(resp)) + ) + + phase_resp[:, i_slice] = phase_resp_crSlice + + return phase_resp @due.dcite(references.GLOVER_2000) def compute_retroicor_regressors( - physio, t_r, n_scans, slice_timings, n_harmonics, card=False, resp=False + physio, + sample_rate, + t_r, + n_scans, + slice_timings, + n_harmonics, + card=False, + resp=False, ): """Compute RETROICOR regressors. @@ -81,12 +128,22 @@ def compute_retroicor_regressors( ---------- physio : array_like 1D array, whether cardiac or respiratory signal. + If cardiac, the array is a set of peaks in seconds. + If respiratory, the array is the actual respiratory signal. + sample_rate : float + Physio sample rate, in Hertz. t_r : float + Imaging sample rate, in seconds. n_scans : int - slice_timings + Number of volumes in the imaging run. + slice_timings : array_like + Slice times, in seconds. n_harmonics : int + ??? card : bool, optional + Whether the physio data correspond to cardiac or repiratory signal. resp : bool, optional + Whether the physio data correspond to cardiac or repiratory signal. Returns ------- @@ -105,12 +162,6 @@ def compute_retroicor_regressors( """ n_slices = np.shape(slice_timings) # number of slices - # if respiration, compute histogram and temporal derivative of respiration signal - if resp: - # TODO: Replace with numpy.histogram - resp_hist, resp_hist_bins = plt.hist(physio, bins=100) - resp_diff = np.diff(physio, n=1) - # initialize output variables retroicor_regressors = [] phase = np.empty((n_scans, n_slices)) @@ -118,12 +169,18 @@ def compute_retroicor_regressors( for i_slice in range(n_slices): # Initialize slice timings for current slice crslice_timings = t_r * np.arange(n_scans) + slice_timings[i_slice] - # Compute physiological phases using the timings of physio events (e.g. peaks) slice sampling times + + # Compute physiological phases using the timings of physio events (e.g. peaks) + # slice sampling times if card: phase[:, i_slice] = compute_phase_card(physio, crslice_timings) if resp: phase[:, i_slice] = compute_phase_resp( - resp_diff, resp_hist, resp_hist_bins, crslice_timings + physio, + sample_rate, + n_scans, + slice_timings, + t_r, ) # Compute retroicor regressors From 58b908c9597c92e325ba60859e4a7e32dfc8fd32 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 25 Feb 2021 11:37:15 -0500 Subject: [PATCH 35/52] Rename function. --- phys2denoise/metrics/retroicor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py index deb57d9..a99d3cc 100644 --- a/phys2denoise/metrics/retroicor.py +++ b/phys2denoise/metrics/retroicor.py @@ -112,7 +112,7 @@ def compute_phase_resp(resp, sample_rate, n_scans, slice_timings, t_r): @due.dcite(references.GLOVER_2000) -def compute_retroicor_regressors( +def retroicor( physio, sample_rate, t_r, From d6e198469055839b00a6c6f76ba04b374b1a33dd Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 25 Feb 2021 11:51:07 -0500 Subject: [PATCH 36/52] Fix function call. --- phys2denoise/metrics/retroicor.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py index a99d3cc..7072501 100644 --- a/phys2denoise/metrics/retroicor.py +++ b/phys2denoise/metrics/retroicor.py @@ -173,7 +173,13 @@ def retroicor( # Compute physiological phases using the timings of physio events (e.g. peaks) # slice sampling times if card: - phase[:, i_slice] = compute_phase_card(physio, crslice_timings) + phase[:, i_slice] = compute_phase_card( + physio, + crslice_timings, + n_scans, + t_r, + ) + if resp: phase[:, i_slice] = compute_phase_resp( physio, From 6e7f7ee1d496056142ced2030a38d63d4430befd Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 18:21:04 +0100 Subject: [PATCH 37/52] Revert "Change call to retroicor into argument selection of type" This reverts commit 5adc4f50a9d8ae479363bb7d9b86bfd93fee1bf9. --- phys2denoise/cli/run.py | 40 +++++++++++++++++++++++------------- phys2denoise/phys2denoise.py | 15 ++++++++++++-- 2 files changed, 39 insertions(+), 16 deletions(-) diff --git a/phys2denoise/cli/run.py b/phys2denoise/cli/run.py index eaea2d1..b6a7248 100644 --- a/phys2denoise/cli/run.py +++ b/phys2denoise/cli/run.py @@ -9,6 +9,22 @@ from phys2denoise.metrics.chest_belt import rpv, rv, rvt, rrf, env +class MetricsArgDict(argparse.Action): + """ + Custom Argparse Action to create a dictionary with the metrics' arguments in parser's output. + + """ + def __call__(self, parser, namespace, values, option_strings): + if not hasattr(namespace, "metrics_arg"): + setattr(namespace, "metrics_arg", dict()) + Keys = ["sample_rate", "peaks", "throughs", "oversampling", "time_length", "onset", + "tr", "window", "lags", "nscans", "nharm"] + Vals = ["None", "None", "None", "50", "None", "0", "None", "6", "None", "1", "None"] + for k, v in zip(Keys, Vals): + getattr(namespace, "metrics_arg")[k] = v + getattr(namespace, "metrics_arg")[self.dest] = values + + def _get_parser(): """ Parse command line inputs for this function. @@ -23,6 +39,7 @@ def _get_parser(): # Argument parser follow template provided by RalphyZ. # https://stackoverflow.com/a/43456577 """ + parser = argparse.ArgumentParser() optional = parser._action_groups.pop() required = parser.add_argument_group("Required Argument") @@ -89,12 +106,19 @@ def _get_parser(): help="Respiratory response function. Requires the following inputs: " "sample-rate, oversampling, time-length, onset and tr.", default=[]) - metrics.add_argument("-rcor", "--retroicor", + metrics.add_argument("-rcard", "--retroicor-card", dest="metrics", action="append_const", const="r_card", help="Computes regressors for cardiac signal. Requires the following " - "inputs: either card or resp, tr, nscans and n_harm.", + "inputs: tr, nscans and n_harm.", + default=[]) + metrics.add_argument("-rresp", "--retroicor-resp", + dest="metrics", + action="append_const", + const="r_resp", + help="Computes regressors for respiratory signal. Requires the following " + "inputs: tr, nscans and n_harm.", default=[]) # Metric arguments metric_arg.add_argument("-sr", "--sample-rate", @@ -161,18 +185,6 @@ def _get_parser(): type=int, help="Number of harmonics.", default=None) - metric_arg.add_argument("-card", "--cardiac", - dest="card", - type=bool, - action="store_true", - help="Compute *cardiac* RETROICOR.", - default=False) - metric_arg.add_argument("-resp", "--resp", - dest="resp", - type=bool, - action="store_true", - help="Compute *respiratory* RETROICOR.", - default=False) # Other optional arguments optional.add_argument("-debug", "--debug", diff --git a/phys2denoise/phys2denoise.py b/phys2denoise/phys2denoise.py index 9a6a9b3..b9cd65e 100644 --- a/phys2denoise/phys2denoise.py +++ b/phys2denoise/phys2denoise.py @@ -172,8 +172,19 @@ def phys2denoise(filename, outdir='.', # Goes through the list of metrics and calls them for metric in metrics: - args = select_input_args(metric, kwargs) - regr[f'{metric}'] = metric(physio, **args) + if metric == 'retroicor_card': + args = select_input_args(compute_retroicor_regressors, kwargs) + args['card'] = True + regr['retroicor_card'] = compute_retroicor_regressors(physio, + **args) + elif metric == 'retroicor_resp': + args = select_input_args(compute_retroicor_regressors, kwargs) + args['resp'] = True + regr['retroicor_resp'] = compute_retroicor_regressors(physio, + **args) + else: + args = select_input_args(metric, kwargs) + regr[f'{metric}'] = metric(physio, **args) #!# Add regressors visualisation From 76fe4371c4b33a58928a60f792d7ee316bb70ea4 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 18:31:14 +0100 Subject: [PATCH 38/52] Better retroicor call --- phys2denoise/phys2denoise.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/phys2denoise/phys2denoise.py b/phys2denoise/phys2denoise.py index b9cd65e..284b36c 100644 --- a/phys2denoise/phys2denoise.py +++ b/phys2denoise/phys2denoise.py @@ -175,16 +175,18 @@ def phys2denoise(filename, outdir='.', if metric == 'retroicor_card': args = select_input_args(compute_retroicor_regressors, kwargs) args['card'] = True - regr['retroicor_card'] = compute_retroicor_regressors(physio, - **args) + retroicor_regrs = compute_retroicor_regressors(physio, **args) + for vslice in range(len(args['slice_timings'])): + regr[f'retroicor_card_slice-{vslice}'] = retroicor_regrs[vslice] elif metric == 'retroicor_resp': args = select_input_args(compute_retroicor_regressors, kwargs) args['resp'] = True - regr['retroicor_resp'] = compute_retroicor_regressors(physio, - **args) + retroicor_regrs = compute_retroicor_regressors(physio, **args) + for vslice in range(len(args['slice_timings'])): + regr[f'retroicor_resp_slice-{vslice}'] = retroicor_regrs[vslice] else: args = select_input_args(metric, kwargs) - regr[f'{metric}'] = metric(physio, **args) + regr[metric.__name__] = metric(physio, **args) #!# Add regressors visualisation From 58f41ba7b58b1c0a6364bd80680e5cab930d65d6 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 25 Feb 2021 12:33:16 -0500 Subject: [PATCH 39/52] Move phase functions to appropriate files. --- phys2denoise/metrics/cardiac.py | 50 +++++++++++++ phys2denoise/metrics/chest_belt.py | 55 ++++++++++++++ phys2denoise/metrics/retroicor.py | 111 ++--------------------------- 3 files changed, 109 insertions(+), 107 deletions(-) diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py index f295848..43c8079 100644 --- a/phys2denoise/metrics/cardiac.py +++ b/phys2denoise/metrics/cardiac.py @@ -59,3 +59,53 @@ def _crf(t): crf_arr = _crf(time_stamps) crf_arr = crf_arr / max(abs(crf_arr)) return crf_arr + + +def cardiac_phase(peaks, slice_timings, n_scans, t_r): + """Calculate cardiac phase from cardiac peaks. + + Assumes that timing of cardiac events are given in same units + as slice timings, for example seconds. + + Parameters + ---------- + peaks : 1D array_like + Cardiac peak times, in seconds. + slice_timings : 1D array_like + Slice times, in seconds. + n_scans : int + Number of volumes in the imaging run. + t_r : float + Sampling rate of the imaging run, in seconds. + + Returns + ------- + phase_card : array_like + Cardiac phase signal. + """ + n_slices = np.shape(slice_timings) + phase_card = np.zeros(n_scans) + + for i_slice in range(n_slices): + # find previous cardiac peaks + previous_card_peaks = np.asarray(np.nonzero(peaks < slice_timings[i_slice])) + if np.size(previous_card_peaks) == 0: + t1 = 0 + else: + last_peak = previous_card_peaks[0][-1] + t1 = peaks[last_peak] + + # find posterior cardiac peaks + next_card_peaks = np.asarray(np.nonzero(peaks > slice_timings[i_slice])) + if np.size(next_card_peaks) == 0: + t2 = n_scans * t_r + else: + next_peak = next_card_peaks[0][0] + t2 = peaks[next_peak] + + # compute cardiac phase + phase_card[i_slice] = (2 * np.math.pi * (slice_timings[i_slice] - t1)) / ( + t2 - t1 + ) + + return phase_card diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py index 9ff8225..7c5ae1c 100644 --- a/phys2denoise/metrics/chest_belt.py +++ b/phys2denoise/metrics/chest_belt.py @@ -1,10 +1,14 @@ """Denoising metrics for chest belt recordings.""" +import matplotlib as mpl import numpy as np import pandas as pd from scipy.ndimage.filters import convolve1d from scipy.signal import detrend, resample from scipy.stats import zscore +mpl.use("TkAgg") +import matplotlib.pyplot as plt + from .. import references from ..due import due from . import utils @@ -215,3 +219,54 @@ def _rrf(t): rrf_arr = _rrf(time_stamps) rrf_arr = rrf_arr / max(abs(rrf_arr)) return rrf_arr + + +def respiratory_phase(resp, sample_rate, n_scans, slice_timings, t_r): + """Calculate respiratory phase from respiratory signal. + + Parameters + ---------- + resp : 1D array_like + Respiratory signal. + sample_rate : float + Sample rate of physio, in Hertz. + n_scans + Number of volumes in the imaging run. + slice_timings + Slice times, in seconds. + t_r + Sample rate of the imaging run, in seconds. + + Returns + ------- + phase_resp : array_like + Respiratory phase signal. + """ + n_slices = np.shape(slice_timings) + phase_resp = np.zeros((n_scans, n_slices)) + + # generate histogram from respiratory signal + # TODO: Replace with numpy.histogram + resp_hist, resp_hist_bins = plt.hist(resp, bins=100) + + # first compute derivative of respiration signal + resp_diff = np.diff(resp, n=1) + + for i_slice in range(n_slices): + # generate slice acquisition timings across all scans + times_crSlice = t_r * np.arange(n_scans) + slice_timings[i_slice] + phase_resp_crSlice = np.zeros(n_scans) + for j_scan in range(n_scans): + iphys = int( + max([1, round(times_crSlice[j_scan] * sample_rate)]) + ) # closest idx in resp waveform + iphys = min([iphys, len(resp_diff)]) # cannot be longer than resp_diff + thisBin = np.argmin(abs(resp[iphys] - resp_hist_bins)) + numerator = np.sum(resp_hist[0:thisBin]) + phase_resp_crSlice[j_scan] = ( + np.math.pi * np.sign(resp_diff[iphys]) * (numerator / len(resp)) + ) + + phase_resp[:, i_slice] = phase_resp_crSlice + + return phase_resp diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py index 7072501..690b797 100644 --- a/phys2denoise/metrics/retroicor.py +++ b/phys2denoise/metrics/retroicor.py @@ -1,116 +1,13 @@ """These functions compute RETROICOR regressors (Glover et al. 2000).""" import numpy as np -import matplotlib as mpl - -mpl.use("TkAgg") -import matplotlib.pyplot as plt +from .cardiac import cardiac_phase +from .chest_belt import respiratory_phase from .. import references from ..due import due -def compute_phase_card(peaks, slice_timings, n_scans, t_r): - """Calculate cardiac phase from cardiac peaks. - - Assumes that timing of cardiac events are given in same units - as slice timings, for example seconds. - - Parameters - ---------- - peaks : 1D array_like - Cardiac peak times, in seconds. - slice_timings : 1D array_like - Slice times, in seconds. - n_scans : int - Number of volumes in the imaging run. - t_r : float - Sampling rate of the imaging run, in seconds. - - Returns - ------- - phase_card : array_like - Cardiac phase signal. - """ - n_slices = np.shape(slice_timings) - phase_card = np.zeros(n_scans) - - for i_slice in range(n_slices): - # find previous cardiac peaks - previous_card_peaks = np.asarray(np.nonzero(peaks < slice_timings[i_slice])) - if np.size(previous_card_peaks) == 0: - t1 = 0 - else: - last_peak = previous_card_peaks[0][-1] - t1 = peaks[last_peak] - - # find posterior cardiac peaks - next_card_peaks = np.asarray(np.nonzero(peaks > slice_timings[i_slice])) - if np.size(next_card_peaks) == 0: - t2 = n_scans * t_r - else: - next_peak = next_card_peaks[0][0] - t2 = peaks[next_peak] - - # compute cardiac phase - phase_card[i_slice] = (2 * np.math.pi * (slice_timings[i_slice] - t1)) / ( - t2 - t1 - ) - - return phase_card - - -def compute_phase_resp(resp, sample_rate, n_scans, slice_timings, t_r): - """Calculate respiratory phase from respiratory signal. - - Parameters - ---------- - resp : 1D array_like - Respiratory signal. - sample_rate : float - Sample rate of physio, in Hertz. - n_scans - Number of volumes in the imaging run. - slice_timings - Slice times, in seconds. - t_r - Sample rate of the imaging run, in seconds. - - Returns - ------- - phase_resp : array_like - Respiratory phase signal. - """ - n_slices = np.shape(slice_timings) - phase_resp = np.zeros((n_scans, n_slices)) - - # generate histogram from respiratory signal - # TODO: Replace with numpy.histogram - resp_hist, resp_hist_bins = plt.hist(resp, bins=100) - - # first compute derivative of respiration signal - resp_diff = np.diff(resp, n=1) - - for i_slice in range(n_slices): - # generate slice acquisition timings across all scans - times_crSlice = t_r * np.arange(n_scans) + slice_timings[i_slice] - phase_resp_crSlice = np.zeros(n_scans) - for j_scan in range(n_scans): - iphys = int( - max([1, round(times_crSlice[j_scan] * sample_rate)]) - ) # closest idx in resp waveform - iphys = min([iphys, len(resp_diff)]) # cannot be longer than resp_diff - thisBin = np.argmin(abs(resp[iphys] - resp_hist_bins)) - numerator = np.sum(resp_hist[0:thisBin]) - phase_resp_crSlice[j_scan] = ( - np.math.pi * np.sign(resp_diff[iphys]) * (numerator / len(resp)) - ) - - phase_resp[:, i_slice] = phase_resp_crSlice - - return phase_resp - - @due.dcite(references.GLOVER_2000) def retroicor( physio, @@ -173,7 +70,7 @@ def retroicor( # Compute physiological phases using the timings of physio events (e.g. peaks) # slice sampling times if card: - phase[:, i_slice] = compute_phase_card( + phase[:, i_slice] = cardiac_phase( physio, crslice_timings, n_scans, @@ -181,7 +78,7 @@ def retroicor( ) if resp: - phase[:, i_slice] = compute_phase_resp( + phase[:, i_slice] = respiratory_phase( physio, sample_rate, n_scans, From e0b4a1166228ba2d8c193f0a6e825fe5288bcd73 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 25 Feb 2021 12:33:25 -0500 Subject: [PATCH 40/52] Fill out init file. --- phys2denoise/metrics/__init__.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/phys2denoise/metrics/__init__.py b/phys2denoise/metrics/__init__.py index e69de29..1337dec 100644 --- a/phys2denoise/metrics/__init__.py +++ b/phys2denoise/metrics/__init__.py @@ -0,0 +1,15 @@ +"""Metrics derived from physiological signals.""" +from .cardiac import crf, cardiac_phase +from .chest_belt import rpv, env, rv, rrf, respiratory_phase +from .retroicor import retroicor + +__all__ = [ + "crf", + "cardiac_phase", + "rpv", + "env", + "rv", + "rrf", + "respiratory_phase", + "retroicor", +] From 6508a5efaee6ceaea3f773fba80d59a042230a74 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 25 Feb 2021 12:46:32 -0500 Subject: [PATCH 41/52] Preallocate RETROICOR contents. --- phys2denoise/metrics/retroicor.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py index 690b797..f4ae089 100644 --- a/phys2denoise/metrics/retroicor.py +++ b/phys2denoise/metrics/retroicor.py @@ -64,6 +64,8 @@ def retroicor( phase = np.empty((n_scans, n_slices)) for i_slice in range(n_slices): + retroicor_regressors.append(np.empty((n_scans, 2 * n_harmonics))) + # Initialize slice timings for current slice crslice_timings = t_r * np.arange(n_scans) + slice_timings[i_slice] From d76735930806d5aef617891bcf73343a38fd9c35 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 18:49:26 +0100 Subject: [PATCH 42/52] Deal better with retroicor outputs --- phys2denoise/phys2denoise.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/phys2denoise/phys2denoise.py b/phys2denoise/phys2denoise.py index 284b36c..a5ff33e 100644 --- a/phys2denoise/phys2denoise.py +++ b/phys2denoise/phys2denoise.py @@ -177,13 +177,19 @@ def phys2denoise(filename, outdir='.', args['card'] = True retroicor_regrs = compute_retroicor_regressors(physio, **args) for vslice in range(len(args['slice_timings'])): - regr[f'retroicor_card_slice-{vslice}'] = retroicor_regrs[vslice] + for harm in range(args['n_harm']): + key = f'rcor-card_s-{vslice}_hrm-{harm}' + regr[f'{key}_cos'] = retroicor_regrs[vslice][:, harm*2] + regr[f'{key}_sin'] = retroicor_regrs[vslice][:, harm*2+1] elif metric == 'retroicor_resp': args = select_input_args(compute_retroicor_regressors, kwargs) args['resp'] = True retroicor_regrs = compute_retroicor_regressors(physio, **args) for vslice in range(len(args['slice_timings'])): - regr[f'retroicor_resp_slice-{vslice}'] = retroicor_regrs[vslice] + for harm in range(args['n_harm']): + key = f'rcor-resp_s-{vslice}_hrm-{harm}' + regr[f'{key}_cos'] = retroicor_regrs[vslice][:, harm*2] + regr[f'{key}_sin'] = retroicor_regrs[vslice][:, harm*2+1] else: args = select_input_args(metric, kwargs) regr[metric.__name__] = metric(physio, **args) From 1ca87800c22cf5666329f028b964dd4047ce0c13 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 18:51:39 +0100 Subject: [PATCH 43/52] Remove unused code --- phys2denoise/cli/run.py | 16 ---------------- phys2denoise/phys2denoise.py | 10 +++++----- 2 files changed, 5 insertions(+), 21 deletions(-) diff --git a/phys2denoise/cli/run.py b/phys2denoise/cli/run.py index b6a7248..4103513 100644 --- a/phys2denoise/cli/run.py +++ b/phys2denoise/cli/run.py @@ -9,22 +9,6 @@ from phys2denoise.metrics.chest_belt import rpv, rv, rvt, rrf, env -class MetricsArgDict(argparse.Action): - """ - Custom Argparse Action to create a dictionary with the metrics' arguments in parser's output. - - """ - def __call__(self, parser, namespace, values, option_strings): - if not hasattr(namespace, "metrics_arg"): - setattr(namespace, "metrics_arg", dict()) - Keys = ["sample_rate", "peaks", "throughs", "oversampling", "time_length", "onset", - "tr", "window", "lags", "nscans", "nharm"] - Vals = ["None", "None", "None", "50", "None", "0", "None", "6", "None", "1", "None"] - for k, v in zip(Keys, Vals): - getattr(namespace, "metrics_arg")[k] = v - getattr(namespace, "metrics_arg")[self.dest] = values - - def _get_parser(): """ Parse command line inputs for this function. diff --git a/phys2denoise/phys2denoise.py b/phys2denoise/phys2denoise.py index a5ff33e..0cc681e 100644 --- a/phys2denoise/phys2denoise.py +++ b/phys2denoise/phys2denoise.py @@ -22,7 +22,7 @@ from phys2denoise.cli.run import _get_parser from phys2denoise.metrics.cardiac import crf from phys2denoise.metrics.chest_belt import rpv, rv, rvt, rrf -from phys2denoise.metrics.retroicor import compute_retroicor_regressors +from phys2denoise.metrics.retroicor import retroicor from . import __version__ from .due import due, Doi @@ -173,18 +173,18 @@ def phys2denoise(filename, outdir='.', # Goes through the list of metrics and calls them for metric in metrics: if metric == 'retroicor_card': - args = select_input_args(compute_retroicor_regressors, kwargs) + args = select_input_args(retroicor, kwargs) args['card'] = True - retroicor_regrs = compute_retroicor_regressors(physio, **args) + retroicor_regrs = retroicor(physio, **args) for vslice in range(len(args['slice_timings'])): for harm in range(args['n_harm']): key = f'rcor-card_s-{vslice}_hrm-{harm}' regr[f'{key}_cos'] = retroicor_regrs[vslice][:, harm*2] regr[f'{key}_sin'] = retroicor_regrs[vslice][:, harm*2+1] elif metric == 'retroicor_resp': - args = select_input_args(compute_retroicor_regressors, kwargs) + args = select_input_args(retroicor, kwargs) args['resp'] = True - retroicor_regrs = compute_retroicor_regressors(physio, **args) + retroicor_regrs = retroicor(physio, **args) for vslice in range(len(args['slice_timings'])): for harm in range(args['n_harm']): key = f'rcor-resp_s-{vslice}_hrm-{harm}' From 8ad11c1f4dde0fd13948b08428971682c6fc09d5 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 17:54:40 +0000 Subject: [PATCH 44/52] Update CHANGELOG.md [skip ci] --- CHANGELOG.md | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 CHANGELOG.md diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..4bc5d0c --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,33 @@ +# 0.1.0 (Thu Feb 25 2021) + +:tada: This release contains work from new contributors! :tada: + +Thanks for all your work! + +:heart: Stefano Moia ([@smoia](https://github.com/smoia)) + +:heart: Inés Chavarría ([@ineschh](https://github.com/ineschh)) + +#### 💥 Breaking Change during development + +- Add main workflow of phys2denoise [#17](https://github.com/physiopy/phys2denoise/pull/17) ([@smoia](https://github.com/smoia) [@ineschh](https://github.com/ineschh)) + +#### ⚠️ Pushed to `master` + +- update all-contributors ([@smoia](https://github.com/smoia)) +- Correct attributes ([@smoia](https://github.com/smoia)) +- Add duecredit to setup ([@smoia](https://github.com/smoia)) +- General configuration of the infrastructure ([@smoia](https://github.com/smoia)) +- Correct development status ([@smoia](https://github.com/smoia)) +- Add infra files ([@smoia](https://github.com/smoia)) +- Initial commit ([@smoia](https://github.com/smoia)) + +#### 🏠 Internal + +- Update PR template to match `phys2bids` [#20](https://github.com/physiopy/phys2denoise/pull/20) ([@smoia](https://github.com/smoia)) +- Add initial infrastructure [#12](https://github.com/physiopy/phys2denoise/pull/12) ([@smoia](https://github.com/smoia)) + +#### Authors: 2 + +- Inés Chavarría ([@ineschh](https://github.com/ineschh)) +- Stefano Moia ([@smoia](https://github.com/smoia)) From d5f51485007e9f126ed3bbca1a5b3dad3b76b9df Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 17:55:23 +0000 Subject: [PATCH 45/52] Update contributors [skip ci] --- .all-contributorsrc | 9 +++++++++ README.md | 1 + 2 files changed, 10 insertions(+) diff --git a/.all-contributorsrc b/.all-contributorsrc index de5f715..3401917 100644 --- a/.all-contributorsrc +++ b/.all-contributorsrc @@ -47,6 +47,15 @@ "ideas", "review" ] + }, + { + "login": "ineschh", + "name": "Inés Chavarría", + "avatar_url": "https://avatars.githubusercontent.com/u/72545702?v=4", + "profile": "https://github.com/ineschh", + "contributions": [ + "infra" + ] } ], "contributorsPerLine": 8, diff --git a/README.md b/README.md index 8742d16..447728a 100644 --- a/README.md +++ b/README.md @@ -31,6 +31,7 @@ Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/d
Ross Markello
🚇
Stefano Moia
🔣 🚇 📆 💻 🤔
Taylor Salo
💻 🤔 👀 +
Inés Chavarría
🚇 From e729bc386bc329fc4d6d3bbe5ab25d347c2155ca Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 25 Feb 2021 13:09:38 -0500 Subject: [PATCH 46/52] Reorg retroicor. --- phys2denoise/metrics/__init__.py | 2 +- phys2denoise/metrics/{retroicor.py => multimodal.py} | 5 +++++ 2 files changed, 6 insertions(+), 1 deletion(-) rename phys2denoise/metrics/{retroicor.py => multimodal.py} (94%) diff --git a/phys2denoise/metrics/__init__.py b/phys2denoise/metrics/__init__.py index 1337dec..8eb43c8 100644 --- a/phys2denoise/metrics/__init__.py +++ b/phys2denoise/metrics/__init__.py @@ -1,7 +1,7 @@ """Metrics derived from physiological signals.""" from .cardiac import crf, cardiac_phase from .chest_belt import rpv, env, rv, rrf, respiratory_phase -from .retroicor import retroicor +from .multimodal import retroicor __all__ = [ "crf", diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/multimodal.py similarity index 94% rename from phys2denoise/metrics/retroicor.py rename to phys2denoise/metrics/multimodal.py index f4ae089..fb95eae 100644 --- a/phys2denoise/metrics/retroicor.py +++ b/phys2denoise/metrics/multimodal.py @@ -48,6 +48,11 @@ def retroicor( phase : array_like 2D array of shape (n_scans, n_slices) + Notes + ----- + RETROICOR regressors should be regressed from the imaging data *before* + any other preprocessing, including slice-timing correction and motion correction. + References ---------- * Glover, G. H., Li, T. Q., & Ress, D. (2000). From 58cb95a099cc049d766f81842df4823cae44c589 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 18:14:12 +0000 Subject: [PATCH 47/52] Update CHANGELOG.md [skip ci] --- CHANGELOG.md | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4bc5d0c..126ffa6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,17 @@ +# 0.2.0 (Thu Feb 25 2021) + +#### 💥 Breaking Change during development + +- Add initial set of metrics [#19](https://github.com/physiopy/phys2denoise/pull/19) ([@CesarCaballeroGaudes](https://github.com/CesarCaballeroGaudes) [@tsalo](https://github.com/tsalo) [@smoia](https://github.com/smoia)) + +#### Authors: 3 + +- Cesar Caballero Gaudes ([@CesarCaballeroGaudes](https://github.com/CesarCaballeroGaudes)) +- Stefano Moia ([@smoia](https://github.com/smoia)) +- Taylor Salo ([@tsalo](https://github.com/tsalo)) + +--- + # 0.1.0 (Thu Feb 25 2021) :tada: This release contains work from new contributors! :tada: From 63abe3791f2c4f6fe3423736f095686dd68227d1 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 18:14:47 +0000 Subject: [PATCH 48/52] Update contributors [skip ci] --- .all-contributorsrc | 3 ++- README.md | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/.all-contributorsrc b/.all-contributorsrc index 3401917..4e8b663 100644 --- a/.all-contributorsrc +++ b/.all-contributorsrc @@ -34,7 +34,8 @@ "infra", "projectManagement", "code", - "ideas" + "ideas", + "doc" ] }, { diff --git a/README.md b/README.md index 447728a..d6b7b14 100644 --- a/README.md +++ b/README.md @@ -29,7 +29,7 @@ Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/d
Cesar Caballero Gaudes
💻 🤔
Ross Markello
🚇 -
Stefano Moia
🔣 🚇 📆 💻 🤔 +
Stefano Moia
🔣 🚇 📆 💻 🤔 📖
Taylor Salo
💻 🤔 👀
Inés Chavarría
🚇 From cfced3b47e7e7e4254e775cdfbc448e4c7cb15b0 Mon Sep 17 00:00:00 2001 From: Stefano Moia Date: Thu, 25 Feb 2021 19:21:25 +0100 Subject: [PATCH 49/52] Fix python version call [skip ci] --- .github/workflows/pythonpublish.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/pythonpublish.yml b/.github/workflows/pythonpublish.yml index 9e24640..4ddc1bd 100644 --- a/.github/workflows/pythonpublish.yml +++ b/.github/workflows/pythonpublish.yml @@ -18,7 +18,7 @@ jobs: - name: Set up Python uses: actions/setup-python@v1 with: - python-version: 3.7 + python-version: '3.7' - name: Install dependencies run: | python -m pip install --upgrade pip @@ -29,4 +29,4 @@ jobs: TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }} run: | python -m python setup.py sdist bdist_wheel - python -m twine upload dist/* \ No newline at end of file + python -m twine upload dist/* From 82c3022dbfc5f003098d31ccf6c0228fa19840b5 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 25 Feb 2021 13:49:47 -0500 Subject: [PATCH 50/52] Add some smoke tests and fix bugs. --- phys2denoise/metrics/cardiac.py | 56 +++++++++++-------- phys2denoise/metrics/chest_belt.py | 7 ++- phys2denoise/tests/test_metrics_cardiac.py | 41 ++++++++++++++ phys2denoise/tests/test_metrics_chest_belt.py | 42 ++++++++++++++ 4 files changed, 119 insertions(+), 27 deletions(-) create mode 100644 phys2denoise/tests/test_metrics_cardiac.py create mode 100644 phys2denoise/tests/test_metrics_chest_belt.py diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py index 43c8079..f52d7b2 100644 --- a/phys2denoise/metrics/cardiac.py +++ b/phys2denoise/metrics/cardiac.py @@ -61,7 +61,7 @@ def _crf(t): return crf_arr -def cardiac_phase(peaks, slice_timings, n_scans, t_r): +def cardiac_phase(peaks, sample_rate, slice_timings, n_scans, t_r): """Calculate cardiac phase from cardiac peaks. Assumes that timing of cardiac events are given in same units @@ -71,6 +71,8 @@ def cardiac_phase(peaks, slice_timings, n_scans, t_r): ---------- peaks : 1D array_like Cardiac peak times, in seconds. + sample_rate : float + Sample rate of physio, in Hertz. slice_timings : 1D array_like Slice times, in seconds. n_scans : int @@ -81,31 +83,37 @@ def cardiac_phase(peaks, slice_timings, n_scans, t_r): Returns ------- phase_card : array_like - Cardiac phase signal. + Cardiac phase signal, of shape (n_scans,) """ - n_slices = np.shape(slice_timings) - phase_card = np.zeros(n_scans) + assert slice_timings.ndim == 1, "Slice times must be a 1D array" + n_slices = np.size(slice_timings) + phase_card = np.zeros((n_scans, n_slices)) + card_peaks_sec = peaks / sample_rate for i_slice in range(n_slices): - # find previous cardiac peaks - previous_card_peaks = np.asarray(np.nonzero(peaks < slice_timings[i_slice])) - if np.size(previous_card_peaks) == 0: - t1 = 0 - else: - last_peak = previous_card_peaks[0][-1] - t1 = peaks[last_peak] - - # find posterior cardiac peaks - next_card_peaks = np.asarray(np.nonzero(peaks > slice_timings[i_slice])) - if np.size(next_card_peaks) == 0: - t2 = n_scans * t_r - else: - next_peak = next_card_peaks[0][0] - t2 = peaks[next_peak] - - # compute cardiac phase - phase_card[i_slice] = (2 * np.math.pi * (slice_timings[i_slice] - t1)) / ( - t2 - t1 - ) + # generate slice acquisition timings across all scans + times_crSlice = t_r * np.arange(n_scans) + slice_timings[i_slice] + phase_card_crSlice = np.zeros(n_scans) + for j_scan in range(n_scans): + previous_card_peaks = np.asarray( + np.nonzero(card_peaks_sec < times_crSlice[j_scan]) + ) + if np.size(previous_card_peaks) == 0: + t1 = 0 + else: + last_peak = previous_card_peaks[0][-1] + t1 = card_peaks_sec[last_peak] + next_card_peaks = np.asarray( + np.nonzero(card_peaks_sec > times_crSlice[j_scan]) + ) + if np.size(next_card_peaks) == 0: + t2 = n_scans * t_r + else: + next_peak = next_card_peaks[0][0] + t2 = card_peaks_sec[next_peak] + phase_card_crSlice[j_scan] = ( + 2 * np.math.pi * (times_crSlice[j_scan] - t1) + ) / (t2 - t1) + phase_card[:, i_slice] = phase_card_crSlice return phase_card diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py index 7c5ae1c..def8f06 100644 --- a/phys2denoise/metrics/chest_belt.py +++ b/phys2denoise/metrics/chest_belt.py @@ -240,14 +240,15 @@ def respiratory_phase(resp, sample_rate, n_scans, slice_timings, t_r): Returns ------- phase_resp : array_like - Respiratory phase signal. + Respiratory phase signal, of shape (n_scans, n_slices). """ - n_slices = np.shape(slice_timings) + assert slice_timings.ndim == 1, "Slice times must be a 1D array" + n_slices = np.size(slice_timings) phase_resp = np.zeros((n_scans, n_slices)) # generate histogram from respiratory signal # TODO: Replace with numpy.histogram - resp_hist, resp_hist_bins = plt.hist(resp, bins=100) + resp_hist, resp_hist_bins, _ = plt.hist(resp, bins=100) # first compute derivative of respiration signal resp_diff = np.diff(resp, n=1) diff --git a/phys2denoise/tests/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py new file mode 100644 index 0000000..4b10ebd --- /dev/null +++ b/phys2denoise/tests/test_metrics_cardiac.py @@ -0,0 +1,41 @@ +"""Tests for phys2denoise.metrics.cardiac.""" +import numpy as np + +from phys2denoise.metrics import cardiac + + +def test_crf_smoke(): + """Basic smoke test for CRF calculation.""" + samplerate = 0.01 # in seconds + oversampling = 20 + time_length = 20 + onset = 0 + tr = 0.72 + crf_arr = cardiac.crf( + samplerate, + oversampling=oversampling, + time_length=time_length, + onset=onset, + tr=tr, + ) + pred_len = np.rint(time_length / (tr / oversampling)) + assert crf_arr.ndim == 1 + assert crf_arr.shape == pred_len + + +def test_cardiac_phase(): + """Basic smoke test for cardiac phase calculation.""" + t_r = 1.0 + n_scans = 200 + sample_rate = 1 / 0.01 + slice_timings = np.linspace(0, t_r, 22)[1:-1] + peaks = np.array([0.534, 0.577, 10.45, 20.66, 50.55, 90.22]) + card_phase = cardiac.cardiac_phase( + peaks, + sample_rate=sample_rate, + slice_timings=slice_timings, + n_scans=n_scans, + t_r=t_r, + ) + assert card_phase.ndim == 2 + assert card_phase.shape == (n_scans, slice_timings.size) diff --git a/phys2denoise/tests/test_metrics_chest_belt.py b/phys2denoise/tests/test_metrics_chest_belt.py new file mode 100644 index 0000000..58a3021 --- /dev/null +++ b/phys2denoise/tests/test_metrics_chest_belt.py @@ -0,0 +1,42 @@ +"""Tests for phys2denoise.metrics.chest_belt.""" +import numpy as np + +from phys2denoise.metrics import chest_belt + + +def test_rrf_smoke(): + """Basic smoke test for RRF calculation.""" + samplerate = 0.01 # in seconds + oversampling = 20 + time_length = 20 + onset = 0 + tr = 0.72 + rrf_arr = chest_belt.rrf( + samplerate, + oversampling=oversampling, + time_length=time_length, + onset=onset, + tr=tr, + ) + pred_len = int(np.rint(time_length / (tr / oversampling))) + assert rrf_arr.ndim == 1 + assert rrf_arr.size == pred_len + + +def test_respiratory_phase(): + """Basic smoke test for respiratory phase calculation.""" + t_r = 1.0 + n_scans = 200 + sample_rate = 1 / 0.01 + slice_timings = np.linspace(0, t_r, 22)[1:-1] + n_samples = int(np.rint((n_scans * t_r) * sample_rate)) + resp = np.random.normal(size=n_samples) + resp_phase = chest_belt.respiratory_phase( + resp, + sample_rate=sample_rate, + slice_timings=slice_timings, + n_scans=n_scans, + t_r=t_r, + ) + assert resp_phase.ndim == 2 + assert resp_phase.shape == (n_scans, slice_timings.size) From 7e8e4c5c52e8edd7a846de60076f7adbd4737625 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 25 Feb 2021 14:11:17 -0500 Subject: [PATCH 51/52] More improvements. --- phys2denoise/metrics/chest_belt.py | 84 ++++++++----------- phys2denoise/tests/test_metrics_cardiac.py | 2 +- phys2denoise/tests/test_metrics_chest_belt.py | 33 +++++++- 3 files changed, 68 insertions(+), 51 deletions(-) diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py index def8f06..2d9459a 100644 --- a/phys2denoise/metrics/chest_belt.py +++ b/phys2denoise/metrics/chest_belt.py @@ -3,7 +3,7 @@ import numpy as np import pandas as pd from scipy.ndimage.filters import convolve1d -from scipy.signal import detrend, resample +from scipy.signal import detrend from scipy.stats import zscore mpl.use("TkAgg") @@ -15,17 +15,18 @@ @due.dcite(references.POWER_2018) -def rpv(belt_ts, window): +def rpv(resp, window=6): """Calculate respiratory pattern variability. Parameters ---------- - belt_ts - window + resp : 1D array_like + window : int Returns ------- - rpv_arr + rpv_val : float + Respiratory pattern variability value. Notes ----- @@ -43,35 +44,29 @@ def rpv(belt_ts, window): 115, pp. 2105-2114, 2018. """ # First, z-score respiratory traces - resp_z = zscore(belt_ts) + resp_z = zscore(resp) # Collect upper envelope rpv_upper_env = utils.rms_envelope_1d(resp_z, window) # Calculate standard deviation - rpv_arr = np.std(rpv_upper_env) - return rpv_arr + rpv_val = np.std(rpv_upper_env) + return rpv_val @due.dcite(references.POWER_2020) -def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)): +def env(resp, samplerate, window=10): """Calculate respiratory pattern variability across a sliding window. Parameters ---------- - belt_ts : (X,) :obj:`numpy.ndarray` + resp : (X,) :obj:`numpy.ndarray` A 1D array with the respiratory belt time series. samplerate : :obj:`float` - Sampling rate for belt_ts, in Hertz. - out_samplerate : :obj:`float` - Sampling rate for the output time series in seconds. - Corresponds to TR in fMRI data. + Sampling rate for resp, in Hertz. window : :obj:`int`, optional - Size of the sliding window, in the same units as out_samplerate. - Default is 6. - lags : (Y,) :obj:`tuple` of :obj:`int`, optional - List of lags to apply to the rv estimate. In the same units as - out_samplerate. + Size of the sliding window, in seconds. + Default is 10. Returns ------- @@ -92,42 +87,39 @@ def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)): young adults scanned at rest, including systematic changes and 'missed' deep breaths," Neuroimage, vol. 204, 2020. """ - window = window * samplerate / out_samplerate + # Convert window to Hertz + window = int(window * samplerate) + # Calculate RPV across a rolling window env_arr = ( - pd.Series(belt_ts).rolling(window=window, center=True).apply(rpv, window=window) + pd.Series(resp) + .rolling(window=window, center=True) + .apply(rpv, args=(window,)) ) env_arr[np.isnan(env_arr)] = 0.0 return env_arr @due.dcite(references.CHANG_GLOVER_2009) -def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)): +def rv(resp, samplerate, window=6, lags=(0,)): """Calculate respiratory variance. Parameters ---------- - belt_ts : (X,) :obj:`numpy.ndarray` + resp : (X,) :obj:`numpy.ndarray` A 1D array with the respiratory belt time series. samplerate : :obj:`float` - Sampling rate for belt_ts, in Hertz. - out_samplerate : :obj:`float` - Sampling rate for the output time series in seconds. - Corresponds to TR in fMRI data. + Sampling rate for resp, in Hertz. window : :obj:`int`, optional - Size of the sliding window, in the same units as out_samplerate. + Size of the sliding window, in seconds. Default is 6. - lags : (Y,) :obj:`tuple` of :obj:`int`, optional - List of lags to apply to the rv estimate. In the same units as - out_samplerate. Returns ------- - rv_out : (Z, 2Y) :obj:`numpy.ndarray` - Respiratory variance values, with lags applied, downsampled to - out_samplerate, convolved with an RRF, and detrended/normalized. - The first Y columns are not convolved with the RRF, while the second Y - columns are. + rv_out : (X, 2) :obj:`numpy.ndarray` + Respiratory variance values. + The first column is raw RV values, after detrending/normalization. + The second column is RV values convolved with the RRF, after detrending/normalization. Notes ----- @@ -145,25 +137,19 @@ def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)): end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage, issue 4, vol. 47, pp. 1381-1393, 2009. """ + # Convert window to Hertz + window = int(window * samplerate) + # Raw respiratory variance - rv_arr = pd.Series(belt_ts).rolling(window=window, center=True).std() + rv_arr = pd.Series(resp).rolling(window=window, center=True).std() rv_arr[np.isnan(rv_arr)] = 0.0 - # Apply lags - n_out_samples = int((belt_ts.shape[0] / samplerate) / out_samplerate) - # convert lags from out_samplerate to samplerate - delays = [abs(int(lag * samplerate)) for lag in lags] - rv_with_lags = utils.apply_lags(rv_arr, lags=delays) - - # Downsample to out_samplerate - rv_with_lags = resample(rv_with_lags, num=n_out_samples, axis=0) - # Convolve with rrf - rrf_arr = rrf(out_samplerate, oversampling=1) - rv_convolved = convolve1d(rv_with_lags, rrf_arr, axis=0) + rrf_arr = rrf(samplerate, oversampling=1) + rv_convolved = convolve1d(rv_arr, rrf_arr, axis=0) # Concatenate the raw and convolved versions - rv_combined = np.hstack((rv_with_lags, rv_convolved)) + rv_combined = np.stack((rv_arr, rv_convolved), axis=-1) # Detrend and normalize rv_combined = rv_combined - np.mean(rv_combined, axis=0) diff --git a/phys2denoise/tests/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py index 4b10ebd..bb10d76 100644 --- a/phys2denoise/tests/test_metrics_cardiac.py +++ b/phys2denoise/tests/test_metrics_cardiac.py @@ -23,7 +23,7 @@ def test_crf_smoke(): assert crf_arr.shape == pred_len -def test_cardiac_phase(): +def test_cardiac_phase_smoke(): """Basic smoke test for cardiac phase calculation.""" t_r = 1.0 n_scans = 200 diff --git a/phys2denoise/tests/test_metrics_chest_belt.py b/phys2denoise/tests/test_metrics_chest_belt.py index 58a3021..f256918 100644 --- a/phys2denoise/tests/test_metrics_chest_belt.py +++ b/phys2denoise/tests/test_metrics_chest_belt.py @@ -23,7 +23,7 @@ def test_rrf_smoke(): assert rrf_arr.size == pred_len -def test_respiratory_phase(): +def test_respiratory_phase_smoke(): """Basic smoke test for respiratory phase calculation.""" t_r = 1.0 n_scans = 200 @@ -40,3 +40,34 @@ def test_respiratory_phase(): ) assert resp_phase.ndim == 2 assert resp_phase.shape == (n_scans, slice_timings.size) + + +def test_rpv_smoke(): + """Basic smoke test for respiratory pattern variability calculation.""" + n_samples = 2000 + resp = np.random.normal(size=n_samples) + window = 50 + rpv_val = chest_belt.rpv(resp, window) + assert isinstance(rpv_val, float) + + +def test_env_smoke(): + """Basic smoke test for ENV calculation.""" + n_samples = 2000 + resp = np.random.normal(size=n_samples) + samplerate = 1 / 0.01 + window = 6 + env_arr = chest_belt.env(resp, samplerate=samplerate, window=window) + assert env_arr.ndim == 1 + assert env_arr.shape == (n_samples,) + + +def test_rv_smoke(): + """Basic smoke test for respiratory variance calculation.""" + n_samples = 2000 + resp = np.random.normal(size=n_samples) + samplerate = 1 / 0.01 + window = 6 + rv_arr = chest_belt.rv(resp, samplerate=samplerate, window=window) + assert rv_arr.ndim == 2 + assert rv_arr.shape == (n_samples, 2) From df0a85bb851cea87d817274e788e755fd584e50e Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 25 Feb 2021 14:12:27 -0500 Subject: [PATCH 52/52] Cleanup. --- phys2denoise/metrics/chest_belt.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py index 2d9459a..e8dafea 100644 --- a/phys2denoise/metrics/chest_belt.py +++ b/phys2denoise/metrics/chest_belt.py @@ -92,9 +92,7 @@ def env(resp, samplerate, window=10): # Calculate RPV across a rolling window env_arr = ( - pd.Series(resp) - .rolling(window=window, center=True) - .apply(rpv, args=(window,)) + pd.Series(resp).rolling(window=window, center=True).apply(rpv, args=(window,)) ) env_arr[np.isnan(env_arr)] = 0.0 return env_arr