Core Physics¶

Foundational physics constants, safe numerical primitives, and integral building blocks shared by both the NLSQ meshgrid path and the CMC element-wise path.

Physics Constants and Validation¶

Physical constants, parameter bounds, and validation for heterodyne model.

class heterodyne.core.physics.ValidationResult[source]

Bases: object

Result of parameter validation with detailed error reporting.

valid: True if all parameters are within bounds.

violations: List of human-readable violation messages.

parameters_checked: Number of parameters validated.

message: Summary message about validation result.

valid: bool

violations: list[str]

parameters_checked: int = 0

message: str = ''

__init__(valid, violations=<factory>, parameters_checked=0, message='')

class heterodyne.core.physics.PhysicsConstants[source]

Bases: object

Physical constants for XPCS scattering analysis.

All values in SI base units unless otherwise noted.

k_B: ClassVar[float] = 1.380649e-23

h: ClassVar[float] = 6.62607015e-34

c: ClassVar[float] = 299792458.0

WAVELENGTH_8KEV: ClassVar[float] = 1.55

WAVELENGTH_10KEV: ClassVar[float] = 1.24

WAVELENGTH_12KEV: ClassVar[float] = 1.0332

Q_MIN_TYPICAL: ClassVar[float] = 0.0001

Q_MAX_TYPICAL: ClassVar[float] = 1.0

TIME_MIN_XPCS: ClassVar[float] = 1e-06

TIME_MAX_XPCS: ClassVar[float] = 1000.0

VELOCITY_MIN: ClassVar[float] = 1e-06

VELOCITY_MAX: ClassVar[float] = 10000.0

FRACTION_MIN: ClassVar[float] = 0.0

FRACTION_MAX: ClassVar[float] = 1.0

EPS: ClassVar[float] = 1e-12

MAX_EXP_ARG: ClassVar[float] = 700.0

MIN_POSITIVE: ClassVar[float] = 1e-100

__init__()

heterodyne.core.physics.get_default_bounds_array()[source]

Get default bounds as arrays in canonical parameter order.

Return type:: tuple[ndarray, ndarray]
Returns:: (lower_bounds, upper_bounds) each of shape (14,)

class heterodyne.core.physics.TransportPhysics[source]

Bases: object

Physical interpretation of transport parameters.

Transport coefficient: J(t) = D0 * t^alpha + offset

Physical regimes based on alpha: - alpha = 1.0: Normal (Brownian) diffusion - alpha < 1.0: Subdiffusion (crowded/constrained) - alpha > 1.0: Superdiffusion (active/directed) - alpha = 2.0: Ballistic motion

NORMAL_DIFFUSION: ClassVar[float] = 1.0

BALLISTIC: ClassVar[float] = 2.0

static interpret_alpha(alpha)[source]

Interpret alpha value physically.

For J(t) = D0 * t^alpha + offset: - alpha ≈ 0: constant transport rate (equilibrium) - alpha ≈ 1: linearly growing transport (normal diffusion) - alpha < 1: sub-linear (subdiffusive) - alpha > 1: super-linear (superdiffusive) - alpha ≈ 2: quadratic (ballistic)

Parameters:: alpha (float) – Transport rate exponent
Return type:: str
Returns:: Physical interpretation string

static diffusion_coefficient(D0, alpha, t=1.0)[source]

Compute effective diffusion coefficient at time t.

For J(t) = D0 * t^alpha, the effective D is: D_eff = dJ/dt = D0 * alpha * t^(alpha-1)

Parameters:

D0 (float) – Transport prefactor
alpha (float) – Transport exponent
t (float) – Time point (default 1.0)

Return type:

float

Returns:

Effective diffusion coefficient

__init__()

heterodyne.core.physics.validate_parameters(params, bounds=None)[source]

Validate parameter values against physical bounds.

Parameters:

params (dict[str, float]) – Dictionary mapping parameter names to values.
bounds (dict[str, tuple[float, float]] | None) – Optional custom bounds. Defaults to PARAMETER_BOUNDS.

Return type:

ValidationResult

Returns:

ValidationResult with violations list if any bounds are exceeded.

Numerical Primitives¶

Safe arithmetic wrappers and integral utilities used throughout the correlation computation pipeline.

Numerically safe mathematical primitives for heterodyne physics.

All functions are designed to be compatible with both NumPy and JAX arrays, avoiding NaN/Inf from edge cases (division by zero, overflow in exp, negative bases in power).

Shared utilities used by both the NLSQ meshgrid path and the CMC element-wise path: - trapezoid_cumsum: O(dt²) cumulative integral - create_time_integral_matrix: N×N from cumsum (NLSQ only) - smooth_abs: gradient-safe |x| for NUTS - compute_transport_rate: J(t) = D0·t^α + offset - compute_velocity_rate: v(t) = v0·t^β + v_offset

