diff --git a/petab/v1/models/bngl_model.py b/petab/v1/models/bngl_model.py index 2f8cef39..3dcf0d0e 100644 --- a/petab/v1/models/bngl_model.py +++ b/petab/v1/models/bngl_model.py @@ -12,10 +12,16 @@ parse/semantic check, no network generation); otherwise the model is assumed valid. -Two things worth knowing if a model doesn't parse the way you expect: +Three things worth knowing if a model doesn't parse the way you expect: * Symbols usable in an observable formula are parameters, observables, and functions -- *not* compartments. +* A parameter whose value is an expression over other parameters + (``kon koff/(Kd*NA*V)``) is evaluated, since a parameters block is + arithmetic over other parameters and needs no reaction network. The + arithmetic follows BNGL rather than Python, so ``^`` is a power, it + groups from the left, and unary minus binds tighter than it does in + Python. See :func:`evaluate_bngl_parameters`. * The reader accepts line continuations (a trailing ``\\``), ``begin species`` as an alias for ``begin seed species``, line labels (both the numeric ``1 L0 1`` and named ``CD14: ...`` forms), and a leading ``$`` @@ -31,11 +37,13 @@ from __future__ import annotations +import math import os import re import shutil import subprocess -from collections.abc import Iterable +import warnings +from collections.abc import Callable, Iterable from dataclasses import dataclass from pathlib import Path @@ -247,6 +255,442 @@ def _compartment_name(line: str) -> str | None: return tokens[0] if tokens else None +# -- parameter expression evaluation ----------------------------------------- +# +# A ``parameters`` block may give a parameter an expression over other +# parameters (``kon koff/(Kd*NA*V)``) rather than a literal, which is +# ordinary BNGL style rather than an edge case: across 303 models drawn from +# the BioNetGen model collections, 1934 of 9323 parameter declarations +# (20.8%) are expression-valued. Resolving them needs no BNG2.pl and no +# network generation, because a parameters block is arithmetic over other +# parameters. +# +# 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})``, the only export path that emits +# numbers instead of echoing the source text. The function table and the +# precedence order are BioNetGen's ``Perl2/Expression.pm`` (``%functions``, +# ``%NARGS``, and the operator list in ``arrayToExpression``): +# +# * ``^`` raises to a power, where Python's is bitwise exclusive-or. +# * ``^`` is *left* associative, so ``2^3^2`` is 64 rather than 512. +# * Unary minus binds *tighter* than ``^``, so ``-2^2`` is ``(-2)^2`` == 4 +# rather than ``-(2^2)`` == -4. This holds for literals, parameters, +# parenthesised groups and function calls alike (``-exp(0)^2`` == 1). +# * The natural logarithm is ``ln``. A bare ``log`` is rejected, as BNG2.pl +# rejects it, so a typo stays an error instead of becoming a plausible +# wrong number. +# * ``rint`` is ``floor(x + 0.5)``, rounding a half upward, where Python's +# ``round`` sends a half to the nearest even number. +# * ``_pi`` and ``_e`` are zero-argument functions, written ``_pi()``. +# * 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. +# +# Expressions are tokenized and parsed rather than handed to ``eval``, which +# would import Python's precedence and operator meanings along with the +# obvious injection problem. + + +class BnglExpressionError(ValueError): + """A parameter expression could not be parsed or evaluated.""" + + +class CircularParameterError(BnglExpressionError): + """A parameter's definition depends on itself, directly or not.""" + + +def _bngl_if(condition: float, then_: float, else_: float) -> float: + """BNGL's ``if``, which selects on ``condition != 0``.""" + return then_ if condition != 0 else else_ + + +#: The built-in functions BNG2.pl accepts, mirroring ``%functions`` in +#: ``Expression.pm``. ``log`` is absent because BNGL has no bare ``log``. +#: ``floor`` and ``ceil`` are absent because ``Expression.pm`` keeps them +#: commented out as unsupported, so BNG2.pl rejects them. ``TFUN`` is absent +#: deliberately: it reads a data file while a simulation runs, so it is not a +#: parameters-block constant. +_FUNCTIONS: dict[str, Callable[..., float]] = { + "_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": _bngl_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_PARAMETER_NAMES = frozenset(_FUNCTIONS) + +# Longest-first, so ``**``, ``>=`` and ``&&`` 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: dict[str, Callable[[float, float], bool]] = { + "<": 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]]: + """``(kind, value)`` tokens for a BNGL expression.""" + 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} " + f"in {text!r}" + ) + pos = match.end() + kind = match.lastgroup + if kind == "space": + continue + value = match.group() + # BNGL writes exponentiation as ^, and BNG2.pl also accepts **. + tokens.append(("op", "^") if value == "**" else (kind, value)) + return tokens + + +class _Parser: + """Recursive-descent parser for the arithmetic sublanguage. + + Precedence, loosest to tightest, is the order of the operator list in + ``arrayToExpression``, 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: Callable[[str], float], + ): + self._tokens = tokens + self._text = text + self._lookup = lookup + self._pos = 0 + + def parse(self) -> float: + """The value of the whole expression.""" + 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: + if self._pos < len(self._tokens): + return self._tokens[self._pos] + return 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: Iterable[str]) -> 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() + seen = repr(found[1]) if found else "end of expression" + raise BnglExpressionError( + f"Expected {value!r} in {self._text!r}, found {seen}" + ) + + 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 an + # operand the way bare Perl would, so ``0||5`` is 1.0. + if op == "&&": + value = float(value != 0 and rhs != 0) + else: + value = float(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. + 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, 2^3^2 == 64. + value = self._parse_unary() + while self._accept("^"): + exponent = self._parse_unary() + try: + value = float(value**exponent) + except (ArithmeticError, TypeError, ValueError) as e: + # 0^-1, an overflow, or a negative base raised to a + # fractional power, which Python answers with a complex. + raise BnglExpressionError( + f"Cannot raise {value!r} to the power {exponent!r} " + f"in {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. + 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 " + f"{self._text!r}: {e}" + ) from e + + +def evaluate_bngl_expression(text: str, symbols: dict[str, float]) -> float: + """Evaluate one BNGL expression against already-resolved ``symbols``. + + :param text: The expression, for example ``koff/(Kd*NA*V)``. + :param symbols: Values for the names the expression refers to. + :returns: The value of the expression. + :raises BnglExpressionError: If it cannot be parsed or evaluated, or + refers to a name ``symbols`` does not define. + """ + + def lookup(name: str) -> float: + try: + return symbols[name] + except KeyError: + raise BnglExpressionError( + f"Unknown parameter {name!r} in {text!r}" + ) from None + + return _Parser(_tokenize(text), text, lookup).parse() + + +def _parameter_resolver( + parameters: dict[str, str], +) -> tuple[Callable[[str], float], dict[str, float]]: + """A memoizing ``lookup(name)`` over a parameters block, and its cache.""" + resolved: dict[str, float] = {} + resolving: list[str] = [] + + def lookup(name: str) -> float: + if name in resolved: + return resolved[name] + if name in resolving: + start = resolving.index(name) + cycle = " -> ".join([*resolving[start:], 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_PARAMETER_NAMES: + raise BnglExpressionError( + f"{name!r} is a BNGL built-in function name and cannot be " + f"used as 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_bngl_parameters( + parameters: dict[str, str], +) -> dict[str, float]: + """Resolve a BNGL parameters block to numbers. + + Values are resolved lazily in dependency order, so a parameter may be + defined before the ones it depends on. BNG2.pl is stricter here, since + it drops a forward-referencing parameter, but accepting the + order-independent form loses no model BNG2.pl would have accepted. + + :param parameters: Parameter name to raw right-hand side, literal or + expression, as :func:`parse_bngl` collects it. + :returns: Parameter name to value. + :raises CircularParameterError: On a definition that depends on itself. + :raises BnglExpressionError: On anything unparseable, or a reference to + a name the block does not define. Use + :func:`evaluate_bngl_parameters_partial` when one bad definition + should not cost the caller the whole block. + """ + lookup, resolved = _parameter_resolver(parameters) + for name in parameters: + lookup(name) + return resolved + + +def evaluate_bngl_parameters_partial( + parameters: dict[str, str], +) -> tuple[dict[str, float], dict[str, str]]: + """Resolve what can be resolved in a parameters block, and report the rest. + + A block is a single namespace, so one unusable definition should cost + the caller that parameter and whatever depends on it, rather than the + entire block. + + :param parameters: Parameter name to raw right-hand side. + :returns: ``(values, errors)``, where ``values`` maps each parameter + that could be computed to its value and ``errors`` maps each one + that could not to the reason. Every parameter appears in exactly + one of the two. + """ + lookup, resolved = _parameter_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 + + class BnglModel(Model): """PEtab wrapper for BNGL models.""" @@ -266,6 +710,9 @@ def __init__( self.model = model self._model_id = model_id + self._resolved_parameters: ( + tuple[dict[str, float], dict[str, str]] | None + ) = None if not is_valid_identifier(self._model_id): raise ValueError( @@ -305,33 +752,55 @@ def model_id(self, model_id): def get_parameter_ids(self) -> Iterable[str]: return list(self.model.parameters) + def _parameter_values(self) -> tuple[dict[str, float], dict[str, str]]: + """``(values, errors)`` for the parameters block, computed once. + + A parameters block is arithmetic over other parameters, so this + needs no BNG2.pl and no network generation. Resolution is partial: + one unusable definition costs that parameter and whatever depends + on it, rather than the whole block. + """ + if self._resolved_parameters is None: + self._resolved_parameters = evaluate_bngl_parameters_partial( + dict(self.model.parameters) + ) + return self._resolved_parameters + def get_parameter_value(self, id_: str) -> float: - try: - rhs = self.model.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 BNGL parameter expression requires BNG2.pl / " - "network generation, which is out of scope for the " - "introspection-only BnglModel." - ) from e + if id_ not in self.model.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.model.parameters[id_]}' that could not be evaluated: " + f"{errors[id_]}" + ) def get_free_parameter_ids_with_values( self, ) -> Iterable[tuple[str, float]]: - out = [] - for name, rhs in self.model.parameters.items(): - try: - out.append((name, float(rhs))) - except ValueError: - # An expression-valued parameter has no introspection-grade - # value; skip it rather than evaluate the expression. - continue - return out + # An expression-valued parameter used to be skipped here, which + # lost it from the PEtab problem with nothing said. They are + # resolved now, and anything still unusable is named in a warning + # rather than disappearing, without taking 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.model.parameters)} parameters could not be " + f"evaluated and are omitted: {detail}", + stacklevel=2, + ) + return [ + (name, values[name]) + for name in self.model.parameters + if name in values + ] def get_valid_parameters_for_parameter_table(self) -> Iterable[str]: # All parameters are allowed in the parameter table. diff --git a/tests/v1/test_model_bngl.py b/tests/v1/test_model_bngl.py index b4641faf..9933a5db 100644 --- a/tests/v1/test_model_bngl.py +++ b/tests/v1/test_model_bngl.py @@ -49,19 +49,22 @@ def test_get_parameter_value_unknown_raises(model): model.get_parameter_value("nope") -def test_expression_valued_parameter_is_not_evaluated(): - # A numeric RHS coerces to float; an expression RHS is confined to - # NotImplementedError rather than evaluated (that needs BNG2.pl). +def test_expression_valued_parameter_is_evaluated(): + # 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 without BNG2.pl. The BNGL arithmetic itself is + # pinned against a real BNG2.pl in test_model_bngl_expressions.py. entities = parse_bngl( "begin parameters\n base 2\n k_on 2*base\nend parameters\n" ) model = BnglModel(entities, model_id="m") assert model.get_parameter_value("base") == 2.0 - with pytest.raises(NotImplementedError): - model.get_parameter_value("k_on") - # The expression-valued parameter is still an enumerated entity, but it - # contributes no introspection-grade value. - assert dict(model.get_free_parameter_ids_with_values()) == {"base": 2.0} + assert model.get_parameter_value("k_on") == 4.0 + assert dict(model.get_free_parameter_ids_with_values()) == { + "base": 2.0, + "k_on": 4.0, + } # -- grammar hardening: block aliases + seed-species "$" clamp --------------- diff --git a/tests/v1/test_model_bngl_expressions.py b/tests/v1/test_model_bngl_expressions.py new file mode 100644 index 00000000..e9addc90 --- /dev/null +++ b/tests/v1/test_model_bngl_expressions.py @@ -0,0 +1,365 @@ +"""BNGL parameter-expression evaluation. + +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 +the evaluator against it with no BNG2.pl needed, so the contract holds in +ordinary continuous integration. :func:`test_table_still_matches_bng2pl` +re-derives the same values from a real BNG2.pl where one is available, so +the table cannot quietly rot if BioNetGen changes. +""" + +import math +import re +import subprocess + +import pytest + +from petab.v1.models.bngl_model import ( + BnglExpressionError, + BnglModel, + CircularParameterError, + _locate_bng2, + evaluate_bngl_expression, + evaluate_bngl_parameters, + evaluate_bngl_parameters_partial, + parse_bngl, +) + +#: ``(expression, the 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), rounding a half up, where Python's round + # sends a half to the nearest even number. + ("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 rather +#: than a plausible wrong number. +BNG_REJECTS = [ + "log(10)", # BNGL's natural log is ln, and there is no bare log + "floor(1.7)", # commented out in Expression.pm as unsupported + "ceil(1.2)", + "_pi", # the constants are functions, written _pi() + "_e", + "foo(1)", + "1/0", + "2 @ 3", + "if(1,5,1/0)", # BNG2.pl evaluates all three arguments +] + + +@pytest.mark.parametrize( + "text, expected", BNG_VERIFIED, ids=[e for e, _ in BNG_VERIFIED] +) +def test_bng_verified_table(text, expected): + """The evaluator reproduces what BNG2.pl computes.""" + assert evaluate_bngl_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_bngl_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"}}) +""" + + +def test_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. + """ + bng2 = _locate_bng2() + if bng2 is None: + pytest.skip("BNG2.pl not available") + + names = {f"p{i}": text for i, (text, _) in enumerate(BNG_VERIFIED)} + block = "\n".join(f" {n} {t}" for n, t in names.items()) + (tmp_path / "probe.bngl").write_text(_PROBE_MODEL.format(block=block)) + + proc = subprocess.run( # noqa: S603 + [bng2, "probe.bngl"], + cwd=tmp_path, + capture_output=True, + text=True, + timeout=300, + check=False, + ) + net = tmp_path / "ev.net" + assert net.exists(), ( + f"BNG2.pl wrote no 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) + ) + + +# -- resolution order and cycles --------------------------------------------- + + +def test_expression_over_other_parameters_resolves(): + text = ( + "begin parameters\n" + " NA 6.022e23\n" + " V 1e-12\n" + " Kd 5.0\n" + " koff 0.1\n" + " kon koff/(Kd*NA*V)\n" + "end parameters\n" + ) + model = BnglModel(parse_bngl(text), model_id="demo") + assert model.get_parameter_value("kon") == pytest.approx( + 0.1 / (5.0 * 6.022e23 * 1e-12) + ) + ids = [name for name, _ in model.get_free_parameter_ids_with_values()] + assert ids == list(model.get_parameter_ids()) + + +def test_declaration_order_does_not_matter(): + """A parameter may be defined before the ones it depends on.""" + assert evaluate_bngl_parameters({"b": "a*2", "a": "3"}) == { + "a": 3.0, + "b": 6.0, + } + + +def test_chained_expression_dependencies_resolve(): + values = evaluate_bngl_parameters({"a": "2", "b": "a*3", "c": "b+a"}) + assert values == {"a": 2.0, "b": 6.0, "c": 8.0} + + +@pytest.mark.parametrize( + "params, target, expected", + [ + # Shapes taken from real BioNetGen models. + ( + { + "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": "20", "Gf_BRAF": "ln(Kd_BRAF)"}, + "Gf_BRAF", + math.log(20), + ), + ( + { + "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_bngl_parameters(params)[target] == pytest.approx(expected) + + +def test_circular_definition_names_the_cycle(): + with pytest.raises(CircularParameterError) as excinfo: + evaluate_bngl_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_bngl_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_bngl_parameters({"_e": "5"}) + + +def test_evaluate_expression_against_known_symbols(): + assert evaluate_bngl_expression("x*2 + y", {"x": 1.5, "y": 1.0}) == 4.0 + + +# -- partial resolution ------------------------------------------------------ + + +def test_partial_resolution_keeps_the_usable_parameters(): + values, errors = evaluate_bngl_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"} + + +def test_partial_resolution_also_drops_dependents_of_a_bad_parameter(): + values, errors = evaluate_bngl_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_bngl_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(): + """A whole-block failure would lose more than the original bug did.""" + text = ( + "begin parameters\n" + " good1 2\n" + " good2 good1*3\n" + " bad not_a_parameter\n" + "end parameters\n" + ) + model = BnglModel(parse_bngl(text), model_id="demo") + with pytest.warns(UserWarning, match="could not be evaluated"): + pairs = dict(model.get_free_parameter_ids_with_values()) + assert pairs == {"good1": 2.0, "good2": 6.0} + + +def test_unevaluable_parameter_surfaces_from_the_model(): + text = "begin parameters\n a b\nend parameters\n" + model = BnglModel(parse_bngl(text), model_id="demo") + with pytest.raises(ValueError, match="could not be evaluated"): + model.get_parameter_value("a") + + +def test_missing_parameter_still_raises_value_error(): + text = "begin parameters\n a 1\nend parameters\n" + model = BnglModel(parse_bngl(text), model_id="demo") + with pytest.raises(ValueError, match="does not exist"): + model.get_parameter_value("nope")