settle.py

"""The stopping rule of the settle loops: when a relaxation run has come to rest.

It is the rule of the hanging-bar settle loop and of the material test's interior
(run_stress_strain).

The Newton solvers are dynamic integrators. To read a STATIC answer from them, the hanging bar
is released under gravity and run until it stops moving. "Stops moving" must not be read off
the kinetic energy of a single frame: a ringing body has zero velocity at the turning point of
every swing, so a kinetic-energy threshold can fire at the extreme of a swing instead of at
rest. The recorded SemiImplicit run of 2026-06-30 stopped exactly that way -- its kinetic
energy fell from 6.7e-5 J to 1.0e-6 J within four frames while the tip was still swinging
through 0.5 mm.

The rule used instead looks at POSITIONS over a trailing window that is longer than one ringing
period (``params.SETTLE_WINDOW_FRAMES``): a run counts as settled once no node coordinate has
moved by more than a small fraction of the current deflection over the whole window. A turning
point cannot pass this test, because the window spans the swing on both sides of it.

A run can also blow up instead of settling -- an explicit integrator past its stable step does.
Warp keeps the positions finite (Newton caps the particle speed), so a NaN check alone misses it:
``diverged`` also calls a run blown up once a node has moved farther than a given limit (the
performance study uses the bar's own length; no stable answer comes near it).

Pure numpy and no engine, so the rules are exercised by tests/test_settle.py on synthetic traces.
"""

from __future__ import annotations

from collections.abc import Sequence

import numpy as np

from common.arrays import Float64Array, FloatArray


def window_spread(window_q: Sequence[FloatArray] | FloatArray) -> float:
    """Largest peak-to-peak excursion of any node coordinate over the window [m].

    Args:
        window_q: W position arrays of shape (n, 3) [m], oldest first (a sequence of arrays or one
            array of shape (W, n, 3)).

    Returns:
        The largest peak-to-peak excursion of any node coordinate over the window [m].
    """
    w: Float64Array = np.asarray(window_q, dtype=np.float64)
    return float(np.ptp(w, axis=0).max())


def is_settled(
    window_q: Sequence[FloatArray] | FloatArray, rest_q: FloatArray, tol_rel: float, tol_abs: float
) -> tuple[bool, float]:
    """Whether the run is at rest over the window, and the spread that decided it.

    Settled means ``spread <= max(tol_rel * u_max, tol_abs)``, where ``u_max`` is the largest
    nodal displacement in the newest frame -- a relative tolerance, so a 40 mm and a 140 mm
    answer are held to the same 0.1 %, with an absolute floor for a body that barely deforms.

    Args:
        window_q: W position arrays of shape (n, 3) [m], oldest first.
        rest_q: Rest positions, shape (n, 3) [m].
        tol_rel: Tolerance as a fraction of the largest nodal displacement in the newest frame.
        tol_abs: Absolute tolerance floor [m].

    Returns:
        (settled, spread): whether the window is at rest, and its ``window_spread`` [m].
    """
    w: Float64Array = np.asarray(window_q, dtype=np.float64)
    spread: float = window_spread(w)
    u_max: float = float(np.abs(w[-1] - np.asarray(rest_q, dtype=np.float64)).max())
    return spread <= max(tol_rel * u_max, tol_abs), spread


def diverged(q: FloatArray, rest_q: FloatArray, limit: float) -> bool:
    """Whether a run has blown up: a non-finite position, or a node farther than ``limit`` from rest [m].

    Args:
        q: Current positions, shape (n, 3) [m].
        rest_q: Rest positions, shape (n, 3) [m].
        limit: Largest displacement from rest a stable answer can have [m].

    Returns:
        ``True`` if any position is NaN / inf or any node is more than ``limit`` from its rest position.
    """
    q = np.asarray(q, dtype=np.float64)
    return not np.isfinite(q).all() or float(np.abs(q - np.asarray(rest_q, dtype=np.float64)).max()) > limit