heterodyne.core.physics_utils.safe_exp(x, limit=500.0)[source]

Exponential with overflow protection.

Clips the argument to [-limit, limit] before computing exp() to avoid Inf outputs. The default limit of 500 gives exp(500) ≈ 1.4e217 which is within float64 range.

Parameters:

x (Array | ndarray) – Input array
limit (float) – Clipping threshold (symmetric)

Return type:

Array

Returns:

exp(clip(x)), same shape as x

heterodyne.core.physics_utils.safe_power(base, exponent)[source]

Power function safe for non-positive bases.

For base ≤ 0, returns 0.0 (the physical limit for t^α transport). For base > 0, returns base^exponent normally.

Parameters:

base (Array | ndarray) – Base array (typically time values)
exponent (float) – Power exponent

Return type:

Array

Returns:

Safe power result, same shape as base

heterodyne.core.physics_utils.safe_divide(numerator, denominator, fill=0.0, min_denom=1e-30)[source]

Division with protection against zero/near-zero denominators.

Parameters:

numerator (Array | ndarray) – Dividend array
denominator (Array | ndarray) – Divisor array
fill (float) – Value to return where denominator is too small
min_denom (float) – Minimum absolute denominator value

Return type:

Array

Returns:

Safe quotient, same shape as inputs

heterodyne.core.physics_utils.safe_log(x, floor=1e-30)[source]

Logarithm with protection against non-positive arguments.

Parameters:

x (Array | ndarray) – Input array
floor (float) – Minimum value before taking log

Return type:

Array

Returns:

log(max(x, floor)), same shape as x

heterodyne.core.physics_utils.safe_sqrt(x)[source]

Square root with protection against negative arguments.

Parameters:: x (Array | ndarray) – Input array
Return type:: Array
Returns:: sqrt(max(x, 0)), same shape as x

heterodyne.core.physics_utils.compute_relative_difference(a, b)[source]

Compute element-wise relative difference |a - b| / max(|a|, |b|, 1e-10).

Useful for comparing correlation matrices or parameter arrays where absolute differences may mislead at different scales.

Parameters:

a (Array | ndarray) – First array
b (Array | ndarray) – Second array

Return type:

Array

Returns:

Relative difference array, values in [0, 2]

heterodyne.core.physics_utils.symmetrize(matrix)[source]

Force a matrix to be exactly symmetric: (M + M^T) / 2.

Parameters:: matrix (Array | ndarray) – Square matrix
Return type:: Array
Returns:: Symmetric matrix

heterodyne.core.physics_utils.smooth_abs(x, eps=1e-12)[source]

Gradient-safe absolute value: sqrt(x² + ε).

jnp.abs(x) has undefined gradient at x=0, which causes NaN in NUTS backpropagation on matrix diagonals where integrals are zero. This smooth approximation matches |x| to O(√ε) and has well-defined gradients everywhere.

Parameters:

x (Array) – Input array.
eps (float) – Smoothing parameter. 1e-12 gives ~1e-6 bias on diagonal.

Return type:

Array

Returns:

Smooth |x|, same shape as x.

heterodyne.core.physics_utils.smooth_clip(x, low, high, sharpness=50.0)[source]

Soft clip to [low, high] with continuous gradient at the boundaries.

