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_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/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_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, diff --git a/tests/property/test_sysid_targeted_search.py b/tests/property/test_sysid_targeted_search.py index 0d1a4602..13ac77e1 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,65 @@ 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: 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] + + 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 + 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 + 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