Skills › Scientific Agent Skills › tellurium › kinetic_experiment.py
kinetic_experiment.py
209 lines · skills/tellurium/scripts/kinetic_experiment.py at 5f2bc03
1"""Simulate local SBML/Antimony models and verify replay of a generated COMBINE archive."""
2from __future__ import annotations
4import argparse
5import csv
6import hashlib
7import importlib.metadata
8import json
9import math
10from pathlib import Path
11import platform
12import re
15def model_document(path, model_format):
16 import libsbml
17 import tellurium as te
18 if model_format not in {"antimony", "sbml"}:
19 raise ValueError("model_format must be antimony or sbml")
20 source = Path(path).read_text(encoding="utf-8")
21 sbml = te.antimonyToSBML(source) if model_format == "antimony" else source
22 document = libsbml.readSBMLFromString(sbml)
23 if document.getModel() is None:
24 raise ValueError("Input has no SBML model")
25 document.checkConsistency()
26 issues = [{"id": document.getError(i).getErrorId(), "severity": document.getError(i).getSeverityAsString(),
27 "message": document.getError(i).getMessage().strip()} for i in range(document.getNumErrors())]
28 if any(issue["severity"] in {"Error", "Fatal"} for issue in issues):
29 raise ValueError("SBML validation errors: " + json.dumps(issues))
30 return document, issues
33def validate_experiment(config, model):
34 allowed = {"end_time", "points", "relative_tolerance", "absolute_tolerance", "species", "scenarios"}
35 if not isinstance(config, dict) or set(config) != allowed:
36 raise ValueError(f"Experiment requires exactly {sorted(allowed)}")
37 if any(rule.isRate() for rule in model.getListOfRules()):
38 raise ValueError("Rate-rule models require a separate validated tolerance workflow; "
39 "RoadRunner 2.10.0 scalar tolerance and state-vector orders can differ")
40 for name in ("end_time", "relative_tolerance", "absolute_tolerance"):
41 value = config[name]
42 if isinstance(value, bool) or not isinstance(value, (int, float)) or not math.isfinite(value) or value <= 0:
43 raise ValueError(f"{name} must be positive and finite")
44 if isinstance(config["points"], bool) or not isinstance(config["points"], int) or config["points"] < 2:
45 raise ValueError("points must be an integer >=2 (including both endpoints)")
46 species = config["species"]
47 if (not isinstance(species, list) or not species
48 or any(not isinstance(s, str) for s in species)
49 or len(set(species)) != len(species) or any(model.getSpecies(s) is None for s in species)):
50 raise ValueError("species must be a nonempty list of unique SBML species IDs")
51 # SED-ML species targets below use SBML concentration semantics.
52 if any(model.getSpecies(s).getHasOnlySubstanceUnits() for s in species):
53 raise ValueError("This archive workflow selects concentrations; amount-only species need an explicit amount workflow")
54 if any(model.getCompartment(model.getSpecies(s).getCompartment()).getSpatialDimensions() == 0 for s in species):
55 raise ValueError("Concentration selections require compartments with nonzero spatial dimensions")
56 scenarios = config["scenarios"]
57 if not isinstance(scenarios, dict) or scenarios.get("baseline") != {}:
58 raise ValueError("scenarios must include an unchanged baseline mapping")
59 for name, changes in scenarios.items():
60 if not isinstance(name, str) or not re.fullmatch(r"[A-Za-z_][A-Za-z0-9_]*", name) or not isinstance(changes, dict):
61 raise ValueError("Scenario names must be SBML-style identifiers and changes must be a mapping")
62 for symbol, value in changes.items():
63 if not isinstance(symbol, str):
64 raise ValueError("Parameter IDs must be strings")
65 parameter = model.getParameter(symbol)
66 if parameter is None or not parameter.getConstant() or model.getInitialAssignmentBySymbol(symbol) is not None:
67 raise ValueError(f"{symbol}: perturbations require a constant global parameter without an initial assignment")
68 if isinstance(value, bool) or not isinstance(value, (int, float)) or not math.isfinite(value):
69 raise ValueError(f"{symbol}: perturbation must be a finite number")
70 return config
73def unique_json_object(pairs):
74 """Do not silently discard a duplicated scientific setting or condition."""
75 result = {}
76 for key, value in pairs:
77 if key in result:
78 raise ValueError(f"Duplicate JSON key: {key}")
79 result[key] = value
80 return result
83def make_sedml(config, names, sbml_namespace):
84 import libsedml
85 namespaces = libsedml.SedNamespaces(1, 3)
86 namespaces.getNamespaces().add(sbml_namespace, "sbml")
87 sed = libsedml.SedDocument(namespaces)
88 simulation = sed.createUniformTimeCourse(); simulation.setId("simulation")
89 simulation.setInitialTime(0.); simulation.setOutputStartTime(0.)
90 simulation.setOutputEndTime(config["end_time"]); simulation.setNumberOfPoints(config["points"] - 1)
91 algorithm = simulation.createAlgorithm(); algorithm.setKisaoID("KISAO:0000019")
92 # RoadRunner's scalar absolute_tolerance scales a state/amount tolerance vector.
93 # KiSAO 571 captures that adjustment-factor meaning; 211 would claim a fixed tolerance.
94 for kisao, value in (("KISAO:0000209", config["relative_tolerance"]),
95 ("KISAO:0000571", config["absolute_tolerance"]),
96 ("KISAO:0000671", "true"), ("KISAO:0000656", "false")):
97 parameter = algorithm.createAlgorithmParameter(); parameter.setKisaoID(kisao); parameter.setValue(str(value))
98 for scenario in names:
99 model = sed.createModel(); model.setId(f"model_{scenario}")
100 model.setSource(f"model_{scenario}.xml"); model.setLanguage("urn:sedml:language:sbml")
101 task = sed.createTask(); task.setId(f"task_{scenario}")
102 task.setModelReference(f"model_{scenario}"); task.setSimulationReference("simulation")
103 report = sed.createReport(); report.setId(f"report_{scenario}")
104 for index, symbol in enumerate(["time"] + config["species"]):
105 generator = sed.createDataGenerator(); generator.setId(f"dg_{scenario}_{index}")
106 variable = generator.createVariable(); variable.setId(f"v_{scenario}_{index}")
107 variable.setTaskReference(f"task_{scenario}")
108 if index == 0:
109 variable.setSymbol("urn:sedml:symbol:time")
110 else:
111 variable.setTarget(f"/sbml:sbml/sbml:model/sbml:listOfSpecies/sbml:species[@id='{symbol}']")
112 generator.setMath(libsedml.parseFormula(variable.getId()))
113 dataset = report.createDataSet(); dataset.setId(f"ds_{scenario}_{index}")
114 dataset.setLabel(symbol); dataset.setDataReference(generator.getId())
115 return sed
118def run(model_path, model_format, experiment_path, output):
119 import libsbml
120 import libsedml
121 import numpy as np
122 import tellurium as te
123 model_path, experiment_path, output = Path(model_path), Path(experiment_path), Path(output)
124 if output.exists():
125 raise FileExistsError(f"Refusing to overwrite {output}")
126 document, issues = model_document(model_path, model_format)
127 config = validate_experiment(json.loads(experiment_path.read_text(encoding="utf-8"),
128 object_pairs_hook=unique_json_object), document.getModel())
129 output.mkdir(parents=True)
130 selections = ["time"] + [f"[{s}]" for s in config["species"]]
131 results, model_files, initial_tolerances = {}, [], {}
132 for scenario, changes in config["scenarios"].items():
133 # Fresh SBML per condition avoids reset/state carryover ambiguity.
134 condition = document.clone()
135 for symbol, value in changes.items():
136 condition.getModel().getParameter(symbol).setValue(value)
137 sbml = libsbml.writeSBMLToString(condition)
138 model_file = output / f"model_{scenario}.xml"; model_file.write_text(sbml); model_files.append(model_file)
139 runner = te.loadSBMLModel(sbml)
140 runner.conservedMoietyAnalysis = False
141 runner.setIntegrator("cvode")
142 integrator = runner.getIntegrator()
143 integrator.relative_tolerance = config["relative_tolerance"]
144 integrator.absolute_tolerance = config["absolute_tolerance"]
145 integrator.stiff = True
146 integrator.variable_step_size = False
147 initial_tolerances[scenario] = {
148 runner.model.getStateVectorId(i): float(value)
149 for i, value in enumerate(integrator.getAbsoluteToleranceVector())}
150 values = np.asarray(runner.simulate(start=0., end=config["end_time"], points=config["points"],
151 selections=selections))
152 if values.shape != (config["points"], len(selections)) or not np.isfinite(values).all():
153 raise RuntimeError(f"{scenario}: invalid simulation output")
154 results[scenario] = values
155 with (output / f"{scenario}.csv").open("w", newline="") as handle:
156 writer = csv.writer(handle); writer.writerow(selections); writer.writerows(values)
157 sed = make_sedml(config, config["scenarios"], document.getModel().getURI())
158 sed_path = output / "experiment.sedml"
159 libsedml.writeSedMLToFile(sed, str(sed_path))
160 parsed_sed = libsedml.readSedMLFromFile(str(sed_path))
161 sed_issues = [{"severity": parsed_sed.getError(i).getSeverityAsString(), "message": parsed_sed.getError(i).getMessage()}
162 for i in range(parsed_sed.getNumErrors())]
163 if any(issue["severity"] in {"Error", "Fatal"} for issue in sed_issues):
164 raise ValueError("Generated SED-ML failed validation")
165 archive = output / "experiment.omex"
166 te.createCombineArchive(str(archive), [str(p) for p in model_files] + [str(sed_path)],
167 [p.name for p in model_files] + [sed_path.name],
168 ["http://identifiers.org/combine.specifications/sbml"] * len(model_files) + ["http://identifiers.org/combine.specifications/sed-ml"],
169 [False] * len(model_files) + [True])
170 # Replay only this generated archive; no arbitrary archive extraction entry point.
171 replay = te.executeCombineArchive(str(archive), createOutputs=False)
172 if len(replay) != 1:
173 raise RuntimeError("Expected one SED-ML document in generated archive")
174 generators = next(iter(replay.values()))["dataGenerators"]
175 replay_errors = {}
176 for scenario, values in results.items():
177 replayed = np.column_stack([generators[f"dg_{scenario}_{i}"].reshape(-1) for i in range(len(selections))])
178 np.testing.assert_allclose(replayed, values, rtol=config["relative_tolerance"] * 10,
179 atol=config["absolute_tolerance"] * 10)
180 replay_errors[scenario] = float(np.max(np.abs(replayed - values)))
181 model = document.getModel()
182 units = {"time": model.getTimeUnits(), "substance": model.getSubstanceUnits(), "volume": model.getVolumeUnits(),
183 "species": {s: libsbml.UnitDefinition.printUnits(model.getSpecies(s).getDerivedUnitDefinition()) for s in config["species"]}}
184 report = {"versions": {**dict(te.getVersionInfo()), "libcombine": importlib.metadata.version("python-libcombine"),
185 "python": platform.python_version(), "numpy": np.__version__}, "integrator": "cvode", "stiff": True, "experiment": config,
186 "absolute_tolerance_semantics": "RoadRunner scalar adjustment factor, not fixed concentration error",
187 "initial_state_absolute_tolerances": initial_tolerances,
188 "source_model_sha256": hashlib.sha256(model_path.read_bytes()).hexdigest(),
189 "experiment_sha256": hashlib.sha256(experiment_path.read_bytes()).hexdigest(),
190 "sbml_validation_findings": issues, "sedml_parse_findings": sed_issues, "units": units,
191 "archive_sha256": hashlib.sha256(archive.read_bytes()).hexdigest(),
192 "archive_replay_max_absolute_difference": replay_errors,
193 "minimum_concentration": {s: float(v[:, 1:].min()) for s, v in results.items()}}
194 (output / "report.json").write_text(json.dumps(report, indent=2) + "\n")
195 return report
198def main():
199 parser = argparse.ArgumentParser(description=__doc__)
200 parser.add_argument("--model", type=Path, required=True)
201 parser.add_argument("--format", choices=["sbml", "antimony"], required=True)
202 parser.add_argument("--experiment", type=Path, required=True)
203 parser.add_argument("--output", type=Path, required=True)
204 args = parser.parse_args()
205 print(json.dumps(run(args.model, args.format, args.experiment, args.output), indent=2))
208if __name__ == "__main__":
209 main()