Linear solves and diagnostics

Solver choice and solve diagnostics for the inference's linear systems.

Every tension and pressure of this package is the solution of a linear least-squares system. Small systems are solved densely (scipy.linalg.lstsq), exactly; large ones sparsely (lsqr, or a sparse direct factorisation), because a dense system grows as the product of its dimensions in memory and as their cube in time: the Young-Dupre system of a 1 866-cell tissue is 49 984 x 11 746, 4.7 GB dense, and at 10 000 cells it would be some 120 GB.

This module holds what the public entry points (infer_tensions*, infer_pressures*, infer_forces) pass down to every solve:

• solver: "auto" (the default) solves a system densely below sparse_threshold unknowns and sparsely at or above it; "dense" and "sparse" force one path whatever the size. Below the threshold, "auto" is bitwise the dense solve. Forcing "dense" on a system whose dense form would not fit in the memory available emits a DenseSolveMemoryWarning with the estimate, before anything is allocated. The tension methods Equilibrium, Variational and ForceBalance* have no sparse path: they are solved densely under "auto", and "sparse" raises ValueError for them. A sparse lsqr solve that does not converge falls back to the dense one with a SparseSolveFallbackWarning (or, when the dense one cannot fit in memory, returns lsqr's last iterate with that warning).
• sparse_threshold: the number of unknowns at which "auto" switches (default foambryo.junction_graph.SPARSE_SOLVE_MIN_UNKNOWNS, 2 000). What counts as "unknowns" is the number of columns of the system the dense path would assemble: the interfaces of one junction-graph component for a tension solve, the materials for a Laplace pressure, the interfaces plus the cells for the variational pressure. One exception: the variational pressure of a cropped mesh is solved densely up to foambryo.pressure_inference.DENSE_CROP_PRESSURE_MAX_ENTRIES matrix entries under "auto", whatever sparse_threshold is.
• diagnostics: False (the default: nothing is computed and each entry point returns its usual values), True (one SolveRecord per solve: solver, shape, rank, residual norm, condition number, time) or "spectrum" (also the extreme singular values of the large, sparse systems, by scipy.sparse.linalg.svds, at a measured extra cost). With diagnostics on, each entry point returns one more value, a SolveDiagnostics, last; save_solve_diagnostics writes it to a JSON file (the scalars) and an NPZ file (the singular values).

They are keyword-only on every entry point that takes them, and an invalid value raises ValueError (validate_solve_options). A lower-level function called directly (for example infer_tensions_young_dupre_local) uses the options in force, the defaults outside any entry point; solve_options puts others in force around such a call.

What the condition numbers are. On a dense least-squares solve, the 2-norm condition number s_max / s_min from the singular values scipy.linalg.lstsq computes anyway (free; inf when the smallest is zero). On an lsqr solve, lsqr's own running estimate of the Frobenius-norm condition number ||A||_F ||A^+||_F (free; an estimate, and an over-estimate of the 2-norm one by up to a factor n). On a sparse direct solve, the 1-norm condition number estimated from the factorisation (Hager-Higham, scipy.sparse.linalg.onenormest; cheap). The source is recorded beside every number: the three are not the same quantity.

How the options travel. A public entry point sets them in a context variable for the duration of its call, and the solve helpers read it, so the options reach every solve without threading an argument through each builder. The context is per thread and per asynchronous task (contextvars); outside any entry point the defaults apply.

DenseSolveMemoryWarning

Bases: UserWarning

A dense solve was forced on a system whose dense form exceeds the memory available.

Emitted before anything is allocated, with the estimate, so that the caller can stop the run (turn it into an error with warnings.simplefilter("error", DenseSolveMemoryWarning)). Only solver="dense" emits it: "auto" never chooses the dense path for a large system.

SolveOptions dataclass

SolveOptions(solver: str = 'auto', sparse_threshold: int | None = None, diagnostics: bool | str = False)

The solver options in force for the solves of one entry-point call.

Attributes:

Name	Type	Description
solver	str	"auto", "dense" or "sparse".
sparse_threshold	int | None	unknowns at which "auto" goes sparse; None means foambryo.junction_graph.SPARSE_SOLVE_MIN_UNKNOWNS, read when the solve runs.
diagnostics	bool | str	False, True or "spectrum".

SolveRecord dataclass

SolveRecord(system: str, solver: str, n_rows: int, n_unknowns: int, rank: int | None = None, residual_norm: float | None = None, condition_number: float | None = None, condition_number_source: str | None = None, seconds: float | None = None, lsqr_istop: int | None = None, lsqr_iterations: int | None = None, spectrum_seconds: float | None = None, spectrum_note: str | None = None, singular_values: NDArray[float64] | None = None)

What one linear solve did. Scalars are JSON-ready; the singular values go to NPZ.

Attributes:

Name	Type	Description
system	str	which system, e.g. "tensions YoungDupre, junction-graph component 0".
solver	str	the solver that produced the solution ("dense lstsq", "sparse lsqr", "dense lstsq after lsqr did not converge", "sparse direct (SuperLU)", "dense inv", ...).
n_rows	int	rows of the system solved (after any reduction, e.g. undetermined cells left out).
n_unknowns	int	columns of the system solved.
rank	int | None	numerical rank, where the solver computes it (dense lstsq: its rank); None on the sparse paths, which do not.
residual_norm	float | None	||b - A x||_2 of the system solved.
condition_number	float | None	see the module docstring for what each source means; None when the path has none.
condition_number_source	str | None	how condition_number was obtained.
seconds	float | None	wall time of the solve, including the diagnostics' own cost.
lsqr_istop,	lsqr_iterations	lsqr's stop code and iteration count (sparse path only).
spectrum_seconds	float | None	time spent on svds for diagnostics="spectrum" (sparse path only).
spectrum_note	str | None	why a requested spectrum is missing or partial, if it is.
singular_values	NDArray[float64] | None	the singular values in decreasing order: all of them on a dense lstsq path; on a sparse path with diagnostics="spectrum", the largest and the smallest SPECTRUM_EXTREME_VALUES, concatenated (a gap, not a full spectrum).

scalars

scalars() -> dict[str, Any]

The record without its singular values, with non-finite floats as strings (JSON-safe).

SolveDiagnostics dataclass

SolveDiagnostics(options: SolveOptions, records: list[SolveRecord] = list())

The records of every solve one entry-point call ran, in the order they ran.

Attributes:

Name	Type	Description
options	SolveOptions	the options in force.
records	list[SolveRecord]	one SolveRecord per solve.

table

table() -> list[dict[str, Any]]

The records' scalars, one dict per solve (what the JSON file holds).

save_solve_diagnostics

save_solve_diagnostics(diagnostics: SolveDiagnostics, path_stem: str | PathLike[str]) -> tuple[Path, Path]

Write solve diagnostics to <path_stem>.json (scalars) and <path_stem>.npz (spectra).

The JSON file holds the options and one object per solve (SolveRecord.scalars, non-finite numbers as the strings "inf"/"nan"). The NPZ file holds, for every solve that has them, its singular values under the key singular_values_<index>, <index> being the solve's position in the JSON list.

Parameters:

Name	Type	Description	Default
diagnostics	SolveDiagnostics	what an entry point returned with diagnostics on.	required
path_stem	str | PathLike[str]	path without extension; parent directories must exist.	required

Returns:

Type	Description
tuple[Path, Path]	tuple[Path, Path]: the JSON path and the NPZ path.

validate_solve_options

validate_solve_options(solver: str, sparse_threshold: int | None, diagnostics: bool | str) -> SolveOptions

Check the three options and return them as a SolveOptions.

Raises:

Type	Description
ValueError	on a solver outside SOLVERS, a diagnostics outside DIAGNOSTICS, or a negative sparse_threshold.

solve_options

solve_options(solver: str = 'auto', sparse_threshold: int | None = None, diagnostics: bool | str = False) -> Iterator[SolveDiagnostics | None]

Put solver options in force for the solves run inside the block.

Used by the public entry points; usable directly around any lower-level call. Yields the SolveDiagnostics the solves inside fill, or None when diagnostics is False.

Raises:

Type	Description
ValueError	see validate_solve_options.

current_solve_options

current_solve_options() -> SolveOptions

The options in force (the defaults outside any entry point).

collecting_diagnostics

collecting_diagnostics() -> bool

Whether the solves now running should build SolveRecords.

spectrum_requested

spectrum_requested() -> bool

Whether diagnostics="spectrum" is in force.

record_solve

record_solve(record: SolveRecord) -> None

Append record to the diagnostics being collected; nothing when none are.

sparse_threshold_in_force

sparse_threshold_in_force() -> int

The sparse_threshold in force, the package default when the caller gave none.

use_sparse_solver

use_sparse_solver(n_unknowns: int) -> bool

Whether the system of n_unknowns unknowns takes the sparse path under the options in force.

require_dense_only_method_is_not_forced_sparse

require_dense_only_method_is_not_forced_sparse(method: str) -> None

Refuse solver="sparse" for a method that has no sparse path, rather than run it densely.

Raises:

Type	Description
ValueError	when solver="sparse" is in force.

available_memory_bytes

available_memory_bytes() -> int | None

The memory a new allocation can use without swapping, as well as this platform tells.

psutil.virtual_memory().available when psutil is installed (it is not a dependency); else MemAvailable from /proc/meminfo on Linux; else the total physical memory (os.sysconf), an upper bound, on platforms that expose no available-memory figure without psutil (macOS). None when nothing is readable.

dense_solve_bytes

dense_solve_bytes(n_rows: int, n_cols: int, copies: int = _DENSE_COPIES) -> int

Estimated bytes of a dense float64 least-squares solve of an n_rows x n_cols system.

warn_if_forced_dense_exceeds_memory

warn_if_forced_dense_exceeds_memory(n_rows: int, n_cols: int, system: str, *, copies: int = _DENSE_COPIES, stacklevel: int = 3) -> None

Emit DenseSolveMemoryWarning when solver="dense" is forced on a system too large for memory.

Call it before allocating the dense matrix. Does nothing under "auto" or "sparse".

Parameters:

Name	Type	Description	Default
n_rows	int	rows of the dense system.	required
n_cols	int	columns of the dense system.	required
system	str	what the system is, for the message.	required
copies	int	how many dense copies of that shape the solve holds at once.	_DENSE_COPIES
stacklevel	int	passed to warnings.warn.	3

condition_from_singular_values

condition_from_singular_values(singular_values: NDArray[float64]) -> float | None

2-norm condition number s_max / s_min (inf if s_min is 0, None if there are none).

extreme_singular_values

extreme_singular_values(operator: spmatrix | sparray, *, symmetric_positive: bool = False) -> tuple[NDArray[np.float64] | None, float, str | None]

The largest and smallest SPECTRUM_EXTREME_VALUES singular values of a sparse operator.

The largest by scipy.sparse.linalg.svds (ARPACK, which="LM"). The smallest first by svds with which="SM"; if ARPACK does not converge within SPECTRUM_MAX_ITERATIONS, by shift-invert at zero (eigsh(..., sigma=0), a sparse LU factorisation) of the operator itself when it is symmetric positive (symmetric_positive: its singular values are its eigenvalues) or of the normal matrix A^T A otherwise. The normal matrix squares the condition number, so those values are accurate only down to about 1e-8 of the largest; the note says which route gave them. Never raises on a failure: what is missing is said in the note.

Returns:

Name	Type	Description
tuple	NDArray[float64] | None	the values in decreasing order (None if none were obtained), the seconds spent, and
	float	a note (None when both ends came from svds).

record_dense_lstsq

record_dense_lstsq(system: str, solver: str, dense: NDArray[float64], vector: NDArray[float64], x: NDArray[float64], residuals: NDArray[float64], rank: int, singular_values: NDArray[float64], start: float, *, lsqr_istop: int | None = None, lsqr_iterations: int | None = None) -> dict[str, Any]

Append the SolveRecord of a dense scipy.linalg.lstsq solve; return its scalar summary.

Reads the solve's own outputs and never alters x. The residual norm is the square root of the residual lstsq returned (a 0-d value for a 1-D right-hand side, shape (1,) from the sparse path) when it is a finite number, else (rank-deficient system) computed as ||vector - dense @ x||.

Parameters:

Name	Type	Description	Default
system	str	which system, see SolveRecord.	required
solver	str	which solver, see SolveRecord.	required
dense	NDArray[float64]	the matrix solved.	required
vector	NDArray[float64]	the right-hand side.	required
x	NDArray[float64]	lstsq's solution.	required
residuals	NDArray[float64]	lstsq's residuals.	required
rank	int	lstsq's rank.	required
singular_values	NDArray[float64]	lstsq's singular values.	required
start	float	time.perf_counter() taken before the solve.	required
lsqr_istop	int | None	lsqr's stop code, when this dense solve is a fallback after lsqr.	None
lsqr_iterations	int | None	lsqr's iterations, likewise.	None

Returns:

Name	Type	Description
dict	dict[str, Any]	n_rows, rank, residual_norm, condition_number, condition_number_source,
	dict[str, Any]	seconds.