Equilibrium Stability

equilibrium_stability

Eager local linear-stability evidence for candidate equilibria.

The diagnostic in this module is deliberately dense, host-only and local. It does not change :func:op_engine.steady_state, certify nonlinear or global stability, or imply that every Array API namespace provides an eigensolver.

EquilibriumDomain

Bases: StrEnum

Whether declared one-sided state constraints affect interpretation.

LocalSpectrum

Bases: StrEnum

Strict classification of an assessed reduced linear spectrum.

LocalStabilityConfig(residual_tolerance, spectral_absolute_tolerance, spectral_relative_tolerance, boundary_tolerance, finite_difference_step=None) dataclass

Explicit numerical thresholds for :func:local_stability.

residual_tolerance and boundary_tolerance may be scalars or one value per state coordinate. Spectral tolerances have inverse-time units; their effective band is absolute + relative * max(abs(eigenvalue)).

__post_init__()

Reject invalid scalar controls before evaluating a model.

Raises:

Type	Description
ValueError	A spectral tolerance is negative, non-finite or not scalar.

Source code in src/op_engine/equilibrium_stability.py

def __post_init__(self) -> None:
    """Reject invalid scalar controls before evaluating a model.

    Raises:
        ValueError: A spectral tolerance is negative, non-finite or not scalar.
    """
    for name in ("spectral_absolute_tolerance", "spectral_relative_tolerance"):
        value = getattr(self, name)
        if (
            not isinstance(value, (int, float))
            or isinstance(value, bool)
            or not math.isfinite(value)
            or value < 0
        ):
            msg = f"{name} must be finite and non-negative"
            raise ValueError(msg)

LocalStabilityResult(status, classification, message, state, time, coordinates, coordinate_system, time_unit, jacobian_source, jacobian_reference, model_digest, solver_converged, residual_tolerance, boundary_tolerance, fixed_mask, nonnegative_mask, domain=EquilibriumDomain.UNSPECIFIED, raw_residual=None, active_residual=None, residual_scaled_rms=None, boundary_coordinates=(), finite_difference_steps=None, jacobian=None, full_eigenvalues=None, active_jacobian=None, active_eigenvalues=None, invariant_rows=None, invariant_rank=0, constraint_rank=0, constraint_svd_tolerance=None, tangent_basis=None, tangent_jacobian=None, tangent_eigenvalues=None, conservation_error=None, conservation_tolerance=None, invariant_flow_error=None, invariant_flow_tolerance=None, expected_conservation_neutral_modes=0, observed_active_neutral_modes=None, effective_spectral_tolerance=None, spectral_abscissa=None, leading_eigenvalue=None, leading_timescale=None, timescale_kind=None) dataclass

Evidence and interpretation from one local linearization.

LocalStabilityStatus

Bases: StrEnum

Numerical status, independent of any spectral classification.

StabilityJacobianSource

Bases: StrEnum

Provenance of the dense full-state Jacobian.

StabilityTimescaleKind

Bases: StrEnum

Interpretation of an e-folding timescale.

SymbolicDerivativeArtifact

Bases: Protocol

Structural surface consumed from op_system.SymbolicDerivatives.

evaluate(t, y, **parameters)

Evaluate derivative matrices in the artifact's declared order.

Source code in src/op_engine/equilibrium_stability.py

def evaluate(self, t: object, y: object, **parameters: object) -> _DerivativeValues:
    """Evaluate derivative matrices in the artifact's declared order."""
    ...

local_stability(rhs, candidate, *, config, t=0.0, coordinates=None, coordinate_system='supplied state', fixed=None, invariants=None, nonnegative=False, derivatives=None, parameters=None, jacobian=None, jacobian_reference=None, time_unit=None)

Assess a candidate equilibrium through an eager NumPy linearization.

The classification applies only to an interior, residual-qualified point and, when constraints are declared, only to their tangent manifold. A boundary status deliberately suppresses an interior stability conclusion.

Parameters:

