homodyne.core¶

The homodyne.core package provides the physical foundation for all XPCS computations: constants, parameter bounds, validation logic, and numerically stable JAX-compiled helper functions shared across the NLSQ and CMC backends.

Note

For physics background on the homodyne model and the g1/g2 correlation functions, see the Theory & Physics section.

homodyne.core.physics¶

PhysicsConstants gathers reference values for synchrotron experiments (wavelengths, q-ranges, diffusion coefficient bounds). ValidationResult captures the outcome of parameter validation, including human-readable violation messages.

validate_parameters() checks that a set of physical parameters falls within the bounds specified by ParameterBounds. It returns a ValidationResult and optionally raises ValueError for critical violations.

Physical Constants and Parameter Validation for Homodyne¶

Centralized physical constants, parameter bounds, and validation functions for homodyne scattering analysis. Provides reference values and constraints based on experimental physics and numerical stability requirements.

This module establishes the physical framework for all model computations and ensures parameter values remain within reasonable bounds for stable numerical computation.

class homodyne.core.physics.ValidationResult[source]

Bases: object

Result of parameter validation with detailed error reporting.

Provides comprehensive information about parameter validation including which parameters violated bounds and by how much.

valid

True if all parameters are within bounds

Type:: bool

violations

List of human-readable violation messages

Type:: list of str

parameters_checked

Number of parameters validated

Type:: int

message

Summary message about validation result

Type:: str

valid: bool

violations: list[str]

parameters_checked: int = 0

message: str = ''

__str__()[source]

String representation for logging.

Return type:: str

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

class homodyne.core.physics.PhysicsConstants[source]

Bases: object

Physical constants and reference values for XPCS analysis.

These values are based on typical synchrotron X-ray scattering experiments and provide reasonable defaults for most analyses.

WAVELENGTH_CU_KA = 1.54

WAVELENGTH_8KEV = 1.55

WAVELENGTH_12KEV = 1.03

WAVELENGTH_15KEV = 0.83

Q_MIN_TYPICAL = 0.0001

Q_MAX_TYPICAL = 1.0

TIME_MIN_XPCS = 1e-06

TIME_MAX_XPCS = 1000.0

DIFFUSION_MIN = 100.0

DIFFUSION_MAX = 100000.0

DIFFUSION_TYPICAL = 100.0

SHEAR_RATE_MIN = 1e-06

SHEAR_RATE_MAX = 0.5

SHEAR_RATE_TYPICAL = 0.01

ANGLE_MIN = -10.0

ANGLE_MAX = 10.0

DIFFUSION_OFFSET_MIN = -100000.0

DIFFUSION_OFFSET_MAX = 100000.0

SHEAR_OFFSET_MIN = -0.1

SHEAR_OFFSET_MAX = 0.1

EPS = 1e-12

MAX_EXP_ARG = 700.0

MIN_POSITIVE = 1e-100

ALPHA_MIN = -2.0

ALPHA_MAX = 2.0

BETA_MIN = -2.0

BETA_MAX = 2.0

homodyne.core.physics.parameter_bounds()[source]

Get standard parameter bounds for all model types.

Return type:: dict[str, list[tuple[float, float]]]
Returns:: Dictionary mapping model types to parameter bounds

homodyne.core.physics.validate_parameters_detailed(params, bounds, param_names=None, tolerance=1e-10)[source]

Validate parameter values against bounds with detailed error reporting.

This is the enhanced validation function that provides comprehensive information about which parameters violated bounds and by how much.

Parameters:

params (ndarray) – Parameter array to validate
bounds (list[tuple[float, float]]) – List of (min, max) tuples for each parameter
param_names (list[str] | None) – Names of parameters for better error messages. If None, uses indices.
tolerance (float) – Tolerance for bounds checking (default: 1e-10)

Returns:

Detailed validation result with violations list

Return type:

ValidationResult

Examples

