diff --git a/.all-contributorsrc b/.all-contributorsrc
index a3e0a5d..431b2ba 100644
--- a/.all-contributorsrc
+++ b/.all-contributorsrc
@@ -29,6 +29,17 @@
"infra"
]
},
+ {
+ "login": "beccaclements99",
+ "name": "Becca Clements",
+ "avatar_url": "https://avatars.githubusercontent.com/u/95257798?v=4",
+ "profile": "https://github.com/beccaclements99",
+ "contributions": [
+ "bug",
+ "code",
+ "userTesting"
+ ]
+ },
{
"login": "avigotsky",
"name": "Andrew Vigotsky",
@@ -38,6 +49,16 @@
"code"
]
},
+ {
+ "login": "ccomellalue",
+ "name": "ccomellalue",
+ "avatar_url": "https://avatars.githubusercontent.com/u/209748720?v=4",
+ "profile": "https://github.com/ccomellalue",
+ "contributions": [
+ "code",
+ "ideas"
+ ]
+ },
{
"login": "merelvdthiel",
"name": "merelvdthiel",
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1f7d1c5..ccdce20 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,38 @@
+# 0.28.3 (Thu Mar 12 2026)
+
+:tada: This release contains work from a new contributor! :tada:
+
+Thank you, Becca Clements ([@beccaclements99](https://github.com/beccaclements99)), for all your work!
+
+#### 🐛 Bug Fix
+
+- Fix usage of previously computed lag maps to run a (new) lagged regression [#154](https://github.com/smoia/phys2cvr/pull/154) ([@beccaclements99](https://github.com/beccaclements99) [@pre-commit-ci[bot]](https://github.com/pre-commit-ci[bot]) [@smoia](https://github.com/smoia))
+- Give an error message if peaks are not provided AND comp_endtidal is True [#156](https://github.com/smoia/phys2cvr/pull/156) ([@beccaclements99](https://github.com/beccaclements99))
+
+#### ⚠️ Pushed to `master`
+
+- tests: Reorganise tests, remove or mark xfails the failing ones ([@smoia](https://github.com/smoia))
+- docs: update licenses ([@smoia](https://github.com/smoia))
+- docs: Fix logos ([@smoia](https://github.com/smoia))
+
+#### ⚠️ Tests
+
+- Fix workflow halting on missing peaks when they are not required (i.e. when skipping end-tidal convolution) [#153](https://github.com/smoia/phys2cvr/pull/153) ([@smoia](https://github.com/smoia) [@beccaclements99](https://github.com/beccaclements99) [@pre-commit-ci[bot]](https://github.com/pre-commit-ci[bot]))
+
+#### 🏠 Internal
+
+- [pre-commit.ci] pre-commit autoupdate [#152](https://github.com/smoia/phys2cvr/pull/152) ([@pre-commit-ci[bot]](https://github.com/pre-commit-ci[bot]))
+- Bump actions/checkout from 5 to 6 [#151](https://github.com/smoia/phys2cvr/pull/151) ([@dependabot[bot]](https://github.com/dependabot[bot]))
+
+#### Authors: 4
+
+- [@dependabot[bot]](https://github.com/dependabot[bot])
+- [@pre-commit-ci[bot]](https://github.com/pre-commit-ci[bot])
+- Becca Clements ([@beccaclements99](https://github.com/beccaclements99))
+- Stefano Moia ([@smoia](https://github.com/smoia))
+
+---
+
# 0.28.2 (Wed Nov 19 2025)
#### 🐛 Bug Fix
diff --git a/README.md b/README.md
index c67265e..01caf92 100644
--- a/README.md
+++ b/README.md
@@ -20,7 +20,7 @@ phys2cvr
[](https://pypi.org/project/phys2cvr/)
-[](#contributors)
+[](#contributors)
A python-based tool to generate regressor for and/or estimate CVR maps and their lag.
@@ -92,7 +92,11 @@ Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/d
Stefano Moia
💻 🤔 🚇 📆 |
Kristina Zvolanek
💻 🐛 🚇 |
+ 
Becca Clements
🐛 💻 📓 |
Andrew Vigotsky
💻 |
+ 
ccomellalue
💻 🤔 |
+
+
merelvdthiel
📖 |
razkin
🎨 📖 ⚠️ |
@@ -108,7 +112,7 @@ Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/d
License
-------
-Copyright 2021-2025, Stefano Moia & phys2cvr contributors.
+Copyright 2021-2026, Stefano Moia & phys2cvr contributors.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
diff --git a/phys2cvr/cli/run.py b/phys2cvr/cli/run.py
index e07b9f6..2597b92 100644
--- a/phys2cvr/cli/run.py
+++ b/phys2cvr/cli/run.py
@@ -701,7 +701,7 @@ def _check_opt_conf(parser):
"""
-Copyright 2021-2025, Stefano Moia & phys2cvr contributors.
+Copyright 2021-2026, Stefano Moia & phys2cvr contributors.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
diff --git a/phys2cvr/io.py b/phys2cvr/io.py
index 60cf05f..4028c79 100644
--- a/phys2cvr/io.py
+++ b/phys2cvr/io.py
@@ -119,6 +119,7 @@ def load_txt(fname, shape=None):
'.txt': ' ',
'.1d': ' ',
'.par': ' ',
+ '': ' ',
}
mtx = np.genfromtxt(fname, delimiter=delimiter_map.get(ext))
@@ -235,10 +236,11 @@ def load_array(fname, shape=''):
"""
_, _, ext = utils.check_ext(EXT_ARRAY, fname, scan=True, remove=True)
- if ext.lower() in EXT_1D:
- return load_txt(fname, shape=shape)
- if ext.lower() in EXT_MAT:
- return load_mat(fname, shape=shape)
+ if ext:
+ if ext.lower() in EXT_1D:
+ return load_txt(fname, shape=shape)
+ if ext.lower() in EXT_MAT:
+ return load_mat(fname, shape=shape)
raise NotImplementedError(
f'{fname} file extension {ext} was not found or is not supported yet'
@@ -268,7 +270,7 @@ def load_physio(fname):
def export_regressor(
- regressors_matrix, ntp, outname, suffix='petco2hrf', ext='.1D', axis=-1
+ regressors_matrix, ntp, outprefix, suffix='petco2hrf', ext='.1D', axis=-1
):
"""
Export generated regressors for fMRI analysis.
@@ -279,7 +281,7 @@ def export_regressor(
The regressors that needs to be exported, in its original sample
ntp : int
The number of fMRI timepoints
- outname : str or path
+ outprefix : str or path
Prefix of the output file - can contain a path.
suffix : str, optional
The suffix of the output file.
@@ -301,7 +303,7 @@ def export_regressor(
regressors_demeaned = regressors_matrix - regressors_matrix.mean(
axis=axis, keepdims=True
)
- np.savetxt(f'{outname}_{suffix}{ext}', regressors_demeaned, fmt='%.6f')
+ np.savetxt(f'{outprefix}_{suffix}{ext}', regressors_demeaned, fmt='%.6f')
return regressors_demeaned
@@ -324,7 +326,7 @@ def export_nifti(data, img, fname):
"""
-Copyright 2021-2025, Stefano Moia & phys2cvr contributors.
+Copyright 2021-2026, Stefano Moia & phys2cvr contributors.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
diff --git a/phys2cvr/regressors.py b/phys2cvr/regressors.py
index 9b42de7..c213666 100644
--- a/phys2cvr/regressors.py
+++ b/phys2cvr/regressors.py
@@ -100,7 +100,7 @@ def compute_petco2hrf(
'Arrays with more than 2 dimensions are not supported.'
)
- if comp_endtidal:
+ if comp_endtidal and pidx is not None:
petco2 = endtidal_interpolation(co2, pidx, axis=-1)
# Plot PetCO2 vs CO2
@@ -111,6 +111,10 @@ def compute_petco2hrf(
# Demean and export
petco2 = petco2 - petco2.mean()
np.savetxt(f'{outprefix}_petco2.1D', petco2, fmt='%.18f')
+ elif comp_endtidal and pidx is None:
+ raise ValueError(
+ 'End-tidal interpolation was requested, but peaks were not provided.'
+ )
else:
LGR.info(
'Skipping End Tidal interpolation of PetCO2 trace (if you provided raw CO2 '
@@ -412,7 +416,7 @@ def create_physio_regressor(
"""
-Copyright 2021-2025, Stefano Moia & phys2cvr contributors.
+Copyright 2021-2026, Stefano Moia & phys2cvr contributors.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
diff --git a/phys2cvr/signal.py b/phys2cvr/signal.py
index b290c85..597d079 100644
--- a/phys2cvr/signal.py
+++ b/phys2cvr/signal.py
@@ -70,7 +70,7 @@ def create_hrf(freq=40):
u = np.arange(0, p[6] / dt + 1, 1) - p[5] / dt
a1 = p[0] / p[2]
- b1 = 1 / p[3]
+ b1 = 1 / p[2]
a2 = p[1] / p[3]
b2 = 1 / p[3]
@@ -325,7 +325,7 @@ def resample_signal_freqs(ts, freq1, freq2, axis=-1):
"""
-Copyright 2021-2025, Stefano Moia & phys2cvr contributors.
+Copyright 2021-2026, Stefano Moia & phys2cvr contributors.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
diff --git a/phys2cvr/stats.py b/phys2cvr/stats.py
index e1fa273..df3f6b6 100644
--- a/phys2cvr/stats.py
+++ b/phys2cvr/stats.py
@@ -365,7 +365,7 @@ def regression(
"""
-Copyright 2021-2025, Stefano Moia & phys2cvr contributors.
+Copyright 2021-2026, Stefano Moia & phys2cvr contributors.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
diff --git a/phys2cvr/tests/test_io.py b/phys2cvr/tests/test_io.py
index f875211..54e5ec3 100644
--- a/phys2cvr/tests/test_io.py
+++ b/phys2cvr/tests/test_io.py
@@ -2,9 +2,14 @@
"""Tests for io."""
import os
+import sys
import nibabel as nib
import numpy as np
+import peakdet as pk
+import pymatreader
+import pytest
+from scipy.io import savemat
from phys2cvr import io
@@ -27,8 +32,78 @@ def test_load_nifti_get_mask(nifti_data, nifti_mask):
assert isinstance(img, nib.nifti1.Nifti1Image)
+@pytest.mark.parametrize(
+ 'extension, delimiter',
+ [
+ ['.csv', ','],
+ ['.csv.gz', ','],
+ ['.tsv', '\t'],
+ ['.tsv.gz', '\t'],
+ ['.txt', ' '],
+ ['.1d', ' '],
+ ['.par', ' '],
+ ['', ' '],
+ ],
+)
+def test_load_txt(extension, delimiter, testdir):
+ """Test load_txt."""
+ a = np.random.randn(5)
+ n = os.join.path(testdir, f'zoe{extension}')
+ np.savetxt(n, a, delimiter=delimiter)
+ b = io.load_txt(n)
+
+ assert (a == b).all()
+ os.remove(n)
+
+
+def test_load_mat():
+ """Test load_mat."""
+ a = np.random.randn(5)
+ b = 'likealeaf'
+ n = 'inthewind'
+
+ savemat(n, {'data': a, 'other': b})
+
+ c = io.load_mat(n)
+
+ assert (a == c).all()
+ os.remove(n)
+
+
+def test_load_array():
+ """Test load_mat."""
+ a = np.random.randn(5)
+ b = 'likealeaf'
+ n = 'inthewind.mat'
+ m = 'zoe.txt'
+ np.savetxt(m, a)
+
+ savemat(n, {'data': a, 'other': b})
+
+ d = io.load_array(n)
+ e = io.load_array(m)
+
+ assert (a == d).all()
+ assert (a == e).all()
+ os.remove(n)
+ os.remove(m)
+
+
+def test_load_physio():
+ d = pk.Physio(np.random.randn(40), fs=10)
+ d = pk.peakfind_physio(d)
+ n = 'kaylee.phys'
+ pk.save_physio(n, d)
+ p = io.load_physio(n)
+
+ assert (p[0] == d.data).all()
+ assert (p[1] == d.peaks).all()
+ assert p[2] == 10
+ os.remove(n)
+
+
def test_export_regressor(testdir):
- petco2hrf_shift = np.random.rand(10)
+ petco2hrf_shift = np.random.randn(10)
ntp = 10
outname = os.path.join(testdir, 'test_regressor')
suffix = 'petco2hrf'
@@ -49,7 +124,7 @@ def test_export_regressor(testdir):
def test_export_nifti(testdir):
# create some test data
- data = np.random.rand(10, 10, 10)
+ data = np.random.randn(10, 10, 10)
affine = np.eye(4)
header = nib.Nifti1Header()
img = nib.Nifti1Image(data, affine, header)
@@ -66,8 +141,39 @@ def test_export_nifti(testdir):
assert np.allclose(out_img.get_fdata(), data, atol=1e-6, rtol=0)
assert np.allclose(out_img.affine, affine, atol=1e-6, rtol=0)
assert out_img.header.__dict__.keys() == header.__dict__.keys()
+ os.remove(fname)
# Eventually check that headers have the same content (although they don't now
# for good reasons!)
# ## Break tests
+def test_break_load_mat():
+ """Break load_mat."""
+ sys.modules['pymatreader'] = None
+ with pytest.raises(ImportError) as errorinfo:
+ io.load_mat('simon')
+ assert 'is required' in str(errorinfo.value)
+ sys.modules['pymatreader'] = pymatreader
+
+ a = 'heart'
+ n = 'ofgold'
+ savemat(n, {'data': a})
+
+ with pytest.raises(EOFError) as errorinfo:
+ io.load_mat(n)
+ assert f'{n} does not seem' in str(errorinfo.value)
+ os.remove(n)
+
+
+def test_break_load_xls():
+ """Break load_xls."""
+ with pytest.raises(NotImplementedError) as errorinfo:
+ io.load_xls('firefly')
+ assert 'loading is not' in str(errorinfo.value)
+
+
+def test_break_load_array():
+ """Break load_xls."""
+ with pytest.raises(NotImplementedError) as errorinfo:
+ io.load_array('firefly.leaf')
+ assert 'file extension' in str(errorinfo.value)
diff --git a/phys2cvr/tests/test_regressors.py b/phys2cvr/tests/test_regressors.py
new file mode 100644
index 0000000..3e66fa6
--- /dev/null
+++ b/phys2cvr/tests/test_regressors.py
@@ -0,0 +1,116 @@
+#!/usr/bin/env python3
+"""Tests for io."""
+
+import os
+from unittest.mock import patch
+
+import numpy as np
+import pytest
+
+from phys2cvr import regressors
+
+# ## Unit tests
+
+
+def test_create_legendre():
+ L = regressors.create_legendre(3, 10)
+ assert L.shape == (10, 4)
+
+ L = regressors.create_legendre(1, 5)
+ assert np.allclose(L[:, 1], np.linspace(-1, 1, 5))
+
+ L = regressors.create_legendre(0, 7)
+ assert L.shape == (7, 1)
+ assert np.all(L[:, 0] == 1)
+
+
+@patch(
+ 'phys2cvr.regressors.endtidal_interpolation', return_value=np.array([1, 2, 3, 4])
+)
+@patch('phys2cvr.regressors.plot_two_timeseries')
+@patch('phys2cvr.regressors.convolve_signal', return_value=np.array([9, 8, 7, 6]))
+def test_compute_petco2hrf(mock_conv, mock_plot, mock_endtidal, testdir):
+ out = os.path.join(testdir, 'Lune')
+ co2 = np.array([1, 2, 3, 4])
+ pidx = np.array([0, 2])
+ r = regressors.compute_petco2hrf(co2, pidx, 1.0, out)
+ assert np.all(r == np.array([9, 8, 7, 6]))
+ os.remove(f'{out}_petco2.1D')
+
+
+def test_compute_petco2hrf_skip(testdir):
+ out = os.path.join(testdir, 'Catherine')
+ r = regressors.compute_petco2hrf(
+ np.array([1, 2, 3]),
+ np.array([0]),
+ 1.0,
+ out,
+ comp_endtidal=False,
+ response_function=None,
+ )
+ assert np.allclose(r, np.array([-1, 0, 1])) # returns petco2 demeaned
+
+
+@patch('phys2cvr.regressors.x_corr', return_value=(None, 3, np.array([1, 2, 3])))
+@patch('phys2cvr.regressors.plot_xcorr')
+def test_compute_bulk_shift(mock_plotx, mock_xcorr, testdir):
+ out = os.path.join(testdir, 'Monoco')
+ func = np.arange(10.0)
+ pet = np.arange(20.0)
+ s = regressors.compute_bulk_shift(func, pet, 1.0, out)
+ assert s == 3
+ assert os.path.exists(f'{out}_optshift.1D')
+ os.remove(f'{out}_optshift.1D')
+
+
+@patch('phys2cvr.regressors.export_regressor', return_value=np.zeros((5, 5)))
+def test_create_fine_shift_regressors(mock_export, testdir):
+ out = os.path.join(testdir, 'Esquie')
+ pet = np.arange(50.0)
+ r = regressors.create_fine_shift_regressors(pet, 3, 2, 1.0, 10, 20, out)
+ assert r.shape == (5, 5)
+
+
+@patch('phys2cvr.regressors.export_regressor', return_value=np.zeros((3, 3)))
+def test_create_fine_shift_regressors_padding(mock_export, testdir):
+ out = os.path.join(testdir, 'Maelle')
+ pet = np.arange(10.0)
+ r = regressors.create_fine_shift_regressors(pet, 9, 4, 1.0, 5, 8, out)
+ assert r.shape == (3, 3)
+
+
+@patch('phys2cvr.regressors.resample_signal_freqs', return_value=np.arange(20.0))
+@patch('phys2cvr.regressors.plot_two_timeseries')
+@patch('phys2cvr.regressors.export_regressor', return_value=np.arange(10.0))
+def test_create_physio_regressor(mock_exp, mock_plot, mock_bs, testdir):
+ out = os.path.join(testdir, 'Sciel')
+ func = np.arange(10.0)
+ pet = np.arange(30.0)
+ d, lags = regressors.create_physio_regressor(
+ func, pet, 1.0, 1.0, out, lag_max=2, skip_xcorr=True
+ )
+ assert d.shape == (10,)
+ assert lags is not None
+
+
+@patch('phys2cvr.regressors.resample_signal_freqs', return_value=np.arange(20.0))
+@patch('phys2cvr.regressors.export_regressor', return_value=np.arange(10.0))
+def test_create_physio_regressor_no_lag_max(mock_exp, mock_bs, testdir):
+ out = os.path.join(testdir, 'Francois')
+ func = np.arange(10.0)
+ pet = np.arange(30.0)
+ d, lags = regressors.create_physio_regressor(
+ func, pet, 1.0, 1.0, out, lag_max=None, skip_xcorr=True
+ )
+ assert d.shape == (10,)
+ assert lags is None
+ os.remove(f'{out}_petco2hrf_vs_avgroi.png')
+
+
+# ## Break tests
+
+
+def test_break_compute_petco2hrf():
+ with pytest.raises(NotImplementedError) as errorinfo:
+ regressors.compute_petco2hrf(np.zeros((2, 2)), np.array([0]), 1.0, 'x')
+ assert 'Arrays with more' in str(errorinfo.value)
diff --git a/phys2cvr/tests/test_signal.py b/phys2cvr/tests/test_signal.py
index 84ae188..0c302c4 100644
--- a/phys2cvr/tests/test_signal.py
+++ b/phys2cvr/tests/test_signal.py
@@ -2,8 +2,11 @@
"""Tests for io."""
import os
+from pathlib import PosixPath
+from unittest.mock import patch
import numpy as np
+import pytest
from phys2cvr import signal
@@ -62,3 +65,99 @@ def test_create_hrf():
expected_peak_index = 25
actual_peak_index = np.argmax(hrf)
assert actual_peak_index == expected_peak_index
+
+
+def test_filter_signal():
+ x = np.sin(np.linspace(0, 10, 200))
+ out = signal.filter_signal(x, tr=1.0)
+ assert out.shape == x.shape
+
+ x = np.random.randn(5, 200)
+ out = signal.filter_signal(x, tr=1.0, axis=-1)
+ assert out.shape == x.shape
+
+
+def test_endtidal_interpolation():
+ sig = np.array([0, 1, 0, 2, 0])
+ peaks = [1, 3]
+ out = signal.endtidal_interpolation(sig, peaks)
+ assert out.shape == sig.shape
+ assert np.allclose(out[[1, 3]], sig[[1, 3]])
+
+ sig = np.arange(10.0)
+ peaks = [8, 2, 5]
+ out = signal.endtidal_interpolation(sig, peaks)
+ assert out.shape == sig.shape
+ assert np.allclose(out[[2, 5, 8]], sig[[2, 5, 8]])
+
+
+def test_convolve_signal(testdir):
+ x = np.ones(20)
+ out = signal.convolve_signal(x, freq=10, response_function='hrf')
+ assert out.ndim == 1
+ assert out.shape[0] > x.shape[0]
+ # Add assertion on first 20 points of hrf at 10 Hz)
+
+ x = np.arange(10.0)
+ out = signal.convolve_signal(x, freq=10, response_function=None)
+ assert np.allclose(out, x)
+
+ x = np.arange(5.0)
+ rf = [1, 0, -1]
+ out = signal.convolve_signal(x, freq=10, response_function=rf)
+ assert out.ndim == 1
+
+ p = os.path.join(testdir, 'rf.1D')
+ np.savetxt(p, np.array([1.0, 2.0, 3.0]))
+ x = np.ones(10)
+ with patch('phys2cvr.io.load_array', return_value=np.array([1.0, 2.0, 3.0])):
+ out = signal.convolve_signal(x, 10, response_function=str(p))
+ assert out.ndim == 1
+
+
+def test_resample_signal_samples_basic():
+ x = np.arange(10.0)
+ out = signal.resample_signal_samples(x, 20)
+ assert out.shape[0] == 20
+
+ x = np.random.randn(5, 10)
+ out = signal.resample_signal_samples(x, 15, axis=-1)
+ assert out.shape == (5, 15)
+
+
+def test_resample_signal_freqs_basic():
+ x = np.arange(10.0)
+ out = signal.resample_signal_freqs(x, 1, 2)
+ assert out.shape[0] > x.shape[0]
+
+ x = np.random.randn(4, 10)
+ out = signal.resample_signal_freqs(x, 1, 3, axis=-1)
+ assert out.shape[1] > 10
+
+
+# ## Break tests
+
+
+def test_break_convolve_signal_invalid_array(testdir):
+ x = np.arange(5.0)
+ rf = ['a', 'b']
+ with pytest.raises(ValueError):
+ signal.convolve_signal(x, freq=10, response_function=rf)
+
+ x = np.zeros((3, 3, 3))
+ with pytest.raises(NotImplementedError):
+ signal.convolve_signal(x, freq=10)
+
+ p = os.path.join(testdir, 'does_not_exist.1D')
+ x = np.ones(10)
+ with pytest.raises(OSError):
+ signal.convolve_signal(x, 10, response_function=PosixPath(p))
+
+ p = os.path.join(testdir, 'does_not_exist.1D')
+ x = np.ones(10)
+ with pytest.raises(NotImplementedError):
+ signal.convolve_signal(x, 10, response_function=str(p))
+
+ x = np.ones(10)
+ with pytest.raises(NotImplementedError):
+ signal.convolve_signal(x, 10, response_function='not_a_rf')
diff --git a/phys2cvr/tests/test_stats.py b/phys2cvr/tests/test_stats.py
index 8fbf474..cd06a8c 100644
--- a/phys2cvr/tests/test_stats.py
+++ b/phys2cvr/tests/test_stats.py
@@ -2,7 +2,160 @@
"""Tests for io."""
import numpy as np
+import pytest
from phys2cvr import stats
+
# ## Unit tests
+def test_x_corr():
+ func = np.array([1.0, -2.0, 3.0])
+ co2 = np.array([0.0, 1.0, -2.0, 3.0, 4.0])
+
+ # Manually compute expected correlations
+ # sco2 windows: [0,1,-2], [1,-2,3], [-2,3,4]
+ # zscores:
+ # func z: [-1.224745, 0, 1.224745]
+ # window z, then dot product / N
+ val, idx, arr = stats.x_corr(func, co2)
+
+ sco2 = np.array([[0.0, 1.0, -2.0], [1.0, -2.0, 3.0], [-2.0, 3.0, 4.0]])
+ # Pre-computed expected values
+ expected = np.corrcoef(func, sco2)[1:, 0]
+
+ assert isinstance(val, float)
+ assert isinstance(idx, (int | np.int64 | np.int32))
+ assert isinstance(arr, np.ndarray)
+ assert arr.shape[0] == co2.size - func.size + 1
+ assert np.allclose(arr, expected, atol=1e-6)
+ assert val == pytest.approx(expected.max())
+ assert idx == expected.argmax()
+
+ func = np.array([1.0, -2.0, 1.0])
+ co2 = np.array([1.0, -1.0, 2.0, -1.0, 1.0])
+ sco2 = np.array([[1, -1, 2], [-1, 2, -1], [2, -1, 1]])
+
+ val, idx, arr = stats.x_corr(func, co2, abs_xcorr=True)
+ expected_corr = np.corrcoef(func, sco2)[1:, 0]
+
+ assert val >= 0
+ assert arr.size == 3
+ assert np.allclose(arr, expected_corr)
+ assert val == pytest.approx(expected.max())
+ assert idx == 1
+
+
+def test_ols():
+ Y = np.array([1, 2, 3, 4], float)
+ X = np.column_stack([np.ones(4), np.arange(4)])
+ betas, t, r = stats.ols(Y, X)
+ assert betas.shape == (2, 1)
+ assert t.shape == (2, 1)
+ assert r.shape == (1,)
+ assert np.isfinite(r).all()
+
+ # Fit Y = b0*1 + b1*x → exact solution
+ X = np.column_stack([np.ones(5), np.arange(5)])
+ Y = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
+
+ betas, tstats, r = stats.ols(Y, X)
+
+ assert np.allclose(betas.flatten(), [0.0, 1.0])
+ assert r == pytest.approx(1.0)
+
+ # t-stats for perfect fit → huge (but finite after clipping)
+ assert tstats[1, 0] > 1000
+
+ X = np.column_stack([np.ones(4), np.arange(4)])
+ Y = np.array([1.0, 2.0, 1.0, 2.0]) # variance around intercept only
+
+ _, _, r = stats.ols(Y, X, r2model='intercept')
+
+ # r^2 intercept = R^2 of X vs intercept-only model
+ # Should be >0 because x explains alternating pattern
+ assert r[0] > 0.2
+
+
+def test_ols_residuals():
+ X = np.column_stack([np.ones(3), np.arange(3)])
+ Y = np.array([1.0, 3.0, 5.0])
+ # Perfect Y = 1 + 2x
+ res = stats.ols(Y, X, residuals=True)
+
+ assert res.shape == (3, 1)
+ assert np.allclose(res.flatten(), 0.0)
+
+
+def test_regression():
+ data = np.random.randn(4, 4, 4, 10)
+ regr = np.random.randn(10, 1)
+ b, t, r = stats.regression(data, regr)
+ assert b.shape == data.shape[:-1]
+ assert t.shape == data.shape[:-1]
+ assert r.shape == data.shape[:-1]
+
+ mask = np.zeros((4, 4, 4))
+ mask[1, 1, 1] = 1
+ b, t, r = stats.regression(data, regr, mask=mask)
+ assert b[1, 1, 1] != 0
+ assert b.sum() != 0
+ assert r.shape == data.shape[:-1]
+
+ data = np.zeros((1, 1, 1, 5))
+ data[0, 0, 0] = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
+
+ regr = np.arange(5)[:, None] # perfect predictor
+
+ b, t, r = stats.regression(data, regr)
+
+ assert b[0, 0, 0] == pytest.approx(1.0)
+ assert r[0, 0, 0] == pytest.approx(1.0)
+ assert t[0, 0, 0] > 1000
+
+ data = np.zeros((2, 2, 2, 4))
+ data[..., :] = np.array([1.0, 2.0, 3.0, 4.0])
+
+ regr = np.arange(4)[:, None]
+
+ mask = np.zeros((2, 2, 2))
+ mask[0, 0, 0] = 1
+
+ b, t, r = stats.regression(data, regr, mask=mask)
+
+ assert b[0, 0, 0] == pytest.approx(1.0)
+ assert r[0, 0, 0] == pytest.approx(1.0)
+ assert b.sum() == b[0, 0, 0]
+ assert r.sum() == r[0, 0, 0]
+
+
+# ## Break tests
+
+
+def test_break_x_corr():
+ func = np.ones(5)
+ co2 = np.ones(6)
+ with pytest.raises(ValueError):
+ stats.x_corr(func, co2, offset=5)
+
+ with pytest.raises(NotImplementedError):
+ stats.x_corr(func, co2, offset=-1)
+
+
+def test_break_ols():
+ Y = np.arange(5)
+ X = np.column_stack([np.ones(5), np.arange(5)])
+ with pytest.raises(ValueError):
+ stats.ols(Y, X, r2model='invalid')
+
+ Y = np.zeros((2, 2, 2))
+ X = np.ones((2, 3))
+ with pytest.raises(NotImplementedError):
+ stats.ols(Y, X)
+
+
+def test_break_regression():
+ data = np.random.randn(4, 4, 4, 10)
+ regr = np.random.randn(10, 1)
+ denoise = np.random.randn(9, 2)
+ with pytest.raises(ValueError):
+ stats.regression(data, regr, denoise_mat=denoise)
diff --git a/phys2cvr/tests/test_utils.py b/phys2cvr/tests/test_utils.py
index ad12bb0..148c15c 100644
--- a/phys2cvr/tests/test_utils.py
+++ b/phys2cvr/tests/test_utils.py
@@ -11,19 +11,18 @@
# ## Unit tests
-@pytest.mark.xfail(strict=False, reason="Something's off with lists")
@pytest.mark.parametrize(
'var, dtype, out',
[
(10, 'int', 10),
(10.0, 'int', 10),
('10', 'int', 10),
- ([10], 'int', [10]),
+ ([10], 'int', 10),
(None, 'int', None),
(10, 'float', 10.0),
(10.0, 'float', 10.0),
('10.0', 'float', 10.0),
- ([10.0], 'float', [10.0]),
+ ([10.0], 'float', 10.0),
(None, 'float', None),
(10, 'str', '10'),
(10.0, 'str', '10.0'),
@@ -108,6 +107,19 @@ def test_break_if_declared_force_type():
assert 'not supported' in str(errorinfo.value)
+# Check this one
+@pytest.mark.parametrize(
+ 'var, dtype',
+ [
+ ([10], 'int'),
+ ([10.0], 'float'),
+ ],
+)
+def test_if_declared_force_type_errors(var, dtype):
+ with pytest.raises(TypeError):
+ utils.if_declared_force_type(var, dtype)
+
+
@pytest.mark.xfail(strict=False, reason='Message changed')
def test_break_check_nifti_dim_missing_dims():
with pytest.raises(ValueError) as errorinfo:
diff --git a/phys2cvr/utils.py b/phys2cvr/utils.py
index 4754e73..8edf75b 100644
--- a/phys2cvr/utils.py
+++ b/phys2cvr/utils.py
@@ -246,7 +246,7 @@ def save_bash_call(fname, outdir):
"""
-Copyright 2021-2025, Stefano Moia & phys2cvr contributors.
+Copyright 2021-2026, Stefano Moia & phys2cvr contributors.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
diff --git a/phys2cvr/viz.py b/phys2cvr/viz.py
index 5c8f4db..980d176 100644
--- a/phys2cvr/viz.py
+++ b/phys2cvr/viz.py
@@ -111,7 +111,7 @@ def plot_xcorr(xcorr, outprefix, freq=None):
"""
-Copyright 2021-2025, Stefano Moia & phys2cvr contributors.
+Copyright 2021-2026, Stefano Moia & phys2cvr contributors.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
diff --git a/phys2cvr/workflows.py b/phys2cvr/workflows.py
index 03a0276..df463f8 100755
--- a/phys2cvr/workflows.py
+++ b/phys2cvr/workflows.py
@@ -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)
@@ -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!')
@@ -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!')
@@ -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:
@@ -427,6 +431,9 @@ def phys2cvr(
'file containing its peaks was provided. '
' Please provide peak file!'
)
+ else:
+ # pidx to None for compatibility
+ pidx = None
if freq is None:
raise NameError(
@@ -462,6 +469,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:
@@ -476,35 +516,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,
@@ -608,28 +619,16 @@ 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)
@@ -637,8 +636,8 @@ def phys2cvr(
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')
@@ -656,6 +655,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)
@@ -762,7 +771,7 @@ def _main(argv=None):
"""
-Copyright 2021-2025, Stefano Moia & phys2cvr contributors.
+Copyright 2021-2026, Stefano Moia & phys2cvr contributors.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.