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_tensionas a hard constraint: minimise||M beta||^2subject toa^T beta = mean_tension(a=1/nion 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 rowsum(gamma) = ni * mean_tension, the convention the angle-reading, cotangent andEquilibriummethods used until 0.6.0 (since 0.7.0 they meet the normalisation exactly, seefoambryo.junction_graph). The normalisation then competes with the physics rows, so the returned mean is only approximatelymean_tension.BayesianMap— Gaussian prior onbeta, 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 oneTensionComputationMethod.ForceBalanceuses.SoftRow: the same condition as one more least-squares row, met only approximately (TensionComputationMethod.ForceBalanceSoftRow).BayesianMap: a Gaussian prior centred onmean_tension, so the mean is deliberately not exactlymean_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 |
cells |
list[int]
|
Cell labels, in the order of |
tensions |
NDArray[float64]
|
Inferred surface tensions, shape |
pressures |
NDArray[float64]
|
Inferred pressures, shape |
tension_std |
NDArray[float64]
|
Per-interface standard deviation, shape |
normalisation |
ForceBalanceNormalisation
|
Which gauge treatment produced this. |
mean_tension_deviation |
float
|
|
residual_norm |
float
|
|
rho |
float
|
|
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]]: |
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 |
1
|
normalisation
|
ForceBalanceNormalisation
|
How the normalisation is imposed.
Defaults to |
Kkt
|
Returns:
| Name | Type | Description |
|---|---|---|
ForceBalanceResult |
ForceBalanceResult
|
Tensions (relative, keyed by interface |
Warns:
| Type | Description |
|---|---|
CutCellPressureWarning
|
on a cropped mesh, naming the cut cells whose pressure is |
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 |
1
|
normalisation
|
ForceBalanceNormalisation
|
How the normalisation is imposed.
Defaults to |
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 |