The junction graph and the normalisation

The junction graph of a mesh, and the tension solves normalised in each of its connected pieces.

The normalisation (the "gauge"). A mesh determines tension ratios only, so every tension solve adds one condition that fixes the overall scale: the tensions are normalised so that their mean is mean_tension, exactly, in each connected piece of the junction graph. This module calls that condition the gauge, and a connected piece of the graph a component.

Why per piece. A Young-Dupre-family solver couples two interfaces only when they share a triple line. Build the graph with interfaces as nodes and an edge wherever two interfaces meet at a triple line (foambryo.validation.mesh_io.list_junctions). If that graph has k connected components, the tension system decomposes into k blocks each determined only up to its own scale, while mean(gamma) = mean_tension supplies a single global constraint. Ordinary least squares then returns the minimum-norm solution, which drives the k - 1 unpinned blocks towards zero -- giving a gauge-matched relative error of exactly 1.0, silently. This was found on 12 of the 72 ground-truth simulator meshes of 91-166 cells described in BENCHMARKS.md.

A degenerate sub-case. An interface that is a fully closed surface (an engulfed cell, a sealed lumen) has no triple line at all, so it is an isolated node in no junction-graph edge and no equation in any builder's physics rows. Its tension is not merely unpinned in scale -- it is wholly undetermined, and this module reports it as such (float("nan")), never as a number, never silently folded into a component.

The design. Infer per component, normalised per component: each block's tensions are internally consistent and comparable only within the block. That is surfaced in the output, not buried in a docstring -- ComponentTensionSolution carries, for every interface, which component it belongs to and which normalisation was applied to that component; isolated/closed interfaces are reported as undetermined, never as a component with a normalisation and never as a number; and any code that would aggregate tensions across components must refuse rather than silently blend two incompatible scales (require_single_component_solution below).

The normalisation is exact (since 0.7.0). The builders still append a row of ones, sum(gamma) = n * f(mean_tension); solve_gauged_least_squares reads its value and replaces the row by solve_with_exact_gauge: a gauge row weighted to the physics rows' median norm (the best-conditioned weight), then an exact rescale (homogeneous rows) or shift (difference rows), so that each component's tensions meet their gauge exactly. Until 0.6.0 the row of ones was a soft least-squares row of weight 1. Because the physics rows are homogeneous, its weight changed only the overall scale of the solution, never a tension ratio, but it missed the requested mean (by ~1e-6 for YoungDupre, up to 80 % for Equilibrium) and raised the condition number by up to three orders of magnitude (2.8e7 against 311 on a 1 866-cell tissue's YoungDupre system), which is what limited the iterative solve's accuracy to ~1e-6 of the tension scale.

Scope: which methods. This module decomposes any solver whose assembled system is [physics rows; one trailing global gauge rowsum(gamma) = n * f(mean_tension)], which is every method in foambryo.tension_inference and foambryo.per_edge_tension_inference except:

• TensionComputationMethod.Equilibrium, whose physics rows are keyed by geometric mesh edges (len(edge.twin) != 1), which include quadrijunction edges that mesh_io.list_junctions does not enumerate (it only finds edges bordering exactly three regions) -- so this module's triple-line-only graph would not necessarily match Equilibrium's actual row coupling. A known limitation: Equilibrium is normalised over the whole mesh at once, never decomposed;
• the ForceBalance* and Variational family, whose operator additionally carries area and volume gradients tied to shared cells and is column-coupled through the pressure block rather than through the triple-junction graph alone -- measured not to fragment the same way, so these are normalised over the whole mesh at once and never decomposed here.

SparseSolveFallbackWarning

Bases: UserWarning

lsqr did not converge on a large linear system.

The dense least-squares solve was run instead; or, when the dense solve would not fit in the memory available, lsqr's last iterate is returned and the message says it is not converged.

MalformedTripleLineError

MalformedTripleLineError(edge: tuple[int, int] | None, triple: tuple[int, int, int] | None, missing_interface: Interface)

Bases: RuntimeError

A triple-line edge names an interface pair absent from the mesh's own interface list.

Structurally, this should not be reachable after mesh_io.list_junctions's stricter admission rule: a genuine triple line's three sub-interfaces are, by construction, its three incident triangles' own region pairs (see that function's docstring). This is a named, catchable safety net standing in for what used to be a bare KeyError naming an internal union-find dict key -- for any future change to list_junctions's guarantee, or any caller that builds junctions/edge_to_triple some other way.

