Skip to content

The force-balance family

Tensions and pressures in one solve, from the unified discrete force-balance operator.

infer_force_balance is the entry point; the ForceBalance* members of foambryo.TensionComputationMethod call it through infer_tensions. Where the eight classical methods in foambryo.tension_inference each impose a local balance written in measured angles (Young-Dupre, Lami, cotangent) or in per-interface area gradients (Variational), this module inverts the global discrete equilibrium condition directly:

E(v) = sum_interfaces gamma_m A_m(v) - sum_cells p_c V_c(v),
equilibrium  <=>  grad_v E = 0  <=>  M beta = 0,     beta = (gamma, p),

with M = [M_gamma | M_p] assembled by foambryo.validation.force_balance.force_balance_operator: row block i is the gradient of E with respect to vertex i, built from the discrete area and volume gradients of its triangles (the documentation's methods note describes it). By the package author's design choice the line-tension block M_t is not used: beta = (gamma, p).

Sign convention. The operator's pressure columns are +dV_c/dv, so the energy above puts pressure in with a minus and the solved pressure-block entries are -p. This module negates them, so ForceBalanceResult.pressures are physical pressures on foambryo's usual convention (exterior medium at 0, higher pressure inside a convex cell). The exterior label is not a column at all — mesh_io.cell_labels drops it — which is what fixes the additive pressure gauge; without that the constant-pressure shift would be an exact null direction of M_p and every solve below would be singular.

The normalisation (the gauge), and the three variants. M beta = 0 is homogeneous: beta = 0 minimises it and only tension ratios are physical, so one scale convention has to be added: the tensions are normalised so that their mean is mean_tension. This module offers three ways of imposing it, ForceBalanceNormalisation:

  • Kkt — mean(gamma) = mean_tension as a hard constraint: minimise ||M beta||^2 subject to a^T beta = mean_tension (a = 1/ni on the surface-tension entries, 0 on the pressures) by solving the KKT system [[2 M^T M, a], [a^T, 0]] @ (beta, lambda) = (0, mean_tension). Exact by construction, and the constraint cannot compete with the physics rows.
  • SoftRow — one extra least-squares row sum(gamma) = ni * mean_tension, the convention the angle-reading, cotangent and Equilibrium methods used until 0.6.0 (since 0.7.0 they meet the normalisation exactly, see foambryo.junction_graph). The normalisation then competes with the physics rows, so the returned mean is only approximately mean_tension.
  • BayesianMap — Gaussian prior on beta, MAP estimate and posterior covariance (map_estimate). The prior anchors the scale softly and shrinks the tensions toward the prior mean, so this variant's gauge is deliberately not exact; report it gauge-matched (the rule every accuracy comparison in this package follows).

All three produce a per-interface standard deviation. For Kkt it is the constrained-GLS covariance with sigma_f estimated from the fit's own residual, so it carries no free hyperparameter and is exactly zero along the fixed gauge direction; see foambryo.validation.force_balance.constrained_covariance for the model and for the caveat that i.i.d. force noise understates the true (spatially correlated) uncertainty.

Reliability on reconstructed meshes -- read every number from this module with this. The operator is inverted on the mesh it is given, and a dw3d-reconstructed mesh is not a force balance for any tension assignment: its rho_free (foambryo.uncertainty) is about fifty times the reference meshes' own value, on every case and independent of which sampler produced it. Measured on the 40 equilibrium cases of the benchmark cohort against their own reference meshes (shipped default reconstruction: dithered point placement, seed 42, min_distance = 3): rho_free median 0.578, range 0.479-0.718, against the references' 0.0118, range 0.00066-0.036 -- a ratio of medians of 49, a median per-case ratio of 50, per-case ratios from 17 to 850 -- and every one of the 40 above the 0.2 alarm threshold; the deterministic variant reads 0.623. The five non-equilibrium cases' references already sit at 0.26 and their reconstructions at 0.59. rho at the Kkt solution equals rho_free by construction, so the solve cannot lower it. Every ForceBalance* tension and pressure inferred on a reconstruction, and every comparison of one reconstruction against another through this family, stands on that floor until the reconstruction itself moves it; the angle-reading family (YoungDupre and relatives) never inverts M and is not bounded this way. Whether the family is repairable on reconstructed meshes is an open scientific question, not addressed here.

ForceBalanceNormalisation

Bases: Enum

How the overall tension scale (the normalisation, or gauge) is fixed when inverting M beta = 0.

  • Kkt: mean(gamma) = mean_tension, exactly, through a KKT system. The default, and the one TensionComputationMethod.ForceBalance uses.
  • SoftRow: the same condition as one more least-squares row, met only approximately (TensionComputationMethod.ForceBalanceSoftRow).
  • BayesianMap: a Gaussian prior centred on mean_tension, so the mean is deliberately not exactly mean_tension (TensionComputationMethod.ForceBalanceMAP).

See the module docstring.

ForceBalanceResult dataclass

ForceBalanceResult(interfaces: list[tuple[int, int]], cells: list[int], tensions: NDArray[float64], pressures: NDArray[float64], tension_std: NDArray[float64], normalisation: ForceBalanceNormalisation, mean_tension_deviation: float, residual_norm: float, rho: float, condition_number: float, sigma_f: float)

Everything one force-balance solve produces, including its own diagnostics.

Attributes:

Name Type Description
interfaces list[tuple[int, int]]

Interface keys, in the order of tensions.

cells list[int]

Cell labels, in the order of pressures; the exterior medium is absent (it is the pressure gauge, fixed at 0).

tensions NDArray[float64]

Inferred surface tensions, shape (ni,).

pressures NDArray[float64]

Inferred pressures, shape (nc,), physical sign, exterior at 0, in the unit of the tensions divided by the mesh's length unit; nan for a cut cell or a cell without an equation on a cropped mesh. A byproduct of the same solve, not a separate inference.

tension_std NDArray[float64]

Per-interface standard deviation, shape (ni,), from the variant's own model (for Kkt, i.i.d. force noise, which understates the uncertainty; see foambryo.uncertainty for measured error bars).

normalisation ForceBalanceNormalisation

Which gauge treatment produced this.

mean_tension_deviation float

|mean(tensions) / mean_tension - 1|. Exactly the quantity the three variants are compared on; ~1e-16 for Kkt by construction.

residual_norm float

||M beta||.

rho float

||M beta|| / ||M_gamma gamma|| at the inferred tensions — the fraction of the surface-tension force this mesh leaves unbalanced under them. For Kkt it is rho_free: the ratio at the minimiser of ||M beta|| under the gauge, which is close to, but not strictly, the lowest ratio any tensions reach (see foambryo.validation.force_balance.effective_surface_tensions).

condition_number float

Condition number of the linear system actually solved (the KKT matrix, the augmented operator, or the posterior precision).

sigma_f float

Force-residual noise scale used for the error bars.

dict_tensions property

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

Map interface key -> inferred tension.

dict_pressures property

dict_pressures: dict[int, float]

Map cell label -> inferred pressure, with the exterior medium at 0.

dict_tension_std property

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

Map interface key -> standard deviation of the inferred tension.

build_force_balance_operator

build_force_balance_operator(mesh: DcelData) -> tuple[NDArray[np.float64], list[tuple[int, int]], list[int]]

Assemble M = [M_gamma | M_p] for a DcelData mesh.

DcelData normalises every triangle so that its stored winding puts the right-hand-rule normal outward from the first (lower) label (dcel._normalize_label_order), which is exactly the convention force_balance_operator assumes for the pressure signs — so the operator can be assembled from mesh.f directly, and gives the same matrix as assembling it from the raw .rec arrays (pinned by test_operator_agrees_between_dcel_and_raw_rec).

Parameters:

Name Type Description Default
mesh DcelData

Mesh to analyse.

required

Returns:

Type Description
tuple[NDArray[float64], list[tuple[int, int]], list[int]]

tuple[NDArray[np.float64], list[tuple[int, int]], list[int]]: (M, interfaces, cells), with M of shape (3 nv, ni + nc) and both key lists sorted.

infer_force_balance

infer_force_balance(mesh: DcelData, mean_tension: float = 1, normalisation: ForceBalanceNormalisation = ForceBalanceNormalisation.Kkt) -> ForceBalanceResult

Infer tensions and pressures by inverting the discrete force-balance operator.

One solve produces both blocks of beta = (gamma, p), so unlike the classical path (tensions from angles, then pressures from a second variational solve) the two are consistent by construction.

Warning -- reconstructed meshes. On dw3d output the returned rho is about fifty times what the same operator reads on the simulator's reference meshes (median 0.578 against 0.0118 on the 40-case equilibrium cohort, all 40 above the 0.2 alarm threshold of InferenceDiagnostics.identifiability_report); no tension assignment balances such a mesh, so the tensions and pressures returned are the least-bad fit to a geometry that is not an equilibrium. Read rho before reading the tensions. See the module docstring.

Cropped meshes (foambryo.crop_faces): the force balances of the vertices on crop faces are left out; the crop interfaces are absent from interfaces; an interface or cell left without an equation is nan; and the pressures of the cut cells are solved for and reported as nan, with a CutCellPressureWarning naming them.

Parameters:

Name Type Description Default
mesh DcelData

Mesh to analyse, in any length unit.

required
mean_tension float

The normalisation: the mean of the returned tensions (exactly for Kkt, approximately for SoftRow, not exactly for BayesianMap). Defaults to 1, matching foambryo.infer_tensions.

1
normalisation ForceBalanceNormalisation

How the normalisation is imposed. Defaults to Kkt (exact by construction).

Kkt

Returns:

Name Type Description
ForceBalanceResult ForceBalanceResult

Tensions (relative, keyed by interface (label_1, label_2), label_1 < label_2), pressures (exterior medium at 0), per-interface standard deviations and diagnostics. There is no base_pressure: add one to every pressure to move the exterior's.

Warns:

Type Description
CutCellPressureWarning

on a cropped mesh, naming the cut cells whose pressure is nan.

infer_tensions_force_balance

infer_tensions_force_balance(mesh: DcelData, mean_tension: float = 1, normalisation: ForceBalanceNormalisation = ForceBalanceNormalisation.Kkt) -> tuple[NDArray[np.float64], dict[tuple[int, int], float], NDArray[np.float64]]

Adapter putting infer_force_balance on the infer_tensions_and_residuals signature.

Parameters:

Name Type Description Default
mesh DcelData

Mesh to analyse.

required
mean_tension float

The normalisation, as for infer_force_balance. Defaults to 1.

1
normalisation ForceBalanceNormalisation

How the normalisation is imposed. Defaults to Kkt.

Kkt

Returns:

Type Description
tuple[NDArray[float64], dict[tuple[int, int], float], NDArray[float64]]

tuple[NDArray[np.float64], dict[tuple[int, int], float], NDArray[np.float64]]: - array of relative tensions; - map interface id (label 1, label 2) -> tension on this interface; - the sum of squared residuals ||M beta||^2, as a length-1 array, which is what the other methods' lstsq returns for an over-determined system.