From 267ae5a0a1185421f3e3e576d616e84fccc04943 Mon Sep 17 00:00:00 2001 From: Bill Hlavacek Date: Tue, 22 Sep 2026 15:55:36 -0600 Subject: [PATCH] Read SBML initial assignments instead of the placeholder attribute they supersede SBML lets a model set an entity's starting value in two places. A parameter carries a value attribute, a compartment a size, and a species an initial amount or concentration, and a listOfInitialAssignments entry may supersede any of them. When it does, the attribute is a placeholder the model never starts from. The import-time scanner read only the attributes, so it reported numbers the model does not start from, and said nothing. Two consumers used them. The measurement layer binds them as constants, so an observable or noise formula naming such a parameter was scored with the placeholder. The exporter turns them into absolute numbers for a relative condition, so conditions.tsv got a wrong value. A model that writes no attribute at all, which is what antimony emits for a derived parameter, failed differently and not much better: the symbol stayed in the namespace with no value behind it, so a formula naming it died at materialize time saying it was neither a simulation output nor a model parameter, which is not true. An entity an initial assignment derives now has no value the model file settles, so the scanner reports none. An assignment that is arithmetic over numbers alone is the value and replaces the attribute. An assignment computed from other entities drops the entity from the value map and the formula namespace and records its right hand side, and the measurement layer inlines that right hand side, which is what it already does for an assignment rule. So an observable over Bertozzi's beta_N, which the model computes as (R0_ * gamma_) / N_ behind a value of 0.0, is scored with the value the model starts from, and it tracks the fit when R0_ is estimated. Inlining is refused when the assignment reads something that moves during the simulation, because the substituted expression is read at a measurement time while the assignment fixed a value at the start. A species is deliberately left alone beyond losing a stale initial: an initial assignment pins the start only, so the species is still a dynamical state and still an output column, and excluding it would break the Boehm fixture, whose observables name two species set this way. Two adjacent defects are fixed in the same pass because they are the same sentence. An assignment rule target that carries a leftover value attribute no longer reports that number either. A prediction noise formula now inlines the same map a measurement formula does, so the two layers no longer disagree about which symbols a model offers. The refusal in mutation_target_value no longer blames BNGL alone. It names the parameter, gives both languages' version of the cause, and says what to write instead. That message was already reachable for an SBML model before this change. Measured on 25 models from the public PEtab benchmark collection, 16 use initial assignments, 31 of them on a parameter or compartment, and 28 of those were reporting a stale value. No model committed to this repository has one on a parameter, so nothing here changes. ADR-0147 records the decision. --- ...easurement-layer-inlines-the-derivation.md | 103 +++++ pybnf/config.py | 62 +-- pybnf/petab/_sbml.py | 356 ++++++++++++++++-- pybnf/petab/conditions.py | 20 +- pybnf/petab/formula.py | 61 +-- pybnf/petab/import_.py | 69 ++-- tests/test_petab_export.py | 44 ++- tests/test_petab_formula.py | 76 ++++ tests/test_petab_sbml_layer.py | 135 ++++++- tests/test_petab_sbml_scanner.py | 300 +++++++++++++++ 10 files changed, 1109 insertions(+), 117 deletions(-) create mode 100644 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 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 000000000..b6b9c917f --- /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 e46db2335..c16e28471 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 f833d59d3..a39c2dc4e 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 1567bb75f..aad30fe51 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 027f5d3b0..0d84a7b83 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 9d9101d73..fc97e54e5 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 81e18f808..b4cd90bba 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 354fe04b3..03e5b5819 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 5d8dfba43..2d1c40fdb 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 4966b3611..1a55be8a9 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