Error bars and diagnostics

Error bars on inferred tensions, with their measured calibration, and identifiability diagnostics.

infer_tensions_with_uncertainty returns the tensions of infer_tensions (it never changes a point estimate) together with an interval and the diagnostics needed to decide whether to trust it; compute_inference_diagnostics returns the diagnostics alone. Neither has a formulation for crop faces yet: both raise ValueError on a cropped mesh.

Large tissues. Both compute the diagnostics from the dense force-balance operator, three rows per mesh vertex and one column per interface and per cell, and its singular value decomposition. Past DENSE_OPERATOR_ENTRY_LIMIT entries (2**31 - 1, the most scipy's LAPACK can index) that decomposition cannot run, so both raise ValueError at once, before building anything; a mesh of the default reconstruction reaches it at about 1 300 cells of real size. Below the limit the cost grows with the operator: on a generated 1 200-cell foam (0.8 of the limit) infer_tensions_with_uncertainty(mesh, 1, YoungDupre) took about 30 minutes and 46 GiB.

Read this before using an error bar from here. No interval of this module is a calibrated error bar on a reconstructed mesh. The first-order interval of the YoungDupre tensions (AnglePropagation) is a lower bound on the error: measured on 2 400 reconstructions with ground truth, one reconstruction's +- std contained the true tension for 26-54 % of interfaces (nominal 68 %), because the error is mostly the reconstruction's bias, which no propagated scatter measures. The calibration figures below are for the simulator's reference meshes.

The goal here is not "add error bars" — it is produce an interval whose error/std is O(1), or state plainly that we cannot. error/std (the ratio of the actual error against ground truth to the reported standard deviation, gauge-matched) is the calibration metric: 1 is calibrated, >>1 is optimistic (dangerous), <<1 is conservative (safe but uninformative). Measured over the 42 equilibrium cases of benchmarking-dataset (foambryo_benchmarks/run_uncertainty_calibration.py):

model	estimator	median error/std	worst case	verdict
ProbeResidual	force balance	0.233	0.93	O(1), and conservative on 42/42
CoherentResidual	force balance	1.33	29.9	O(1) at the median, not uniformly
EdgeBootstrap	YoungDupreLocal	0.27	1.57	~4x conservative (measured on all 47
cases)
BayesianMap	force balance	0.21	0.61	O(1) by choice of prior, not by measurement
AnglePropagation	YoungDupre	0.015	0.078	~65x conservative here; too narrow on reconstructions
ClusterRobust	force balance	204	1.5e4	optimistic; not recommended
ClusterBootstrap	force balance	216	3e4	optimistic; not recommended
IidGls	force balance	246	5.4e3	optimistic; not recommended (the naive posterior)
LineBootstrap	YoungDupre	1.3e-05	—	~1e5 conservative, weight-confounded; diagnostic only

AnglePropagation's row is on the reference meshes, whose angles are nearly exact. On 2 400 reconstructions of five tissue populations with ground truth, one reconstruction's +- std contains the true tension for only 26-54 % of interfaces (nominal 68 %), +- 2 std for 48-79 % (nominal 95 %): it is not a calibrated interval there, and no single factor makes it one.

So that goal is met: ProbeResidual is an interval with error/std = O(1), and it errs on the safe side on every one of the 42 cases.

Why the naive posterior fails by 250x, and why the fix is one factor. The i.i.d. model says the per-component vertex force residual is N(0, sigma_f^2) with sigma_f estimated as ||M beta|| / sqrt(dof). On a 11 000-vertex mesh that is ~34 000 rows, so the estimator averages the noise down by sqrt(34 000) ~ 185 and claims to know each tension to ~1e-4. It does not: the residual was measured to be systematic — the meshes are the equilibrium of the generator's volume-constraint projection, not of M beta = 0 — so it does not average down at all. Measured, the implied n_eff (the number of independent observations the rows are actually worth, dof / (error/std)^2) has median 0.57 on the 42 and is below 1 on 29 of them: the ~33 000 rows are worth less than one independent observation.

CoherentResidual replaces the probabilistic assumption with a deterministic one, and it is the same linear algebra with sigma_f -> ||M beta||: treat the residual as one coherent perturbation of measured norm and unknown direction. The worst-case-direction first-order bound on the estimator is then

|d gamma_i| <= sigma * sqrt(Sigma_unit_ii),

with Sigma_unit the constrained-GLS shape matrix of foambryo.validation.force_balance.constrained_covariance at sigma_f = 1. That is the n_eff = 1 extreme, and it moves the median calibration from 246 to 1.33 — the whole 250x gap is that one sqrt(dof) factor, nothing else. But it is not uniform: 10 of the 42 cases still come out above 10, worst 29.9.

ProbeResidual closes that, and the mechanism is measured rather than guessed. The remaining failure is self-fitting: ||M beta_hat|| is the residual at the estimate, and beta_hat was chosen to minimise it. The residual at the truth is a median 5.4x and up to 87x larger, and that ratio is what predicts the failures (Spearman 0.80). Worse, the bound fails worst on the best meshes — a mesh close to M's equilibrium lets the fit drive ||M beta_hat|| to almost nothing while the true error stays put (Spearman against rho_free: -0.81). ProbeResidual therefore takes sigma from the residual at an independent tension estimate (PROBE_METHOD, YoungDupre, which never touches M). Measured: the probe residual reproduces the residual at the ground truth to 0.1 % on every case, so this is the oracle bound without an oracle — and probe error can only inflate sigma, so it fails safe. See _probe_residual_norm.

Why resampling does not rescue the i.i.d. model. ClusterRobust (a cluster-robust sandwich) and ClusterBootstrap (resampling whole cells) both model correlation; neither can model a bias. Measured, they move 246 to 204 and 216 — a factor 1.2, not 250. That is the honest reason to prefer the deterministic bound: the dominant term is not variance, so there is nothing for a resampling scheme to find.

Why ForceBalanceMAP's apparently-good 0.21 is not an answer. Its interval's width is set by MAP_PRIOR_TENSION_STD, a hyperparameter documented as "deliberately broad" and never derived. Sweeping it over 0.03-1.0 moves error/std over 2.93-0.061, a log-log slope of -1.03: the interval is the prior. Any calibration it shows was chosen. Both CoherentResidual and ProbeResidual have no such knob — their scale is a measured residual norm.

The negative control behaves differently, and that is a check on the mechanism. On the five non-equilibrium cases excluded from every accuracy claim, the residual is genuinely large and spread over interface interiors, so the implied n_eff is 370 rather than 0.57, IidGls improves to 9.0, and ProbeResidual/CoherentResidual become ~20x conservative. The models fail in the safe direction exactly where the mesh is bad — which is the behaviour one wants and not something that was fitted.

Why the angle family is conservative rather than optimistic. YoungDupre reads triple-line angles and never touches M. Propagating the measured along-line angle spread through the operator over-predicts its error by ~2 orders, because (as measured) the dataset's angle deviations from Neumann lie nearly orthogonal to the Young-Dupre operator's column space. So its interval is safe and nearly useless, and its point estimate (2.3e-05) is far better than any interval we can defend around it.

One caveat, and it is load-bearing. YoungDupre's accuracy is weighting-critical: a factor-~4 change in relative inter-line weights costs two orders of magnitude. LineBootstrap resamples whole triple lines, which is such a reweighting (n copies of a line of weight w act as one row of weight w sqrt(n)), so its width is dominated by weight perturbation rather than by angle noise. It is reported for that reason, not despite it. Nothing here changes the junction-length weighting of any shipped method — every point estimate returned by this module is bit-identical to the corresponding infer_tensions call.

UncertaintyModel

Bases: Enum

How the interval around an inferred tension is produced.

See the module docstring for the measured error/std of each. ProbeResidual is the one recommended for reporting — O(1) and conservative on every one of the 42 equilibrium cases. The others are kept because a calibration claim that cannot be reproduced against its alternatives is not a calibration claim, and because two of them (IidGls, ClusterRobust) are the models a reader would otherwise reach for.

Attributes:

Name	Type	Description
ProbeResidual		CoherentResidual, but scaled by the residual at an independent tension estimate instead of at the fitted one, which removes the self-fitting. Force-balance / Variational family. Default, and the only model measured to be conservative on every one of the 42 cases.
CoherentResidual		Deterministic worst-case-direction bound, sigma_f = ||M beta|| (n_eff = 1). Force-balance / Variational family. O(1) at the median but not uniformly — see the module docstring.
IidGls		the naive constrained-GLS posterior, sigma_f = ||M beta|| / sqrt(dof). Retained for comparison and not for reporting — ~250x optimistic.
ClusterRobust		Cluster-robust (sandwich) covariance, clusters = cells.
ClusterBootstrap		Nonparametric bootstrap resampling whole cells.
AnglePropagation		First-order propagation of the measured along-line angle spread through the Young-Dupre solve. Angle family (YoungDupre), and its default. Not a calibrated interval on reconstructions: a lower bound on the error (one reconstruction's +- std held the true tension for 26-54 % of interfaces, nominal 68 %).
LineBootstrap		Nonparametric bootstrap resampling whole triple lines. Angle family; weight-confounded, see the module docstring.

InferenceDiagnostics dataclass

InferenceDiagnostics(n_interfaces: int, n_cells: int, n_rows: int, rho: float | None, rho_free: float, sigma_gauge: float, sigma_soft: float, soft_mode_fractions: tuple[float, float, float], condition_number: float)

Per-mesh diagnostics saying whether any tension inference on it can be trusted.

These are properties of the mesh and the operator, not of a method, which is what makes them usable on dw3d output where no ground truth exists. rho_free is the one to read first: it was measured that the force-balance / Variational family's error is limited by it (median rho_free 0.012 against a median error 0.021 on the 42), while YoungDupre's is not.

Attributes:

Name	Type	Description
n_interfaces	int	Number of interfaces (tension unknowns).
n_cells	int	Number of cells (pressure unknowns; the exterior is the gauge).
n_rows	int	Rows of M = 3 * n_vertices.
rho	float | None	||P_perp M_gamma gamma|| / ||M_gamma gamma|| at the tensions passed in, pressures marginalised — the fraction of the surface-tension force this mesh leaves unbalanced at those tensions. None if no tensions were given. Gauge-invariant.
rho_free	float	The same quantity with the tensions themselves free — the distance of this mesh from any equilibrium of M. It is the ratio at the minimiser of ||M beta|| under the mean-tension gauge, so rho at given tensions can fall slightly below it, but not by more than the gap between that minimiser and the minimum of the ratio. Gauge-invariant. ~1e-10 for a mesh that minimises M; 0.0007-0.036 on the 42 simulation meshes; 0.26-0.28 on the five non-equilibrium cases excluded from every accuracy claim. On dw3d-reconstructed meshes of the same cases it is about 50x the reference value: median 0.578 (range 0.48-0.72) on the 40-case equilibrium cohort, every case above the 0.2 alarm below, so the force-balance family's numbers on reconstructions stand on a floor no gauge choice can lower.
sigma_gauge	float	Smallest singular value of the column-normalised M — the overall-scale gauge direction. Should be near zero; if it is not, the operator has no scale gauge and something is wrong.
sigma_soft	float	Second-smallest singular value: the worst-determined physical direction. 1 / sigma_soft is the noise amplification (1.2-5.5 on the 42).
soft_mode_fractions	tuple[float, float, float]	How much of that soft direction sits in the (surface tension, line tension, pressure) blocks. The line entry is always 0 here — the harness operator carries no line-tension block, by the package author's design choice.
condition_number	float	cond(M), the operator actually inverted.

identifiability_report

identifiability_report() -> str

Human-readable identifiability and equilibrium summary.

Returns:

Name	Type	Description
str	str	A short multi-line report, safe to print next to an inferred tension table.

TensionUncertaintyResult dataclass

TensionUncertaintyResult(interfaces: list[tuple[int, int]], tensions: NDArray[float64], tension_std: NDArray[float64], method: TensionComputationMethod, model: UncertaintyModel, sigma_scale: float, residual_norm: float, n_rows: int, diagnostics: InferenceDiagnostics, calibration_note: str = '')

An inferred tension field with an interval whose calibration is stated, not assumed.

Attributes:

Name	Type	Description
interfaces	list[tuple[int, int]]	Interface keys, in the order of tensions.
tensions	NDArray[float64]	Inferred tensions — bit-identical to the corresponding infer_tensions call. This module never changes a point estimate.
tension_std	NDArray[float64]	Per-interface interval half-width, in the unit of tensions. For ProbeResidual and CoherentResidual this is a worst-case-direction bound, not a standard deviation; for AnglePropagation a first-order standard error that is a lower bound on the error of a reconstruction, not a calibrated interval; read calibration_note.
method	TensionComputationMethod	Which estimator produced tensions.
model	UncertaintyModel	Which model produced tension_std.
sigma_scale	float	The residual scale the interval was built on.
residual_norm	float	||M beta|| (force-balance family) or ||M gamma - b|| (angle family) at the returned estimate.
n_rows	int	Rows of the operator inverted.
diagnostics	InferenceDiagnostics	Mesh-level identifiability and equilibrium.
calibration_note	str	The measured error/std of this model on the 42 equilibrium cases, and what it means. Always non-empty: an uncalibrated interval must not leave this module without saying so.

dict_tensions property

dict_tensions: dict[tuple[int, int], float]

Map interface key -> inferred tension.

dict_tension_std property

dict_tension_std: dict[tuple[int, int], float]

Map interface key -> interval half-width on that tension.

compute_inference_diagnostics

compute_inference_diagnostics(mesh: DcelData, dict_tensions: dict[tuple[int, int], float] | None = None) -> InferenceDiagnostics

Equilibrium and identifiability diagnostics for a mesh, readable without ground truth.

The rho / sigma_soft diagnostics of the benchmark harness, as a public function, so they can be computed on a mesh with no ground truth (i.e. on dw3d output). Everything here is a property of the mesh and of the force-balance operator; nothing is fitted to a ground-truth tension.

Parameters:

Name	Type	Description	Default
mesh	DcelData	Mesh to analyse, in any length unit.	required
dict_tensions	dict[tuple[int, int], float] | None	Tensions at which to evaluate rho, one per interface (label_1, label_2), label_1 < label_2, on any scale (rho is scale-invariant). Defaults to None, in which case rho is None and only rho_free (which needs no tensions) is reported.	None

Returns:

Name	Type	Description
InferenceDiagnostics	InferenceDiagnostics	See that class; call identifiability_report() for text.

Raises:

Type	Description
ValueError	on a mesh with crop faces, for which these diagnostics have no formulation yet (foambryo.crop_faces), or on a mesh whose dense force-balance operator has more than DENSE_OPERATOR_ENTRY_LIMIT entries (about 1 300 cells of real size; see the module docstring), raised before the operator is built.
KeyError	if dict_tensions misses an interface of the mesh.

infer_tensions_with_uncertainty

infer_tensions_with_uncertainty(mesh: DcelData, mean_tension: float = 1, method: TensionComputationMethod = TensionComputationMethod.ForceBalance, model: UncertaintyModel | None = None, n_resamples: int = DEFAULT_N_RESAMPLES, seed: int = 0) -> TensionUncertaintyResult

Infer tensions and an interval whose calibration has been measured and is reported.

A companion to foambryo.infer_tensions: the point estimate is that of the corresponding infer_tensions call -- bit-identical for ForceBalance and YoungDupre (pinned by tests/test_uncertainty.py), to rounding (within 1e-12 on benchmark case 000) for Variational, whose interval is built on the ForceBalance operator. The YoungDupre tensions are normalised over the whole mesh at once, so they match infer_tensions on a mesh whose junction graph is in one piece. The returned calibration_note states the measured error/std of the interval rather than leaving the reader to assume it is 1.

The interval is not a calibrated error bar on a reconstructed mesh. The first-order AnglePropagation interval of YoungDupre is a lower bound on the error: one reconstruction's +- std held the true tension for 26-54 % of interfaces (nominal 68 %) and +- 2 std for 48-79 % (nominal 95 %), measured on 2 400 reconstructions with ground truth. The force-balance family's error is floored by rho_free on reconstructed meshes (see the module docstring and InferenceDiagnostics).

Only the estimators whose intervals have been calibrated are accepted. That is deliberate: an interval for a method nobody measured would be exactly the failure mode this module exists to close.

Parameters:

Name	Type	Description	Default
mesh	DcelData	Mesh to analyse, in any length unit; not cropped.	required
mean_tension	float	The normalisation: the mean of the returned tensions, met exactly. Defaults to 1.	1
method	TensionComputationMethod	The estimator. One of ForceBalance, Variational (both the force-balance family, same estimator up to the solve route) or YoungDupre (the angle family). Defaults to ForceBalance.	ForceBalance
model	UncertaintyModel | None	How to build the interval. Defaults to None, which selects ProbeResidual for the force-balance family and AnglePropagation for the angle family — in both cases the best-calibrated model measured for that estimator.	None
n_resamples	int	Bootstrap draws, for the bootstrap models (ClusterBootstrap, LineBootstrap). Defaults to DEFAULT_N_RESAMPLES (200).	DEFAULT_N_RESAMPLES
seed	int	Bootstrap RNG seed, recorded so any number is reproducible. Defaults to 0.	0

Returns:

Name	Type	Description
TensionUncertaintyResult	TensionUncertaintyResult	Tensions (keyed by interface (label_1, label_2), label_1 < label_2), the per-interface interval half-width tension_std in the same unit, the mesh's diagnostics and the calibration note.

Raises:

Type	Description
NotImplementedError	If method has no calibrated interval (any method but ForceBalance, Variational and YoungDupre).
ValueError	If model does not apply to method's estimator family, if the mesh has crop faces, for which no interval is calibrated yet (foambryo.crop_faces), or if its dense force-balance operator has more than DENSE_OPERATOR_ENTRY_LIMIT entries (about 1 300 cells of real size, for every method and model, because the diagnostics are always computed; see the module docstring), raised before anything is built.

error_over_std

error_over_std(result: TensionUncertaintyResult, dict_truth_tensions: dict[tuple[int, int], float]) -> dict[str, float]

The calibration metric used throughout this module: |error| / std, gauge-matched, per interface.

Gauge matching is mandatory (the rule every accuracy comparison in this package follows): only tension ratios are physical, so the inferred tensions are rescaled by the single factor making their mean match the truth's, and the same factor is applied to the interval — an interval is in the same units as the quantity it bounds, so rescaling one without the other would be a units error dressed up as a calibration result.

Parameters:

Name	Type	Description	Default
result	TensionUncertaintyResult	The result to score.	required
dict_truth_tensions	dict[tuple[int, int], float]	Ground-truth tension per interface.	required

Returns:

Type	Description
dict[str, float]	dict[str, float]: {"median", "p95", "max", "n", "gauge_matching_scale", "median_std", "median_abs_error"}. Values are nan where the interval is zero.