Name	Type	Description	Default
rhs	StabilityRhs	Dense NumPy right-hand side in the supplied coordinates.	required
candidate	NumpyArray | SteadyStateResult	NumPy state or NumPy-backed steady-state result.	required
config	LocalStabilityConfig	Explicit residual, spectral, boundary and difference controls.	required
t	float	Frozen evaluation time. A time-dependent RHS is only diagnosed there.	0.0
coordinates	Sequence[str] | None	Ordered state labels. An op_system artifact supplies these.	None
coordinate_system	str	Human-readable coordinate-system declaration.	'supplied state'
fixed	Sequence[int] | NumpyArray | None	Coordinates held fixed, as indices or a boolean mask.	None
invariants	ArrayLike | None	Declared invariant rows in the supplied coordinates.	None
nonnegative	bool | Sequence[bool] | NumpyArray	Coordinates with a one-sided non-negative domain.	False
derivatives	SymbolicDerivativeArtifact | None	Optional op_system symbolic derivative artifact.	None
parameters	Mapping[str, object] | None	Exact runtime bindings for derivatives.	None
jacobian	StabilityJacobian | None	Optional dense full-state Jacobian callback.	None
jacobian_reference	str | None	Required provenance label for jacobian.	None
time_unit	str | None	Unit attached to any reported e-folding timescale.	None

Returns:

Type	Description
LocalStabilityResult	Partial or complete numerical evidence with an explicit status.

Raises:

Type	Description
TypeError	The state, solver flag or configuration is not eager NumPy.
ValueError	Input declarations, shapes, dtypes or controls conflict.

Source code in src/op_engine/equilibrium_stability.py

