diff --git a/CHANGELOG.md b/CHANGELOG.md index 931bd66..ace26cf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,34 @@ and this project adheres to [Semantic Versioning](https://semver.org/). ### Added +- `qsp-codegen verify` subcommand: end-to-end self-test that generates C++, + scaffolds and builds a minimal `qsp_sim`, and compares its trajectories + against a MATLAB SimBiology reference (PASS/FAIL). The CLI now uses + `generate`/`verify` subcommands; the legacy `qsp-codegen --sbml … --out-dir …` + form still routes to `generate`. +- `generate` now also emits a complete, ready-to-run `param_all.xml` (ICs + + parameters from the SBML), so the generated `qsp_sim` runs with no manual + `` wrapping or `qsp-refresh-param-xml` merge step. +- SBML `` (user lambda) inlining at call sites. +- Codegen-time validation of generated C++: catches `AUX_VAR_*` + use-before-definition (a dependency-ordering bug class) with a located, + human-readable error instead of an opaque downstream C++ compiler error. +- Tests for the MathML→C++ converter, generated-code validation, + `functionDefinition` inlining, and SBML event parsing. + +### Fixed + +- MathML→C++ converter now handles SimBiology's `nthroot` (and a + broader set of ``-exported operators: `log2`/`log10`, `ceiling`, and the + hyperbolic / inverse-trig families). Previously emitted a `/* unknown op */` + comment that crashed the Jacobian `sympify`. Unrecognized operators now + raise a clear, named error. +- Dependency-ordering bug for concentration assignment-rule species in a + dynamic-volume compartment (e.g. a drug concentration in a growing tumor + compartment): the compartment-volume temporary is now seeded into the + init / `update_y_other` emission closures, fixing an undefined `AUX_VAR_*` + reference in the generated C++. + - `qsp_sim` runtime flags `--time-unit days|seconds` and `--time-factor ` to override the compile-time time-scaling default per invocation. Lets a single binary integrate both unit-annotated SBML (default, runs in SI diff --git a/README.md b/README.md index e7e6eba..2cc5571 100644 --- a/README.md +++ b/README.md @@ -7,36 +7,87 @@ consumed by a CVODEBase-backed QSP simulator: - `QSP_enum.h`, `ODE_system.h`, `ODE_system.cpp` - `QSPParam.h`, `QSPParam.cpp` -- `qsp_params_xml_snippet.xml` (merged into consumer `param_all.xml` via - `qsp-refresh-param-xml`) +- `param_all.xml` — a complete, ready-to-run parameter file (initial + conditions + parameters taken straight from the SBML). Run the generated + `qsp_sim` against it directly — no hand-editing, no merge step. +- `qsp_params_xml_snippet.xml` — the bare `` block, for consumers that + merge it into a maintained `param_all.xml` via `qsp-refresh-param-xml`. + +## Prerequisites + +- **Python ≥ 3.10** with `sympy` (for the analytical Jacobian) and `numpy`. +- **A C++17 compiler + CMake ≥ 3.18** to build the generated simulator. The + first configure fetches and compiles SUNDIALS and yaml-cpp via CMake + `FetchContent` (one-time, a few minutes). +- **MATLAB + SimBiology** only if you want to (a) export SBML from a + SimBiology model or (b) run `qsp-codegen verify` against a SimBiology + reference. Pure SBML-in → C++-out needs neither. ## Install ```bash -pip install ~/Projects/qsp-codegen +pip install ~/Projects/qsp-codegen # or: uv pip install -e ~/Projects/qsp-codegen ``` -## Usage +## Quick start (SBML in → trajectories out) ```bash -qsp-codegen --sbml path/to/PDAC_model.sbml --out-dir path/to/ode/ +# 1. Generate C++ + a ready-to-run param file from a SimBiology SBML export. +qsp-codegen generate --sbml MyModel.sbml --out-dir build/qsp/ode -qsp-refresh-param-xml \ - --snippet path/to/ode/qsp_params_xml_snippet.xml \ - --xml path/to/param_all.xml \ - --xml path/to/param_all_test.xml +# 2. Build a minimal qsp_sim against the generated ODE (driver + ODE, default +# no-op init hook). See "Bundled C++ driver" below for the CMakeLists; or +# let `qsp-codegen verify` scaffold and build it for you. +cmake -S build/sim -B build/sim/build -DCMAKE_BUILD_TYPE=Release +cmake --build build/sim/build --target qsp_sim -j4 + +# 3. Run it directly against the emitted param_all.xml — no hand-editing. +build/sim/build/qsp_sim build/qsp/ode/param_all.xml out.csv 365 4.0 ``` -Run whenever the QSP model structure changes (species/parameters/reactions). -Not needed for parameter-value tweaks. +The legacy form `qsp-codegen --sbml … --out-dir …` (no `generate` word) still +works. Re-run codegen whenever the model *structure* changes +(species/parameters/reactions); not needed for parameter-value tweaks. + +## Self-test: does the C++ match SimBiology? + +`qsp-codegen verify` runs the whole pipeline end-to-end — codegen, scaffold + +build a minimal `qsp_sim`, then compare its trajectories against a MATLAB +SimBiology reference over the same window — and reports PASS/FAIL: -## Scope +```bash +qsp-codegen verify \ + --sbml MyModel.sbml \ + --matlab-dir /path/to/model/repo \ # has startup.m + the model script + --matlab-script build_my_model \ # builds `model`; must NOT `clear` + --stop-time 365 +``` -- Parses SBML Level 2 v4 (SimBiology export dialect). -- Derives analytical Jacobian via sympy + CSE when sympy is available; - falls back to numerical Jacobian otherwise. -- Consumer-side invariants (sync checks, ABM-specific param codegen) are - *not* in scope — they live in the consumer repo. +## Supported SBML features + +Parses **SBML Level 2 v4** (the SimBiology export dialect): + +- **Rate-law math**: the full SBML L2 operator set, including SimBiology's + ``-exported named operators (`max`, `min`, `nthroot`, …), `` with + ``, and the trig / hyperbolic / log families. +- **User functions**: SBML `` lambdas are inlined at their + call sites. +- **Rules**: assignment (`repeatedAssignment`) and initial assignments, + including dynamic compartment volumes (e.g. a growing tumor compartment). +- **Events**: single-comparison triggers (`lt`/`leq`/`gt`/`geq`) with event + assignments, mapped to CVODE root functions. *Not yet supported* (fail + loudly with a clear message): event ``s and compound (`and`/`or`) + triggers. +- **Jacobian**: analytical via sympy + CSE when sympy is present; numerical + fallback otherwise. +- **Units**: converts `` to SI for integration (see + *Time units* below). + +Anything outside this surface raises a **located, human-readable error** at +codegen time (naming the operator/function/event) rather than emitting broken +C++ — and the generator self-checks its output for undefined temporaries +before writing. Consumer-side invariants (sync checks, ABM-specific param +codegen) are out of scope — they live in the consumer repo. ## Bundled C++ driver (`qsp_sim_core`) diff --git a/cpp/include/qsp_sim_core/model_hooks.h b/cpp/include/qsp_sim_core/model_hooks.h index e217cbb..bb535bf 100644 --- a/cpp/include/qsp_sim_core/model_hooks.h +++ b/cpp/include/qsp_sim_core/model_hooks.h @@ -21,7 +21,7 @@ #include #include -namespace CancerVCT { +namespace qsp_sim_core { class ODE_system; @@ -73,6 +73,6 @@ struct EvolveResult { // success=true, t_diagnosis_days=0 — i.e. no pre-scenario evolve. EvolveResult evolve_to_diagnosis(ODE_system& ode, const EvolveOpts& opts); -} // namespace CancerVCT +} // namespace qsp_sim_core #endif diff --git a/cpp/include/qsp_sim_core/trajectory_writer.h b/cpp/include/qsp_sim_core/trajectory_writer.h index f140664..e1bcea6 100644 --- a/cpp/include/qsp_sim_core/trajectory_writer.h +++ b/cpp/include/qsp_sim_core/trajectory_writer.h @@ -44,7 +44,7 @@ #include #include -namespace CancerVCT { +namespace qsp_sim_core { class TrajectoryWriter { public: @@ -128,6 +128,6 @@ class TrajectoryWriter { uint64_t n_times_ = 0; }; -} // namespace CancerVCT +} // namespace qsp_sim_core #endif // QSP_SIM_CORE_TRAJECTORY_WRITER_H diff --git a/cpp/src/default_hooks.cpp b/cpp/src/default_hooks.cpp index aa78652..b048f76 100644 --- a/cpp/src/default_hooks.cpp +++ b/cpp/src/default_hooks.cpp @@ -10,7 +10,7 @@ #include "qsp_sim_core/model_hooks.h" -namespace CancerVCT { +namespace qsp_sim_core { EvolveResult evolve_to_diagnosis(ODE_system& /*ode*/, const EvolveOpts& /*opts*/) { EvolveResult r; @@ -20,4 +20,4 @@ EvolveResult evolve_to_diagnosis(ODE_system& /*ode*/, const EvolveOpts& /*opts*/ return r; } -} // namespace CancerVCT +} // namespace qsp_sim_core diff --git a/cpp/src/qsp_sim_main.cpp b/cpp/src/qsp_sim_main.cpp index 2be08f2..dc03332 100644 --- a/cpp/src/qsp_sim_main.cpp +++ b/cpp/src/qsp_sim_main.cpp @@ -74,7 +74,7 @@ #include "qsp_sim_core/model_hooks.h" #include "qsp_sim_core/trajectory_writer.h" -using namespace CancerVCT; +using namespace qsp_sim_core; namespace { diff --git a/src/qsp_codegen/codegen.py b/src/qsp_codegen/codegen.py index dbd1232..9dcef2f 100644 --- a/src/qsp_codegen/codegen.py +++ b/src/qsp_codegen/codegen.py @@ -51,11 +51,14 @@ def __init__(self, sbml_path: str): self.initial_assignments: List[dict] = [] self.events: List[dict] = [] self.unit_defs: Dict[str, float] = {} # unit_id → SI factor + # SBML id → {"bvars": [...], "body": } + self.function_defs: Dict[str, dict] = {} self._parse_units() self._parse_compartments() self._parse_parameters() # before species (species may ref param units) self._parse_species() + self._parse_function_definitions() # before reactions/rules (they may call them) self._parse_reactions() self._parse_rules() self._parse_initial_assignments() @@ -317,6 +320,63 @@ def _parse_initial_assignments(self): # --- MathML → Infix Converter ---------------------------------------- + def _parse_function_definitions(self): + """Parse lambdas so calls can be inlined. + + SBML stores a user function as ``x ... + ``. We keep the bound-variable names and the body + MathML element; :meth:`_inline_function_def` substitutes call-site + arguments for the bvars and converts the body in place. + """ + for fd in self.model.findall(f".//{SBML_NS}functionDefinition"): + fid = fd.get("id") + math = fd.find(f"{MATH_NS}math") + if not fid or math is None: + continue + lam = math.find(f"{MATH_NS}lambda") + if lam is None: + continue + bvars, body = [], None + for child in lam: + if child.tag.replace(MATH_NS, "") == "bvar": + ci = child.find(f"{MATH_NS}ci") + if ci is not None and ci.text: + bvars.append(ci.text.strip()) + else: + body = child # last non-bvar child is the function body + if body is not None: + self.function_defs[fid] = {"bvars": bvars, "body": body} + + def _inline_function_def(self, name: str, arg_nodes) -> str: + """Inline a call to an SBML as infix C++. + + Converts each actual argument, then converts the lambda body with the + bound variables temporarily resolving to those argument expressions. + Reentrant: nested/self-referential calls save and restore the shadowed + id_to_name entries. + """ + fd = self.function_defs[name] + bvars, body = fd["bvars"], fd["body"] + if len(arg_nodes) != len(bvars): + raise ValueError( + f"SBML functionDefinition '{name}' expects {len(bvars)} " + f"argument(s), called with {len(arg_nodes)}." + ) + arg_infix = [f"({self._mathml_to_infix(a)})" for a in arg_nodes] + missing = object() + saved = {} + for bv, ai in zip(bvars, arg_infix): + saved[bv] = self.id_to_name.get(bv, missing) + self.id_to_name[bv] = ai + try: + return f"({self._mathml_to_infix(body)})" + finally: + for bv, old in saved.items(): + if old is missing: + self.id_to_name.pop(bv, None) + else: + self.id_to_name[bv] = old + def _mathml_to_infix(self, node) -> str: """Convert a MathML or tree to infix C++ string. @@ -358,13 +418,37 @@ def _mathml_to_infix(self, node) -> str: op_node = children[0] op_tag = op_node.tag.replace(MATH_NS, "") - # SimBiology exports max/min as max instead of + # SimBiology exports several operators as name rather + # than a native MathML tag (e.g. max, nthroot). + # Treat any naming a known math function as that operator and + # dispatch it like a native tag. Keep this set in sync with the + # functions _apply_op (and the branch) can emit. if op_tag == "ci": func_name = (op_node.text or "").strip() - _CI_FUNCTIONS = {"max", "min", "abs", "floor", "ceil", - "exp", "log", "ln", "sqrt", "power"} + # User-defined SBML lambda: inline it. + if func_name in self.function_defs: + return self._inline_function_def(func_name, children[1:]) + _CI_FUNCTIONS = { + "max", "min", "abs", "floor", "ceil", "ceiling", + "exp", "ln", "log", "log2", "log10", "sqrt", + "power", "nthroot", + "sin", "cos", "tan", "sinh", "cosh", "tanh", + "asin", "acos", "atan", + } if func_name in _CI_FUNCTIONS: op_tag = func_name + else: + # A operator that is neither a known math function nor + # a parsed . Fail loudly with context + # rather than emitting a comment that later crashes the + # Jacobian sympify with no breadcrumb. + raise ValueError( + f"Unsupported function '{func_name}' applied in a " + f"MathML expression: not a known SBML L2 math operator " + f"(see _mathml_to_infix) and not a " + f"in this model. If it is a standard function, add it " + f"to _CI_FUNCTIONS and _apply_op." + ) # Handle specially: extract child if op_tag == "root": @@ -448,6 +532,9 @@ def _apply_op(self, op: str, args: List[str]) -> str: "log10": "std::log10", "exp": "std::exp", "sqrt": "std::sqrt", "abs": "std::abs", "floor": "std::floor", "ceiling": "std::ceil", "sin": "std::sin", "cos": "std::cos", "tan": "std::tan", + "sinh": "std::sinh", "cosh": "std::cosh", "tanh": "std::tanh", + "asin": "std::asin", "acos": "std::acos", "atan": "std::atan", + "ceil": "std::ceil", } if op in func_map: return f"{func_map[op]}({', '.join(args)})" @@ -475,7 +562,18 @@ def _apply_op(self, op: str, args: List[str]) -> str: return f"std::pow({args[1]}, 1.0 / {args[0]})" return f"std::sqrt({args[0]})" - return f"/* unknown op: {op} */({', '.join(args)})" + if op == "nthroot": + # SimBiology's nthroot: nthroot(x, n) == x^(1/n). + # Args are (radicand, degree) — the reverse of 's order. + if len(args) == 2: + return f"std::pow({args[0]}, 1.0 / ({args[1]}))" + return f"std::sqrt({args[0]})" + + raise ValueError( + f"Unsupported MathML operator '{op}'. Add a handler in _apply_op " + f"(and list it in _mathml_to_infix's set if SimBiology " + f"exports it as {op})." + ) def _convert_piecewise(self, node) -> str: """Convert MathML to C++ ternary.""" @@ -1230,7 +1328,7 @@ def gen_enum_h(sbml: SBMLModel) -> str: "", "#include ", "", - "namespace CancerVCT{", + "namespace qsp_sim_core{", "", "// QSP Species Enum (ODE state vector indices)", "enum QSPSpeciesEnum", @@ -1273,7 +1371,7 @@ def gen_enum_h(sbml: SBMLModel) -> str: lines.append(f"QSP_FILE_PARAM_COUNT") lines.append("};") lines.append("") - lines.append("} // namespace CancerVCT") + lines.append("} // namespace qsp_sim_core") return "\n".join(lines) + "\n" @@ -1299,8 +1397,8 @@ def gen_ode_h(jac_nnz: int = 0) -> str: ' CVLsJacFn getJacobianFn() const override { return &ODE_system::jac; }\n' ) - return ('''#ifndef __CancerVCT_ODE__ -#define __CancerVCT_ODE__ + return ('''#ifndef __qsp_sim_core_ODE__ +#define __qsp_sim_core_ODE__ // Auto-generated by qsp_codegen.py from SBML — do not edit manually @@ -1310,7 +1408,7 @@ def gen_ode_h(jac_nnz: int = 0) -> str: #include #include -namespace CancerVCT{ +namespace qsp_sim_core{ class ODE_system : public CVODEBase @@ -1395,7 +1493,7 @@ class ODE_system : _class_parameter[i] = v; } -} // namespace CancerVCT +} // namespace qsp_sim_core #endif ''') @@ -1411,7 +1509,7 @@ def gen_ode_cpp(sbml: SBMLModel, mapping: dict) -> str: lines.append('#define PARAM(x) _class_parameter[x]') lines.append('#define PFILE(x) param.getVal(x)') lines.append('') - lines.append('namespace CancerVCT{') + lines.append('namespace qsp_sim_core{') lines.append('#define QSP_W ODE_system::_QSP_weight') lines.append('') lines.append('bool ODE_system::use_steady_state = false;') @@ -1978,6 +2076,18 @@ def collect_ar_deps(vn): collect_ar_deps(other_vn) for r, sp in ar_species: collect_ar_deps(r["variable_name"]) + # A concentration rule-species writes back as + # `_species_var[SP] = AUX_VAR * AUX_VAR_{comp}` when its + # compartment volume is itself an assignment rule (dynamic + # volume, e.g. a growing tumor compartment V_T). That multiplier + # is introduced by the writeback, not by the rule expression, so + # collect_ar_deps over the expression alone misses it and + # AUX_VAR_{comp} is emitted undefined. Seed it explicitly, the + # same way the output generator seeds dynamic compartment volumes. + if not sp.get("has_only_substance_units", False): + comp_name = sp["compartment"] + if comp_name in rule_names_ar: + collect_ar_deps(comp_name) needed_list = [r for r in sbml.assignment_rules if r["variable_name"] in needed_rules] ordered_needed = order_rules(needed_list, needed_rules, sbml=sbml) @@ -2119,7 +2229,7 @@ def collect_ar_deps(vn): # via CVodeSetJacFn when built with KLU support; dense builds ignore it. if gen_ode_cpp._jacobian_info is not None: lines.append(gen_jacobian_cpp(sbml, gen_ode_cpp._jacobian_info)) - lines.append('} // namespace CancerVCT') + lines.append('} // namespace qsp_sim_core') return "\n".join(lines) + "\n" @@ -2129,14 +2239,14 @@ def collect_ar_deps(vn): def gen_qsp_param_h() -> str: - return '''#ifndef __CancerVCT_QSPParam__ -#define __CancerVCT_QSPParam__ + return '''#ifndef __qsp_sim_core_QSPParam__ +#define __qsp_sim_core_QSPParam__ // Auto-generated by qsp_codegen.py from SBML — do not edit manually #include "qsp_sim_core/ParamBase.h" -namespace CancerVCT{ +namespace qsp_sim_core{ class QSPParam : public SP_QSP_IO::ParamBase { @@ -2157,7 +2267,7 @@ class QSPParam : public SP_QSP_IO::ParamBase static const char* _xml_paths[]; }; -} // namespace CancerVCT +} // namespace qsp_sim_core #endif ''' @@ -2172,7 +2282,7 @@ def gen_qsp_param_cpp(sbml: SBMLModel) -> str: lines.append('#include ') lines.append('#include ') lines.append('') - lines.append('namespace CancerVCT{') + lines.append('namespace qsp_sim_core{') lines.append('') # XML paths: Compartment ICs, Species ICs, Model Parameters @@ -2214,7 +2324,7 @@ def gen_qsp_param_cpp(sbml: SBMLModel) -> str: lines.append(' std::cout << _xml_paths[i] << " = " << _param[i] << std::endl;') lines.append('}') lines.append('') - lines.append('} // namespace CancerVCT') + lines.append('} // namespace qsp_sim_core') return "\n".join(lines) + "\n" @@ -2265,6 +2375,83 @@ def gen_xml_snippet(sbml: SBMLModel) -> str: # Main # ========================================================================= +_AUX_DECL_RE = re.compile(r"\brealtype\s+(AUX_VAR_\w+)\s*=") +_AUX_USE_RE = re.compile(r"\bAUX_VAR_\w+\b") + + +def validate_generated_cpp(cpp: str, filename: str = "ODE_system.cpp") -> None: + """Catch use-before-definition of ``AUX_VAR_*`` temporaries at codegen time. + + The generator emits assignment-rule / compartment-volume temporaries as + function-local ``realtype AUX_VAR_x = ...;``. A dependency-ordering bug can + emit a *use* of ``AUX_VAR_x`` in a function/block where it was never + declared (e.g. a concentration rule-species whose dynamic compartment + volume multiplier was not seeded into the block's dependency closure). + That otherwise surfaces only as an opaque C++ compiler error pointing at + machine-generated line numbers; here we fail at generation time with a + located, human-readable message. + + Scoping model: a stack of brace scopes (the file's ``namespace`` is the + outermost). An ``AUX_VAR`` is in scope if declared in the current or any + enclosing scope, and — within a scope — a declaration must precede its + uses (C++ rule). Each function body is its own scope, so the same temp may + legitimately be redeclared across functions. + """ + scopes: List[set] = [set()] + offenders: List[Tuple[int, str]] = [] + for lineno, line in enumerate(cpp.splitlines(), 1): + # Process the line as an ordered sequence of brace tokens and + # ';'-separated statements (one-statement-per-line is the usual + # generated form, but handling multiple keeps this robust). Within a + # statement, uses are checked against the current scope chain *before* + # the statement's own declaration is registered, so same-statement + # `realtype AUX_VAR_x = ` and same-line `... ; ... AUX_VAR_x` + # both resolve correctly. + for seg in re.split(r"([{};])", line): + if seg == "{": + scopes.append(set()) + elif seg == "}": + if len(scopes) > 1: + scopes.pop() + elif seg in (";", ""): + continue + else: + decl = _AUX_DECL_RE.search(seg) + # On a declaration the LHS name is being declared, not used — + # only scan the RHS (after '=') for uses. + check_region = seg[decl.end():] if decl else seg + in_scope = set().union(*scopes) + for m in _AUX_USE_RE.finditer(check_region): + if m.group(0) not in in_scope: + offenders.append((lineno, m.group(0))) + if decl is not None: + scopes[-1].add(decl.group(1)) + + if offenders: + first = offenders[0] + raise ValueError( + f"Generated {filename} references {len(offenders)} undefined " + f"AUX_VAR temporary(ies) — a codegen dependency-ordering bug. " + f"First: '{first[1]}' used before declaration at line {first[0]}. " + f"This means an assignment-rule/compartment-volume dependency was " + f"not seeded into a block's emission closure." + ) + + +def wrap_param_xml(snippet: str) -> str: + """Wrap a codegen ```` snippet into a complete, runnable param_all.xml. + + The snippet emitted alongside the C++ is a bare ``...`` block + (for consumers that merge it into a maintained file via + qsp-refresh-param-xml). Wrapping it in ```` yields a file the + generated ``qsp_sim`` can read directly — no hand-editing, no merge step. + """ + s = snippet.strip() + if not (s.startswith("") and s.endswith("")): + raise ValueError("snippet is not a bare block") + return '\n\n' + s + "\n\n" + + def generate(sbml_path: str, out_dir: str) -> Dict[str, str]: """Run codegen for the given SBML, writing all outputs under ``out_dir``. @@ -2300,15 +2487,25 @@ def generate(sbml_path: str, out_dir: str) -> Dict[str, str]: jac_nnz = 0 print("\nGenerating C++ files...") + xml_snippet = gen_xml_snippet(sbml) files = { "QSP_enum.h": gen_enum_h(sbml), "ODE_system.h": gen_ode_h(jac_nnz=jac_nnz), "ODE_system.cpp": gen_ode_cpp(sbml, mapping), "QSPParam.h": gen_qsp_param_h(), "QSPParam.cpp": gen_qsp_param_cpp(sbml), - "qsp_params_xml_snippet.xml": gen_xml_snippet(sbml), + "qsp_params_xml_snippet.xml": xml_snippet, + # Complete, ready-to-run parameter file (ICs + params straight from the + # SBML). Lets users run the generated qsp_sim immediately — no manual + # wrapping, no qsp-refresh-param-xml merge. The bare snippet + # above is still emitted for the merge-into-maintained-file workflow. + "param_all.xml": wrap_param_xml(xml_snippet), } + # Fail fast on generated-code dependency-ordering bugs (clear message at + # codegen time, vs. an opaque C++ compiler error on generated lines). + validate_generated_cpp(files["ODE_system.cpp"]) + os.makedirs(out_dir, exist_ok=True) for fname, content in files.items(): path = os.path.join(out_dir, fname) @@ -2320,23 +2517,44 @@ def generate(sbml_path: str, out_dir: str) -> Dict[str, str]: return files +def _handle_generate(args) -> int: + generate(args.sbml, args.out_dir) + return 0 + + def main(argv: Optional[List[str]] = None) -> int: + if argv is None: + argv = sys.argv[1:] + # Back-compat: the original CLI was `qsp-codegen --sbml X --out-dir Y` with + # no subcommand. If the first token is an option (not -h/--help), assume the + # `generate` subcommand so existing callers (Makefiles, scripts) keep working. + if argv and argv[0].startswith("-") and argv[0] not in ("-h", "--help"): + argv = ["generate", *argv] + parser = argparse.ArgumentParser( description="SBML → C++ CVODE ODE code generator.", ) - parser.add_argument( - "--sbml", - required=True, + sub = parser.add_subparsers(dest="command") + + gen = sub.add_parser("generate", help="Generate C++ ODE sources from SBML.") + gen.add_argument( + "--sbml", required=True, help="Path to SBML Level 2 v4 file (e.g. PDAC_model.sbml).", ) - parser.add_argument( - "--out-dir", - required=True, + gen.add_argument( + "--out-dir", required=True, help="Directory where generated C++ sources are written.", ) + gen.set_defaults(_handler=_handle_generate) + + from .verify import add_subparser as _add_verify + _add_verify(sub) + args = parser.parse_args(argv) - generate(args.sbml, args.out_dir) - return 0 + if not getattr(args, "command", None): + parser.print_help() + return 1 + return args._handler(args) if __name__ == "__main__": diff --git a/src/qsp_codegen/verify.py b/src/qsp_codegen/verify.py new file mode 100644 index 0000000..14c049b --- /dev/null +++ b/src/qsp_codegen/verify.py @@ -0,0 +1,176 @@ +"""``qsp-codegen verify`` — end-to-end self-test for a generated model. + +Answers the first question any new user has: *does the generated C++ actually +reproduce my SimBiology model?* In one command it: + + 1. generates C++ from the SBML (reusing :func:`qsp_codegen.codegen.generate`), + 2. scaffolds a minimal ``qsp_sim`` CMake project (driver + generated ODE + + the library's default no-op model-init hook) and builds it, + 3. assembles a ready-to-run ``param_all.xml`` from the codegen snippet, + 4. runs the compiled binary and a MATLAB SimBiology reference over the same + window, and compares every matched species. + +The build step needs a C++ toolchain + CMake (the first configure fetches and +compiles SUNDIALS/yaml-cpp); the reference step needs MATLAB on the path (or +``--matlab``). The MATLAB model script must build a ``model`` variable and must +NOT start with ``clear`` (it is ``run`` inside the parity harness's workspace). +""" +from __future__ import annotations + +import argparse +import subprocess +import sys +from pathlib import Path +from typing import Optional + +from .codegen import generate + +# Minimal, model-agnostic consumer CMake project. Mirrors the thin +# pdac-build/cpp/sim consumer but with NO model-specific init hook: the +# library's default no-op evolve_to_diagnosis (default_hooks.cpp) is linked, +# which is exactly right for a default-IC parity. ``${ODE_DIR}`` is filled in. +_CMAKE_TEMPLATE = """cmake_minimum_required(VERSION 3.18) +project(qsp_sim CXX) + +if(NOT DEFINED QSP_SIM_CORE_PREFIX) + find_package(Python3 REQUIRED COMPONENTS Interpreter) + execute_process( + COMMAND "${Python3_EXECUTABLE}" -m qsp_codegen.cmake --prefix + OUTPUT_VARIABLE QSP_SIM_CORE_PREFIX + OUTPUT_STRIP_TRAILING_WHITESPACE + RESULT_VARIABLE _qsp_rc) + if(NOT _qsp_rc EQUAL 0) + message(FATAL_ERROR "Could not locate qsp_sim_core (qsp-codegen wheel).") + endif() +endif() +list(APPEND CMAKE_PREFIX_PATH "${QSP_SIM_CORE_PREFIX}") +find_package(qsp_sim_core CONFIG REQUIRED) + +add_executable(qsp_sim + ${QSP_SIM_CORE_DRIVER_SOURCE} + @ODE_DIR@/ODE_system.cpp + @ODE_DIR@/QSPParam.cpp +) +target_include_directories(qsp_sim PRIVATE @ODE_DIR@) +target_link_libraries(qsp_sim PRIVATE qsp_sim_core::qsp_sim_core) +""" + + +def assemble_param_xml(snippet_path: Path, out_path: Path) -> Path: + """Wrap a codegen ```` snippet into a complete ``param_all.xml``. + + ``generate()`` now emits ``param_all.xml`` directly; this remains for + callers that only have the snippet on hand. + """ + from .codegen import wrap_param_xml + + out_path.write_text(wrap_param_xml(snippet_path.read_text())) + return out_path + + +def _build_qsp_sim(ode_dir: Path, work_dir: Path, python_exe: str) -> Path: + """Scaffold + build a minimal qsp_sim against the generated ODE.""" + sim_dir = work_dir / "sim" + sim_dir.mkdir(parents=True, exist_ok=True) + cmake_txt = _CMAKE_TEMPLATE.replace("@ODE_DIR@", str(ode_dir.resolve())) + (sim_dir / "CMakeLists.txt").write_text(cmake_txt) + build_dir = sim_dir / "build" + + print(f" configuring + building qsp_sim in {build_dir} ...") + subprocess.run( + ["cmake", "-S", str(sim_dir), "-B", str(build_dir), + "-DCMAKE_BUILD_TYPE=Release", f"-DPython3_EXECUTABLE={python_exe}"], + check=True, capture_output=True, text=True, + ) + subprocess.run( + ["cmake", "--build", str(build_dir), "--target", "qsp_sim", "-j", "4"], + check=True, capture_output=True, text=True, + ) + binary = build_dir / "qsp_sim" + if not binary.exists(): + raise RuntimeError(f"build reported success but no binary at {binary}") + return binary + + +def run_verify( + sbml: Path, + matlab_dir: Path, + matlab_script: str, + work_dir: Path, + stop_time: float = 365.0, + rtol: float = 0.05, + atol: float = 1e-6, + matlab_binary: str = "matlab", + python_exe: Optional[str] = None, +) -> bool: + """Run the full codegen → build → parity self-test. Returns pass/fail.""" + # Imported lazily so `generate`/validation work without numpy/matlab present. + from .parity import compare, run_cpp_trajectories, run_matlab_trajectories + + python_exe = python_exe or sys.executable + work_dir.mkdir(parents=True, exist_ok=True) + ode_dir = work_dir / "qsp" / "ode" + + print(f"[1/4] codegen {sbml} -> {ode_dir}") + generate(str(sbml), str(ode_dir)) + + print("[2/4] build qsp_sim") + qsp_sim = _build_qsp_sim(ode_dir, work_dir, python_exe) + + param_xml = assemble_param_xml( + ode_dir / "qsp_params_xml_snippet.xml", work_dir / "param_all.xml" + ) + + cpp_csv = work_dir / "cpp.csv" + matlab_csv = work_dir / "matlab.csv" + print(f"[3/4] run C++ + MATLAB ({stop_time}d, default ICs, grid-pinned)") + run_cpp_trajectories( + qsp_sim=qsp_sim, param_xml=param_xml, out_csv=cpp_csv, + t_end_days=stop_time, min_cadence_hours=4.0, + ) + # Pin MATLAB to the C++ grid so the compare is row-aligned (no interp + # artifacts on stiff early transients). + import numpy as np + cpp_times = np.loadtxt(cpp_csv, delimiter=",", skiprows=1, usecols=0) + times_csv = work_dir / "cpp_times.csv" + np.savetxt(times_csv, cpp_times) + run_matlab_trajectories( + matlab_model_dir=matlab_dir, matlab_model_script=matlab_script, + sbml_path=sbml, param_xml=param_xml, out_csv=matlab_csv, + stop_time=stop_time, output_times_csv=times_csv, matlab_binary=matlab_binary, + ) + + print(f"[4/4] compare (rtol={rtol}, atol={atol})") + passed, report = compare(str(matlab_csv), str(cpp_csv), rtol=rtol, atol=atol) + print(report) + print("\nVERIFY:", "PASS ✅" if passed else "FAIL ❌") + return passed + + +def add_subparser(subparsers) -> None: + """Register the ``verify`` subcommand on a codegen argparse subparsers.""" + p = subparsers.add_parser( + "verify", + help="Codegen + build + C++↔MATLAB parity self-test for an SBML model.", + ) + p.add_argument("--sbml", required=True, type=Path, help="SBML model file.") + p.add_argument("--matlab-dir", required=True, type=Path, + help="Consumer repo root (has startup.m + the model script).") + p.add_argument("--matlab-script", required=True, + help="Bare script name that builds `model` (no `.m`, no `clear`).") + p.add_argument("--work-dir", type=Path, default=Path("qsp_verify_out"), + help="Scratch dir for generated code, build, and CSVs.") + p.add_argument("--stop-time", type=float, default=365.0, help="Sim days.") + p.add_argument("--rtol", type=float, default=0.05) + p.add_argument("--atol", type=float, default=1e-6) + p.add_argument("--matlab", default="matlab", help="MATLAB binary path.") + p.set_defaults(_handler=_handle) + + +def _handle(args) -> int: + ok = run_verify( + sbml=args.sbml, matlab_dir=args.matlab_dir, matlab_script=args.matlab_script, + work_dir=args.work_dir, stop_time=args.stop_time, rtol=args.rtol, + atol=args.atol, matlab_binary=args.matlab, + ) + return 0 if ok else 1 diff --git a/tests/test_events.py b/tests/test_events.py new file mode 100644 index 0000000..6d83871 --- /dev/null +++ b/tests/test_events.py @@ -0,0 +1,79 @@ +"""Tests for SBML parsing: supported single-comparison triggers + +assignments, and the loud failures for delays / compound triggers. + +End-to-end event behavior is additionally covered by the TNBC parity run +(3 events, including a `V_T.C < 0.5*cell -> V_T.C = 0.01*cell` reset), which +matches the SimBiology reference to within tolerance. +""" +import xml.etree.ElementTree as ET + +import pytest + +from qsp_codegen.codegen import SBMLModel, SBML_NS, MATH_NS + +SNS = SBML_NS.strip("{}") +MNS = MATH_NS.strip("{}") + + +def _parse_events(event_xml: str): + doc = ( + f'' + f'' + f"{event_xml}" + f"" + ) + root = ET.fromstring(doc) + m = object.__new__(SBMLModel) + m.model = root.find(f"{SBML_NS}model") + m.id_to_name = {} + m.function_defs = {} + m.events = [] + m._parse_events() + return m.events + + +def _trigger(op="lt", delay=""): + return ( + f'' + f'<{op}/>C0.5' + f"{delay}" + f'' + f'0.01' + f"" + ) + + +def test_single_comparison_event_parsed(): + events = _parse_events(_trigger("lt")) + assert len(events) == 1 + ev = events[0] + assert ev["trigger_op"] == "lt" + assert ev["trigger_left"] == "C" + assert ev["trigger_right"] == "0.5" + assert ev["assignments"][0]["variable_id"] == "C" + assert "0.01" in ev["assignments"][0]["expression"] + + +@pytest.mark.parametrize("op", ["lt", "leq", "gt", "geq"]) +def test_supported_comparison_ops(op): + assert _parse_events(_trigger(op))[0]["trigger_op"] == op + + +def test_delay_raises(): + delay = f'1' + with pytest.raises(NotImplementedError, match="delay"): + _parse_events(_trigger("lt", delay=delay)) + + +def test_compound_trigger_raises(): + ev = ( + f'' + f"C0.5" + f"C0.1" + f'' + f'' + f'0.01' + f"" + ) + with pytest.raises(NotImplementedError, match="unsupported trigger op"): + _parse_events(ev) diff --git a/tests/test_function_defs.py b/tests/test_function_defs.py new file mode 100644 index 0000000..162c51e --- /dev/null +++ b/tests/test_function_defs.py @@ -0,0 +1,55 @@ +"""Tests for SBML inlining in the MathML->C++ converter.""" +import xml.etree.ElementTree as ET + +import pytest + +from qsp_codegen.codegen import SBMLModel + +MATH = "http://www.w3.org/1998/Math/MathML" + + +def _model_with_hill(): + """Bare SBMLModel with a hand-built hill(x, k) = x / (x + k) lambda.""" + m = object.__new__(SBMLModel) + m.id_to_name = {} + body = ET.fromstring( + f'x' + f"xk" + ) + m.function_defs = {"hill": {"bvars": ["x", "k"], "body": body}} + return m + + +def _call(model, mathml): + return model._mathml_to_infix(ET.fromstring(mathml)) + + +def test_function_def_is_inlined(): + m = _model_with_hill() + out = _call( + m, + f'hill' + f'S2', + ) + # Arguments substituted in; bound variable names must not leak. + assert "S" in out and "2.0" in out and "/" in out + assert "x" not in out and "k" not in out + + +def test_function_def_arg_count_mismatch_raises(): + m = _model_with_hill() + with pytest.raises(ValueError, match="expects 2"): + _call(m, f'hillS') + + +def test_nested_function_call_inlines(): + # hill(hill(S, 1), 2) — exercises reentrant bvar save/restore. + m = _model_with_hill() + out = _call( + m, + f'hill' + f'hillS1' + f'2', + ) + assert "S" in out and "1.0" in out and "2.0" in out + assert "x" not in out and "k" not in out diff --git a/tests/test_mathml_ops.py b/tests/test_mathml_ops.py new file mode 100644 index 0000000..9197de9 --- /dev/null +++ b/tests/test_mathml_ops.py @@ -0,0 +1,63 @@ +"""Unit tests for the MathML -> C++ infix converter (SBMLModel._mathml_to_infix). + +Guards the operator coverage that SimBiology exercises via SBML export. The +motivating regression: SimBiology exports `nthroot` (and max/min) as +``name`` rather than a native MathML tag; an unsupported name used to +fall through to a ``/* unknown op */`` comment that only blew up much later in +the Jacobian ``sympify``. The converter now dispatches the full L2 math set and +fails loudly on anything it genuinely cannot translate. +""" +import xml.etree.ElementTree as ET + +import pytest + +from qsp_codegen.codegen import SBMLModel + +MATH = "http://www.w3.org/1998/Math/MathML" + + +def _convert(mathml: str) -> str: + # _mathml_to_infix only touches self.id_to_name; build a bare instance. + model = object.__new__(SBMLModel) + model.id_to_name = {} + model.function_defs = {} + return model._mathml_to_infix(ET.fromstring(mathml)) + + +def _apply(op_xml: str, *operand_xml: str) -> str: + operands = "".join(operand_xml) + return _convert(f'{op_xml}{operands}') + + +X = "x" +N3 = '3' +N2 = '2' + + +def test_nthroot_ci_becomes_pow_reciprocal(): + # nthroot(x, 3) == x^(1/3) — the TNBC vasculature K^(2/3) term. + out = _apply("nthroot", X, N3) + assert "std::pow(x, 1.0 / (3.0))" in out + + +def test_power_tag(): + assert _apply("", X, N2) == "std::pow(x, 2.0)" + + +def test_root_with_degree(): + out = _apply("", f"{N3}", X) + assert "std::pow(x, 1.0 / 3.0)" in out + + +def test_ci_exported_max(): + # SimBiology exports max as max. + assert "std::max" in _apply("max", X, '0') + + +def test_hyperbolic_ci(): + assert "std::tanh(x)" == _apply("tanh", X) + + +def test_unknown_ci_function_raises_with_context(): + with pytest.raises(ValueError, match="mysteryFunc"): + _apply("mysteryFunc", X) diff --git a/tests/test_smoke.py b/tests/test_smoke.py index 38df42e..92a6107 100644 --- a/tests/test_smoke.py +++ b/tests/test_smoke.py @@ -18,6 +18,7 @@ "QSPParam.h", "QSPParam.cpp", "qsp_params_xml_snippet.xml", + "param_all.xml", # complete, ready-to-run param file (wraps the snippet) } @@ -30,8 +31,25 @@ def test_generate_pdac(tmp_path): assert p.exists() and p.stat().st_size > 0 -def test_cli_requires_sbml_and_out(): +def test_cli_generate_requires_sbml_and_out(): from qsp_codegen.codegen import main + # The `generate` subcommand still enforces its required args. with pytest.raises(SystemExit): - main([]) + main(["generate"]) + + +def test_cli_no_subcommand_prints_help_and_returns_nonzero(): + from qsp_codegen.codegen import main + + assert main([]) == 1 + + +def test_cli_legacy_bare_options_still_route_to_generate(tmp_path): + # Back-compat: `qsp-codegen --sbml ... --out-dir ...` (no subcommand word) + # must still reach `generate`. With a non-existent SBML it should fail-fast, + # not silently no-op — assert it raises rather than returning 0. + from qsp_codegen.codegen import main + + with pytest.raises(Exception): + main(["--sbml", str(tmp_path / "nope.sbml"), "--out-dir", str(tmp_path)]) diff --git a/tests/test_validate_generated.py b/tests/test_validate_generated.py new file mode 100644 index 0000000..f0242f6 --- /dev/null +++ b/tests/test_validate_generated.py @@ -0,0 +1,60 @@ +"""Tests for validate_generated_cpp — the codegen-time guard against +use-before-definition of AUX_VAR_* temporaries in generated C++. + +Regression context: a concentration rule-species in a dynamic-volume +compartment emitted `AUX_VAR_V_T` undefined in the init/update_y_other blocks, +surfacing only as an opaque C++ compiler error on machine-generated lines. This +guard turns that class of dependency-ordering bug into a clear codegen-time error. +""" +import pytest + +from qsp_codegen.codegen import validate_generated_cpp as validate + + +def test_clean_passes(): + cpp = ( + "namespace X {\n" + "void f(){\n" + " realtype AUX_VAR_V_T = 1.0;\n" + " realtype c = AUX_VAR_V_T * 2.0;\n" + "}\n" + "}\n" + ) + validate(cpp) # no raise + + +def test_redeclaration_across_functions_ok(): + # Same temp legitimately redeclared in separate function scopes. + cpp = ( + "namespace X {\n" + "void f(){ realtype AUX_VAR_V_T = 1.0; realtype z = AUX_VAR_V_T; }\n" + "void g(){ realtype AUX_VAR_V_T = 2.0; realtype w = AUX_VAR_V_T; }\n" + "}\n" + ) + validate(cpp) # no raise + + +def test_use_before_decl_multiline_raises(): + cpp = ( + "namespace X {\n" + "void f(){\n" + " realtype y = AUX_VAR_V_T * 2.0;\n" # used here + " realtype AUX_VAR_V_T = 3.0;\n" # declared after + "}\n" + "}\n" + ) + with pytest.raises(ValueError, match="AUX_VAR_V_T"): + validate(cpp) + + +def test_use_in_function_without_decl_raises(): + # The exact shape of the real bug: declared in one function, used in + # another that never declares it. + cpp = ( + "namespace X {\n" + "void f(){ realtype AUX_VAR_V_T = 1.0; realtype a = AUX_VAR_V_T; }\n" + "void g(){ realtype b = AUX_VAR_V_T * 2.0; }\n" # no decl in g + "}\n" + ) + with pytest.raises(ValueError, match="line"): + validate(cpp)