expconv

expconv

Exponential convolution functions for pharmacokinetic modeling.

This module implements efficient recursive formulas for convolving signals with exponential and multi-exponential functions. These are the core building blocks for compartmental pharmacokinetic models.

GPU/CPU agnostic using the xp array module pattern. NO scipy dependency - see XP Compatibility Requirements in plan.md.

References

Flouri D, Lesnic D, Mayrovitz HN (2016). Numerical solution of the convolution integral equations in pharmacokinetics. Comput Methods Biomech Biomed Engin. 19(13):1359-1367.

dcmri library: https://github.com/dcmri/dcmri

expconv

expconv(f, T, t)

Convolve a signal with exponential decay function(s).

Computes the convolution: (f * exp(-t/T))(t) = integral_0^t f(u) * exp(-(t-u)/T) du

using an efficient recursive formula that avoids explicit numerical integration. This is the fundamental operation for compartmental pharmacokinetic models.

Handles both single-voxel (scalar T) and batch (array T) cases. When T is an array, every voxel is processed in a single pass with the loop running over time points — efficient on CPU and GPU.

PARAMETER	DESCRIPTION
f	Input signal (e.g., arterial input function). Shape (n_times,) or (n_times, 1) or (n_times, n_voxels). TYPE: ndarray
T	Time constant(s) of the exponential decay in seconds. Scalar for a single convolution, or array of shape (n_voxels,) for batch convolution. TYPE: float or ndarray
t	Time points in seconds. Shape (n_times,) or (n_times, 1). Must be monotonically increasing. TYPE: ndarray

RETURNS	DESCRIPTION
ndarray	Convolution result. Shape (n_times,) for scalar T, (n_times, n_voxels) for array T.

Notes

The recursive formula from Flouri et al. (2016) is used:

E[i] = E[i-1] * exp(-dt/T) + integral_{t[i-1]}^{t[i]} f(u) * exp(-(t[i]-u)/T) du

where the integral within each interval is evaluated analytically assuming piecewise-linear interpolation of f.

This is O(n) in computation time, compared to O(n^2) for naive numerical integration.

Examples:

>>> import numpy as np
>>> from osipy.common.convolution import expconv
>>> t = np.linspace(0, 10, 101)
>>> aif = np.exp(-t / 2)  # Input function
>>> T = 5.0  # 5 second time constant
>>> result = expconv(aif, T, t)

References

.. [1] Flouri D, Lesnic D, Mayrovitz HN (2016). Comput Methods Biomech Biomed Engin. 19(13):1359-1367.

biexpconv

biexpconv(f, T1, T2, t)

Convolve a signal with a bi-exponential function.

Computes the convolution of f with the bi-exponential function: h(t) = (exp(-t/T1) - exp(-t/T2)) / (T1 - T2)

This is useful for two-compartment exchange models.

PARAMETER	DESCRIPTION
f	Input signal. Shape: (n_times,). TYPE: ndarray
T1	First time constant in seconds. Must be positive. TYPE: float
T2	Second time constant in seconds. Must be positive and different from T1. TYPE: float
t	Time points in seconds. Shape: (n_times,). TYPE: ndarray

RETURNS	DESCRIPTION
ndarray	Convolution result. Shape: (n_times,).

Notes

The bi-exponential convolution is computed as: biexpconv(f, T1, T2, t) = (expconv(f, T1, t) - expconv(f, T2, t)) / (T1 - T2)

For T1 -> T2, uses the limiting form involving the derivative.

Examples:

>>> import numpy as np
>>> from osipy.common.convolution import biexpconv
>>> t = np.linspace(0, 10, 101)
>>> aif = np.exp(-t / 2)
>>> result = biexpconv(aif, 3.0, 5.0, t)

References

.. [1] dcmri library: https://github.com/dcmri/dcmri

nexpconv

nexpconv(f, T, n_exp, t)

Convolve a signal with an n-exponential (gamma variate) function.

Computes the convolution of f with the n-exponential function: h(t) = (t/T)^(n-1) * exp(-t/T) / (T * (n-1)!)

This represents a chain of n identical compartments, each with time constant T. Also known as the gamma variate function.

PARAMETER	DESCRIPTION
f	Input signal. Shape: (n_times,). TYPE: ndarray
T	Time constant of each compartment in seconds. Must be positive. TYPE: float
n_exp	Number of exponentials (compartments). Must be >= 1. TYPE: int
t	Time points in seconds. Shape: (n_times,). TYPE: ndarray

RETURNS	DESCRIPTION
ndarray	Convolution result. Shape: (n_times,).

Notes

For n=1, this is equivalent to expconv. For n>1, the function is computed recursively: nexpconv(f, T, n, t) = expconv(nexpconv(f, T, n-1, t), T, t) / T

For large n, the gamma variate approaches a Gaussian centered at t = nT with standard deviation sqrt(n)T.

For n > 20, numerical overflow may occur in the factorial term. In this case, a Gaussian approximation is used.

Examples:

>>> import numpy as np
>>> from osipy.common.convolution import nexpconv
>>> t = np.linspace(0, 20, 201)
>>> aif = np.exp(-t / 2)
>>> result = nexpconv(aif, 2.0, 3, t)  # 3-compartment chain

References

.. [1] dcmri library: https://github.com/dcmri/dcmri