diff --git a/docs/adr/0147-an-sbml-entity-an-initialassignment-derives-has-no-value-the-model-file-settles-so-the-scanner-drops-the-declared-attribute-and-the-measurement-layer-inlines-the-derivation.md b/docs/adr/0147-an-sbml-entity-an-initialassignment-derives-has-no-value-the-model-file-settles-so-the-scanner-drops-the-declared-attribute-and-the-measurement-layer-inlines-the-derivation.md new file mode 100644 index 00000000..b6b9c917 --- /dev/null +++ b/docs/adr/0147-an-sbml-entity-an-initialassignment-derives-has-no-value-the-model-file-settles-so-the-scanner-drops-the-declared-attribute-and-the-measurement-layer-inlines-the-derivation.md @@ -0,0 +1,103 @@ +# An SBML entity an initialAssignment derives has no value the model file settles, so the scanner drops the declared attribute and the measurement layer inlines the derivation (issue #795) + +## Status + +Accepted and implemented (2026-09-22), for every edition. It changes the numbers PyBNF reads +out of an SBML or antimony model only when that model uses an `initialAssignment`, which no +model committed to this repository does on a parameter. + +## The problem + +SBML lets a model give an entity its starting value in two places. A parameter carries a +`value` attribute, a compartment a `size`, and a species an `initialAmount` or an +`initialConcentration`. A `` entry may then supersede any of them, +and when it does, the attribute is a placeholder the model never starts from. Tool-exported +SBML writes both routinely, and antimony writes the assignment with no attribute at all, so +`k_derived = k_base + 1` becomes `` plus an +assignment. + +The import-time scanner, `pybnf/petab/_sbml.py`, read only the attributes. Two consumers used +the result, and both were wrong in the same way: + +- `Configuration._model_expression_namespace` puts the scanned values into the measurement + layer's constant snapshot, so an observable or noise formula naming such a parameter was + scored with the placeholder. On Bertozzi_PNAS2020, where `beta_N` carries `value="0.0"` and + an assignment computing it from three parameters, that is a fit scored with a rate constant + of zero, silently. +- The exporter's `_numeric_nominal` feeds `conditions.mutation_target_value`, which turns a + relative condition into an absolute number, so a `condition: kon * 2` on such a parameter + wrote a wrong number into `conditions.tsv`. + +Where no attribute is written the failure was different and not much better. The symbol stayed +in the namespace with no value behind it, so a formula naming it reached the branch in +`pybnf/measurement/base.py` that reports the symbol as "neither a simulation-output column nor +a fit/model parameter", a message that is wrong and that sits under a comment saying validation +should have made it unreachable. + +Measured on the public PEtab benchmark collection, 16 of 25 models use initial assignments, 31 +of them on a parameter or compartment, and 28 of those still handed out a stale attribute. + +## The decision + +An entity whose initial value an assignment derives has **no value the model file settles**, so +the scanner reports none, and the measurement layer inlines the derivation instead of binding a +number. + +Concretely, for a parameter or a compartment: + +- An assignment that is arithmetic over numbers alone **is** the value. It replaces the + attribute. `_evaluate_mathml_number` walks the MathML and refuses a ``, which makes + "evaluates to a number" and "reads no other entity" one test rather than two. +- An assignment computed from other entities means the entity leaves `parameter_values` and + leaves `namespace_symbols`, and its serialized right-hand side is recorded in + `derived_initial_values`. A formula naming it is rewritten to the entities it is computed + from, which is what `#465` already does for an `assignmentRule` target, through the same + inliner and the same error shape. + +A **species** is deliberately different. An initial assignment pins t=0 only, so the species is +still a dynamical state and still a simulation-output column. It keeps its place in the +namespace and is never inlined; it only loses a stale declared initial. Treating it like a rule +target would break a committed fixture on the first run, since Boehm's observables name +`STAT5A` and `STAT5B`, both set by initial assignment. + +Inlining is gated on soundness. The substituted expression is read at a measurement time, while +the assignment fixed a value from its inputs at t=0, so the two agree only when every input +holds still for the whole simulation. The scanner therefore refuses to inline an assignment +that reads a species, a `constant="false"` parameter, an assignment-rule or rate-rule target, +an event-assignment target, or any symbol an algebraic rule touches. A refused assignment is +recorded with no expression and a clause naming the offender, which the loader composes into +one error. + +## Why not the alternatives + +Dropping the value and leaving the symbol in the namespace with a better error message was +rejected because it turns 28 measured silent wrong numbers into 28 refusals when a correct +answer is available, and because the better message would live in a branch the code says is +unreachable. + +Dropping the value and excluding the symbol with no inlining was rejected for the same first +reason. Both alternatives also leave a user whose model simulates correctly unable to write an +observable over it without rewriting the expression by hand. + +Evaluating a derived parameter from the file's own values, which is the obvious third option, +is the defect this repository has just argued against upstream. A PEtab parameter table may +override or estimate what the value is computed from, so a number settled from file defaults is +handed to a simulator as a constant and goes stale. The same rule now holds on both sides: +petab's `SbmlModel` reports an initial assignment only when it reduces to a number, and +`PEtab-dev/libpetab-python#517` applies it to BioNetGen parameter expressions. + +## Consequences + +The import-time view now agrees with the runtime. ADR-0094 already recomputes a parameter an +initial assignment fixes whenever a dependency moves, and because `MeasurementModel.materialize` +resolves a symbol from the parameter set before the constant snapshot, an inlined dependency +that the fit estimates tracks the fit. An observable over Bertozzi's `beta_N` is scored with +the value the model starts from, and it moves when `R0_` is estimated. + +Two adjacent defects are fixed in the same pass, because they are the same sentence. An +`assignmentRule` target that carries a vestigial `value` attribute no longer reports that +number either, and `sigma = prediction_formula` now inlines the same map the measurement layer +does, so the two layers no longer disagree about which symbols a model offers. + +The refusal in `conditions.mutation_target_value` no longer blames BNGL alone. It names the +parameter, states both languages' version of the same cause, and says what to write instead. diff --git a/pybnf/config.py b/pybnf/config.py index e46db233..c16e2847 100644 --- a/pybnf/config.py +++ b/pybnf/config.py @@ -2505,9 +2505,9 @@ def _load_measurement_models(self): ed, 2, "the 'observable: , formula: ' measurement-model syntax") from .measurement import MeasurementLayer, MeasurementModel, PerMeasurementModel - from .petab.formula import compile_petab_formula, inline_assignment_rules + from .petab.formula import compile_petab_formula, inline_derived_symbols - namespace, constants, assignment_rules = self._model_expression_namespace() + namespace, constants, derived = self._model_expression_namespace() free_names = {v.name for v in self.variables} # Free parameters resolve from the PSet at eval time, not from the constant snapshot. constants = {n: val for n, val in constants.items() if n not in free_names} @@ -2521,18 +2521,20 @@ def _load_measurement_models(self): models = [] # constant measurement models -> pre-materialized layer per_measurement = {} # row-varying ones -> bound per data point in the objective for obs_id, formula in specs: - # An SBML assignment-rule variable is declared in the model (so it would pass the - # namespace check if it were not excluded) but is algebraically computed -- never a - # simulation-output column and value-less -- so it cannot be resolved as a symbol at - # fit time. Inline it instead: substitute the rule's RHS down to the species the rule - # is defined over (recursively), so `observable: Epo_cells, formula: Epo_cells` just - # works -- the D2D convenience observable IS an assignment rule (#465, the option-2 - # successor to #464's reconstruct-from-species rejection). A formula naming no rule - # variable is returned verbatim; an unresolvable rule (untranslatable MathML / a - # cyclic dependency) raises a pointed error here at load, not late in materialize. - if assignment_rules: - formula = inline_assignment_rules( - formula, assignment_rules, observable_id=obs_id) + # A derived SBML entity is declared in the model (so it would pass the namespace check + # if it were not excluded) but has no value of its own, so it cannot be resolved as a + # symbol at fit time. Inline it instead: substitute its defining RHS down to the + # species/parameters it is computed from (recursively), so + # `observable: Epo_cells, formula: Epo_cells` just works -- the D2D convenience + # observable IS an assignment rule (#465, the option-2 successor to #464's + # reconstruct-from-species rejection), and an observable over a parameter an + # initialAssignment derives scores the value the model actually starts from instead of + # a stale declared attribute (#795). A formula naming no derived entity is returned + # verbatim; one that cannot be resolved (untranslatable MathML, a cyclic dependency, or + # an initial assignment over something that moves) raises a pointed error here at + # load, not late in materialize. + if derived: + formula = inline_derived_symbols(formula, derived, observable_id=obs_id) # A surviving per-measurement placeholder marks a row-varying scale/offset (ADR-0045): # its token differs per data row, so it is bound per data point, NOT pre-materialized. # Admit the placeholder symbol(s) to the allowed set so the rest validates fail-fast. @@ -2608,11 +2610,18 @@ def _load_prediction_noise(self): ed = edition.resolve_edition(self.config.get('edition')) edition.require_edition( ed, 2, "a prediction-dependent 'sigma = prediction_formula ' noise source") - from .petab.formula import compile_petab_formula - namespace, _constants, _rules = self._model_expression_namespace() + from .petab.formula import compile_petab_formula, inline_derived_symbols + namespace, _constants, derived = self._model_expression_namespace() free_names = {v.name for v in self.variables} allowed = namespace | free_names for label, src in sources: + # A prediction noise formula reads the same model namespace as a measurement formula, + # so a derived SBML entity has to be inlined here too. Without this the sibling layers + # disagree: the same symbol works in `observable:` and is rejected as unknown in + # `sigma = prediction_formula` (#465 for a rule target, #795 for an initialAssignment). + if derived: + src.formula = inline_derived_symbols( + src.formula, derived, observable_id=label) # Validate every symbol is a model entity or a declared free parameter, and adopt the # compiler's canonical ordering (the callable is rebuilt lazily worker-side, dropped # across pickling -- this is the validation pass, ADR-0036 §5). @@ -2664,18 +2673,19 @@ def _model_expression_namespace(self): formula over a ``.ant`` model's species must not be rejected as "not a known model entity" just because the namespace was built by parsing the file as BNGL (#463). - Returns ``(namespace_symbols, constants, assignment_rules)``, where ``assignment_rules`` - maps an SBML assignment-rule variable (declared as a parameter, but algebraically - computed -- never a simulation-output column and value-less, so it cannot be resolved as - a symbol at ``materialize``) to its rule's RHS as a PEtab-math infix string (``None`` if - the rule's MathML was not translatable). Such a variable is **excluded** from - ``namespace_symbols``; ``assignment_rules`` lets the loader **inline** a formula that - references one down to the species the rule is defined over (#465).""" + Returns ``(namespace_symbols, constants, derived)``, where ``derived`` maps every SBML + entity the model file defines in terms of others to a ``DerivedSymbol``: an + assignment-rule variable, which is algebraically computed at every step (#465), and a + parameter or compartment an ``initialAssignment`` derives, whose declared attribute is a + placeholder the model never starts from (#795). Neither can be resolved as a symbol at + ``materialize``, so both are **excluded** from ``namespace_symbols``, and ``derived`` lets + the loader **inline** a formula that references one down to the entities it is computed + from.""" from .petab._bngl import parse_model as parse_bngl from .petab._sbml import parse_model as parse_sbml namespace = set() constants = {} - assignment_rules = {} + derived = {} for mf in self.config['models']: text = Path(self._absolute(mf)).read_text(encoding='utf-8', errors='replace') if mf.endswith('.xml') or mf.endswith('.ant'): @@ -2689,7 +2699,7 @@ def _model_expression_namespace(self): ent = parse_sbml(text) namespace |= ent.namespace_symbols constants.update(ent.constants) - assignment_rules.update(ent.assignment_rules) + derived.update(ent.derived_symbols) else: # .bngl -- the BNGL ParamList (parameters u observables u functions) ent = parse_bngl(text) namespace |= (set(ent.parameters) | set(ent.observable_names) @@ -2699,7 +2709,7 @@ def _model_expression_namespace(self): constants[name] = float(rhs) except (TypeError, ValueError): pass # an expression-valued parameter is not a numeric constant - return namespace, constants, assignment_rules + return namespace, constants, derived def _load_simulators(self): diff --git a/pybnf/petab/_sbml.py b/pybnf/petab/_sbml.py index f833d59d..a39c2dc4 100644 --- a/pybnf/petab/_sbml.py +++ b/pybnf/petab/_sbml.py @@ -13,10 +13,61 @@ *direct children of the model element* and never descends into reactions. SBML uses XML namespaces (``{http://www.sbml.org/...}species``), so every tag is matched by its *local* name. The model file itself is carried **verbatim** (ADR-0036); this reads it, never edits it. + +A ```` entry supersedes the value an entity declares as an attribute, +so it is read here too: the assignment settles the value when it is arithmetic over numbers, and +otherwise the model file settles no value at all and the declared attribute is a placeholder the +model never starts from (#795). """ +import math import xml.etree.ElementTree as ET from dataclasses import dataclass +from typing import NamedTuple + + +class DerivedSymbol(NamedTuple): + """One entity the model file defines in terms of others, and how to resolve it. + + ``kind`` is ``'assignment_rule'`` (#465) or ``'initial_assignment'`` (#795). ``expression`` + is the defining RHS as a PEtab-math infix string, or ``None`` when the definition cannot be + trusted as a substitution -- then ``refusal`` is the clause naming why, which the measurement + layer composes into its error. + """ + + kind: str + expression: str + refusal: str + + +#: An ``assignmentRule`` whose MathML this stdlib serializer does not translate (#465). +_R_RULE_UNTRANSLATABLE = ("its defining assignment rule uses a construct this stdlib reader does " + "not translate, such as a piecewise or a relational operator") +#: The same for an ``initialAssignment`` (#795). ``time`` is listed because a `csymbol` for it is +#: the construct an author is most likely to reach for in an initial assignment. +_R_INITIAL_UNTRANSLATABLE = ("its defining initial assignment uses a construct this stdlib reader " + "does not translate, such as a piecewise, a relational operator, or " + "the symbol for time") + + +def _r_time_varying(offenders): + """The refusal for an ``initialAssignment`` computed from something that moves (#795). + + An initial assignment fixes a value from its inputs' *initial* values. Substituting it into a + measurement formula would read those inputs at the measurement time instead, so the + substitution is only sound when every input is constant for the whole simulation. + """ + quoted = [f"'{name}'" for name in offenders] + if len(quoted) == 1: + named = quoted[0] + elif len(quoted) == 2: + named = f'{quoted[0]} and {quoted[1]}' + else: + named = ', '.join(quoted[:-1]) + f' and {quoted[-1]}' + changes = 'whose value changes' if len(quoted) == 1 else 'whose values change' + return (f'its defining initial assignment is computed from {named}, {changes} during the ' + f'simulation, and an initial assignment fixes a value from the initial values of its ' + f'inputs rather than from their values at a measurement time') @dataclass(frozen=True) @@ -25,8 +76,8 @@ class SbmlEntities: ``species_names`` are the floating/boundary species (the trajectory's output columns at run time). ``parameter_names`` are the **global** parameters; ``compartment_names`` the - compartments. ``species_initial`` maps a species id to its initial amount/concentration, - and ``parameter_values`` maps a global parameter or compartment id to its numeric value -- + compartments. ``species_initial`` maps a species id to the initial value the model file + settles, and ``parameter_values`` maps a global parameter or compartment id to its value -- the fixed-constant snapshot a :class:`~pybnf.measurement.MeasurementModel` resolves a non-column, non-PSet symbol against (ADR-0036 §4). @@ -41,6 +92,21 @@ class SbmlEntities: rule into a formula that references it -- ``observable: Epo_cells, formula: Epo_cells`` just works, resolving down to the species the rule is computed from (#465, the option-2 successor to #464's reconstruct-from-species rejection). + + ``derived_initial_values`` is the same idea for an ``initialAssignment``, which SBML lets + supersede a parameter's ``value``, a compartment's ``size``, and a species' initial amount or + concentration (#795). When the assignment is arithmetic over numbers alone the scan evaluates + it and that number *is* the entity's value. When it is computed from other entities the model + file settles no value at all -- the declared attribute is a placeholder the model never starts + from -- so the id is dropped from ``parameter_values`` and recorded here with its RHS, exactly + as a rule target is. The id keeps its place in ``parameter_names``/``compartment_names``, + because the scan stays faithful to the file. + + A **species** whose initial value an assignment derives is the one case that works + differently, and the difference is the point: an ``initialAssignment`` pins t=0 only, so the + species is still a genuine dynamical state and still a simulation-output column. It stays in + :attr:`namespace_symbols`, it is never recorded in ``derived_initial_values``, and it simply + loses its stale entry in ``species_initial``. """ text: str @@ -50,17 +116,39 @@ class SbmlEntities: species_initial: dict # 'S1' -> 10.0 parameter_values: dict # 'k1' -> 0.5, 'cell' -> 1.0 (params u compartment sizes) assignment_rules: dict # 'Epo_cells' -> 'Epo_EpoRi + dEpoi' (RHS infix; None if untranslatable) (#465) + derived_initial_values: dict = None # 'beta_N' -> '(R0_ * gamma_) / N_' (None if not inlinable) (#795) + derived_refusals: dict = None # 'beta_N' -> the clause naming why its RHS is None (#795) @property def namespace_symbols(self): """The symbols an ``observableFormula`` may reference: species u parameters u compartments (the SBML analogue of the BNGL ``ParamList``, ADR-0026/0036), **minus** any assignment-rule variable -- which is declared as a parameter but is not resolvable - at ``materialize`` *as a symbol* (not an output column, no fixed value). A formula naming - one is resolved by inlining the rule's RHS down to species (#465), not by binding the - symbol, so the target stays out of the namespace.""" + at ``materialize`` *as a symbol* (not an output column, no fixed value) -- and minus any + parameter or compartment an ``initialAssignment`` derives, for the same reason (#795). A + formula naming either is resolved by inlining the definition (#465/#795), not by binding + the symbol, so the target stays out of the namespace. A *species* with an initial + assignment is **not** subtracted: it is still an output column.""" return ((self.species_names | self.parameter_names | self.compartment_names) - - set(self.assignment_rules)) + - set(self.assignment_rules) - set(self.derived_initial_values or {})) + + @property + def derived_symbols(self): + """Every entity the model file defines in terms of others, id -> :class:`DerivedSymbol`. + + The single map the measurement and import layers inline through. An ``assignmentRule`` + target (#465) and an ``initialAssignment``-derived parameter or compartment (#795) are the + same kind of entity to a formula, so they get one map, one namespace subtraction and one + error shape. A rule wins the merge: SBML forbids a symbol having both, and a file that + writes both anyway is governed by the rule at t=0 as well.""" + refusals = self.derived_refusals or {} + out = {name: DerivedSymbol('initial_assignment', rhs, + None if rhs is not None else refusals.get(name)) + for name, rhs in (self.derived_initial_values or {}).items()} + out.update({name: DerivedSymbol('assignment_rule', rhs, + None if rhs is not None else _R_RULE_UNTRANSLATABLE) + for name, rhs in self.assignment_rules.items()}) + return out @property def constants(self): @@ -81,6 +169,12 @@ def parse_model(text): species, species_initial = {}, {} parameters, compartments = {}, {} assignment_rules = {} + # The entities whose value moves during a simulation, collected so an initialAssignment + # computed from one of them is never inlined into a measurement formula (#795). A + # `constant="false"` parameter is included for the same reason: SBML lets it be changed by a + # rule or an event, so its initial value is not its value at a measurement time. + rate_rule_targets, algebraic_symbols, event_targets, non_constant = set(), set(), set(), set() + initial_elems = {} for container in list(model): ctag = _local(container.tag) if ctag == 'listOfSpecies': @@ -97,11 +191,34 @@ def parse_model(text): pid = e.get('id') if pid: parameters[pid] = _float_or_none(e.get('value')) + if not _is_constant(e): + non_constant.add(pid) elif ctag == 'listOfCompartments': for e in _children(container, 'compartment'): cid = e.get('id') if cid: compartments[cid] = _float_or_none(e.get('size')) + if not _is_constant(e): + non_constant.add(cid) + elif ctag == 'listOfInitialAssignments': + # An supersedes X's declared value/size/initial + # amount, so the attribute read above is a placeholder the model never starts from. + # The elements are kept and settled after the loop, so the result does not depend on + # the order the listOf* containers appear in (#795). + for e in _children(container, 'initialAssignment'): + sym = e.get('symbol') + if sym: + initial_elems[sym] = e + elif ctag == 'listOfEvents': + # The one place this scan reads below a direct child of , and only to learn + # which entities an event assigns to. An event's trigger/delay/priority math is not + # read. + for ev in _children(container, 'event'): + for lst in _children(ev, 'listOfEventAssignments'): + for ea in _children(lst, 'eventAssignment'): + var = ea.get('variable') + if var: + event_targets.add(var) elif ctag == 'listOfRules': # An makes X an algebraically-computed entity, not a # simulation output -- record its RHS (serialized to PEtab-math infix) so it is @@ -114,10 +231,59 @@ def parse_model(text): for e in _children(container, 'assignmentRule'): var = e.get('variable') if var: - assignment_rules[var] = _assignment_rule_formula(e) - + assignment_rules[var] = _math_formula(e) + for e in _children(container, 'rateRule'): + var = e.get('variable') + if var: + rate_rule_targets.add(var) + for e in _children(container, 'algebraicRule'): + # An algebraic rule determines one of the symbols in its math and the scan + # cannot tell which, so none of them can be trusted to hold still. + algebraic_symbols |= _mathml_identifiers(_expression_node(e)) + + # Settle the initial assignments. One that is arithmetic over numbers alone IS the entity's + # value. One computed from other entities means the model file settles no value at all, so + # the entity is dropped and its definition recorded for the measurement layer to inline -- + # but only when every entity it reads holds still for the whole simulation, because the + # substitution is read at a measurement time and an initial assignment is not (#795). + time_varying = (set(species) | set(assignment_rules) | rate_rule_targets + | algebraic_symbols | event_targets | non_constant) + invariant = (set(parameters) | set(compartments)) - time_varying + derived_initial_values, derived_refusals = {}, {} + for sym, elem in initial_elems.items(): + value = _evaluate_mathml_number(elem) + if sym in species: + # A species stays an output column either way: an initial assignment pins t=0, it + # does not make the species algebraic. It only loses a stale declared initial. + if value is None: + species_initial.pop(sym, None) + else: + species_initial[sym] = value + continue + if sym not in parameters and sym not in compartments: + continue # a stoichiometry target, or an id this scan does not collect + table = parameters if sym in parameters else compartments + if value is not None: + table[sym] = value + continue + table[sym] = None + rhs = _math_formula(elem) + deps = _mathml_identifiers(_expression_node(elem)) + if rhs is None: + derived_initial_values[sym] = None + derived_refusals[sym] = _R_INITIAL_UNTRANSLATABLE + elif sym not in invariant or not deps <= invariant: + derived_initial_values[sym] = None + derived_refusals[sym] = _r_time_varying(sorted((deps | {sym}) - invariant)) + else: + derived_initial_values[sym] = rhs + + # An entity the model file defines in terms of others never reports a number, whichever + # construct defines it. A rule target that also carries a vestigial `value` attribute is the + # same defect as the initial-assignment case: the rule overwrites it at t=0 anyway. + derived_all = set(assignment_rules) | set(derived_initial_values) parameter_values = {k: v for k, v in {**parameters, **compartments}.items() - if v is not None} + if v is not None and k not in derived_all} return SbmlEntities( text=text, species_names=frozenset(species), @@ -126,6 +292,8 @@ def parse_model(text): species_initial=species_initial, parameter_values=parameter_values, assignment_rules=assignment_rules, + derived_initial_values=derived_initial_values, + derived_refusals=derived_refusals, ) @@ -157,35 +325,148 @@ class _UnsupportedMathML(Exception): # ``^`` for exponentiation). ``minus`` is special-cased (unary negation vs binary subtraction). _MATHML_NARY = {'plus': ' + ', 'times': ' * '} _MATHML_BINARY = {'divide': ' / ', 'power': ' ^ '} +# MathML function element -> the python callable that evaluates it (#795). The serializer's +# accepted function set is derived from this one below, so the printed and the evaluated readings +# of a model can never drift apart. +_EVAL_FUNCS = {'exp': math.exp, 'ln': math.log, 'sqrt': math.sqrt, 'abs': abs, + 'sin': math.sin, 'cos': math.cos, 'tan': math.tan} # MathML function element -> a PEtab-math function call ``name(arg, ...)``. A conservative set # the petab grammar parses unambiguously; anything else is _UnsupportedMathML (-> None RHS), so # an exotic rule defers to a clear error instead of risking a wrong inline (#465 / ADR-0035). -_MATHML_FUNCS = frozenset({'exp', 'ln', 'sqrt', 'abs', 'sin', 'cos', 'tan'}) +_MATHML_FUNCS = frozenset(_EVAL_FUNCS) # The operator applications that bind looser than a surrounding operator and so are wrapped in # parens when used as an operand (a function call / atom is already self-delimiting). _MATHML_OPERATORS = frozenset(_MATHML_NARY) | frozenset(_MATHML_BINARY) | {'minus'} -def _assignment_rule_formula(rule_elem): - """An ````'s defining ```` serialized to a PEtab-math infix string, - or ``None`` if its MathML uses a construct :func:`_serialize_mathml` does not translate. +def _expression_node(elem): + """The expression node inside ``elem``'s ```` child, or ``None``. + + Skips a MathML ````/```` sibling of the expression. Shared by the + serializer and the evaluator so both read the same node of an ```` (#465) or + an ```` (#795). + """ + math_elem = _find_child(elem, 'math') + if math_elem is None: + return None + for child in math_elem: + if _local(child.tag) == 'annotation': + continue + return child + return None + + +def _mathml_identifiers(node): + """Every ```` name under ``node``, empty when ``node`` is ``None``. + + What an expression reads, used to decide whether an initial assignment may be inlined into a + measurement formula (#795). + """ + if node is None: + return set() + names = set() + if _local(node.tag) == 'ci': + text = (node.text or '').strip() + if text: + names.add(text) + for child in node: + names |= _mathml_identifiers(child) + return names + + +def _is_constant(elem): + """Whether a ````/```` is declared constant. + + The attribute is required in SBML L3 and optional in L2, where it defaults to true. + """ + return (elem.get('constant') or 'true').strip() in ('true', '1') + + +def _evaluate_mathml_number(elem): + """``elem``'s defining ```` as a float, or ``None`` when it is not arithmetic over + numbers alone. + + This is how the scan tells a self-contained definition from a derived one (#795). Walking the + tree and refusing a ```` makes "evaluates to a number" and "reads no other entity" the + same test, and it is why the check is not "serialize it and try ``float``": the serializer + renders ``13`` as ``(1 / 3)``, which ``float`` rejects, and + ``float`` accepts ``inf``, which a model's ```` could spell. + """ + node = _expression_node(elem) + if node is None: + return None + try: + value = _eval_mathml(node) + except _UnsupportedMathML: + return None + return value if math.isfinite(value) else None + + +def _eval_mathml(node): + """A content-MathML expression ``node`` -> a float, raising :class:`_UnsupportedMathML` on + anything that is not arithmetic over numeric literals -- an identifier included.""" + tag = _local(node.tag) + if tag == 'cn': + return _cn_value(node) + if tag == 'apply': + return _eval_apply(node) + raise _UnsupportedMathML # a lands here: an identifier is not a number + + +def _eval_apply(node): + """An ```` evaluated over the same grammar :func:`_serialize_apply` prints.""" + children = list(node) + if not children: + raise _UnsupportedMathML + op = _local(children[0].tag) + operands = children[1:] + try: + if op in ('plus', 'times'): + if not operands: + raise _UnsupportedMathML + total = 0.0 if op == 'plus' else 1.0 + for operand in operands: + if op == 'plus': + total += _eval_mathml(operand) + else: + total *= _eval_mathml(operand) + return total + if op == 'minus': + if len(operands) == 1: + return -_eval_mathml(operands[0]) + if len(operands) == 2: + return _eval_mathml(operands[0]) - _eval_mathml(operands[1]) + raise _UnsupportedMathML + if op in ('divide', 'power'): + if len(operands) != 2: + raise _UnsupportedMathML + left, right = _eval_mathml(operands[0]), _eval_mathml(operands[1]) + return left / right if op == 'divide' else left ** right + if op in _EVAL_FUNCS: + return float(_EVAL_FUNCS[op](*[_eval_mathml(a) for a in operands])) + except (ArithmeticError, ValueError, TypeError) as e: + raise _UnsupportedMathML from e + raise _UnsupportedMathML + + +def _math_formula(elem): + """An element's defining ```` serialized to a PEtab-math infix string, or + ``None`` if its MathML uses a construct :func:`_serialize_mathml` does not translate. + + Two callers: an ```` (#465) and an ```` (#795). The stdlib (libsbml-free) counterpart of a MathML pretty-printer: enough of content MathML to carry the algebraic convenience observables SBML authors actually write (the D2D ``Epo_cells := Epo_EpoRi + dEpoi``) into the measurement layer, where the loader inlines it down to species (#465). The result feeds the PEtab-math parser + a round-trip self-check at inline time, so any serialization defect is caught loudly there, never silently scored.""" - math = _find_child(rule_elem, 'math') - if math is None: + node = _expression_node(elem) + if node is None: + return None + try: + return _serialize_mathml(node) + except _UnsupportedMathML: return None - for child in math: - if _local(child.tag) == 'annotation': - continue # skip a MathML / sibling of the expression - try: - return _serialize_mathml(child) - except _UnsupportedMathML: - return None - return None def _serialize_mathml(node): @@ -241,9 +522,12 @@ def _operand(node): return text -def _serialize_cn(node): - """A MathML ```` numeric literal -> its infix spelling, honoring ``e-notation`` (a - ````-split mantissa/exponent) and ``rational`` (a ````-split numerator/denom).""" +def _cn_parts(node): + """A MathML ````'s ``(type, [number parts])``, honoring the ```` split. + + One reading of a literal for both the serializer and the evaluator, so ``13`` cannot be + printed as one thing and evaluated as another. + """ ctype = (node.get('type') or 'real').strip() nums = [] if node.text and node.text.strip(): @@ -253,6 +537,13 @@ def _serialize_cn(node): nums.append(child.tail.strip()) if not nums: raise _UnsupportedMathML + return ctype, nums + + +def _serialize_cn(node): + """A MathML ```` numeric literal -> its infix spelling, honoring ``e-notation`` (a + ````-split mantissa/exponent) and ``rational`` (a ````-split numerator/denom).""" + ctype, nums = _cn_parts(node) if ctype == 'e-notation' and len(nums) == 2: return f'{nums[0]}e{nums[1]}' if ctype == 'rational' and len(nums) == 2: @@ -260,6 +551,19 @@ def _serialize_cn(node): return nums[0] +def _cn_value(node): + """A MathML ```` numeric literal -> its float value (#795).""" + ctype, nums = _cn_parts(node) + try: + if ctype == 'e-notation' and len(nums) == 2: + return float(f'{nums[0]}e{nums[1]}') + if ctype == 'rational' and len(nums) == 2: + return float(nums[0]) / float(nums[1]) + return float(nums[0]) + except (ValueError, ArithmeticError) as e: + raise _UnsupportedMathML from e + + def _species_initial(species_elem): """A species' initial value (``initialAmount`` or ``initialConcentration``), or None.""" for attr in ('initialAmount', 'initialConcentration'): diff --git a/pybnf/petab/conditions.py b/pybnf/petab/conditions.py index 1567bb75..aad30fe5 100644 --- a/pybnf/petab/conditions.py +++ b/pybnf/petab/conditions.py @@ -146,7 +146,7 @@ def _species_target_value(pattern, op, val): # Asset: one mutation's operator -> a PEtab targetValue string # --------------------------------------------------------------------------- -def mutation_target_value(op, val, *, nominal=None, surrogate=None): +def mutation_target_value(op, val, *, nominal=None, surrogate=None, target=None): """Map one PyBNF mutation `` `` to a PEtab ``targetValue`` string. An absolute set (``=``) is the bare number, regardless of target kind. A relative op @@ -181,11 +181,19 @@ def mutation_target_value(op, val, *, nominal=None, surrogate=None): if surrogate is not None: return f'{surrogate} {op} {num(val)}' if nominal is None: + named = f"the fixed parameter '{target}'" if target else 'a fixed parameter' + fix = (f"Write this condition as an absolute set ('= '), or declare '{target}' as " + f"a fit parameter so the condition can be expressed relative to its surrogate" + if target else + "Write this condition as an absolute set ('= ')") raise NotImplementedError( - f"A relative mutation ('{op}' {num(val)}) of a fixed parameter needs the " - f"parameter's numeric nominal value, but it has a non-numeric (expression) " - f"value in the model. Evaluating a BNGL parameter expression is " - f"simulation-grade work, out of scope for the exporter (ADR-0027).") + f"A relative mutation ('{op}' {num(val)}) of {named} needs that parameter's numeric " + f"nominal value, but the model file does not settle one. In a BNGL model its value is " + f"an expression over other parameters. In an SBML model an assignment rule computes it " + f"every step, or an initial assignment derives it at the start of the simulation. " + f"PyBNF does not report a value the model file alone cannot settle, because a PEtab " + f"parameter table may override or estimate the entities it is computed from, and the " + f"reported number would then be stale. {fix} (ADR-0027, #465, #795).") if op == '*': return num(nominal * val) if op == '/': @@ -234,7 +242,7 @@ def _condition_rows_for(cid, perturbations, surrogate, nominal_of, species_id_of cid, species_id_of[var], _species_target_value(var, op, val))) continue rows.append(PetabConditionRow( - cid, var, mutation_target_value(op, val, nominal=nominal_of(var)))) + cid, var, mutation_target_value(op, val, nominal=nominal_of(var), target=var))) return rows diff --git a/pybnf/petab/formula.py b/pybnf/petab/formula.py index 027f5d3b..0d84a7b8 100644 --- a/pybnf/petab/formula.py +++ b/pybnf/petab/formula.py @@ -182,9 +182,9 @@ def substitute_placeholders(formula, substitutions): return petab_math -def inline_assignment_rules(formula, rules, *, observable_id=None): - """Inline the SBML assignment-rule variables a measurement ``observableFormula`` references - down to the species/parameters their rules are defined over (#465, ADR-0036). +def inline_derived_symbols(formula, derived, *, observable_id=None): + """Inline the derived SBML entities a measurement ``observableFormula`` references down to the + species/parameters they are defined over (#465 and #795, ADR-0036). An SBML ``assignmentRule variable="X"`` makes ``X`` an algebraic function of other model entities, recomputed every step -- never a simulation-output column and value-less, so the @@ -197,48 +197,65 @@ def inline_assignment_rules(formula, rules, *, observable_id=None): :func:`compile_petab_formula` + the chain-rule gradient (:meth:`MeasurementModel.\ prediction_sensitivity`) handle unchanged: the rule's species sensitivities flow in automatically. - ``rules`` maps every assignment-rule target id to its RHS as a PEtab-math infix string (the - stdlib ``_sbml`` scanner's serialization), or to ``None`` when the rule's MathML used a - construct the scanner could not translate. A formula that references **no** rule variable is - returned **verbatim** (the common case never reaches sympy, so it stays byte-stable). The - substitution + serialization go through ``sympy`` (never a string tokenizer -- ADR-0033), + The second construct is an ``initialAssignment`` on a parameter or compartment (#795). SBML + lets one supersede the declared ``value``/``size``, so the attribute is a placeholder and the + entity has no value of its own. Inlining is what makes such a model work rather than merely + fail politely: an observable over Bertozzi's ``beta_N``, defined as ``(R0_ * gamma_) / N_``, + resolves to an expression over the three parameters, and if the fit estimates one of them the + measurement tracks the fit, because :meth:`MeasurementModel.materialize` resolves a symbol + from the parameter set before the constant snapshot. An initial assignment is only inlinable + when every entity it reads holds still for the whole simulation, since the substituted + expression is read at a measurement time while the assignment fixed a value at t=0; the + scanner applies that gate and records a non-inlinable definition as ``None``. + + ``derived`` maps every derived entity id to a ``DerivedSymbol`` carrying its ``kind``, its RHS + as a PEtab-math infix string (the stdlib ``_sbml`` scanner's serialization) or ``None``, and + the clause naming why a ``None`` cannot be inlined. A formula that references **no** derived + entity is returned **verbatim** (the common case never reaches sympy, so it stays byte-stable). + The substitution + serialization go through ``sympy`` (never a string tokenizer -- ADR-0033), guarded by the same numeric round-trip self-check as the exporter. Raises ``PybnfError`` on a missing ``petab`` extra, an unparseable formula/RHS, a reference to - a rule whose MathML was not translatable (``None`` RHS), a cyclic rule dependency, or a + a definition the scanner could not carry (``None`` RHS), a cyclic dependency, or a substitution that fails its round-trip self-check. """ - if not rules: + if not derived: return formula sympify_petab = _require_petab_math() expr = _parse(sympify_petab, formula, source='observableFormula') - referenced = {str(s) for s in expr.free_symbols} & set(rules) + referenced = {str(s) for s in expr.free_symbols} & set(derived) if not referenced: - return formula # no assignment-rule variable in the formula -> carry it verbatim + return formula # no derived entity in the formula -> carry it verbatim where = f"Measurement model '{observable_id}': " if observable_id else 'Measurement model: ' - resolved = {} # rule name -> its fully rule-free sympy expr (memoized) + resolved = {} # name -> its fully derived-free sympy expr (memoized) def resolve(name, path): if name in resolved: return resolved[name] + symbol = derived.get(name) + noun = ('initial assignment' if getattr(symbol, 'kind', None) == 'initial_assignment' + else 'assignment rule') if name in path: chain = ' -> '.join(path[path.index(name):] + (name,)) raise PybnfError( - f"{where}the SBML assignment rule for '{name}' is cyclic ({chain}), so it cannot " - f"be resolved to a closed-form observable. (ADR-0036, #465.)") - rhs_infix = rules.get(name) + f"{where}the SBML {noun} for '{name}' is cyclic ({chain}), so it cannot " + f"be resolved to a closed-form observable. (ADR-0036, #465 and #795.)") + rhs_infix = symbol.expression if symbol is not None else None if rhs_infix is None: + reason = getattr(symbol, 'refusal', None) or ( + f'its defining {noun} uses a construct this reader does not translate') raise PybnfError( - f"{where}the observableFormula references the SBML assignment-rule variable " - f"'{name}', whose defining math uses a construct PyBNF's measurement layer cannot " - f"inline (e.g. a piecewise/relational rule). Write the observable directly over " - f"the species/parameters it is computed from instead. (ADR-0036, #465.)") - rhs = _parse(sympify_petab, rhs_infix, source=f"assignment rule for {name!r}") + f"{where}the observableFormula references '{name}', which the SBML model file " + f"defines in terms of other model entities and gives no value of its own. PyBNF " + f"cannot inline that definition here because {reason}. Write the observable " + f"directly over the species and parameters it is computed from instead. " + f"(ADR-0036, #465 and #795.)") + rhs = _parse(sympify_petab, rhs_infix, source=f"{noun} for {name!r}") # Substitute by the ACTUAL free-symbol objects of this RHS parse (petab tags symbols with # assumptions, so a plain sp.Symbol(name) would be a distinct object subs/diff ignore -- # the same gotcha inline_constants / substitute_placeholders guard). sub = {s: resolve(str(s), path + (name,)) - for s in rhs.free_symbols if str(s) in rules} + for s in rhs.free_symbols if str(s) in derived} resolved[name] = rhs.subs(sub) if sub else rhs return resolved[name] diff --git a/pybnf/petab/import_.py b/pybnf/petab/import_.py index 9d9101d7..fc97e54e 100644 --- a/pybnf/petab/import_.py +++ b/pybnf/petab/import_.py @@ -288,7 +288,7 @@ def import_job(problem_yaml_path, out_dir, job_type='de', method='ode', # post-sim observation layer, never a model-file edit (ADR-0036). model_texts = {} # location -> verbatim text namespaces, entity_name_sets = [], [] - assignment_rules = {} # SBML assignmentRule target -> RHS infix (inlined below, #493) + derived = {} # SBML entity defined by others -> DerivedSymbol (inlined below, #493/#795) for m in models: loc, lang = m['location'], (m['language'] or 'bngl').lower() text = (base / loc).read_text(encoding='utf-8', errors='replace') @@ -296,7 +296,7 @@ def import_job(problem_yaml_path, out_dir, job_type='de', method='ode', ns, ents, rules = _model_namespace(text, lang) namespaces.append(ns) entity_name_sets.append(ents) - assignment_rules.update(rules) + derived.update(rules) namespace = set().union(*namespaces) entity_names = set().union(*entity_name_sets) @@ -305,7 +305,7 @@ def import_job(problem_yaml_path, out_dir, job_type='de', method='ode', # observableFormulas (ADR-0036: emitted as conf `observable: ... formula:` lines). observable_id_to_column, measurement_models = _observable_id_to_column( observable_rows, namespace, entity_names, fixed_params, obs_params, free_names, - row_varying_obs_params, assignment_rules) + row_varying_obs_params, derived) # Pre-equilibrated dose-response reconstruction (ADR-0062): pull out the two-period scan groups # (a -inf pre-equilibration period + a per-dose measurement period) FIRST, so the plain @@ -504,21 +504,21 @@ def _model_namespace(model_text, language): """The model's expression namespace + entity name set + assignment rules, per language (ADR-0036). - Returns ``(namespace_symbols, entity_names, assignment_rules)``: ``namespace_symbols`` are + Returns ``(namespace_symbols, entity_names, derived)``: ``namespace_symbols`` are the names an ``observableFormula`` may reference (the BNGL ``ParamList`` -- parameters u observables u functions; or SBML species u parameters u compartments -- ADR-0026/0036); ``entity_names`` is the broader declared-name set used for the shadow check (a measurement - model's id must not collide with a model output column); ``assignment_rules`` maps each SBML - ``assignmentRule`` target id to its rule's RHS as a PEtab-math infix string (``None`` if the - MathML was not translatable), the map the importer **inlines** so a formula naming a rule - variable resolves down to the species/parameters the rule is computed from (#493, the import - peer of the config-load inlining -- #465). It is ``{}`` for a BNGL model (no assignment - rules). Read from the model text directly with the stdlib scanners (``_bngl`` / ``_sbml``), - simulator-free. + model's id must not collide with a model output column); ``derived`` maps each SBML entity + the model file defines in terms of others to a ``DerivedSymbol`` -- an ``assignmentRule`` + target (#465/#493) and a parameter or compartment an ``initialAssignment`` derives (#795) -- + the map the importer **inlines** so a formula naming one resolves down to the + species/parameters it is computed from, the import peer of the config-load inlining. It is + ``{}`` for a BNGL model. Read from the model text directly with the stdlib scanners + (``_bngl`` / ``_sbml``), simulator-free. """ if language == 'sbml': ent = parse_sbml_model(model_text) - return ent.namespace_symbols, set(ent.namespace_symbols), ent.assignment_rules + return ent.namespace_symbols, set(ent.namespace_symbols), ent.derived_symbols ent = parse_bngl_model(model_text) namespace = (set(ent.parameters) | set(ent.observable_names) | set(ent.function_names)) @@ -527,7 +527,7 @@ def _model_namespace(model_text, language): return namespace, entity_names, {} -def _shared_bare_entities(observable_rows, namespace, assignment_rules, row_varying_obs_params): +def _shared_bare_entities(observable_rows, namespace, derived, row_varying_obs_params): """Model entities named as a bare ``observableFormula`` by more than one observable (#503). A bare-name observable (no observableParameters placeholder, formula an identifier in the @@ -543,10 +543,10 @@ def _shared_bare_entities(observable_rows, namespace, assignment_rules, row_vary counts = Counter() for row in observable_rows: raw = (row.observable_formula or '').strip() - if assignment_rules: - from .formula import inline_assignment_rules - raw = inline_assignment_rules( - raw, assignment_rules, observable_id=row.observable_id) + if derived: + from .formula import inline_derived_symbols + raw = inline_derived_symbols( + raw, derived, observable_id=row.observable_id) if _PLACEHOLDER.search(raw) or row.observable_id in row_varying_obs_params: continue if _IDENTIFIER.match(raw) and raw in namespace: @@ -556,7 +556,7 @@ def _shared_bare_entities(observable_rows, namespace, assignment_rules, row_vary def _observable_id_to_column(observable_rows, namespace, entity_names, fixed_params, obs_params, free_names, row_varying_obs_params=(), - assignment_rules=None): + derived=None): """Map each ``observableId`` to the model column it measures, recording a measurement model for any expression ``observableFormula`` (ADR-0036). Iteration order = table order, which fixes the wide-data column order on the measurement pivot. @@ -593,18 +593,19 @@ def _observable_id_to_column(observable_rows, namespace, entity_names, fixed_par validate against the model namespace. A *constant*-per-observable placeholder is substituted away as in Phase 1; an unresolved (neither constant nor row-varying) placeholder still raises. - ``assignment_rules`` (#493) is the SBML ``assignmentRule`` map (target id -> rule RHS as - PEtab-math infix). An assignment-rule variable is a *derived* model output -- declared as a - parameter/species but algebraically computed every step, so it is never a simulation-output - column and cannot be resolved *as a symbol* (that is exactly why it is absent from - ``namespace``). Each observableFormula is therefore **inlined** first: every referenced rule - variable is replaced by its rule's RHS (recursively) so the formula reduces to species / - parameters the layer can evaluate -- the SBML analogue of a BNGL global function in an - observableFormula, which PyBNF already accepts. A formula naming no rule variable is returned - verbatim (the bare-name common case stays dependency-free); this mirrors the config-load - inlining (#465). + ``derived`` (#493/#795) maps an SBML entity the model file defines in terms of others to a + ``DerivedSymbol``. Such an entity is declared as a parameter/species but has no value of its + own -- an ``assignmentRule`` target is recomputed every step, and a parameter an + ``initialAssignment`` derives has only a placeholder attribute -- so it is never a + simulation-output column and cannot be resolved *as a symbol* (that is exactly why it is + absent from ``namespace``). Each observableFormula is therefore **inlined** first: every + referenced entity is replaced by its defining RHS (recursively) so the formula reduces to + species and parameters the layer can evaluate -- the SBML analogue of a BNGL global function + in an observableFormula, which PyBNF already accepts. A formula naming no derived entity is + returned verbatim (the bare-name common case stays dependency-free); this mirrors the + config-load inlining (#465/#795). """ - assignment_rules = assignment_rules or {} + derived = derived or {} # Fixed PEtab parameters that are NOT model entities are inlined as literals; one that # IS a model entity stays a symbol (it resolves as a model constant at eval time). inline = {n: v for n, v in fixed_params.items() if n not in namespace} @@ -620,7 +621,7 @@ def _observable_id_to_column(observable_rows, namespace, entity_names, fixed_par # columns below so their per-observable noise_model overrides key by distinct ids, not the # one shared column (#503, ADR-0077). A uniquely-targeted entity keeps the bare-name path. shared_entities = _shared_bare_entities( - observable_rows, namespace, assignment_rules, row_varying_obs_params) + observable_rows, namespace, derived, row_varying_obs_params) for row in observable_rows: raw = (row.observable_formula or '').strip() # Inline any SBML assignment-rule variable the formula names down to the species / @@ -630,10 +631,10 @@ def _observable_id_to_column(observable_rows, namespace, entity_names, fixed_par # the bare-name common case never reaches the translator (dependency-free, byte-stable); # the inlining leaves any observableParameters placeholder untouched (rules are model # MathML, never placeholders), so the placeholder handling below is unchanged. - if assignment_rules: - from .formula import inline_assignment_rules - raw = inline_assignment_rules( - raw, assignment_rules, observable_id=row.observable_id) + if derived: + from .formula import inline_derived_symbols + raw = inline_derived_symbols( + raw, derived, observable_id=row.observable_id) had_placeholder = bool(_PLACEHOLDER.search(raw)) row_varying = row.observable_id in row_varying_obs_params if row_varying: diff --git a/tests/test_petab_export.py b/tests/test_petab_export.py index 81e18f80..b4cd90bb 100644 --- a/tests/test_petab_export.py +++ b/tests/test_petab_export.py @@ -2738,8 +2738,50 @@ def test_relative_op_on_fixed_target_is_precomputed(self, op, val, nominal, expe assert mutation_target_value(op, val, nominal=nominal) == expected def test_relative_op_on_expression_nominal_raises(self): - with pytest.raises(NotImplementedError): + with pytest.raises(NotImplementedError) as excinfo: mutation_target_value('*', 2.0, nominal=None) + message = str(excinfo.value) + # The refusal covers both languages: a BNGL parameter whose value is an expression, and + # an SBML parameter a rule or an initial assignment defines (#795). It used to name BNGL + # alone, though an SBML model could already reach it. + assert 'BNGL' in message and 'SBML' in message + assert 'initial assignment' in message + + def test_the_refusal_names_the_parameter_when_the_caller_knows_it(self): + with pytest.raises(NotImplementedError, match="'beta_N'"): + mutation_target_value('*', 2.0, nominal=None, target='beta_N') + + def test_a_relative_condition_on_a_derived_sbml_parameter_refuses(self, tmp_path): + # End to end through the exporter's own reader: the model file settles no value for + # beta_N, so the exporter refuses rather than writing the placeholder attribute times two + # into conditions.tsv, which is what issue #795 reported. + from pybnf.petab.conditions import mutation_target_value as target_value + from pybnf.petab.export import _numeric_nominal, _read_model + xml = tmp_path / 'derived.xml' + xml.write_text( + ''' + + + + + + + + + + + koffKd + + + + + +''') + view = _read_model(xml.read_text(), xml, 'sbml') + assert _numeric_nominal(view, 'koff') == 1.0 # a literal still has a nominal + assert _numeric_nominal(view, 'kon') is None # the placeholder value="0" is gone + with pytest.raises(NotImplementedError, match="'kon'"): + target_value('*', 2.0, nominal=_numeric_nominal(view, 'kon'), target='kon') def test_surrogate_marker_is_double_underscore(self): assert surrogate_name('v1') == 'v1__REF' diff --git a/tests/test_petab_formula.py b/tests/test_petab_formula.py index 354fe04b..03e5b581 100644 --- a/tests/test_petab_formula.py +++ b/tests/test_petab_formula.py @@ -34,6 +34,7 @@ """ import builtins +import re import sys import pytest @@ -395,3 +396,78 @@ def test_formula_free_symbols_lists_sorted_names(self): from pybnf.petab.formula import formula_free_symbols assert formula_free_symbols('0.1 + 0.05*slope + base') == ['base', 'slope'] assert formula_free_symbols('0.5') == [] # a pure constant has no free symbols + + +class TestDerivedSymbolInlining: + """Inlining a derived SBML entity into a measurement formula (#465 for an assignment rule, + #795 for a parameter an initialAssignment derives).""" + + @staticmethod + def _derived(**kwargs): + from pybnf.petab._sbml import DerivedSymbol + return {name: DerivedSymbol('initial_assignment', expr, None) + for name, expr in kwargs.items()} + + def test_a_derived_parameter_inlines_to_its_definition(self): + pytest.importorskip('petab') + from pybnf.petab.formula import inline_derived_symbols + # Bertozzi's shape: the model file gives beta_N no value of its own, so scoring an + # observable over it used to use the placeholder attribute. + out = inline_derived_symbols('beta_N * I', + self._derived(beta_N='(R0_ * gamma_) / N_'), + observable_id='rate') + assert 'beta_N' not in out + assert {'R0_', 'gamma_', 'N_', 'I'} <= set(re.findall(r'[A-Za-z_]\w*', out)) + + def test_an_alias_inlines(self): + pytest.importorskip('petab') + from pybnf.petab.formula import inline_derived_symbols + # Laske's shape: 27 parameters of the form ModelValue_82 = D_rib. + assert inline_derived_symbols('ModelValue_82', + self._derived(ModelValue_82='D_rib')) == 'D_rib' + + def test_a_chain_resolves_through(self): + pytest.importorskip('petab') + from pybnf.petab.formula import inline_derived_symbols + out = inline_derived_symbols('chain', self._derived(chain='stale * 2', stale='k + 1')) + assert 'stale' not in out and 'chain' not in out + assert 'k' in out + + def test_a_formula_naming_nothing_derived_is_returned_verbatim(self): + pytest.importorskip('petab') + from pybnf.petab.formula import inline_derived_symbols + assert inline_derived_symbols('A + B', self._derived(x='k + 1')) == 'A + B' + + def test_a_non_inlinable_initial_assignment_names_the_real_cause(self): + pytest.importorskip('petab') + from pybnf.petab._sbml import DerivedSymbol, _r_time_varying + from pybnf.petab.formula import inline_derived_symbols + from pybnf.printing import PybnfError + derived = {'over_species': DerivedSymbol('initial_assignment', None, + _r_time_varying(['A']))} + with pytest.raises(PybnfError) as excinfo: + inline_derived_symbols('over_species * 2', derived, observable_id='rate') + message = str(excinfo.value) + assert "Measurement model 'rate'" in message + assert 'gives no value of its own' in message + assert "'A'" in message and 'changes during the simulation' in message + assert '#795' in message + + def test_an_untranslatable_assignment_rule_still_says_assignment_rule(self): + pytest.importorskip('petab') + from pybnf.petab._sbml import DerivedSymbol, _R_RULE_UNTRANSLATABLE + from pybnf.petab.formula import inline_derived_symbols + from pybnf.printing import PybnfError + derived = {'ruled': DerivedSymbol('assignment_rule', None, _R_RULE_UNTRANSLATABLE)} + with pytest.raises(PybnfError, match='assignment rule'): + inline_derived_symbols('ruled', derived, observable_id='rate') + + def test_a_cyclic_definition_names_the_construct_and_the_chain(self): + pytest.importorskip('petab') + from pybnf.petab.formula import inline_derived_symbols + from pybnf.printing import PybnfError + with pytest.raises(PybnfError) as excinfo: + inline_derived_symbols('a', self._derived(a='b', b='a'), observable_id='rate') + message = str(excinfo.value) + assert 'initial assignment' in message + assert 'a -> b -> a' in message diff --git a/tests/test_petab_sbml_layer.py b/tests/test_petab_sbml_layer.py index 5d8dfba4..2d1c40fd 100644 --- a/tests/test_petab_sbml_layer.py +++ b/tests/test_petab_sbml_layer.py @@ -737,7 +737,7 @@ class TestAssignmentRuleObservableMaterializesLikeSpecies: def test_inlined_rule_observable_equals_hand_species_column(self, tmp_path): pytest.importorskip('petab') from pybnf.petab._sbml import parse_model as parse_sbml - from pybnf.petab.formula import inline_assignment_rules + from pybnf.petab.formula import inline_derived_symbols from pybnf.pset import FreeParameter, PSet, SbmlModelNoTimeout, TimeCourse # RULE_SBML is DECAY_SBML + ``total := A + B``; simulate it on RoadRunner. @@ -746,7 +746,7 @@ def test_inlined_rule_observable_equals_hand_species_column(self, tmp_path): namespace = set(ent.namespace_symbols) # Inline ``total`` exactly as config._load_measurement_models does, then materialize it # alongside the hand-reconstructed species formula on the same trace. - inlined = inline_assignment_rules('total', ent.assignment_rules, observable_id='from_rule') + inlined = inline_derived_symbols('total', ent.derived_symbols, observable_id='from_rule') assert inlined == 'A + B' # the rule resolved to its species assert 'total' not in namespace # the target is never bound as a symbol @@ -885,3 +885,134 @@ def test_layer_agrees_across_roadrunner_and_bngsim_on_boehm(self, tmp_path): for oid in _BOEHM_OBS: on_rr_grid = np.interp(rr['time'], bg['time'], bg[oid]) np.testing.assert_allclose(on_rr_grid, rr[oid], rtol=1e-3, atol=1e-3) + + +# An SBML model whose rate constant is derived: `beta_N` carries a placeholder value="0.0" and +# an initialAssignment computing it from three parameters. This is Bertozzi_PNAS2020's shape, +# which the scan of the public PEtab benchmark collection found wrong today (#795). +DERIVED_SBML = """ + + + + + + + + + + + + + + + + + + + R0_gamma_ + N_ + + + + + + +""" + + +class TestDerivedParameterObservable: + """An observable over a parameter an ``initialAssignment`` derives (#795). + + The scanner drops the placeholder attribute, so the value is not silently wrong; the loader + inlines the definition, so the observable is not merely refused either.""" + + def _conf(self, tmp_path, formula, model_text=DERIVED_SBML, extra='', + objective='objective = chi_sq'): + import os + import textwrap + from pybnf import config as config_mod + from pybnf.parse import ploop + (tmp_path / 'derived.xml').write_text(model_text) + (tmp_path / 'tc.exp').write_text('# time\trate\trate_SD\n0\t1.0\t0.1\n1\t1.0\t0.1\n') + conf_text = textwrap.dedent(f"""\ + edition = 2 + job_type = trf + {objective} + sbml_backend = roadrunner + model: derived.xml + observable: rate, formula: {formula} + experiment: tc, data: tc.exp + uniform_var = R0_ 1 5 + population_size = 4 + max_iterations = 1 + verbosity = 0 + """) + extra + home = os.getcwd() + os.chdir(tmp_path) + try: + return config_mod.Configuration(ploop(conf_text.splitlines(keepends=True))) + finally: + os.chdir(home) + + def test_the_observable_inlines_to_the_entities_it_is_computed_from(self, tmp_path): + pytest.importorskip('petab') + conf = self._conf(tmp_path, 'beta_N * I') + mm = conf.obj.measurement.models[0] + assert 'beta_N' not in mm.formula # the placeholder value is never bound + assert {'R0_', 'gamma_', 'N_', 'I'} <= mm.allowed_symbols + assert 'beta_N' not in mm.allowed_symbols + + def test_the_materialized_column_is_the_value_the_model_starts_from(self, tmp_path): + # Before the fix the scanner reported beta_N = 0.0, so this column was all zeros. + pytest.importorskip('petab') + import numpy as np + from pybnf.data import Data + conf = self._conf(tmp_path, 'beta_N * I') + mm = conf.obj.measurement.models[0] + data = Data() + data.cols = {'time': 0, 'I': 1} + data.data = np.array([[0., 2.], [1., 4.]]) + got = mm.materialize(data, pset_values={'R0_': 3.0}) + np.testing.assert_allclose(got, (3.0 * 0.5 / 10.0) * np.array([2., 4.])) + assert not np.allclose(got, 0.0) + + def test_the_column_tracks_a_fitted_dependency(self, tmp_path): + # The property that makes the import-time view agree with the runtime (ADR-0094): R0_ is + # a free parameter, so the derived rate constant moves with the fit rather than staying + # at whatever the model file implied. + pytest.importorskip('petab') + import numpy as np + from pybnf.data import Data + conf = self._conf(tmp_path, 'beta_N * I') + mm = conf.obj.measurement.models[0] + data = Data() + data.cols = {'time': 0, 'I': 1} + data.data = np.array([[0., 2.]]) + low = mm.materialize(data, pset_values={'R0_': 1.0}) + high = mm.materialize(data, pset_values={'R0_': 4.0}) + np.testing.assert_allclose(high / low, 4.0) + + def test_an_assignment_over_something_that_moves_refuses_at_load(self, tmp_path): + pytest.importorskip('petab') + from pybnf.printing import PybnfError + moving = DERIVED_SBML.replace('N_\n ', + 'I\n ') + assert 'I' in moving # the assignment really reads the species now + with pytest.raises(PybnfError) as excinfo: + self._conf(tmp_path, 'beta_N * I', model_text=moving) + message = str(excinfo.value) + assert 'initial assignment' in message + assert 'changes during the simulation' in message + + def test_a_prediction_noise_formula_may_name_a_derived_parameter(self, tmp_path): + # The sibling layer reads the same model namespace, so it has to inline the same way. + # Without this a derived parameter works in `observable:` and is rejected as an unknown + # symbol in `sigma = prediction_formula`, which is a #465 gap the initial-assignment work + # would otherwise have widened. + pytest.importorskip('petab') + conf = self._conf( + tmp_path, 'I', + objective='noise_model = normal, sigma = prediction_formula beta_N * I') + formulas = [src.formula for _label, src in conf._prediction_noise_sources()] + assert formulas and all('beta_N' not in f for f in formulas) + assert any('R0_' in f for f in formulas) diff --git a/tests/test_petab_sbml_scanner.py b/tests/test_petab_sbml_scanner.py index 4966b361..1a55be8a 100644 --- a/tests/test_petab_sbml_scanner.py +++ b/tests/test_petab_sbml_scanner.py @@ -240,3 +240,303 @@ def test_untranslatable_construct_degrades_to_none(self): ent = parse_model(doc) assert ent.assignment_rules == {'x': None} assert 'x' not in ent.namespace_symbols + + +# A model exercising every shape at once (#795). SBML lets an initial +# assignment supersede a parameter's value, a compartment's size, and a species' initial +# amount/concentration, so the declared attribute is a placeholder the model never starts from. +# settled -- arithmetic over numbers alone, so the assignment IS the value +# stale -- computed from a constant parameter: no value, inlinable +# valueless -- the same, written the way antimony emits it (no value attribute at all) +# chain -- computed from another derived parameter +# dcomp -- a compartment, to show size is superseded the same way +# over_* -- one per kind of entity that moves during a simulation, none of them inlinable +# piecewise -- MathML this stdlib reader does not translate +# A -- a species by assignment, which stays an output column +SBML_INITIAL = """ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + k + + + k + + + + + + + 1 + + + + + + + + 23 + + + + + k1 + + + + + k1 + + + + + stale2 + + + + + cell4 + + + + A + + + moving + + + ruled + + + rated + + + evented + + + + 1 + + + + + 221 + + + + + k5 + + + + + +""" + + +class TestInitialAssignmentValues: + """An supersedes the declared attribute, so the attribute is a + placeholder the model never starts from (#795).""" + + def test_self_contained_assignment_supersedes_the_attribute(self): + ent = parse_model(SBML_INITIAL) + assert ent.parameter_values['settled'] == 6.0 # 2 * 3, not the value="0" + assert 'settled' in ent.namespace_symbols # it has a value, so it binds + assert 'settled' not in ent.derived_initial_values + + def test_a_derived_parameter_reports_no_value(self): + # The wrong-number symptom: the scanner used to hand out value="99" while the model + # starts from k + 1 == 5. + ent = parse_model(SBML_INITIAL) + assert 'stale' not in ent.parameter_values + assert 'stale' not in ent.constants + + def test_a_derived_parameter_is_recorded_with_its_expression(self): + ent = parse_model(SBML_INITIAL) + assert ent.derived_initial_values['stale'] == 'k + 1' + + def test_a_derived_parameter_leaves_the_namespace_but_not_the_declaration(self): + ent = parse_model(SBML_INITIAL) + assert 'stale' not in ent.namespace_symbols # not resolvable as a symbol + assert 'stale' in ent.parameter_names # the scan stays faithful to the file + assert 'k' in ent.namespace_symbols + + def test_a_value_less_derived_parameter_behaves_identically(self): + # What antimony emits for `k_derived = k_base + 1`: no value attribute at all. This is + # the shape that used to reach the measurement layer's "should be unreachable" branch. + ent = parse_model(SBML_INITIAL) + assert ent.derived_initial_values['valueless'] == 'k + 1' + assert 'valueless' not in ent.parameter_values + assert 'valueless' not in ent.namespace_symbols + + def test_a_derived_compartment_loses_its_size(self): + ent = parse_model(SBML_INITIAL) + assert 'dcomp' not in ent.parameter_values + assert ent.derived_initial_values['dcomp'] == 'cell * 4' + assert 'dcomp' in ent.compartment_names + + def test_a_chain_of_initial_assignments_is_recorded(self): + # `chain` reads `stale`, which is itself derived. Both are recorded, and the inliner + # resolves the chain. + ent = parse_model(SBML_INITIAL) + assert ent.derived_initial_values['chain'] == 'stale * 2' + + def test_block_order_does_not_matter(self): + # The initial assignments are settled after the whole container loop, so a document that + # declares them before the parameters reads the same. + reordered = SBML_INITIAL.replace('', '\n' + ' \n' + ' \n' + ' 23\n' + ' \n' + ' \n' + ' \n ', 1) + assert parse_model(reordered).parameter_values['settled'] == 6.0 + + def test_a_model_without_initial_assignments_has_none(self): + assert parse_model(SBML_L3).derived_initial_values == {} + + +class TestInitialAssignmentSoundnessGate: + """An initial assignment fixes a value from its inputs' *initial* values, so it may only be + inlined into a measurement formula when every input holds still (#795).""" + + def test_an_assignment_over_a_moving_entity_is_not_inlinable(self): + ent = parse_model(SBML_INITIAL) + for name, offender in (('over_species', 'A'), ('over_moving', 'moving'), + ('over_ruled', 'ruled'), ('over_rated', 'rated'), + ('over_evented', 'evented')): + assert ent.derived_initial_values[name] is None, name + assert offender in ent.derived_refusals[name], name + assert 'changes during the simulation' in ent.derived_refusals[name] + + def test_untranslatable_math_is_recorded_as_none(self): + ent = parse_model(SBML_INITIAL) + assert ent.derived_initial_values['piecewise'] is None + assert 'does not translate' in ent.derived_refusals['piecewise'] + + def test_an_algebraic_rule_makes_its_symbols_untrusted(self): + algebraic = SBML_INITIAL.replace( + '', + ' \n' + ' k\n' + ' \n ', 1) + ent = parse_model(algebraic) + assert ent.derived_initial_values['stale'] is None + assert "'k'" in ent.derived_refusals['stale'] + + +class TestInitialAssignmentSpecies: + """A species with an initial assignment is still a dynamical state and still an output + column, so it keeps its place in the namespace and only loses a stale declared initial.""" + + def test_a_species_stays_in_the_namespace(self): + ent = parse_model(SBML_INITIAL) + assert 'B' in ent.namespace_symbols + assert 'B' in ent.species_names + assert 'B' not in ent.derived_initial_values + + def test_a_derived_species_loses_its_stale_initial(self): + ent = parse_model(SBML_INITIAL) + assert 'B' not in ent.species_initial # the file says k * 5, not the declared 3 + + def test_a_self_contained_species_assignment_is_evaluated(self): + ent = parse_model(SBML_INITIAL) + assert ent.species_initial['A'] == 42.0 # 2 * 21, not initialConcentration="1" + + def test_boehm_species_lose_their_placeholder_initials(self): + # The only committed model with this shape: STAT5A and STAT5B carry + # initialConcentration="1" while the assignments set them from 207.6 * ratio. + from pathlib import Path + text = (Path(__file__).parent / 'petab_fixtures' / 'boehm_v2' + / 'model_Boehm_JProteomeRes2014.xml').read_text() + ent = parse_model(text) + assert 'STAT5A' not in ent.species_initial + assert 'STAT5B' not in ent.species_initial + assert 'STAT5A' in ent.namespace_symbols # still an output column + assert ent.species_initial['pApB'] == 0.0 # a literal initial is untouched + + +class TestDerivedSymbolMap: + """The one map the measurement and import layers inline through.""" + + def test_each_kind_is_labelled(self): + symbols = parse_model(SBML_INITIAL).derived_symbols + assert symbols['stale'].kind == 'initial_assignment' + assert symbols['ruled'].kind == 'assignment_rule' + + def test_a_refusal_travels_with_the_symbol(self): + symbols = parse_model(SBML_INITIAL).derived_symbols + assert symbols['stale'].refusal is None + assert 'changes during the simulation' in symbols['over_species'].refusal + + def test_an_assignment_rule_target_reports_no_value_either(self): + # The sibling of the same defect: a rule target that carries a vestigial value attribute + # still had that number reported, though the rule overwrites it at t=0 anyway. + vestigial = SBML_RULES.replace('', + '') + assert 'value="7"' in vestigial # the fixture really changed + ent = parse_model(vestigial) + assert 'ratio' not in ent.parameter_values + assert 'ratio' not in ent.namespace_symbols + + +class TestInitialAssignmentEvaluation: + """The numeric reading of MathML, which is how a self-contained assignment is told from a + derived one.""" + + def _value(self, math_xml, attr='value="0"'): + doc = SBML_INITIAL.replace( + ' \n' + ' \n' + ' 23\n' + ' \n' + ' \n', + f' \n' + f' {math_xml}\n' + f' \n', 1) + return parse_model(doc).parameter_values.get('settled') + + def test_e_notation_literal(self): + assert self._value('13') == 1000.0 + + def test_rational_literal(self): + assert self._value('14') == 0.25 + + def test_unary_minus_and_power(self): + assert self._value('2') == -2.0 + assert self._value('23') == 8.0 + + def test_a_function_call(self): + assert self._value('0') == 1.0 + + def test_division_by_zero_is_not_a_number(self): + # Not an exception: the entity simply has no value the file settles, and it is recorded + # as derived like any other. + assert self._value('10') is None