>>> params = np.array([100.0, -1.5, 10.0])
>>> bounds = [(1.0, 1000.0), (-2.0, 2.0), (0.0, 100.0)]
>>> result = validate_parameters_detailed(params, bounds, ["D0", "alpha", "D_offset"])
>>> if not result.valid:
...     print(result.violations)

homodyne.core.physics.validate_parameters(params, bounds, tolerance=1e-10)[source]

Validate parameter values against bounds with tolerance.

This is the legacy function that returns just a boolean. For detailed validation, use validate_parameters_detailed().

Parameters:

params (ndarray) – Parameter array to validate
bounds (list[tuple[float, float]]) – List of (min, max) tuples for each parameter
tolerance (float) – Tolerance for bounds checking

Returns:

True if all parameters are within bounds, False otherwise

Return type:

bool

homodyne.core.physics.clip_parameters(params, bounds)[source]

Clip parameters to stay within bounds.

Parameters:

params (ndarray) – Parameter array to clip
bounds (list[tuple[float, float]]) – List of (min, max) tuples for each parameter

Return type:

ndarray

Returns:

Clipped parameter array

homodyne.core.physics.get_default_parameters(model_type)[source]

Get sensible default parameters for a model type.

Parameters:: model_type (str) – “diffusion”, “shear”, or “combined”
Return type:: ndarray
Returns:: Array of default parameter values

homodyne.core.physics.validate_experimental_setup(q, L, wavelength=None)[source]

Validate experimental setup parameters for physical reasonableness.

Parameters:

q (float) – Scattering wave vector magnitude (Å⁻¹)
L (float) – Sample-detector distance (mm)
wavelength (float | None) – X-ray wavelength (Å), optional

Return type:

bool

Returns:

True if setup is physically reasonable

homodyne.core.physics.estimate_correlation_time(D0, alpha, q)[source]

Estimate characteristic correlation time for diffusion process.

For normal diffusion (alpha=0): τ ≈ 1/(q²D₀) For anomalous diffusion: scaling is more complex

Parameters:

D0 (float) – Reference diffusion coefficient (Å²/s)
alpha (float) – Diffusion exponent
q (float) – Scattering wave vector (Å⁻¹)

Return type:

float

Returns:

Estimated correlation time (seconds)

homodyne.core.physics.get_parameter_info(model_type)[source]

Get comprehensive parameter information for a model type.

Parameters:: model_type (str) – “diffusion”, “shear”, or “combined”
Return type:: dict[str, Any]
Returns:: Dictionary with parameter names, bounds, defaults, and descriptions

homodyne.core.physics_utils¶

Shared utility functions used by both the NLSQ (meshgrid) and CMC (element-wise) computational backends. These were consolidated from jax_backend.py, physics_nlsq.py, and physics_cmc.py to eliminate code duplication and ensure consistent numerical behaviour.

Key functions:

safe_exp — JIT-compiled overflow-protected exponential. Clips input to [-max_val, max_val] before evaluating, preventing NaN propagation.
safe_sinc — Numerically stable unnormalized sinc function with Taylor expansion near zero. Uses EPS = 1e-12 as the zero threshold.
calculate_diffusion_coefficient / _calculate_diffusion_coefficient_impl_jax — Time-dependent diffusion \(D(t) = D_0 \cdot t^\alpha + D_\text{offset}\).
calculate_shear_rate / calculate_shear_rate_cmc — Time-dependent shear rate \(\dot\gamma(t) = \dot\gamma_0 \cdot t^\beta + \dot\gamma_\text{offset}\).

Gradient-safe numerical floors

All positivity floors in physics_utils use jnp.where(x > eps, x, eps) rather than jnp.maximum(x, eps). This preserves non-zero gradients when the input is below the floor — critical for NLSQ Jacobian computation and NUTS leapfrog integration. Using jnp.maximum would silently zero the gradient, causing the optimizer to stall near the boundary.

Warning

safe_exp and safe_sinc are decorated with @jit. Calling them with Python scalars on first use will trigger XLA compilation. Subsequent calls with arrays of the same shape are fast.

