From a0a97c10e124280655efff7a3ecf98bcf97c0477 Mon Sep 17 00:00:00 2001 From: smoia Date: Wed, 23 Feb 2022 21:49:20 +0100 Subject: [PATCH] Add legacy option for pythonic vs non-pythonic regression ranges. --- phys2cvr/cli/run.py | 31 ++++++++++++++++++++----------- phys2cvr/phys2cvr.py | 44 ++++++++++++++++++++++++++++---------------- phys2cvr/stats.py | 27 +++++++++++++++++---------- 3 files changed, 65 insertions(+), 37 deletions(-) diff --git a/phys2cvr/cli/run.py b/phys2cvr/cli/run.py index 0937a3a..466ce59 100644 --- a/phys2cvr/cli/run.py +++ b/phys2cvr/cli/run.py @@ -280,9 +280,9 @@ def _get_parser(): help=('Maximum lag to consider during lag regression ' 'in seconds. The same lag will be considered in ' 'both directions.\n' - 'Remember that being this python, the upper limit ' - 'is excluded from the computation, i.e. 9 is ' - '[-9, +8.7] or [-9, +9).'), + 'Despite the code being python, the upper limit ' + 'is included in the computation. E.g. -lm 9 ' + '-ls .3 means [-9, +9] (61 regressors).'), default=None) opt_lreg.add_argument('-ls', '--lag-step', dest='lag_step', @@ -290,6 +290,13 @@ def _get_parser(): help=('Lag step to consider during lagged regression ' 'in seconds. Default is 0.3 seconds.'), default=None) + opt_lreg.add_argument('--legacy', + dest='legacy', + action='store_true', + help=('Use pythonic ranges, i.e. the upper limit ' + 'is excluded from the computation.\nE.g. -lm 9 ' + '-ls .3 means [-9, +8.7] or [-9, +9) (60 regressors).'), + default=False) opt_regr = parser.add_argument_group('Optional Arguments to re-run a lagged ' 'regression (also useful to use a lag estimation ' @@ -321,8 +328,9 @@ def _get_parser(): 'as used in:\n' 'S. Moia, et al., \'ICA-based denoising strategies in ' 'breath-hold induced cerebrovascular reactivity mapping ' - 'with multi echo BOLD fMRI\' (2021), NeuroImage\n' - 'Same as setting --lag-max 9 --lag-step 0.3 --r2full'), + 'with multi echo BOLD fMRI\' (2021), NeuroImage.\n' + 'Same as setting --lag-max 9 --lag-step 0.3 ' + '--legacy --r2full'), default=None) opt_conf.add_argument('--brightspin-clinical', dest='workflow_config', @@ -338,16 +346,16 @@ def _get_parser(): help=('Estimate CVR using the average timeseries in the ' '0.02-0.04 frequency spectrum, as used in:\n' 'P. Liu, et al., \'Cerebrovascular reactivity ' - 'mapping without gas challenges\' (2017), NeuroImage\n' - 'Same as setting --apply-filter -hf -lf ' + 'mapping without gas challenges\' (2017), NeuroImage.\n' + 'Same as setting --apply-filter -hf 0.04 -lf 0.02 ' '-skip_conv -skip_lagreg -co2 \'\' '), default=None) opt_conf.add_argument('--baltimore-lag', dest='workflow_config', action='store_const', const='baltimore-lag', - help=('Like "baltimore", but use a L-GLM instead\n' - 'Same as setting --apply-filter -hf -lf ' + help=('Like "baltimore", but use a L-GLM instead.\n' + 'Same as setting --apply-filter -hf 0.04 -lf 0.02 ' '-skip_conv -co2 \'\' '), default=None) @@ -399,6 +407,7 @@ def _check_opt_conf(parser): parser.lagged_regression = True parser.run_conv = True parser.apply_filter = False + parser.legacy = True parser.r2model = 'full' elif parser.workflow_config == 'brightspin-clinical': parser.lag_max = 20 @@ -412,14 +421,14 @@ def _check_opt_conf(parser): parser.apply_filter = True parser.lowcut = 0.02 parser.highcut = 0.04 - parser.fname_co2 = '' + parser.fname_co2 = None parser.lagged_regression = False elif parser.workflow_config == 'baltimore-lag': parser.run_conv = False parser.apply_filter = True parser.lowcut = 0.02 parser.highcut = 0.04 - parser.fname_co2 = '' + parser.fname_co2 = None parser.lagged_regression = True else: raise NotImplementedError(f'{parser.workflow_config} is not configured. ' diff --git a/phys2cvr/phys2cvr.py b/phys2cvr/phys2cvr.py index c3762a5..a4f1407 100755 --- a/phys2cvr/phys2cvr.py +++ b/phys2cvr/phys2cvr.py @@ -51,7 +51,7 @@ def phys2cvr(fname_func, fname_co2=None, fname_pidx=None, fname_roi=None, fname_ outdir=None, freq=None, tr=None, trial_len=None, n_trials=None, highcut=None, lowcut=None, apply_filter=False, run_regression=False, lagged_regression=True, r2model='full', lag_max=None, - lag_step=None, l_degree=0, denoise_matrix=[], scale_factor=None, + lag_step=None, legacy=False, l_degree=0, denoise_matrix=[], scale_factor=None, lag_map=None, regr_dir=None, run_conv=True, quiet=False, debug=False): """ Run main workflow of phys2cvr. @@ -129,11 +129,16 @@ def phys2cvr(fname_func, fname_co2=None, fname_pidx=None, fname_roi=None, fname_ Default: 'full' lag_max : int or float, optional Limits (both positive and negative) of the temporal area to explore, - expressed in seconds (i.e. ±9 seconds). + expressed in seconds (e.g. ±9 seconds). Caution: this is not a pythonic + range, but a real range, i.e. the upper limit is included (e.g. [-9, +9]). Default: None lag_step : int or float, optional Step of the lag to take into account in seconds. Default: None + legacy : bool, optional + If True, use pythonic ranges when creating the regressors, i.e. exclude + the upper range (e.g. [-9, +9) ). + Default: False l_degree : int, optional Only used if performing the regression step. Highest order of the Legendre polynomial to add to the denoising matrix. @@ -363,7 +368,7 @@ def phys2cvr(fname_func, fname_co2=None, fname_pidx=None, fname_roi=None, fname_ if regr_dir is None: regr, regr_shifts = stats.get_regr(func_avg, petco2hrf, tr, freq, outname, lag_max, trial_len, n_trials, - '.1D', lagged_regression) + '.1D', lagged_regression, legacy) elif run_regression: try: regr = np.genfromtxt(f'{outname}_petco2hrf.1D') @@ -372,7 +377,8 @@ def phys2cvr(fname_func, fname_co2=None, fname_pidx=None, fname_roi=None, fname_ 'Estimating it.') regr, regr_shifts = stats.get_regr(func_avg, petco2hrf, tr, freq, outname, lag_max, trial_len, - n_trials, '.1D') + n_trials, '.1D', lagged_regression, + legacy) # Run internal regression if required and possible! if func_is_nifti and run_regression: @@ -430,12 +436,10 @@ def phys2cvr(fname_func, fname_co2=None, fname_pidx=None, fname_roi=None, fname_ elif lag_map is not None: LGR.info(f'Running lagged CVR estimation with lag map {lag_map}! ' '(might take a while...)') - - nrep = int(lag_max * freq * 2) - if lag_step: - step = int(lag_step * freq) + if legacy: + nrep = int(lag_max * freq * 2) else: - step = 1 + nrep = int(lag_max * freq * 2) + 1 if regr_dir: outprefix = os.path.join(regr_dir, os.path.split(outname)[1]) @@ -460,6 +464,8 @@ def phys2cvr(fname_func, fname_co2=None, fname_pidx=None, fname_roi=None, fname_ else: LGR.warning(f'Forcing delta lag to be {lag_step}') + step = int(lag_step * freq) + if lag_max is None: lag_max = np.abs(lag_list).max() LGR.warning(f'phys2cvr detected a max lag of {lag_max} seconds') @@ -476,7 +482,7 @@ def phys2cvr(fname_func, fname_co2=None, fname_pidx=None, fname_roi=None, fname_ for i in lag_idx_list: LGR.info(f'Perform L-GLM for lag {lag_list[i]} ({i+1} of ' - f'{nrep // step})') + f'{len(lag_idx_list)}') try: regr = regr_shifts[:, (i*step)] except NameError: @@ -490,13 +496,19 @@ def phys2cvr(fname_func, fname_co2=None, fname_pidx=None, fname_roi=None, fname_ x1D) else: + # Check the number of repetitions first + if lag_step: + step = int(lag_step * freq) + else: + step = 1 + lag_range = list(range(0, nrep, step)) # Prepare empty matrices - r_square_all = np.empty(list(func.shape[:3]) + [nrep // step]) - beta_all = np.empty(list(func.shape[:3]) + [nrep // step]) - tstat_all = np.empty(list(func.shape[:3]) + [nrep // step]) + r_square_all = np.empty(list(func.shape[:3]) + [len(lag_range)]) + beta_all = np.empty(list(func.shape[:3]) + [len(lag_range)]) + tstat_all = np.empty(list(func.shape[:3]) + [len(lag_range)]) - for n, i in enumerate(range(0, nrep, step)): - LGR.info(f'Perform L-GLM number {n+1} of {nrep // step}') + for n, i in enumerate(lag_range): + LGR.info(f'Perform L-GLM number {n+1} of {len(lag_range)}') try: regr = regr_shifts[:, i] LGR.debug(f'Using shift {i} from matrix in memory: {regr}') @@ -517,7 +529,7 @@ def phys2cvr(fname_func, fname_co2=None, fname_pidx=None, fname_roi=None, fname_ if debug: LGR.debug('Export all betas, tstats, and R^2 volumes.') newdim_all = deepcopy(img.header['dim']) - newdim_all[0], newdim_all[4] = 4, int(nrep // step) + newdim_all[0], newdim_all[4] = 4, int(len(lag_range)) oimg_all = deepcopy(img) oimg_all.header['dim'] = newdim_all io.export_nifti(r_square_all, oimg_all, f'{fname_out_func}_r_square_all') diff --git a/phys2cvr/stats.py b/phys2cvr/stats.py index d13aae5..ecfa75e 100644 --- a/phys2cvr/stats.py +++ b/phys2cvr/stats.py @@ -74,7 +74,8 @@ def x_corr(func, co2, lastrep, firstrep=0, offset=0): def get_regr(func_avg, petco2hrf, tr, freq, outname, lag_max=None, - trial_len=None, n_trials=None, ext='.1D', lagged_regression=True): + trial_len=None, n_trials=None, ext='.1D', lagged_regression=True, + legacy=False): """ Create regressor(s) of interest for nifti GLM. @@ -106,7 +107,10 @@ def get_regr(func_avg, petco2hrf, tr, freq, outname, lag_max=None, Extension to be used for the exported regressors. lagged_regression : bool, optional Estimate regressors for each possible lag of `petco2hrf`. - If True, the maximum number of regressors will be `(freq*lag_max*2)-1` + If True, the maximum number of regressors will be `(freq*lag_max*2)+1` + legacy : bool, optional + If True, exclude the upper lag limit from the regression estimation. + If True, the maximum number of regressors will be `(freq*lag_max*2)` Returns ------- @@ -203,24 +207,27 @@ def get_regr(func_avg, petco2hrf, tr, freq, outname, lag_max=None, os.makedirs(os.path.join(os.path.split(outname)[0], 'regr'), exist_ok=True) # Set num of fine shifts: 9 seconds is a bit more than physiologically feasible - nrep = int(lag_max * freq) - - petco2hrf_shifts = np.empty([func_len, nrep*2]) + negrep = int(lag_max * freq) + if legacy: + posrep = negrep + else: + posrep = negrep + 1 + petco2hrf_shifts = np.empty([func_len, negrep+posrep]) # Padding regressor for shift, and padding optshift too - if (optshift - nrep) < 0: - lpad = nrep - optshift + if (optshift - negrep) < 0: + lpad = negrep - optshift else: lpad = 0 - if (optshift + nrep + len_upd) > len(petco2hrf): - rpad = (optshift + nrep + len_upd) - len(petco2hrf) + if (optshift + posrep + len_upd) > len(petco2hrf): + rpad = (optshift + posrep + len_upd) - len(petco2hrf) else: rpad = 0 petco2hrf_padded = np.pad(petco2hrf, (int(lpad), int(rpad)), 'mean') - for n, i in enumerate(range(-nrep, nrep)): + for n, i in enumerate(range(-negrep, posrep)): petco2hrf_lagged = petco2hrf_padded[optshift+lpad-i:optshift+lpad-i+len_upd] petco2hrf_shifts[:, n] = io.export_regressor(regr_t, petco2hrf_lagged, func_t, outprefix,