From 5cae9532b648829159b73107eb8e60e0c47f3162 Mon Sep 17 00:00:00 2001 From: Arpit Jain Date: Sat, 29 Aug 2026 14:46:21 -0400 Subject: [PATCH 1/2] Evaluate BNGL parameter expressions instead of dropping them A BNGL parameters block may define a parameter as an expression over other parameters, which is ordinary style rather than an edge case: across the model corpora on hand, 1839 of 8822 parameter declarations (20.8%) are expression-valued, in 156 of 297 files. BnglModel could not evaluate any of them. get_parameter_value raised NotImplementedError, and get_free_parameter_ids_with_values swallowed the ValueError and skipped the parameter, so building a PEtab problem lost parameters with nothing said to the user. Resolving one needs no BNG2.pl and no network generation, because a parameters block is arithmetic over other parameters. Add a small evaluator that walks the definitions in dependency order, so declaration order does not matter, and reports a circular definition by naming the cycle rather than recursing until the stack gives out. The sublanguage is BNGL's, so the expressions are tokenized and parsed rather than passed to eval, which would quietly import Python's meanings: ^ is exponentiation and not bitwise xor, ln is the natural logarithm while log10 and log2 are separate names, and division is floating point. A bare log is rejected rather than assumed to mean ln, so a typo is an error rather than a plausible wrong number. _bngl_expr is stdlib-only and imports nothing from the rest of PyBNF, so it can travel with the adapter to libpetab-python (#591, #420 Step B). This supersedes the confinement the previous behaviour pinned, so the test that asserted the NotImplementedError now asserts the computed value. Fixes #666 Signed-off-by: Arpit Jain --- pybnf/petab/_bngl_expr.py | 296 ++++++++++++++++++++++++++++++++++ pybnf/petab/bngl_model.py | 43 +++-- tests/test_petab_bngl_expr.py | 145 +++++++++++++++++ tests/test_petab_export.py | 12 +- 4 files changed, 476 insertions(+), 20 deletions(-) create mode 100644 pybnf/petab/_bngl_expr.py create mode 100644 tests/test_petab_bngl_expr.py diff --git a/pybnf/petab/_bngl_expr.py b/pybnf/petab/_bngl_expr.py new file mode 100644 index 000000000..16c361dd8 --- /dev/null +++ b/pybnf/petab/_bngl_expr.py @@ -0,0 +1,296 @@ +"""Evaluator for the BNGL parameters-block expression sublanguage. + +A BNGL ``parameters`` block may give a parameter an expression over other +parameters rather than a literal (``kon koff/(Kd*NA*V)``), which is ordinary +style rather than an edge case. Network generation is not needed to resolve +one: a parameters block is arithmetic over other parameters, so the values can +be computed by walking the definitions in dependency order. + +The sublanguage is BNGL's, not Python's, and the two disagree in ways that are +silent rather than loud: + +* ``^`` raises to a power. In Python it is bitwise exclusive-or, so passing an + expression through unchanged would compute a different number rather than + fail. +* The natural logarithm is ``ln``. Python's ``log`` is the natural logarithm + too, but BNGL also has ``log10`` and ``log2``, so the names cannot be mapped + across by position. +* Division is floating point throughout. + +Expressions are therefore tokenized and parsed here rather than handed to +``eval``, which would import Python's precedence and operator meanings along +with the obvious injection problem. + +This module is deliberately self-contained: stdlib only, and no imports from +the rest of PyBNF, so that it can move to ``libpetab-python`` alongside the +BNGL model adapter it serves (see #591, #420 Step B). +""" + +from __future__ import annotations + +import math +import re + +__all__ = [ + "BnglExpressionError", + "CircularParameterError", + "evaluate_expression", + "evaluate_parameters", +] + + +class BnglExpressionError(ValueError): + """A parameter expression could not be parsed or evaluated.""" + + +class CircularParameterError(BnglExpressionError): + """A parameter's definition depends on itself, directly or through others.""" + + +# BNGL's named constants, as BNG2.pl's Expression.pm exposes them. +_CONSTANTS = { + "_pi": math.pi, + "_e": math.e, +} + +# Functions BNGL accepts in a parameter expression. ``ln`` is the natural +# logarithm; ``log`` is deliberately absent, because BNG2.pl does not accept a +# bare ``log`` and silently treating it as ``ln`` would let a typo produce a +# plausible wrong number. +_FUNCTIONS = { + "exp": math.exp, + "ln": math.log, + "log10": math.log10, + "log2": math.log2, + "sqrt": math.sqrt, + "abs": abs, + "sin": math.sin, + "cos": math.cos, + "tan": math.tan, + "asin": math.asin, + "acos": math.acos, + "atan": math.atan, + "sinh": math.sinh, + "cosh": math.cosh, + "tanh": math.tanh, + "floor": math.floor, + "ceil": math.ceil, + "rint": lambda x: float(round(x)), + "min": min, + "max": max, +} + +_TOKEN_RE = re.compile( + r""" + (?P\d+\.\d*(?:[eE][+-]?\d+)? + |\.\d+(?:[eE][+-]?\d+)? + |\d+(?:[eE][+-]?\d+)?) + | (?P[A-Za-z_]\w*) + | (?P\*\*|[-+*/^(),]) + | (?P\s+) + """, + re.VERBOSE, +) + + +def _tokenize(text: str) -> list[tuple[str, str]]: + tokens: list[tuple[str, str]] = [] + pos = 0 + while pos < len(text): + match = _TOKEN_RE.match(text, pos) + if match is None: + raise BnglExpressionError( + f"Unexpected character {text[pos]!r} at position {pos} in {text!r}" + ) + pos = match.end() + kind = match.lastgroup + if kind == "space": + continue + value = match.group() + # BNGL writes exponentiation as ^; accept ** as well, since BNG2.pl does. + tokens.append(("op", "^") if value == "**" else (kind, value)) + return tokens + + +class _Parser: + """Recursive-descent parser for the arithmetic sublanguage. + + Precedence, loosest to tightest: ``+ -``, then ``* /``, then unary ``-``, + then ``^`` (right associative). ``^`` binding tighter than unary minus is + what makes ``-2^2`` come out as ``-4``. + """ + + def __init__(self, tokens: list[tuple[str, str]], text: str, lookup): + self._tokens = tokens + self._text = text + self._lookup = lookup + self._pos = 0 + + def parse(self) -> float: + value = self._parse_sum() + if self._pos != len(self._tokens): + raise BnglExpressionError( + f"Unexpected trailing input in {self._text!r} at token " + f"{self._tokens[self._pos][1]!r}" + ) + return value + + def _peek(self) -> tuple[str, str] | None: + return self._tokens[self._pos] if self._pos < len(self._tokens) else None + + def _accept(self, value: str) -> bool: + token = self._peek() + if token is not None and token[0] == "op" and token[1] == value: + self._pos += 1 + return True + return False + + def _expect(self, value: str) -> None: + if not self._accept(value): + found = self._peek() + raise BnglExpressionError( + f"Expected {value!r} in {self._text!r}, found " + + (repr(found[1]) if found else "end of expression") + ) + + def _parse_sum(self) -> float: + value = self._parse_product() + while True: + if self._accept("+"): + value += self._parse_product() + elif self._accept("-"): + value -= self._parse_product() + else: + return value + + def _parse_product(self) -> float: + value = self._parse_unary() + while True: + if self._accept("*"): + value *= self._parse_unary() + elif self._accept("/"): + divisor = self._parse_unary() + if divisor == 0: + raise BnglExpressionError(f"Division by zero in {self._text!r}") + # True division throughout: BNGL has no integer division, and + # Python's / on two ints would still be float, but being + # explicit keeps that from depending on operand types. + value = float(value) / float(divisor) + else: + return value + + def _parse_unary(self) -> float: + if self._accept("-"): + return -self._parse_unary() + if self._accept("+"): + return self._parse_unary() + return self._parse_power() + + def _parse_power(self) -> float: + base = self._parse_atom() + if self._accept("^"): + # Right associative, and the exponent may itself be signed. + return base ** self._parse_unary() + return base + + def _parse_atom(self) -> float: + token = self._peek() + if token is None: + raise BnglExpressionError(f"Expression ended unexpectedly: {self._text!r}") + kind, value = token + + if kind == "number": + self._pos += 1 + return float(value) + + if kind == "op" and value == "(": + self._pos += 1 + inner = self._parse_sum() + self._expect(")") + return inner + + if kind == "name": + self._pos += 1 + if self._accept("("): + args = [self._parse_sum()] + while self._accept(","): + args.append(self._parse_sum()) + self._expect(")") + return self._call(value, args) + if value in _CONSTANTS: + return _CONSTANTS[value] + return self._lookup(value) + + raise BnglExpressionError(f"Unexpected token {value!r} in {self._text!r}") + + def _call(self, name: str, args: list[float]) -> float: + try: + func = _FUNCTIONS[name] + except KeyError: + raise BnglExpressionError( + f"Unknown function {name!r} in {self._text!r}" + ) from None + try: + return float(func(*args)) + except TypeError as e: + raise BnglExpressionError( + f"Wrong number of arguments to {name!r} in {self._text!r}" + ) from e + except ValueError as e: + raise BnglExpressionError( + f"{name}() is undefined for its argument in {self._text!r}: {e}" + ) from e + + +def evaluate_expression(text: str, symbols: dict[str, float]) -> float: + """Evaluate one BNGL expression against already-resolved ``symbols``.""" + return _Parser(_tokenize(text), text, lambda n: _resolve_known(n, symbols, text)).parse() + + +def _resolve_known(name: str, symbols: dict[str, float], text: str) -> float: + try: + return symbols[name] + except KeyError: + raise BnglExpressionError( + f"Unknown parameter {name!r} in {text!r}" + ) from None + + +def evaluate_parameters(parameters: dict[str, str]) -> dict[str, float]: + """Resolve a BNGL parameters block to numbers. + + ``parameters`` maps a parameter name to its raw right-hand side, literal or + expression, as :func:`pybnf.petab._bngl.parse_model` collects it. Values are + resolved lazily in dependency order, so declaration order does not matter, + which matches BNG2.pl. + + Raises :class:`CircularParameterError` on a definition that depends on + itself, and :class:`BnglExpressionError` on anything unparseable or on a + reference to a name the block does not define. + """ + resolved: dict[str, float] = {} + resolving: list[str] = [] + + def lookup(name: str) -> float: + if name in resolved: + return resolved[name] + if name in resolving: + cycle = " -> ".join([*resolving[resolving.index(name):], name]) + raise CircularParameterError( + f"Parameter {name!r} is defined in terms of itself: {cycle}" + ) + if name not in parameters: + raise BnglExpressionError(f"Unknown parameter {name!r}") + resolving.append(name) + try: + value = _Parser( + _tokenize(parameters[name]), parameters[name], lookup + ).parse() + finally: + resolving.pop() + resolved[name] = value + return value + + for name in parameters: + lookup(name) + return resolved diff --git a/pybnf/petab/bngl_model.py b/pybnf/petab/bngl_model.py index 843cb1389..797e9a3bb 100644 --- a/pybnf/petab/bngl_model.py +++ b/pybnf/petab/bngl_model.py @@ -29,6 +29,7 @@ from petab.v1.models.model import Model from ._bngl import parse_model +from ._bngl_expr import BnglExpressionError, evaluate_parameters #: BNGL model type, as used in a PEtab v2 yaml file as ``language``. MODEL_TYPE_BNGL = 'bngl' @@ -44,6 +45,7 @@ def __init__(self, entities, model_id, path=None): self._entities = entities self._model_id = model_id self._path = Path(path) if path is not None else None + self._resolved_parameters = None @staticmethod def from_file(filepath_or_buffer, model_id=None, base_path=None): @@ -68,28 +70,35 @@ def model_id(self): def get_parameter_ids(self): return list(self._entities.parameters) + def _parameter_values(self): + """Every parameter resolved to a number, computed once and cached. + + A parameters block is arithmetic over other parameters, so this needs + no BNG2.pl and no network generation; see :mod:`pybnf.petab._bngl_expr`. + """ + if self._resolved_parameters is None: + self._resolved_parameters = evaluate_parameters( + dict(self._entities.parameters) + ) + return self._resolved_parameters + def get_parameter_value(self, id_): + if id_ not in self._entities.parameters: + raise ValueError(f"Parameter {id_} does not exist.") try: - rhs = self._entities.parameters[id_] - except KeyError as e: - raise ValueError(f"Parameter {id_} does not exist.") from e - try: - return float(rhs) - except ValueError as e: - raise NotImplementedError( - f"Parameter '{id_}' has an expression value '{rhs}'. Evaluating a " - f"BNGL parameter expression needs BNG2.pl/network generation, which " - f"is out of scope for the validation-grade BnglModel (ADR-0026)." + return self._parameter_values()[id_] + except BnglExpressionError as e: + raise ValueError( + f"Parameter '{id_}' has an expression value " + f"'{self._entities.parameters[id_]}' that could not be evaluated: {e}" ) from e def get_free_parameter_ids_with_values(self): - out = [] - for name, rhs in self._entities.parameters.items(): - try: - out.append((name, float(rhs))) - except ValueError: - continue # an expression-valued parameter has no validation-grade value - return out + # Expression-valued parameters used to be skipped here, which lost them + # from the PEtab problem with nothing said. They are resolved now; a + # block that still cannot be evaluated raises rather than going quiet. + values = self._parameter_values() + return [(name, values[name]) for name in self._entities.parameters] def get_valid_parameters_for_parameter_table(self): return list(self._entities.parameters) diff --git a/tests/test_petab_bngl_expr.py b/tests/test_petab_bngl_expr.py new file mode 100644 index 000000000..d3e32dcf9 --- /dev/null +++ b/tests/test_petab_bngl_expr.py @@ -0,0 +1,145 @@ +"""BNGL parameter-expression evaluation (issue #666).""" + +import math + +import pytest + +from pybnf.petab._bngl import parse_model +from pybnf.petab._bngl_expr import ( + BnglExpressionError, + CircularParameterError, + evaluate_expression, + evaluate_parameters, +) +from pybnf.petab.bngl_model import BnglModel + + +def test_expression_valued_parameter_is_resolved(): + """The case from issue #666: kon is an expression over other parameters.""" + text = """ +begin model +begin parameters + NA 6.022e23 + V 1e-12 + Kd 5.0 + koff 0.1 + kon koff/(Kd*NA*V) +end parameters +end model +""" + m = BnglModel(parse_model(text), model_id='demo') + + assert m.get_parameter_value('kon') == pytest.approx(0.1 / (5.0 * 6.022e23 * 1e-12)) + + # The parameter used to be dropped from this list entirely. + ids = [name for name, _ in m.get_free_parameter_ids_with_values()] + assert ids == list(m.get_parameter_ids()) + assert 'kon' in ids + + +def test_caret_is_exponentiation_not_xor(): + """BNGL's ^ raises to a power; Python's is bitwise xor, and 2^3 there is 1.""" + assert evaluate_parameters({'a': '2^3'})['a'] == 8.0 + + +def test_exponentiation_binds_tighter_than_unary_minus(): + assert evaluate_parameters({'a': '-2^2'})['a'] == -4.0 + + +def test_exponentiation_is_right_associative(): + assert evaluate_parameters({'a': '2^3^2'})['a'] == 512.0 + + +def test_ln_is_the_natural_logarithm_and_bare_log_is_rejected(): + assert evaluate_parameters({'a': 'ln(_e)'})['a'] == pytest.approx(1.0) + assert evaluate_parameters({'a': 'log10(1000)'})['a'] == pytest.approx(3.0) + + # BNG2.pl has no bare log(); treating it as ln would turn a typo into a + # plausible wrong number rather than an error. + with pytest.raises(BnglExpressionError, match='log'): + evaluate_parameters({'a': 'log(10)'}) + + +def test_division_is_floating_point(): + assert evaluate_parameters({'a': '1/2'})['a'] == 0.5 + + +def test_declaration_order_does_not_matter(): + """BNG2.pl resolves by dependency, not by position in the block.""" + assert evaluate_parameters({'b': 'a*2', 'a': '3'}) == {'a': 3.0, 'b': 6.0} + + +@pytest.mark.parametrize( + 'params, target, expected', + [ + # Shapes taken from the survey in issue #666. + ({'kp18': '2', 'km18': '1', 'kp19': '3', 'km19': '1', 'kp22': '4', + 'km22': '2', 'kp20': '5', 'km20': '1', + 'loop3': '(kp18/km18)*(kp19/km19)/((kp22/km22)*(kp20/km20))'}, + 'loop3', (2 / 1) * (3 / 1) / ((4 / 2) * (5 / 1))), + ({'p_RM_AC': '7', 'p_RM_A': 'p_RM_AC'}, 'p_RM_A', 7.0), + ({'lifetime': '4', 'gamma_R': '1/lifetime'}, 'gamma_R', 0.25), + ({'krZapTcr': '3', 'krZapCd3e': '10*krZapTcr'}, 'krZapCd3e', 30.0), + ({'Kd_BRAF_RAFi2': '20', 'Gf_BRAF_RAFi2': 'ln(Kd_BRAF_RAFi2)'}, + 'Gf_BRAF_RAFi2', math.log(20)), + ], +) +def test_real_world_expression_shapes(params, target, expected): + assert evaluate_parameters(params)[target] == pytest.approx(expected) + + +def test_chained_expression_dependencies_resolve(): + """A parameter may depend on another that is itself an expression.""" + values = evaluate_parameters({'a': '2', 'b': 'a*3', 'c': 'b+a'}) + assert values == {'a': 2.0, 'b': 6.0, 'c': 8.0} + + +def test_circular_definition_names_the_cycle(): + with pytest.raises(CircularParameterError) as excinfo: + evaluate_parameters({'a': 'b', 'b': 'a'}) + assert 'a -> b -> a' in str(excinfo.value) + + +def test_self_referential_definition_is_reported(): + with pytest.raises(CircularParameterError): + evaluate_parameters({'a': 'a+1'}) + + +@pytest.mark.parametrize( + 'rhs', + ['b', '2 +', 'foo(1)', '1/0', '2 @ 3'], +) +def test_unusable_expressions_raise_rather_than_go_quiet(rhs): + with pytest.raises(BnglExpressionError): + evaluate_parameters({'a': rhs}) + + +def test_unevaluable_parameter_surfaces_from_the_model(): + """The adapter reports the failure instead of dropping the parameter.""" + text = """ +begin model +begin parameters + a b +end parameters +end model +""" + m = BnglModel(parse_model(text), model_id='demo') + with pytest.raises(ValueError, match='could not be evaluated'): + m.get_parameter_value('a') + + +def test_missing_parameter_still_raises_value_error(): + text = """ +begin model +begin parameters + a 1 +end parameters +end model +""" + m = BnglModel(parse_model(text), model_id='demo') + with pytest.raises(ValueError, match='does not exist'): + m.get_parameter_value('nope') + + +def test_evaluate_expression_against_known_symbols(): + assert evaluate_expression('x*2 + y', {'x': 1.5, 'y': 1.0}) == 4.0 diff --git a/tests/test_petab_export.py b/tests/test_petab_export.py index a96407de4..03f8e6ca1 100644 --- a/tests/test_petab_export.py +++ b/tests/test_petab_export.py @@ -2141,7 +2141,11 @@ def test_is_state_variable_is_seed_species_only(self, model): assert not model.is_state_variable('v1') # a parameter is not a species assert not model.is_state_variable('x') # nor is an observable - def test_expression_valued_parameter_is_not_evaluated(self): + def test_expression_valued_parameter_is_evaluated(self): + # Superseded #666: an expression RHS used to raise NotImplementedError, + # and get_free_parameter_ids_with_values dropped the parameter without + # saying so. A parameters block is arithmetic over other parameters, so + # it is resolved here without BNG2.pl; see pybnf.petab._bngl_expr. pytest.importorskip('petab') from pybnf.petab._bngl import parse_model from pybnf.petab.bngl_model import BnglModel @@ -2149,8 +2153,10 @@ def test_expression_valued_parameter_is_not_evaluated(self): "begin parameters\n base 2\n k_on 2*base\nend parameters\n") model = BnglModel(ent, model_id='m') assert model.get_parameter_value('base') == 2.0 # numeric RHS -> float - with pytest.raises(NotImplementedError): # expression RHS -> confined - model.get_parameter_value('k_on') + assert model.get_parameter_value('k_on') == 4.0 # expression RHS -> resolved + assert dict(model.get_free_parameter_ids_with_values()) == { + 'base': 2.0, 'k_on': 4.0, + } # -- is_valid: both contract paths pinned (#437) -------------------------- # The contract (ADR-0026): shell to `BNG2.pl --check` when a BNG2.pl is From 2b59d5edc4fb1cd8456374f9495495a562a7ed84 Mon Sep 17 00:00:00 2001 From: Bill Hlavacek Date: Sat, 29 Aug 2026 15:41:33 -0600 Subject: [PATCH 2/2] Pin the BNGL expression evaluator to BNG2.pl, and stop one bad parameter costing the block The evaluator's semantics were derived from the issue text rather than from a run, and BNG2.pl 2.9.3 disagrees in three ways that produce a wrong number silently -- the failure class #666 set out to fix. Checked by putting each expression through writeNET({evaluate_expressions=>1}), the only export path that emits numbers instead of echoing the source: -2^2 BNG2.pl 4 was -4 unary minus binds tighter than ^ 2^3^2 BNG2.pl 64 was 512 ^ is left associative, not right rint(2.5) BNG2.pl 3 was 2 rint is floor(x+0.5), not round-half-even The tests asserted the first two, so they would have kept anyone from noticing. Also aligns the built-in table with Expression.pm: _pi/_e are zero-argument functions (_pi(), not _pi); floor/ceil are not BNGL (commented out upstream as unsupported by muParser); asinh/acosh/atanh/sum/avg were missing; and the comparison and logical operators plus if() were unlexable, which is real BNGL that BNG2.pl evaluates. TFUN is deliberately still absent -- it reads a data file at simulation time, so it is not a parameter-block constant. Resolution is now partial. get_free_parameter_ids_with_values() raised for the whole block on a single unusable definition, which lost more than the original bug did: antigen_pulses_harmon2017_simplified.bngl went from 9 parameters to an exception. It now returns everything it can resolve and warns, naming what it could not. Across the 303 corpus models this moves parameter coverage from 79.3% to 99.0% (+1845), with no model raising; the 89 residuals are all PyBNF __FREE placeholders, which correctly have no value until a fit substitutes one. BNG_VERIFIED is the contract, checked from both sides: pinned offline so it holds in ordinary CI, and re-derived from a real BNG2.pl in one run under the existing `bionetgen` marker so it cannot quietly rot. 84 of 86 differential probes agree; the two that do not are documented and both in the permissive direction (BNG2.pl drops a forward-referencing parameter, and accepts a trailing operator). ADR-0026 is annotated where it says an expression RHS is out of scope. --- ...del-language-runtime-registered-adapter.md | 18 +- pybnf/petab/_bngl_expr.py | 229 +++++++++++--- pybnf/petab/bngl_model.py | 45 ++- tests/test_petab_bngl_expr.py | 291 +++++++++++++++--- tests/test_petab_export.py | 5 +- 5 files changed, 488 insertions(+), 100 deletions(-) diff --git a/docs/adr/0026-bngl-is-a-first-class-petab-model-language-runtime-registered-adapter.md b/docs/adr/0026-bngl-is-a-first-class-petab-model-language-runtime-registered-adapter.md index 3786f6d25..199a785e3 100644 --- a/docs/adr/0026-bngl-is-a-first-class-petab-model-language-runtime-registered-adapter.md +++ b/docs/adr/0026-bngl-is-a-first-class-petab-model-language-runtime-registered-adapter.md @@ -65,6 +65,13 @@ in BNGL: | `get_parameter_ids()` | `parameters` names | `v1, v2, v3` | — | | `get_parameter_value(id)` | numeric RHS → `float`; **expression RHS → `NotImplementedError`**; unknown id → `ValueError` | `v1=5` | — | | `get_free_parameter_ids_with_values()` | `parameters` with numeric RHS | `(v1,5),(v2,5),(v3,5)` | — | + +> **Superseded by issue #666 (PR #673):** both rows now *evaluate* an expression RHS. +> A parameters block is arithmetic over other parameters, so resolving it needs no +> network generation — see `pybnf.petab._bngl_expr`. `get_parameter_value` returns the +> number (or a `ValueError` naming why it could not be computed), and +> `get_free_parameter_ids_with_values()` returns every parameter it can resolve, +> warning about the rest instead of dropping them in silence. | `has_entity_with_id(id)` | **params ∪ observables ∪ global functions ∪ molecule types ∪ compartments ∪ seed species** | `True` for `v1`/`x`/`y`/`counter`; `False` for `obs_x` | `CheckObservablesDoNotShadowModelEntities`, `CheckValidParameterInConditionOrParameterTable` | | `get_valid_parameters_for_parameter_table()` | `parameters` names | `v1, v2, v3` | `CheckAllParametersPresentInParameterTable` (`allowed`) | | `get_valid_ids_for_condition_table()` | parameters ∪ compartments | `v1, v2, v3` | `CheckValidConditionTargets`, `CheckValidParameterInConditionOrParameterTable` | @@ -144,6 +151,8 @@ Grammar points the enrichment must honor (verified against the EBNF + BNG2.pl): bare-number RHS, raises `NotImplementedError` for an **expression RHS** (`k_on 2*base_rate` — evaluating an expression tree is the simulation-grade work scoped out), and `ValueError` for an unknown id (the ABC's documented exception). + **Superseded by issue #666 (PR #673):** the expression RHS is evaluated now; of this + sentence only the unknown-id `ValueError` survives. - **Observables** are `("Molecules"|"Species"|"Counter") WS Name WS Pattern` — name is the *second* token, guarded by a leading observable keyword. - **Global functions** are `Name "(" [args] ")" (WS|"=") MathExpression` — `y()=…`, `f(x)=…`. @@ -213,7 +222,8 @@ environments and de-risked on the acceptance artifact: `obs_x`/unknowns; `symbol_allowed_in_observable_formula` True for `x`,`y`,`v1` and False for `obs_x`/`func_y`/compartments/unknowns; `is_state_variable` True for the seed species, False for `v1`; `get_parameter_value` returns the nominal, raises `ValueError` on unknown and - `NotImplementedError` on an expression-valued parameter. + `NotImplementedError` on an expression-valued parameter + (**superseded by issue #666 (PR #673):** the expression is evaluated instead). - **Unit-test `register_bngl()`** idempotency (two calls; `sbml`/`pysb` still route to the originals) and that the `_bngl.parse_model` refactor leaves the exporter's 22 tests green. - **Make the oracle run in CI**: add `petab` to the `setup-pybnf` composite action so @@ -235,6 +245,12 @@ functions, seed species, optional compartments; the demo (`parabola.bngl`) and i - **Expression-valued parameter *values*** (`k_on 2*base_rate`): the id is enumerated (it is a model entity), but `get_parameter_value` raises `NotImplementedError` rather than evaluating an expression tree. Confined, not silent. + **Superseded by issue #666 (PR #673).** The confinement was the wrong call, and it was not + really confined: `get_free_parameter_ids_with_values()` *dropped* the parameter in silence, + and 20.8% of parameter declarations across our model corpora (1934 of 9323) are + expression-valued. A parameters block needs no network generation to resolve, so this is + now in scope; the BNGL semantics are pinned against a real BNG2.pl in + `tests/test_petab_bngl_expr.py`. - **The full generated-species list** (network gen). `is_state_variable` answers at the seed-species grain validation needs; it does not enumerate the reaction network. For the demo (no conditions) it is never consequential. Revisited when the conditions/experiments export diff --git a/pybnf/petab/_bngl_expr.py b/pybnf/petab/_bngl_expr.py index 16c361dd8..57ddfcb2a 100644 --- a/pybnf/petab/_bngl_expr.py +++ b/pybnf/petab/_bngl_expr.py @@ -7,14 +7,30 @@ be computed by walking the definitions in dependency order. The sublanguage is BNGL's, not Python's, and the two disagree in ways that are -silent rather than loud: +silent rather than loud. Every rule below was checked against BNG2.pl 2.9.3 by +running the expression through ``writeNET({evaluate_expressions=>1})``, which is +the only export path that emits *numbers* instead of echoing the source text; +the function table and precedence order are ``Perl2/Expression.pm`` (``%functions`` +at l.53, ``%NARGS`` at l.245, and the precedence list in ``arrayToExpression`` +at l.2036): * ``^`` raises to a power. In Python it is bitwise exclusive-or, so passing an expression through unchanged would compute a different number rather than fail. +* ``^`` is **left** associative: ``2^3^2`` is 64, not 512. +* Unary minus binds **tighter** than ``^``, so ``-2^2`` is ``(-2)^2`` == 4, not + ``-(2^2)`` == -4. This holds uniformly for literals, parameters, parenthesised + groups and function calls (``-exp(0)^2`` == 1). * The natural logarithm is ``ln``. Python's ``log`` is the natural logarithm too, but BNGL also has ``log10`` and ``log2``, so the names cannot be mapped - across by position. + across by position. A bare ``log`` is rejected, as BNG2.pl rejects it. +* ``rint`` is ``floor(x + 0.5)`` -- round half *up*, not Python's round-half-to- + even. ``rint(2.5)`` is 3 and ``rint(0.5)`` is 1. +* ``_pi`` and ``_e`` are zero-argument *functions*, written ``_pi()``/``_e()``. + Bare ``_pi`` is not a name BNG2.pl resolves. +* Comparison and logical operators yield 1.0/0.0, and ``if(cond, a, b)`` selects + on ``cond != 0``. BNG2.pl evaluates all three arguments before selecting, so + ``if(1, 5, 1/0)`` is an error there and here. * Division is floating point throughout. Expressions are therefore tokenized and parsed here rather than handed to @@ -36,6 +52,7 @@ "CircularParameterError", "evaluate_expression", "evaluate_parameters", + "evaluate_parameters_partial", ] @@ -47,23 +64,29 @@ class CircularParameterError(BnglExpressionError): """A parameter's definition depends on itself, directly or through others.""" -# BNGL's named constants, as BNG2.pl's Expression.pm exposes them. -_CONSTANTS = { - "_pi": math.pi, - "_e": math.e, -} +def _if(cond, then_, else_): + # BNG2.pl's built-in is a plain Perl sub, so all three arguments are already + # evaluated by the time it chooses; it does not short-circuit. Taking floats + # here reproduces that -- `if(1, 5, 1/0)` fails in both. + return then_ if cond != 0 else else_ + -# Functions BNGL accepts in a parameter expression. ``ln`` is the natural -# logarithm; ``log`` is deliberately absent, because BNG2.pl does not accept a -# bare ``log`` and silently treating it as ``ln`` would let a typo produce a -# plausible wrong number. +# The built-in functions BNG2.pl accepts, mirroring %functions in Expression.pm. +# `log` is absent because BNG2.pl has no bare `log` and silently treating it as +# `ln` would let a typo produce a plausible wrong number. `floor` and `ceil` are +# absent because Expression.pm keeps them commented out ("not supported by +# muParser"); BNG2.pl rejects them. `TFUN` is deliberately not implemented: it +# reads a data file at simulation time, so it is not a parameter-block constant. _FUNCTIONS = { + "_pi": lambda: math.pi, + "_e": lambda: math.e, "exp": math.exp, "ln": math.log, "log10": math.log10, "log2": math.log2, "sqrt": math.sqrt, "abs": abs, + "rint": lambda x: float(math.floor(x + 0.5)), "sin": math.sin, "cos": math.cos, "tan": math.tan, @@ -73,25 +96,44 @@ class CircularParameterError(BnglExpressionError): "sinh": math.sinh, "cosh": math.cosh, "tanh": math.tanh, - "floor": math.floor, - "ceil": math.ceil, - "rint": lambda x: float(round(x)), + "asinh": math.asinh, + "acosh": math.acosh, + "atanh": math.atanh, + "if": _if, "min": min, "max": max, + "sum": lambda *a: math.fsum(a), + "avg": lambda *a: math.fsum(a) / len(a), } +#: Names BNG2.pl refuses to accept as a parameter name ("Cannot use built-in +#: function name '_pi' as a parameter"). +RESERVED_NAMES = frozenset(_FUNCTIONS) + +# Longest-first, so `**`, `>=`, `&&` and friends are not split into single +# characters. `~=` is BNG2.pl's alias for `!=`. _TOKEN_RE = re.compile( r""" (?P\d+\.\d*(?:[eE][+-]?\d+)? |\.\d+(?:[eE][+-]?\d+)? |\d+(?:[eE][+-]?\d+)?) | (?P[A-Za-z_]\w*) - | (?P\*\*|[-+*/^(),]) + | (?P\*\*|&&|\|\||<=|>=|==|!=|~=|[-+*/^(),<>]) | (?P\s+) """, re.VERBOSE, ) +_COMPARISONS = { + "<": lambda a, b: a < b, + ">": lambda a, b: a > b, + "<=": lambda a, b: a <= b, + ">=": lambda a, b: a >= b, + "==": lambda a, b: a == b, + "!=": lambda a, b: a != b, + "~=": lambda a, b: a != b, +} + def _tokenize(text: str) -> list[tuple[str, str]]: tokens: list[tuple[str, str]] = [] @@ -115,9 +157,14 @@ def _tokenize(text: str) -> list[tuple[str, str]]: class _Parser: """Recursive-descent parser for the arithmetic sublanguage. - Precedence, loosest to tightest: ``+ -``, then ``* /``, then unary ``-``, - then ``^`` (right associative). ``^`` binding tighter than unary minus is - what makes ``-2^2`` come out as ``-4``. + Precedence, loosest to tightest -- the order of ``arrayToExpression``'s + operator list in Expression.pm, which folds each level left to right: + + ``&& ||`` < ``< > <= >= == != ~=`` < ``+ -`` < ``* /`` < + unary ``- +`` < ``^`` + + Unary minus sitting *below* ``^`` is what makes ``-2^2`` come out as 4, and + the left fold is what makes ``2^3^2`` come out as 64. """ def __init__(self, tokens: list[tuple[str, str]], text: str, lookup): @@ -127,7 +174,7 @@ def __init__(self, tokens: list[tuple[str, str]], text: str, lookup): self._pos = 0 def parse(self) -> float: - value = self._parse_sum() + value = self._parse_logical() if self._pos != len(self._tokens): raise BnglExpressionError( f"Unexpected trailing input in {self._text!r} at token " @@ -145,6 +192,13 @@ def _accept(self, value: str) -> bool: return True return False + def _accept_any(self, values) -> str | None: + token = self._peek() + if token is not None and token[0] == "op" and token[1] in values: + self._pos += 1 + return token[1] + return None + def _expect(self, value: str) -> None: if not self._accept(value): found = self._peek() @@ -153,6 +207,26 @@ def _expect(self, value: str) -> None: + (repr(found[1]) if found else "end of expression") ) + def _parse_logical(self) -> float: + value = self._parse_comparison() + while True: + op = self._accept_any(("&&", "||")) + if op is None: + return value + rhs = self._parse_comparison() + # BNG2.pl normalises these to 1/0 rather than returning the operand + # the way bare Perl `||` would: `0||5` is 1.0, not 5. + value = float((value != 0 and rhs != 0) if op == "&&" + else (value != 0 or rhs != 0)) + + def _parse_comparison(self) -> float: + value = self._parse_sum() + while True: + op = self._accept_any(_COMPARISONS) + if op is None: + return value + value = float(_COMPARISONS[op](value, self._parse_sum())) + def _parse_sum(self) -> float: value = self._parse_product() while True: @@ -164,12 +238,12 @@ def _parse_sum(self) -> float: return value def _parse_product(self) -> float: - value = self._parse_unary() + value = self._parse_power() while True: if self._accept("*"): - value *= self._parse_unary() + value *= self._parse_power() elif self._accept("/"): - divisor = self._parse_unary() + divisor = self._parse_power() if divisor == 0: raise BnglExpressionError(f"Division by zero in {self._text!r}") # True division throughout: BNGL has no integer division, and @@ -179,19 +253,34 @@ def _parse_product(self) -> float: else: return value + def _parse_power(self) -> float: + # Left associative, and a signed operand belongs to the base rather than + # to the whole power: BNG2.pl gives -2^2 == 4 and 2^3^2 == 64. + value = self._parse_unary() + while self._accept("^"): + exponent = self._parse_unary() + try: + value = float(value ** exponent) + except (ArithmeticError, ValueError) as e: + # 0^-1, an overflow, or a negative base raised to a fractional + # power (which Python answers with a complex number). + raise BnglExpressionError( + f"Cannot raise {value!r} to the power {exponent!r} in " + f"{self._text!r}: {e}" + ) from e + except TypeError as e: # complex result from a negative fractional base + raise BnglExpressionError( + f"Cannot raise {value!r} to the power {exponent!r} in " + f"{self._text!r}" + ) from e + return value + def _parse_unary(self) -> float: if self._accept("-"): return -self._parse_unary() if self._accept("+"): return self._parse_unary() - return self._parse_power() - - def _parse_power(self) -> float: - base = self._parse_atom() - if self._accept("^"): - # Right associative, and the exponent may itself be signed. - return base ** self._parse_unary() - return base + return self._parse_atom() def _parse_atom(self) -> float: token = self._peek() @@ -205,20 +294,21 @@ def _parse_atom(self) -> float: if kind == "op" and value == "(": self._pos += 1 - inner = self._parse_sum() + inner = self._parse_logical() self._expect(")") return inner if kind == "name": self._pos += 1 if self._accept("("): - args = [self._parse_sum()] - while self._accept(","): - args.append(self._parse_sum()) + # `_pi()` and `_e()` take no arguments, so an empty list is legal. + args = [] + if self._peek() != ("op", ")"): + args.append(self._parse_logical()) + while self._accept(","): + args.append(self._parse_logical()) self._expect(")") return self._call(value, args) - if value in _CONSTANTS: - return _CONSTANTS[value] return self._lookup(value) raise BnglExpressionError(f"Unexpected token {value!r} in {self._text!r}") @@ -236,6 +326,10 @@ def _call(self, name: str, args: list[float]) -> float: raise BnglExpressionError( f"Wrong number of arguments to {name!r} in {self._text!r}" ) from e + except ArithmeticError as e: + raise BnglExpressionError( + f"{name}() could not be evaluated in {self._text!r}: {e}" + ) from e except ValueError as e: raise BnglExpressionError( f"{name}() is undefined for its argument in {self._text!r}: {e}" @@ -256,18 +350,8 @@ def _resolve_known(name: str, symbols: dict[str, float], text: str) -> float: ) from None -def evaluate_parameters(parameters: dict[str, str]) -> dict[str, float]: - """Resolve a BNGL parameters block to numbers. - - ``parameters`` maps a parameter name to its raw right-hand side, literal or - expression, as :func:`pybnf.petab._bngl.parse_model` collects it. Values are - resolved lazily in dependency order, so declaration order does not matter, - which matches BNG2.pl. - - Raises :class:`CircularParameterError` on a definition that depends on - itself, and :class:`BnglExpressionError` on anything unparseable or on a - reference to a name the block does not define. - """ +def _resolver(parameters: dict[str, str]): + """A memoizing ``lookup(name) -> float`` over a parameters block.""" resolved: dict[str, float] = {} resolving: list[str] = [] @@ -281,6 +365,11 @@ def lookup(name: str) -> float: ) if name not in parameters: raise BnglExpressionError(f"Unknown parameter {name!r}") + if name in RESERVED_NAMES: + raise BnglExpressionError( + f"{name!r} is a BNGL built-in function name and cannot be used as " + f"a parameter name" + ) resolving.append(name) try: value = _Parser( @@ -291,6 +380,50 @@ def lookup(name: str) -> float: resolved[name] = value return value + return lookup, resolved + + +def evaluate_parameters(parameters: dict[str, str]) -> dict[str, float]: + """Resolve a BNGL parameters block to numbers. + + ``parameters`` maps a parameter name to its raw right-hand side, literal or + expression, as :func:`pybnf.petab._bngl.parse_model` collects it. Values are + resolved lazily in dependency order, so a parameter may be defined before + the ones it depends on. (BNG2.pl itself is stricter here -- it drops a + forward-referencing parameter -- but accepting the order-independent form + costs nothing and loses no model BNG2.pl would have accepted.) + + Raises :class:`CircularParameterError` on a definition that depends on + itself, and :class:`BnglExpressionError` on anything unparseable or on a + reference to a name the block does not define. Use + :func:`evaluate_parameters_partial` when one bad definition should not cost + the caller the whole block. + """ + lookup, resolved = _resolver(parameters) for name in parameters: lookup(name) return resolved + + +def evaluate_parameters_partial( + parameters: dict[str, str], +) -> tuple[dict[str, float], dict[str, str]]: + """Resolve what can be resolved, and report the rest. + + Returns ``(values, errors)``: ``values`` maps each parameter that could be + computed to its number, and ``errors`` maps each one that could not to the + reason. Every parameter appears in exactly one of the two. + + A block is a single namespace, so one unusable definition should cost the + caller that parameter and whatever depends on it -- not the entire block. + """ + lookup, resolved = _resolver(parameters) + errors: dict[str, str] = {} + for name in parameters: + if name in resolved: + continue + try: + lookup(name) + except BnglExpressionError as e: + errors[name] = str(e) + return resolved, errors diff --git a/pybnf/petab/bngl_model.py b/pybnf/petab/bngl_model.py index 797e9a3bb..88e4e992c 100644 --- a/pybnf/petab/bngl_model.py +++ b/pybnf/petab/bngl_model.py @@ -24,12 +24,13 @@ import os import shutil import subprocess +import warnings from pathlib import Path from petab.v1.models.model import Model from ._bngl import parse_model -from ._bngl_expr import BnglExpressionError, evaluate_parameters +from ._bngl_expr import evaluate_parameters_partial #: BNGL model type, as used in a PEtab v2 yaml file as ``language``. MODEL_TYPE_BNGL = 'bngl' @@ -71,13 +72,18 @@ def get_parameter_ids(self): return list(self._entities.parameters) def _parameter_values(self): - """Every parameter resolved to a number, computed once and cached. + """``(values, errors)`` for the parameters block, computed once and cached. A parameters block is arithmetic over other parameters, so this needs no BNG2.pl and no network generation; see :mod:`pybnf.petab._bngl_expr`. + Resolution is *partial*: one unusable definition costs that parameter + and whatever depends on it, not the whole block. A real model can carry + a construct we do not evaluate (an NFsim ``TFUN``, say), and taking the + other 16 parameters down with it would be a worse failure than the one + #666 set out to fix. """ if self._resolved_parameters is None: - self._resolved_parameters = evaluate_parameters( + self._resolved_parameters = evaluate_parameters_partial( dict(self._entities.parameters) ) return self._resolved_parameters @@ -85,20 +91,31 @@ def _parameter_values(self): def get_parameter_value(self, id_): if id_ not in self._entities.parameters: raise ValueError(f"Parameter {id_} does not exist.") - try: - return self._parameter_values()[id_] - except BnglExpressionError as e: - raise ValueError( - f"Parameter '{id_}' has an expression value " - f"'{self._entities.parameters[id_]}' that could not be evaluated: {e}" - ) from e + values, errors = self._parameter_values() + if id_ in values: + return values[id_] + raise ValueError( + f"Parameter '{id_}' has an expression value " + f"'{self._entities.parameters[id_]}' that could not be evaluated: " + f"{errors[id_]}" + ) def get_free_parameter_ids_with_values(self): # Expression-valued parameters used to be skipped here, which lost them - # from the PEtab problem with nothing said. They are resolved now; a - # block that still cannot be evaluated raises rather than going quiet. - values = self._parameter_values() - return [(name, values[name]) for name in self._entities.parameters] + # from the PEtab problem with nothing said (#666). They are resolved + # now, and anything still unusable is named in a warning rather than + # disappearing -- but it no longer takes the rest of the block with it. + values, errors = self._parameter_values() + if errors: + detail = '; '.join(f'{name} ({errors[name]})' for name in sorted(errors)) + warnings.warn( + f"Model {self._model_id!r}: {len(errors)} of " + f"{len(self._entities.parameters)} parameters could not be evaluated " + f"and are omitted from the PEtab problem: {detail}", + stacklevel=2, + ) + return [(name, values[name]) + for name in self._entities.parameters if name in values] def get_valid_parameters_for_parameter_table(self): return list(self._entities.parameters) diff --git a/tests/test_petab_bngl_expr.py b/tests/test_petab_bngl_expr.py index d3e32dcf9..abf0a3b8e 100644 --- a/tests/test_petab_bngl_expr.py +++ b/tests/test_petab_bngl_expr.py @@ -1,6 +1,22 @@ -"""BNGL parameter-expression evaluation (issue #666).""" +"""BNGL parameter-expression evaluation (issue #666). + +The semantics here are not guesses: :data:`BNG_VERIFIED` is a table of +expressions with the value BNG2.pl 2.9.3 actually computes for them, obtained by +running each through ``writeNET({evaluate_expressions=>1})`` -- the only export +path that emits numbers rather than echoing the source text. + +The table is checked from both sides. :func:`test_bng_verified_table` pins our +evaluator against it with no BNG2.pl needed, so the contract holds in ordinary +CI; :func:`test_bng_verified_table_still_matches_bng2pl` re-derives the same +values from a real BNG2.pl when one is on PATH, so the table cannot quietly rot +if BioNetGen changes. Anything that disagrees is a bug in one of the two, which +is the point. +""" import math +import re +import shutil +import subprocess import pytest @@ -10,9 +26,173 @@ CircularParameterError, evaluate_expression, evaluate_parameters, + evaluate_parameters_partial, ) from pybnf.petab.bngl_model import BnglModel +#: ``(expression, value BNG2.pl computes)``. Self-contained, so each one can be +#: dropped straight into a parameters block. +BNG_VERIFIED = [ + # -- operators --------------------------------------------------------- + ('2^3', 8.0), + ('2**3', 8.0), # BNG2.pl accepts ** as a synonym for ^ + ('1/2', 0.5), # float division, never integer + ('8/4/2', 1.0), # / is left associative + ('1-2-3', -4.0), # - is left associative + ('1+2*3', 7.0), + ('(1+2)*3', 9.0), + ('2*-3', -6.0), + # Unary minus binds TIGHTER than ^, so this is (-2)^2, not -(2^2). + ('-2^2', 4.0), + ('-2^3', -8.0), + ('3*-2^2', 12.0), + ('-(2^2)', -4.0), # explicit parens do give -(2^2) + ('0-2^2', -4.0), # binary minus is looser, as usual + ('-exp(0)^2', 1.0), # the rule covers function calls too + # ^ is LEFT associative: (2^3)^2, not 2^(3^2). + ('2^3^2', 64.0), + ('2^2^3', 64.0), + ('4^0.5^2', 4.0), + ('2^(3^2)', 512.0), + ('2^-2', 0.25), + ('2^-2^2', 0.0625), + # -- comparison and logical, which yield 1.0/0.0 ----------------------- + ('1<2', 1.0), + ('2<1', 0.0), + ('1==1', 1.0), + ('1!=1', 0.0), + ('1~=2', 1.0), # ~= is BNG2.pl's alias for != + ('2&&3', 1.0), # normalised, unlike Perl's own && + ('0&&3', 0.0), + ('0||5', 1.0), # 1.0, not 5 + ('1+2>2', 1.0), # + binds tighter than > + ('1<2&&2<3', 1.0), # comparison binds tighter than && + ('if(1,5,7)', 5.0), + ('if(0,5,7)', 7.0), + ('if(2>1,5,7)', 5.0), + # -- functions --------------------------------------------------------- + ('_pi()', math.pi), # zero-argument functions, not bare names + ('_e()', math.e), + ('ln(_e())', 1.0), + ('exp(1)', math.e), + ('log10(1000)', 3.0), + ('log2(8)', 3.0), + ('sqrt(4)', 2.0), + ('abs(-3)', 3.0), + ('sin(1)', math.sin(1)), + ('cos(1)', math.cos(1)), + ('tan(1)', math.tan(1)), + ('asin(0.5)', math.asin(0.5)), + ('acos(0.5)', math.acos(0.5)), + ('atan(0.5)', math.atan(0.5)), + ('sinh(1)', math.sinh(1)), + ('cosh(1)', math.cosh(1)), + ('tanh(1)', math.tanh(1)), + ('asinh(1)', math.asinh(1)), + ('acosh(2)', math.acosh(2)), + ('atanh(0.5)', math.atanh(0.5)), + ('min(1,2)', 1.0), + ('min(3,1,2)', 1.0), # min/max/sum/avg are variadic + ('max(1,2)', 2.0), + ('sum(1,2,3,4)', 10.0), + ('avg(2,4)', 3.0), + # rint is floor(x + 0.5) -- round half UP, not Python's round-half-even. + ('rint(0.5)', 1.0), + ('rint(1.5)', 2.0), + ('rint(2.5)', 3.0), + ('rint(-0.5)', 0.0), + ('rint(-2.5)', -2.0), +] + +#: Expressions BNG2.pl refuses. Rejecting them keeps a typo an error instead of +#: a plausible wrong number. +BNG_REJECTS = [ + 'log(10)', # BNGL's natural log is ln; there is no bare log + 'floor(1.7)', # commented out in Expression.pm ("not supported by muParser") + 'ceil(1.2)', + '_pi', # the constants are functions: _pi(), not _pi + '_e', + 'foo(1)', + '1/0', + '2 @ 3', + 'if(1,5,1/0)', # BNG2.pl evaluates all three args, so this dies there too +] + + +@pytest.mark.parametrize('text, expected', BNG_VERIFIED, ids=[e for e, _ in BNG_VERIFIED]) +def test_bng_verified_table(text, expected): + """Our evaluator reproduces what BNG2.pl computes. No BNG2.pl needed.""" + assert evaluate_parameters({'z': text})['z'] == pytest.approx(expected) + + +@pytest.mark.parametrize('text', BNG_REJECTS) +def test_bng_rejected_expressions_are_rejected_here_too(text): + with pytest.raises(BnglExpressionError): + evaluate_parameters({'z': text}) + + +# -- the differential against a real BNG2.pl --------------------------------- + +_NET_PARAM = re.compile(r'^\s*\d+\s+(\w+)\s+(\S+)') + +_PROBE_MODEL = """\ +begin model +begin parameters +{block} +end parameters +begin molecule types + A() + B() +end molecule types +begin seed species + A() 1 +end seed species +begin reaction rules + A() -> B() 1 +end reaction rules +end model +generate_network({{overwrite=>1}}) +writeNET({{evaluate_expressions=>1,prefix=>"ev"}}) +""" + + +@pytest.mark.bionetgen +def test_bng_verified_table_still_matches_bng2pl(tmp_path): + """Re-derive :data:`BNG_VERIFIED` from BNG2.pl itself, in one run. + + Every expression goes into a single parameters block, so this costs one + BNG2.pl invocation rather than one per case. + """ + names = {f'p{i}': text for i, (text, _) in enumerate(BNG_VERIFIED)} + block = '\n'.join(f' {name} {text}' for name, text in names.items()) + (tmp_path / 'probe.bngl').write_text(_PROBE_MODEL.format(block=block)) + + proc = subprocess.run([shutil.which('BNG2.pl'), 'probe.bngl'], check=False, + cwd=tmp_path, capture_output=True, text=True, timeout=300) + net = tmp_path / 'ev.net' + assert net.exists(), f'BNG2.pl did not write a network:\n{proc.stdout}\n{proc.stderr}' + + computed, in_block = {}, False + for line in net.read_text().splitlines(): + if line.strip().startswith('begin parameters'): + in_block = True + continue + if line.strip().startswith('end parameters'): + break + if in_block: + m = _NET_PARAM.match(line.split('#')[0]) + if m: + computed[m.group(1)] = float(m.group(2)) + + mismatched = [] + for i, (text, expected) in enumerate(BNG_VERIFIED): + actual = computed.get(f'p{i}') + if actual is None or actual != pytest.approx(expected): + mismatched.append(f' {text!r}: table says {expected!r}, BNG2.pl says {actual!r}') + assert not mismatched, 'BNG2.pl disagrees with BNG_VERIFIED:\n' + '\n'.join(mismatched) + + +# -- the issue's own case ---------------------------------------------------- def test_expression_valued_parameter_is_resolved(): """The case from issue #666: kon is an expression over other parameters.""" @@ -37,35 +217,12 @@ def test_expression_valued_parameter_is_resolved(): assert 'kon' in ids -def test_caret_is_exponentiation_not_xor(): - """BNGL's ^ raises to a power; Python's is bitwise xor, and 2^3 there is 1.""" - assert evaluate_parameters({'a': '2^3'})['a'] == 8.0 - - -def test_exponentiation_binds_tighter_than_unary_minus(): - assert evaluate_parameters({'a': '-2^2'})['a'] == -4.0 - - -def test_exponentiation_is_right_associative(): - assert evaluate_parameters({'a': '2^3^2'})['a'] == 512.0 - - -def test_ln_is_the_natural_logarithm_and_bare_log_is_rejected(): - assert evaluate_parameters({'a': 'ln(_e)'})['a'] == pytest.approx(1.0) - assert evaluate_parameters({'a': 'log10(1000)'})['a'] == pytest.approx(3.0) - - # BNG2.pl has no bare log(); treating it as ln would turn a typo into a - # plausible wrong number rather than an error. - with pytest.raises(BnglExpressionError, match='log'): - evaluate_parameters({'a': 'log(10)'}) - - -def test_division_is_floating_point(): - assert evaluate_parameters({'a': '1/2'})['a'] == 0.5 - - def test_declaration_order_does_not_matter(): - """BNG2.pl resolves by dependency, not by position in the block.""" + """Lazy resolution by dependency, so a parameter may precede its inputs. + + BNG2.pl is stricter (it drops a forward-referencing parameter), but being + permissive here loses no model it would have accepted. + """ assert evaluate_parameters({'b': 'a*2', 'a': '3'}) == {'a': 3.0, 'b': 6.0} @@ -82,6 +239,10 @@ def test_declaration_order_does_not_matter(): ({'krZapTcr': '3', 'krZapCd3e': '10*krZapTcr'}, 'krZapCd3e', 30.0), ({'Kd_BRAF_RAFi2': '20', 'Gf_BRAF_RAFi2': 'ln(Kd_BRAF_RAFi2)'}, 'Gf_BRAF_RAFi2', math.log(20)), + # The real if() from blbr_heterogeneity_goldstein1980, which the first + # cut of this evaluator could not lex at all. + ({'LT': '3', 'RT': '1', 'excess_ratio': '1', + 'use_excess': 'if(LT/(RT+0.01)>=excess_ratio,1,0)'}, 'use_excess', 1.0), ], ) def test_real_world_expression_shapes(params, target, expected): @@ -105,15 +266,63 @@ def test_self_referential_definition_is_reported(): evaluate_parameters({'a': 'a+1'}) -@pytest.mark.parametrize( - 'rhs', - ['b', '2 +', 'foo(1)', '1/0', '2 @ 3'], -) +def test_builtin_name_is_rejected_as_a_parameter_name(): + """BNG2.pl: "Cannot use built-in function name '_e' as a parameter".""" + with pytest.raises(BnglExpressionError, match='built-in'): + evaluate_parameters({'_e': '5'}) + + +@pytest.mark.parametrize('rhs', ['b', '2 +', 'foo(1)', '1/0', '2 @ 3']) def test_unusable_expressions_raise_rather_than_go_quiet(rhs): with pytest.raises(BnglExpressionError): evaluate_parameters({'a': rhs}) +def test_evaluate_expression_against_known_symbols(): + assert evaluate_expression('x*2 + y', {'x': 1.5, 'y': 1.0}) == 4.0 + + +# -- partial resolution: one bad definition must not cost the whole block ---- + +def test_partial_resolution_keeps_the_usable_parameters(): + values, errors = evaluate_parameters_partial( + {'a': '2', 'b': 'a*3', 'bad': 'nosuch', 'c': '4'}) + assert values == {'a': 2.0, 'b': 6.0, 'c': 4.0} + assert set(errors) == {'bad'} + assert 'nosuch' in errors['bad'] + + +def test_partial_resolution_also_drops_dependents_of_a_bad_parameter(): + values, errors = evaluate_parameters_partial( + {'bad': 'nosuch', 'downstream': 'bad*2', 'fine': '1'}) + assert values == {'fine': 1.0} + assert set(errors) == {'bad', 'downstream'} + + +def test_every_parameter_is_either_resolved_or_reported(): + params = {'a': '1', 'b': 'a+1', 'c': 'oops', 'd': 'c*2'} + values, errors = evaluate_parameters_partial(params) + assert set(values) | set(errors) == set(params) + assert not (set(values) & set(errors)) + + +def test_one_unevaluable_parameter_does_not_take_down_the_model(): + """Regression: a whole-block abort loses more than the original bug did.""" + text = """ +begin model +begin parameters + good1 2 + good2 good1*3 + bad k__FREE +end parameters +end model +""" + m = BnglModel(parse_model(text), model_id='demo') + with pytest.warns(UserWarning, match='could not be evaluated'): + pairs = dict(m.get_free_parameter_ids_with_values()) + assert pairs == {'good1': 2.0, 'good2': 6.0} + + def test_unevaluable_parameter_surfaces_from_the_model(): """The adapter reports the failure instead of dropping the parameter.""" text = """ @@ -141,5 +350,17 @@ def test_missing_parameter_still_raises_value_error(): m.get_parameter_value('nope') -def test_evaluate_expression_against_known_symbols(): - assert evaluate_expression('x*2 + y', {'x': 1.5, 'y': 1.0}) == 4.0 +def test_fully_resolvable_model_warns_about_nothing(): + text = """ +begin model +begin parameters + a 2 + b a*3 +end parameters +end model +""" + m = BnglModel(parse_model(text), model_id='demo') + import warnings as _w + with _w.catch_warnings(): + _w.simplefilter('error') + assert dict(m.get_free_parameter_ids_with_values()) == {'a': 2.0, 'b': 6.0} diff --git a/tests/test_petab_export.py b/tests/test_petab_export.py index 03f8e6ca1..f885684b6 100644 --- a/tests/test_petab_export.py +++ b/tests/test_petab_export.py @@ -2142,10 +2142,11 @@ def test_is_state_variable_is_seed_species_only(self, model): assert not model.is_state_variable('x') # nor is an observable def test_expression_valued_parameter_is_evaluated(self): - # Superseded #666: an expression RHS used to raise NotImplementedError, + # Superseded by #666: an expression RHS used to raise NotImplementedError, # and get_free_parameter_ids_with_values dropped the parameter without # saying so. A parameters block is arithmetic over other parameters, so - # it is resolved here without BNG2.pl; see pybnf.petab._bngl_expr. + # it is resolved here without BNG2.pl; see pybnf.petab._bngl_expr, whose + # semantics are pinned against a real BNG2.pl in test_petab_bngl_expr.py. pytest.importorskip('petab') from pybnf.petab._bngl import parse_model from pybnf.petab.bngl_model import BnglModel