Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Include upper boundary of regression ranges (and add a legacy option for pythonic ranges). #39

Merged
merged 1 commit into from
Feb 23, 2022
Merged

Include upper boundary of regression ranges (and add a legacy option for pythonic ranges). #39

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 20 additions & 11 deletions phys2cvr/cli/run.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -280,16 +280,23 @@ 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',
type=float,
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 '
Expand Down Expand Up @@ -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',
Expand All @@ -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)

Expand Down Expand Up @@ -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
Expand All @@ -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. '
Expand Down
44 changes: 28 additions & 16 deletions phys2cvr/phys2cvr.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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')
Expand All @@ -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:
Expand Down Expand Up @@ -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])
Expand All @@ -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')
Expand All @@ -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:
Expand All @@ -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}')
Expand All @@ -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')
Expand Down
27 changes: 17 additions & 10 deletions phys2cvr/stats.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down Expand Up @@ -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
-------
Expand Down Expand Up @@ -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,
Expand Down