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 3786f6d2..199a785e 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 new file mode 100644 index 00000000..57ddfcb2 --- /dev/null +++ b/pybnf/petab/_bngl_expr.py @@ -0,0 +1,429 @@ +"""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. 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. 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 +``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", + "evaluate_parameters_partial", +] + + +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.""" + + +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_ + + +# 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, + "asin": math.asin, + "acos": math.acos, + "atan": math.atan, + "sinh": math.sinh, + "cosh": math.cosh, + "tanh": math.tanh, + "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\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]] = [] + 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 -- 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): + self._tokens = tokens + self._text = text + self._lookup = lookup + self._pos = 0 + + def parse(self) -> float: + value = self._parse_logical() + 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 _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() + raise BnglExpressionError( + f"Expected {value!r} in {self._text!r}, found " + + (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: + 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_power() + while True: + if self._accept("*"): + value *= self._parse_power() + elif self._accept("/"): + 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 + # 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_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_atom() + + 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_logical() + self._expect(")") + return inner + + if kind == "name": + self._pos += 1 + if self._accept("("): + # `_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) + 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 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}" + ) 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 _resolver(parameters: dict[str, str]): + """A memoizing ``lookup(name) -> float`` over a parameters block.""" + 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}") + 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( + _tokenize(parameters[name]), parameters[name], lookup + ).parse() + finally: + resolving.pop() + 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 843cb138..88e4e992 100644 --- a/pybnf/petab/bngl_model.py +++ b/pybnf/petab/bngl_model.py @@ -24,11 +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 evaluate_parameters_partial #: BNGL model type, as used in a PEtab v2 yaml file as ``language``. MODEL_TYPE_BNGL = 'bngl' @@ -44,6 +46,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 +71,51 @@ def model_id(self): def get_parameter_ids(self): return list(self._entities.parameters) + def _parameter_values(self): + """``(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_partial( + dict(self._entities.parameters) + ) + return self._resolved_parameters + def get_parameter_value(self, id_): - 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)." - ) from e + if id_ not in self._entities.parameters: + raise ValueError(f"Parameter {id_} does not exist.") + 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): - 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 (#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 new file mode 100644 index 00000000..abf0a3b8 --- /dev/null +++ b/tests/test_petab_bngl_expr.py @@ -0,0 +1,366 @@ +"""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 + +from pybnf.petab._bngl import parse_model +from pybnf.petab._bngl_expr import ( + BnglExpressionError, + 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.""" + 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_declaration_order_does_not_matter(): + """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} + + +@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)), + # 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): + 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'}) + + +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 = """ +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_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 a96407de..f885684b 100644 --- a/tests/test_petab_export.py +++ b/tests/test_petab_export.py @@ -2141,7 +2141,12 @@ 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 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, 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 @@ -2149,8 +2154,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