From fed42833ba0cea04761bd524c1736145676b6272 Mon Sep 17 00:00:00 2001 From: Bill Hlavacek Date: Mon, 27 Jul 2026 23:14:44 -0600 Subject: [PATCH] fix(codegen): bound the sensitivity df/dp derivation with a #95-style budget (#90) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit #95/#187 budget the build-time symbolic derivation of the analytical Jacobian, so a model that does not derive in time falls back to the finite-difference Jacobian instead of hanging the load. The df/dp path #55 added (#65/#66/#67/#68) does the same kind of sympy work and had none: `_codegen.py` contained zero references to `deadline`. Its worst site runs one `sp.diff` per (distinct rate law, parameter it reads) pair — the product of rate-law size and parameter count, the axis #95 found super-linear — so a genome-scale Functional model with a long `sensitivity_params` list presented as a build that HANGS rather than one that declines and says why. The derivation is now bounded by `BNGSIM_SENS_DERIV_BUDGET_S`, resolved by `_sens_derivation_budget_s` off the same base, slope and override grammar as the Jacobian's (one implementation, so the two cannot drift) but read independently, so `BNGSIM_JAC_DERIV_BUDGET_S=inf` does not silently uncap it. The deadline is threaded to every sympy site on the build — the Functional rate laws, both flavours of derived-parameter chain rule, the model path and the .net text path — and checked on entry to each rate law as well as before each partial, so overshoot is bounded to one (law, parameter) pair. Expiry declines through the existing `_warn_functional_sens_rhs_refused`, naming how far it got and the override. It also covers the derivation the issue's three-site table does not list: the `J*yS` half re-derives df/dx through `_functional_jacobian_groups`, which is the same math `attach_functional_jacobian` ran at load under #95 but with no clock of its own — and a model whose load-time attach was itself cut off has no earlier bound to inherit, so bounding only df/dp left the identical hang one call later. The #151 native saturable emitters run no sympy and are untouched. One difference from the Jacobian budget is deliberate: this one never becomes unbounded by size. `_FD_NONVIABLE_SPECIES` exists because past that scale an FD Jacobian does not converge — there is nothing to fall back to. Declining a sensitivity RHS only hands the columns to CVODES' difference quotient, which is what every Functional model used before #55 and is correct at every scale. Not cheap (9-37x per column, and at tight rtol the DQ's ~sqrt(rtol) accuracy collapses the step size outright), which is why the decline warns and names the knob. Corpus A/B over all 585 .net models, default budget vs unbounded: 544 emit an analytic sens RHS in both arms, 0 source hashes differ, total derivation 7.8s with the slowest model at 0.46s — 46x headroom under the 20s base. An explicit override also joins the .net path's in-process memo key and its on-disk model_hash, which key on model content rather than on the generated C. --- .github/workflows/mir.yml | 4 + .github/workflows/windows-tail.yml | 4 + CHANGELOG.md | 56 +++ docs/user-guide/solvers.md | 26 ++ python/bngsim/_codegen.py | 239 +++++++++++- python/bngsim/_jacobian.py | 104 +++++- python/tests/test_codegen_sens_budget.py | 447 +++++++++++++++++++++++ 7 files changed, 842 insertions(+), 38 deletions(-) create mode 100644 python/tests/test_codegen_sens_budget.py diff --git a/.github/workflows/mir.yml b/.github/workflows/mir.yml index 40da236..034fe48 100644 --- a/.github/workflows/mir.yml +++ b/.github/workflows/mir.yml @@ -78,6 +78,9 @@ on: - "python/bngsim/_codegen.py" - "python/bngsim/_simulator.py" - "python/bngsim/_switch_sensitivity.py" + # In _CODEGEN_SOURCE_MODULES alongside the two above, i.e. it decides what + # gets emitted (the derivation budgets and every symbolic derivative). + - "python/bngsim/_jacobian.py" - "python/bngsim/__init__.py" - "python/tests/test_mir_equivalence.py" - "python/tests/test_codegen*.py" @@ -153,6 +156,7 @@ jobs: python/tests/test_codegen_functional_sens_rhs.py \ python/tests/test_codegen_switch_condition_sens.py \ python/tests/test_codegen_mm_sens.py \ + python/tests/test_codegen_sens_budget.py \ python/tests/test_codegen_determinism.py \ -rA --tb=short ${{ github.event.inputs.pytest_args }} \ 2>&1 | tee mir-${{ matrix.os }}.log diff --git a/.github/workflows/windows-tail.yml b/.github/workflows/windows-tail.yml index e37c1fe..22b00c7 100644 --- a/.github/workflows/windows-tail.yml +++ b/.github/workflows/windows-tail.yml @@ -42,6 +42,9 @@ on: branches: [main] paths: &tail_paths - "python/bngsim/_codegen.py" + # In _CODEGEN_SOURCE_MODULES alongside _codegen.py, i.e. it decides what + # gets emitted (the derivation budgets and every symbolic derivative). + - "python/bngsim/_jacobian.py" - "python/bngsim/_simulator.py" - "python/bngsim/_version.py" - "python/bngsim/_net_reader.py" @@ -128,6 +131,7 @@ jobs: python/tests/test_codegen_chunking.py \ python/tests/test_codegen_jacobian.py \ python/tests/test_model_codegen_sensitivity.py \ + python/tests/test_codegen_sens_budget.py \ python/tests/test_ssa_propensity_export.py \ -rA --tb=long ${{ inputs.pytest_args }} \ 2>&1 | tee tail-windows.log diff --git a/CHANGELOG.md b/CHANGELOG.md index d01756f..0bef798 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -216,6 +216,62 @@ in `CMakeLists.txt`) is derived from it. re-bootstrap after a pin bump fail on step one. ### Added +- **A build-time derivation budget for the sensitivity `∂f/∂p` path (issue #90), + the last unaddressed "still has to bail" item of #55.** #95/#187 bound the + symbolic derivation of the analytical *Jacobian*, so a model that does not + derive in time falls back to the finite-difference Jacobian instead of hanging + the load. The `∂f/∂p` path added by #55 (#65/#66/#67/#68) does the same kind of + sympy work and had no budget at all: `python/bngsim/_codegen.py` contained zero + references to `deadline`, and its three `sp.diff` sites were unguarded. The one + that scales badly runs **one `sp.diff` per (distinct rate law, parameter it + reads) pair** — the product of rate-law size and parameter count, exactly the + axis #95 found super-linear — so a genome-scale Functional model with a long + `sensitivity_params` list would present as a build that *hangs* rather than one + that declines and says why (and to a `rr_parity`-style harness, which times + build and solve together, as an ODE timeout). + + The derivation is now bounded by `BNGSIM_SENS_DERIV_BUDGET_S`. It shares the + Jacobian budget's base, per-species slope and override grammar — one + implementation, so the two policies cannot drift — but is resolved + independently, so `BNGSIM_JAC_DERIV_BUDGET_S=inf` (the documented genome-scale + workaround for keeping an analytical Jacobian) does not silently uncap it. + The deadline is threaded down to every `sp.diff` on the path — the Functional + rate laws and both flavours of derived-parameter chain rule, on the model path + and the `.net` text path alike — and checked on entry to each rate law as well + as before each partial, so overshoot is bounded to one (law, parameter) pair + rather than one whole model. Expiry declines through the existing + `_warn_functional_sens_rhs_refused`, landing beside every other decline reason + with a warning naming how far the derivation got and the override to raise. + + It also covers the *other* derivation on the same build, which the issue's + three-site table does not list: `ySdot = J·yS + ∂f/∂p`, and the `J·yS` half + re-derives ∂f/∂x through `_functional_jacobian_groups`. That is the same math + `attach_functional_jacobian` already ran at load under the #95 budget, but a + re-derivation with no clock of its own — and on a model whose load-time attach + was itself cut off there is no earlier bound to inherit, so bounding only ∂f/∂p + would have left the identical hang reachable one call later. The #151 native + saturable emitters run no sympy and are unaffected; only their fallback takes + the deadline. + + **One difference from the Jacobian budget is deliberate: this one never becomes + unbounded by size.** `_FD_NONVIABLE_SPECIES` exists because past that scale an + FD Jacobian does not converge — there is nothing to fall back *to*, so the + analytical Jacobian is mandatory. Declining a sensitivity RHS only hands the + columns to CVODES' internal difference quotient, which is what every Functional + model used before #55 and which is correct at every scale. Correct is not cheap + (measured 9–37x per column, and on a stiff model at a tight `rtol` the DQ's + `~sqrt(rtol)` accuracy collapses the step size outright), which is why the + decline is a warning that names the knob rather than a silent downgrade. + + Corpus A/B over all 585 `.net` models, default budget vs unbounded: **544 models + emit an analytic sensitivity RHS in both arms, 0 source hashes differ**, total + derivation 7.6 s with the slowest single model at 0.44 s — 46x headroom under + the 20 s base, so nothing real is near the cut-off. An explicit override also + joins the `.net` path's in-process memo key and its on-disk `model_hash` + (that path keys on the model's *content*, not on the generated C, so without it + a build made under a tight budget would be served back to one made without it — + the same trap #67's A/B hatch had to sidestep). + - **Per-model composition counts in the ODE Jacobian characterization report (issue #42).** `jacobian_characterization.py` recorded rich Jacobian-side structure (`N`, `rank`, `n_reactions`, `density`, `stiffness_ratio_*`) but diff --git a/docs/user-guide/solvers.md b/docs/user-guide/solvers.md index 416f458..6d5074d 100644 --- a/docs/user-guide/solvers.md +++ b/docs/user-guide/solvers.md @@ -78,6 +78,32 @@ treated as mandatory and is always derived to completion. Override the budget wi `BNGSIM_JAC_DERIV_BUDGET_S` (seconds, or `inf`/`none`/`0` to disable it entirely — the manual genome-scale escape hatch). +**Sensitivity derivation budget (GH #90)**: a sensitivity run derives a second +set of terms at build time — both halves of the compiled sensitivity RHS's +`ySdot = J·yS + ∂f/∂p` — and that derivation has its own budget, +`BNGSIM_SENS_DERIV_BUDGET_S`. It +shares the Jacobian budget's base and its per-species scaling, and takes the same +values (seconds, or `inf`/`none`/`0` for unbounded), but the two are independent: +raising one does not raise the other. + +One difference is deliberate. The sensitivity budget **never becomes unbounded**, +at any model size. Where a dropped analytical Jacobian falls back to something +that may not converge at all, a declined sensitivity RHS falls back to CVODES' +internal difference quotient, which is correct at every scale — so there is never +a reason to let this derivation run without a bound, and a genome-scale model with +Functional rate laws gets a decline it can read rather than a build that appears +to hang. + +The fallback is correct but not cheap: measured 9–37x per sensitivity column, and +on a stiff model at a tight `rtol` the difference quotient's own `~sqrt(rtol)` +accuracy collapses the step size (it can exhaust `max_steps`). So the decline is +logged as a warning naming this variable — if a model you care about reports it, +raise the budget rather than accept the slower path: + +```bash +BNGSIM_SENS_DERIV_BUDGET_S=inf python fit.py +``` + ## Logging ```python diff --git a/python/bngsim/_codegen.py b/python/bngsim/_codegen.py index ada0951..c0a0596 100644 --- a/python/bngsim/_codegen.py +++ b/python/bngsim/_codegen.py @@ -1414,6 +1414,79 @@ def _inline_derived_param_refs( return s +def _check_derivation_deadline(deadline: float | None) -> None: + """Raise :class:`_DerivationBudgetExceeded` if the build-time symbolic + derivation has run past ``deadline`` (GH #90). + + ``deadline`` is a ``time.perf_counter()`` stamp, or ``None`` for unbounded — + which is what every caller outside the sensitivity build passes, so this is a + no-op for them. One spelling for every check site so the sensitivity budget + cannot drift from the Jacobian's, whose per-observable check + (``_jacobian.differentiate_rate_law``) this mirrors. + """ + if deadline is not None and time.perf_counter() > deadline: + from bngsim._jacobian import _DerivationBudgetExceeded + + raise _DerivationBudgetExceeded + + +def _sens_derivation_deadline(n_species: int) -> float | None: + """Absolute ``time.perf_counter()`` deadline for one sensitivity-RHS build, or + ``None`` when the budget is disabled (GH #90). + + Resolved once per :func:`generate_sens_from_model` / :func:`generate_sens_rhs_c` + call and threaded down, so every ``sp.diff`` on the ∂f/∂p path shares a single + wall-clock bound instead of each site getting its own. + """ + from bngsim._jacobian import _sens_derivation_budget_s + + budget = _sens_derivation_budget_s(n_species=n_species) + return None if budget is None else time.perf_counter() + budget + + +def _sens_budget_cache_tag() -> str: + """Cache-key fragment for an explicitly overridden sensitivity derivation + budget (GH #90), or ``""`` when the env var is unset. + + The budget decides whether a model gets an analytic sensitivity RHS at all, so + it belongs in the key of any cache that is not content-addressed on the + generated source — the ``.net`` path's in-process memo and its on-disk + ``model_hash``, which both key on the .net's *content*. Without it a build made + under a deliberately tight budget would be served back to one made without it, + the same trap ``functional_sens_rhs_enabled`` already sidesteps (GH #67). + Empty when unset, so the default key — and every ``.so`` already cached — is + byte-identical to before. + + This does not make a *default*-budget expiry cache-safe: the budget is + wall-clock, so a model that derives near the limit can emit on one run and + decline on the next, and whichever came first is what the .net path cached. + Raising the budget is the fix, and doing so lands in a fresh namespace. + """ + from bngsim._jacobian import _SENS_BUDGET_ENV + + raw = os.environ.get(_SENS_BUDGET_ENV) + return "" if raw is None else f":sens_budget={raw.strip().lower()}" + + +def _sens_budget_decline_reason(n_species: int, progress: str) -> str: + """The decline reason for a build-time ∂f/∂p derivation that ran past its + budget (GH #90), phrased for :func:`_warn_functional_sens_rhs_refused`. + + A budget expiry is reported through that same channel as every other decline + rather than one of its own, and names both how far the derivation got and the + override — the alternative to declining is a build that appears to hang, which + is the whole point of the budget. + """ + from bngsim._jacobian import _SENS_BUDGET_ENV, _sens_derivation_budget_s + + budget = _sens_derivation_budget_s(n_species=n_species) + return ( + f"the build-time ∂f/∂p derivation exceeded its {budget:g}s budget " + f"({progress}); set {_SENS_BUDGET_ENV} to raise or disable it (seconds, or " + "inf/none/0 for unbounded)" + ) + + def _compute_derived_param_jacobian( expr: str, primary_param_names: set, @@ -1473,6 +1546,7 @@ def _derived_param_jacobian_checked( primary_param_names: set, param_idx: dict, derived_exprs: dict[str, str] | None = None, + deadline: float | None = None, ) -> tuple[dict[str, str] | None, str | None]: """:func:`_compute_derived_param_jacobian`, plus the reason it gave up. @@ -1485,10 +1559,20 @@ def _derived_param_jacobian_checked( an analytic sensitivity RHS must refuse the whole RHS on a failure (falling back to CVODES' correct-but-slower internal difference quotient) rather than ship a gradient component that is confidently, exactly wrong. + + ``deadline`` (GH #90) is a ``time.perf_counter()`` stamp bounding the caller's + whole build-time derivation. It is checked once on entry (parsing an + expression is itself unbounded work) and again before each ``sp.diff``, so a + single pathological expression overshoots by at most one partial. Expiry + raises :class:`bngsim._jacobian._DerivationBudgetExceeded` rather than + returning a reason: it is a property of the *build*, not of this expression, + and it must unwind past the per-expression caches to decline the whole model. + Callers that pass no deadline (the default) never see it. """ s = expr.strip() if not s: return None, None + _check_derivation_deadline(deadline) try: import sympy as sp from sympy.parsing.sympy_parser import parse_expr @@ -1567,6 +1651,7 @@ def _derived_param_jacobian_checked( result: dict[str, str] = {} for p_name in referenced: + _check_derivation_deadline(deadline) deriv = sp.diff(sym_expr, sym_map[sym_name_of[p_name]]) if deriv == 0: continue @@ -2284,6 +2369,10 @@ def generate_sens_rhs_c(net_path: str) -> str | None: For Functional/MM: returns None (fall back to CVODES internal FD). + The only symbolic work here is the derived-rate-constant chain rule below; it + shares the GH #90 build-time budget with the model path, so a .net carrying + enough ``# ConstantExpression`` rate constants declines rather than hangs. + Parameters ---------- net_path : str @@ -2305,6 +2394,11 @@ def generate_sens_rhs_c(net_path: str) -> str | None: n_sp = len(species) n_params = len(params) + # GH #90: one deadline for this build's symbolic work, resolved before it. + from bngsim._jacobian import _DerivationBudgetExceeded + + deadline = _sens_derivation_deadline(n_sp) + # Build name→index maps param_idx = {name: i for i, (_, name, _, _) in enumerate(params)} func_names = {name for _, name, _ in functions} @@ -2339,9 +2433,21 @@ def generate_sens_rhs_c(net_path: str) -> str | None: # standard deviation) has no bearing here. if is_const or name not in rate_const_names: continue - jac, reason = _derived_param_jacobian_checked( - expr, primary_param_names, param_idx, derived_exprs=derived_exprs - ) + try: + jac, reason = _derived_param_jacobian_checked( + expr, + primary_param_names, + param_idx, + derived_exprs=derived_exprs, + deadline=deadline, + ) + except _DerivationBudgetExceeded: + # GH #90: decline to CVODES' difference quotient rather than let the + # chain-rule derivation run unbounded (see generate_sens_from_model). + _warn_functional_sens_rhs_refused( + _sens_budget_decline_reason(n_sp, f"deriving the rate constant {name} = {expr!r}") + ) + return None if reason is not None: # Issue #56: emitting the RHS without this chain rule would report # ∂/∂primary as exactly zero. Refuse the analytic RHS instead so the @@ -5052,7 +5158,9 @@ def _mm_jacobian_groups(plan_mm, add) -> list[list[str]] | None: return groups -def _functional_jacobian_groups(core, data, add) -> tuple[list[list[str]], ...] | None: +def _functional_jacobian_groups( + core, data, add, deadline: float | None = None +) -> tuple[list[list[str]], ...] | None: """Reconstruct every Functional reaction's Jacobian contribution as balanced ``{ … }`` C line-groups, scattered through the caller's ``add``. @@ -5074,6 +5182,13 @@ def _functional_jacobian_groups(core, data, add) -> tuple[list[list[str]], ...] rule and the per-observable product rule follow ``set_functional_jacobian`` / ``scatter_functional_observable_terms``, with the derivative math from the native saturable C emitters (GH #151). + + ``deadline`` (GH #90) bounds the SymPy fallback of that derivative math. This + is a *re*-derivation — ``attach_functional_jacobian`` already ran it at load + under the #95 budget — but a re-derivation that ignores its own clock, and on + a model whose load-time attach was itself cut off there is no earlier bound to + inherit. The sensitivity RHS passes its build's deadline; the Jacobian emitter + passes ``None`` and is unchanged. """ from collections import Counter @@ -5187,7 +5302,13 @@ def _scatter_existing(sp_j, i, coeff, live_idx, sdiv, rhs): continue affected = _net_affected(rxn) terms = build_per_species_c( - rxn["rate_expr"], func_map, obs_groups, species_meta, constants, resolve_symbol + rxn["rate_expr"], + func_map, + obs_groups, + species_meta, + constants, + resolve_symbol, + deadline, ) if terms is None: return None @@ -5252,7 +5373,7 @@ def _scatter_existing(sp_j, i, coeff, live_idx, sdiv, rhs): affected = _net_affected(rxn) rxn_idx = int(rxn["rxn_idx"]) od = differentiate_rate_law_c( - rxn["rate_expr"], func_map, ctx_obs_names, constants, resolve_symbol + rxn["rate_expr"], func_map, ctx_obs_names, constants, resolve_symbol, deadline ) if od is None: return None @@ -5808,6 +5929,9 @@ class _FunctionalDfdpScope(NamedTuple): ``switch_scope`` is GH #68's gate context: non-``None`` only for a model that has a condition to gate *and* whose clock/parameter view could be assembled. ``None`` keeps the pre-#68 behaviour — every condition declines the model. + + ``deadline`` is GH #90's build-time derivation bound, shared by every rate law + of the model (see :func:`_sens_derivation_deadline`); ``None`` is unbounded. """ func_map: dict[str, str] @@ -5817,6 +5941,7 @@ class _FunctionalDfdpScope(NamedTuple): primary_param_names: set[str] derived_exprs: dict[str, str] switch_scope: SwitchConditionScope | None = None + deadline: float | None = None def _functional_rate_law_partials( @@ -5826,8 +5951,12 @@ def _functional_rate_law_partials( Returns ``([(param_idx, c_expr)], None)`` — possibly empty, which is the *success* case for a rate law with no parameter dependence at all — or - ``(None, reason)`` naming what blocked it. Never raises, and never returns a - partial list: a rate law is either fully differentiated or refused (#56). + ``(None, reason)`` naming what blocked it. Never returns a partial list: a + rate law is either fully differentiated or refused (#56). The one exception + it raises is ``scope.deadline``'s :class:`_DerivationBudgetExceeded` (GH #90), + which is deliberately not a ``reason``: a reason is memoized per rate-law text + by the caller, and a wall-clock expiry is a property of the build, not of the + expression. The derivative is taken w.r.t. every parameter symbol surviving inlining, *including* a derived (ConstantExpression) one, whose own column is that @@ -5850,6 +5979,12 @@ def _functional_rate_law_partials( except ImportError: return None, "sympy is not installed" + # GH #90: bound the whole build, but check on entry as well as per-parameter — + # inlining, the construct scan and ``_exprtk_to_sympy`` are themselves + # unbounded work on a large enough law, so a check only at the ``sp.diff`` + # loop would let a single rate law overshoot arbitrarily. + _check_derivation_deadline(scope.deadline) + inlined = _inline_functions(rate_expr, scope.func_map) if inlined is None: return None, "function references form a cycle or nest deeper than 64 levels" @@ -5910,6 +6045,7 @@ def resolve_symbol(name: str) -> str | None: terms: list[tuple[int, str]] = [] for a in sorted(free & set(scope.param_of_alias)): pname = scope.param_of_alias[a] + _check_derivation_deadline(scope.deadline) deriv = sp.diff(sym_expr, sp.Symbol(a)) if deriv == 0: continue @@ -5933,6 +6069,7 @@ def resolve_symbol(name: str) -> str | None: scope.primary_param_names, scope.param_idx_by_name, derived_exprs=scope.derived_exprs, + deadline=scope.deadline, ) if why is not None: return None, ( @@ -5946,7 +6083,9 @@ def resolve_symbol(name: str) -> str | None: return terms, None -def _functional_dfdp_terms(core, data) -> tuple[dict[int, list[tuple[int, str]]], str | None]: +def _functional_dfdp_terms( + core, data, deadline: float | None = None +) -> tuple[dict[int, list[tuple[int, str]]], str | None]: """Differentiate every Functional rate law w.r.t. every parameter it reads. Returns ``({reaction_idx: [(param_idx, c_expr_for_∂func/∂p)]}, None)``, or @@ -5978,13 +6117,19 @@ def _functional_dfdp_terms(core, data) -> tuple[dict[int, list[tuple[int, str]]] nor ``time`` (this is what catches ``rateOf``, whose accessor is an evaluator variable and not a model parameter), a derivative sympy cannot render as C, and a derived parameter whose own Jacobian was lost. + + ``deadline`` (GH #90) bounds the sympy work: this loop runs one ``sp.diff`` + per (distinct rate law, parameter it reads) pair, the one axis on this path + that grows super-linearly with model size, so an unbudgeted genome-scale + Functional model would appear to hang the build rather than decline. Expiry + declines like any other reason, naming how far it got. """ params = data["parameters"] observables = data["observables"] functions = data["functions"] reactions = data["reactions"] - from bngsim._jacobian import has_condition_construct + from bngsim._jacobian import _DerivationBudgetExceeded, has_condition_construct ctx = core.functional_jacobian_context() frxn_by_idx = {int(r["rxn_idx"]): r for r in (ctx.get("functional_reactions") or [])} @@ -6055,6 +6200,7 @@ def _alias(n: str) -> str: primary_param_names=primary_param_names, derived_exprs=derived_exprs, switch_scope=switch_scope, + deadline=deadline, ) # A rule-generated network reuses one rate-law expression across many @@ -6095,7 +6241,19 @@ def _decline(reason: str) -> tuple[dict[int, list[tuple[int, str]]], str]: rate_expr = frxn["rate_expr"] hit = cache.get(rate_expr) if hit is None: - hit = _functional_rate_law_partials(rate_expr, scope) + try: + hit = _functional_rate_law_partials(rate_expr, scope) + except _DerivationBudgetExceeded: + # GH #90. The deadline is checked on entry to the differentiation + # as well as per-parameter, so this doubles as the between-laws + # check — and it costs nothing on a cache hit, which is what a + # rule-generated network is almost entirely made of. + return _decline( + _sens_budget_decline_reason( + len(data["species"]), + f"after {len(cache)} distinct rate law(s), at {label}", + ) + ) cache[rate_expr] = hit terms, why = hit if terms is None: @@ -6156,6 +6314,12 @@ def generate_sens_from_model(model, *, functional: bool = False) -> str | None: ``(∂p_d/∂primary_k) * sf * ∏y^m`` to the sensitivity of every primary parameter that appears in ``expr``, exactly mirroring the .net path's ``derived_expansion`` machinery. + + Every symbolic derivation this triggers — the Functional ∂func/∂p and both + flavours of derived-parameter chain rule — shares one wall-clock budget + (GH #90), so a model whose ``sp.diff`` work does not finish in time declines + with a warning and falls back to CVODES' internal difference quotient instead + of hanging the build. """ core = model._core if hasattr(model, "_core") else model data = core.codegen_data() @@ -6167,6 +6331,11 @@ def generate_sens_from_model(model, *, functional: bool = False) -> str | None: n_sp = len(species) n_params = len(params) + # GH #90: one deadline for the whole build, resolved before any sympy runs. + from bngsim._jacobian import _DerivationBudgetExceeded + + deadline = _sens_derivation_deadline(n_sp) + # Bail if any reaction is non-Elementary — analytical sens RHS is only # defined for k * sf * ∏y^m kinetics. Same constraint as # generate_sens_rhs_c (line 762-765) for the .net path. @@ -6193,7 +6362,7 @@ def generate_sens_from_model(model, *, functional: bool = False) -> str | None: functional_terms: dict[int, list[tuple[int, str]]] = {} functional_jacv_groups: list[list[str]] = [] if functional: - functional_terms, decline = _functional_dfdp_terms(core, data) + functional_terms, decline = _functional_dfdp_terms(core, data, deadline) if decline is not None: return None if functional_terms: @@ -6203,7 +6372,19 @@ def generate_sens_from_model(model, *, functional: bool = False) -> str | None: # A decline here is the #151 emitters' ("this derivative is not # representable in C"), which the ∂f/∂p pass cannot see: it # differentiates w.r.t. parameters, this one w.r.t. species. - groups = _functional_jacobian_groups(core, data, _jacv_add) + # GH #90: this half runs sympy too (the #151 native emitters cover the + # saturable family, everything else falls through to it), so it shares + # the build's deadline rather than being the one unbounded derivation + # left on the path. + try: + groups = _functional_jacobian_groups(core, data, _jacv_add, deadline) + except _DerivationBudgetExceeded: + _warn_functional_sens_rhs_refused( + _sens_budget_decline_reason( + n_sp, "deriving J*yS over the Functional reactions" + ) + ) + return None if groups is None: _warn_functional_sens_rhs_refused( "a Functional rate law's derivative with respect to the species it " @@ -6268,9 +6449,26 @@ def generate_sens_from_model(model, *, functional: bool = False) -> str | None: expr = p.get("expression", "") if not expr: continue - jac, reason = _derived_param_jacobian_checked( - expr, primary_param_names, param_idx_by_name, derived_exprs=derived_exprs - ) + try: + jac, reason = _derived_param_jacobian_checked( + expr, + primary_param_names, + param_idx_by_name, + derived_exprs=derived_exprs, + deadline=deadline, + ) + except _DerivationBudgetExceeded: + # GH #90: same budget as the Functional pass above, and the same + # outcome — decline to CVODES' difference quotient rather than let a + # model with thousands of derived rate constants hang the build. This + # loop is reached by Elementary-only models too, which have no other + # sympy on this path. + _warn_functional_sens_rhs_refused( + _sens_budget_decline_reason( + n_sp, f"deriving the rate constant {p['name']} = {expr!r}" + ) + ) + return None if reason is not None: # Issue #56 — see generate_sens_rhs_c: a dropped chain rule here # reads downstream as a hard zero, so refuse the analytic RHS and @@ -7347,7 +7545,7 @@ def prepare_ssa_propensity_lib(model, *, force_recompile: bool = False) -> str | # in place — so every flag is part of the key, and an entry for one combination # must never satisfy another. _PREPARE_CODEGEN_MEMO: dict[ - tuple[str, bool, bool, bool, bool], tuple[Path, tuple[tuple[str, int], ...], str] + tuple[str, bool, bool, bool, bool, str], tuple[Path, tuple[tuple[str, int], ...], str] ] = {} _PREPARE_CODEGEN_MEMO_LOCK = threading.Lock() @@ -7455,13 +7653,17 @@ def prepare_codegen(net_path: str, model=None, emit_jac: bool = True) -> Path: want_jac, want_outputs, want_output_sens = _codegen_emit_flags(model, emit_jac) # The GH #67 hatch is process-scoped, not file-scoped, so it belongs in the # in-process memo key as well as the on-disk one below — a test that flips - # it mid-process must not be handed the other variant's .so. + # it mid-process must not be handed the other variant's .so. GH #90's + # derivation-budget override is process-scoped in exactly the same way and + # decides exactly the same thing (whether the analytic sens RHS is emitted), + # so it rides along in both keys. memo_key = ( net_key, want_jac, want_outputs, want_output_sens, functional_sens_rhs_enabled(), + _sens_budget_cache_tag(), ) # Fast path (T2): an unchanged .net (and its .tfun deps) resolves to the @@ -7500,6 +7702,7 @@ def prepare_codegen(net_path: str, model=None, emit_jac: bool = True) -> Path: # so the default key — and every .so already in the cache — is unchanged. if not functional_sens_rhs_enabled(): suffix += ":no_functional_sens" + suffix += _sens_budget_cache_tag() if suffix: model_hash = hashlib.sha256((model_hash + suffix).encode()).hexdigest()[:16] diff --git a/python/bngsim/_jacobian.py b/python/bngsim/_jacobian.py index e887619..c1e352c 100644 --- a/python/bngsim/_jacobian.py +++ b/python/bngsim/_jacobian.py @@ -1096,14 +1096,19 @@ def _native_per_species_terms(rate_expr, func_map, obs_groups, species_amount, c def build_per_species_c( - rate_expr, func_map, obs_groups, species_amount, constant_names, resolve_symbol + rate_expr, func_map, obs_groups, species_amount, constant_names, resolve_symbol, deadline=None ) -> list[tuple[int, str]] | None: """Codegen per-species path: ``[(species_idx0, c_str)]`` or ``None``. Tries the native saturable family first (no SymPy), then the SymPy path (:func:`build_per_species_sympy` + :func:`sympy_to_c`). ``[]`` is a *success* (constant-rate ⇒ zero column). A native success whose C emission fails (an - unresolvable symbol) falls through to SymPy rather than declining outright.""" + unresolvable symbol) falls through to SymPy rather than declining outright. + + ``deadline`` (GH #90) is forwarded to the SymPy fallback only — the native + path runs no SymPy and needs no bound. Passed by the sensitivity RHS, whose + ``J·yS`` half re-derives this on its own build; ``None`` (the Jacobian + emitter's call) is the unbudgeted behaviour.""" from bngsim import _saturable_jacobian as _sat native = _sat.build_per_species_native( @@ -1122,7 +1127,7 @@ def build_per_species_c( return out terms = build_per_species_sympy( - rate_expr, func_map, obs_groups, species_amount, constant_names + rate_expr, func_map, obs_groups, species_amount, constant_names, deadline ) if terms is None: return None @@ -1136,10 +1141,13 @@ def build_per_species_c( def differentiate_rate_law_c( - rate_expr, func_map, observable_names, constant_names, resolve_symbol + rate_expr, func_map, observable_names, constant_names, resolve_symbol, deadline=None ) -> list[tuple[str, str]] | None: """Codegen per-observable path: ordered ``[(observable_name, c_str)]`` or - ``None``. Native saturable family first (no SymPy), then SymPy.""" + ``None``. Native saturable family first (no SymPy), then SymPy. + + ``deadline`` (GH #90) is forwarded to the SymPy fallback only; see + :func:`build_per_species_c`.""" from bngsim import _saturable_jacobian as _sat nd = _sat.differentiate_rate_law_native(rate_expr, func_map, observable_names, constant_names) @@ -1155,7 +1163,7 @@ def differentiate_rate_law_c( if emitted: return out - dd = differentiate_rate_law(rate_expr, func_map, observable_names, constant_names) + dd = differentiate_rate_law(rate_expr, func_map, observable_names, constant_names, deadline) if dd is None: return None out2: list[tuple[str, str]] = [] @@ -1264,10 +1272,58 @@ def differentiate_rate_law_c( _FD_COSTLY_SPECIES = 2000 +# GH #90: the sensitivity ∂f/∂p derivation (``_codegen._functional_dfdp_terms`` +# and the derived-parameter chain rules on the same build) gets its own budget, +# resolved by _sens_derivation_budget_s below. Same base, same slope, same +# override grammar as the Jacobian's — deliberately, so the two policies cannot +# drift — but its own env var, and one substantive difference: it never goes +# unbounded by size. +# +# The difference is about the FALLBACK, not the derivation. _FD_NONVIABLE_SPECIES +# exists because past that size there is nothing to fall back TO: an FD Jacobian +# needs ~n_species RHS evals per Newton setup and simply does not converge, so +# cutting the derivation off breaks the solve outright and the analytical Jacobian +# is mandatory. The sensitivity path has no such cliff — declining hands the +# columns to CVODES' internal difference quotient, which is what every Functional +# model used before #55 and is correct at every scale (measured 9-37x more +# expensive per column, not divergent). So the sensitivity budget stays finite +# everywhere: a build that would otherwise appear to HANG instead declines, says +# so, and solves. Anyone who would rather wait sets BNGSIM_SENS_DERIV_BUDGET_S. +_JAC_BUDGET_ENV = "BNGSIM_JAC_DERIV_BUDGET_S" +_SENS_BUDGET_ENV = "BNGSIM_SENS_DERIV_BUDGET_S" + + class _DerivationBudgetExceeded(Exception): """Internal signal: the build-time symbolic derivation passed its wall-clock budget. Caught by :func:`attach_functional_jacobian`, which logs the fallback - and leaves the model on the finite-difference Jacobian (GH #95).""" + and leaves the model on the finite-difference Jacobian (GH #95), and by + ``_codegen.generate_sens_from_model`` / ``generate_sens_rhs_c``, which decline + the analytic sensitivity RHS and leave the model on CVODES' internal + difference quotient (GH #90).""" + + +def _budget_env_override(env_var: str, default: float | None) -> float | None: + """Apply an explicit ``env_var`` budget over a size-derived ``default``. + + The override grammar, shared by every derivation budget: an absolute number of + seconds, or ``inf``/``none``/``off``/``0`` for unbounded (the pre-#95 and + documented genome-scale workaround). A non-positive or non-finite value also + disables the budget; a value that does not parse falls through to ``default``, + so a typo degrades to the policy rather than to no budget at all. + """ + raw = os.environ.get(env_var) + if raw is None: + return default + raw = raw.strip().lower() + if raw in ("inf", "none", "off", "0"): + return None + try: + val = float(raw) + except ValueError: + return default + if val <= 0 or val != val or val == float("inf"): + return None + return val def _derivation_budget_s(n_species: int = 0) -> float | None: @@ -1288,19 +1344,27 @@ def _derivation_budget_s(n_species: int = 0) -> float | None: else: default = max(_DEFAULT_DERIVATION_BUDGET_S, _BUDGET_PER_SPECIES_S * n_species) - raw = os.environ.get("BNGSIM_JAC_DERIV_BUDGET_S") - if raw is None: - return default - raw = raw.strip().lower() - if raw in ("inf", "none", "off", "0"): - return None - try: - val = float(raw) - except ValueError: - return default - if val <= 0 or val != val or val == float("inf"): - return None - return val + return _budget_env_override(_JAC_BUDGET_ENV, default) + + +def _sens_derivation_budget_s(n_species: int = 0) -> float | None: + """Resolve the build-time budget for the *sensitivity* ∂f/∂p derivation in + seconds, or ``None`` for unbounded (GH #90). + + The sensitivity counterpart of the Jacobian's :func:`_derivation_budget_s`, + sharing its base, its per-species slope and its override grammar — but read + from ``BNGSIM_SENS_DERIV_BUDGET_S``, and **never unbounded by species count**. + See the block comment above ``_SENS_BUDGET_ENV`` for why + ``_FD_NONVIABLE_SPECIES`` has no counterpart here: the sensitivity fallback + (CVODES' internal difference quotient) stays viable at every model size, so + there is never a reason to let this derivation run without a bound. + + The two budgets are independent: setting ``BNGSIM_JAC_DERIV_BUDGET_S=inf`` to + keep a genome-scale model's analytical Jacobian does not also uncap this one, + because they buy different things and fall back to different paths. + """ + default = max(_DEFAULT_DERIVATION_BUDGET_S, _BUDGET_PER_SPECIES_S * n_species) + return _budget_env_override(_SENS_BUDGET_ENV, default) def eager_jacobian_requested(defer_jacobian: bool | None = None) -> bool: diff --git a/python/tests/test_codegen_sens_budget.py b/python/tests/test_codegen_sens_budget.py new file mode 100644 index 0000000..5cf6f57 --- /dev/null +++ b/python/tests/test_codegen_sens_budget.py @@ -0,0 +1,447 @@ +"""GH #90 — a build-time derivation budget for the sensitivity ∂f/∂p path. + +#95/#187 bounded the symbolic derivation of the analytical *Jacobian*: a model +that does not derive in time falls back to the finite-difference Jacobian instead +of hanging the load. The ∂f/∂p path added by #55 (#65/#66/#67/#68) does the same +kind of sympy work — one ``sp.diff`` per (distinct rate law, parameter it reads) +pair, the one axis here that grows with the *product* of model size and parameter +count — and had no budget at all. The failure mode that removes is a build that +appears to hang rather than one that declines and says why. + +Two things are under test, and the second is why the first is safe: + +* **The policy.** The sensitivity budget shares the Jacobian's base, slope and + override grammar (one spelling, so they cannot drift) but has its own env var + and is **never unbounded by size**. That last difference is the design decision + worth pinning: ``_FD_NONVIABLE_SPECIES`` exists because past that size an FD + Jacobian does not converge, so there is nothing to fall back *to*; the + sensitivity fallback is CVODES' internal difference quotient, which is what + every Functional model used before #55 and stays viable at every scale. + +* **The plumbing.** The deadline is checked *during* the derivation (on entry to + each rate law, before each ``sp.diff``, and inside the derived-parameter chain + rule), not merely before it — otherwise the budget bounds nothing on the single + pathological model it exists for. An expiry declines through the same channel + as every other decline, names the override, and leaves a model that still + solves. + +The one thing deliberately *not* budgeted here is ``_derived_expr_partials_numeric`` +(the #43 IC-seed path). It runs at solve setup rather than on this build, and a +partial result there is a silently wrong ∂x₀/∂p rather than a slower-but-correct +fallback, so it needs its own all-or-nothing story and its own issue. +""" + +from __future__ import annotations + +import logging +import math + +import bngsim +import numpy as np +import pytest +from bngsim import _codegen as cg +from bngsim._jacobian import ( + _BUDGET_PER_SPECIES_S, + _DEFAULT_DERIVATION_BUDGET_S, + _FD_NONVIABLE_SPECIES, + _derivation_budget_s, + _DerivationBudgetExceeded, + _sens_derivation_budget_s, +) + +pytest.importorskip("sympy") + +_SENS_ENV = "BNGSIM_SENS_DERIV_BUDGET_S" +_JAC_ENV = "BNGSIM_JAC_DERIV_BUDGET_S" + + +def _has_cc() -> bool: + try: + cg._find_c_compiler() + return True + except Exception: + return False + + +requires_cc = pytest.mark.skipif(not _has_cc(), reason="no C compiler available") + +# GH #85, open and pre-existing on main: under the MIR JIT backend a Functional +# model constructed *with* ``sensitivity_params`` does not compile — c2mir rejects +# the GH #198 ``bngsim_codegen_output_sens`` block that ``cc`` accepts. Only the +# one test here that builds such a Simulator is blocked; every source-level test +# runs on both backends. Spelled the same way (and retiring the same day) as the +# markers in test_codegen_functional_sens_rhs.py and +# test_codegen_switch_condition_sens.py. +blocked_on_mir_jit = pytest.mark.xfail( + cg._codegen_jit_backend() == "mir", + reason="GH #85: c2mir cannot compile bngsim_codegen_output_sens", + raises=RuntimeError, + strict=True, +) + + +@pytest.fixture(autouse=True) +def _clear_budget_env(monkeypatch): + """Both budgets resolve from their env var first, so every case here starts + from the size-derived default unless it sets one explicitly.""" + monkeypatch.delenv(_SENS_ENV, raising=False) + monkeypatch.delenv(_JAC_ENV, raising=False) + + +# ─── fixtures ────────────────────────────────────────────────────────────── +# +# Same shapes as test_codegen_functional_sens_rhs.py: a Functional law (which +# reaches _functional_rate_law_partials) and an Elementary model with a derived +# rate constant (which reaches only the chain-rule loop, the *other* budgeted +# site and the one an all-Elementary model can still hang on). + +SIR = """\ +begin parameters + 1 S0 2e7 # Constant + 2 I0 1 # Constant + 3 beta 1/S0 # ConstantExpression + 4 gamma 1/7 # Constant +end parameters +begin functions + 1 betaI() beta*I +end functions +begin species + 1 person(state~S) S0 + 2 person(state~I) I0 + 3 person(state~R) 0 +end species +begin reactions + 1 1 2 betaI #_R1 + 2 2 3 gamma #_R2 +end reactions +begin groups + 1 S 1 + 2 I 2 + 3 R 3 +end groups +""" + +# Smooth saturating algebra, and — unlike SIR — a model whose difference-quotient +# fallback actually resolves, so the cost of a decline can be measured against the +# analytic answer rather than against a step-budget failure. +HILL = """\ +begin parameters + 1 kmax 3.5 # Constant + 2 Km 4.0 # Constant + 3 n 2.0 # Constant + 4 kdeg 0.2 # Constant +end parameters +begin functions + 1 vsat() kmax*Atot^n/(Km^n + Atot^n) +end functions +begin species + 1 A() 6.0 + 2 B() 1.0 +end species +begin reactions + 1 1 2 vsat #_R1 + 2 2 1 kdeg #_R2 +end reactions +begin groups + 1 Atot 1 + 2 Btot 2 +end groups +""" + +# ``sin`` puts the ∂f/∂x half (J·yS) on the SymPy fallback: the #151 native +# saturable emitters cover mass-action, Hill and exponential shapes without ever +# importing SymPy, so a fixture from that family would exercise no derivation to +# budget. test_codegen_determinism.py needs ``sin()`` for the same reason. +TRANSCENDENTAL = """\ +begin parameters + 1 kmax 3.5 # Constant + 2 Km 4.0 # Constant + 3 kdeg 0.2 # Constant +end parameters +begin functions + 1 vosc() kmax*sin(Atot/Km) +end functions +begin species + 1 A() 6.0 + 2 B() 1.0 +end species +begin reactions + 1 1 2 vosc #_R1 + 2 2 1 kdeg #_R2 +end reactions +begin groups + 1 Atot 1 + 2 Btot 2 +end groups +""" + +ELEMENTARY = """\ +begin parameters + 1 k1 0.3 # Constant + 2 scale 2.0 # Constant + 3 k2 k1*scale # ConstantExpression +end parameters +begin species + 1 A() 10.0 + 2 B() 0.0 +end species +begin reactions + 1 1 2 k1 #_R1 + 2 2 1 k2 #_R2 +end reactions +begin groups + 1 Atot 1 + 2 Btot 2 +end groups +""" + + +def _model(tmp_path, text, name="m.net"): + net = tmp_path / name + net.write_text(text) + return bngsim.Model.from_net(net) + + +def _net(tmp_path, text, name="m.net"): + net = tmp_path / name + net.write_text(text) + return str(net) + + +# ─── the policy ──────────────────────────────────────────────────────────── + + +class TestPolicy: + @pytest.mark.parametrize("n_species", [0, 1, 295, 1000]) + def test_small_models_share_the_jacobian_base(self, n_species): + """Below the scaling knee the sensitivity budget is the #95 base, exactly + as the Jacobian's is — the two policies are one policy with two dials.""" + assert _sens_derivation_budget_s(n_species=n_species) == _DEFAULT_DERIVATION_BUDGET_S + assert _sens_derivation_budget_s(n_species=n_species) == _derivation_budget_s( + n_species=n_species + ) + + def test_the_budget_scales_with_species_count(self): + n = int(4 * _DEFAULT_DERIVATION_BUDGET_S / _BUDGET_PER_SPECIES_S) + assert n < _FD_NONVIABLE_SPECIES + assert _sens_derivation_budget_s(n_species=n) == pytest.approx(_BUDGET_PER_SPECIES_S * n) + + def test_it_stays_finite_where_the_jacobian_goes_unbounded(self): + """The one substantive difference between the two budgets, and the reason + it exists: past ``_FD_NONVIABLE_SPECIES`` an FD Jacobian does not converge, + so the Jacobian derivation is mandatory and runs to completion. Declining a + sensitivity RHS only hands the columns to CVODES' difference quotient, + which is correct at any size — so there is never a reason to let *this* + derivation run unbounded, and a genome-scale model gets a decline instead + of a build that appears to hang.""" + n = _FD_NONVIABLE_SPECIES * 4 + assert _derivation_budget_s(n_species=n) is None + budget = _sens_derivation_budget_s(n_species=n) + assert budget is not None + assert math.isfinite(budget) + # Generous, though: ~10x the observed derivation rate, so a model that + # scales like a real one is not cut off for being large. + assert budget == pytest.approx(_BUDGET_PER_SPECIES_S * n) + + @pytest.mark.parametrize("raw", ["inf", "none", "off", "0", "-1", " INF ", "nan"]) + def test_the_override_can_disable_it(self, monkeypatch, raw): + monkeypatch.setenv(_SENS_ENV, raw) + assert _sens_derivation_budget_s(n_species=10) is None + + def test_the_override_wins_over_the_size_policy(self, monkeypatch): + monkeypatch.setenv(_SENS_ENV, "2.5") + assert _sens_derivation_budget_s(n_species=_FD_NONVIABLE_SPECIES * 4) == 2.5 + + def test_a_malformed_override_degrades_to_the_policy(self, monkeypatch): + """A typo must not read as "no budget at all" — that is the failure this + whole issue is about.""" + monkeypatch.setenv(_SENS_ENV, "twenty") + assert _sens_derivation_budget_s(n_species=10) == _DEFAULT_DERIVATION_BUDGET_S + + def test_the_two_budgets_are_independent(self, monkeypatch): + """``BNGSIM_JAC_DERIV_BUDGET_S=inf`` is the documented way to keep a + genome-scale model's analytical Jacobian. It buys a different thing and + falls back to a different path, so it must not silently uncap this one — + which would restore exactly the hang #90 removes.""" + monkeypatch.setenv(_JAC_ENV, "inf") + assert _derivation_budget_s(n_species=10) is None + assert _sens_derivation_budget_s(n_species=10) == _DEFAULT_DERIVATION_BUDGET_S + + monkeypatch.setenv(_SENS_ENV, "inf") + monkeypatch.setenv(_JAC_ENV, "3") + assert _derivation_budget_s(n_species=10) == 3.0 + assert _sens_derivation_budget_s(n_species=10) is None + + +# ─── the cache key ───────────────────────────────────────────────────────── + + +class TestCacheKey: + """The ``.net`` path keys its ``.so`` on the model's *content*, not on the + generated C, so anything that changes whether the sens RHS is emitted has to + reach the key — the trap #67's A/B hatch already had to sidestep.""" + + def test_unset_leaves_every_existing_key_untouched(self): + assert cg._sens_budget_cache_tag() == "" + + def test_each_override_gets_its_own_namespace(self, monkeypatch): + monkeypatch.setenv(_SENS_ENV, "1e-9") + tight = cg._sens_budget_cache_tag() + monkeypatch.setenv(_SENS_ENV, "inf") + loose = cg._sens_budget_cache_tag() + assert tight and loose and tight != loose + + def test_the_tag_is_normalized(self, monkeypatch): + monkeypatch.setenv(_SENS_ENV, " INF ") + assert cg._sens_budget_cache_tag() == ":sens_budget=inf" + + @requires_cc + def test_a_tight_budget_is_not_served_a_cached_analytic_so(self, tmp_path, monkeypatch): + """End to end through the real cache: emit under the default budget, then + re-request the same .net under a budget that cannot be met. Sharing a key + would hand back the analytic .so and make the decline invisible.""" + model = _model(tmp_path, SIR) + net = str(tmp_path / "m.net") + analytic = cg.prepare_codegen(net, model=model) + monkeypatch.setenv(_SENS_ENV, "1e-9") + declined = cg.prepare_codegen(net, model=model) + assert analytic != declined + + +# ─── the deadline is checked *during* the derivation ─────────────────────── + + +class TestDeadlinePlumbing: + def test_an_expired_deadline_stops_the_derived_chain_rule(self): + """``_derived_param_jacobian_checked`` raises rather than returning a + reason: a wall-clock expiry is a property of the build, not of this + expression, and the caller memoizes reasons per expression.""" + args = ("k1*scale", {"k1", "scale"}, {"k1": 0, "scale": 1}) + jac, reason = cg._derived_param_jacobian_checked(*args) + assert reason is None and jac + with pytest.raises(_DerivationBudgetExceeded): + cg._derived_param_jacobian_checked(*args, deadline=cg.time.perf_counter() - 1.0) + + def test_no_deadline_is_the_unbudgeted_path(self): + """Every caller outside the sensitivity build passes ``None``, so the + default must be a no-op rather than a budget of zero.""" + cg._check_derivation_deadline(None) # must not raise + + def test_the_jacobian_half_of_the_sens_rhs_is_bounded_too(self, tmp_path): + """``ySdot = J·yS + ∂f/∂p`` has two derivations, and #90's issue text names + only the second. The first is a *re*-derivation of what + ``attach_functional_jacobian`` ran at load — but one that ignores its own + clock, and on a model whose load-time attach was itself cut off by the #95 + budget there is no earlier bound to inherit. Bounding only ∂f/∂p would + leave the same hang reachable one call later.""" + model = _model(tmp_path, TRANSCENDENTAL) + data = model._core.codegen_data() + assert cg._functional_jacobian_groups(model._core, data, cg._jacv_add) is not None + with pytest.raises(_DerivationBudgetExceeded): + cg._functional_jacobian_groups( + model._core, data, cg._jacv_add, cg.time.perf_counter() - 1.0 + ) + + def test_the_budget_is_enforced_mid_derivation(self, tmp_path, monkeypatch): + """The point of checking per rate law and per parameter rather than once + up front. A clock that advances a second per reading expires a 2 s budget + only *after* the derivation has started, so a build that declines here + cannot have been stopped by a pre-flight check.""" + model = _model(tmp_path, SIR) + ticks = iter(range(0, 10_000)) + monkeypatch.setattr(cg.time, "perf_counter", lambda: float(next(ticks))) + monkeypatch.setenv(_SENS_ENV, "2") + assert cg.generate_sens_from_model(model, functional=True) is None + + def test_a_generous_budget_on_the_same_clock_still_emits(self, tmp_path, monkeypatch): + """The control for the case above: same fake clock, same model, a budget + the derivation fits inside. Without this, the test above would also pass + if the deadline were simply always expired.""" + model = _model(tmp_path, SIR) + ticks = iter(range(0, 10_000)) + monkeypatch.setattr(cg.time, "perf_counter", lambda: float(next(ticks))) + monkeypatch.setenv(_SENS_ENV, "1e6") + assert cg.generate_sens_from_model(model, functional=True) is not None + + +# ─── what a caller sees ──────────────────────────────────────────────────── + + +class TestDecline: + def test_the_default_budget_changes_nothing(self, tmp_path): + """The regression guard for every model that already had an analytic sens + RHS: the corpus derives in 2.6 s total, so nothing real is near the + budget and emission must be untouched.""" + _src, has_sens = cg.generate_combined_from_model(_model(tmp_path, SIR)) + assert has_sens is True + + def test_a_functional_model_declines_when_the_budget_expires(self, tmp_path, monkeypatch): + monkeypatch.setenv(_SENS_ENV, "1e-9") + _src, has_sens = cg.generate_combined_from_model(_model(tmp_path, SIR)) + assert has_sens is False + + def test_an_elementary_model_declines_too(self, tmp_path, monkeypatch): + """The chain-rule loop is the *only* sympy an all-Elementary model runs on + this path, and a model with thousands of ``# ConstantExpression`` rate + constants can hang on it without ever touching a Functional law.""" + model = _model(tmp_path, ELEMENTARY) + assert cg.generate_sens_from_model(model) is not None + monkeypatch.setenv(_SENS_ENV, "1e-9") + assert cg.generate_sens_from_model(model) is None + + def test_the_net_text_path_declines_too(self, tmp_path, monkeypatch): + """``generate_sens_rhs_c`` reads the .net as text and never sees a model, + so it is a separate entry point with its own build — and its own + derived-rate-constant derivation to bound.""" + net = _net(tmp_path, ELEMENTARY) + assert cg.generate_sens_rhs_c(net) is not None + monkeypatch.setenv(_SENS_ENV, "1e-9") + assert cg.generate_sens_rhs_c(net) is None + + def test_the_escape_hatch_restores_the_unbudgeted_path(self, tmp_path, monkeypatch): + monkeypatch.setenv(_SENS_ENV, "inf") + _src, has_sens = cg.generate_combined_from_model(_model(tmp_path, SIR)) + assert has_sens is True + + def test_the_decline_is_loud_and_names_the_override(self, tmp_path, monkeypatch, caplog): + """A silent decline is a 9-37x slowdown nobody can attribute. It routes + through the same warning as every other decline reason (#56's precedent), + says how far the derivation got, and names the knob.""" + monkeypatch.setenv(_SENS_ENV, "1e-9") + with caplog.at_level(logging.WARNING, logger="bngsim"): + cg.generate_combined_from_model(_model(tmp_path, SIR)) + msg = "\n".join(r.getMessage() for r in caplog.records) + assert _SENS_ENV in msg + assert "∂f/∂p" in msg and "budget" in msg + assert "difference quotient" in msg + + +class TestFallbackStillSolves: + """What a decline actually costs the caller: the columns are still computed, + by CVODES' internal difference quotient, and they agree with the analytic + ones. + + Not a claim that the fallback is *cheap*. #55 measured 9-37x per column, and + on a stiff model at a tight tolerance the DQ's ~sqrt(rtol) accuracy collapses + the step size outright — the SIR fixture above burns its entire step budget + and returns zeros, which is why it is not the model used here. That is the + honest cost of this budget, and the reason the warning names the override: + declining beats hanging, but it is not free. + """ + + @blocked_on_mir_jit + def test_a_declined_model_still_returns_sensitivities(self, tmp_path, monkeypatch): + analytic = np.asarray(_run(_model(tmp_path, HILL), ["kmax", "Km"]).sensitivities) + + monkeypatch.setenv(_SENS_ENV, "1e-9") + declined = np.asarray(_run(_model(tmp_path, HILL, "m2.net"), ["kmax", "Km"]).sensitivities) + + assert np.all(np.isfinite(declined)) + assert np.any(declined != 0.0) + scale = max(float(np.max(np.abs(analytic))), 1e-300) + np.testing.assert_allclose(declined, analytic, rtol=1e-4, atol=1e-5 * scale) + + +def _run(model, params, t_end=40.0, n=21): + sim = bngsim.Simulator(model, method="ode", sensitivity_params=list(params)) + return sim.run(t_span=(0.0, t_end), n_points=n)