def local_stability(  # noqa: C901,PLR0911,PLR0912,PLR0913,PLR0914,PLR0915
    rhs: StabilityRhs,
    candidate: NumpyArray | SteadyStateResult,
    *,
    config: LocalStabilityConfig,
    t: float = 0.0,
    coordinates: Sequence[str] | None = None,
    coordinate_system: str = "supplied state",
    fixed: Sequence[int] | NumpyArray | None = None,
    invariants: ArrayLike | None = None,
    nonnegative: bool | Sequence[bool] | NumpyArray = False,
    derivatives: SymbolicDerivativeArtifact | None = None,
    parameters: Mapping[str, object] | None = None,
    jacobian: StabilityJacobian | None = None,
    jacobian_reference: str | None = None,
    time_unit: str | None = None,
) -> LocalStabilityResult:
    """Assess a candidate equilibrium through an eager NumPy linearization.

    The classification applies only to an interior, residual-qualified point
    and, when constraints are declared, only to their tangent manifold.  A
    boundary status deliberately suppresses an interior stability conclusion.

    Args:
        rhs: Dense NumPy right-hand side in the supplied coordinates.
        candidate: NumPy state or NumPy-backed steady-state result.
        config: Explicit residual, spectral, boundary and difference controls.
        t: Frozen evaluation time. A time-dependent RHS is only diagnosed there.
        coordinates: Ordered state labels. An op_system artifact supplies these.
        coordinate_system: Human-readable coordinate-system declaration.
        fixed: Coordinates held fixed, as indices or a boolean mask.
        invariants: Declared invariant rows in the supplied coordinates.
        nonnegative: Coordinates with a one-sided non-negative domain.
        derivatives: Optional op_system symbolic derivative artifact.
        parameters: Exact runtime bindings for ``derivatives``.
        jacobian: Optional dense full-state Jacobian callback.
        jacobian_reference: Required provenance label for ``jacobian``.
        time_unit: Unit attached to any reported e-folding timescale.

    Returns:
        Partial or complete numerical evidence with an explicit status.

    Raises:
        TypeError: The state, solver flag or configuration is not eager NumPy.
        ValueError: Input declarations, shapes, dtypes or controls conflict.
    """
    if not isinstance(config, LocalStabilityConfig):
        msg = "config must be an explicit LocalStabilityConfig"
        raise TypeError(msg)
    solver_converged: bool | None = None
    state_value: object = candidate
    if isinstance(candidate, SteadyStateResult):
        state_value = candidate.state
    if not isinstance(state_value, np.ndarray):
        msg = "local_stability requires an eager NumPy state; transfer explicitly"
        raise TypeError(msg)
    if isinstance(candidate, SteadyStateResult):
        if not isinstance(candidate.converged, np.ndarray):
            msg = "local_stability requires an eager NumPy convergence flag"
            raise TypeError(msg)
        if candidate.converged.shape != ():
            msg = "SteadyStateResult.converged must be a scalar array"
            raise ValueError(msg)
        solver_converged = bool(candidate.converged.item())
    state = state_value
    if state.ndim != 1 or state.size < 1:
        msg = f"candidate must be a non-empty vector; got shape {state.shape}"
        raise ValueError(msg)
    if state.dtype not in {np.dtype(np.float32), np.dtype(np.float64)}:
        msg = "candidate must have NumPy float32 or float64 dtype"
        raise TypeError(msg)
    if not isinstance(t, (int, float)) or isinstance(t, bool) or not math.isfinite(t):
        msg = "t must be a finite real scalar"
        raise ValueError(msg)
    if not isinstance(coordinate_system, str) or not coordinate_system:
        msg = "coordinate_system must be a non-empty string"
        raise ValueError(msg)
    if time_unit is not None and (not isinstance(time_unit, str) or not time_unit):
        msg = "time_unit must be None or a non-empty string"
        raise ValueError(msg)
    n = state.size
    residual_tolerance = _vector_control(
        config.residual_tolerance, n, "residual_tolerance", positive=True
    )
    boundary_tolerance = _vector_control(
        config.boundary_tolerance, n, "boundary_tolerance", positive=False
    )
    fixed_mask = _selection_mask(fixed, n, "fixed")
    nonnegative_mask = _nonnegative_mask(nonnegative, n)
    invariant_rows = _invariant_rows(invariants, n)
    source, names, bound_parameters, model_digest = _source_and_coordinates(
        derivatives,
        jacobian,
        parameters,
        coordinates,
        n,
        jacobian_reference,
    )

    evidence: dict[str, Any] = {
        "state": state.copy(),
        "time": float(t),
        "coordinates": names,
        "coordinate_system": coordinate_system,
        "time_unit": time_unit,
        "jacobian_source": source,
        "jacobian_reference": jacobian_reference,
        "model_digest": model_digest,
        "solver_converged": solver_converged,
        "residual_tolerance": residual_tolerance,
        "boundary_tolerance": boundary_tolerance,
        "fixed_mask": fixed_mask,
        "nonnegative_mask": nonnegative_mask,
        "invariant_rows": invariant_rows,
    }

    def result(
        status: LocalStabilityStatus,
        message: str,
        classification: LocalSpectrum | None = None,
        **updates: object,
    ) -> LocalStabilityResult:
        return LocalStabilityResult(
            status=status,
            classification=classification,
            message=message,
            **(evidence | updates),
        )

    if not np.isfinite(state).all():
        return result(LocalStabilityStatus.NONFINITE_STATE, "candidate is non-finite")
    infeasible = nonnegative_mask & (state < -boundary_tolerance)
    if np.any(infeasible):
        return result(
            LocalStabilityStatus.INFEASIBLE,
            "candidate violates a declared non-negative coordinate",
            domain=EquilibriumDomain.INFEASIBLE,
        )
    try:
        raw_residual = _strict_numpy_vector(rhs(float(t), state), state, "rhs")
    except (
        ArithmeticError,
        AttributeError,
        LookupError,
        RuntimeError,
        TypeError,
        ValueError,
    ) as error:
        return result(
            LocalStabilityStatus.RESIDUAL_FAILURE,
            f"RHS evaluation failed: {type(error).__name__}: {error}",
        )
    evidence["raw_residual"] = raw_residual
    active = ~fixed_mask
    active_residual = raw_residual[active]
    evidence["active_residual"] = active_residual
    if not np.isfinite(raw_residual).all():
        return result(
            LocalStabilityStatus.NONFINITE_RESIDUAL,
            "candidate residual is non-finite",
        )
    if active_residual.size:
        ratios = active_residual / residual_tolerance[active]
        residual_scaled_rms = float(np.sqrt(np.mean(ratios * ratios)))
    else:
        residual_scaled_rms = 0.0
    evidence["residual_scaled_rms"] = residual_scaled_rms
    if solver_converged is False:
        return result(
            LocalStabilityStatus.SOLVE_FAILURE,
            "input SteadyStateResult did not converge",
        )
    if residual_scaled_rms > 1:
        return result(
            LocalStabilityStatus.NOT_EQUILIBRIUM,
            "active residual exceeds its declared tolerance",
        )

    constrained = nonnegative_mask & active
    boundary = constrained & (np.abs(state) <= boundary_tolerance)
    boundary_coordinates = tuple(
        name for name, hit in zip(names, boundary, strict=True) if hit
    )
    evidence["boundary_coordinates"] = boundary_coordinates
    domain = (
        EquilibriumDomain.BOUNDARY
        if np.any(boundary)
        else EquilibriumDomain.INTERIOR
        if np.any(constrained)
        else EquilibriumDomain.UNSPECIFIED
    )
    evidence["domain"] = domain
    if (
        domain is EquilibriumDomain.BOUNDARY
        and source is StabilityJacobianSource.FINITE_DIFFERENCE
    ):
        return result(
            LocalStabilityStatus.BOUNDARY,
            "central finite differences are not evaluated across a declared boundary",
        )

    try:
        finite_difference_steps = None
        if derivatives is not None:
            values = derivatives.evaluate(float(t), state, **bound_parameters)
            matrix = _strict_numpy_matrix(
                values.state, state, "symbolic derivative artifact"
            )
        elif jacobian is not None:
            matrix = _strict_numpy_matrix(jacobian(float(t), state), state, "jacobian")
        else:
            matrix, finite_difference_steps = _finite_difference_jacobian(
                rhs, float(t), state, config.finite_difference_step
            )
    except (
        ArithmeticError,
        AttributeError,
        LookupError,
        RuntimeError,
        TypeError,
        ValueError,
    ) as error:
        return result(
            LocalStabilityStatus.JACOBIAN_FAILURE,
            f"Jacobian evaluation failed: {type(error).__name__}: {error}",
        )
    evidence["finite_difference_steps"] = finite_difference_steps
    evidence["jacobian"] = matrix
    if not np.isfinite(matrix).all():
        return result(
            LocalStabilityStatus.NONFINITE_JACOBIAN,
            "Jacobian contains non-finite entries",
        )

    free = np.flatnonzero(active)
    active_jacobian = matrix[np.ix_(free, free)]
    evidence["active_jacobian"] = active_jacobian
    try:
        full_eigenvalues = np.linalg.eigvals(matrix)
        active_eigenvalues = np.linalg.eigvals(active_jacobian)
    except np.linalg.LinAlgError as error:
        return result(
            LocalStabilityStatus.EIGENSOLVE_FAILED,
            f"full or active eigensolve failed: {error}",
        )
    evidence["full_eigenvalues"] = full_eigenvalues
    evidence["active_eigenvalues"] = active_eigenvalues
    if (
        not np.isfinite(full_eigenvalues).all()
        or not np.isfinite(active_eigenvalues).all()
    ):
        return result(
            LocalStabilityStatus.NONFINITE_EIGENVALUES,
            "full or active eigenspectrum is non-finite",
        )

    try:
        tangent_basis, constraint_rank, constraint_svd_tolerance = _constraint_basis(
            fixed_mask, invariant_rows
        )
    except (np.linalg.LinAlgError, ValueError) as error:
        return result(
            LocalStabilityStatus.INCONSISTENT_INVARIANTS,
            f"constraint reduction failed: {error}",
        )
    evidence["constraint_rank"] = constraint_rank
    evidence["constraint_svd_tolerance"] = constraint_svd_tolerance
    evidence["tangent_basis"] = tangent_basis

    invariant_rank = invariant_rows.shape[0]
    evidence["invariant_rank"] = invariant_rank
    evidence["expected_conservation_neutral_modes"] = invariant_rank
    if tangent_basis.shape[1] == 0:
        evidence["observed_active_neutral_modes"] = 0
        return result(
            LocalStabilityStatus.NO_DYNAMIC_DIRECTIONS,
            "declared constraints leave no dynamic tangent directions",
        )
    active_spectral_tolerance = (
        config.spectral_absolute_tolerance
        + config.spectral_relative_tolerance * float(np.max(np.abs(active_eigenvalues)))
    )
    evidence["observed_active_neutral_modes"] = int(
        np.count_nonzero(np.abs(active_eigenvalues) <= active_spectral_tolerance)
    )
    if invariant_rank:
        invariant_free = invariant_rows[:, active]
        try:
            active_jacobian_scale = float(np.linalg.norm(active_jacobian, ord=2))
            conservation_error = float(
                np.linalg.norm(invariant_free @ active_jacobian, ord=2)
            )
        except np.linalg.LinAlgError as error:
            return result(
                LocalStabilityStatus.INCONSISTENT_INVARIANTS,
                f"invariant norm evaluation failed: {error}",
            )
        conservation_tolerance = (
            config.spectral_absolute_tolerance
            + config.spectral_relative_tolerance * active_jacobian_scale
        )
        invariant_flow_error = float(np.max(np.abs(invariant_free @ active_residual)))
        invariant_flow_tolerance = float(
            np.max(np.abs(invariant_free) @ residual_tolerance[active])
        )
        evidence.update({
            "conservation_error": conservation_error,
            "conservation_tolerance": conservation_tolerance,
            "invariant_flow_error": invariant_flow_error,
            "invariant_flow_tolerance": invariant_flow_tolerance,
        })
        if conservation_error > conservation_tolerance:
            return result(
                LocalStabilityStatus.INCONSISTENT_INVARIANTS,
                "declared invariant rows do not annihilate the local Jacobian",
            )

    tangent_jacobian = tangent_basis.T @ matrix.astype(np.float64) @ tangent_basis
    evidence["tangent_jacobian"] = tangent_jacobian
    try:
        tangent_eigenvalues = np.linalg.eigvals(tangent_jacobian)
    except np.linalg.LinAlgError as error:
        return result(
            LocalStabilityStatus.EIGENSOLVE_FAILED,
            f"tangent eigensolve failed: {error}",
        )
    evidence["tangent_eigenvalues"] = tangent_eigenvalues
    if not np.isfinite(tangent_eigenvalues).all():
        return result(
            LocalStabilityStatus.NONFINITE_EIGENVALUES,
            "tangent eigenspectrum is non-finite",
        )
    effective_tolerance = (
        config.spectral_absolute_tolerance
        + config.spectral_relative_tolerance
        * float(np.max(np.abs(tangent_eigenvalues)))
    )
    real_parts = np.real(tangent_eigenvalues)
    spectral_abscissa = float(np.max(real_parts))
    evidence["effective_spectral_tolerance"] = effective_tolerance
    evidence["spectral_abscissa"] = spectral_abscissa
    if domain is EquilibriumDomain.BOUNDARY:
        return result(
            LocalStabilityStatus.BOUNDARY,
            "active boundary requires a one-sided tangent-cone analysis",
        )
    if np.any(real_parts > effective_tolerance):
        classification = LocalSpectrum.UNSTABLE
        leading_index = int(np.argmax(real_parts))
        leading_eigenvalue = complex(tangent_eigenvalues[leading_index])
        evidence.update({
            "leading_eigenvalue": leading_eigenvalue,
            "leading_timescale": 1 / spectral_abscissa,
            "timescale_kind": StabilityTimescaleKind.GROWTH,
        })
    elif np.all(real_parts < -effective_tolerance):
        classification = LocalSpectrum.ASYMPTOTICALLY_STABLE
        leading_index = int(np.argmax(real_parts))
        leading_eigenvalue = complex(tangent_eigenvalues[leading_index])
        evidence.update({
            "leading_eigenvalue": leading_eigenvalue,
            "leading_timescale": -1 / spectral_abscissa,
            "timescale_kind": StabilityTimescaleKind.DECAY,
        })
    elif np.all(np.abs(real_parts) <= effective_tolerance):
        classification = LocalSpectrum.NEUTRAL_SPECTRUM
    else:
        classification = LocalSpectrum.NONHYPERBOLIC
    return result(
        LocalStabilityStatus.ANALYZED,
        (
            "classification is local to the declared constraint manifold"
            if invariant_rank or np.any(fixed_mask)
            else "classification is a local interior linearization"
        ),
        classification,
    )