diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 4a2ca8a..eae9e39 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,9 +1,9 @@ repos: - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.4.8 # Use the latest version + rev: v0.4.8 # Use the latest version hooks: - id: ruff - args: [--fix] # Optional: to enable lint fixes + args: [--fix] # Optional: to enable lint fixes - id: ruff-format - repo: https://github.com/pre-commit/pre-commit-hooks rev: v4.6.0 diff --git a/docs/examples/tissue/plot_two_cxm.py b/docs/examples/tissue/plot_two_cxm.py new file mode 100644 index 0000000..3e0639b --- /dev/null +++ b/docs/examples/tissue/plot_two_cxm.py @@ -0,0 +1,51 @@ +""" +============= +The Two Compartment Exchange Model +============= + +Simulating tissue concentrations from two compartment models with different settings. +""" + +import matplotlib.pyplot as plt + +# %% +# Import necessary packages +import numpy as np +import osipi + +# %% +# Generate Parker AIF with default settings. + +# Define time points in units of seconds - in this case we use a time +# resolution of 1 sec and a total duration of 6 minutes. +t = np.arange(0, 5 * 60, 1, dtype=float) + +# Create an AIF with default settings +ca = osipi.aif_parker(t) + +# %% +# Plot the tissue concentrations for an extracellular volume fraction +# of 0.2, plasma volume fraction of 0.1, permeability serface area of 5 ml/min +# and flow rate of 10 ml/min +PS = 15 # Permeability surface area product in ml/min +Fp = [5, 25] # Flow rate in ml/min +ve = 0.1 # Extracellular volume fraction +vp = [0.1, 0.02] # Plasma volume fraction + +ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp[0], PS=PS, ve=ve, vp=vp[0]) +plt.plot(t, ct, "b-", label=f" Fp = {Fp[0]},PS = {PS}, ve = {ve}, vp = {vp[0]}") + +ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp[1], PS=PS, ve=ve, vp=vp[0]) +plt.plot(t, ct, "r-", label=f" Fp = {Fp[1]},PS = {PS}, ve = {ve}, vp = {vp[0]}") + +ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp[0], PS=PS, ve=ve, vp=vp[1]) +plt.plot(t, ct, "g-", label=f" Fp = {Fp[0]},PS = {PS}, ve = {ve}, vp = {vp[1]}") + +ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp[1], PS=PS, ve=ve, vp=vp[1]) +plt.plot(t, ct, "y-", label=f" Fp = {Fp[1]},PS = {PS}, ve = {ve}, vp = {vp[1]}") + + +plt.xlabel("Time (sec)") +plt.ylabel("Tissue concentration (mM)") +plt.legend() +plt.show() diff --git a/docs/references/models/tissue_models/2cxm.md b/docs/references/models/tissue_models/2cxm.md new file mode 100644 index 0000000..68f2f7f --- /dev/null +++ b/docs/references/models/tissue_models/2cxm.md @@ -0,0 +1,17 @@ +# osipi.Two Compartment Exchange Model + +::: osipi.two_compartment_exchange_model + + +## Example using `osipi.two_compartment_exchange_model` + +
+ The Two Compartment Exchange Model +
+ + + The Two Compartment Exchange Model + + +
+
diff --git a/docs/references/models/tissue_models/index.md b/docs/references/models/tissue_models/index.md index 298b575..5584f43 100644 --- a/docs/references/models/tissue_models/index.md +++ b/docs/references/models/tissue_models/index.md @@ -1,3 +1,4 @@ - [tofts](tofts.md) - [extended_tofts](extended_tofts.md) +- [two_cxm](2cxm.md) diff --git a/mkdocs.yml b/mkdocs.yml index a742a1c..183e016 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -75,6 +75,7 @@ nav: - references/models/tissue_models/index.md - osipi.tofts: references/models/tissue_models/tofts.md - osipi.extended_tofts: references/models/tissue_models/extended_tofts.md + - osipi.two_compartment_exchange: references/models/tissue_models/2cxm.md - Examples: generated/gallery diff --git a/src/osipi/__init__.py b/src/osipi/__init__.py index 454bb7e..ffe955b 100644 --- a/src/osipi/__init__.py +++ b/src/osipi/__init__.py @@ -7,7 +7,8 @@ from ._tissue import ( tofts, - extended_tofts + extended_tofts, + two_compartment_exchange_model ) from ._signal import ( diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index 372b41c..dbc3e9f 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -2,6 +2,7 @@ import numpy as np from scipy.interpolate import interp1d +from scipy.signal import convolve from ._convolution import exp_conv @@ -73,6 +74,22 @@ def tofts( >>> plt.plot(t, ca, "r", t, ct, "b") """ + + # Input validation + if t.size != ca.size: + raise ValueError(f"t and ca must have the same size, got {t.size} and {ca.size}") + + for param, name in [(Ktrans, "Ktrans"), (ve, "ve"), (Ta, "Ta")]: + if not (np.isscalar(param) and np.isreal(param)): + raise TypeError(f"{name} must be a numeric scalar value, got {type(param).__name__}") + + if Ktrans < 0: + raise ValueError(f"Ktrans must be non-negative, got {Ktrans}") + if not (0 <= ve <= 1): + raise ValueError(f"ve must be in range [0, 1], got {ve}") + if Ta < 0: + raise ValueError(f"Ta must be non-negative, got {Ta}") + if not np.allclose(np.diff(t), np.diff(t)[0]): warnings.warn( ("Non-uniform time spacing detected. Time array may be" " resampled."), @@ -237,6 +254,24 @@ def extended_tofts( >>> plt.plot(t, ca, "r", t, ct, "b") """ + # Input validation + if t.size != ca.size: + raise ValueError(f"t and ca must have the same size, got {t.size} and {ca.size}") + + for param, name in [(Ktrans, "Ktrans"), (ve, "ve"), (vp, "vp")]: + if not (np.isscalar(param) and np.isreal(param)): + raise TypeError(f"{name} must be a numeric scalar value, got {type(param).__name__}") + + if not (0 <= ve <= 1): + raise ValueError(f"ve must be in range [0, 1], got {ve}") + if not (0 <= vp <= 1): + raise ValueError(f"vp must be in range [0, 1], got {vp}") + if not (0 <= ve + vp <= 1): + raise ValueError(f"Sum of ve ({ve}) and vp ({vp}) exceeds 1") + if Ta < 0: + raise ValueError(f"Ta must be non-negative, got {Ta}") + if Ktrans < 0: + raise ValueError(f"Ktrans must be non-negative, got {Ktrans}") if not np.allclose(np.diff(t), np.diff(t)[0]): warnings.warn( @@ -332,3 +367,165 @@ def extended_tofts( ct = ct_func(t) return ct + + +def two_compartment_exchange_model( + t: np.ndarray, + ca: np.ndarray, + Fp: float, + PS: float, + ve: float, + vp: float, + Ta: float = 30.0, +) -> np.ndarray: + """Two compartment exchange model + + Tracer flows from the AIF to the blood plasma compartment; two-way leakage + between the plasma and extracellular compartments(EES) is permitted. + + Args: + t: 1D array of times(s). [OSIPI code is Q.GE1.004] + ca: Arterial concentrations in mM for each time point in t. [OSIPI code is Q.IC1.001.[a]] + Fp: Blood plasma flow rate into a unit tissue volume in ml/min. [OSIPI code is Q.PH1.002] + PS: Permeability surface area product in ml/min. [OSIPI code is Q.PH1.004] + ve: Extracellular volume fraction. [OSIPI code Q.PH1.001.[e]] + vp: Plasma volume fraction. [OSIPI code Q.PH1.001.[p]] + Ta: Arterial delay time, i.e., + difference in onset time between tissue curve and AIF in units of sec. [OSIPI code Q.PH1.007] + + Returns: + Ct: Tissue concentrations in mM + + See Also: + `extended_tofts` + + References: + - Lexicon url: https://osipi.github.io/OSIPI_CAPLEX/perfusionModels/#indicator-kinetic-models + - Lexicon code: M.IC1.009 + - OSIPI name: Two Compartment Exchange Model + - Adapted from contributions by: MJT_UoEdinburgh_UK + + Example: + Create an array of time points covering 6 min in steps of 1 sec, + calculate the Parker AIF at these time points, calculate tissue concentrations + using the Two Compartment Exchange model and plot the results. + + >>> import matplotlib.pyplot as plt + >>> import osipi + + Calculate AIF + + >>> t = np.arange(0, 6 * 60, 0.1) + >>> ca = osipi.aif_parker(t) + + Plot the tissue concentrations for an extracellular volume fraction + of 0.2, plasma volume fraction of 0.1, permeability serface area of 5 ml/min + and flow rate of 10 ml/min + + >>> PS = 5 # Permeability surface area product in ml/min + >>> Fp = 10 # Flow rate in ml/min + >>> ve = 0.2 # Extracellular volume fraction + >>> vp = 0.1 # Plasma volume fraction + >>> ct = osipi.two_compartment_exchange_model(t, ca, Fp, PS, ve, vp) + >>> plt.plot(t, ca, "r", t, ct, "g") + + """ + # Input validation + if t.size != ca.size: + raise ValueError(f"t and ca must have the same size, got {t.size} and {ca.size}") + + for param, name in [(Fp, "Fp"), (PS, "PS"), (ve, "ve"), (vp, "vp")]: + if not (np.isscalar(param) and np.isreal(param)): + raise TypeError(f"{name} must be a numeric scalar value, got {type(param).__name__}") + + if Fp <= 0: + raise ValueError(f"Fp must be positive, got {Fp}") + if PS < 0: + raise ValueError(f"PS must be non-negative, got {PS}") + if not (0 <= ve <= 1): + raise ValueError(f"ve must be in range [0, 1], got {ve}") + if not (0 <= vp <= 1): + raise ValueError(f"vp must be in range [0, 1], got {vp}") + if not (0 <= ve + vp <= 1): + raise ValueError(f"Sum of ve ({ve}) and vp ({vp}) exceeds 1") + + if vp == 0: + E = 1 - np.exp(-PS / Fp) + Ktrans = E * Fp + return tofts(t, ca, Ktrans, ve, Ta, discretization_method="conv") + + if not np.allclose(np.diff(t), np.diff(t)[0]): + warnings.warn( + ("Non-uniform time spacing detected. Time array may be" " resampled."), + stacklevel=2, + ) + + # Convert units + fp_per_s = Fp / (60.0 * 100.0) + ps_per_s = PS / (60.0 * 100.0) + + # Calculate the impulse response function + v = ve + vp + + # Mean transit time + T = v / fp_per_s + tc = vp / fp_per_s + te = ve / ps_per_s + + upsample_factor = 1 + n = t.size + n_upsample = (n - 1) * upsample_factor + 1 + t_upsample = np.linspace(t[0], t[-1], n_upsample) + tau_upsample = t_upsample - t[0] + + sig_p = ((T + te) + np.sqrt((T + te) ** 2 - 4 * tc * te)) / (2 * tc * te) + sig_n = ((T + te) - np.sqrt((T + te) ** 2 - 4 * tc * te)) / (2 * tc * te) + + # Calculate the impulse response function for the plasma compartment and EES + + irf_cp = ( + vp + * sig_p + * sig_n + * ( + (1 - te * sig_n) * np.exp(-tau_upsample * sig_n) + + (te * sig_p - 1.0) * np.exp(-tau_upsample * sig_p) + ) + / (sig_p - sig_n) + ) + + irf_ce = ( + ve + * sig_p + * sig_n + * (np.exp(-tau_upsample * sig_n) - np.exp(-tau_upsample * sig_p)) + / (sig_p - sig_n) + ) + + irf_cp[[0]] /= 2 + irf_ce[[0]] /= 2 + + dt = np.min(np.diff(t)) / upsample_factor + + if Ta != 0: + f = interp1d( + t, + ca, + kind="linear", + bounds_error=False, + fill_value=0, + ) + ca = (t > Ta) * f(t - Ta) + + # get concentration in plasma and EES + Cp = dt * convolve(ca, irf_cp, mode="full", method="auto")[: len(t)] + Ce = dt * convolve(ca, irf_ce, mode="full", method="auto")[: len(t)] + + t_upsample = np.linspace(t[0], t[-1], n_upsample) + + Cp = np.interp(t, t_upsample, Cp) + Ce = np.interp(t, t_upsample, Ce) + + # get tissue concentration + Ct = Cp + Ce + return Ct diff --git a/tests/test_tissue.py b/tests/test_tissue.py index 304cba9..3e5fb5f 100755 --- a/tests/test_tissue.py +++ b/tests/test_tissue.py @@ -2,130 +2,240 @@ import numpy as np import osipi +import pytest -def test_tissue_tofts(): - """1. +@pytest.fixture +def time_array(): + return np.arange(0, 6 * 60, 1) - Basic operation of the function - test that the peak tissue concentration is less than the peak - AIF - """ +@pytest.fixture +def aif(time_array): + return osipi.aif_parker(time_array) - t = np.linspace(0, 6 * 60, 360) - ca = osipi.aif_parker(t) - ct = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2) - assert np.round(np.max(ct)) < np.round(np.max(ca)) - # 2. Basic operation of the function - test with non-uniform spacing of - # time array - t = np.geomspace(1, 6 * 60 + 1, num=360) - 1 - ca = osipi.aif_parker(t) - ct = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2) - assert np.round(np.max(ct)) < np.round(np.max(ca)) +# Tests for the Tofts model - # 3. The offset option - test that the tissue concentration is shifted - # from the AIF by the specified offset time - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2, Ta=60.0) - assert (np.min(np.where(ct > 0.0)) - np.min(np.where(ca > 0.0)) - 1) * 1 == 60.0 - # 4. Test that the discretization options give almost the same result - - # time step must be very small +@pytest.mark.parametrize( + "t", + [np.arange(0, 6 * 60, 1), np.geomspace(1, 6 * 60 + 1, num=360) - 1], + ids=["uniform", "non-uniform"], +) +def test_tofts_peak_concentration(t, aif): + """Test peak tissue concentration is less than peak AIF""" + ct = osipi.tofts(t, aif, Ktrans=0.6, ve=0.2) + assert np.round(np.max(ct)) < np.round(np.max(aif)) + + +def test_tofts_time_offset(time_array, aif): + """Test tissue concentration is shifted by specified offset time""" + ct = osipi.tofts(time_array, aif, Ktrans=0.6, ve=0.2, Ta=60.0) + offset_diff = (np.min(np.where(ct > 0.0)) - np.min(np.where(aif > 0.0)) - 1) * 1 + assert offset_diff == 60.0 + + +def test_tofts_discretization_agreement(): + """Test different discretization methods give similar results with small time steps""" t = np.arange(0, 6 * 60, 0.01) ca = osipi.aif_parker(t) ct_conv = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2) ct_exp = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2, discretization_method="exp") assert np.allclose(ct_conv, ct_exp, rtol=1e-4, atol=1e-3) - # 5. Test that the ratio of the area under the ct and ca curves is - # approximately the extracellular volume - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct_conv = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2) - ct_exp = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2, discretization_method="exp") - assert math.isclose(np.trapz(ct_conv, t) / np.trapz(ca, t), 0.2, abs_tol=1e-1) - assert math.isclose(np.trapz(ct_exp, t) / np.trapz(ca, t), 0.2, abs_tol=1e-1) - # 6. Test specific use cases - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct_conv = osipi.tofts(t, ca, Ktrans=0, ve=0.2) - assert np.count_nonzero(ct_conv) == 0 - - ct_exp = osipi.tofts(t, ca, Ktrans=0, ve=0.2, discretization_method="exp") - assert np.count_nonzero(ct_exp) == 0 - - ct_conv = osipi.tofts(t, ca, Ktrans=0.6, ve=0) - assert np.count_nonzero(ct_conv) == 0 - - ct_exp = osipi.tofts(t, ca, Ktrans=0.6, ve=0, discretization_method="exp") - assert np.count_nonzero(ct_exp) == 0 - - -def test_tissue_extended_tofts(): - # 1. Basic operation of the function - test that the peak tissue - # concentration is less than the peak AIF - t = np.linspace(0, 6 * 60, 360) +@pytest.mark.parametrize("method", [None, "exp"], ids=["conv", "exp"]) +def test_tofts_volume_ratio(time_array, aif, method): + """Test ratio of AUC matches extracellular volume""" + ct = osipi.tofts(time_array, aif, Ktrans=0.6, ve=0.2, discretization_method=method) + ratio = np.trapz(ct, time_array) / np.trapz(aif, time_array) + assert math.isclose(ratio, 0.2, abs_tol=1e-1) + + +@pytest.mark.parametrize("Ktrans, ve", [(0, 0.2), (0.6, 0)], ids=["Ktrans=0", "ve=0"]) +@pytest.mark.parametrize("method", [None, "exp"], ids=["conv", "exp"]) +def test_tofts_specific_cases(time_array, aif, Ktrans, ve, method): + """Test edge cases with parameter combinations""" + ct = osipi.tofts(time_array, aif, Ktrans=Ktrans, ve=ve, discretization_method=method) + assert np.count_nonzero(ct) == 0 + + +@pytest.mark.parametrize( + "param, value, match", + [ + ("Ktrans", -0.1, "non-negative"), + ("Ktrans", "invalid", "numeric scalar"), + ("ve", -0.2, "in range \[0, 1\]"), + ("ve", [0.2], "numeric scalar"), + ("Ta", -5.0, "non-negative"), + ], + ids=["Ktrans_negative", "Ktrans_type", "ve_negative", "ve_type", "Ta_negative"], +) +def test_tofts_invalid_parameters(time_array, aif, param, value, match): + """Test invalid parameter values/type for all model parameters""" + with pytest.raises((ValueError, TypeError), match=match): + # Construct kwargs with one invalid parameter at a time + kwargs = {"Ktrans": 0.6, "ve": 0.2, "Ta": 30.0} + kwargs[param] = value + osipi.tofts(time_array, aif, **kwargs) + + +# Tests for the Extended Tofts model + + +@pytest.mark.parametrize( + "t", + [np.linspace(0, 6 * 60, 360), np.geomspace(1, 6 * 60 + 1, num=360) - 1], + ids=["uniform", "non-uniform"], +) +def test_extended_tofts_peak_concentration(t): + """Test peak tissue concentration is less than peak AIF for different time arrays""" ca = osipi.aif_parker(t) ct = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3) assert np.round(np.max(ct)) < np.round(np.max(ca)) - # 2. Basic operation of the function - test with non-uniform spacing of - # time array - t = np.geomspace(1, 6 * 60 + 1, num=360) - 1 - ca = osipi.aif_parker(t) - ct = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3) - assert np.round(np.max(ct)) < np.round(np.max(ca)) - # 3. The offset option - test that the tissue concentration is shifted - # from the AIF by the specified offset time - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3, Ta=60.0) - assert (np.min(np.where(ct > 0.0)) - np.min(np.where(ca > 0.0)) - 1) * 1 == 60.0 +def test_extended_tofts_time_offset(time_array, aif): + """Test tissue concentration is shifted by specified offset time""" + ct = osipi.extended_tofts(time_array, aif, Ktrans=0.6, ve=0.2, vp=0.3, Ta=60.0) + offset_diff = (np.min(np.where(ct > 0.0)) - np.min(np.where(aif > 0.0)) - 1) * 1 + assert offset_diff == 60.0 - # 4. Test that the discretization options give almost the same result - - # time step must be very small + +def test_extended_tofts_discretization_agreement(): + """Test different discretization methods give similar results with small time steps""" t = np.arange(0, 6 * 60, 0.01) ca = osipi.aif_parker(t) ct_conv = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3) ct_exp = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3, discretization_method="exp") assert np.allclose(ct_conv, ct_exp, rtol=1e-4, atol=1e-3) - # 5. Test that the ratio of the area under the ct and ca curves is - # approximately the extracellular volume plus the plasma volume - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct_conv = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3) - ct_exp = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3, discretization_method="exp") - assert math.isclose( - np.trapz(ct_conv, t) / np.trapz(ca, t), - 0.2 + 0.3, - abs_tol=1e-1, +@pytest.mark.parametrize("method", [None, "exp"], ids=["conv", "exp"]) +def test_extended_tofts_volume_ratio(time_array, aif, method): + """Test ratio of AUC matches ve + vp for different discretization methods""" + ct = osipi.extended_tofts( + time_array, aif, Ktrans=0.6, ve=0.2, vp=0.3, discretization_method=method ) - assert math.isclose(np.trapz(ct_exp, t) / np.trapz(ca, t), 0.2 + 0.3, abs_tol=1e-1) + ratio = np.trapz(ct, time_array) / np.trapz(aif, time_array) + assert math.isclose(ratio, 0.2 + 0.3, abs_tol=1e-1) - # 6. Test specific use cases - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct_conv = osipi.extended_tofts(t, ca, Ktrans=0, ve=0.2, vp=0.3) - assert np.allclose(ct_conv, ca * 0.3, rtol=1e-4, atol=1e-3) - ct_exp = osipi.extended_tofts(t, ca, Ktrans=0, ve=0.2, vp=0.3, discretization_method="exp") - assert np.allclose(ct_conv, ca * 0.3, rtol=1e-4, atol=1e-3) +@pytest.mark.parametrize("Ktrans, ve", [(0, 0.2), (0.6, 0)], ids=["Ktrans=0", "ve=0"]) +@pytest.mark.parametrize("method", [None, "exp"], ids=["conv", "exp"]) +def test_extended_tofts_specific_cases(time_array, aif, Ktrans, ve, method): + """Test edge cases with parameter combinations""" + ct = osipi.extended_tofts( + time_array, aif, Ktrans=Ktrans, ve=ve, vp=0.3, discretization_method=method + ) + expected = aif * 0.3 + assert np.allclose(ct, expected, rtol=1e-4, atol=1e-3) + + +@pytest.mark.parametrize( + "param, value, match", + [ + ("Ktrans", -0.1, "non-negative"), + ("Ktrans", "invalid", "numeric scalar"), + ("ve", -0.2, "in range \[0, 1\]"), + ("ve", [0.2], "numeric scalar"), + ("vp", -0.2, "in range \[0, 1\]"), + ("vp", [0.2], "numeric scalar"), + ("vp", 1.1, "in range \[0, 1\]"), + ("vp", [0.2], "numeric scalar"), + ("Ta", -5.0, "non-negative"), + ], + ids=[ + "Ktrans_negative", + "Ktrans_type", + "ve_negative", + "ve_type", + "vp_negative", + "vp_type", + "vp_range", + "vp_type", + "Ta_negative", + ], +) +def test_extended_tofts_invalid_parameters(time_array, aif, param, value, match): + """Test invalid parameter values/type for all model parameters""" + with pytest.raises((ValueError, TypeError), match=match): + # Construct kwargs with one invalid parameter at a time + kwargs = {"Ktrans": 0.6, "ve": 0.2, "vp": 0.3, "Ta": 30.0} + kwargs[param] = value + osipi.extended_tofts(time_array, aif, **kwargs) + + +# Tests for the 2CXM model + + +@pytest.mark.parametrize( + "time_points", + [ + (np.linspace(0, 6 * 60, 360), "uniform"), + (np.geomspace(1, 6 * 60 + 1, num=360) - 1, "non-uniform"), + ], +) +def test_2cxm_basic_operation(time_points): + """Test that the peak tissue concentration is less than the peak AIF for different time arrays""" + t, _ = time_points + ca = osipi.aif_parker(t) + ct = osipi.two_compartment_exchange_model(t, ca, Fp=10, PS=5, ve=0.2, vp=0.3) + assert np.round(np.max(ct)) < np.round(np.max(ca)) - ct_conv = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0, vp=0.3) - assert np.allclose(ct_conv, ca * 0.3, rtol=1e-4, atol=1e-3) - ct_exp = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0, vp=0.3, discretization_method="exp") - assert np.allclose(ct_conv, ca * 0.3, rtol=1e-4, atol=1e-3) +def test_2cxm_time_offset(time_array, aif): + """Test that the tissue concentration is shifted from the AIF by the specified offset time""" + ct = osipi.two_compartment_exchange_model(time_array, aif, Fp=10, PS=5, ve=0.2, vp=0.3, Ta=60.0) + assert (np.min(np.where(ct > 0.0)) - np.min(np.where(aif > 0.0)) - 1) * 1 == 60.0 + + +def test_2cxm_zero_vp_matches_tofts(time_array, aif): + """Test that 2CXM with vp=0 behaves like the Tofts model""" + ct = osipi.two_compartment_exchange_model(time_array, aif, Fp=10, PS=5, ve=0.2, vp=0) + ct_tofts = osipi.tofts(time_array, aif, Ktrans=3.93, ve=0.2) + assert np.allclose(ct, ct_tofts, rtol=1e-4, atol=1e-3) + + +@pytest.mark.parametrize( + "param, value, match", + [ + ("Fp", [10], "numeric scalar"), # Non-scalar Fp + ("PS", "five", "numeric scalar"), # Non-numeric PS + ("ve", np.array([0.2]), "numeric scalar"), # Array ve + ("vp", None, "numeric scalar"), # None value + ("Fp", -5, "positive"), # Negative Fp + ("PS", -1, "non-negative"), # Negative PS + ("ve", -0.1, "in range \[0, 1\]"), # ve < 0 + ("vp", 1.1, "in range \[0, 1\]"), # vp > 1 + (("ve", "vp"), (0.8, 0.3), "Sum of ve \(0.8\) and vp \(0.3\) exceeds 1"), + ], + ids=[ + "Fp_type", + "PS_type", + "ve_type", + "vp_type", + "Fp_negative", + "PS_negative", + "ve_range", + "vp_range", + "ve_vp_sum", + ], +) +def test_2cxm_invalid_inputs(time_array, aif, param, value, match): + """Test invalid parameter values/type for all model parameters""" + with pytest.raises((ValueError, TypeError), match=match): + # Construct kwargs with default valid values + kwargs = {"Fp": 10, "PS": 5, "ve": 0.2, "vp": 0.1} + # Handle tuple parameters for combined validation case + if isinstance(param, tuple): + for p, v in zip(param, value): + kwargs[p] = v + else: + kwargs[param] = value + osipi.two_compartment_exchange_model(time_array, aif, **kwargs) if __name__ == "__main__": - test_tissue_tofts() - test_tissue_extended_tofts() - - print("All tissue concentration model tests passed!!") + pytest.main([__file__])