From 42e6a7e93aa899fdaeaf2a245608f3f0e3ba612a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matthias=20K=C3=B6nig?= Date: Tue, 2 Apr 2024 21:06:39 +0200 Subject: [PATCH] starting simulator update --- setup.cfg | 44 ++-- src/sbmlsim/combine/sedml/runner.py | 1 - src/sbmlsim/examples/datagenerator_example.py | 4 +- src/sbmlsim/examples/example_parallel.py | 190 ------------------ .../examples/experiments/covid/simulate.py | 1 - .../experiments/curve_types/experiment.py | 4 +- src/sbmlsim/examples/experiments/demo/demo.py | 4 +- .../initial_assignment/initial_assignment.py | 4 +- .../experiments/midazolam/simulate.py | 4 +- .../repressilator/repressilator_scans.py | 4 +- .../examples/multiprocessing/ray_example.py | 55 ----- src/sbmlsim/experiment/experiment.py | 2 +- src/sbmlsim/experiment/runner.py | 12 +- src/sbmlsim/model/model.py | 122 +---------- src/sbmlsim/model/model_roadrunner.py | 86 ++++---- src/sbmlsim/simulator/__init__.py | 1 - src/sbmlsim/simulator/simulation.py | 135 ------------- src/sbmlsim/simulator/simulation_ray.py | 172 ---------------- src/sbmlsim/simulator/simulation_serial.py | 188 ++++++++++++----- src/sbmlsim/test/test_simulator.py | 27 --- 20 files changed, 222 insertions(+), 838 deletions(-) delete mode 100644 src/sbmlsim/examples/example_parallel.py delete mode 100644 src/sbmlsim/examples/multiprocessing/ray_example.py delete mode 100644 src/sbmlsim/simulator/simulation.py delete mode 100644 src/sbmlsim/simulator/simulation_ray.py diff --git a/setup.cfg b/setup.cfg index 803d2761..d1dcb986 100644 --- a/setup.cfg +++ b/setup.cfg @@ -46,15 +46,15 @@ keywords = [options] zip_safe = True -python_requires = >=3.9 +python_requires = >=3.11 install_requires = sbmlutils>=0.8.7 - numpy==1.21.3 # see https://github.com/sys-bio/roadrunner/issues/963 - libroadrunner==2.1.3 # see https://github.com/sys-bio/roadrunner/issues/963 - scipy==1.10.1 - pint==0.19.2 + numpy>=1.26.4 + libroadrunner>=2.6.0 + scipy>=1.12.0 + pint>=0.23 - # -> sbmlutils (no version pin) + # from sbmlutils (no version pin) python-libsbml rich depinfo @@ -66,28 +66,22 @@ install_requires = matplotlib # dependencies - petab>=0.1.25 - # fix for windows access violations - ray[default]>=2.9.3 - grpcio==1.49.1 - pydantic>=2.6.3 - tables>=3.8.0 - xarray>=2023.6.0 - bottleneck>=1.3.4 + petab>=0.2.9 + pydantic>=2.6.4 + tables>=3.9.2 + xarray>=2024.3.0 + bottleneck>=1.3.8 + psutil>=5.9.8 + setproctitle>=1.3.3 + + plotly>=5.20.0 + altair>=5.3.0 + seaborn>=0.13.2 + xmltodict>=0.13.0 + pyDOE>=0.3.8 python-libsedml>=2.0.32 python-libnuml>=1.1.6 - # biosimulators-utils>=0.1.171 - - psutil>=5.9.0 - setproctitle>=1.2.2 - - plotly>=5.6.0 - altair>=4.2.0 - seaborn>=0.11.2 - - xmltodict>=0.12.0 - pyDOE>=0.3.8 tests_require = tox>=3.24.3 pytest>=7.0.1 diff --git a/src/sbmlsim/combine/sedml/runner.py b/src/sbmlsim/combine/sedml/runner.py index 1ac56e14..2d6b4380 100644 --- a/src/sbmlsim/combine/sedml/runner.py +++ b/src/sbmlsim/combine/sedml/runner.py @@ -10,7 +10,6 @@ from sbmlsim.combine.sedml.parser import SEDMLParser from sbmlsim.experiment import ExperimentRunner, SimulationExperiment from sbmlsim.simulator import SimulatorSerial -from sbmlsim.simulator.simulation_ray import SimulatorParallel def sedmltojson(sedml_path: Path) -> None: diff --git a/src/sbmlsim/examples/datagenerator_example.py b/src/sbmlsim/examples/datagenerator_example.py index 85594d61..58f2ca35 100644 --- a/src/sbmlsim/examples/datagenerator_example.py +++ b/src/sbmlsim/examples/datagenerator_example.py @@ -4,13 +4,13 @@ from sbmlsim.result import XResult from sbmlsim.result.datagenerator import DataGeneratorIndexingFunction from sbmlsim.simulation import Dimension, ScanSim, Timecourse, TimecourseSim -from sbmlsim.simulator.simulation_ray import SimulatorParallel +from sbmlsim.simulator.simulation_serial import SimulatorSerial from sbmlsim.test import MODEL_MIDAZOLAM def example_scan() -> XResult: """Run scan and return results.""" - simulator = SimulatorParallel(model=MODEL_MIDAZOLAM) + simulator = SimulatorSerial(model=MODEL_MIDAZOLAM) Q_ = simulator.Q_ scan = ScanSim( diff --git a/src/sbmlsim/examples/example_parallel.py b/src/sbmlsim/examples/example_parallel.py deleted file mode 100644 index c031f417..00000000 --- a/src/sbmlsim/examples/example_parallel.py +++ /dev/null @@ -1,190 +0,0 @@ -""" -Parallel execution of timecourses -""" -import tempfile -import time - -import numpy as np -import roadrunner - -from sbmlsim.model import RoadrunnerSBMLModel -from sbmlsim.simulation import Dimension, ScanSim, Timecourse, TimecourseSim -from sbmlsim.simulator import SimulatorSerial -from sbmlsim.simulator.simulation_ray import ( - SimulatorActor, - SimulatorParallel, - cpu_count, - ray, -) -from sbmlsim.test import MODEL_GLCWB, MODEL_REPRESSILATOR -from sbmlsim.units import Units, UnitsInformation - - -def example_single_actor(): - """Creates a single stateful simulator actor and executes timecourse. - - Actors should never created manually ! - Use the SimulatorParallel to run parallel simulations. - - :return: - """ - # create state file - r = roadrunner.RoadRunner(str(MODEL_REPRESSILATOR)) - RoadrunnerSBMLModel.set_timecourse_selections(r) - uinfo = UnitsInformation.from_sbml(MODEL_REPRESSILATOR) - - f_state = tempfile.NamedTemporaryFile(suffix=".dat") - r.saveState(f_state.name) - - # Create single actor process - sa = SimulatorActor.remote(f_state.name) - - # run simulation - tcsim = TimecourseSim( - [ - Timecourse(start=0, end=100, steps=100), - Timecourse(start=0, end=100, steps=100, changes={"X": 10, "Y": 20}), - ] - ) - - tcsim.normalize(uinfo=uinfo) - tc_id = sa._timecourse.remote(tcsim) - results = ray.get(tc_id) - return results - - -def example_multiple_actors(): - """Multiple independent simulator actors. - - Actors should never be created manually, use the SimulatorParallel! - """ - # create state file - r = roadrunner.RoadRunner(str(MODEL_REPRESSILATOR)) - RoadrunnerSBMLModel.set_timecourse_selections(r) - uinfo = UnitsInformation.from_sbml(MODEL_REPRESSILATOR) - - f_state = tempfile.NamedTemporaryFile(suffix=".dat") - r.saveState(f_state.name) - - # create ten Simulators. - simulators = [SimulatorActor.remote(f_state.name) for _ in range(16)] - - # define timecourse - tcsim = TimecourseSim( - [ - Timecourse(start=0, end=100, steps=100), - Timecourse(start=0, end=100, steps=100, changes={"X": 10, "Y": 20}), - ] - ) - tcsim.normalize(uinfo=uinfo) - - # run simulation on simulators - tc_ids = [s._timecourse.remote(tcsim) for s in simulators] - # collect results - results = ray.get(tc_ids) - return results - - -def example_parallel_timecourse(nsim=40, actor_count=15): - """Execute multiple simulations with model in parallel. - - :param nsim: number of simulations - :return: - """ - # collect all simulation definitions (see also the ensemble functions) - - scan = ScanSim( - simulation=TimecourseSim( - [ - Timecourse( - start=0, - end=100, - steps=100, - changes={ - "IVDOSE_som": 0.0, # [mg] - "PODOSE_som": 0.0, # [mg] - "Ri_som": 10.0e-6, # [mg/min] - }, - ) - ] - ), - dimensions=[Dimension("dim1", index=np.arange(nsim))], - ) - - def message(info, time): - print(f"{info:<10}: {time:4.3f}") - - # load model once for caching (fair comparison) - r = roadrunner.RoadRunner(str(MODEL_GLCWB)) - # simulator definitions - simulator_defs = [ - { - "key": "parallel", - "simulator": SimulatorParallel, - "kwargs": {"actor_count": actor_count}, - }, - {"key": "serial", "simulator": SimulatorSerial, "kwargs": {}}, - ] - print("-" * 80) - print(f"Run '{nsim}' simulations") - print("-" * 80) - sim_info = [] - for sim_def in simulator_defs: - key = sim_def["key"] - print("***", key, "***") - Simulator = sim_def["simulator"] - kwargs = sim_def["kwargs"] - - # run simulation (with model reading) - start_time = time.time() - # create a simulator with 16 parallel actors - simulator = Simulator(model=MODEL_GLCWB, **kwargs) - load_time = time.time() - start_time - message(f"load", load_time) - - start_time = time.time() - results = simulator.run_scan(scan=scan) - sim_time = time.time() - start_time - total_time = load_time + sim_time - message("simulate", sim_time) - message("total", total_time) - assert len(results.coords["dim1"]) == nsim - - # run parallel simulation (without model reading) - start_time = time.time() - results = simulator.run_scan(scan=scan) - repeat_time = time.time() - start_time - message(f"repeat", repeat_time) - assert len(results.coords["dim1"]) == nsim - - actor_count = kwargs.get("actor_count", 1) - times = { - "load": load_time, - "simulate": sim_time, - "total": total_time, - "repeat": repeat_time, - } - sim_info.extend( - [ - { - "key": key, - "nsim": nsim, - "actor_count": actor_count, - "time_type": k, - "time": v, - } - for (k, v) in times.items() - ] - ) - - print("-" * 80) - - return sim_info - - -if __name__ == "__main__": - # example_single_actor() - # example_multiple_actors() - - sim_info = example_parallel_timecourse(nsim=100, actor_count=15) - ray.timeline(filename="timeline.json") diff --git a/src/sbmlsim/examples/experiments/covid/simulate.py b/src/sbmlsim/examples/experiments/covid/simulate.py index 62a2fc71..9b9d0139 100644 --- a/src/sbmlsim/examples/experiments/covid/simulate.py +++ b/src/sbmlsim/examples/experiments/covid/simulate.py @@ -24,7 +24,6 @@ def run_covid_examples(output_path: Path) -> None: output_path=output_path / "sbmlsim", data_path=output_path, base_path=output_path, - parallel=True, ) for experiment in experiments: diff --git a/src/sbmlsim/examples/experiments/curve_types/experiment.py b/src/sbmlsim/examples/experiments/curve_types/experiment.py index 81d1cb53..0f359c35 100644 --- a/src/sbmlsim/examples/experiments/curve_types/experiment.py +++ b/src/sbmlsim/examples/experiments/curve_types/experiment.py @@ -10,7 +10,7 @@ from sbmlsim.plot import Axis, Figure, Plot from sbmlsim.result.report import Report from sbmlsim.simulation import AbstractSim, Timecourse, TimecourseSim -from sbmlsim.simulator.simulation_ray import SimulatorParallel, SimulatorSerial +from sbmlsim.simulator.simulation_serial import SimulatorSerial from sbmlsim.task import Task from sbmlsim.test import MODEL_REPRESSILATOR @@ -95,7 +95,7 @@ def run_curve_types_experiments(output_path: Path) -> Path: runner = ExperimentRunner( CurveTypesExperiment, - simulator=SimulatorParallel(), + simulator=SimulatorSerial(), data_path=data_path, base_path=base_path, ) diff --git a/src/sbmlsim/examples/experiments/demo/demo.py b/src/sbmlsim/examples/experiments/demo/demo.py index 61de108a..2c3760d0 100644 --- a/src/sbmlsim/examples/experiments/demo/demo.py +++ b/src/sbmlsim/examples/experiments/demo/demo.py @@ -20,7 +20,7 @@ TimecourseSim, ) from sbmlsim.simulation.sensitivity import ModelSensitivity, SensitivityType -from sbmlsim.simulator.simulation_ray import SimulatorParallel, SimulatorSerial +from sbmlsim.simulator.simulation_serial import SimulatorSerial from sbmlsim.task import Task from sbmlsim.test import MODEL_DEMO @@ -106,7 +106,7 @@ def run_demo_experiments(output_path: Path) -> None: runner = ExperimentRunner( DemoExperiment, - simulator=SimulatorParallel(), + simulator=SimulatorSerial(), data_path=data_path, base_path=base_path, ) diff --git a/src/sbmlsim/examples/experiments/initial_assignment/initial_assignment.py b/src/sbmlsim/examples/experiments/initial_assignment/initial_assignment.py index 2036302f..e55e254b 100644 --- a/src/sbmlsim/examples/experiments/initial_assignment/initial_assignment.py +++ b/src/sbmlsim/examples/experiments/initial_assignment/initial_assignment.py @@ -8,7 +8,7 @@ from sbmlsim.model import AbstractModel, RoadrunnerSBMLModel from sbmlsim.plot import Axis, Figure from sbmlsim.simulation import AbstractSim, Timecourse, TimecourseSim -from sbmlsim.simulator.simulation_ray import SimulatorParallel +from sbmlsim.simulator.simulation_serial import SimulatorSerial from sbmlsim.task import Task @@ -119,7 +119,7 @@ def run(output_path): runner = ExperimentRunner( AssignmentExperiment, - simulator=SimulatorParallel(), + simulator=SimulatorSerial(), base_path=base_path, data_path=base_path, ) diff --git a/src/sbmlsim/examples/experiments/midazolam/simulate.py b/src/sbmlsim/examples/experiments/midazolam/simulate.py index 81c8773b..96603b7f 100644 --- a/src/sbmlsim/examples/experiments/midazolam/simulate.py +++ b/src/sbmlsim/examples/experiments/midazolam/simulate.py @@ -12,7 +12,7 @@ from sbmlsim.examples.experiments.midazolam.experiments.mandema1992 import Mandema1992 from sbmlsim.experiment import ExperimentRunner, SimulationExperiment from sbmlsim.report.experiment_report import ExperimentReport, ReportResults -from sbmlsim.simulator.simulation_ray import SimulatorParallel +from sbmlsim.simulator.simulation_serial import SimulatorSerial def run_midazolam_experiments(output_path: Path) -> None: @@ -23,7 +23,7 @@ def run_midazolam_experiments(output_path: Path) -> None: Mandema1992, Kupferschmidt1995, ], - simulator=SimulatorParallel(), + simulator=SimulatorSerial(), base_path=base_path, data_path=base_path / "data", ) diff --git a/src/sbmlsim/examples/experiments/repressilator/repressilator_scans.py b/src/sbmlsim/examples/experiments/repressilator/repressilator_scans.py index e798a462..97d2e804 100644 --- a/src/sbmlsim/examples/experiments/repressilator/repressilator_scans.py +++ b/src/sbmlsim/examples/experiments/repressilator/repressilator_scans.py @@ -17,7 +17,7 @@ Timecourse, TimecourseSim, ) -from sbmlsim.simulator.simulation_ray import SimulatorParallel, SimulatorSerial +from sbmlsim.simulator.simulation_serial import SimulatorSerial from sbmlsim.task import Task from sbmlsim.test import MODEL_REPRESSILATOR @@ -219,7 +219,7 @@ def run_repressilator_experiments(output_path: Path) -> Path: base_path = Path(__file__).parent data_path = base_path - for simulator in [SimulatorSerial(), SimulatorParallel()]: + for simulator in [SimulatorSerial()]: runner = ExperimentRunner( [RepressilatorScanExperiment], simulator=simulator, diff --git a/src/sbmlsim/examples/multiprocessing/ray_example.py b/src/sbmlsim/examples/multiprocessing/ray_example.py deleted file mode 100644 index 933a4dee..00000000 --- a/src/sbmlsim/examples/multiprocessing/ray_example.py +++ /dev/null @@ -1,55 +0,0 @@ -"""Example parallel simulation with ray. - -requirements - roadrunner - ray - pandas -""" -import time - -import pandas as pd -import ray -import roadrunner - - -# start ray -ray.init(ignore_reinit_error=True) - - -@ray.remote -class SimulatorActorPath(object): - """Ray actor to execute simulations.""" - - def __init__(self, r: roadrunner.RoadRunner): - self.r: roadrunner.RoadRunner = r - - def simulate(self, size=1): - """Simulate.""" - print("Start simulations") - ts = time.time() - results = [] - for _ in range(size): - self.r.resetAll() - s = self.r.simulate(0, 100, steps=100) - # create numpy array (which can be pickled), better use shared memory - df = pd.DataFrame(s, columns=s.colnames) - results.append(df) - te = time.time() - print("Finished '{}' simulations: {:2.2f} ms".format(size, (te - ts) * 1000)) - return results - - -if __name__ == "__main__": - actor_count = 10 # cores to run this on - - rr = roadrunner.RoadRunner("icg_body_flat.xml") - simulators = [SimulatorActorPath.remote(rr) for _ in range(actor_count)] - - # run simulations - sim_per_actor = 10 - tc_ids = [] - for simulator in simulators: - tcs_id = simulator.simulate.remote(size=sim_per_actor) - tc_ids.append(tcs_id) - - results = ray.get(tc_ids) diff --git a/src/sbmlsim/experiment/experiment.py b/src/sbmlsim/experiment/experiment.py index 9eb07c74..595b4e26 100644 --- a/src/sbmlsim/experiment/experiment.py +++ b/src/sbmlsim/experiment/experiment.py @@ -451,7 +451,7 @@ def _run_tasks(self, simulator, reduced_selections: bool = True): # load model in simulator model: AbstractModel = self._models[model_id] - logger.info("set model") + logger.info(f"set model: {type(model)}: {model=}") simulator.set_model(model=model) logger.info("set selections") diff --git a/src/sbmlsim/experiment/runner.py b/src/sbmlsim/experiment/runner.py index 8824541f..ce8e02b1 100644 --- a/src/sbmlsim/experiment/runner.py +++ b/src/sbmlsim/experiment/runner.py @@ -17,7 +17,6 @@ from sbmlsim.model import RoadrunnerSBMLModel from sbmlsim.report.experiment_report import ExperimentReport, ReportResults from sbmlsim.simulator import SimulatorSerial -from sbmlsim.simulator.simulation_ray import SimulatorParallel from sbmlsim.units import UnitRegistry, UnitsInformation from sbmlsim.utils import timeit @@ -104,15 +103,14 @@ def initialize(self, experiment_classes: Union[List[Type[SimulationExperiment]], for model_id, source in experiment.models().items(): if source not in self.models: # not cashed yet, cash the model for lookup - self.models[source] = RoadrunnerSBMLModel( - source=source, ureg=self.ureg + self.models[source] = RoadrunnerSBMLModel.from_abstract_model( + abstract_model=source, ureg=self.ureg ) _models[model_id] = self.models[source] # set resolved models in experiment experiment._models = _models # only after model loading the unit registry is filled - experiment.initialize() self.experiments[experiment.sid] = experiment @@ -130,7 +128,8 @@ def run_experiments( output_path.mkdir(parents=True) exp_results = [] - for sid, experiment in self.experiments.items(): # type: SimulationExperiment + experiment: SimulationExperiment + for sid, experiment in self.experiments.items(): logger.info(f"Running SimulationExperiment: {sid}") # ExperimentResult used to create report @@ -151,12 +150,11 @@ def run_experiments( output_path: Path, base_path: Path = None, data_path: Union[List[Path], Tuple[Path], Optional[Path]] = None, - parallel: bool = True, ) -> Path: """Run simulation experiments.""" if not isinstance(experiments, (list, tuple)): experiments = [experiments] - simulator = SimulatorParallel() if parallel else SimulatorSerial() + simulator = SimulatorSerial() runner = ExperimentRunner( experiments, diff --git a/src/sbmlsim/model/model.py b/src/sbmlsim/model/model.py index 05917e3a..f9573b40 100644 --- a/src/sbmlsim/model/model.py +++ b/src/sbmlsim/model/model.py @@ -18,7 +18,7 @@ logger = log.get_logger(__name__) -class AbstractModel(object): +class AbstractModel: """Abstract base class to store a model in sbmlsim. Depending on the model language different subclasses are implemented. @@ -99,123 +99,3 @@ def to_dict(self): "changes": self.changes, } return d - - ''' - # SED-ML HANDLING - - def apply_change(self, target, value): - """Applies change to model""" - if target.startswith("/"): - # xpath expression - target = AbstractModel._resolve_xpath(self._model, target) - - - def apply_model_changes(self, changes): - """Applies dictionary of model changes.""" - for key, value in self.changes.items(): - AbstractModel.apply_change(target=key, value=value) - - - @staticmethod - def _target_from_xpath(model: 'AbstractModel', xpath: str): - """ Resolve the target from the xpath expression. - - A single target in the model corresponding to the modelId is resolved. - Currently, the model is not used for xpath resolution. - - :param xpath: xpath expression. - :type xpath: str - :param modelId: id of model in which xpath should be resolved - :type modelId: str - :return: single target of xpath expression - :rtype: Target (namedtuple: id type) - """ - # TODO: via better xpath expression - # get type from the SBML document for the given id. - # The xpath expression can be very general and does not need to contain the full - # xml path - # For instance: - # /sbml:sbml/sbml:model/descendant::*[@id='S1'] - # has to resolve to species. - # TODO: figure out concentration or amount (from SBML document) - # FIXME: getting of sids, pids not very robust, handle more cases (rules, reactions, ...) - - Target = namedtuple('Target', 'id type') - - def getId(xpath): - xpath = xpath.replace('"', "'") - match = re.findall(r"id='(.*?)'", xpath) - if (match is None) or (len(match) is 0): - logger.warn("Xpath could not be resolved: {}".format(xpath)) - return match[0] - - # parameter value change - if ("model" in xpath) and ("parameter" in xpath): - return Target(getId(xpath), 'parameter') - # species concentration change - elif ("model" in xpath) and ("species" in xpath): - return Target(getId(xpath), 'concentration') - # other - elif ("model" in xpath) and ("id" in xpath): - return Target(getId(xpath), 'other') - # cannot be parsed - else: - raise ValueError("Unsupported target in xpath: {}".format(xpath)) - - @staticmethod - def set_xpath_value(xpath: str, value: float, model): - """ Creates python line for given xpath target and value. - :param xpath: - :type xpath: - :param value: - :type value: - :return: - :rtype: - """ - target = SEDMLParser._resolve_xpath(xpath) - if target: - if target.type == "concentration": - # initial concentration - expr = f'init([{target.id}])' - elif target.type == "amount": - # initial amount - expr = f'init({target.id})' - else: - # other (parameter, flux, ...) - expr = target.id - print(f"{expr} = {value}") - model[expr] = value - else: - logger.error(f"Unsupported target xpath: {xpath}") - - @staticmethod - def selectionFromVariable(var, model): - """ Resolves the selection for the given variable. - - First checks if the variable is a symbol and returns the symbol. - If no symbol is set the xpath of the target is resolved - and used in the selection - - :param var: variable to resolve - :type var: SedVariable - :return: a single selection - :rtype: Selection (namedtuple: id type) - """ - Selection = namedtuple('Selection', 'id type') - - # parse symbol expression - if var.isSetSymbol(): - cvs = var.getSymbol() - astr = cvs.rsplit("symbol:") - sid = astr[1] - return Selection(sid, 'symbol') - # use xpath - elif var.isSetTarget(): - xpath = var.getTarget() - target = SEDMLParser._resolveXPath(xpath, model) - return Selection(target.id, target.type) - - else: - warnings.warn(f"Unrecognized Selection in variable: {var}") - return None - ''' diff --git a/src/sbmlsim/model/model_roadrunner.py b/src/sbmlsim/model/model_roadrunner.py index 2a9c9f65..6d8da532 100644 --- a/src/sbmlsim/model/model_roadrunner.py +++ b/src/sbmlsim/model/model_roadrunner.py @@ -30,7 +30,7 @@ class RoadrunnerSBMLModel(AbstractModel): def __init__( self, - source: Union[str, Path, AbstractModel], + source: Union[str, Path], base_path: Path = None, changes: Dict = None, sid: str = None, @@ -39,66 +39,63 @@ def __init__( ureg: UnitRegistry = None, settings: Dict = None, ): - logger.debug(f"source: {type(source)}, {source}") - if isinstance(source, AbstractModel): - logger.debug("RoadrunnerSBMLModel from AbstractModel") - super(RoadrunnerSBMLModel, self).__init__( - source=source.source, - language_type=source.language_type, - changes=source.changes, - sid=source.sid, - name=source.name, - base_path=source.base_path, - selections=selections, - ) - else: - logger.info("RoadrunnerSBMLModel from source") - super(RoadrunnerSBMLModel, self).__init__( - source=source, - language_type=AbstractModel.LanguageType.SBML, - changes=changes, - sid=sid, - name=name, - base_path=base_path, - selections=selections, - ) + super(RoadrunnerSBMLModel, self).__init__( + source=source, + language_type=AbstractModel.LanguageType.SBML, + changes=changes, + sid=sid, + name=name, + base_path=base_path, + selections=selections, + ) + + # check SBML if self.language_type != AbstractModel.LanguageType.SBML: - raise ValueError( - f"{self.__class__.__name__} only supports " - f"language_type '{AbstractModel.LanguageType.SBML}'." - ) + raise ValueError(f"language_type not supported '{self.language_type}'.") - # load the model - self.state_path = self.get_state_path() - logger.debug(f"Load model from state: {self.state_path}") - self._model = self.load_roadrunner_model( + # load model + self.state_path: Path = self.get_state_path() + self.r: Optional[roadrunner.RoadRunner] = self.load_roadrunner_model( source=self.source, state_path=self.state_path ) + # set selections self.selections = self.set_timecourse_selections( - self._model, selections=self.selections + self.r, selections=self.selections ) # set integrator settings - if settings is not None: - RoadrunnerSBMLModel.set_integrator_settings(self._model, **settings) - - self.uinfo = self.parse_units(ureg) + if settings: + RoadrunnerSBMLModel.set_integrator_settings(self.r, **settings) # normalize model changes + self.uinfo = self.parse_units(ureg) self.normalize(uinfo=self.uinfo) - logger.debug(f"model.changes: {self.changes}") - @property def Q_(self) -> Quantity: """Quantity to create quantities for model changes.""" return self.uinfo.ureg.Quantity - @property - def r(self) -> roadrunner.RoadRunner: - """Roadrunner instance.""" - return self._model + @staticmethod + def from_abstract_model( + abstract_model: AbstractModel, + selections: List[str] = None, + ureg: UnitRegistry = None, + settings: Dict = None + ): + """Create from AbstractModel.""" + logger.debug("RoadrunnerSBMLModel from AbstractModel") + return RoadrunnerSBMLModel( + source=abstract_model.source.source, + changes=abstract_model.changes, + sid=abstract_model.sid, + name=abstract_model.name, + base_path=abstract_model.base_path, + selections=selections, + ureg=ureg, + settings=settings, + ) def get_state_path(self) -> Optional[Path]: """Get path of the state file. @@ -107,8 +104,6 @@ def get_state_path(self) -> Optional[Path]: """ if self.source.is_path(): md5 = md5_for_path(self.source.path) - # FIXME: get unique hash for library version - return Path(f"{self.source.path}_rr{roadrunner.__version__}_{md5}.state") else: return None @@ -140,6 +135,7 @@ def load_roadrunner_model( if state_path: r.saveState(str(state_path)) logger.info(f"Save state: '{state_path}'") + elif source.is_content(): r = roadrunner.RoadRunner(str(source.content)) diff --git a/src/sbmlsim/simulator/__init__.py b/src/sbmlsim/simulator/__init__.py index a1aa685d..b68fa2e8 100644 --- a/src/sbmlsim/simulator/__init__.py +++ b/src/sbmlsim/simulator/__init__.py @@ -1,3 +1,2 @@ """Package for simulator.""" -from .simulation import SimulatorWorker from .simulation_serial import SimulatorSerial diff --git a/src/sbmlsim/simulator/simulation.py b/src/sbmlsim/simulator/simulation.py deleted file mode 100644 index c27b6341..00000000 --- a/src/sbmlsim/simulator/simulation.py +++ /dev/null @@ -1,135 +0,0 @@ -"""Classes for running simulations with SBML models.""" -from typing import Optional - -import pandas as pd -import roadrunner -from sbmlutils import log - -from sbmlsim.model import ModelChange -from sbmlsim.simulation import Timecourse, TimecourseSim - - -logger = log.get_logger(__name__) - - -class SimulatorWorker: - """Worker running simulations. - - Implements the timecourse simulation once which can be reused by - the different simulators. - """ - - def __init__(self): - """Init placeholder.""" - self.r: Optional[roadrunner.RoadRunner] = None - - def _timecourse(self, simulation: TimecourseSim) -> pd.DataFrame: - """Timecourse simulation. - - Requires for all timecourse definitions in the timecourse simulation - to be unit normalized. The changes have no units any more - for parallel simulations. - You should never call this function directly! - - :param simulation: Simulation definition(s) - :return: DataFrame with results - """ - if isinstance(simulation, Timecourse): - simulation = TimecourseSim(timecourses=[simulation]) - - if simulation.reset: - self.r.resetToOrigin() - - frames = [] - t_offset = simulation.time_offset - for k, tc in enumerate(simulation.timecourses): - - if k == 0 and tc.model_changes: - # [1] apply model changes of first simulation - logger.debug("Applying model changes") - for key, item in tc.model_changes.items(): - if key.startswith("init"): - logger.error( - f"Initial model changes should be provided " - f"without 'init': '{key} = {item}'" - ) - # FIXME: implement model changes via init - # init_key = f"init({key})" - init_key = key - try: - value = item.magnitude - except AttributeError: - value = item - - try: - self.r[init_key] = value - except RuntimeError: - logger.error(f"roadrunner RuntimeError: '{init_key} = {item}'") - # boundary condition=true species, trying direct fallback - # see https://github.com/sys-bio/roadrunner/issues/711 - init_key = key - self.r[key] = value - - logger.debug(f"\t{init_key} = {item}") - - # [2] re-evaluate initial assignments - # https://github.com/sys-bio/roadrunner/issues/710 - # logger.debug("Reevaluate initial conditions") - # FIXME/TODO: support initial model changes - # self.r.resetAll() - # self.r.reset(SelectionRecord.DEPENDENT_FLOATING_AMOUNT) - # self.r.reset(SelectionRecord.DEPENDENT_INITIAL_GLOBAL_PARAMETER) - - # [3] apply model manipulations - # model manipulations are applied to model - if len(tc.model_manipulations) > 0: - # FIXME: update to support roadrunner model changes - for key, value in tc.model_changes.items(): - if key == ModelChange.CLAMP_SPECIES: - for sid, formula in value.items(): - ModelChange.clamp_species(self.r, sid, formula) - else: - raise ValueError( - f"Unsupported model change: " - f"'{key}': {value}. Supported changes are: " - f"['{ModelChange.CLAMP_SPECIES}']" - ) - - # [4] apply changes - if tc.changes: - logger.debug("Applying simulation changes") - for key, item in tc.changes.items(): - - # FIXME: handle concentrations/amounts/default - # TODO: Figure out the hasOnlySubstanceUnit flag! (roadrunner) - # r: roadrunner.ExecutableModel = self.r - - try: - self.r[key] = float(item.magnitude) - except AttributeError: - self.r[key] = float(item) - logger.debug(f"\t{key} = {item}") - - # debug model state - # FIXME: report issue - # sbml_str = self.r.getCurrentSBML() - # with open("/home/mkoenig/git/pkdb_models/pkdb_models/models/dextromethorphan/results/debug/test.xml", "w") as f_out: - # f_out.write(sbml_str) - - # run simulation - integrator = self.r.integrator - # FIXME: support simulation by times - if integrator.getValue("variable_step_size"): - s = self.r.simulate(start=tc.start, end=tc.end) - else: - s = self.r.simulate(start=tc.start, end=tc.end, steps=tc.steps) - - df = pd.DataFrame(s, columns=s.colnames) - df.time = df.time + t_offset - - if not tc.discard: - # discard timecourses (pre-simulation) - t_offset += tc.end - frames.append(df) - - return pd.concat(frames, sort=False) diff --git a/src/sbmlsim/simulator/simulation_ray.py b/src/sbmlsim/simulator/simulation_ray.py deleted file mode 100644 index 7427dfc3..00000000 --- a/src/sbmlsim/simulator/simulation_ray.py +++ /dev/null @@ -1,172 +0,0 @@ -"""Parallel simulation using ray.""" -from pathlib import Path -from typing import Iterator, List - -import numpy as np -import pandas as pd -import psutil -import ray -import roadrunner -from sbmlutils import log -import tempfile - -from sbmlsim.simulation import TimecourseSim -from sbmlsim.simulator import SimulatorSerial -from sbmlsim.simulator.simulation import SimulatorWorker - - -logger = log.get_logger(__name__) - -# start ray -ray.init(ignore_reinit_error=True) - - -def cpu_count() -> int: - """Get physical CPU count.""" - return psutil.cpu_count(logical=False) - - -@ray.remote -class SimulatorActor(SimulatorWorker): - """Ray actor to execute simulations. - - An actor is essentially a stateful worker - """ - - def __init__(self, path_state=None): - """State contains model, integrator settings and selections.""" - self.r = None - if path_state is not None: - self.set_model(path_state) - - def set_model(self, state: bytes) -> None: - """Set model using the loaded state - - Faster to load the state. - """ - self.r: roadrunner.RoadRunner = roadrunner.RoadRunner() - if state is not None: - # write state to temporary file for reading - with tempfile.NamedTemporaryFile("wb") as f_temp: - f_temp.write(state) - filename = f_temp.name - self.r.loadState(str(filename)) - - def set_timecourse_selections(self, selections: Iterator[str]): - """Set the timecourse selections.""" - logger.info(f"'set_timecourse_selections':{selections}") - try: - if selections is None: - r_model: roadrunner.ExecutableModel = self.r.model - self.r.timeCourseSelections = ( - ["time"] - + r_model.getFloatingSpeciesIds() - + r_model.getBoundarySpeciesIds() - + r_model.getGlobalParameterIds() - + r_model.getReactionIds() - + r_model.getCompartmentIds() - ) - self.r.timeCourseSelections += [ - f"[{key}]" - for key in ( - r_model.getFloatingSpeciesIds() - + r_model.getBoundarySpeciesIds() - ) - ] - else: - self.r.timeCourseSelections = selections - except RuntimeError as err: - logger.error(f"{err}") - raise (err) - - def work(self, simulations): - """Run a bunch of simulations on a single worker.""" - results = [] - for tc_sim in simulations: - results.append(self._timecourse(tc_sim)) - return results - - -class SimulatorParallel(SimulatorSerial): - """Parallel simulator. - - The parallel simulator is a subclass of the SimulatorSerial reusing the - logic for running simulations. - """ - - def __init__(self, model=None, **kwargs): - """Initialize parallel simulator with multiple workers. - - :param model: model source or model - :param actor_count: int, - """ - if "actor_count" in kwargs: - self.actor_count = kwargs.pop("actor_count") - else: - - # FIXME: get virtual cores - self.actor_count = max(cpu_count() - 1, 1) - logger.info(f"Using '{self.actor_count}' cpu/core for parallel simulation.") - - # Create actors once - logger.debug(f"Creating '{self.actor_count}' SimulationActors") - self.simulators = [SimulatorActor.remote() for _ in range(self.actor_count)] - - super(SimulatorParallel, self).__init__(model=None, **kwargs) - if model is not None: - self.set_model(model=model) - - def set_model(self, model): - """Set model.""" - super(SimulatorParallel, self).set_model(model) - if model: - if not self.model.state_path: - raise ValueError("State path does not exist.") - - # read state only once - with open(self.model.state_path, "rb") as f_state: - state = f_state.read() - for simulator in self.simulators: - simulator.set_model.remote(state) - self.set_timecourse_selections(self.r.selections) - - # FIXME: set integrator settings - - def set_timecourse_selections(self, selections: Iterator[str]): - """Set the timecourse selections.""" - for simulator in self.simulators: - simulator.set_timecourse_selections.remote(selections) - - def _timecourses(self, simulations: List[TimecourseSim]) -> List[pd.DataFrame]: - """Run all simulations with given model and collect the results. - - :param simulations: List[TimecourseSim] - :return: Result - """ - # Strip units for parallel simulations - for sim in simulations: - sim.strip_units() - - # Split simulations in chunks for actors - # !simulation have to stay in same order to reconstruct dimensions! - chunk_indices = np.array_split(np.arange(len(simulations)), self.actor_count) - chunks = [[] for _ in range(self.actor_count)] - for k, indices in enumerate(chunk_indices): - for index in indices: - chunks[k].append(simulations[index]) - - tc_ids = [] - for k, simulator in enumerate(self.simulators): - tcs_id = simulator.work.remote(chunks[k]) - tc_ids.append(tcs_id) - - results = ray.get(tc_ids) - # flatten list of lists [[df, df], [df, df], ...] - # indices = [k for sublist in chunks_indices for k in sublist] - return [df for sublist in results for df in sublist] - - @staticmethod - def _create_chunks(item, size: int): - """Yield successive sized chunks from item.""" - for i in range(0, len(item), size): - yield item[i : i + size] diff --git a/src/sbmlsim/simulator/simulation_serial.py b/src/sbmlsim/simulator/simulation_serial.py index b528ef81..2096b506 100644 --- a/src/sbmlsim/simulator/simulation_serial.py +++ b/src/sbmlsim/simulator/simulation_serial.py @@ -1,22 +1,21 @@ """Serial simulator.""" -from typing import List - +from typing import List, Optional, Union +from pathlib import Path +import roadrunner import pandas as pd from pint import Quantity -from roadrunner import roadrunner from sbmlutils import log -from sbmlsim.model import AbstractModel, RoadrunnerSBMLModel +from sbmlsim.model import AbstractModel, RoadrunnerSBMLModel, ModelChange from sbmlsim.result import XResult -from sbmlsim.simulation import ScanSim, TimecourseSim -from sbmlsim.simulator.simulation import SimulatorWorker +from sbmlsim.simulation import ScanSim, Timecourse, TimecourseSim from sbmlsim.units import UnitsInformation logger = log.get_logger(__name__) -class SimulatorSerial(SimulatorWorker): +class SimulatorSerial: """Serial simulator using a single core. A single simulator can run many different models. @@ -24,59 +23,52 @@ class SimulatorSerial(SimulatorWorker): cores. """ - def __init__(self, model=None, **kwargs): + def __init__(self, model: Union[str|Path|RoadrunnerSBMLModel|AbstractModel] = None, **kwargs): """Initialize serial simulator. :param model: Path to model or model :param kwargs: integrator settings """ - # default settings + self.r: Optional[roadrunner.RoadRunner] = None + self.model: Optional[RoadrunnerSBMLModel] = None + + # integrator settings self.integrator_settings = { "absolute_tolerance": 1e-10, "relative_tolerance": 1e-10, + **kwargs } - self.integrator_settings.update(kwargs) - self.set_model(model) - # self.model: Optional[AbstractModel, RoadrunnerSBMLModel] = None - def set_model(self, model): - """Set model for simulator and updates the integrator settings. - """ - logger.debug("SimulatorSerial.set_model") - if model is None: - self.model = None - else: - # if isinstance(model, AbstractModel): - # self.model = model - # else: - # handle path, urn, ... - - # FIXME: this is probably the issue - self.model = RoadrunnerSBMLModel( - source=model, settings=self.integrator_settings - ) + # set model + self.set_model(model) + def set_model(self, model: Union[str|Path|RoadrunnerSBMLModel|AbstractModel]): + """Set model for simulator and updates the integrator settings.""" + logger.info("SimulatorSerial.set_model") + self.model = None + if model is not None: + if isinstance(model, RoadrunnerSBMLModel): + self.model = model + elif isinstance(model, AbstractModel): + self.model = RoadrunnerSBMLModel.from_abstract_model( + abstract_model=model + ) + elif isinstance(model, (str, Path)): + self.model = RoadrunnerSBMLModel( + source=model, + ) + + self.r = model.r self.set_integrator_settings(**self.integrator_settings) def set_integrator_settings(self, **kwargs): """Set settings in the integrator.""" - if isinstance(self.model, RoadrunnerSBMLModel): - RoadrunnerSBMLModel.set_integrator_settings(self.model.r, **kwargs) - else: - logger.warning( - "Integrator settings can only be set on RoadrunnerSBMLModel." - ) + RoadrunnerSBMLModel.set_integrator_settings(self.r, **kwargs) def set_timecourse_selections(self, selections): """Set timecourse selection in model.""" - logger.debug(f"set_timecourse_selections: {selections}") RoadrunnerSBMLModel.set_timecourse_selections(self.r, selections=selections) - @property - def r(self) -> roadrunner.ExecutableModel: - """Get the RoadRunner model.""" - return self.model._model - @property def uinfo(self) -> UnitsInformation: """Get model unit information.""" @@ -111,9 +103,115 @@ def run_scan(self, scan: ScanSim) -> XResult: return XResult.from_dfs(dfs=dfs, scan=scan, uinfo=self.uinfo) def _timecourses(self, simulations: List[TimecourseSim]) -> List[pd.DataFrame]: - if len(simulations) > 1: - logger.warning( - "Use of SimulatorSerial to run multiple timecourses. " - "Use SimulatorParallel instead." - ) return [self._timecourse(sim) for sim in simulations] + + def _timecourse(self, simulation: TimecourseSim) -> pd.DataFrame: + """Timecourse simulation. + + Requires for all timecourse definitions in the timecourse simulation + to be unit normalized. The changes have no units any more + for parallel simulations. + You should never call this function directly! + + :param simulation: Simulation definition(s) + :return: DataFrame with results + """ + if isinstance(simulation, Timecourse): + simulation = TimecourseSim(timecourses=[simulation]) + + if simulation.reset: + self.r.resetToOrigin() + + frames = [] + t_offset = simulation.time_offset + for k, tc in enumerate(simulation.timecourses): + + if k == 0 and tc.model_changes: + # [1] apply model changes of first simulation + logger.debug("Applying model changes") + for key, item in tc.model_changes.items(): + if key.startswith("init"): + logger.error( + f"Initial model changes should be provided " + f"without 'init': '{key} = {item}'" + ) + # FIXME: implement model changes via init + # init_key = f"init({key})" + init_key = key + try: + value = item.magnitude + except AttributeError: + value = item + + try: + self.r[init_key] = value + except RuntimeError: + logger.error(f"roadrunner RuntimeError: '{init_key} = {item}'") + # boundary condition=true species, trying direct fallback + # see https://github.com/sys-bio/roadrunner/issues/711 + init_key = key + self.r[key] = value + + logger.debug(f"\t{init_key} = {item}") + + # [2] re-evaluate initial assignments + # https://github.com/sys-bio/roadrunner/issues/710 + # logger.debug("Reevaluate initial conditions") + # FIXME/TODO: support initial model changes + # self.r.resetAll() + # self.r.reset(SelectionRecord.DEPENDENT_FLOATING_AMOUNT) + # self.r.reset(SelectionRecord.DEPENDENT_INITIAL_GLOBAL_PARAMETER) + + # [3] apply model manipulations + # model manipulations are applied to model + if len(tc.model_manipulations) > 0: + # FIXME: update to support roadrunner model changes + for key, value in tc.model_changes.items(): + if key == ModelChange.CLAMP_SPECIES: + for sid, formula in value.items(): + ModelChange.clamp_species(self.r, sid, formula) + else: + raise ValueError( + f"Unsupported model change: " + f"'{key}': {value}. Supported changes are: " + f"['{ModelChange.CLAMP_SPECIES}']" + ) + + # [4] apply changes + if tc.changes: + logger.debug("Applying simulation changes") + for key, item in tc.changes.items(): + + # FIXME: handle concentrations/amounts/default + # TODO: Figure out the hasOnlySubstanceUnit flag! (roadrunner) + # r: roadrunner.ExecutableModel = self.r + + try: + self.r[key] = float(item.magnitude) + except AttributeError: + self.r[key] = float(item) + logger.debug(f"\t{key} = {item}") + + # debug model state + # FIXME: report issue + # sbml_str = self.r.getCurrentSBML() + # with open("/home/mkoenig/git/pkdb_models/pkdb_models/models/dextromethorphan/results/debug/test.xml", "w") as f_out: + # f_out.write(sbml_str) + + # run simulation + integrator = self.r.integrator + # FIXME: support simulation by times + if integrator.getValue("variable_step_size"): + s = self.r.simulate(start=tc.start, end=tc.end) + else: + s = self.r.simulate(start=tc.start, end=tc.end, steps=tc.steps) + + df = pd.DataFrame(s, columns=s.colnames) + df.time = df.time + t_offset + + if not tc.discard: + # discard timecourses (pre-simulation) + t_offset += tc.end + frames.append(df) + + return pd.concat(frames, sort=False) diff --git a/src/sbmlsim/test/test_simulator.py b/src/sbmlsim/test/test_simulator.py index e92cd8c7..dde0d721 100644 --- a/src/sbmlsim/test/test_simulator.py +++ b/src/sbmlsim/test/test_simulator.py @@ -3,7 +3,6 @@ from sbmlsim.model import RoadrunnerSBMLModel from sbmlsim.simulator import SimulatorSerial -from sbmlsim.simulator.simulation_ray import SimulatorParallel from sbmlsim.test import MODEL_REPRESSILATOR @@ -37,29 +36,3 @@ def test_tolerances_serial2(): ) assert isinstance(simulator, SimulatorSerial) _tolerance_test(r=simulator.model.r, abs_tol=abs_tol, rel_tol=rel_tol) - - -def test_tolerances_parallel(): - abs_tol = 1e-14 - rel_tol = 1e-14 - - simulator = SimulatorParallel( - model=MODEL_REPRESSILATOR, - absolute_tolerance=abs_tol, - relative_tolerance=rel_tol, - ) - assert isinstance(simulator, SimulatorParallel) - _tolerance_test(r=simulator.model.r, abs_tol=abs_tol, rel_tol=rel_tol) - - -def test_tolerances_parallel2(): - abs_tol = 1e-14 - rel_tol = 1e-14 - - simulator = SimulatorParallel( - model=RoadrunnerSBMLModel(source=MODEL_REPRESSILATOR), - absolute_tolerance=abs_tol, - relative_tolerance=rel_tol, - ) - assert isinstance(simulator, SimulatorParallel) - _tolerance_test(r=simulator.model.r, abs_tol=abs_tol, rel_tol=rel_tol)