Acts as the identity in the interior and softplus-smoothed at the boundaries. Use this for physical bounds (e.g. sample fraction in [0, 1]) where a hard jnp.clip would zero the gradient and stall NUTS leapfrog integration or NLSQ Jacobian descent (CLAUDE.md rule #7).

The boundary smear scales as 1/sharpness — at the default value sharpness=50 the boundary lands within ~``(high-low)/50`` of the target (≈2% of the range), with monotonic identity in the interior. Raise sharpness for a tighter approximation at the cost of gradient magnitude near the boundary; lower it for stronger regularisation.

Parameters:

x (Array) – Input array (any shape).
low (float) – Lower physical bound (inclusive in the limit).
high (float) – Upper physical bound (inclusive in the limit).
sharpness (float) – Softplus sharpness; default 50 gives ~2% boundary smear.

Return type:

Array

Returns:

Smoothly bounded array, asymptotically in (low, high), with well-defined gradients everywhere.

heterodyne.core.physics_utils.trapezoid_cumsum(f, dt)[source]

Trapezoidal cumulative integral with O(dt²) accuracy.

Computes cumsum[0] = 0, cumsum[k] = Σ_{i=0}^{k-1} (f[i]+f[i+1])/2 × dt.

This matches homodyne’s trapezoid_cumsum pattern. The dt factor is included in the returned values (unlike homodyne which factors it out into the wavevector prefactor).

Parameters:

f (Array) – Function values at uniformly spaced time points, shape (N,).
dt (float | Array) – Time step.

Return type:

Array

Returns:

Cumulative integral, shape (N,). cumsum[0] = 0 always.

heterodyne.core.physics_utils.create_time_integral_matrix(cumsum_values)[source]

Build N×N integral matrix from cumulative sums (NLSQ meshgrid path).

M[i,j] = cumsum[j] - cumsum[i] (signed difference).

For transport integrals, call smooth_abs on the result to get direction-independent decay. For velocity integrals, use the signed result directly (it feeds into cos(q cos(φ) ∫v dt)).

Parameters:: cumsum_values (Array) – Cumulative integral, shape (N,).
Return type:: Array
Returns:: Signed integral matrix, shape (N, N).

heterodyne.core.physics_utils.compute_transport_rate(t, D0, alpha, offset)[source]

Transport rate function J(t) = D0·t^α + offset.

Shared by both NLSQ and CMC paths — the rate function is the same, only the integral evaluation strategy differs.

Parameters:

t (Array) – Time array, shape (N,).
D0 (float | Array) – Transport prefactor (Å²/s^α).
alpha (float | Array) – Transport exponent (dimensionless).
offset (float | Array) – Constant rate offset (Å²/s).

Return type:

Array

Returns:

Rate values, shape (N,), floored at 0.

heterodyne.core.physics_utils.compute_velocity_rate(t, v0, beta, v_offset)[source]

Velocity rate function v(t) = v0·t^β + v_offset.

Unlike transport rate, the velocity is NOT floored at 0 because the velocity integral enters as cos(q·cos(φ)·∫v dt) which is naturally bounded.

Parameters:

t (Array) – Time array, shape (N,).
v0 (float | Array) – Velocity prefactor (Å/s^β).
beta (float | Array) – Velocity exponent (dimensionless).
v_offset (float | Array) – Constant velocity offset (Å/s).

Return type:

Array

Returns:

Velocity values, shape (N,).

heterodyne.core.physics_utils.safe_sinc(x)[source]

Unnormalized sinc function sin(x)/x, safe at x=0.

Returns 1.0 at x=0 (the mathematical limit).

Parameters:: x (Array) – Input array (radians, unnormalized).
Return type:: Array
Returns:: sin(x)/x with sinc(0) = 1.

Physics Factors¶

Pre-computed physics factor tables for correlation model evaluation.

Pre-computed physics factors for efficient correlation computation.

class heterodyne.core.physics_factors.PhysicsFactors[source]

Bases: object

Pre-computed factors that don’t depend on fit parameters.

These are computed once from experimental setup and reused across all optimization iterations for efficiency.

t: Array

q: float

q_squared: float

dt: float

n_times: int

phi_angle: float

__post_init__()[source]

Validate factors.

Return type:: None

property time_extent: float: Total time span.

get_q_cosine(phi0=0.0)[source]

Get q * cos(phi_total) for cross-term phase.

Parameters:: phi0 (float) – Additional angle from fit parameters
Return type:: Array
Returns:: q * cos(phi_angle + phi0) as JAX scalar

__init__(t, q, q_squared, dt, n_times, phi_angle)

heterodyne.core.physics_factors.create_physics_factors(n_times, dt, q, phi_angle=0.0, t_start=0.0)[source]

Create physics factors from experimental parameters.

Parameters:

n_times (int) – Number of time points
dt (float) – Time step
q (float) – Scattering wavevector magnitude
phi_angle (float) – Detector phi angle (degrees)
t_start (float) – Starting time (default 0)

Return type:

PhysicsFactors

Returns:

PhysicsFactors instance

heterodyne.core.physics_factors.create_physics_factors_from_config(config)[source]

Create physics factors from configuration dictionary.

Reads from analyzer_parameters (canonical) with fallback to legacy temporal/scattering top-level sections for backwards compatibility.

Parameters:: config (dict) – Configuration with analyzer_parameters or legacy temporal/scattering sections.
Return type:: PhysicsFactors
Returns:: PhysicsFactors instance

class heterodyne.core.physics_factors.CachedMatrices[source]

Bases: object

Cached matrices that depend only on time grid.

These are expensive to recompute and don’t change during fitting.

time_diff: Array

mean_time: Array

triu_indices: tuple[Array, Array]

tril_indices: tuple[Array, Array]

__init__(time_diff, mean_time, triu_indices, tril_indices)

heterodyne.core.physics_factors.create_cached_matrices(factors)[source]

Create cached matrices from physics factors.

Parameters:: factors (PhysicsFactors) – PhysicsFactors instance
Return type:: CachedMatrices
Returns:: CachedMatrices instance