Merge pull request #1 from smoia/fix/pidxcall

Fix lack of provided peaks when not necessary

Author: beccaclements99

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
 
 

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 '

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]
 

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)

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)