Shared Physics Utility Functions for Homodyne¶

This module provides common utility functions and physics helpers used by both NLSQ (meshgrid) and CMC (element-wise) computational backends.

These functions were consolidated from: - jax_backend.py - physics_nlsq.py - physics_cmc.py

to eliminate code duplication and ensure consistent behavior across backends.

Key Functions: - safe_len: JAX-safe length function for scalars and arrays - safe_exp: Overflow-protected exponential - safe_sinc: Numerically stable unnormalized sinc function - _calculate_diffusion_coefficient_impl_jax: Time-dependent diffusion D(t) - _calculate_shear_rate_impl_jax: Time-dependent shear rate γ̇(t) - _create_time_integral_matrix_impl_jax: Trapezoidal cumulative integral matrix

homodyne.core.physics_utils.safe_len(obj)[source]

JAX-safe length function that handles scalars, arrays, and JAX objects.

Parameters:: obj (object) – Any object that might have a length or shape
Returns:: Length of the object, or 1 for scalars
Return type:: int

homodyne.core.physics_utils.safe_exp(x, max_val=700.0)[source]

Safe exponential to prevent overflow.

Parameters:

x (Array) – Input array
max_val (float) – Maximum absolute value to clip to (default 700.0)

Return type:

Array

Returns:

exp(clip(x, -max_val, max_val))

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

Safe UNNORMALIZED sinc function: sin(x) / x (NOT sin(πx) / (πx)).

This matches the reference implementation which uses sin(arg) / arg directly. The phase argument already includes all necessary scaling factors.

P2-4: Uses a Taylor expansion near zero (1 - x²/6 + x⁴/120) for smooth gradient continuity. The old hard switch from sin(x)/x to 1.0 at |x|=EPS created a gradient discontinuity that caused spurious NUTS rejections near gamma_dot_t0 ≈ 0.

Parameters:: x (Array) – Input array
Return type:: Array
Returns:: sin(x)/x for |x| >= 1e-4, Taylor approximation for |x| < 1e-4

homodyne.core.physics_utils.calculate_diffusion_coefficient(time_array, D0, alpha, D_offset)[source]

Calculate time-dependent diffusion coefficient using discrete evaluation.

Follows reference v1 implementation: D_t[i] = D0 * (time_array[i] ** alpha) + D_offset Physical constraint: D(t) should be positive and finite

Parameters:

time_array (Array) – Array of time points
D0 (float) – Diffusion coefficient amplitude
alpha (float) – Anomalous diffusion exponent
D_offset (float) – Baseline diffusion offset

Return type:

Array

Returns:

D(t) evaluated at each time point with physical bounds applied

homodyne.core.physics_utils.calculate_shear_rate(time_array, gamma_dot_0, beta, gamma_dot_offset)[source]

Calculate time-dependent shear rate using discrete evaluation.

Follows reference v1 implementation: γ̇_t[i] = γ̇₀ * (time_array[i] ** β) + γ̇_offset

Parameters:

time_array (Array) – Array of time points
gamma_dot_0 (float) – Shear rate amplitude
beta (float) – Shear rate exponent
gamma_dot_offset (float) – Baseline shear rate offset

Return type:

Array

Returns:

γ̇(t) evaluated at each time point

homodyne.core.physics_utils.calculate_shear_rate_cmc(time_array, gamma_dot_0, beta, gamma_dot_offset)[source]

Calculate time-dependent shear rate for CMC (element-wise) computations.

This variant includes an additional safety check for consecutive zeros in CMC element-wise data where dt could be zero.

Parameters:

time_array (Array) – Array of time points
gamma_dot_0 (float) – Shear rate amplitude
beta (float) – Shear rate exponent
gamma_dot_offset (float) – Baseline shear rate offset

Return type:

Array

Returns:

γ̇(t) evaluated at each time point