Attributes:

Name	Type	Description
edge		The mesh edge (u, v), u < v, that synthesized triple -- None if no edge in edge_to_triple could be matched to it.
triple		The region-triple (a, b, c) the edge synthesized.
missing_interface		The interface pair from triple that has no key in the union-find parent map, i.e. does not occur as any triangle's actual region pair.

Build the error from the offending edge, its synthesized triple, and the missing interface.

JunctionGraphPartition dataclass

JunctionGraphPartition(interfaces: tuple[Interface, ...], components: tuple[tuple[Interface, ...], ...], component_of: dict[Interface, int], isolated: tuple[Interface, ...])

Connected-component partition of a mesh's interfaces by its triple-junction graph.

Attributes:

Name	Type	Description
interfaces	tuple[Interface, ...]	Every interface present on the mesh (mesh_io.list_interfaces).
components	tuple[tuple[Interface, ...], ...]	One tuple of interface keys per connected component, ordered largest first (ties broken by the component's smallest key); a mesh with a fully connected junction graph has exactly one component containing every touched interface.
component_of	dict[Interface, int]	Map interface -> component index, for every interface except the isolated ones (which are never assigned a component -- see isolated).
isolated	tuple[Interface, ...]	Interfaces touching no triple line at all (closed surfaces). Disjoint from every entry of components.

n_components property

n_components: int

Number of proper (gaugeable) components -- excludes isolated interfaces.

fragmented property

fragmented: bool

True iff a single global mean(gamma) = mean_tension normalisation no longer pins one scale.

This is the case whenever there is more than one component, or at least one isolated interface -- an isolated interface is already a failure of the single-gauge assumption even if the rest of the mesh's junction graph is fully connected.

JunctionGraphGaugeWarning

Bases: UserWarning

Emitted when a mesh's triple-junction graph is in several pieces or has an undetermined interface.

Fires whenever JunctionGraphPartition.fragmented is true: more than one component (a connected piece of the junction graph, normalised on its own), at least one undetermined interface (no surviving triple line, e.g. a fully closed surface), or both. Tensions in different components sit on different, arbitrary relative scales -- comparing or averaging them directly is meaningless -- and an undetermined interface is reported as nan, never a number. Emitted by infer_tensions, infer_tensions_and_residuals, infer_forces and infer_tensions_over_seeds, whose flat maps hold the tensions of every component side by side; use foambryo.tension_inference.infer_tensions_by_junction_component to get each interface's component explicitly instead of blending them.

ComponentTensionSolution dataclass

ComponentTensionSolution(tensions: NDArray[float64], dict_tensions: dict[Interface, float], component_of: dict[Interface, int | None], gauge_notes: dict[int, str], undetermined: tuple[Interface, ...], lstsq_residuals: NDArray[float64], fragmented: bool, n_components: int, solver_notes: dict[int, dict[str, object]] = dict())

One inference's tensions, decomposed by the junction graph.

When fragmented is False this is the solve of the whole mesh with one normalisation, mean(gamma) = mean_tension met exactly. When fragmented is True, each connected component of the interfaces' triple-junction graph was solved with its own mean(gamma) = mean_tension normalisation, met exactly; a component's tensions are therefore comparable only to each other, never across components, and gauge_notes says so per component. Isolated interfaces (no triple line at all) are undetermined: dict_tensions maps them to float("nan"), and component_of maps them to None.

Attributes:

Name	Type	Description
tensions	NDArray[float64]	(n_interfaces,), aligned with the keys the caller passed in; float("nan") at every undetermined interface.
dict_tensions	dict[Interface, float]	The same, keyed by interface. float("nan") for undetermined interfaces -- never a number.
component_of	dict[Interface, int | None]	Interface -> component index, or None for an undetermined interface.
gauge_notes	dict[int, str]	Component index -> a one-line description of the gauge applied to it.
undetermined	tuple[Interface, ...]	Interfaces with no triple line at all (closed surfaces).
lstsq_residuals	NDArray[float64]	Concatenation of each component's squared physics residual (a length-1 array per component; the normalisation, met exactly, adds none). When fragmented is False, the single whole-mesh solve's.
fragmented	bool	Whether the junction graph required per-component gauging.
n_components	int	Number of gauged components (1 when not fragmented).
solver_notes	dict[int, dict[str, object]]	Component index -> which least-squares solver ran on it and, for the sparse solver, its stop code, iteration count, residual norm and condition estimate (_solve_least_squares). Component 0 is the whole mesh when not fragmented. Empty for solutions not built by solve_gauged_least_squares, and for LogLami, which rebuilds its solution from the one this function returned.

GaugeNotExactWarning

Bases: UserWarning

A gauge could not be made exact: the solution is the well-scaled soft row's.

Raised only for a system whose physics rows are neither homogeneous (zero right-hand side) nor differences (each row sums to zero): no method of this package builds one.

warn_on_junction_graph_fragmentation

warn_on_junction_graph_fragmentation(fragmented: bool, n_components: int, isolated: tuple[Interface, ...], stacklevel: int = 3) -> None

Emit JunctionGraphGaugeWarning when a mesh's junction-graph gauging was fragmented.

Called only from the entry points that drop the underlying gauge structure on the floor -- mirrors second_order_angle_reading.warn_on_partial_coverage's precedent: the report/solution-returning entry point (infer_tensions_by_junction_component) hands the caller the ComponentTensionSolution this warning would have described, and stays silent; infer_tensions_and_residuals (and, through it, infer_tensions) does not, so it calls this instead.

Takes the three scalars directly, rather than a JunctionGraphPartition, so a caller that already solved the system (and therefore already has a ComponentTensionSolution or YoungDupreLocalResult carrying these same three fields) does not need to recompute partition_from_mesh a second time just to decide whether to warn.

Parameters:

Name	Type	Description	Default
fragmented	bool	Whether the junction graph required per-component gauging.	required
n_components	int	Number of gauged components (1 when not fragmented).	required
isolated	tuple[Interface, ...]	Interfaces with no triple line at all.	required
stacklevel	int	passed through to warnings.warn, so the warning points at the caller's own line rather than at this module. Defaults to 3.	3

partition_by_junction_graph

partition_by_junction_graph(f: NDArray[int64], regions: NDArray[int64]) -> JunctionGraphPartition

Connected components of the triple-junction graph, from raw mesh arrays.

One node per interface (mesh_io.list_interfaces); an edge between two interfaces whenever they co-occur in a triple line (mesh_io.list_junctions gives each triple line's three constituent interfaces, all C(3,2) = 3 pairs among them become edges).

Parameters:

Name	Type	Description	Default
f	NDArray[int64]	Mesh triangles (vertex indices), shape (nf, 3).	required
regions	NDArray[int64]	Per-triangle region pair, shape (nf, 2).	required

Returns:

Name	Type	Description
JunctionGraphPartition	JunctionGraphPartition	The partition. A "healthy" mesh (every interface reachable from every other through some chain of shared triple lines, no closed interface) returns exactly one component holding every interface and an empty isolated.

Raises:

Type	Description
MalformedTripleLineError	If a triple line names an interface pair that is not present in mesh_io.list_interfaces(regions). mesh_io.list_junctions's admission rule (three interfaces on exactly three incident triangles) should make this unreachable in practice; this is the defensive, named replacement for what used to be a bare KeyError.

partition_over_surviving_junctions

partition_over_surviving_junctions(f: NDArray[int64], regions: NDArray[int64], surviving: Collection[tuple[int, int, int]]) -> JunctionGraphPartition

:func:partition_by_junction_graph, restricted to the triple lines still carrying a row.

An estimator that declines some of its observations no longer couples the interfaces those observations coupled, so the gauge structure has to be read off the surviving junction lines rather than off the mesh. An interface left touching no surviving line is then isolated — exactly what this module already reports as float("nan") for a fully closed surface — and a survivor graph that falls into several pieces is gauged per component by :func:solve_gauged_least_squares, as a fragmented mesh already is.

surviving only ever removes lines: a triple that this mesh's own mesh_io.list_junctions does not admit as a genuine triple line is absent from the graph here exactly as it is absent from :func:partition_by_junction_graph's, so an estimator whose rows outrun that admission rule is treated identically with and without refusal.

Parameters:

Name	Type	Description	Default
f	NDArray[int64]	Mesh triangles (vertex indices), shape (nf, 3).	required
regions	NDArray[int64]	Per-triangle region pair, shape (nf, 2).	required
surviving	Collection[tuple[int, int, int]]	sorted region-triples whose row survived.	required

Returns:

Name	Type	Description
JunctionGraphPartition	JunctionGraphPartition	the partition of the same interfaces by the surviving lines alone. Passing every triple the mesh has reproduces :func:partition_by_junction_graph exactly.

partition_from_mesh

partition_from_mesh(mesh: object) -> JunctionGraphPartition

Convenience: partition a DcelData mesh directly (mesh.f's vertex/label columns).

Parameters:

Name	Type	Description	Default
mesh	object	A foambryo.dcel.DcelData (typed as object to avoid a hard import cycle; only mesh.f, an (nf, 5) array of [v0, v1, v2, label_a, label_b], is used).	required

Returns:

Name	Type	Description
JunctionGraphPartition	JunctionGraphPartition	See partition_by_junction_graph.

gauge_row_weight

gauge_row_weight(physics: NDArray[float64] | csr_matrix, n_unknowns: int) -> float

Weight of the gauge row: its norm is the median non-zero physics-row norm (to a power of two).

The row is weight * 1^T, norm weight * sqrt(n). Matching it to the physics rows' typical norm puts the system at the bottom of the condition number's V in the row weight (measured on a 1 866-cell tissue's Young-Dupre system: 311 against 2.8e7 for the unit row). 1.0 when no physics row has an entry.

solve_with_exact_gauge

solve_with_exact_gauge(physics: NDArray[float64] | csr_matrix, physics_b: NDArray[float64], gauge_total: float, *, system: str = 'tension system', use_sparse: bool | None = None, scale_columns: bool = False) -> tuple[NDArray[np.float64], NDArray[np.float64], dict[str, object]]

Least squares on physics @ x = physics_b with the gauge sum(x) = gauge_total met exactly.

The gauge row is appended with the weight gauge_row_weight (its norm equal to the median physics-row norm, where the system is best conditioned), the system is solved by _solve_least_squares (dense or sparse as the options say), and the gauge is then made exact without moving the solution off the physics rows' least-squares set:

• homogeneous physics rows (physics_b == 0: every Young-Dupre-family, cotangent, Lami, per-edge and Equilibrium system): the solution of a soft gauge row of any weight is the exact-gauge solution times one factor (Sherman-Morrison; the two were measured to agree to 1e-13), so it is rescaled to sum(x) = gauge_total. Tension ratios are those of the exact constrained least-squares solution;
• difference rows (each physics row sums to zero, e.g. LogLami in log space): adding a constant to every unknown leaves every physics residual unchanged, so the solution is shifted to meet the gauge;
• otherwise (no builder of this package): the soft row's solution is returned with a GaugeNotExactWarning.

Parameters:

Name	Type	Description	Default
physics	NDArray[float64] | csr_matrix	(n_rows, n) physics rows, dense or CSR, without any gauge row.	required
physics_b	NDArray[float64]	(n_rows,) their right-hand side.	required
gauge_total	float	the value sum(x) must take (n * mean_tension for the tension methods).	required
system	str	label, for messages and diagnostics. Keyword-only.	'tension system'
use_sparse	bool | None	force a path, see _solve_least_squares. Keyword-only.	None
scale_columns	bool	see _solve_least_squares. Keyword-only.	False

Returns:

Name	Type	Description
tuple	NDArray[float64]	(x, residuals, record). residuals is the squared 2-norm of the physics rows'
	NDArray[float64]	residual, physics @ x - physics_b, as a length-1 array: the gauge row is met exactly and
	dict[str, object]	contributes nothing. record is _solve_least_squares's, plus gauge (how the gauge was
	tuple[NDArray[float64], NDArray[float64], dict[str, object]]	made exact) and gauge_row_weight.

Nothing is raised beyond _solve_least_squares's exceptions; a zero-sum solution of homogeneous rows (no scale to rescale) is returned as solved, with a GaugeNotExactWarning.

solve_gauged_least_squares

solve_gauged_least_squares(matrix_m: NDArray[float64] | csr_matrix, vector_b: NDArray[float64], keys: Sequence[Interface], mean_tension: float, partition: JunctionGraphPartition, *, system: str = 'tension system') -> ComponentTensionSolution

Solve matrix_m @ x = vector_b, gauging per junction-graph component if fragmented.

matrix_m must be exactly the shape every builder in foambryo.tension_inference produces: (n_rows, len(keys)), whose last row is the single global gauge sum(gamma) = vector_b[-1] (vector_b[-1] is len(keys) * mean_tension for every method except LogLami, which works in log-space -- this function reads vector_b[-1] directly rather than recomputing it, so it is agnostic to which convention the caller used). The row itself is not used: each component is solved by solve_with_exact_gauge, which meets its share of the gauge exactly (sum(gamma) over the component = vector_b[-1] * n_component / n, i.e. mean(gamma) = mean_tension). The residuals returned are the physics rows' (the gauge contributes none).

Each component's solve is dense or sparse as foambryo.linear_solves's options in force say (by default: dense below SPARSE_SOLVE_MIN_UNKNOWNS interfaces, lsqr at or above), see _solve_least_squares.

Parameters:

Name	Type	Description	Default
matrix_m	NDArray[float64] | csr_matrix	Assembled system, gauge row last. A sparse system is never densified as a whole: only the components that take the dense path are.	required
vector_b	NDArray[float64]	Right-hand side, gauge value last.	required
keys	Sequence[Interface]	Interface key per column of matrix_m, in order.	required
mean_tension	float	The gauge value the caller requested (only used for the human-readable gauge_notes; the actual constraint value is read from vector_b[-1]).	required
partition	JunctionGraphPartition	This mesh's junction-graph partition. Must cover exactly keys (checked).	required
system	str	what the system is, for messages and diagnostics. Keyword-only.	'tension system'

Returns:

Name	Type	Description
ComponentTensionSolution	ComponentTensionSolution	See its docstring.

Raises:

Type	Description
ValueError	If partition's interfaces do not exactly match keys, or (defensive, fragmented case only) if a physics row's nonzero columns span more than one component -- meaning the caller's row coupling does not match the triple-junction graph, and this decomposition must not be trusted for that builder.

require_single_component_solution

require_single_component_solution(solution: ComponentTensionSolution) -> int

Refuse to let a caller silently aggregate a fragmented solution across components.

Call this at the top of any function that reduces one mesh's tensions to a single number (a mean, a single gauge-matching scale, a pooled error). It raises rather than returning a number computed by blending two components that are not on the same scale.

Parameters:

Name	Type	Description	Default
solution	ComponentTensionSolution	A solve from solve_gauged_least_squares.	required

Returns:

Name	Type	Description
int	int	0, the sole component index, when solution is not fragmented.

Raises:

Type	Description
ValueError	If solution.fragmented is True -- i.e. more than one gauge component, or at least one undetermined interface, is present.