Evaluate observable rules in SI, not in declared units

src/qsp_codegen/observables.py

@@ -108,27 +108,39 @@ 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, ...]:
         """Requested symbols that are rules rather than raw states."""
         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),
         "    }",
         "",
     ]

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