Skip to content

Code reorganisation to use lag maps #2

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

Merged
beccaclements99 merged 2 commits into from
Mar 12, 2026
Merged

Code reorganisation to use lag maps #2

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6792125
Add info on loading data
smoia Mar 12, 2026
dbf88a4
Move lag map reading and lag_step and lag_max estimation out of regr_…
smoia Mar 12, 2026
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
102 changes: 54 additions & 48 deletions phys2cvr/workflows.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -311,6 +311,7 @@ def phys2cvr(
f'R^2 model {r2model} not supported. Supported models are {stats.R2MODEL}'
)

LGR.info('Load functional data')
if func_is_1d:
if tr:
func_avg = io.load_array(fname_func)
Expand Down Expand Up @@ -345,6 +346,7 @@ def phys2cvr(

# Read mask (and mask func) if provided
if fname_mask:
LGR.info('Load mask to restrict operations on functional data')
_, mask, _ = io.load_nifti_get_mask(fname_mask, is_mask=True)
if func.shape[:3] != mask.shape:
raise ValueError(f'{fname_mask} and {fname_func} have different sizes!')
Expand All @@ -364,6 +366,7 @@ def phys2cvr(

# Read roi if provided
if fname_roi:
LGR.info('Load ROI to obtain a reference from functional data')
_, roi, _ = io.load_nifti_get_mask(fname_roi, is_mask=True)
if func.shape[:3] != roi.shape:
raise ValueError(f'{fname_roi} and {fname_func} have different sizes!')
Expand All @@ -382,6 +385,7 @@ def phys2cvr(
LGR.info(f'Obtaining filtered average signal in {roiref}')
func_avg = signal.filter_signal(func_avg, tr, lowcut, highcut, butter_order)

LGR.info('Load physiological data')
if fname_co2 is None:
LGR.info(f'Computing "CVR" (approximation) maps using {fname_func} only')
if func_is_1d:
Expand Down Expand Up @@ -462,6 +466,39 @@ def phys2cvr(
else:
petco2hrf = co2

# If the user provided a lag map, read it and extract information from it
if lag_map:
LGR.info('Load lag map')
lag, _, _ = io.load_nifti_get_mask(lag_map)
if func.shape[:3] != lag.shape:
raise ValueError(f'{lag_map} and {fname_func} have different sizes!')

# Read lag_step and lag_max from file (or try to)
lag = lag * mask

lag_list = np.unique(lag[mask])

if lag_step is None:
# np.unique sorts results already
lag_step = np.unique(np.round(lag_list[1:] - lag_list[:-1], 3))
if lag_step.size > 1:
# Check if extra steps found are multiple of the first within numerical tolerance
if not np.isclose(np.mod(lag_step, lag_step[0]), 0).all():
raise ValueError(
f'phys2cvr found different delta lags in {lag_map}: {lag_step}'
)

lag_step = lag_step[0]
LGR.warning(f'phys2cvr detected a delta lag of {lag_step} seconds')
else:
LGR.warning(f'Forcing delta lag to be {lag_step}')

if lag_max is None:
lag_max = np.abs(lag_list).max()
LGR.warning(f'phys2cvr detected a max lag of {lag_max} seconds')
else:
LGR.warning(f'Forcing max lag to be {lag_max}')

# If a regressor dir is specified, try load the data,
# If failing or otherwise, compute the regressors.
if run_regression and regr_dir is not None:
Expand All @@ -476,35 +513,6 @@ def phys2cvr(
regr_dir = None

if regr_dir is None:
if lag_map:
lag, _, _ = io.load_nifti_get_mask(lag_map)
if func.shape[:3] != lag.shape:
raise ValueError(f'{lag_map} and {fname_func} have different sizes!')

# Read lag_step and lag_max from file (or try to)
lag = lag * mask

lag_list = np.unique(lag[mask > 0])

if lag_step is None:
lag_step = np.unique(np.round(lag_list[1:] - lag_list[:-1], 3))
if lag_step.size > 1:
raise ValueError(
f'phys2cvr found different delta lags in {lag_map}'
)
else:
lag_step = lag_step[0]
LGR.warning(f'phys2cvr detected a delta lag of {lag_step} seconds')
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')
else:
LGR.warning(f'Forcing max lag to be {lag_max}')
regr, regr_shifts = create_physio_regressor(
func_avg,
petco2hrf,
Expand Down Expand Up @@ -608,37 +616,25 @@ def phys2cvr(
r_square, oimg, f'{fname_out_func}_r_square_simple{fname_ext}'
)

if (
lagged_regression
and regr_shifts is not None
and ((lag_max and lag_step) or lag_map)
):
if lag_max:
LGR.info(
f'Running lagged CVR estimation with max lag = {lag_max}! '
'(might take a while...)'
)
elif lag_map is not None:
if lagged_regression and regr_shifts is not None and (lag_max and lag_step):
# If user specified a lag map, run regression based on it (see "Load lag map")
if lag_map is not None:
LGR.info(
f'Running lagged CVR estimation with lag map {lag_map}! '
'(might take a while...)'
)
if legacy:
nrep = int(lag_max * freq * 2)
else:
nrep = int(lag_max * freq * 2) + 1

# If user specified a lag map, use that one to regress things
if lag_map:
step = int(lag_step * freq)

lag_idx = np.round((lag + lag_max) * freq / step).astype(int)
lag_idx_list = np.unique(lag_idx)

# Prepare empty matrices
beta = np.empty_like(lag, dtype='float32')
tstat = np.empty_like(lag, dtype='float32')

for index, i in enumerate(lag_idx_list):
LGR.info(f'Perform L-GLM number {index + 1} of {len(lag_idx_list)}')
for n, i in enumerate(lag_idx_list):
LGR.info(f'Perform L-GLM number {n + 1} of {len(lag_idx_list)}')
regr = regr_shifts[(i * step), :, np.newaxis]

x1D = os.path.join(outdir, 'mat', f'mat_{i:04g}.1D')
Expand All @@ -656,6 +652,16 @@ def phys2cvr(
)

else:
LGR.info(
f'Running lagged CVR estimation with max lag = {lag_max}! '
'(might take a while...)'
)

if legacy:
nrep = int(lag_max * freq * 2)
else:
nrep = int(lag_max * freq * 2) + 1

# Check the number of repetitions first
if lag_step:
step = int(lag_step * freq)
Expand Down