Read SBML initial assignments instead of the placeholder attribute they supersede

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 `<listOfInitialAssignments>` 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 `<parameter id="k_derived" constant="true"/>` 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 `<ci>`, 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.

pybnf/config.py

@@ -2505,9 +2505,9 @@ def _load_measurement_models(self):
             ed, 2, "the 'observable: <id>, formula: <expr>' 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 <expr>' 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):