From 585729111e1d09618c6e8c196cc6f61b869ab72a Mon Sep 17 00:00:00 2001 From: Joel Eliason Date: Tue, 4 Aug 2026 17:56:49 -0600 Subject: [PATCH] fix: evaluate observable rules in SI, not in declared units A rule is written in the model's declared units and is dimensionally consistent only after conversion. V_T adds V_Tmin in mL to a cell count times vol_cell in micrometre^3; SimBiology reconciles that by tracking units, and the C++ emitter reconciles it by scaling every parameter on load. The Python emitter read the raw numbers, so it was summing incompatible quantities. On the PDAC model all 13 constants in the live closure carry a factor, spanning 1.66e-24 to 6.02e5. V_T came out as 2.16e13 rather than 33 mL, which drove stromal_fraction to exactly 1.0 for every patient: a calibration target that could never move and would read as total misfit. Inputs convert in, results convert back out, so the module takes and returns the model's own units, which is what the simulator writes and what a target reads. The factors are emitted as a visible SI_FACTOR table rather than folded into the arithmetic. Also stamps the emitter's version and a hash of its source into the header. The version alone does not move under an editable install, and an artifact that does not rebuild when its generator changes is the same drift this module exists to prevent. --- src/qsp_codegen/observables.py | 79 ++++++++++++++++++++++++++++++---- tests/test_observables.py | 68 +++++++++++++++++++++++++++-- 2 files changed, 134 insertions(+), 13 deletions(-) diff --git a/src/qsp_codegen/observables.py b/src/qsp_codegen/observables.py index cb97b6f..edabc8a 100644 --- a/src/qsp_codegen/observables.py +++ b/src/qsp_codegen/observables.py @@ -108,6 +108,7 @@ class Closure: constants: Dict[str, float] # fixed parameters and compartment sizes order: Tuple[str, ...] # rules to evaluate, dependencies first expressions: Dict[str, str] # rule name -> infix, as the C++ emitter sees it + si: Dict[str, float] # symbol -> SI factor, native * si = SI @property def derived(self) -> Tuple[str, ...]: @@ -115,20 +116,31 @@ def derived(self) -> Tuple[str, ...]: return tuple(s for s in self.requested if s in self.expressions) -def _classify(sbml) -> Tuple[Dict[str, str], Set[str], Dict[str, float]]: - """Rules, species, and everything with a fixed value. +def _classify(sbml) -> Tuple[Dict[str, str], Set[str], Dict[str, float], + Dict[str, float]]: + """Rules, species, everything with a fixed value, and the SI factors. A compartment carrying a rule is derived, not a constant: ``V_T`` is computed from the species it holds, and its listed size is only an initial value. + + The factors are what make the arithmetic mean anything. A rule is written in + the model's declared units and is dimensionally consistent only after + conversion: ``V_T`` adds ``V_Tmin`` in mL to a cell count times ``vol_cell`` + in micrometre^3, which SimBiology reconciles by tracking units and the C++ + emitter reconciles by scaling every parameter on load. Reading the raw + numbers instead silently evaluates a sum of incompatible quantities. """ rules = {r["variable_name"]: _rewrite_calls(r["expression"]) for r in sbml.assignment_rules} species = {s["name"] for s in sbml.species} fixed: Dict[str, float] = {p["name"]: p["value"] for p in sbml.parameters} + si: Dict[str, float] = {} + for entry in (*sbml.parameters, *sbml.species, *sbml.compartments): + si[entry["name"]] = float(sbml.get_si_factor(entry.get("units"))) for c in sbml.compartments: if c["name"] not in rules: fixed[c["name"]] = c["size"] - return rules, species, fixed + return rules, species, fixed, si def resolve(sbml, symbols: Iterable[str]) -> Closure: @@ -137,7 +149,7 @@ def resolve(sbml, symbols: Iterable[str]) -> Closure: Raises :class:`UnresolvedSymbol` rather than dropping anything: a symbol that silently vanishes is a readout that silently stops being computed. """ - rules, species, fixed = _classify(sbml) + rules, species, fixed, si = _classify(sbml) requested = tuple(sorted(set(symbols))) states: Set[str] = set() @@ -175,12 +187,23 @@ def visit(name: str) -> None: + ", ".join(sorted(unresolved)) ) + named = (*requested, *states, *constants, *order) + # Only the inputs are required to declare units. A rule variable that is not + # itself a declared component carries none, and 1.0 is then right rather than + # a guess; an undeclared *input* would silently rescale the whole expression. + missing = sorted(n for n in (*states, *constants) if n not in si) + if missing: + raise UnresolvedSymbol( + f"no SI factor for {len(missing)} input(s): {missing[:5]}. Without one " + f"the emitted module would evaluate the rules in mixed units." + ) return Closure( requested=requested, states=tuple(sorted(states)), constants=dict(sorted(constants.items())), order=tuple(order), expressions={n: rules[n] for n in order}, + si={n: si.get(n, 1.0) for n in named}, ) @@ -214,6 +237,21 @@ def _translate(expression: str, rename: Dict[str, str]) -> str: return _TOKEN.sub(lambda m: rename.get(m.group(0), m.group(0)), expression) +def _emitter_stamp() -> str: + """Version plus a hash of this file, so an editable install is still detectable. + + The version alone does not move when the emitter is edited in place, and a + generated artifact that does not rebuild when its generator changes is the + same drift this module exists to prevent. + """ + import hashlib + + from . import __version__ + + digest = hashlib.sha256(Path(__file__).read_bytes()).hexdigest()[:12] + return f"{__version__}+{digest}" + + def emit(closure: Closure, *, source: str, function_name: str = "observables") -> str: """Python source for a module computing ``closure.requested`` from raw states. @@ -228,7 +266,10 @@ def emit(closure: Closure, *, source: str, function_name: str = "observables") - body = {name: _translate(closure.expressions[name], rename) for name in closure.order} lines = [ - f'"""Observables generated by qsp-codegen from {source}. Do not edit."""', + f'"""Observables generated by qsp-codegen from {source}. Do not edit.', + "", + f"emitter: {_emitter_stamp()}", + '"""', "", ] if any("_maximum(" in line for line in body.values()): @@ -239,6 +280,7 @@ def emit(closure: Closure, *, source: str, function_name: str = "observables") - " return (a + b + abs(a - b)) / 2", "", ] + si = closure.si lines += [ "STATES = (", *(f" {name!r}," for name in closure.states), @@ -248,25 +290,44 @@ def emit(closure: Closure, *, source: str, function_name: str = "observables") - *(f" {name!r}: {value!r}," for name, value in closure.constants.items()), "}", "", + "#: native * SI_FACTOR = SI. The rules are dimensionally consistent only in", + "#: SI, so inputs are converted in and results converted back out. Both the", + "#: arguments and the return value are therefore in the model's own units,", + "#: which is what the simulator writes and what a calibration target reads.", + "SI_FACTOR = {", + *(f" {name!r}: {si[name]!r}," + for name in (*closure.states, *closure.constants, *closure.requested) + if name in si), + "}", + "", "OBSERVABLES = (", *(f" {name!r}," for name in closure.requested), ")", "", "", f"def {function_name}(states, constants=CONSTANTS):", - f' """Every symbol in OBSERVABLES, from a mapping over STATES."""', + f' """Every symbol in OBSERVABLES, from a mapping over STATES.', + "", + " Native units in, native units out.", + ' """', ] for name in closure.states: - lines.append(f" {rename[name]} = states[{name!r}]") + factor = si.get(name, 1.0) + scale = "" if factor == 1.0 else f" * {factor!r}" + lines.append(f" {rename[name]} = states[{name!r}]{scale}") for name in closure.constants: - lines.append(f" {rename[name]} = constants[{name!r}]") + factor = si.get(name, 1.0) + scale = "" if factor == 1.0 else f" * {factor!r}" + lines.append(f" {rename[name]} = constants[{name!r}]{scale}") lines.append("") for name in closure.order: lines.append(f" {rename[name]} = {body[name]}") lines += [ "", " return {", - *(f" {name!r}: {rename[name]}," for name in closure.requested), + *(f" {name!r}: {rename[name]}" + + ("," if si.get(name, 1.0) == 1.0 else f" / {si[name]!r},") + for name in closure.requested), " }", "", ] diff --git a/tests/test_observables.py b/tests/test_observables.py index 30533bb..b2507fc 100644 --- a/tests/test_observables.py +++ b/tests/test_observables.py @@ -22,12 +22,21 @@ class _Model: """The four attributes resolve() reads. SBML parsing is covered elsewhere.""" - def __init__(self, rules=(), species=(), parameters=(), compartments=()): + def __init__(self, rules=(), species=(), parameters=(), compartments=(), + si=None): + # Each symbol is its own unit id, so a test can give one symbol a factor + # without inventing a unit system. self.assignment_rules = [{"variable_name": n, "expression": e} for n, e in rules] - self.species = [{"name": n} for n in species] - self.parameters = [{"name": n, "value": v} for n, v in parameters] - self.compartments = [{"name": n, "size": s} for n, s in compartments] + self.species = [{"name": n, "units": n} for n in species] + self.parameters = [{"name": n, "value": v, "units": n} + for n, v in parameters] + self.compartments = [{"name": n, "size": s, "units": n} + for n, s in compartments] + self._si = dict(si or {}) + + def get_si_factor(self, unit_id): + return self._si.get(unit_id, 1.0) MODEL = _Model( @@ -167,6 +176,44 @@ def test_a_translated_call_does_not_leak_into_the_closure(self): m = _Model(rules=[("x", "std::max(a, b)")], species=["a", "b"]) assert set(resolve(m, ["x"]).states) == {"a", "b"} + def test_a_rule_is_evaluated_in_si_not_in_declared_units(self): + """The PDAC V_T bug: mL plus micrometre^3 is only a sum after conversion.""" + m = _Model(rules=[("v", "vmin + a * vol")], species=["a"], + parameters=[("vmin", 1.0), ("vol", 2.0)], + si={"vol": 1e-12, "a": 1.0, "vmin": 1.0}) + # a=3 -> SI 1.0 + 3 * 2e-12; reading the raw numbers would give 7.0. + got = _run(emit(resolve(m, ["v"]), source="t.sbml"), {"a": 3.0})["v"] + assert got == pytest.approx(1.0 + 6e-12) + + def test_a_result_comes_back_in_the_units_it_was_asked_in(self): + """Native in, native out, so a calibration target reads what it expects. + + The rule variable is declared as a component, which is where its units + live and what SBML requires of it. + """ + m = _Model(rules=[("v", "a")], species=["a"], parameters=[("v", 0.0)], + si={"a": 1e-6, "v": 1e-6}) + assert _run(emit(resolve(m, ["v"]), source="t.sbml"), {"a": 5.0})["v"] \ + == pytest.approx(5.0) + + def test_a_requested_state_round_trips(self): + m = _Model(rules=[("v", "a")], species=["a"], parameters=[("v", 0.0)], + si={"a": 1e-6, "v": 1e-6}) + assert _run(emit(resolve(m, ["a", "v"]), source="t.sbml"), + {"a": 5.0})["a"] == pytest.approx(5.0) + + def test_an_input_without_units_is_refused(self): + """1.0 for an input silently rescales the expression; for an output it does not.""" + + class _NoUnits(_Model): + def get_si_factor(self, unit_id): + raise KeyError(unit_id) + + m = _NoUnits(rules=[("v", "a")], species=["a"]) + m.species = [{"name": "a"}] # no units key at all + with pytest.raises((UnresolvedSymbol, KeyError)): + resolve(m, ["v"]) + def test_colliding_names_raise(self): m = _Model(rules=[("s", "V_T.C + V_T_C")], species=["V_T.C", "V_T_C"]) with pytest.raises(NameCollision, match="sanitise"): @@ -219,3 +266,16 @@ def test_round_trip_from_real_sbml(tmp_path): out = _run(emit(closure, source="m.sbml"), {"V_T.CD8": 3.0, "V_T.CD8_exh": 4.0}) assert out["CD8_total"] == 3.0 + 4.0 * 10.0 + + +class TestStamp: + def test_the_header_carries_an_emitter_stamp(self): + m = _Model(rules=[("v", "a")], species=["a"], parameters=[("v", 0.0)]) + assert "emitter: " in emit(resolve(m, ["v"]), source="t.sbml") + + def test_the_stamp_moves_when_the_emitter_does(self): + """A version alone does not move on an editable install; the hash does.""" + from qsp_codegen.observables import _emitter_stamp + + assert _emitter_stamp().split("+")[1] != "" + assert len(_emitter_stamp().split("+")[1]) == 12