StressPy v0.1: Core Developer Guide

stresspy is an early reference implementation of the deterministic numerical core of the Geometric Stress Criterion (GSC). It decomposes a model–data discrepancy into components that are locally accessible and inaccessible through variation of the model's adjustable parameters.

Version 0.1 is intentionally small. It does not fit models, infer the meaning of uploaded data, construct ODE Jacobians, or perform bootstrap calibration. Those capabilities can be added around the core after the central calculation has been independently tested.

1. Mathematical Convention

At a specified parameter point $\hat{\theta}$, let

$$r = y - f(\hat{\theta})$$

be the model–data discrepancy and let

$$J = \left.\frac{\partial f}{\partial x}\right\vert{}_{\hat{\theta}}$$

be the Jacobian with respect to the chosen parameter coordinates $x$.

After applying an observation-space whitening transformation $L$, GSC uses

$$r_W = Lr, \qquad J_W = LJ$$

If $U_r$ contains the retained left singular vectors of $J_W$, then

$$r_{\parallel,W} = U_r U_r^\top r_W, \qquad r_{\perp,W} = r_W - r_{\parallel,W}$$

The reported stresses are

$$S_{\mathrm{total}} = \Vert{}r_W\Vert{}2^2, \qquad S{\parallel} = \Vert{}r_{\parallel,W}\Vert{}2^2, \qquad S{\perp} = \Vert{}r_{\perp,W}\Vert{}_2^2$$

with normal fraction

$$F_{\perp} = \frac{S_{\perp}}{S_{\mathrm{total}}}$$

The minimum-norm local repair is

$$\Delta x = V_r \Sigma_r^{-1} U_r^\top r_W$$

Under the local linear approximation,

$$r - J \Delta x = r_\perp$$

The repair is minimum-norm only in the coordinates used to construct $J$. Its magnitude is therefore coordinate-dependent. Local tangent accessibility also does not imply that the corresponding finite nonlinear step is practical.

2. What the Core Returns

decompose() returns a GSCResult object containing:

Note: When the total discrepancy is exactly zero, normal_fraction is reported as nan, because the ratio is mathematically undefined.

3. Basic Usage

import numpy as np
from stresspy import decompose

J = np.array([
    [1.0, 0.0],
    [1.0, 1.0],
    [0.0, 1.0],
])

predicted = np.array([1.0, 2.0, 1.0])
observed = np.array([1.1, 1.9, 0.9])
sigma = np.array([0.1, 0.1, 0.1])

result = decompose(
    observed=observed,
    predicted=predicted,
    jacobian=J,
    sigma=sigma,
    rank_rtol=1e-8,
)

print(result.rank)
print(result.total_stress)
print(result.tangent_stress)
print(result.normal_stress)
print(result.normal_fraction_pct)
print(result.repair_vector)