StressPy

Core numerical routines for the Geometric Stress Criterion (GSC).

stresspy decomposes a model--data discrepancy into a component aligned with the model's local parameter-accessible tangent space and a component normal to that space. It also reports numerical rank, singular values, the normal fraction, conditioning and a minimum-norm local repair vector.

What is StressPy?

When a model disagrees with experimental data, conventional goodness-of-fit measures quantify the size of the discrepancy but do not reveal its geometric character. StressPy applies the Geometric Stress Criterion (GSC) to separate it into two components:

• Tangent stress — locally parameter-accessible discrepancy: The part aligned with changes the model can produce, to first order, by varying its parameters near the current point. Local accessibility does not guarantee a practical finite nonlinear repair.
• Normal stress — locally inaccessible discrepancy: The part orthogonal to the retained local parameter-response directions in the selected metric. It cannot be removed to first order at the current point. It can reflect noise as well as model discrepancy; it is not proof of global inadequacy.

StressPy complements goodness-of-fit measures such as RMSE and chi-squared error, and model-selection criteria such as AIC and BIC. These quantify fit or balance fit against complexity; StressPy asks how much of the discrepancy is aligned with locally parameter-accessible response directions.

Version 0.0.8 adds nonlinear repair evaluation and parameter, observation, rank, weighting and joint-condition diagnostics. Calibration against an explicit noise model, introduced in 0.0.7, remains available. StressPy remains an early reference implementation: it constructs Jacobians and supports bootstrap orchestration, but requires a user-supplied fitting function and does not solve ODEs on the user's behalf.

Installation

pip install stresspy

Quick start

import numpy as np
from stresspy import evaluate_gsc

def model(parameters, x):
    intercept, slope = parameters
    return intercept + slope * x

x = np.array([-1.0, 0.0, 1.0, 2.0])
parameters = np.array([1.0, 0.5])
observed = np.array([0.45, 1.10, 1.70, 2.35])
sigma = np.full(observed.size, 0.10)

result = evaluate_gsc(
    model_func=model,
    parameters=parameters,
    observed=observed,
    model_args=(x,),
    sigma=sigma,
)

print("Predictions:", result.predicted)
print("Jacobian:\n", result.jacobian)
print("Tangent stress:", result.tangent_stress)
print("Normal stress:", result.normal_stress)
print("Normal fraction (%):", result.normal_fraction_pct)
print("Numerical rank:", result.rank)
print("Repair vector:", result.repair_vector)

evaluate_gsc evaluates the model, constructs its Jacobian and performs the decomposition in one call. Its default Jacobian method is dependency-free forward finite differencing.

Jacobian construction

The model function must accept the parameter vector as its first argument and return one finite prediction per observation. Additional inputs can be supplied through model_args and model_kwargs.

Forward finite differences are the default:

result = evaluate_gsc(
    model,
    parameters,
    observed,
    model_args=(x,),
    jacobian_method="forward",
)

Central finite differences require twice as many perturbed model evaluations but commonly improve derivative accuracy:

result = evaluate_gsc(
    model,
    parameters,
    observed,
    model_args=(x,),
    jacobian_method="central",
)

StressPy chooses parameter-scaled finite-difference steps from machine precision. A positive scalar or one step per parameter can instead be supplied with step.

Finite differences assume that model outputs are deterministic and locally smooth at the supplied parameter point. Discontinuities, solver failures, stochastic simulations and poorly scaled parameters can make a numerical Jacobian unreliable. Important analyses should be repeated with alternative step sizes or central differences as a sensitivity check.

The Jacobian and repair vector use exactly the coordinates supplied in parameters. To work in log-parameter coordinates, pass log parameters to a model wrapper that exponentiates them before evaluating the underlying model.

The adapters can also be used independently:

from stresspy import finite_difference_jacobian

jacobian = finite_difference_jacobian(
    model,
    parameters,
    method="central",
    model_args=(x,),
)

Optional JAX automatic differentiation

Install the optional dependency with:

pip install "stresspy[jax]"

Then use a JAX-traceable model written with jax.numpy operations:

result = evaluate_gsc(
    jax_model,
    parameters,
    observed,
    jacobian_method="jax",
)

StressPy never silently substitutes finite differences when JAX is explicitly requested. An informative error is raised if JAX is unavailable or the model cannot be differentiated by JAX.

Analysis from a precomputed residual and Jacobian

When predictions and the Jacobian have already been calculated, use analyze:

from stresspy import analyze

residual = observed - predicted
result = analyze(residual, jacobian, sigma=sigma)

Weighting

An unweighted Euclidean analysis requires no additional argument:

result = analyze(residual, jacobian)

Independent observational standard deviations can be supplied with sigma:

result = analyze(residual, jacobian, sigma=sigma)

An optional absolute or quantile-based lower floor can prevent extremely small standard deviations from dominating the observation metric:

absolute_floor = analyze(
    residual,
    jacobian,
    sigma=sigma,
    sigma_floor=0.05,
)

quantile_floor = analyze(
    residual,
    jacobian,
    sigma=sigma,
    sigma_floor_quantile=0.10,
)

Positive diagonal precision weights may be supplied directly. They define the metric sum(weights * residual**2) and are equivalent to sigma = 1 / sqrt(weights):

weights = 1.0 / sigma**2
result = analyze(residual, jacobian, weights=weights)