homodyne.core.physics_utils.create_time_integral_matrix(time_dependent_array)[source]

Create time integral matrix using trapezoidal numerical integration.

Computes the full N x N matrix of pairwise trapezoidal integral differences via broadcasting. The dt scaling happens in wavevector_q_squared_half_dt, NOT in this cumsum.

Algorithm:

Trapezoidal integration: cumsum[i] = Sum(k=0 to i-1) 0.5 * (f[k] + f[k+1])
Compute full difference matrix: matrix[i,j] = smooth_abs(cumsum[i] - cumsum[j])
The dt factor is applied via wavevector_q_squared_half_dt = 0.5 * q^2 * dt

This gives: matrix[i,j] = number of integration steps. Actual integral: dt * matrix[i,j] approximates the integral from 0 to abs(ti-tj) of f(t’) dt’

Benefits over simple cumsum: - Reduces oscillations from discretization by ~50% - Second-order accuracy (O(dt^2)) vs. first-order (O(dt)) - Eliminates checkerboard artifacts in diagonal-corrected results

Parameters:: time_dependent_array (Array) – f(t) evaluated at discrete time points
Return type:: Array
Returns:: Time integral matrix (in units of integration steps)

homodyne.core.physics_utils.trapezoid_cumsum(values)[source]

Cumulative trapezoid integral without dt scaling (dt is applied outside).

Returns cumsum so that cumsum[j] - cumsum[i] equals the trapezoidal sum over all intervals between indices i and j. The caller applies a smooth absolute value to that difference when mapping each (t1, t2) pair, keeping gradients well-behaved at zero-length intervals.

This is used by the CMC element-wise computations.

Parameters:: values (Array) – 1D array of values to integrate
Return type:: Array
Returns:: Cumulative trapezoidal sums

Usage Example¶

import jax.numpy as jnp
from homodyne.core.physics_utils import safe_exp, safe_sinc

# Safe exponential — no overflow for large inputs
x = jnp.array([0.0, 100.0, 1000.0])
result = safe_exp(x)          # last element clamped, not NaN

# Numerically stable sinc near zero
phi = jnp.linspace(0, 0.1, 10)
s = safe_sinc(phi)            # uses Taylor series for |phi| < EPS

homodyne.core.scaling_utils¶

Quantile-based estimation of per-angle contrast and offset parameters. These utilities implement the physics relation:

\[C_2 = \text{contrast} \times g_1^2 + \text{offset}\]

At large time lags \(g_1^2 \to 0\), so \(C_2 \to \text{offset}\) (the floor). At small time lags \(g_1^2 \approx 1\), so \(C_2 \approx \text{contrast} + \text{offset}\) (the ceiling).

Both NLSQAdapter (anti-degeneracy layer) and fit_mcmc_jax() (per_angle_mode="auto") call estimate_per_angle_scaling() to obtain initialisation values that prevent parameter absorption degeneracy.

Shared scaling utilities for per-angle contrast/offset estimation.

This module provides unified quantile-based estimation of contrast and offset parameters that can be used by both NLSQ and CMC optimization backends.

The physics basis:

C2 = contrast × g1² + offset

At large time lags, g1² → 0, so C2 → offset (the “floor”)
At small time lags, g1² ≈ 1, so C2 ≈ contrast + offset (the “ceiling”)

Version: 2.18.0

homodyne.core.scaling_utils.estimate_contrast_offset_from_quantiles(c2_data, delta_t, contrast_bounds=(0.0, 1.0), offset_bounds=(0.5, 1.5), lag_floor_quantile=0.80, lag_ceiling_quantile=0.20, value_quantile_low=0.10, value_quantile_high=0.90)[source]

Estimate contrast and offset from C2 data using physics-informed quantile analysis.

Uses the correlation decay structure: C2 = contrast × g1² + offset - At large time lags, g1² → 0, so C2 → offset (the “floor”) - At small time lags, g1² ≈ 1, so C2 ≈ contrast + offset (the “ceiling”)

