"""
Temporal dynamics functions for time-resolved spectroscopy.
Function Conventions
--------------------
Use CamelCase naming (UpperCamelCase or lowerCamelCase) for function names.
**Dynamics Functions:**
Signature: func(t, par1, par2, ..., t0)
- t: Time axis (numpy array)
- par1, par2, ...: Function-specific parameters
- t0: Time zero (function starts at this time)
- Returns: f(t) = 0 for t < t0, dynamics for t >= t0
**Convolution Kernels:**
Signature: funcCONV(t, par1, par2, ...)
- t: Time differences centered at zero. Kernels must be elementwise in
this argument: the kernel-matrix convolution evaluates them on a 2D
matrix of time differences (see utils.arrays.conv_matrix_operator)
- par1, par2, ...: Kernel parameters
- Returns: Kernel values (unnormalized; the convolution row-normalizes)
Every kernel additionally requires a private edge-mass companion
``_<name>_edge_mass(dt_left, dt_right, **params)`` registered in
``CONV_EDGE_MASS``. It returns the exact exterior masses of the kernel
body (same normalization) beyond the data window — the analytic
integrals over ``(-inf, t[0]]`` and ``[t[-1], inf)`` used for
edge-value padding. Use cancellation-safe tail forms (erfc, direct
exp, clip) rather than CDF differences, and validate the parameters
(strictly positive, finite) — companions are the runtime backstop for
expression-driven kernel parameters that bypass model-load bound
checks. A kernel without a companion is rejected at model validation
(mcp) and scheduling (GIR).
**Time Zero Convention:**
All dynamics functions are zero before t0 and activate at t >= t0.
This reflects physical causality: response begins after excitation.
Exception: erfFun is a smoothed step, so it is nonzero (approaching 0)
shortly before t0; it crosses A/2 at t0.
**Offsets and Baselines:**
Dynamics components combine by addition, so constant offsets are their
own component rather than a parameter of every function:
- stepFun: sharp onset to a constant value A at t0
- erfFun: Gaussian-broadened onset, rises from 0 to A around t0
For example, a decay to a nonzero plateau is expFun + stepFun sharing t0.
**Time Resolution:**
Functions inherit time axis from Dynamics model. Consider:
- Time step size relative to dynamics (dt << tau)
- Time range coverage (include full decay/rise)
- Kernel width appropriate for convolution
Parameter Naming
----------------
Common parameter names:
- A: Amplitude (change in signal)
- tau: Time constant (decay/rise time, 1/e point)
- t0: Time zero (start of dynamics)
- f: Frequency (for oscillations)
- phi: Phase (for oscillations)
- SD: Standard deviation (for Gaussian kernels)
- W: FWHM (for Lorentzian kernels)
Adding New Functions
--------------------
To add a new dynamics or convolution function:
1. Implement following conventions above
2. Ensure f(t<t0) = 0 for dynamics functions
3. Keep convolution kernels elementwise in their first argument
4. Test with realistic time-resolved data
"""
from collections.abc import Callable
import numpy as np
from scipy.special import erf, erfc
#
[docs]
def none(t: np.ndarray) -> np.ndarray:
"""
Placeholder function to define empty subcycles in a mcp.Dynamics model.
Used to define empty subcycles in multi-cycle Dynamics models without
adding any time-dependent behavior. This allows subcycle numbering to
work correctly when some subcycles should have no dynamics.
Usage (in model YAML file)::
model_sub2:
none: {}
Parameters
----------
t : ndarray
Time axis (not used)
Returns
-------
ndarray
Array of zeros with same shape as t
"""
# This function should never actually be called:
# It is caught in mcp.Model.combine() and skipped entirely.
return np.zeros_like(t)
#
[docs]
def stepFun(t: np.ndarray, A: float, t0: float) -> np.ndarray:
"""
Step function (constant offset switching on at t0).
The causal offset primitive: dynamics components combine by addition,
so add stepFun to model baselines or plateaus (e.g. expFun + stepFun
sharing t0 gives a decay to a nonzero plateau). For a Gaussian-broadened
onset use erfFun instead.
Parameters
----------
t : ndarray
Time axis
A : float
Step height (constant value for t >= t0)
t0 : float
Time zero (onset of the step)
Returns
-------
ndarray
Step function: 0 for t<t0, A for t>=t0
"""
return np.where(t < t0, 0.0, A)
#
[docs]
def linFun(t: np.ndarray, m: float, t0: float) -> np.ndarray:
"""
Linear dynamics (constant rate of change).
Parameters
----------
t : ndarray
Time axis
m : float
Slope (rate of change). Units: [signal units]/[time units]
- m > 0: Linear increase
- m < 0: Linear decrease
t0 : float
Time zero (start of linear change)
Returns
-------
ndarray
Linear function: 0 for t<t0, m*(t-t0) for t>=t0
"""
return np.where(t < t0, 0.0, m * (t - t0))
#
[docs]
def expFun(t: np.ndarray, A: float, tau: float, t0: float) -> np.ndarray:
"""
Exponential decay or rise dynamics.
Parameters
----------
t : ndarray
Time axis
A : float
Amplitude (initial change at t0).
- A > 0: Jumps to A at t0, decays toward 0
- A < 0: Jumps to ``-|A|`` at t0, rises toward 0
tau : float
Time constant (1/e time). Units: [time units]
At t = t0 + tau, signal changes by factor of e (≈2.718)
t0 : float
Time zero (start of exponential)
Returns
-------
ndarray
Exponential: 0 for t<t0, A*exp(-(t-t0)/tau) for t>=t0
"""
return np.where(t < t0, 0.0, A * np.exp(-1 / tau * (t - t0)))
#
[docs]
def sinFun(t: np.ndarray, A: float, f: float, phi: float, t0: float) -> np.ndarray:
"""
Sinusoidal oscillations (coherent dynamics).
Parameters
----------
t : ndarray
Time axis
A : float
Oscillation amplitude (peak-to-peak = 2A)
f : float
Frequency in [1/time units]
Period = 1/f
phi : float
Phase offset in radians
- phi = 0: Sine starts at zero
- phi = π/2: Starts at maximum (cosine)
- phi = π: Starts at zero (negative slope)
t0 : float
Time zero (start of oscillation)
Returns
-------
ndarray
Sinusoid: 0 for t<t0, A*sin(2πf(t-t0)+phi) for t>=t0
Oscillates around 0; add stepFun to shift the center line.
"""
return np.where(t < t0, 0.0, A * np.sin(2 * np.pi * f * (t - t0) + phi))
#
[docs]
def sinDivX(t: np.ndarray, A: float, f: float, t0: float) -> np.ndarray:
"""
Damped sinc function: sin(x)/x oscillation.
Parameters
----------
t : ndarray
Time axis
A : float
Amplitude scaling factor
f : float
Frequency in [1/time units]
t0 : float
Time zero (start of oscillation)
Returns
-------
ndarray
Sinc oscillation: 0 for t<t0, A*sin(2πf(t-t0))/(2πf(t-t0)) for t>=t0
"""
# np.sinc(u) = sin(pi*u)/(pi*u), so u=2*f*(t-t0) gives sin(2*pi*f*dt)/(2*pi*f*dt)
return np.where(t < t0, 0.0, A * np.sinc(2 * f * (t - t0)))
#
[docs]
def erfFun(t: np.ndarray, A: float, SD: float, t0: float) -> np.ndarray:
"""
Error function rise (step with Gaussian broadening).
erfFun ≈ stepFun ⊗ Gaussian(SD)
As a smoothed step this is the one dynamics function that is nonzero
(approaching 0) shortly before t0; it crosses A/2 at t0 and rises to A.
Parameters
----------
t : ndarray
Time axis
A : float
Amplitude (final value, asymptote as t → ∞)
SD : float
Standard deviation of Gaussian broadening (rise time ~2.355*SD)
Smaller SD → sharper rise
t0 : float
Center of rise (50% point)
Returns
-------
ndarray
Error function: A/2 * (1 + erf((t-t0)/(SD*√2)))
"""
return np.asarray(A / 2 * (1 + erf((t - t0) / (SD * np.sqrt(2)))))
#
[docs]
def sqrtFun(t: np.ndarray, A: float, t0: float) -> np.ndarray:
"""
Square root rise (diffusion dynamics).
Parameters
----------
t : ndarray
Time axis
A : float
Amplitude scaling factor
t0 : float
Time zero (start of diffusion)
Returns
-------
ndarray
Square root rise: 0 for t<t0, A*√(t-t0) for t>=t0
"""
# numpy array .clip sets all t<t0 to zero
return np.asarray(A * np.sqrt((t - t0).clip(0)))
#
# convolution kernels
# elementwise in the time-difference argument (evaluated on the dt
# matrix of the kernel-matrix convolution operator)
#
#
[docs]
def gaussCONV(x: np.ndarray, SD: float) -> np.ndarray:
"""
Gaussian convolution kernel (instrumental response function).
Parameters
----------
x : ndarray
Time differences (centered at 0)
SD : float
Standard deviation (Gaussian width).
FWHM = 2.355 * SD = 2*√(2ln2) * SD
Returns
-------
ndarray
Gaussian kernel (unnormalized, will be normalized in convolution)
"""
return np.exp(-1 / 2 * (x / SD) ** 2)
#
[docs]
def expSymCONV(x: np.ndarray, tau: float) -> np.ndarray:
"""
Symmetric exponential kernel (double exponential).
Exponential decay in both directions from center:
``exp(-|x|/tau)``
Parameters
----------
x : ndarray
Time axis (centered at 0)
tau : float
Decay time constant
Returns
-------
ndarray
Symmetric exponential kernel
"""
return np.asarray(np.exp(-1 / tau * np.abs(x)))
#
[docs]
def expDecayCONV(x: np.ndarray, tau: float) -> np.ndarray:
"""
Causal exponential kernel (one-sided decay).
Parameters
----------
x : ndarray
Time axis (centered at 0)
tau : float
Decay time constant
Returns
-------
ndarray
One-sided exponential: 0 for x<0, exp(-x/tau) for x≥0
"""
return np.where(x < 0, 0.0, expSymCONV(x, tau))
#
[docs]
def expRiseCONV(x: np.ndarray, tau: float) -> np.ndarray:
"""
Anti-causal exponential rise kernel (mirror of expDecayCONV).
The kernel is nonzero only for x <= 0, so the convolved response at
time t draws from the signal at later times: it rises before the
excitation and saturates at t0.
Parameters
----------
x : ndarray
Time axis (centered at 0)
tau : float
Rise time constant
Returns
-------
ndarray
One-sided exponential: exp(x/tau) for x≤0, 0 for x>0
"""
return np.where(x > 0, 0.0, expSymCONV(x, tau))
#
[docs]
def boxCONV(x: np.ndarray, width: float) -> np.ndarray:
"""
Box (rectangular) convolution kernel.
Parameters
----------
x : ndarray
Time axis (centered at 0)
width : float
Width of rectangular window
Returns
-------
ndarray
Rectangular function: 1 inside width, 0 outside (hard edges)
"""
return np.where(np.abs(x) <= width / 2, 1.0, 0.0)
#
# edge-mass companions
# exact exterior masses of each kernel body (same normalization as the
# body) for the edge-value padding of the kernel-matrix convolution:
# M_L = integral over (-inf, t[0]] = upper tail at dt_left = t - t[0],
# M_R = integral over [t[-1], inf) = lower tail at dt_right = t - t[-1].
# Companions validate their parameters (strictly positive, finite):
# they run on every evaluation, so they are the backstop for
# expression-driven kernel parameters that bypass the model-load bound
# checks (e.g. boxCONV width = 0 would otherwise silently become the
# identity operator). Private (excluded from registry discovery);
# dispatched via CONV_EDGE_MASS below.
#
#
def _validate_kernel_par(name: str, value: float) -> None:
"""Reject nonpositive or non-finite kernel parameters at evaluation."""
if not (np.isfinite(value) and value > 0):
raise ValueError(
f"Convolution kernel parameter '{name}' must be strictly "
f"positive and finite, got {value}."
)
#
def _gaussCONV_edge_mass(
dt_left: np.ndarray, dt_right: np.ndarray, SD: float
) -> tuple[np.ndarray, np.ndarray]:
"""Exterior masses of the peak-1 Gaussian body (erfc tail forms)."""
_validate_kernel_par("SD", SD)
scale = SD * np.sqrt(np.pi / 2)
M_L = scale * erfc(dt_left / (np.sqrt(2) * SD))
M_R = scale * erfc(-dt_right / (np.sqrt(2) * SD))
return np.asarray(M_L), np.asarray(M_R)
#
def _expSymCONV_edge_mass(
dt_left: np.ndarray, dt_right: np.ndarray, tau: float
) -> tuple[np.ndarray, np.ndarray]:
"""Exterior masses of exp(-|x|/tau); dt_left >= 0 and dt_right <= 0."""
_validate_kernel_par("tau", tau)
M_L = tau * np.exp(-dt_left / tau)
M_R = tau * np.exp(dt_right / tau)
return np.asarray(M_L), np.asarray(M_R)
#
def _expDecayCONV_edge_mass(
dt_left: np.ndarray, dt_right: np.ndarray, tau: float
) -> tuple[np.ndarray, np.ndarray]:
"""Exterior masses of the causal kernel: zero body for x < 0."""
_validate_kernel_par("tau", tau)
M_L = tau * np.exp(-dt_left / tau)
M_R = np.zeros_like(dt_right)
return np.asarray(M_L), M_R
#
def _expRiseCONV_edge_mass(
dt_left: np.ndarray, dt_right: np.ndarray, tau: float
) -> tuple[np.ndarray, np.ndarray]:
"""Exterior masses of the anti-causal kernel: zero body for x > 0."""
_validate_kernel_par("tau", tau)
M_L = np.zeros_like(dt_left)
M_R = tau * np.exp(dt_right / tau)
return M_L, np.asarray(M_R)
#
def _boxCONV_edge_mass(
dt_left: np.ndarray, dt_right: np.ndarray, width: float
) -> tuple[np.ndarray, np.ndarray]:
"""Exterior masses of the unit box on [-width/2, width/2]."""
_validate_kernel_par("width", width)
M_L = np.clip(width / 2 - dt_left, 0.0, width)
M_R = np.clip(dt_right + width / 2, 0.0, width)
return np.asarray(M_L), np.asarray(M_R)
# kernel name -> edge-mass companion; single source of truth consumed by
# both evaluation paths (mcp.Component.convolve and the GIR evaluator)
CONV_EDGE_MASS: dict[str, Callable] = {
"gaussCONV": _gaussCONV_edge_mass,
"expSymCONV": _expSymCONV_edge_mass,
"expDecayCONV": _expDecayCONV_edge_mass,
"expRiseCONV": _expRiseCONV_edge_mass,
"boxCONV": _boxCONV_edge_mass,
}