diff --git a/CHANGELOG.md b/CHANGELOG.md index b4a109cb0..96bd83f55 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -197,6 +197,54 @@ in `CMakeLists.txt`) is derived from it. ulp. SED-ML export refuses a per-species `atol` outright rather than writing one entry as `KISAO:0000211`, which would describe a different run. +### Changed +- **The GH #88 periodic step bound is now derived only where the discontinuity + roots provably cannot reach (issue #274).** #262 measured that no corpus model's + answer the bound changes; this acts on it. The bound and the GH #72/#231/#259 + roots address the same hazard — an integrator step spanning a schedule edge — + and where a root already forces the stop, `max_step` is pure cost: it shortens + every step over the whole horizon, not just the ones near an edge. + + The test is value-position reachability from the ODE RHS through the assignment + rules, pruning every piecewise condition that IS a registered root + (`_periodic_disc_escapes_roots`). A `floor`/`ceiling`/`modulo` influences the + RHS two ways and only one of them is rooted: through a piecewise **condition** + (`exposure = piecewise(D, frac < w, 0)` — a root on that condition forces the + stop, so the bound adds nothing), or as a **value** (`dose = D * frac`, or + `k = k0 * i` with `i = floor(time/24)` — the sawtooth jump is an RHS + discontinuity in its own right, at an instant no relational brackets). Only the + second keeps the bound. + + On the 25 corpus models that carried one: kept on 10, dropped on 15, cutting + 1,015,702 internal steps to 519,472 (**-49%**) with no answer moving — + `MODEL0406553884` -60%, `MODEL0406793751` -44%, `MODEL0847869198` -40%, + `MODEL1708310001` -29%. That last one is the check that matters: it is the model + #259 gave the roots to, with an exact segmented oracle (953.07) confirming they + resolve its schedule, and the predicate drops its bound without being told. + + **Not a root count.** The tempting predicate — "apply the bound only when the + model registered no roots" — is unsound, and measurably so: three corpus models + (`BIOMD0000000577`, `BIOMD0000000589`, `MODEL1006230027`) are rooted yet carry a + periodic disc node reaching the RHS ungated. `BIOMD0000000589` is the clearest, + with 15 of its 16 schedule conditions of the form `i == 0`, `i == 1`, … on + `i := floor(time/24)`: **equalities**, which the inequality-only emitters cannot + root, on a step function that genuinely moves with time. + + A new fixture makes that concrete, and is the first thing in the suite to + witness the bound's *necessity* since #259 removed the old one — the reason this + is a narrowing rather than a retirement. Its pulse is gated by an equality on a + floor while an unrelated `time < 5` threshold is rooted, so the model has roots + and still needs the bound: with it, `y(10) = 36.788` against an exact 36.788; + without, 182.212 — 395% high, identically at rtol 1e-9 and 1e-11. Tol-stably + wrong is the signature of a pulse never sampled at all. A root-count predicate + returns 182.212 here. + + `Model._periodic_disc_max_step` is `None` for the 15 dropped models, so anything + reading it now reads "the bound this model needs" rather than "a periodic + schedule was detected". `_periodic_time_disc_max_step` keeps its pre-#274 + behaviour when called without the registered-condition set, so a caller that has + not run the root scan cannot accidentally drop a bound it has no basis to drop. + ### Fixed - **A parameter that a same-named function shadows was still listed as a knob (issues #256 and #266 — the same defect reported from the SBML side and the diff --git a/python/bngsim/_sbml_loader.py b/python/bngsim/_sbml_loader.py index 41507ef0c..037f1b5c9 100644 --- a/python/bngsim/_sbml_loader.py +++ b/python/bngsim/_sbml_loader.py @@ -1809,7 +1809,113 @@ def _replace_exprtk_symbols(expr: str, replacements: dict[str, str]) -> str: return pat.sub(lambda m: f"({replacements[m.group(1)]})", expr) -def _periodic_time_disc_max_step(sbml_model, func_defs, base_ctx): +def _piecewise_branches(node) -> tuple[list, list]: + """``(condition_children, value_children)`` of a MathML ``piecewise``. + + Each ```` contributes a value at an even index and its condition at + the next odd one; a trailing ```` is a value with no condition. + """ + nc = node.getNumChildren() + conds, vals = [], [] + i = 0 + while i < nc: + if i + 1 < nc: + vals.append(node.getChild(i)) + conds.append(node.getChild(i + 1)) + i += 2 + else: + vals.append(node.getChild(i)) + i += 1 + return conds, vals + + +def _periodic_disc_escapes_roots( + rhs_maths: list, + asg_math: dict, + td: set, + func_defs: dict, + registered_conditions: set, +) -> bool: + """Does a periodic discontinuity reach the RHS by a path no root stops at? + + GH #274. The GH #88 step bound and the GH #72/#231/#259 roots address the + same hazard — an integrator step spanning a schedule edge — by different + means, and where a root already forces the stop the bound is pure cost + (``max_step`` shortens every step over the whole horizon, not just the ones + near an edge: up to 2.5x the internal steps on the corpus, for an answer + that does not move). So the bound is kept only for what the roots do not + reach, which is what this decides. + + A time-dependent ``floor``/``ceiling``/``modulo`` influences the RHS two + ways, and only one of them is rooted: + + * through a piecewise CONDITION — ``exposure = piecewise(D, frac < w, 0)``. + A registered root on that condition forces a stop exactly at the edge, so + the branch cannot be stepped over and the bound adds nothing. That path is + pruned here. + * as a VALUE — ``dose = D * frac``, or ``k = k0 * i`` with + ``i = floor(time/24)``. The sawtooth jump is then an RHS discontinuity in + its own right, at an instant no relational brackets. Rooting some other + condition in the same model does not help, and the bound is what keeps a + step from spanning it. + + So this is value-position reachability from the ODE RHS, followed through + the assignment rules and pruned at every rooted condition; ``True`` means + some periodic disc node survives that pruning. + + A note for anyone extending this: it is deliberately ONE pass that carries + context down, rather than two passes compared afterwards. libsbml returns a + *fresh* Python wrapper from every ``getChild()`` call, so ``id(node)`` is not + stable across traversals and any "was this node also seen over there" test + silently answers no — which reads as full root coverage and drops the bound + everywhere. + """ + disc_types = { + libsbml.AST_FUNCTION_FLOOR, + libsbml.AST_FUNCTION_CEILING, + libsbml.AST_FUNCTION_REM, + } + + def _scan(node, names_out: list) -> bool: + """True as soon as an unrooted time-dependent disc node is reached.""" + stack = [node] + while stack: + n = stack.pop() + ntype = n.getType() + if ntype in disc_types and (_ast_references_time(n) or (_ast_name_set(n) & td)): + return True + if ntype == libsbml.AST_FUNCTION_PIECEWISE: + conds, vals = _piecewise_branches(n) + stack.extend(vals) + for cond in conds: + emitted: list[str] = [] + _collect_time_discontinuity_conditions(cond, func_defs, emitted, None, td) + # Emitted AND registered ⇒ the integrator stops at this edge. + if emitted and all(e in registered_conditions for e in emitted): + continue + stack.append(cond) + continue + if ntype == libsbml.AST_NAME: + names_out.append(n.getName()) + stack.extend(n.getChild(k) for k in range(n.getNumChildren())) + return False + + pending: list[str] = [] + for math in rhs_maths: + if _scan(math, pending): + return True + seen: set[str] = set() + while pending: + name = pending.pop() + if name in seen: + continue + seen.add(name) + if name in asg_math and _scan(asg_math[name], pending): + return True + return False + + +def _periodic_time_disc_max_step(sbml_model, func_defs, base_ctx, registered_conditions=None): """Recommend a ``max_step_size`` that resolves periodic floor/modulo dosing discontinuities in the ODE RHS, or ``None`` when the model has none. @@ -1820,6 +1926,14 @@ def _periodic_time_disc_max_step(sbml_model, func_defs, base_ctx): signature — the integer value of each floor/ceil/modulo and the active branch of each piecewise — which changes ONLY at a true edge (continuous values like `rem_time` move every sample and would mask the spacing). + + ``registered_conditions`` is the set of discontinuity conditions the caller + just registered as CVODE roots. Given it, GH #274 returns ``None`` for a + schedule those roots already cover, since a root forces the same stop the + bound would and costs nothing between edges — see + :func:`_periodic_disc_escapes_roots`. Omitting it keeps the pre-#274 + behaviour: a caller that has not run the root scan has no basis for that + judgement, and the safe default is to bound. """ global _DISC_FUNCS if _DISC_FUNCS is None: @@ -1893,6 +2007,17 @@ def _periodic_time_disc_max_step(sbml_model, func_defs, base_ctx): if not has_periodic: return None + # GH #274 — the roots have overtaken the bound on most of the shapes it was + # written for. Derive one only where they provably cannot reach; see + # :func:`_periodic_disc_escapes_roots`. Checked before the numeric edge scan + # below, which is the expensive part and is pure waste when the answer is + # None. ``registered_conditions=None`` means the caller did not run the root + # scan and cannot make this judgement, so the bound stands. + if registered_conditions is not None and not _periodic_disc_escapes_roots( + rhs_maths, asg_math, td, func_defs, registered_conditions + ): + return None + # Order the time-dependent assignment rules by dependency so each can be # evaluated as a pure function of t into a fresh ctx. td_rules: list[tuple[str, libsbml.ASTNode]] = [] @@ -6496,11 +6621,16 @@ def _unalias(expr_str: str) -> str: seen_disc.add(cond) builder.add_discontinuity_trigger(cond) - # Periodic floor()/modulo dosing schedules (GH #88) can't be resolved by the - # monotonic roots above — register a step-size bound that keeps the adaptive + # Periodic floor()/modulo dosing schedules (GH #88) that the monotonic roots + # above can't resolve — register a step-size bound that keeps the adaptive # integrator from stepping over a narrow dose pulse. None for any model with - # no time-dependent floor/ceil/modulo feeding the RHS (the common case). - periodic_disc_max_step = _periodic_time_disc_max_step(sbml_model, func_defs, eval_ctx) + # no time-dependent floor/ceil/modulo feeding the RHS (the common case), and + # since GH #274 also None when every such node's influence on the RHS is + # already gated by one of the roots just registered — hence `seen_disc`, + # which is what makes that judgement possible rather than a guess. + periodic_disc_max_step = _periodic_time_disc_max_step( + sbml_model, func_defs, eval_ctx, seen_disc + ) # ── 10.6. rateOf csymbol support (GH #106) ───────────────────────── # Enable the live-derivative path iff the model actually *references* diff --git a/python/tests/test_sbml_periodic_floor_dosing.py b/python/tests/test_sbml_periodic_floor_dosing.py index ca7adb49d..d5badd2bf 100644 --- a/python/tests/test_sbml_periodic_floor_dosing.py +++ b/python/tests/test_sbml_periodic_floor_dosing.py @@ -41,12 +41,32 @@ it does bind, costs up to 2.5x the steps (MODEL0406553884: 439,367 vs 176,919) for no change in the answer. -So item 3 is now the opposite assertion, on the class GH #262 flagged as the one -roots provably cannot reach — a model carrying a bound and NO discontinuity roots -at all — plus one rooted model where the bound demonstrably binds. Both say the -same thing: the bound changes step selection and nothing else. Whether it should -then be narrowed or retired is a separate call; these pin the measurement so it -cannot rot silently the way the necessity claim did. +GH #274 acted on that: the bound is now derived ONLY where the roots provably +cannot reach. The test is value-position reachability from the ODE RHS through +the assignment rules, pruning every piecewise condition that IS a registered +root (:func:`_sbml_loader._periodic_disc_escapes_roots`). A floor/modulo whose +only influence is a rooted condition needs no bound — the root already forces +the stop; one whose *value* reaches the RHS (``dose = D * frac``) does, because +its jump is an RHS discontinuity no relational brackets. On the corpus that +keeps the bound on 10 of the 25 and drops it on 15, cutting 1,015,702 internal +steps to 519,472 with no answer moving. + +Which means most of this module now exercises the NO-bound path, and the +closed-form agreement below is what says dropping it there was safe. Two +fixtures carry the other side: + + * ``_EQ_GATED_PULSE`` — the necessity witness, and the reason "retire" would + have been wrong. Its pulse is gated by an EQUALITY on a floor + (``sub - 20*cyc == 0``, after BIOMD0000000589's ``i == 0`` cycle-index + selector), which the inequality-only emitters cannot root, while an + unrelated ``time < 5`` threshold IS rooted. So the model has roots and still + needs the bound: with it, y(10) = 36.788 against an exact 36.788; without, + 182.212 — 395% high, at BOTH tolerances, which is the step-over signature + (it never sees the pulses, so tightening rtol cannot help). A predicate that + merely counted roots would drop the bound here and return the wrong number. + * ``BIOMD0000000312`` — a real model in the kept set, where the bound binds + and disabling it changes nothing. The keep set is deliberately conservative; + this pins that. Oracle is closed-form. For ``dy/dt = (g - kd·dose(t))·y`` with ``dose = D`` while ``frac(time) = time - floor(time) < w`` (else 0) and constant ``kd``, each unit @@ -163,24 +183,147 @@ def _floor_dose_sbml(width: float) -> str: """ +# ── GH #274: a pulse the roots cannot reach ───────────────────────────────── +# Same dy/dt and the same closed form, but the dose window is selected by an +# EQUALITY on a floor rather than an inequality on a fraction. The root scan +# emits conditions for relationals (`<`, `<=`, `>`, `>=`); an `==` on a step +# function is not one, so this edge gets no root however many the model has — +# and it has one, from the deliberately unrelated `time < 5` threshold driving +# `z`. That combination (rooted model, unrooted pulse) is what separates a +# root-COUNT predicate from a per-path one. +# +# cyc = floor(time) integer cycle +# sub = floor(time * 20) 1/20-of-a-cycle index +# dose = D while sub - 20*cyc == 0 (the first 0.05 of every cycle) +EQ_K, EQ_W, EQ_D = 20, 0.05, 4.0 + +_EQ_GATED_PULSE = f""" + + + + + + + + + + + + + + + + + + + + + + + t + + + + + + + + t + {EQ_K} + + + + + + + + + D + + + sub + {EQ_K}cyc + + 0 + + + 0 + + + + + + + 1 + + t + 5 + + 0 + + + + + + + gkddose + y + + + + + other + + + +""" + + def _floor_oracle(n: int) -> float: """Closed-form y(n) at integer n for the periodic floor-dose model.""" return Y0 * math.exp(n * (GVAL - KDVAL * DVAL * WVAL)) +def _eq_gated_oracle(n: int) -> float: + """The same closed form for the equality-gated fixture's own D and w.""" + return Y0 * math.exp(n * (GVAL - KDVAL * EQ_D * EQ_W)) + + def _y(result): names = list(result.species_names) return np.asarray(result.species)[:, names.index("y")] def test_loader_derives_step_bound_below_window(): - model = bngsim.Model.from_sbml_string(_floor_dose_sbml(WVAL)) + """The derivation arithmetic, on a schedule that still gets a bound (#274). + + Measured on the equality-gated fixture rather than ``_floor_dose_sbml``, + whose pulse condition is a rootable inequality and so no longer earns one. + """ + model = bngsim.Model.from_sbml_string(_EQ_GATED_PULSE) ms = model._periodic_disc_max_step assert ms is not None - # Bound must keep a step from spanning the W-wide window; the loader targets - # window/3, so it sits comfortably below W. - assert 0.0 < ms < WVAL - assert ms == pytest.approx(WVAL / 3.0, rel=1e-3) + # Bound must keep a step from spanning the window; the loader targets + # window/3, so it sits comfortably below it. + assert 0.0 < ms < EQ_W + assert ms == pytest.approx(EQ_W / 3.0, rel=1e-3) + + +def test_rooted_pulse_condition_gets_no_bound(): + """``_floor_dose_sbml``'s edge is an inequality, so #72's scan roots it and + #274 declines the bound: the root already forces the stop it would buy. + + The closed-form tests below run this same fixture unbounded and still hit + the oracle, which is what makes the decline safe rather than merely cheap. + """ + model = bngsim.Model.from_sbml_string(_floor_dose_sbml(WVAL)) + assert model._core.n_discontinuity_triggers >= 1 + assert model._periodic_disc_max_step is None def test_time_threshold_schedule_gets_no_periodic_bound(): @@ -191,11 +334,64 @@ def test_time_threshold_schedule_gets_no_periodic_bound(): def test_clone_preserves_step_bound(): - model = bngsim.Model.from_sbml_string(_floor_dose_sbml(WVAL)) + # On a fixture that actually carries a bound — cloning `None` to `None` + # would pass without testing anything. + model = bngsim.Model.from_sbml_string(_EQ_GATED_PULSE) + assert model._periodic_disc_max_step is not None clone = model.clone() assert clone._periodic_disc_max_step == model._periodic_disc_max_step +def test_equality_gated_pulse_still_needs_the_bound(): + """GH #274's necessity witness — and why the bound is narrowed, not retired. + + The pulse edge is ``sub - 20*cyc == 0``. An equality on a step function is + not something the root scan emits, so nothing brackets this edge, and the + step bound is the only thing standing between the integrator and a 0.05-wide + window. The model nevertheless HAS a root (the unrelated ``time < 5`` + threshold driving ``z``), which is exactly the shape a root-count predicate + gets wrong. + + Disabling the bound does not degrade the answer, it destroys it: 182.2 + against an exact 36.8, identically at rtol 1e-9 and 1e-11. Tol-stably wrong + is the signature of a pulse that is never sampled at all. + """ + model = bngsim.Model.from_sbml_string(_EQ_GATED_PULSE) + assert model._core.n_discontinuity_triggers >= 1, "the unrelated time<5 root" + assert model._periodic_disc_max_step is not None, ( + "a rooted model can still have an unrooted pulse — this is the case a " + "root-count predicate would drop the bound on" + ) + + exact = _eq_gated_oracle(10) + for rtol in (1e-9, 1e-11): + r = bngsim.Simulator(bngsim.Model.from_sbml_string(_EQ_GATED_PULSE), method="ode").run( + t_span=(0.0, 10.0), + n_points=11, + rtol=rtol, + atol=1e-12, + max_steps=10_000_000, + timeout=120, + ) + got = _y(r)[-1] + assert got == pytest.approx(exact, rel=1e-3), f"bound on, rtol={rtol}: {got}" + + off = bngsim.Simulator(bngsim.Model.from_sbml_string(_EQ_GATED_PULSE), method="ode").run( + t_span=(0.0, 10.0), + n_points=11, + rtol=rtol, + atol=1e-12, + max_steps=10_000_000, + max_step=-1, + timeout=120, + ) + stepped_over = _y(off)[-1] + assert stepped_over > 3.0 * exact, ( + f"rtol={rtol}: expected the unbounded arm to jump the pulses " + f"(~{Y0 * math.exp(10 * GVAL):.1f}), got {stepped_over}" + ) + + @pytest.mark.parametrize("jac", ["analytical", "fd"]) @pytest.mark.parametrize("rtol", [1e-9, 1e-11]) def test_periodic_pulse_resolved_matches_closed_form(jac, rtol): @@ -247,7 +443,10 @@ def test_model1708310001_converges_to_segmented_oracle(jac): cycle arithmetic) drives y; resolved, bngsim reaches the exact segmented answer y(100)=953.07 tol-stably for either Jacobian (pre-fix: 1603 analytical / 1578 fd at the sweep tol, only converging to 953 as tol→1e-11).""" - assert bngsim.Model.from_sbml(_MODEL1708)._periodic_disc_max_step is not None + # GH #274 — this model no longer gets a bound at all: every disc node in its + # cycle arithmetic reaches the RHS only through conditions #259 roots. The + # oracle agreement below is now the evidence that the drop was safe. + assert bngsim.Model.from_sbml(_MODEL1708)._periodic_disc_max_step is None for rtol in (1e-9, 1e-11): # fresh state each run (run() does not reset a reused model) m = bngsim.Model.from_sbml(_MODEL1708) @@ -279,8 +478,9 @@ def test_model1708310001_roots_resolve_the_schedule_without_the_bound(): So this asserts what is true now — the roots alone resolve the schedule — and it still fails loudly if they regress, because the step-over answer is - 68% high. The bound itself is still derived (the test above pins that) and - still on by default; whether any model still *needs* it is GH #262. + 68% high. Since GH #274 the ``max_step=-1`` here is redundant rather than + contrarian: this model gets no bound to disable. Keeping it explicit means + the test still pins root coverage if the derivation ever comes back. """ for rtol in (1e-9, 1e-11): m = bngsim.Model.from_sbml(_MODEL1708) @@ -334,10 +534,15 @@ def _final_state_rel_diff(a, b) -> float: return float(np.max(np.abs(a - b) / denom)) -def _two_arm(path: str, t_end: float, n_points: int): - """(bound-on, bound-off) x (rtol, rtol/100) final states and step counts.""" +def _two_arm(path: str, t_end: float, n_points: int, on_max_step: float | None = None): + """(bound-on, bound-off) x (rtol, rtol/100) final states and step counts. + + ``on_max_step`` forces a specific ceiling for the "on" arm. Needed for a + model GH #274 declines to bound, where the default arm is already unbounded + and the comparison would otherwise be against itself. + """ out = {} - for arm, max_step in (("on", None), ("off", -1.0)): + for arm, max_step in (("on", on_max_step), ("off", -1.0)): for tag, rtol in (("loose", 1e-8), ("tight", 1e-10)): kw = dict( t_span=(0.0, t_end), @@ -361,17 +566,20 @@ def _two_arm(path: str, t_end: float, n_points: int): @pytest.mark.skipif( not os.path.exists(_rr_model("BIOMD0000000312")), reason="BIOMD0000000312 SBML not present" ) -def test_bound_unnecessary_even_where_no_root_can_reach(): - """The one class GH #262 kept the bound for turns out not to need it either. - - #262's caveat was that "the bound also covers schedules whose thresholds the - loader cannot emit as an ExprTk condition at all", so a model with a bound - and *zero* discontinuity roots is where the bound is the only mechanism in - play — nothing else can be resolving its pulses. BIOMD0000000312 is exactly - that, and its bound is not inert (the two arms take different numbers of - steps, which is what keeps this test from being vacuous). Disabling it - changes the answer by ~1e-10, four orders below either arm's own tolerance - stability. +def test_kept_bound_is_conservative_on_a_zero_root_model(): + """A model GH #274 KEEPS the bound for, where measurement says it need not. + + BIOMD0000000312 has a bound and *zero* discontinuity roots, so nothing else + can be resolving its schedule and #274's predicate keeps the bound — as it + must, since no root gates anything here. Empirically the bound is not needed: + it binds (the arms take different step counts, which is what keeps this test + from being vacuous) and disabling it moves the answer by ~1e-10, four orders + below either arm's own tolerance stability. + + Both halves are the point. The predicate is deliberately conservative — it + asks whether a root *provably* covers the edge, not whether the model happens + to integrate fine without one — and this pins the gap between the two so a + future tightening has a witness to argue against. """ path = _rr_model("BIOMD0000000312") model = bngsim.Model.from_sbml(path) @@ -397,28 +605,37 @@ def test_bound_unnecessary_even_where_no_root_can_reach(): @pytest.mark.skipif( not os.path.exists(_rr_model("MODEL0847869198")), reason="MODEL0847869198 SBML not present" ) -def test_bound_costs_steps_and_buys_nothing_on_a_rooted_model(): - """Same verdict on the rooted side, and here the price is visible. - - MODEL0847869198 carries both a bound and four discontinuity roots. The bound - binds hard — it forces >20% more internal steps — and buys an answer that - matches the unbounded one to ~1e-9. This is the shape of the cost #262 - weighed: ``max_step`` is a blunt instrument that shortens every step over the - whole horizon, not just the ones near a pulse edge. +def test_dropped_bound_matches_a_finely_bounded_run(): + """A model GH #274 drops the bound on, checked against a forced ceiling. + + MODEL0847869198 has four discontinuity roots, and every disc node in its + schedule reaches the RHS only through conditions those roots cover — so the + loader no longer bounds it. The claim that has to hold is not "the old bound + was harmless" but the stronger "running unbounded equals running with a step + ceiling fine enough to resolve anything": 0.05 here, well under the 0.333 the + pre-#274 loader derived, so it cannot be the ceiling that is too coarse to + disagree. + + The price the drop recovers is visible in the step counts — the forced + ceiling costs >20% more internal steps for an answer that matches to ~1e-9. + ``max_step`` shortens every step over the whole horizon, not just the ones + near a pulse edge. """ path = _rr_model("MODEL0847869198") model = bngsim.Model.from_sbml(path) - assert model._periodic_disc_max_step is not None assert model._core.n_discontinuity_triggers > 0 + assert model._periodic_disc_max_step is None, ( + "expected #274 to drop this model's bound — its schedule is fully rooted" + ) - arms = _two_arm(path, t_end=100.0, n_points=101) + arms = _two_arm(path, t_end=100.0, n_points=101, on_max_step=0.05) steps_on, steps_off = arms["on_tight"][1], arms["off_tight"][1] assert steps_on > 1.2 * steps_off, ( - f"expected the bound to cost steps ({steps_on} with vs {steps_off} without); " - "if it no longer binds, this model has stopped being a cost witness" + f"expected the forced ceiling to cost steps ({steps_on} with vs " + f"{steps_off} without); if it no longer binds, lower it" ) for arm in ("on", "off"): stability = _final_state_rel_diff(arms[f"{arm}_loose"][0], arms[f"{arm}_tight"][0]) assert stability < 1e-5, f"{arm} arm is not tol-stable ({stability:.2e})" diff = _final_state_rel_diff(arms["on_tight"][0], arms["off_tight"][0]) - assert diff < 1e-6, f"the bound changes this model's answer by {diff:.2e}" + assert diff < 1e-6, f"unbounded disagrees with a fine ceiling by {diff:.2e}"