SPP Analysis API

The Sequence of Physical Processes (SPP) tooling in RheoJAX provides time-domain LAOS analysis without relying on Fourier/Chebyshev expansions. It is built from three primary objects:

• rheojax.transforms.spp_decomposer.SPPDecomposer — per-cycle LAOS decomposition into cage modulus, static/dynamic yield stresses, power-law flow, and nonlinearity metrics.

• rheojax.models.spp.spp_yield_stress.SPPYieldStress — a fit-ready model (NLSQ and NumPyro NUTS) that parameterizes SPP yield behavior across amplitudes or steady shear.

• rheojax.pipeline.workflows.SPPAmplitudeSweepPipeline — convenience pipeline for amplitude-sweep LAOS workflows.

For conceptual background see Sequence of Physical Processes (SPP).

SPPDecomposer

class rheojax.transforms.spp_decomposer.SPPDecomposer(omega, gamma_0, n_harmonics=39, yield_tolerance=0.02, start_cycle=0, end_cycle=None, use_numerical_method=False, step_size=8, num_mode=2, wrap_strain_rate=True)

Bases: BaseTransform

SPP decomposition transform for LAOS stress analysis.

Applies the Sequence of Physical Processes (SPP) framework to decompose LAOS stress signals and extract nonlinear viscoelastic parameters.

The transform requires oscillatory shear data with known frequency and strain amplitude. It computes:

1. Elastic/viscous stress decomposition

2. Yield stress extraction (static and dynamic)

3. Power-law flow parameters

4. Lissajous-Bowditch metrics

5. Harmonic decomposition

Parameters:

• omega (float) – Angular frequency ω (rad/s) of the oscillation

• gamma_0 (float) – Strain amplitude γ_0 (dimensionless)

• n_harmonics (int) – Number of odd harmonics to extract for stress (default: 39 per MATLAB SPPplus)

• yield_tolerance (float) – Fractional tolerance for yield point detection (default: 0.02)

• start_cycle (int) – First cycle to analyze (0-indexed). Use to skip startup transients. Default: 0 (start from beginning).

• end_cycle (int | None) – Last cycle to analyze (exclusive). None means use all available cycles. Default: None.

• use_numerical_method (bool) – If True, use MATLAB-compatible numerical differentiation for raw data. If False (default), use Fourier-based decomposition.

• step_size (int) – Step size k for numerical differentiation (only used if use_numerical_method=True). Default: 8 to mirror SPPplus v2.1.

• num_mode (int) – Numerical differentiation mode (1 = edge-aware, 2 = periodic/looped), matching SPPplus num_mode. Only used when use_numerical_method=True.

omega

Angular frequency

Type:

float

gamma_0

Strain amplitude

Type:

float

gamma_dot_0

Strain rate amplitude (ω * γ_0)

Type:

float

n_harmonics

Number of harmonics for decomposition

Type:

int

start_cycle

First cycle to analyze

Type:

int

end_cycle

Last cycle to analyze

Type:

int or None

use_numerical_method

Whether using numerical differentiation

Type:

bool

results_

Dictionary of computed SPP metrics (after transform)

Type:

dict

Examples

Basic usage with RheoData:

>>> from rheojax.core.data import RheoData
>>> from rheojax.transforms.spp_decomposer import SPPDecomposer
>>>
>>> # LAOS stress-strain data
>>> omega = 1.0  # rad/s
>>> gamma_0 = 1.0  # strain amplitude
>>> t = jnp.linspace(0, 2*jnp.pi, 1000)
>>> strain = gamma_0 * jnp.sin(omega * t)
>>> stress = 100.0 * strain + 20.0 * jnp.sin(3 * omega * t)  # With 3rd harmonic
>>>
>>> data = RheoData(
...     x=t,
...     y=stress,
...     domain='time',
...     metadata={
...         'test_mode': 'oscillation',
...         'omega': omega,
...         'gamma_0': gamma_0,
...         'strain': strain,
...     }
... )
>>>
>>> # Apply SPP decomposition
>>> decomposer = SPPDecomposer(omega=omega, gamma_0=gamma_0)
>>> result = decomposer.transform(data)
>>>
>>> # Access metrics
>>> print(f"Static yield stress: {decomposer.results_['sigma_sy']:.2f} Pa")
>>> print(f"Dynamic yield stress: {decomposer.results_['sigma_dy']:.2f} Pa")

Notes

• Input data must be time-domain stress signal

• Strain data must be available in metadata[‘strain’] or computed from ω, γ_0

• Output includes both decomposed waveforms and extracted parameters

__init__(omega, gamma_0, n_harmonics=39, yield_tolerance=0.02, start_cycle=0, end_cycle=None, use_numerical_method=False, step_size=8, num_mode=2, wrap_strain_rate=True)

Initialize SPP decomposer transform.

Parameters:

• omega (float) – Angular frequency (rad/s)

• gamma_0 (float) – Strain amplitude (dimensionless)

• n_harmonics (int) – Number of odd harmonics to extract (default: 39)

• yield_tolerance (float) – Tolerance for yield point detection (default: 0.02)

• start_cycle (int) – First cycle to analyze, 0-indexed (default: 0)

• end_cycle (int | None) – Last cycle to analyze, exclusive (default: None, use all)

• use_numerical_method (bool) – Use MATLAB-compatible numerical differentiation (default: False)

• step_size (int) – Step size k for numerical differentiation (default: 8)

• num_mode (int) – Numerical differentiation mode (1=edge-aware, 2=periodic). Default: 2.

• wrap_strain_rate (bool) – If True, infer strain rate with periodic wrapping when missing (default: True)

get_results()

Get computed SPP analysis results.

Returns:

Dictionary containing all SPP metrics: - sigma_sy: Static yield stress (Pa) - sigma_dy: Dynamic yield stress (Pa) - K: Consistency index (Pa·s^n) - n_power_law: Power-law exponent - harmonic_amplitudes: Array of harmonic amplitudes - harmonic_phases: Array of harmonic phases - I3_I1_ratio: Third harmonic nonlinearity ratio - G_L, G_M: Large and minimum strain moduli (Pa) - eta_L, eta_M: Large and minimum rate viscosities (Pa·s) - S_factor: Stiffening ratio - T_factor: Thickening ratio - G_cage: Time-resolved cage modulus (array) - sigma_elastic: Elastic stress contribution (array) - sigma_viscous: Viscous stress contribution (array)

Return type:

dict

Raises:

RuntimeError – If transform has not been applied yet

Examples

>>> decomposer = SPPDecomposer(omega=1.0, gamma_0=1.0)
>>> _ = decomposer.transform(data)
>>> results = decomposer.get_results()
>>> print(f"I3/I1 = {results['I3_I1_ratio']:.4f}")

get_yield_stresses()

Get static and dynamic yield stresses.

Returns:

(sigma_sy, sigma_dy) in Pa

Return type:

tuple[float, float]

Raises:

RuntimeError – If transform has not been applied yet

get_nonlinearity_metrics()

Get nonlinearity quantification metrics.

Returns:

Dictionary with: - I3_I1_ratio: Third harmonic ratio (FT-rheology) - S_factor: Strain stiffening ratio - T_factor: Shear thickening ratio

Return type:

dict

Raises:

RuntimeError – If transform has not been applied yet

__repr__()

String representation of transform.

Return type:

str

SPPYieldStress Model

class rheojax.models.spp.spp_yield_stress.SPPYieldStress

Bases: BaseModel

SPP-based yield stress model for LAOS analysis.

This model extracts yield stress parameters from amplitude-sweep LAOS data using the SPP framework. It parameterizes the nonlinear response in terms of physically meaningful quantities: cage modulus, static/dynamic yield stresses, and post-yield flow parameters.

The model supports both OSCILLATION mode (LAOS amplitude sweeps) and ROTATION mode (steady shear flow curves) for comprehensive yield stress characterization.

Parameters:

• G_cage (float) – Apparent cage modulus (Pa), elastic stiffness of the cage structure

• sigma_sy_scale (float) – Static yield stress scale factor (Pa)

• sigma_sy_exp (float) – Static yield stress exponent for amplitude scaling

• sigma_dy_scale (float) – Dynamic yield stress scale factor (Pa)

• sigma_dy_exp (float) – Dynamic yield stress exponent for amplitude scaling

• eta_inf (float) – Infinite-shear viscosity (Pa·s)

• n_power_law (float) – Power-law flow index (dimensionless)

• noise (float) – Observation noise scale (Pa), for Bayesian inference

• Equations (Constitutive)

• ----------------------