Parameters:

c2_data (ndarray) – C2 correlation values (1D array).
delta_t (ndarray) – Time lag values abs(t1 - t2) (same shape as c2_data).
contrast_bounds (tuple[float, float]) – Valid bounds for contrast parameter.
offset_bounds (tuple[float, float]) – Valid bounds for offset parameter.
lag_floor_quantile (float) – Quantile threshold for “large lag” region (default: 0.80 = top 20% of lags).
lag_ceiling_quantile (float) – Quantile threshold for “small lag” region (default: 0.20 = bottom 20% of lags).
value_quantile_low (float) – Quantile for robust floor estimation (default: 0.10).
value_quantile_high (float) – Quantile for robust ceiling estimation (default: 0.90).

Returns:

(contrast_est, offset_est) - Estimated values clipped to bounds.

Return type:

tuple[float, float]

Notes

The estimation is robust to outliers by using quantiles instead of min/max. The lag-based segmentation ensures we’re sampling from the appropriate regions of the correlation decay curve.

homodyne.core.scaling_utils.estimate_per_angle_scaling(c2_data, t1, t2, phi_indices, n_phi, contrast_bounds, offset_bounds, log=None)[source]

Estimate contrast and offset initial values for each phi angle.

This is the unified implementation used by both NLSQ and CMC backends.

Optimization (v2.9.1): Uses vectorized grouped operations instead of sequential loop over angles. Provides 3-5x speedup for typical datasets with 20+ phi angles.

Parameters:

c2_data (ndarray) – Pooled C2 correlation values.
t1 (ndarray) – Pooled first time coordinates.
t2 (ndarray) – Pooled second time coordinates.
phi_indices (ndarray) – Index mapping each data point to its phi angle (0 to n_phi-1).
n_phi (int) – Number of unique phi angles.
contrast_bounds (tuple[float, float]) – Valid bounds for contrast.
offset_bounds (tuple[float, float]) – Valid bounds for offset.
log (Logger | LoggerAdapter[Logger] | None) – Logger for diagnostic messages.

Returns:

Dictionary with keys ‘contrast_0’, ‘offset_0’, ‘contrast_1’, ‘offset_1’, etc.

Return type:

dict[str, float]

homodyne.core.scaling_utils.compute_averaged_scaling(c2_data, t1, t2, phi_indices, n_phi, contrast_bounds, offset_bounds, log=None)[source]

Compute averaged contrast and offset for constant mode.

This function estimates per-angle contrast/offset using quantile analysis, then averages them to produce single values for constant mode optimization.

Parameters:

c2_data (ndarray) – Pooled C2 correlation values.
t1 (ndarray) – Pooled first time coordinates.
t2 (ndarray) – Pooled second time coordinates.
phi_indices (ndarray) – Index mapping each data point to its phi angle (0 to n_phi-1).
n_phi (int) – Number of unique phi angles.
contrast_bounds (tuple[float, float]) – Valid bounds for contrast.
offset_bounds (tuple[float, float]) – Valid bounds for offset.
log (Logger | LoggerAdapter[Logger] | None) – Logger for diagnostic messages.

Returns:

(contrast_avg, offset_avg, contrast_per_angle, offset_per_angle) - contrast_avg: Averaged contrast for constant mode - offset_avg: Averaged offset for constant mode - contrast_per_angle: Per-angle estimates (for diagnostics) - offset_per_angle: Per-angle estimates (for diagnostics)

Return type:

tuple[float, float, ndarray, ndarray]

Usage Example¶

import numpy as np
from homodyne.core.scaling_utils import estimate_per_angle_scaling

# c2_data: shape (N,), delta_t: shape (N,)
c2_data = np.random.uniform(1.0, 1.5, size=5000)
delta_t = np.abs(np.random.randn(5000))

contrast, offset = estimate_per_angle_scaling(c2_data, delta_t)
print(f"contrast={contrast:.3f}, offset={offset:.3f}")