For correlated observations, supply a positive-definite covariance matrix:

result = analyze(residual, jacobian, covariance=covariance)

Supply only one of sigma, weights, covariance or whitener. Weighting is part of the geometry: different defensible metrics can produce different tangent--normal decompositions and should be reported explicitly.

The discrepancy convention is

[ r = y - f(\hat{\theta}). ]

With observation-space whitening matrix (L), StressPy forms (r_W=Lr) and (J_W=LJ). If (U_r) contains the retained left singular vectors of (J_W), then

[ r_{\parallel,W}=U_rU_r^\top r_W, \qquad r_{\perp,W}=r_W-r_{\parallel,W}. ]

The squared norms give total, tangent and normal stress. The minimum-norm local repair is calculated in the parameter coordinates represented by the supplied Jacobian. Consequently, repair magnitude is coordinate-dependent, and local tangent accessibility does not guarantee a practical finite nonlinear repair.

Principal functions

• monte_carlo_calibration: fixed-geometry calibration against explicit null noise.
• parametric_bootstrap: same-data calibration with user-supplied refitting and Jacobian recomputation for every simulated dataset.
• evaluate_gsc: evaluate a Python model, construct its Jacobian and perform the complete GSC decomposition.
• finite_difference_jacobian: dependency-free forward or central numerical differentiation.
• jax_jacobian: optional forward-mode automatic differentiation using JAX.
• analyze: recommended high-level analysis from a residual and Jacobian, including common weighting and uncertainty-floor options.
• decompose: single tangent--normal decomposition with optional uncertainty or covariance weighting.
• decompose_blocks: joint interrogation of multiple independent observation blocks sharing the same parameter coordinates.
• floor_sigma: explicit uncertainty-floor preprocessing.
• jacobian_to_log_coordinates: conversion of selected Jacobian columns to log-parameter coordinates.

analyze and decompose return an immutable GSCResult. evaluate_gsc returns its subclass GSCEvaluationResult, which adds the parameter vector, observations, predictions, constructed Jacobian, Jacobian method and numerical steps while preserving direct access to every geometric result field.

Interpretation

Normal stress measures discrepancy outside the retained local Jacobian column space in the selected observation metric. It is a local geometric diagnostic, not by itself a calibrated hypothesis test. Conclusions can depend on the chosen weighting, parameter point and singular-value threshold.

Calibration (new in 0.0.8)

Using the result and sigma from the quick start above:

from stresspy import monte_carlo_calibration

calibration = monte_carlo_calibration(
    result.residual,
    result.jacobian,
    noise_sigma=sigma,                 # Assumed sampling noise
    analysis_kwargs={"sigma": sigma}, # Chosen observation metric
    n_resamples=999,
    rng=20260905,
)
print("Normal stress:", calibration.observed_stress)
print("Reference mean:", calibration.reference_mean)
print("95th reference percentile:", calibration.reference_quantiles[0.95])
print("Upper-tail p-value:", calibration.p_value)
print("Monte Carlo tail-probability interval:", calibration.tail_probability_interval)

This holds the Jacobian and weighting fixed. It asks whether the observed normal stress is unusually large under that conditional noise-only reference. It is not a nonlinear refitting test, nor a probability that the model is wrong. The 95th reference percentile is not a confidence bound for structural error.

The sampling noise must always be explicit: provide one of noise_sigma, noise_covariance, or noise_sampler(rng). Analysis weights and uncertainty floors are not automatically treated as a generative noise model.

parametric_bootstrap instead generates data from the fitted null model, refits each dataset using a callback and rebuilds its Jacobian. See CALIBRATION.md and examples/StressPy_calibration.ipynb in the source distribution for the complete contract, assumptions and a nonlinear example.

Both methods report (exceedances + 1) / (n_resamples + 1), counting ties in the upper tail. This finite-simulation correction is not a multiple-testing correction and does not make plug-in bootstrap calibration exact. Reference draws and numerical ranks are retained for inspection. Failed replicates stop the calculation rather than being silently discarded.

Licence and commercial use

StressPy is available under the PolyForm Noncommercial License 1.0.0. It may be used, studied, modified and redistributed for permitted non-commercial purposes under those terms. Commercial use requires separate written permission from the copyright holder.

Citation

If StressPy contributes to academic work, please cite the software and the associated GSC publication when available. Citation metadata is provided in CITATION.cff.

Suggested software citation:

James, D. (2026). StressPy: Core numerical routines for the Geometric Stress Criterion (Version 0.0.8) [Computer software]. https://pypi.org/project/stresspy/

New in 0.0.8

StressPy now tests whether a local repair improves the actual nonlinear model, reports parameter accessibility and observation/group stress contributions, and examines rank/weighting sensitivity. Joint-condition diagnostics compare shared and separate local repairs. Optional bounded repair is available via pip install "stresspy[optimize]".

from stresspy import evaluate_repair, observation_breakdown
repair = evaluate_repair(model_func, parameters, observed)
print(repair.best_alpha, repair.best_stress)
# With a supplied prediction Jacobian J:
parts = observation_breakdown(observed - model_func(parameters), J)

See DIAGNOSTICS.md and examples/StressPy_diagnostics.ipynb in the source release for a complete workflow and interpretation limits. These are local geometric diagnostics; normal stress alone does not prove global model failure.