• γ_0) (For OSCILLATION (LAOS at amplitude) – σ_sy(γ_0) = sigma_sy_scale * γ_0^sigma_sy_exp σ_dy(γ_0) = sigma_dy_scale * γ_0^sigma_dy_exp σ_max(γ_0) = G_cage * γ_0 + η_inf * ω * γ_0

• γ̇) (For ROTATION (steady shear at rate) – σ(γ̇) = σ_dy + η_inf * γ̇^n_power_law (Herschel-Bulkley-like)

• Modes (Test)

• ----------

• OSCILLATION (-)

• ROTATION (-)

Examples

Fit to amplitude sweep data:

>>> from rheojax.models import SPPYieldStress
>>> from rheojax.transforms import SPPDecomposer
>>>
>>> # Amplitude sweep data (multiple γ_0 values)
>>> gamma_0_values = jnp.array([0.01, 0.1, 0.5, 1.0, 2.0, 5.0])
>>> sigma_sy_data = jnp.array([1.0, 10.0, 40.0, 80.0, 140.0, 280.0])
>>>
>>> model = SPPYieldStress()
>>> model.fit(gamma_0_values, sigma_sy_data, test_mode='oscillation')
>>> print(model)

Bayesian inference:

>>> result = model.fit_bayesian(gamma_0_values, sigma_sy_data,
...                             test_mode='oscillation')
>>> print(result.summary)

__init__()

Initialize SPP yield stress model.

model_function(X, params, test_mode=None, **kwargs)

Model function for predictions and Bayesian inference.

This is the NumPyro-facing model used by BayesianMixin. It maps inputs and parameter vectors to stress predictions.

Parameters:

• X (Array) – Independent variable - OSCILLATION: strain amplitudes (gamma_0) - ROTATION: shear rates (gamma_dot)

• params (Array) –

  [G_cage, sigma_sy_scale, sigma_sy_exp, sigma_dy_scale,

  sigma_dy_exp, eta_inf, n_power_law, noise]

• test_mode (TestModeEnum | None) – Mode selector; defaults to the model’s stored test_mode.

Returns:

Predicted stress values.

Return type:

Array

bayesian_prior_factory(name, lower, upper)

Create NumPyro-friendly priors for model parameters.

Uses physically-motivated prior distributions: - LogNormal for scale parameters (positive, often log-distributed) - Beta for bounded exponents - HalfCauchy for noise scale

Parameters:

• name (str) – Parameter name

• lower (float | None) – Lower bound

• upper (float | None) – Upper bound

Returns:

NumPyro distribution, or None to use default uniform

Return type:

Distribution | None

predict_amplitude_sweep(gamma_0_array, omega=1.0, yield_type='static')

Predict full amplitude sweep response.

Computes both static and dynamic yield stresses across amplitudes.

Parameters:

• gamma_0_array (ndarray) – Strain amplitudes to evaluate

• omega (float) – Angular frequency (rad/s)

• yield_type (str) – ‘static’, ‘dynamic’, or ‘both’

Returns:

Dictionary with: - gamma_0: strain amplitudes - sigma_sy: static yield stresses (if requested) - sigma_dy: dynamic yield stresses (if requested)

Return type:

dict

predict_flow_curve(gamma_dot_array)

Predict steady shear flow curve.

Parameters:

gamma_dot_array (ndarray) – Shear rates to evaluate

Returns:

Dictionary with: - gamma_dot: shear rates - sigma: shear stress - eta_app: apparent viscosity

Return type:

dict

predict_rheo(rheo_data, test_mode=None)

Predict rheological response for RheoData.

Parameters:

• rheo_data (RheoData) – Input rheological data

• test_mode (TestModeEnum | None) – Test mode override

Returns:

Predicted response

Return type:

RheoData

__repr__()

String representation.

Return type:

str

SPP Amplitude Sweep Pipeline

class rheojax.pipeline.workflows.SPPAmplitudeSweepPipeline(omega=1.0, n_harmonics=39, step_size=8, num_mode=2, wrap_strain_rate=True, use_numerical_method=None)

Bases: Pipeline

Pipeline for SPP analysis of amplitude sweep LAOS data.

This pipeline performs SPP (Sequence of Physical Processes) analysis on amplitude sweep LAOS data to extract yield stress parameters and nonlinear viscoelastic metrics.

Workflow:

1. Load amplitude sweep data (multiple γ_0 values)

2. Apply SPP decomposition at each amplitude

3. Extract yield stresses (static and dynamic)

4. Fit power-law scaling to yield stress vs amplitude

5. Optionally fit Bayesian SPPYieldStress model

omega

Angular frequency of oscillation (rad/s)

results

Dictionary of SPP metrics per amplitude

model

Fitted SPPYieldStress model (after fit_model)

Example

>>> pipeline = SPPAmplitudeSweepPipeline(omega=1.0)
>>> pipeline.run(amplitude_data_list)
>>> pipeline.fit_model(bayesian=True)
>>> print(pipeline.get_yield_stresses())

__init__(omega=1.0, n_harmonics=39, step_size=8, num_mode=2, wrap_strain_rate=True, use_numerical_method=None)

Initialize SPP amplitude sweep pipeline.

Parameters:

• omega (float) – Angular frequency in rad/s (default: 1.0)

• n_harmonics (int) – Number of harmonics for SPP decomposition (default: 39)

• step_size (int) – Differentiation step size k (default: 8, Rogers parity)

• num_mode (int) – Numerical differentiation mode (default: 2 periodic)

• wrap_strain_rate (bool) – Whether to use wrapped differentiation when rate missing

• use_numerical_method (bool | None) – Force numerical path; None keeps default from transform

run(stress_data, gamma_0_values=None)

Execute SPP analysis on amplitude sweep data.

Parameters:

• stress_data (list[RheoData]) – List of RheoData objects, one per amplitude

• gamma_0_values (list[float] | None) – Strain amplitudes corresponding to each dataset. If None, extracted from RheoData metadata.

Return type:

SPPAmplitudeSweepPipeline

Returns:

self for method chaining

Raises:

ValueError – If gamma_0_values not provided and not in metadata

fit_model(bayesian=False, yield_type='static', **fit_kwargs)

Fit SPPYieldStress model to extracted yield stresses.

Parameters:

• bayesian (bool) – Whether to use Bayesian inference (default: False)

• yield_type (str) – Which yield stress to fit (‘static’ or ‘dynamic’)

• **fit_kwargs – Additional arguments passed to fit or fit_bayesian

Return type:

SPPAmplitudeSweepPipeline

Returns:

self for method chaining

get_yield_stresses()

Get extracted yield stresses from amplitude sweep.

Returns:

• gamma_0: strain amplitudes

• sigma_sy: static yield stresses

• sigma_dy: dynamic yield stresses

Return type:

dict[str, ndarray]

get_amplitude_results(gamma_0)

Get full SPP results for a specific amplitude.

Parameters:

gamma_0 (float) – Strain amplitude to retrieve

Return type:

dict

Returns:

Dictionary of SPP metrics for that amplitude

Raises:

KeyError – If amplitude not in results

get_model()

Get fitted SPPYieldStress model.

Return type:

Any

Returns:

Fitted model or None if not fitted

get_nonlinearity_metrics()

Get nonlinearity metrics (I3/I1, S, T) for each amplitude.

Return type:

dict[float, dict]

Returns:

Dictionary mapping gamma_0 to nonlinearity metrics

Helper Functions

These JAX-first kernels are re-exported for completeness; see rheojax.utils.spp_kernels for details.

JAX-compatible SPP (Sequence of Physical Processes) kernel functions.

This module provides efficient, JAX-compatible implementations of SPP analysis kernel functions for LAOS (Large Amplitude Oscillatory Shear) rheology. SPP analysis enables cycle-by-cycle decomposition of nonlinear stress responses into elastic and viscous contributions, extracting physically meaningful yield parameters.

The SPP framework was developed by Rogers (2012, 2017) and provides: - Time-resolved apparent cage modulus G’_cage(t) - Static and dynamic yield stress extraction - Phase reconstruction from harmonic decomposition - Lissajous-Bowditch plot metrics - Frenet-Serret frame analysis (T, N, B vectors) - Moduli rate calculations (Ġ’, Ġ’’, G_speed)

Key Functions

• apparent_cage_modulus: Time-resolved elastic modulus from stress/strain

• static_yield_stress: Yield stress at strain reversal (strain = ±gamma0)

• dynamic_yield_stress: Yield stress at rate reversal (strain rate = 0)

• harmonic_reconstruction: Stress reconstruction from Fourier components

• harmonic_reconstruction_full: Full Fourier with phase alignment (MATLAB-compatible)

• spp_fourier_analysis: Complete Fourier-based SPP analysis with analytical derivatives

• lissajous_metrics: Bowditch diagram derived quantities (S, T ratios)

• zero_crossing_detection: Robust strain/rate zero-crossing finder

• frenet_serret_frame: Compute T, N, B trajectory vectors

• moduli_rates: Compute Ġ’(t), Ġ’’(t), G_speed, δ̇(t)

• yield_from_displacement_stress: SPP-based yield stress extraction

Physical Interpretation

SPP analysis provides a phenomenological interpretation of LAOS behavior: - G’_cage(t): Instantaneous elastic modulus reflecting cage structure - static_yield: Static yield stress (stress at max strain, cage breakage threshold) - dynamic_yield: Dynamic yield stress (stress at zero rate, flow cessation threshold) - Power-law regime: Post-yield flow characterized by sigma ~ strain_rate^n - Frenet-Serret frame: Geometric analysis of (γ, γ̇/ω, σ) trajectory

References

• S.A. Rogers et al., “A sequence of physical processes determined and quantified in large-amplitude oscillatory shear (LAOS): Application to theoretical nonlinear models”, J. Rheol. 56(1), 2012

• S.A. Rogers, “In search of physical meaning: defining transient parameters for nonlinear viscoelasticity”, Rheol. Acta 56, 2017

• G.J. Donley et al., “Time-resolved dynamics of the yielding transition in soft materials”, J. Non-Newton. Fluid Mech. 264, 2019

rheojax.utils.spp_kernels.apparent_cage_modulus(stress, strain, strain_amplitude)

Compute time-resolved apparent cage modulus.

Apparent cage modulus is the instantaneous elastic response, normalized by strain amplitude: G_cage(t) = stress(t) / gamma0 * sign(strain(t)).

Reference: Rogers et al. (2012) J. Rheol. 56(1). Eq. (1): G’_cage(t) = sigma(t) / gamma_0 * sign(gamma(t))

Parameters:

• stress (Array) – Time-resolved stress signal (Pa)

• strain (Array) – Time-resolved strain signal (dimensionless)

• strain_amplitude (float) – Maximum strain amplitude gamma0 (dimensionless)

Returns:

Apparent cage modulus (Pa)

Return type:

Array

Examples

>>> import jax.numpy as jnp
>>> from rheojax.utils.spp_kernels import apparent_cage_modulus
>>>
>>> # Sinusoidal LAOS data
>>> t = jnp.linspace(0, 2*jnp.pi, 1000)
>>> gamma_0 = 1.0
>>> gamma = gamma_0 * jnp.sin(t)
>>> sigma = 100.0 * jnp.sin(t) + 10.0 * jnp.sin(3*t)  # With 3rd harmonic
>>>
>>> G_cage = apparent_cage_modulus(sigma, gamma, gamma_0)

Notes

• G_cage is constant for a purely linear sinusoidal material

• Deviations from constant indicate nonlinearity

• Sign(γ) ensures correct sign during negative strain half-cycle

rheojax.utils.spp_kernels.static_yield_stress(stress, strain, strain_amplitude, tolerance=0.02)

Approximate static yield stress near strain reversal.

Samples near the strain extrema (abs(strain) close to strain_amplitude) and returns the average absolute stress.

Return type:

float

rheojax.utils.spp_kernels.dynamic_yield_stress(stress, strain_rate, rate_amplitude, tolerance=0.02)

Approximate dynamic yield stress near zero strain rate.

Selects samples where abs(strain_rate) is small, averages abs(stress), and returns that average as the dynamic yield estimate.

Return type:

float

rheojax.utils.spp_kernels.harmonic_reconstruction(stress, omega, n_harmonics=39, dt=None)

Reconstruct stress signal from harmonic components (Fourier decomposition).

Extracts odd harmonic amplitudes and phases from LAOS stress signal, enabling reconstruction and harmonic ratio analysis.

Parameters:

• stress (Array) – Time-resolved stress signal σ(t) (Pa)

• omega (float) – Fundamental angular frequency ω (rad/s)

• n_harmonics (int) – Number of odd harmonics to extract (default: 5, gives 1ω, 3ω, 5ω, 7ω, 9ω)

• dt (float | None) – Time step. If None, assumes stress spans exactly one period.

Return type:

tuple[Array, Array, Array]

Returns:

• amplitudes (Array) – Harmonic amplitudes [A_1, A_3, A_5, …] (Pa)

• phases (Array) – Harmonic phases [φ_1, φ_3, φ_5, …] (radians)

• reconstructed (Array) – Reconstructed stress from harmonics (Pa)

Examples

>>> import jax.numpy as jnp
>>> from rheojax.utils.spp_kernels import harmonic_reconstruction
>>>
>>> omega = 1.0
>>> t = jnp.linspace(0, 2*jnp.pi, 1000)
>>> sigma = 100.0 * jnp.sin(t) + 20.0 * jnp.sin(3*t + 0.1)
>>>
>>> amps, phases, reconstructed = harmonic_reconstruction(sigma, omega)
>>> # amps[0] ≈ 100.0, amps[1] ≈ 20.0
>>> # phases[0] ≈ 0, phases[1] ≈ 0.1

Notes

• Only odd harmonics (1, 3, 5, …) are physically relevant in LAOS

• Even harmonics indicate asymmetric response (wall slip, etc.)

• I_n/I_1 ratio quantifies nonlinearity strength

rheojax.utils.spp_kernels.power_law_fit(stress, strain_rate, threshold_fraction=0.1)

Log-log fit of sigma = K * abs(strain_rate) ** n over the flowing region.

Returns (K, n, r_squared).

Return type:

tuple[float, float, float]

rheojax.utils.spp_kernels.lissajous_metrics(stress, strain, strain_rate, strain_amplitude, rate_amplitude)

Compute Lissajous-Bowditch diagram derived metrics.

Extracts nonlinearity measures from Lissajous plots including S-factor (stiffening ratio) and T-factor (thickening ratio).

Parameters:

• stress (Array) – Time-resolved stress signal σ(t) (Pa)

• strain (Array) – Time-resolved strain signal γ(t) (dimensionless)

• strain_rate (Array) – Time-resolved strain rate signal γ̇(t) (1/s)

• strain_amplitude (float) – Maximum strain amplitude γ_0 (dimensionless)

• rate_amplitude (float) – Maximum strain rate amplitude γ̇_0 = ω * γ_0 (1/s)

Returns:

Dictionary containing: - G_L: Large-strain modulus (tangent at γ = γ_0) - G_M: Minimum-strain modulus (tangent at γ = 0) - eta_L: Large-rate viscosity (tangent at γ̇ = γ̇_0) - eta_M: Minimum-rate viscosity (tangent at γ̇ = 0) - S_factor: Stiffening ratio (G_L - G_M) / G_L - T_factor: Thickening ratio (η_L - η_M) / η_L

Return type:

dict

Examples

>>> import jax.numpy as jnp
>>> from rheojax.utils.spp_kernels import lissajous_metrics
>>>
>>> omega = 1.0
>>> t = jnp.linspace(0, 2*jnp.pi, 1000)
>>> gamma_0 = 1.0
>>> gamma = gamma_0 * jnp.sin(omega * t)
>>> gamma_dot = gamma_0 * omega * jnp.cos(omega * t)
>>> sigma = 100.0 * gamma + 10.0 * gamma_dot  # Linear viscoelastic
>>>
>>> metrics = lissajous_metrics(sigma, gamma, gamma_dot, gamma_0, gamma_0 * omega)
>>> # S_factor ≈ 0 (linear), T_factor ≈ 0 (linear)

Notes

• S > 0: strain stiffening, S < 0: strain softening

• T > 0: shear thickening, T < 0: shear thinning

• For linear viscoelastic: S = T = 0

rheojax.utils.spp_kernels.zero_crossing_indices(signal)

Find indices of zero-crossings in a signal (robust implementation).

Uses linear interpolation to find precise crossing locations, handling noise-induced multiple crossings via hysteresis filtering.

Parameters:

signal (Array) – Input signal to analyze for zero-crossings

Returns:

Boolean mask of zero-crossing locations (True at crossings)

Return type:

Array

Examples

>>> import jax.numpy as jnp
>>> from rheojax.utils.spp_kernels import zero_crossing_indices
>>>
>>> signal = jnp.sin(jnp.linspace(0, 4*jnp.pi, 100))
>>> crossings = zero_crossing_indices(signal)
>>> # crossings is True at indices where sin crosses zero

Notes

• Returns mask of same length as signal

• Crossing detected when sign(s[i]) != sign(s[i+1])

• Edge cases (exact zeros) handled properly

rheojax.utils.spp_kernels.harmonic_truncation_robustness(amplitudes, n_harmonics_original, n_harmonics_truncated)

Compute truncation error metric for harmonic decomposition.

Quantifies how much signal energy is lost when truncating to fewer harmonics, enabling assessment of reconstruction quality.

Parameters:

• amplitudes (Array) – Full set of harmonic amplitudes [A_1, A_3, A_5, …]

• n_harmonics_original (int) – Original number of harmonics

• n_harmonics_truncated (int) – Number of harmonics to keep after truncation

Returns:

Fraction of total energy retained after truncation (0 to 1)

Return type:

float

Examples

>>> import jax.numpy as jnp
>>> from rheojax.utils.spp_kernels import harmonic_truncation_robustness
>>>
>>> # Fundamental dominant with small 3rd harmonic
>>> amps = jnp.array([100.0, 10.0, 2.0, 0.5, 0.1])
>>> robustness = harmonic_truncation_robustness(amps, 5, 2)
>>> # robustness ≈ 0.99 (most energy in first 2 harmonics)

Notes

• Value near 1.0 indicates safe truncation

• Value < 0.95 suggests significant information loss

• Useful for adaptive harmonic selection

rheojax.utils.spp_kernels.spp_stress_decomposition(stress, strain, strain_rate, strain_amplitude, rate_amplitude)

Decompose total stress into elastic and viscous contributions via linear projection.

This is a Cho-style orthogonal projection (Cho et al. 2005), NOT the full SPP Frenet-Serret decomposition. It projects stress onto normalized strain and strain-rate directions and distributes the residual symmetrically. For the full SPP decomposition (which includes the displacement stress sigma_d via the osculating plane), use spp_fourier_analysis() instead.

Separates σ(t) = σ’(t) + σ’’(t) where: - σ’(t): Elastic (in-phase with strain) component - σ’’(t): Viscous (in-phase with strain rate) component

Parameters:

• stress (Array) – Time-resolved stress signal σ(t) (Pa)

• strain (Array) – Time-resolved strain signal γ(t) (dimensionless)

• strain_rate (Array) – Time-resolved strain rate signal γ̇(t) (1/s)

• strain_amplitude (float) – Maximum strain amplitude γ_0 (dimensionless)

• rate_amplitude (float) – Maximum strain rate amplitude γ̇_0 = ω * γ_0 (1/s)

Return type:

tuple[Array, Array]

Returns:

• sigma_elastic (Array) – Elastic stress contribution σ’(t) (Pa)

• sigma_viscous (Array) – Viscous stress contribution σ’’(t) (Pa)

Examples

>>> import jax.numpy as jnp
>>> from rheojax.utils.spp_kernels import spp_stress_decomposition
>>>
>>> omega = 1.0
>>> t = jnp.linspace(0, 2*jnp.pi, 1000)
>>> gamma_0 = 1.0
>>> gamma = gamma_0 * jnp.sin(omega * t)
>>> gamma_dot = gamma_0 * omega * jnp.cos(omega * t)
>>> # G' = 100 Pa, G'' = 50 Pa
>>> sigma = 100.0 * gamma + 50.0 * gamma_dot / omega
>>>
>>> sigma_e, sigma_v = spp_stress_decomposition(
...     sigma, gamma, gamma_dot, gamma_0, gamma_0 * omega
... )
>>> # sigma_e ≈ 100 * gamma (elastic)
>>> # sigma_v ≈ 50 * gamma_dot / omega (viscous)

Notes

• This is a linear projection decomposition, exact for sinusoidal signals. For nonlinear responses, the residual (displacement stress) is split equally between elastic and viscous components.

• For the full SPP decomposition that separately tracks the displacement stress, use spp_fourier_analysis().

• For linear viscoelastic: σ_e = G’ * γ, σ_v = G’’ * γ / ω

• Decomposition satisfies σ = σ_e + σ_v at all times

rheojax.utils.spp_kernels.numerical_derivative(signal, dt, order=1, step_size=1)

Alias for numerical_derivative_4th_order(). Kept for backwards compatibility and MATLAB SPPplus naming parity.

Return type:

Array

rheojax.utils.spp_kernels.numerical_derivative_4th_order(signal, dt, order=1, step_size=1)

Compute numerical derivatives using 4th-order finite differences (MATLAB SPPplus compatible).

Implements the EXACT finite-difference schemes from SPPplus_numerical_v2.m: - 4th-order centered differences in the interior - Forward differences at the beginning boundary - Backward differences at the ending boundary

This matches MATLAB’s “standard” differentiation mode (num_mode=1).

Parameters:

• signal (Array) – Input signal to differentiate (1D array)

• dt (float) – Time step between samples (s)

• order (int) – Derivative order: 1, 2, or 3 (default: 1)

• step_size (int) – Step size k for stencil (default: 1, larger = more smoothing)

Returns:

Numerical derivative of same length as input (4th-order accurate in interior)

Return type:

Array

Notes

MATLAB SPPplus_numerical_v2.m stencils (mode 1):

First derivative (interior, 4th order):

rd = (-f[p+2k] + 8*f[p+k] - 8*f[p-k] + f[p-2k]) / (12*k*dt)

Second derivative (interior, 4th order):

rdd = (-f[p+2k] + 16*f[p+k] - 30*f[p] + 16*f[p-k] - f[p-2k]) / (12*(k*dt)^2)

Third derivative (interior, 4th order):

rddd = (-f[p+3k] + 8*f[p+2k] - 13*f[p+k] + 13*f[p-k] - 8*f[p-2k] + f[p-3k]) / (8*(k*dt)^3)

Forward/backward stencils at boundaries use 2nd-order accurate formulas.

rheojax.utils.spp_kernels.numerical_derivative_periodic(signal, dt, step_size=1)

Compute 1st, 2nd, and 3rd derivatives assuming periodic signal (MATLAB “looped” mode).

For LAOS data where the signal is periodic (steady-state oscillation), this uses centered differences everywhere by wrapping around at boundaries. Matches MATLAB SPPplus_numerical_v2.m “looped” differentiation mode.

Parameters:

• signal (Array) – Periodic input signal (e.g., one or more complete LAOS cycles)

• dt (float) – Time step between samples (s)

• step_size (int) – Step size k for stencil (default: 1)

Return type:

tuple[Array, Array, Array]

Returns:

• d1 (Array) – First derivative

• d2 (Array) – Second derivative

• d3 (Array) – Third derivative

Examples

>>> import jax.numpy as jnp
>>> from rheojax.utils.spp_kernels import numerical_derivative_periodic
>>>
>>> # One complete period of sine wave
>>> omega = 1.0
>>> t = jnp.linspace(0, 2*jnp.pi/omega, 1000, endpoint=False)
>>> dt = t[1] - t[0]
>>> signal = jnp.sin(omega * t)
>>>
>>> d1, d2, d3 = numerical_derivative_periodic(signal, dt)
>>> # d1 ≈ omega * cos(omega*t)
>>> # d2 ≈ -omega^2 * sin(omega*t)
>>> # d3 ≈ -omega^3 * cos(omega*t)

Notes

• Assumes signal represents complete periods (periodic boundary)

• Uses higher-order centered differences for accuracy

• More accurate than standard differentiation for periodic LAOS data

rheojax.utils.spp_kernels.spp_numerical_analysis(strain, stress, omega, dt, step_size=8, num_mode=2)

Perform full SPP analysis using numerical differentiation (MATLAB-compatible).

Implements the numerical SPP workflow from SPPplus_numerical_v2.m: 1. Compute strain rate from strain (or use provided) 2. Compute derivatives of [strain, rate, stress] trajectory 3. Calculate instantaneous G'_t and G''_t via cross-product formula 4. Extract tan(δ)_t, phase angle, and displacement stress

Parameters:

• strain (Array) – Strain signal γ(t) (dimensionless)

• stress (Array) – Stress signal σ(t) (Pa)

• omega (float | Array) – Angular frequency ω (rad/s). Can be scalar or per-sample array.

• dt (float) – Time step between samples (s)

• step_size (int) – Finite difference step size k (default: 8 for Rogers parity)

• num_mode (int) – 1 = edge-aware (forward/backward + centered); 2 = periodic/looped (default).

Returns:

Dictionary containing: - Gp_t: Instantaneous storage modulus G’(t) (Pa) - Gpp_t: Instantaneous loss modulus G’’(t) (Pa) - G_star_t: Instantaneous complex modulus |G*(t)| (Pa) - tan_delta_t: Instantaneous tan(δ)(t) - delta_t: Instantaneous phase angle δ(t) (radians) - disp_stress: Displacement stress (Pa) - eq_strain_est: Equivalent strain estimate

Return type:

dict

Examples

>>> import jax.numpy as jnp
>>> from rheojax.utils.spp_kernels import spp_numerical_analysis
>>>
>>> omega = 1.0
>>> t = jnp.linspace(0, 2*jnp.pi, 1000)
>>> dt = t[1] - t[0]
>>> gamma_0 = 1.0
>>> strain = gamma_0 * jnp.sin(omega * t)
>>> # Linear viscoelastic response
>>> stress = 100.0 * strain + 50.0 * gamma_0 * omega * jnp.cos(omega * t)
>>>
>>> result = spp_numerical_analysis(strain, stress, omega, dt)
>>> # result['Gp_t'] ≈ 100.0 (constant for linear material)

Notes

• Matches MATLAB SPPplus cross-product formulation

• G'_t = -rd_x_rdd[:,0] / rd_x_rdd[:,2]

• G''_t = -rd_x_rdd[:,1] / rd_x_rdd[:,2]

• Works directly with raw experimental data (no Fourier decomposition)

rheojax.utils.spp_kernels.spp_fourier_analysis(strain, stress, omega, dt, n_harmonics=39, n_cycles=1)

Complete SPP analysis using Fourier-based analytical derivatives (MATLAB-compatible).

Implements the full workflow from SPPplus_fourier_v2.m: 1. FFT strain and stress signals 2. Compute phase offset and rotate coefficients 3. Compute derivatives ANALYTICALLY from Fourier coefficients 4. Calculate G’(t), G’’(t) via cross-product formula 5. Extract all SPP metrics including moduli rates and Frenet-Serret frame

This is more accurate than numerical differentiation for noisy data.

Parameters:

• strain (Array) – Strain signal γ(t) (dimensionless)

• stress (Array) – Stress signal σ(t) (Pa)

• omega (float) – Angular frequency ω (rad/s)

• dt (float) – Time step (s)

• n_harmonics (int) – Number of odd harmonics for reconstruction (default: 15)

• n_cycles (int) – Number of complete cycles in data (default: 1)

Returns:

Dictionary containing all SPP metrics: - Gp_t: Instantaneous G’(t) (Pa) - Gpp_t: Instantaneous G’’(t) (Pa) - G_star_t: Complex modulus |G*(t)| (Pa) - tan_delta_t: Loss tangent tan(δ)(t) - delta_t: Phase angle δ(t) (radians) - disp_stress: Displacement stress (Pa) - eq_strain_est: Equivalent strain estimate - Gp_t_dot: Time derivative of G’(t) (Pa/s) - Gpp_t_dot: Time derivative of G’’(t) (Pa/s) - G_speed: Moduli rate magnitude (Pa/s) - delta_t_dot: Phase angle rate (rad/s) - T_vec, N_vec, B_vec: Frenet-Serret frame vectors - strain_recon, stress_recon: Reconstructed waveforms - Delta: Phase offset used

Return type:

dict

Notes

ANALYTICAL derivatives from Fourier series:

f(t) = Σ [An*cos(nωt) + Bn*sin(nωt)] f’(t) = Σ [-nω*An*sin(nωt) + nω*Bn*cos(nωt)] f’’(t) = Σ [-n²ω²*An*cos(nωt) - n²ω²*Bn*sin(nωt)] f’’’(t) = Σ [n³ω³*An*sin(nωt) - n³ω³*Bn*cos(nωt)]

rheojax.utils.spp_kernels.compute_phase_offset(strain, omega, dt, n_cycles=1)

Compute phase offset Delta for aligning strain to start at zero crossing.

This matches MATLAB SPPplus_fourier_v2.m phase offset calculation:

Delta = atan(An1_n(p+1)/Bn1_n(p+1)) if Bn1_n(p+1) < 0: Delta = Delta + pi

Parameters:

• strain (Array) – Strain signal γ(t)

• omega (float) – Angular frequency ω (rad/s)

• dt (float) – Time step (s)

• n_cycles (int) – Number of complete cycles in data (default: 1)

Returns:

Phase offset Delta (radians) to align strain reference

Return type:

float

rheojax.utils.spp_kernels.harmonic_reconstruction_full(strain, strain_rate, stress, omega, n_harmonics=39, n_cycles=1, W_int=None)

Full Fourier-based harmonic reconstruction with phase alignment (MATLAB-compatible).

Implements the complete workflow from SPPplus_fourier_v2.m: 1. FFT all three waveforms (strain, rate, stress) 2. Compute phase offset Delta from strain fundamental 3. Rotate all Fourier coefficients to align with phase reference 4. Reconstruct aligned waveforms

Parameters:

• strain (Array) – Strain signal γ(t) (dimensionless)

• strain_rate (Array) – Strain rate signal γ̇(t) (1/s) - will be normalized by omega

• stress (Array) – Stress signal σ(t) (Pa)

• omega (float) – Angular frequency ω (rad/s)

• n_harmonics (int) – Number of odd harmonics for stress reconstruction (default: 15)

• n_cycles (int) – Number of complete cycles in data (default: 1)

Returns:

Dictionary containing: - Delta: Phase offset (radians) - An_strain, Bn_strain: Aligned strain Fourier coefficients - An_rate, Bn_rate: Aligned rate Fourier coefficients - An_stress, Bn_stress: Aligned stress Fourier coefficients - strain_recon: Reconstructed strain - rate_recon: Reconstructed rate/omega - stress_recon: Reconstructed stress - time_new: Phase-aligned time array

Return type:

dict

Notes

This function matches MATLAB SPPplus_fourier_v2.m coefficient rotation:

An_n[nn+1] = An1_n[nn+1]*cos(Delta/p*nn) - Bn1_n[nn+1]*sin(Delta/p*nn) Bn_n[nn+1] = Bn1_n[nn+1]*cos(Delta/p*nn) + An1_n[nn+1]*sin(Delta/p*nn)

rheojax.utils.spp_kernels.frenet_serret_frame(rd, rdd)

Compute the Frenet-Serret frame (T, N, B) for a 3D trajectory.

The Frenet-Serret frame provides a local coordinate system along the (γ, γ̇/ω, σ) trajectory, useful for understanding the geometry of the nonlinear response.

Parameters:

• rd (Array) – First derivative of response wave [d(γ)/dt, d(γ̇/ω)/dt, d(σ)/dt] Shape: (n_points, 3)

• rdd (Array) – Second derivative of response wave Shape: (n_points, 3)

Return type:

tuple[Array, Array, Array, Array, Array]

Returns:

• T_vec (Array) – Tangent vector (unit vector in direction of motion)

• N_vec (Array) – Principal normal vector (direction of curvature)

• B_vec (Array) – Binormal vector (T × N)

• curvature (Array) – Local curvature κ = |rd × rdd| / |rd|³

• torsion (Array) – Local torsion τ (requires third derivative, returns zeros)

Notes

Formulas (matching MATLAB SPPplus):

T = rd / |rd|
N = -(rd × (rd × rdd)) / (|rd| × |rd × rdd|)
B = (rd × rdd) / |rd × rdd|
κ = |rd × rdd| / |rd|³

rheojax.utils.spp_kernels.yield_from_displacement_stress(disp_stress, strain, strain_rate, Gp_t, delta_t, strain_amplitude, rate_amplitude)

Extract yield stresses from displacement stress curve (SPP methodology).

This implements the Donley et al. (2019) framework for yield stress extraction: - σ_sy (static yield): From displacement stress at G’(t) → 0 transition - σ_dy (dynamic yield): From displacement stress at δ(t) → π/2 transition

This is more physically meaningful than simple geometric extraction.

Parameters:

• disp_stress (Array) – Displacement stress σ_disp = σ - (G’·γ + G’’·γ̇/ω) (Pa)

• strain (Array) – Strain signal γ(t) (dimensionless)

• strain_rate (Array) – Strain rate signal γ̇(t)/ω (normalized, dimensionless)

• Gp_t (Array) – Instantaneous storage modulus G’(t) (Pa)

• delta_t (Array) – Instantaneous phase angle δ(t) (radians)

• strain_amplitude (float) – Maximum strain amplitude γ_0 (dimensionless)

• rate_amplitude (float) – Maximum strain rate amplitude γ̇_0 = ω * γ_0 (1/s)

Returns:

Dictionary containing: - sigma_sy: Static yield stress (Pa) - from G’(t) minima - sigma_dy: Dynamic yield stress (Pa) - from δ(t) → π/2 - yield_strain_sy: Strain at static yield - yield_strain_dy: Strain at dynamic yield - yield_indices_sy: Indices of static yield points - yield_indices_dy: Indices of dynamic yield points - sigma_sy_disp: Static yield from displacement stress peak - sigma_dy_disp: Dynamic yield from displacement stress at zero rate

Return type:

dict

Notes

The SPP framework defines yield stresses based on the displacement stress: - Static yield occurs when the cage structure breaks (G’(t) → 0) - Dynamic yield occurs when flow ceases (δ(t) → π/2)

This differs from simple geometric extraction (stress at extrema) and provides a more physically meaningful interpretation of the yielding transition.

References

G.J. Donley et al., “Time-resolved dynamics of the yielding transition in soft materials”, J. Non-Newton. Fluid Mech. 264, 2019

rheojax.utils.spp_kernels.calculate_loop_energy(stress, strain)

Calculate energy metrics from stress-strain loops via integration.

Computes the dissipated energy (area of the hysteresis loop) and elastic energy metrics, ensuring parity with MATLAB SPP analysis.

Parameters:

• stress (Array) – Stress signal σ(t) (Pa)

• strain (Array) – Strain signal γ(t) (dimensionless)

Returns:

Dictionary containing: - dissipated_energy: Energy dissipated per cycle (Pa) [Loop Area] - dissipated_energy_density: Dissipated energy / π (Pa) - elastic_energy: Maximum stored elastic energy (Pa) [Estimated]

Return type:

dict

Notes

• Dissipated Energy (E_d): Defined as the area enclosed by the Lissajous loop of stress vs strain: E_d = ∮ σ dγ. This represents the viscous energy loss per unit volume per cycle.

• Elastic Energy (E_s): Estimated as the energy stored at maximum strain, approximated by 0.5 * σ(γ_max) * γ_max. This corresponds to the potential energy stored in the elastic structure at peak deformation.

• Dissipated Energy Density: E_d / π. This is a normalized metric often used to compare dissipation across different amplitudes.

The integration uses the trapezoidal rule on the time-ordered data points, which correctly computes the signed area of the closed loop.

rheojax.utils.spp_kernels.differentiate_rate_from_strain(strain, dt, step_size=8, looped=True)

Compute strain rate from strain via numerical differentiation.

Provides a wrapped (periodic) 8-point stencil path to mirror the MATLAB/Rogers SPPplus implementation, while keeping the prior finite difference fallback for non-periodic data.

Parameters:

• strain (Array) – Strain signal γ(t) (dimensionless)

• dt (float) – Time step (s)

• step_size (int) – Finite difference step size k (default: 8, Rogers parity)

• looped (bool) – If True, use periodic derivative (wrapped); otherwise edge-aware.

Returns:

Strain rate γ̇(t) (1/s)

Return type:

Array

Notes

• looped=True + step_size=8 matches SPPplus v2.1 wrapped 8-point rate inference when the rate column is absent.

• looped=False preserves the previous 4th-order finite-difference path.

rheojax.utils.spp_kernels.convert_units(data, from_unit, to_unit)

Convert data between common rheological units.

Parameters:

• data (Array) – Input data array

• from_unit (str) – Source unit (e.g., ‘percent’, ‘mPa’, ‘rad’, ‘deg’)

• to_unit (str) – Target unit (e.g., ‘fraction’, ‘Pa’, ‘rad’, ‘deg’)

Returns:

Converted data

Return type:

Array

Examples

>>> strain_fraction = convert_units(strain_percent, 'percent', 'fraction')
>>> stress_Pa = convert_units(stress_mPa, 'mPa', 'Pa')