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