From 2e6f3ac58b91fa533a0dc185bde9156887b0d6ef Mon Sep 17 00:00:00 2001 From: Nicholas Ehsan Roy Date: Wed, 7 Oct 2026 12:03:49 +0200 Subject: [PATCH 1/3] fix(sysid): fit_lm reads a small jump of the residual from a long rejected step The floor rule's one-sided test divides a rejected candidate's departure from the linearisation by the mirror image's departure plus the linearisation's own change over the candidate. That change grows with the candidate's length and a jump does not, so a small jump (a late bounce that moves by one time step changes few samples) read from a long candidate fell under the 2**10 that reads a jump: two fits of the stock ball ended converged=True at ratios of 341 and 456, where the loss still fell along the piece (found by the slow random hunt on jaxlib 0.10.2; the same on 0.11.0). A candidate above 2**5 and not above 2**10 is now read a second time: its step is halved down to the one float spacing across which the residual departs most, and that departure is compared with the same mirror-image departure plus the linearisation's change over one spacing, against 2**8. Measured: 1.8e3 and 2.0e3 at the two endings, 3.4e4 or more at the other jump endings of 4,500 drawn fits; at most 52 on smooth endings, none of which is asked (their first reading is at most 21). The mirror image is kept from the whole candidate because one pair of neighbouring points is no measure of rounding: with it taken across one spacing too, smooth floors read 4e3 to 1e6. MADD-ANO-216 (never released) gains the residual case; SYS-080 quotes the second reading. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkCde7g23gTziUjYAvnKD --- CHANGELOG.md | 1 + docs/validation/known_anomalies.yaml | 37 ++++- docs/validation/sysid_fmu_claims.yaml | 15 +- src/maddening/sysid.py | 138 ++++++++++++++++-- .../test_sysid_non_differentiable_residual.py | 137 ++++++++++++++++- tests/property/test_sysid_targeted_search.py | 84 +++++++++-- 6 files changed, 376 insertions(+), 36 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0b77685e..ecb949a9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -418,6 +418,7 @@ guidance; the itemized changes follow. The `[verify]` extra now only pulls `hypothesis`. ### Fixed +- **`fit_lm` no longer reports `converged=True` beside a small jump of the residual read from a long rejected step** (MADD-ANO-216's residual case, never released): the one-sided test divided the jump by the linearisation's change over the whole step, so a late bounce 22 float spacings from the iterate read 341 and 456 against a threshold of `2**10` on the stock ball. A rejected step reading above `2**5` is now read again across the one float spacing where the residual departs most (a few residual evaluations, at the end of the run only); the run ends `converged=False` with the existing warning. Action: none. - **A coupling group under `acceleration="iqn-ils"` or `"iqn-imvj"` returns from `step()` when its iterate leaves float range** (MADD-ANO-233; 0.3.0 and 0.3.1 under `solver="ift"` with `"iqn-imvj"`): the secant least-squares handed LAPACK's SVD a matrix holding an `inf`, and the step never returned. It now runs to `max_iterations` and reports `converged=False`, `residual=inf`, as the other accelerations do; steps of a finite run are unchanged, bit for bit. Action: none; on 0.3.x use `solver="fori"` or `"aitken"`. - **`gradient_relative_error_bound` covers the constants one pass resolves** (MADD-ANO-234, never released): with a constant in it whose tangent at the returned iterate is rounding (the centre or curve of a nonlinearity evaluated on its centre) it read 1.09 for a relative error of 1.134, usable. Such a constant, and a gain whose whole value moves the pass by less than the residual's float floor, is no longer in the bound; the docstring says which constants are. Action: none. - **`GraphManager.step`, `run` and `run_adaptive` no longer store a step that reshapes a state leaf** (`MADD-ANO-220`, in 0.1.0 to 0.3.1): a list given for a scalar constant (`BallNode(initial_velocity=[1.0, 2.0])`), or an external input of another shape at a later step, broadcast a scalar leaf, and the state's checkpoint then did not load after `reset_state()`. The step is refused by name (compared once per trace, host-side: no compiled program changes), as `run_scan` always did; the FMU sidecar and bridge refuse it too. diff --git a/docs/validation/known_anomalies.yaml b/docs/validation/known_anomalies.yaml index 3bb3d7ff..aec3ed28 100644 --- a/docs/validation/known_anomalies.yaml +++ b/docs/validation/known_anomalies.yaml @@ -14716,6 +14716,32 @@ anomalies: the 90 random starts converged above the minimum, and a targeted search of 800 drawn fits of the ball finds no fit converged where the loss still falls along the piece it is on. + + A residual case, found by the slow random hunt on 2026-10-07 (CPU, + jaxlib 0.10.2; the same to the last digit on 0.11.0) and fixed + before any release: a small jump read from a long candidate. The + ratio divides the jump by the linearisation's change over the + candidate, which grows with its length, and a late bounce that moves + by one time step changes few samples. Two fits of the ball (the + hunt's, and one of 4,500 uniformly drawn) ended `converged=True` at + a loss of 0.006 with ratios of 456 and 341 -- a jump of 0.015 in the + residual's norm, 22 float spacings of the elasticity from the + iterate, read from a candidate 35 spacings long -- where the loss + still fell along the piece by 11.8 and 5.6 times what its rounding + explains. A candidate above `2**5` and not above `2**10` is now + read a second time: its step is halved down to the one float spacing + of the coordinates across which the residual departs most from its + linearisation, and it is a jump where that departure is `2**8` times + its mirror image's (the whole candidate's, as before) plus the + linearisation's change over one spacing. Measured: 2.0e3 and 1.8e3 + at those two endings and 3.4e4 or more at the other 261 jump endings + of the 4,500 draws; at most 52 on smooth fits stopped by the floor + rule (140 endings, none of which is asked: their first reading is + at most 21). After the fix: no fit converged + where the loss still falls along its piece among the 4,500 draws or + in the hunt on three seeds (800 examples each), on either jaxlib; + of the 250 endings that existing per-push and slow tests put to + the floor rule (140 of smooth fits), none changed. severity: "major" safety_relevance: "context_dependent" safety_relevance_rationale: > @@ -14743,9 +14769,11 @@ anomalies: test passes because it is stationary on its own piece; and a fit converged at the interior minimum of a smooth piece that is not the lowest (a local minimum, with the lower piece as near as 1e-3 of a - parameter on the ball). The threshold is measured, not derived: a - jump smaller than `2**10` times the linear change of a step of a few - float spacings reads as rounding. The root cause, the missing + parameter on the ball). The thresholds are measured, not derived: + a jump reads as rounding where it is smaller than `2**5` times the + linear change of the rejected step that crossed it, or than `2**8` + times what the residual departs from its linearisation by at that + step's mirror image. The root cause, the missing event-time derivative, is MADD-ANO-021 and stays open. verification: - "tests/core/test_sysid_non_differentiable_residual.py::test_fit_lm_is_converged_on_the_ball_only_at_the_minimum" @@ -14755,6 +14783,9 @@ anomalies: - "tests/core/test_sysid_non_differentiable_residual.py::test_an_iteration_whose_damped_candidates_never_moved_is_asked_again" - "tests/property/test_sysid_targeted_search.py::test_no_wrong_fit_among_the_fixed_draws" - "tests/property/test_sysid_targeted_search.py::test_fit_lm_on_the_ball_is_not_converged_where_the_loss_falls_along_its_piece" + - "tests/property/test_sysid_targeted_search.py::test_fit_lm_reads_a_small_jump_of_the_ball_from_a_long_candidate" + - "tests/core/test_sysid_non_differentiable_residual.py::test_a_small_jump_inside_a_long_candidate_is_read_across_one_spacing" + - "tests/core/test_sysid_non_differentiable_residual.py::test_a_smooth_residual_read_a_second_time_is_still_not_a_jump" github_issue: null - anomaly_id: "MADD-ANO-217" diff --git a/docs/validation/sysid_fmu_claims.yaml b/docs/validation/sysid_fmu_claims.yaml index b8186593..b75ae58f 100644 --- a/docs/validation/sysid_fmu_claims.yaml +++ b/docs/validation/sysid_fmu_claims.yaml @@ -2624,7 +2624,10 @@ claims: linearisation on the candidate's side than on the other side and than the linearisation's own change, the rejection is a jump and not rounding ... The run then ends converged=False with a RuntimeWarning - naming a residual that is not differentiable"; an iteration "whose + naming a residual that is not differentiable"; "a candidate at more + than 2**5 is read a second time ... and it is a jump where that + departure is 2**8 times the other side's and the linearisation's + change over one spacing"; an iteration "whose ladder began above lam0 while the undamped step was rejected ... Where the undamped step promised, on the linear model, 2**10 times more than the short candidates lost to rounding -- or every candidate @@ -2649,7 +2652,8 @@ claims: an event, a contact or a valve whose timing depends on a trained parameter, MADD-ANO-021): converged=False and the warning where a rejected candidate within 2**-10 of a parameter shows the jump against - its mirror image (the stock TableNode -> BallNode graph, float32; that + its mirror image, over its whole length or across one float spacing + inside it (the stock TableNode -> BallNode graph, float32; that graph does not run under x64); not detected, and not claimed: a kink (a continuous residual whose slope jumps), a jump whose mirror image the bounds do not allow, a jump beside an iterate the proposal test @@ -2679,6 +2683,13 @@ claims: - tests/core/test_sysid_non_differentiable_residual.py::test_a_candidate_with_no_mirror_image_or_no_finite_residual_says_nothing - tests/core/test_sysid_non_differentiable_residual.py::test_an_iteration_whose_damped_candidates_never_moved_is_asked_again - tests/core/test_sysid_non_differentiable_residual.py::test_candidates_that_never_move_are_asked_again_once_only + - tests/core/test_sysid_non_differentiable_residual.py::test_a_small_jump_inside_a_long_candidate_is_read_across_one_spacing + - tests/core/test_sysid_non_differentiable_residual.py::test_a_jump_read_the_second_time_is_the_one_reported + - tests/core/test_sysid_non_differentiable_residual.py::test_a_small_jump_met_by_one_on_the_other_side_is_not_one_sided + - tests/core/test_sysid_non_differentiable_residual.py::test_a_smooth_residual_read_a_second_time_is_still_not_a_jump + - tests/core/test_sysid_non_differentiable_residual.py::test_a_candidate_outside_the_two_thresholds_is_read_once + - tests/core/test_sysid_non_differentiable_residual.py::test_a_second_reading_through_a_residual_that_is_not_finite_says_nothing + - tests/property/test_sysid_targeted_search.py::test_fit_lm_reads_a_small_jump_of_the_ball_from_a_long_candidate - tests/property/test_sysid_targeted_search.py::test_no_wrong_fit_among_the_fixed_draws - tests/property/test_sysid_targeted_search.py::test_fit_lm_on_the_ball_is_not_converged_where_the_loss_falls_along_its_piece status: verified diff --git a/src/maddening/sysid.py b/src/maddening/sysid.py index 5c097cf7..06ab28e2 100644 --- a/src/maddening/sysid.py +++ b/src/maddening/sysid.py @@ -3923,9 +3923,45 @@ def _check_adam_hyper(n_iter, lr, tol, betas, eps, notify_every) -> tuple[int, i #: coordinate's float spacing moves its value by 6e-4, at 21 -- and on the #: stock bouncing ball, whose bounce moves by one time step, 2.2e4 to #: 4.4e5 over 55 endings. ``2**10`` sits between, a factor of 48 above -#: the one and 21 below the other. +#: the one and 21 below the other. A wider draw (4,500 fits of the ball, +#: jaxlib 0.10.2 and 0.11.0) reads 1.9e3 to 5.8e5 at 261 jump endings and +#: 341 at one more: a small jump reads low over a whole candidate, so a +#: ratio above ``2**5`` is read a second time (:data:`_JUMP_SUSPECT`, +#: :data:`_JUMP_ACROSS`). _JUMP_EXCESS = 2.0 ** 10 +#: The ratio above which :func:`_one_sided_excess` reads a rejected +#: candidate a second time, across one spacing of the coordinates +#: (:func:`_excess_across_a_spacing`). ``2**5`` is above every ratio +#: measured on a smooth ending (:data:`_JUMP_EXCESS`: at most 21), so +#: those are read exactly as before. +_JUMP_SUSPECT = 2.0 ** 5 + +#: The second reading at which a candidate between the two thresholds +#: above is a jump. The ratio of a whole candidate understates a small +#: jump: it divides the jump by the model's change over the candidate's +#: length, and a late bounce that moves by one time step changes few +#: samples. The second reading divides by the model's change over one +#: spacing instead; the measure of the rounding is the same. Measured +#: (CPU, jaxlib 0.10.2 and 0.11.0, the same on both) on the stock ball: +#: two fits ended ``converged=True`` at ratios of 341 and 456 -- a jump of +#: 0.015 in the residual's norm, 22 float spacings of the elasticity from +#: the iterate, read from a candidate 35 spacings long -- where the second +#: reading is 1.8e3 and 2.0e3; it is 3.4e4 or more at the other 261 jump +#: endings of 4,500 drawn fits, and at most 5.2 at the six that ended on a +#: piece's own minimum. On smooth residuals stopped by the floor rule +#: (140 endings of the per-push and slow sysid tests, float32 and x64, +#: none of them asked: their first reading is under ``2**5``) it would be +#: at most 13, but for the wide ``logit`` fit of :data:`_JUMP_EXCESS` at +#: 52. ``2**8`` sits between, a factor of 4.9 above the one and 6.9 below +#: the other. +_JUMP_ACROSS = 2.0 ** 8 + +#: Halvings :func:`_excess_across_a_spacing` may take to bring a +#: candidate's step down to one spacing of the coordinates: more than the +#: bits of a float64 significand, so the cap is never what stops it. +_JUMP_BISECTIONS = 64 + def _norm64(x) -> float: """``||x||`` in float64 without squaring a value outside its range.""" @@ -3961,9 +3997,58 @@ def _gain_in_the_gap(r, J, undamped, theta, rises) -> bool: and promised > 0.0) -def _one_sided_excess(theta, r, J, rejected, mirror) -> tuple[float, float]: - """``(excess, move)``: how one-sided the residual is around ``theta``, - read from candidates a Levenberg-Marquardt iteration rejected. +def _excess_across_a_spacing(th, r64, J64, cand, cand_r, behind, mirror, dtype) -> float: + """The one-sided ratio of :func:`_one_sided_excess` for the rejected + ``cand``, with the residual's departure and the model's change taken + across one spacing of the coordinates inside the step from ``th`` + instead of over its whole length. + + The step is halved, each time keeping the half over which the residual + departs further from the linear model, until its two ends are + neighbours on the coordinates' grid. The departure between those two + is then compared with ``behind`` -- the departure at the mirror image + of the whole candidate, which :func:`_one_sided_excess` measured -- + plus the model's change over the one spacing. + + A jump the candidate crossed is between the two ends, whole: the + numerator is what it was for the whole candidate, and the model's + change, which diluted it, is a spacing's worth. Of a differentiable + residual one spacing holds only a part of the candidate's departure, + and its mirror image's departure is of the same size: the reading + is of order one however the residual curves. ``behind`` is kept from + the whole candidate, not measured again across one spacing, because + one pair of neighbours is no measure of the rounding: their residuals + often agree exactly (measured: readings of 4e3 to 1e6 at smooth + floors with the mirror image taken across one spacing too). + + ``0.0`` where a residual is not finite. + """ + lo, lo_r = th, r64 + hi, hi_r = np.asarray(cand, dtype=np.float64), np.asarray(cand_r, dtype=np.float64) + for _ in range(_JUMP_BISECTIONS): + mid, mid_r = mirror(jnp.asarray(0.5 * (lo + hi), dtype=dtype)) + mid = np.asarray(mid, dtype=np.float64) + mid_r = np.asarray(mid_r, dtype=np.float64) + if np.array_equal(mid, lo) or np.array_equal(mid, hi): + break + near = _norm64(mid_r - lo_r - J64 @ (mid - lo)) + far = _norm64(hi_r - mid_r - J64 @ (hi - mid)) + if far > near: + lo, lo_r = mid, mid_r + else: + hi, hi_r = mid, mid_r + predicted = J64 @ (hi - lo) + ahead = _norm64(hi_r - lo_r - predicted) + scale = behind + _norm64(predicted) + if not (np.isfinite(ahead) and np.isfinite(scale)) or ahead == 0.0: + return 0.0 + return float(ahead / scale) if scale > 0.0 else float(np.inf) + + +def _one_sided_excess(theta, r, J, rejected, mirror) -> tuple[float, float, bool]: + """``(excess, move, jumped)``: how one-sided the residual is around + ``theta``, read from candidates a Levenberg-Marquardt iteration + rejected, and whether that is a jump. For a rejected candidate ``theta + d`` the residual's departure from the linear model, ``||r(theta + d) - r - J d||``, is compared with the @@ -3983,6 +4068,17 @@ def _one_sided_excess(theta, r, J, rejected, mirror) -> tuple[float, float]: jump over those. No scale has to be assumed for the rounding: the mirror image measures it on the problem itself. + A ratio above :data:`_JUMP_EXCESS` is a jump. The model's change + grows with the candidate and a jump does not, so a small jump read + from a long candidate gives a ratio between the two kinds: a candidate + above :data:`_JUMP_SUSPECT` and not above :data:`_JUMP_EXCESS` is read + a second time, across one spacing of the coordinates + (:func:`_excess_across_a_spacing`), and is a jump where that reading + is above :data:`_JUMP_ACROSS`. ``excess`` and ``move`` are those of + the candidate with the largest reading among the jumps, or among all + of them where none is one; the second reading is reported only for a + candidate it made a jump. + ``rejected`` holds ``(candidate, residual, move)``; ``mirror(th)`` returns the coordinates ``th`` is evaluated at (projected onto the bounds, on the leaves' grid) and the residual there. A candidate whose @@ -3992,7 +4088,7 @@ def _one_sided_excess(theta, r, J, rejected, mirror) -> tuple[float, float]: th = np.asarray(theta, dtype=np.float64) r64 = np.asarray(r, dtype=np.float64) J64 = np.asarray(J, dtype=np.float64) - excess, at_move = 0.0, 0.0 + excess, at_move, jumped = 0.0, 0.0, False for cand, cand_r, move in rejected: d = np.asarray(cand, dtype=np.float64) - th if not d.any(): @@ -4008,9 +4104,15 @@ def _one_sided_excess(theta, r, J, rejected, mirror) -> tuple[float, float]: if not (np.isfinite(ahead) and np.isfinite(scale)) or ahead == 0.0: continue ratio = ahead / scale if scale > 0.0 else np.inf - if ratio > excess: - excess, at_move = float(ratio), float(move) - return excess, at_move + is_jump = ratio > _JUMP_EXCESS + if _JUMP_SUSPECT < ratio <= _JUMP_EXCESS: + across = _excess_across_a_spacing( + th, r64, J64, cand, cand_r, behind, mirror, theta.dtype) + if across > _JUMP_ACROSS: + ratio, is_jump = across, True + if (is_jump, ratio) > (jumped, excess): + excess, at_move, jumped = float(ratio), float(move), bool(is_jump) + return excess, at_move, jumped @@ -5802,11 +5904,21 @@ def fit_lm( linearisation's own change, the rejection is a jump and not rounding, which is alike on both sides. The run then ends ``converged=False`` with a :class:`RuntimeWarning` naming a residual that is not - differentiable; ``params`` is the lowest loss it found. Measured, in + differentiable; ``params`` is the lowest loss it found. A small jump + read from a long candidate falls short of that (the linearisation's + change grows with the step, the jump does not), so a candidate at more + than ``2**5`` is read a second time: its step is halved down to the + one float spacing of the coordinates across which the residual departs + most (a residual evaluation a halving), and it is a jump where that + departure is ``2**8`` times the other side's and the linearisation's + change over one spacing. Measured, in float32 and under x64: at most 21 on fits of smooth residuals stopped by the floor rule (the spring over its transforms, units and noise levels; ``HeartPumpNode`` at noise 0.02 and 1), and 2.2e4 to 4.4e5 on - the ball (float32). Not detected: a *kink* -- a continuous residual whose slope + the ball (float32); the second reading, at most 52 on the former and + 1.8e3 or more on the ball's jumps, two of which read 341 and 456 the + first time. + Not detected: a *kink* -- a continuous residual whose slope jumps -- in whose crease a run ends ``converged=True``, at a point no step lowers the loss from but whose ``J`` is one-sided; a jump whose mirror image the bounds do not allow; and a jump beside an iterate the @@ -6312,7 +6424,7 @@ def _gauss_newton_stationary(th, before, r, J, loss_at_th=None) -> bool: # shorter only because it is damped more, so it is not tested. accepted = False stationary = False - jump = (0.0, 0.0) + jump = (0.0, 0.0, False) verdict["representable"] = True # A loss of exactly 0.0 from a residual that is not: the framed # sum's unframing underflowed float64 (a float64 residual below @@ -6467,10 +6579,10 @@ def _gauss_newton_stationary(th, before, r, J, loss_at_th=None) -> bool: # loss that can be far above the minimum. The rejected # candidates tell the two apart (:func:`_one_sided_excess`). jump = _one_sided_excess(theta, r, J, rejected, _mirrored) - stationary = not jump[0] > _JUMP_EXCESS + stationary = not jump[2] if again: continue - if jump[0] > _JUMP_EXCESS: + if jump[2]: warnings.warn(_JUMP_WARNING.format(move=jump[1], excess=jump[0]), RuntimeWarning, stacklevel=2) break diff --git a/tests/core/test_sysid_non_differentiable_residual.py b/tests/core/test_sysid_non_differentiable_residual.py index f397fcff..449a2e25 100644 --- a/tests/core/test_sysid_non_differentiable_residual.py +++ b/tests/core/test_sysid_non_differentiable_residual.py @@ -219,17 +219,20 @@ def watched(*args): for k in ("stiffness", "damping")]) assert np.ptp(np.asarray(answers), axis=0).max() < 1e-3 assert read, "premise: at least one of the fits was stopped by the floor rule" - # Of order one; the threshold is 2**10. - assert max(read) < sysid._JUMP_EXCESS / 2.0 ** 4 # noqa: SLF001 + # Of order one; the threshold is 2**10, and nothing is read a second + # time below 2**5. + assert max(read) < sysid._JUMP_SUSPECT # noqa: SLF001 # --------------------------------------------------------------------------- # The one-sided test itself # --------------------------------------------------------------------------- -def _excess(fn, theta, steps, lo=-np.inf, hi=np.inf): +def _reading(fn, theta, steps, lo=-np.inf, hi=np.inf, asked=None): """``_one_sided_excess`` of ``fn`` at ``theta`` for candidates - ``theta + step``, with its Jacobian by ``jacfwd``.""" + ``theta + step``, with its Jacobian by ``jacfwd``: ``(excess, move, + jumped)``. ``asked`` collects the points it evaluated (mirror images + and halvings).""" theta = jnp.asarray(theta, jnp.float64 if jax.config.jax_enable_x64 else jnp.float32) r, J = fn(theta), jax.jacfwd(fn)(theta) rejected = [] @@ -239,11 +242,21 @@ def _excess(fn, theta, steps, lo=-np.inf, hi=np.inf): def mirror(th): th = jnp.clip(th, lo, hi) + if asked is not None: + asked.append(np.asarray(th)) return th, fn(th) return sysid._one_sided_excess(theta, r, J, rejected, mirror) # noqa: SLF001 +def _excess(fn, theta, steps, **kwargs): + """``(excess, move)`` of :func:`_reading`, the verdict being the + threshold's: a jump exactly where ``excess`` is over ``2**10``.""" + excess, move, jumped = _reading(fn, theta, steps, **kwargs) + assert jumped == (excess > sysid._JUMP_EXCESS) # noqa: SLF001 + return excess, move + + def _smooth(th): return jnp.stack([jnp.sin(th[0]) + th[1] ** 2, jnp.exp(th[0] * th[1]), th[0] - th[1]]) @@ -285,6 +298,122 @@ def to_nan(th): assert _excess(to_nan, [1.0, 0.5], [[-1e-4, 0.0]])[0] == 0.0 +# --------------------------------------------------------------------------- +# A small jump read from a long candidate: the second reading +# --------------------------------------------------------------------------- + +#: A step of 0.02 in the second entry, read from a candidate of 2e-4: the +#: jump over the linear change of the whole candidate (2.8e-4) is about 70, +#: between the two thresholds of the first reading. +_SMALL_JUMP, _LONG_STEP = 0.02, [2e-4, 0.0] + + +def _small_jump(at): + def fn(th): + return _smooth(th) + (jnp.where(th[0] > at, _SMALL_JUMP, 0.0) + * jnp.asarray([0.0, 1.0, 0.0])) + return fn + + +def _steep(th): + # Smooth, with a slope that grows by exp(6.5) over a step of 2e-4. + return jnp.stack([jnp.exp(3.25e4 * (th[0] - 1.0)), th[1]]) + + +def _first_reading(monkeypatch, fn, steps, **kwargs): + """The reading with no second one: ``(excess, move, jumped)``, and the + premise that the largest ratio is between the two thresholds.""" + with monkeypatch.context() as patch: + patch.setattr(sysid, "_JUMP_SUSPECT", np.inf) + whole = _reading(fn, [1.0, 0.5], steps, **kwargs) + assert sysid._JUMP_EXCESS > whole[0] > sysid._JUMP_SUSPECT # noqa: SLF001 + assert not whole[2] + return whole + + +# The edge beside the iterate, and inside the candidate's step (the fit of +# the ball that showed this ended 22 float spacings from its edge). +@pytest.mark.parametrize("at", [1.0, 1.00005, 1.00019]) +def test_a_small_jump_inside_a_long_candidate_is_read_across_one_spacing(at, monkeypatch): + whole = _first_reading(monkeypatch, _small_jump(at), [_LONG_STEP]) + excess, move, jumped = _reading(_small_jump(at), [1.0, 0.5], [_LONG_STEP]) + assert jumped and move == pytest.approx(2e-4) + # The jump over the mirror image's departure and one float32 spacing's + # linear change: tens of thousands. + assert excess > 2.0 ** 4 * sysid._JUMP_ACROSS > whole[0] # noqa: SLF001 + + +def test_a_jump_read_the_second_time_is_the_one_reported(monkeypatch): + """Beside a candidate whose first reading is larger and which is no + jump, whichever of the two comes first.""" + def fn(th): + # Along the first coordinate the small jump, met on the other side + # by one 400 times smaller; along the second a smooth residual + # whose slope grows by exp(9) over a step of 1e-4. + met = jnp.where(2.0 - th[0] > 1.00005, _SMALL_JUMP / 400.0, 0.0) + steep = jnp.exp(9e4 * (th[1] - 0.5)) + return _small_jump(1.00005)(th) + jnp.stack([met, 0.0 * met, steep]) + + jump, smooth = _LONG_STEP, [0.0, 1e-4] + only_jump = _reading(fn, [1.0, 0.5], [jump]) + only_smooth = _reading(fn, [1.0, 0.5], [smooth]) + assert only_jump[2] and not only_smooth[2] + # Premise: the smooth candidate's reading is the larger number. + assert sysid._JUMP_EXCESS > only_smooth[0] > only_jump[0] > sysid._JUMP_ACROSS # noqa: SLF001 + assert only_jump[1] == pytest.approx(2e-4) and only_smooth[1] == pytest.approx(1e-4) + for steps in ([jump, smooth], [smooth, jump]): + assert _reading(fn, [1.0, 0.5], steps) == only_jump + + +def test_a_small_jump_met_by_one_on_the_other_side_is_not_one_sided(monkeypatch): + """A second, smaller jump where the mirror image of the candidate + falls: the second reading is the ratio of the two jumps, under its + threshold, and the candidate keeps its first.""" + def fn(th): + return _small_jump(1.00005)(th) + (jnp.where(2.0 - th[0] > 1.00005, 3e-4, 0.0) + * jnp.asarray([1.0, 0.0, 0.0])) + + asked = [] + reading = _reading(fn, [1.0, 0.5], [_LONG_STEP], asked=asked) + assert len(asked) > 4, "premise: the candidate was read a second time" + assert reading == _first_reading(monkeypatch, fn, [_LONG_STEP]) + + +def test_a_smooth_residual_read_a_second_time_is_still_not_a_jump(monkeypatch): + """A residual curved enough over the candidate to be read again: one + spacing holds a small part of the departure, and the mirror image's is + of its size.""" + asked, across = [], [] + real = sysid._excess_across_a_spacing # noqa: SLF001 + monkeypatch.setattr(sysid, "_excess_across_a_spacing", + lambda *args: across.append(real(*args)) or across[-1]) + reading = _reading(_steep, [1.0, 0.5], [_LONG_STEP], asked=asked) + assert len(asked) > 4 and len(across) == 1 + assert 0.0 < across[0] < 1.0 + assert reading == _first_reading(monkeypatch, _steep, [_LONG_STEP]) + + +def test_a_candidate_outside_the_two_thresholds_is_read_once(): + """A ratio of order one costs the mirror image and nothing more, and + one over the upper threshold needs no second reading.""" + for fn, steps in ((_smooth, [[1e-4, 0.0], [-2e-4, 1e-4], [3e-5, 3e-5]]), + (_jumping, [[1e-4, 0.0]])): + asked = [] + _excess(fn, [1.0, 0.5], steps, asked=asked) + assert len(asked) == len(steps) + + +def test_a_second_reading_through_a_residual_that_is_not_finite_says_nothing(monkeypatch): + def fn(th): + gap = (th[0] > 1.00004) & (th[0] < 1.00016) + return _small_jump(1.00005)(th) + jnp.where(gap, jnp.nan, 0.0) + + asked = [] + reading = _reading(fn, [1.0, 0.5], [_LONG_STEP], asked=asked) + assert len(asked) > 4, "premise: the candidate was read a second time" + assert reading == _first_reading(monkeypatch, fn, [_LONG_STEP]) + + # --------------------------------------------------------------------------- # A ladder of candidates that never moved is not a ladder of rejections # --------------------------------------------------------------------------- diff --git a/tests/property/test_sysid_targeted_search.py b/tests/property/test_sysid_targeted_search.py index 0d1a4602..6715db19 100644 --- a/tests/property/test_sysid_targeted_search.py +++ b/tests/property/test_sysid_targeted_search.py @@ -615,24 +615,31 @@ def residual(p): return _BALL["problem"] +def _fit_of_the_ball(case: BallCase, key="ball"): + """``fit_lm`` on ``case`` (inside ``precision(False)``): the problem, + the result and the messages of the warnings it raised.""" + problem = _ball_problem() + start = jax.tree.map(lambda x: x, problem.params) + + def put(node, leaf, value): + start["nodes"][node][leaf] = jnp.asarray(value, start["nodes"][node][leaf].dtype) + + put("record", "elasticity", case.elasticity_true) + put("record", "gravity", -9.81 * case.gravity_true) + put("ball", "elasticity", case.elasticity_start) + put("ball", "gravity", -9.81 * case.gravity_start) + with _shared_programs(key), warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + res = fit_lm(problem.gm, problem.residual, params=start) + return problem, res, [str(w.message) for w in caught] + + def converged_beside_a_lower_loss(case: BallCase): """(d) ``converged=True`` on a residual with jumps is still a point no nearby one lowers the loss from by more than rounding explains.""" key = "ball" with precision(False): - problem = _ball_problem() - start = jax.tree.map(lambda x: x, problem.params) - - def put(node, leaf, value): - start["nodes"][node][leaf] = jnp.asarray(value, start["nodes"][node][leaf].dtype) - - put("record", "elasticity", case.elasticity_true) - put("record", "gravity", -9.81 * case.gravity_true) - put("ball", "elasticity", case.elasticity_start) - put("ball", "gravity", -9.81 * case.gravity_start) - with _shared_programs(key), warnings.catch_warnings(): - warnings.simplefilter("ignore") - res = fit_lm(problem.gm, problem.residual, params=start) + problem, res, _ = _fit_of_the_ball(case, key) details = (bool(res.converged), float(res.best_loss), int(res.n_iter)) if not res.converged: return 0.0, details @@ -700,10 +707,21 @@ def test_no_wrong_fit_found_by_the_search(name): #: too short to gain more than the loss's rounding, with the undamped step #: across the next jump; the floor rule now asks such an iterate again from #: the starting damping (a score of 30 before; ``converged=True`` at a loss -#: of 0.021). +#: of 0.021). The last two were found by the slow hunt on jaxlib 0.10.2 +#: (the same on 0.11.0) and among 4,500 uniform draws: a *small* jump, 22 +#: float spacings of the elasticity from the iterate, read from a rejected +#: candidate 35 spacings long, where the jump over the whole candidate's +#: linear change is 456 and 341 -- under the ``2**10`` that reads a jump; +#: such a candidate is now read a second time, across one spacing (scores +#: of 11.8 and 5.6 before; ``converged=True`` at a loss of 0.006). +_READ_ACROSS_ONE_SPACING = [ + (0.6133024381551229, 0.7699776319888987, 0.5, 1.0), + (0.5903806445545127, 0.7114224216650756, 0.719718582123716, 1.1640531624602226), +] _CONVERGED_WHERE_THE_LOSS_FELL = [ (0.7, 1.0, 0.625, 1.333521432163324), (0.7, 1.2429505644303036, 0.625, 1.0), + *_READ_ACROSS_ONE_SPACING, ] @@ -713,6 +731,44 @@ def test_fit_lm_on_the_ball_is_not_converged_where_the_loss_falls_along_its_piec assert score <= BESIDE, (score, details) +@pytest.mark.parametrize("cell", _READ_ACROSS_ONE_SPACING) +def test_fit_lm_reads_a_small_jump_of_the_ball_from_a_long_candidate(cell, monkeypatch): + """What the two endings are, not only that the score passes them: the + whole candidate reads between the two thresholds, the second reading + across one spacing is over its own, and the run ends + ``converged=False`` with the warning of a residual that is not + differentiable. Without the second reading it is ``converged=True`` + where the loss still falls along the piece.""" + real, second = sysid._excess_across_a_spacing, [] # noqa: SLF001 + + def watched(*args): + second.append(real(*args)) + return second[-1] + + monkeypatch.setattr(sysid, "_excess_across_a_spacing", watched) + with precision(False): + _, res, messages = _fit_of_the_ball(BallCase(*cell)) + assert not res.converged + assert sum("not differentiable" in text for text in messages) == 1 + assert len(second) == 1 and second[0] > 2.0 ** 2 * sysid._JUMP_ACROSS # noqa: SLF001 + assert f"{second[0]:.1e} times further" in "".join(messages) + # The second reading is put to its own threshold: one between that and + # the first reading's ends the run the same way. + between = 2.0 * sysid._JUMP_ACROSS # noqa: SLF001 + assert between < sysid._JUMP_EXCESS # noqa: SLF001 + monkeypatch.setattr(sysid, "_excess_across_a_spacing", lambda *args: between) + with precision(False): + _, capped, messages = _fit_of_the_ball(BallCase(*cell)) + assert not capped.converged and capped.best_loss == res.best_loss + assert sum("not differentiable" in text for text in messages) == 1 + assert f"{between:.1e} times further" in "".join(messages) + # Premise: the reading of the whole candidate alone leaves it converged, + # at the same point, with the score over its threshold. + monkeypatch.setattr(sysid, "_JUMP_SUSPECT", np.inf) + score, details = converged_beside_a_lower_loss(BallCase(*cell)) + assert details[0] and details[1] == res.best_loss and score > BESIDE, (score, details) + + #: Fits that stopped with iterations left, ``converged=False``, the damping #: on the upper end of a wide range (a ``logit`` one with the edge warning) #: and the stiffness two to three times its truth -- at no optimum: the From 8879b311a4319ba3ae15f612f0be7e97a148422f Mon Sep 17 00:00:00 2001 From: Nicholas Ehsan Roy Date: Wed, 7 Oct 2026 12:07:30 +0200 Subject: [PATCH 2/3] test(sysid): the spec strategy draws no logit bound a float32 leaf cannot have spec_and_value drew ParamSpec(bounds=(1.0027701900272462e-38, 1.0), transform='logit'): a subnormal float32 lower bound, which ParamSpec refuses by name for a float32 leaf. The strategy predates that refusal, so two properties of TestBoundsAndTransforms failed at random. A subnormal logit bound is now moved to 0.0 (the bound the message recommends) rather than rejected; the drawn spec joins the refusal's own test, and the accepted side of the boundary is an explicit example. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkCde7g23gTziUjYAvnKD --- ..._inside_bounds_near_the_smallest_normal.py | 11 ++++ tests/property/test_sysid_contract.py | 57 ++++++++++++++++++- 2 files changed, 67 insertions(+), 1 deletion(-) diff --git a/tests/core/test_constrain_lands_inside_bounds_near_the_smallest_normal.py b/tests/core/test_constrain_lands_inside_bounds_near_the_smallest_normal.py index 5f6c796b..5284d2df 100644 --- a/tests/core/test_constrain_lands_inside_bounds_near_the_smallest_normal.py +++ b/tests/core/test_constrain_lands_inside_bounds_near_the_smallest_normal.py @@ -154,6 +154,17 @@ def test_a_logit_interval_narrower_than_four_smallest_normals_is_refused_by_name spec.to_constrained(jnp.asarray(F(0.0))) with pytest.raises(ValueError, match="cannot map .* is a subnormal"): spec.to_unconstrained(jnp.asarray(F(0.5))) + # The same shape as a property test drew it (float32: 0.85 of the + # smallest normal). Under x64 it is an ordinary bound, and mapped. + drawn = ParamSpec(bounds=(1.0027701900272462e-38, 1.0), transform="logit") + if np.dtype(dtype) == np.float32: + for refused in (lambda: drawn.to_constrained(jnp.asarray(F(0.0))), + lambda: drawn.to_unconstrained(jnp.asarray(F(0.5))), + lambda: drawn.check(jnp.asarray(F(0.5)))): + with pytest.raises(ValueError, match="cannot map .* is a subnormal"): + refused() + else: + drawn.check(drawn.to_constrained(jnp.asarray(F(0.0)))) for lo, hi in ((0.0, tiny), (0.0, 3.5 * tiny), (-tiny, tiny), (tiny, 2 * tiny)): spec = ParamSpec(bounds=(lo, hi), transform="logit") with pytest.raises(ValueError, match="cannot map .* width .* smallest normals"): diff --git a/tests/property/test_sysid_contract.py b/tests/property/test_sysid_contract.py index 75752cc3..cdf0e9cf 100644 --- a/tests/property/test_sysid_contract.py +++ b/tests/property/test_sysid_contract.py @@ -641,9 +641,56 @@ def body(s, _): # --------------------------------------------------------------------------- +#: The smallest normal float32: ``spec_and_value`` values are float32 leaves. +_F32_TINY = float(np.finfo(np.float32).tiny) +#: Drawn at random, and failing two tests of ``TestBoundsAndTransforms``, in +#: local runs: a ``logit`` lower bound that is a subnormal float32. +_SUBNORMAL_LOGIT_LOWER_BOUND = 1.0027701900272462e-38 + + +def _logit_bound_a_float32_leaf_can_have(bound: float) -> float: + """``bound``, or ``0.0`` where it is a subnormal float32 number. + + A ``logit`` spec with such a bound cannot map a float32 leaf and says + so by name (``ParamSpec._require_representable``: the transform's + arithmetic reads a subnormal operand as zero, so every value would be + measured from another bound than the one declared). That refusal is + held by ``tests/core/test_constrain_lands_inside_bounds_near_the_smallest_normal.py::test_a_logit_interval_narrower_than_four_smallest_normals_is_refused_by_name``; + the properties here are about specs the library maps. A draw of one + is moved to the bound the message recommends rather than rejected: + nothing is thrown away, and the neighbour is a spec that has to work. + Bounds no float32 holds exactly, and subnormal bounds of the other + transforms, are left as drawn: those are legal. + """ + with np.errstate(under="ignore"): + rounded = abs(float(np.float32(bound))) + return 0.0 if 0.0 < rounded < _F32_TINY else bound + + +def test_the_spec_strategy_moves_only_a_subnormal_logit_bound(): + move = _logit_bound_a_float32_leaf_can_have + assert _SUBNORMAL_LOGIT_LOWER_BOUND < _F32_TINY + for subnormal in (_SUBNORMAL_LOGIT_LOWER_BOUND, -_SUBNORMAL_LOGIT_LOWER_BOUND, + 1.4e-45, -1e-40, float(np.nextafter(np.float32(_F32_TINY), np.float32(0)))): + assert move(subnormal) == 0.0 + # Zero, the smallest normal either side, a number below every float32 + # (it rounds to zero, and the library reads it so), ordinary bounds and + # one no float32 holds exactly: as drawn. + for kept in (0.0, _F32_TINY, -_F32_TINY, 1e-60, 1e-30, -20.0, 1.6977594293117515): + assert move(kept) == kept + # What it guards: the drawn spec is refused for a float32 leaf, and its + # neighbour is mapped. + with pytest.raises(ValueError, match="cannot map .* is a subnormal"): + ParamSpec(bounds=(_SUBNORMAL_LOGIT_LOWER_BOUND, 1.0), + transform="logit").to_constrained(jnp.float32(0.0)) + spec = ParamSpec(bounds=(move(_SUBNORMAL_LOGIT_LOWER_BOUND), 1.0), transform="logit") + spec.check(spec.to_constrained(jnp.float32(0.0))) + + @st.composite def spec_and_value(draw): - """A constructible ``ParamSpec`` and a value strictly inside its bounds. + """A constructible ``ParamSpec`` that maps a float32 leaf, and a value + strictly inside its bounds. The intervals are kept wide (at least 1e-2, and wide relative to their endpoints) because ``to_constrained`` documents that an @@ -656,6 +703,9 @@ def spec_and_value(draw): lo = draw(_finite(-20.0, 20.0)) span = draw(_finite(1e-2, 40.0)) if transform == "logit": + # ``lo + span`` is never subnormal: ``span`` is at least 1e-2, and + # a float64 sum of that size is zero or far above 1e-38. + lo = _logit_bound_a_float32_leaf_can_have(lo) bounds = (lo, lo + span) value = lo + span * draw(_finite(0.05, 0.95)) elif transform == "log": @@ -705,6 +755,10 @@ def test_a_transform_refuses_the_bounds_it_cannot_work_with(self, data): ParamSpec(transform="softplus") @given(spec_value=spec_and_value()) + # The two sides of the bound the strategy moves: zero, and the + # smallest normal float32 (a subnormal one is refused by name). + @example(spec_value=(ParamSpec(bounds=(0.0, 1.0), transform="logit"), 0.5)) + @example(spec_value=(ParamSpec(bounds=(_F32_TINY, 1.0), transform="logit"), 0.5)) @settings(max_examples=EXAMPLES_CHEAP, deadline=None) def test_constrain_inverts_unconstrain_inside_the_bounds(self, spec_value): spec, value = spec_value @@ -726,6 +780,7 @@ def test_constrain_inverts_unconstrain_inside_the_bounds(self, spec_value): @example(spec_value=(ParamSpec(bounds=(1.4e-45, None), transform=None), 1.0), u=-3.0) @example(spec_value=(ParamSpec(bounds=(-1e-40, 1e-40), transform=None), 0.0), u=7.0) @example(spec_value=(ParamSpec(bounds=(0.0, 7.5e-37), transform="logit"), 1e-37), u=-800.0) + @example(spec_value=(ParamSpec(bounds=(_F32_TINY, 1.0), transform="logit"), 0.5), u=-800.0) @settings(max_examples=EXAMPLES_CHEAP, deadline=None) def test_constrain_lands_inside_the_bounds_from_any_coordinate( self, spec_value, u, From 09e10a11e950b94e2cc987971b0273a81c5e9319 Mon Sep 17 00:00:00 2001 From: Nicholas Ehsan Roy Date: Wed, 7 Oct 2026 12:52:20 +0200 Subject: [PATCH 3/3] test(sysid): the ball pin asserts the verdict, not how many candidates were read twice On jaxlib 0.11.2 the ball's trajectory rounds differently and the ending of the first cell has two candidates between the thresholds (second readings 4.4e3 and 2.1e3) where 0.10.2 and 0.11.0 have one (2.0e3). The test counted them. It now asserts what the fix guarantees on every version: converged=False with one warning; the first reading alone finds no jump and is in the band read again; the verdict is a second reading over its threshold, and the warning prints it. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkCde7g23gTziUjYAvnKD --- tests/property/test_sysid_targeted_search.py | 47 ++++++++++++++------ 1 file changed, 34 insertions(+), 13 deletions(-) diff --git a/tests/property/test_sysid_targeted_search.py b/tests/property/test_sysid_targeted_search.py index 6715db19..13ac77e1 100644 --- a/tests/property/test_sysid_targeted_search.py +++ b/tests/property/test_sysid_targeted_search.py @@ -733,25 +733,46 @@ def test_fit_lm_on_the_ball_is_not_converged_where_the_loss_falls_along_its_piec @pytest.mark.parametrize("cell", _READ_ACROSS_ONE_SPACING) def test_fit_lm_reads_a_small_jump_of_the_ball_from_a_long_candidate(cell, monkeypatch): - """What the two endings are, not only that the score passes them: the - whole candidate reads between the two thresholds, the second reading - across one spacing is over its own, and the run ends - ``converged=False`` with the warning of a residual that is not - differentiable. Without the second reading it is ``converged=True`` - where the loss still falls along the piece.""" - real, second = sysid._excess_across_a_spacing, [] # noqa: SLF001 - - def watched(*args): - second.append(real(*args)) + """What the two endings are, not only that the score passes them: no + whole candidate reads as a jump, the largest reads between the two + thresholds, a second reading across one spacing is over its own, and + the run ends ``converged=False`` with the warning of a residual that + is not differentiable. Without the second reading it is + ``converged=True`` where the loss still falls along the piece. + + How many candidates are read a second time, and their readings, are + the trajectory's rounding and differ between jaxlib versions (one at + 2.0e3 on 0.10.2 and 0.11.0, two at 4.4e3 and 2.1e3 on 0.11.2 for the + first cell), so neither is asserted.""" + across, excess = sysid._excess_across_a_spacing, sysid._one_sided_excess # noqa: SLF001 + second, verdicts, first = [], [], [] + + def watched_across(*args): + second.append(across(*args)) return second[-1] - monkeypatch.setattr(sysid, "_excess_across_a_spacing", watched) + def watched_excess(*args): + verdicts.append(excess(*args)) + with monkeypatch.context() as patch: + patch.setattr(sysid, "_JUMP_SUSPECT", np.inf) + first.append(excess(*args)) + return verdicts[-1] + + monkeypatch.setattr(sysid, "_excess_across_a_spacing", watched_across) + monkeypatch.setattr(sysid, "_one_sided_excess", watched_excess) with precision(False): _, res, messages = _fit_of_the_ball(BallCase(*cell)) assert not res.converged assert sum("not differentiable" in text for text in messages) == 1 - assert len(second) == 1 and second[0] > 2.0 ** 2 * sysid._JUMP_ACROSS # noqa: SLF001 - assert f"{second[0]:.1e} times further" in "".join(messages) + reading, _, jumped = verdicts[-1] + whole, _, whole_jumped = first[-1] + # The first reading alone finds no jump, and is in the band that is + # read again; the verdict is a second reading, over its threshold. + assert not whole_jumped + assert sysid._JUMP_SUSPECT < whole <= sysid._JUMP_EXCESS # noqa: SLF001 + assert jumped and reading in second and reading > sysid._JUMP_ACROSS # noqa: SLF001 + assert f"{reading:.1e} times further" in "".join(messages) + monkeypatch.setattr(sysid, "_one_sided_excess", excess) # The second reading is put to its own threshold: one between that and # the first reading's ends the run the same way. between = 2.0 * sysid._JUMP_ACROSS # noqa: SLF001