diff --git a/phys2cvr/io.py b/phys2cvr/io.py index 60cf05f..2b72c03 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 diff --git a/phys2cvr/regressors.py b/phys2cvr/regressors.py index 9b42de7..ca680e1 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 ' diff --git a/phys2cvr/signal.py b/phys2cvr/signal.py index b290c85..4c30b5e 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] 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/workflows.py b/phys2cvr/workflows.py index 79f1286..28d733b 100755 --- a/phys2cvr/workflows.py +++ b/phys2cvr/workflows.py @@ -427,6 +427,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( @@ -451,7 +454,7 @@ def phys2cvr( _, basename_co2, _ = utils.check_ext( io.EXT_ALL, os.path.basename(fname_co2), scan=True, remove=True - )[1] + ) outprefix = os.path.join(outdir, basename_co2)