Skip to content

Commit

Permalink
feat: separated task preprocessing from simulation
Browse files Browse the repository at this point in the history
@jonrkarr
jonrkarr committed Sep 17, 2021
1 parent 54d42b4 commit 0a9a181
Showing 5 changed files with 354 additions and 120 deletions.
2 changes: 1 addition & 1 deletion Dockerfile
View file
Original file line number Diff line number Diff line change
@@ -14,7 +14,7 @@
# Base OS
FROM python:3.9-slim-buster

ARG VERSION="0.1.30"
ARG VERSION="0.1.31"
ARG SIMULATOR_VERSION="4.34.251"

# metadata
2 changes: 1 addition & 1 deletion biosimulators_copasi/_version.py
View file
Original file line number Diff line number Diff line change
@@ -1 +1 @@
__version__ = '0.1.30'
__version__ = '0.1.31'
329 changes: 213 additions & 116 deletions biosimulators_copasi/core.py
View file
Original file line number Diff line number Diff line change
@@ -12,15 +12,15 @@
from biosimulators_utils.log.data_model import CombineArchiveLog, TaskLog, StandardOutputErrorCapturerLevel # noqa: F401
from biosimulators_utils.viz.data_model import VizFormat # noqa: F401
from biosimulators_utils.report.data_model import ReportFormat, VariableResults, SedDocumentResults # noqa: F401
from biosimulators_utils.sedml.data_model import (Task, ModelLanguage, UniformTimeCourseSimulation, # noqa: F401
from biosimulators_utils.sedml.data_model import (Task, ModelLanguage, ModelAttributeChange, UniformTimeCourseSimulation, # noqa: F401
Variable, Symbol)
from biosimulators_utils.sedml import validation
from biosimulators_utils.sedml.exec import exec_sed_doc as base_exec_sed_doc
from biosimulators_utils.simulator.utils import get_algorithm_substitution_policy
from biosimulators_utils.utils.core import raise_errors_warnings
from biosimulators_utils.warnings import warn, BioSimulatorsWarning
from kisao.data_model import AlgorithmSubstitutionPolicy, ALGORITHM_SUBSTITUTION_POLICY_LEVELS
from .data_model import KISAO_ALGORITHMS_MAP
from .data_model import KISAO_ALGORITHMS_MAP, Units
from .utils import (get_algorithm_id, set_algorithm_parameter_value,
get_copasi_model_object_by_sbml_id, get_copasi_model_obj_sbml_ids)
import COPASI
@@ -134,23 +134,171 @@ def exec_sed_task(task, variables, preprocessed_task=None, log=None, config=None
model = task.model
sim = task.simulation

# initialize COPASI task
copasi_model = preprocessed_task['model']['model']

# modify model
if model.changes:
raise_errors_warnings(validation.validate_model_change_types(model.changes, (ModelAttributeChange,)),
error_summary='Changes for model `{}` are not supported.'.format(model.id))
model_change_obj_map = preprocessed_task['model']['model_change_obj_map']
changed_objects = COPASI.ObjectStdVector()
for change in model.changes:
model_obj_set_func, ref = model_change_obj_map[change.target]
new_value = float(change.new_value)
model_obj_set_func(new_value)
changed_objects.push_back(ref)

copasi_model.compileIfNecessary()
copasi_model.updateInitialValues(changed_objects)

# initialize simulation
copasi_model.setInitialTime(sim.initial_time)
copasi_model.forceCompile()

copasi_task = preprocessed_task['task']
copasi_problem = copasi_task.getProblem()
copasi_problem.setOutputStartTime(sim.output_start_time)
copasi_problem.setDuration(sim.output_end_time - sim.initial_time)
if sim.output_end_time == sim.output_start_time:
if sim.output_start_time == sim.initial_time:
step_number = sim.number_of_points
else:
raise NotImplementedError('Output end time must be greater than the output start time.')
else:
step_number = (
sim.number_of_points
* (sim.output_end_time - sim.initial_time)
/ (sim.output_end_time - sim.output_start_time)
)
if step_number != math.floor(step_number):
raise NotImplementedError('Time course must specify an integer number of time points')
else:
step_number = int(step_number)
copasi_problem.setStepNumber(step_number)

# setup data handler
copasi_data_handler = COPASI.CDataHandler()
variable_common_name_map = preprocessed_task['model']['variable_common_name_map']
for variable in variables:
common_name = variable_common_name_map[(variable.target, variable.symbol)]
copasi_data_handler.addDuringName(common_name)
if not copasi_task.initializeRawWithOutputHandler(COPASI.CCopasiTask.OUTPUT_DURING, copasi_data_handler):
raise RuntimeError("Output handler could not be initialized:\n\n {}".format(
get_copasi_error_message(sim.algorithm.kisao_id).replace('\n', "\n ")))

# Execute simulation
result = copasi_task.processRaw(True)
warning_details = copasi_task.getProcessWarning()
if warning_details:
alg_kisao_id = preprocessed_task['simulation']['algorithm_kisao_id']
warn(get_copasi_error_message(alg_kisao_id, warning_details), BioSimulatorsWarning)
if not result:
alg_kisao_id = preprocessed_task['simulation']['algorithm_kisao_id']
error_details = copasi_task.getProcessError()
raise RuntimeError(get_copasi_error_message(alg_kisao_id, error_details))

# collect simulation predictions
number_of_recorded_points = copasi_data_handler.getNumRowsDuring()

if (
variables
and number_of_recorded_points != (sim.number_of_points + 1)
and (sim.output_end_time != sim.output_start_time or sim.output_start_time != sim.initial_time)
):
raise RuntimeError('Simulation produced {} rather than {} time points'.format(
number_of_recorded_points, sim.number_of_points)
) # pragma: no cover # unreachable because COPASI produces the correct number of outputs

variable_results = VariableResults()
for variable in variables:
variable_results[variable.id] = numpy.full((number_of_recorded_points,), numpy.nan)

for i_step in range(number_of_recorded_points):
step_values = copasi_data_handler.getNthRow(i_step)
for variable, value in zip(variables, step_values):
variable_results[variable.id][i_step] = value

if sim.output_end_time == sim.output_start_time and sim.output_start_time == sim.initial_time:
for variable in variables:
variable_results[variable.id] = numpy.concatenate((
variable_results[variable.id][0:1],
numpy.full((sim.number_of_points,), variable_results[variable.id][1]),
))

copasi_data_handler.cleanup()
copasi_data_handler.close()

# log action
if config.LOG:
log.algorithm = preprocessed_task['simulation']['algorithm_kisao_id']
log.simulator_details = {
'methodName': preprocessed_task['simulation']['method_name'],
'methodCode': preprocessed_task['simulation']['algorithm_copasi_id'],
'parameters': preprocessed_task['simulation']['method_parameters'],
}

# return results and log
return variable_results, log


def get_copasi_error_message(algorithm_kisao_id, details=None):
""" Get an error message from COPASI
Args:
algorithm_kisao_id (:obj:`str`): KiSAO id of algorithm of failed simulation
details (:obj:`str`, optional): details of of error
Returns:
:obj:`str`: COPASI error message
"""
error_msg = 'Simulation with algorithm {} ({}) failed'.format(
algorithm_kisao_id, KISAO_ALGORITHMS_MAP.get(algorithm_kisao_id, {}).get('name', 'N/A'))
if not details:
details = COPASI.CCopasiMessage.getLastMessage().getText()
if details:
details = '\n'.join(line[min(2, len(line)):] for line in details.split('\n')
if not (line.startswith('>') and line.endswith('<')))
error_msg += ':\n\n ' + details.replace('\n', '\n ')
return error_msg


def preprocess_sed_task(task, variables, config=None):
""" Preprocess a SED task, including its possible model changes and variables. This is useful for avoiding
repeatedly initializing tasks on repeated calls of :obj:`exec_sed_task`.
Args:
task (:obj:`Task`): task
variables (:obj:`list` of :obj:`Variable`): variables that should be recorded
config (:obj:`Config`, optional): BioSimulators common configuration
Returns:
:obj:`object`: preprocessed information about the task
"""
config = config or get_config()

model = task.model
sim = task.simulation

if config.VALIDATE_SEDML:
raise_errors_warnings(validation.validate_task(task),
error_summary='Task `{}` is invalid.'.format(task.id))
raise_errors_warnings(validation.validate_model_language(model.language, ModelLanguage.SBML),
error_summary='Language for model `{}` is not supported.'.format(model.id))
raise_errors_warnings(validation.validate_model_change_types(model.changes, ()),
raise_errors_warnings(validation.validate_model_change_types(model.changes, (ModelAttributeChange,)),
error_summary='Changes for model `{}` are not supported.'.format(model.id))
raise_errors_warnings(*validation.validate_model_changes(task.model),
raise_errors_warnings(*validation.validate_model_changes(model),
error_summary='Changes for model `{}` are invalid.'.format(model.id))
raise_errors_warnings(validation.validate_simulation_type(sim, (UniformTimeCourseSimulation, )),
error_summary='{} `{}` is not supported.'.format(sim.__class__.__name__, sim.id))
raise_errors_warnings(*validation.validate_simulation(sim),
error_summary='Simulation `{}` is invalid.'.format(sim.id))
raise_errors_warnings(*validation.validate_data_generator_variables(variables),
error_summary='Data generator variables for task `{}` are invalid.'.format(task.id))

model_etree = lxml.etree.parse(model.source)
target_x_paths_ids = validation.validate_target_xpaths(variables, model_etree, attr='id')
model_change_target_sbml_id_map = validation.validate_target_xpaths(model.changes, model_etree, attr='id')
variable_target_sbml_id_map = validation.validate_target_xpaths(variables, model_etree, attr='id')

if config.VALIDATE_SEDML_MODELS:
raise_errors_warnings(*validation.validate_model(model, [], working_dir='.'),
@@ -208,11 +356,52 @@ def exec_sed_task(task, variables, preprocessed_task=None, log=None, config=None
else:
raise

# validate model changes
model_change_obj_map = {}

invalid_changes = []

units = KISAO_ALGORITHMS_MAP[exec_alg_kisao_id]['default_units']
for change in model.changes:
target_sbml_id = model_change_target_sbml_id_map[change.target]
copasi_model_obj = get_copasi_model_object_by_sbml_id(copasi_model, target_sbml_id, units)
if copasi_model_obj is None:
invalid_changes.append(change.target)
else:
model_obj_parent = copasi_model_obj.getObjectParent()

if isinstance(model_obj_parent, COPASI.CCompartment):
set_func = model_obj_parent.setInitialValue
ref = model_obj_parent.getInitialValueReference()

elif isinstance(model_obj_parent, COPASI.CModelValue):
set_func = model_obj_parent.setInitialValue
ref = model_obj_parent.getInitialValueReference()

elif isinstance(model_obj_parent, COPASI.CMetab):
if units == Units.discrete:
set_func = model_obj_parent.setInitialValue
ref = model_obj_parent.getInitialValueReference()
else:
set_func = model_obj_parent.setInitialConcentration
ref = model_obj_parent.getInitialConcentrationReference()

model_change_obj_map[change.target] = (set_func, ref)

if invalid_changes:
raise ValueError(''.join([
'The following change targets are invalid:\n - {}\n\n'.format(
'\n - '.join(sorted(invalid_changes)),
),
'Targets must have one of the following ids:\n - {}'.format(
'\n - '.join(sorted(get_copasi_model_obj_sbml_ids(copasi_model))),
),
]))

# set up observables of the task
# due to a COPASI bug, :obj:`COPASI.CCopasiTask.initializeRawWithOutputHandler` must be called after
# :obj:`COPASI.CCopasiTask.setMethodType`

copasi_data_handler = COPASI.CDataHandler()
variable_common_name_map = {}

invalid_symbols = []
invalid_targets = []
@@ -226,7 +415,7 @@ def exec_sed_task(task, variables, preprocessed_task=None, log=None, config=None
invalid_symbols.append(variable.symbol)

else:
target_sbml_id = target_x_paths_ids[variable.target]
target_sbml_id = variable_target_sbml_id_map[variable.target]

copasi_model_obj = get_copasi_model_object_by_sbml_id(copasi_model, target_sbml_id,
KISAO_ALGORITHMS_MAP[exec_alg_kisao_id]['default_units'])
@@ -236,7 +425,7 @@ def exec_sed_task(task, variables, preprocessed_task=None, log=None, config=None
copasi_model_obj_common_name = copasi_model_obj.getCN().getString()

if copasi_model_obj_common_name is not None:
copasi_data_handler.addDuringName(COPASI.CRegisteredCommonName(copasi_model_obj_common_name))
variable_common_name_map[(variable.target, variable.symbol)] = COPASI.CRegisteredCommonName(copasi_model_obj_common_name)

if invalid_symbols:
raise NotImplementedError("".join([
@@ -256,118 +445,26 @@ def exec_sed_task(task, variables, preprocessed_task=None, log=None, config=None
),
]))

if not copasi_task.initializeRawWithOutputHandler(COPASI.CCopasiTask.OUTPUT_DURING, copasi_data_handler):
raise RuntimeError("Output handler could not be initialized:\n\n {}".format(
get_copasi_error_message(sim.algorithm.kisao_id).replace('\n', "\n ")))

# Execute simulation
copasi_task.setScheduled(True)

copasi_problem = copasi_task.getProblem()
copasi_model.setInitialTime(sim.initial_time)
copasi_model.forceCompile()
copasi_problem.setOutputStartTime(sim.output_start_time)
copasi_problem.setDuration(sim.output_end_time - sim.initial_time)
if sim.output_end_time == sim.output_start_time:
if sim.output_start_time == sim.initial_time:
step_number = sim.number_of_points
else:
raise NotImplementedError('Output end time must be greater than the output start time.')
else:
step_number = (
sim.number_of_points
* (sim.output_end_time - sim.initial_time)
/ (sim.output_end_time - sim.output_start_time)
)
if step_number != math.floor(step_number):
raise NotImplementedError('Time course must specify an integer number of time points')
else:
step_number = int(step_number)
copasi_problem.setStepNumber(step_number)
copasi_problem.setTimeSeriesRequested(False)
copasi_problem.setAutomaticStepSize(False)
copasi_problem.setOutputEvent(False)

result = copasi_task.processRaw(True)
warning_details = copasi_task.getProcessWarning()
if warning_details:
warn(get_copasi_error_message(exec_alg_kisao_id, warning_details), BioSimulatorsWarning)
if not result:
error_details = copasi_task.getProcessError()
raise RuntimeError(get_copasi_error_message(exec_alg_kisao_id, error_details))

# collect simulation predictions
number_of_recorded_points = copasi_data_handler.getNumRowsDuring()

if (
variables
and number_of_recorded_points != (sim.number_of_points + 1)
and (sim.output_end_time != sim.output_start_time or sim.output_start_time != sim.initial_time)
):
raise RuntimeError('Simulation produced {} rather than {} time points'.format(
number_of_recorded_points, sim.number_of_points)
) # pragma: no cover # unreachable because COPASI produces the correct number of outputs

variable_results = VariableResults()
for variable in variables:
variable_results[variable.id] = numpy.full((number_of_recorded_points,), numpy.nan)

for i_step in range(number_of_recorded_points):
step_values = copasi_data_handler.getNthRow(i_step)
for variable, value in zip(variables, step_values):
variable_results[variable.id][i_step] = value

if sim.output_end_time == sim.output_start_time and sim.output_start_time == sim.initial_time:
for variable in variables:
variable_results[variable.id] = numpy.concatenate((
variable_results[variable.id][0:1],
numpy.full((sim.number_of_points,), variable_results[variable.id][1]),
))

# log action
if config.LOG:
log.algorithm = exec_alg_kisao_id
log.simulator_details = {
'methodName': KISAO_ALGORITHMS_MAP[exec_alg_kisao_id]['id'],
'methodCode': alg_copasi_id,
'parameters': method_parameters,
}

# return results and log
return variable_results, log


def get_copasi_error_message(algorithm_kisao_id, details=None):
""" Get an error message from COPASI
Args:
algorithm_kisao_id (:obj:`str`): KiSAO id of algorithm of failed simulation
details (:obj:`str`, optional): details of of error
Returns:
:obj:`str`: COPASI error message
"""
error_msg = 'Simulation with algorithm {} ({}) failed'.format(
algorithm_kisao_id, KISAO_ALGORITHMS_MAP.get(algorithm_kisao_id, {}).get('name', 'N/A'))
if not details:
details = COPASI.CCopasiMessage.getLastMessage().getText()
if details:
details = '\n'.join(line[min(2, len(line)):] for line in details.split('\n')
if not (line.startswith('>') and line.endswith('<')))
error_msg += ':\n\n ' + details.replace('\n', '\n ')
return error_msg


def preprocess_sed_task(task, variables, config=None):
""" Preprocess a SED task, including its possible model changes and variables. This is useful for avoiding
repeatedly initializing tasks on repeated calls of :obj:`exec_sed_task`.
Args:
task (:obj:`Task`): task
variables (:obj:`list` of :obj:`Variable`): variables that should be recorded
config (:obj:`Config`, optional): BioSimulators common configuration
Returns:
:obj:`object`: preprocessed information about the task
"""
pass
# return preprocessed info
return {
'task': copasi_task,
'model': {
'model': copasi_model,
'model_change_obj_map': model_change_obj_map,
'variable_common_name_map': variable_common_name_map,
},
'simulation': {
'algorithm_kisao_id': exec_alg_kisao_id,
'algorithm_copasi_id': alg_copasi_id,
'method_name': KISAO_ALGORITHMS_MAP[exec_alg_kisao_id]['id'],
'method_parameters': method_parameters,
},
}
2 changes: 1 addition & 1 deletion requirements.txt
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
biosimulators_utils[logging,sbml] >= 0.1.116
biosimulators_utils[logging,sbml] >= 0.1.121
kisao
lxml
numpy
139 changes: 138 additions & 1 deletion tests/test_core_main.py
View file
Original file line number Diff line number Diff line change
@@ -8,7 +8,7 @@
"""

from biosimulators_copasi import __main__
from biosimulators_copasi.core import exec_sed_task, exec_sedml_docs_in_combine_archive
from biosimulators_copasi.core import exec_sed_task, exec_sedml_docs_in_combine_archive, preprocess_sed_task
from biosimulators_copasi.data_model import KISAO_ALGORITHMS_MAP
from biosimulators_utils.archive.io import ArchiveReader
from biosimulators_utils.combine import data_model as combine_data_model
@@ -341,6 +341,143 @@ def test_hybrid_rk45_partitioning(self):
)
self.assertFalse(numpy.any(numpy.isnan(results['Mdm2'])))

def test_exec_sed_task_with_changes(self):
task = sedml_data_model.Task(
model=sedml_data_model.Model(
source=os.path.join(os.path.dirname(__file__), 'fixtures', 'model.xml'),
language=sedml_data_model.ModelLanguage.SBML.value,
changes=[],
),
simulation=sedml_data_model.UniformTimeCourseSimulation(
algorithm=sedml_data_model.Algorithm(
kisao_id='KISAO_0000560',
),
initial_time=0.,
output_start_time=0.,
output_end_time=10.,
number_of_points=10,
),
)
model = task.model
sim = task.simulation

variable_ids = ['EmptySet', 'A', 'C', 'DA', 'DAp', "DR", "DRp", "MA", "MR", "R"]
variables = []
for variable_id in variable_ids:
model.changes.append(sedml_data_model.ModelAttributeChange(
id=variable_id,
target="/sbml:sbml/sbml:model/sbml:listOfSpecies/sbml:species[@id='{}']/@initialConcentration".format(variable_id),
target_namespaces=self.NAMESPACES_L2V4,
new_value=None,
))
variables.append(sedml_data_model.Variable(
id=variable_id,
target="/sbml:sbml/sbml:model/sbml:listOfSpecies/sbml:species[@id='{}']".format(variable_id),
target_namespaces=self.NAMESPACES_L2V4,
task=task,
))
preprocessed_task = preprocess_sed_task(task, variables)

model.changes = []
results, _ = exec_sed_task(task, variables, preprocessed_task=preprocessed_task)
with self.assertRaises(AssertionError):
for variable_id in variable_ids:
self.assertEqual(
results[variable_id][0:int(sim.number_of_points / 2 + 1)].shape,
results[variable_id][-int(sim.number_of_points / 2 + 1):].shape,
)
numpy.testing.assert_allclose(
results[variable_id][0:int(sim.number_of_points / 2 + 1)],
results[variable_id][-int(sim.number_of_points / 2 + 1):])

sim.output_end_time = sim.output_end_time / 2
sim.number_of_points = int(sim.number_of_points / 2)
results2, _ = exec_sed_task(task, variables, preprocessed_task=preprocessed_task)

for variable_id in variable_ids:
numpy.testing.assert_allclose(results2[variable_id], results[variable_id][0:sim.number_of_points + 1])

for variable_id in variable_ids:
model.changes.append(sedml_data_model.ModelAttributeChange(
id=variable_id,
target="/sbml:sbml/sbml:model/sbml:listOfSpecies/sbml:species[@id='{}']/@initialConcentration".format(variable_id),
target_namespaces=self.NAMESPACES_L2V4,
new_value=results2[variable_id][-1],
))

results3, _ = exec_sed_task(task, variables, preprocessed_task=preprocessed_task)

with self.assertRaises(AssertionError):
for variable_id in variable_ids:
numpy.testing.assert_allclose(results3[variable_id], results2[variable_id])

for variable_id in variable_ids:
numpy.testing.assert_allclose(results3[variable_id], results[variable_id][-(sim.number_of_points + 1):], rtol=5e-6)

task.model.changes = [
sedml_data_model.ModelAttributeChange(
id="model_change",
target="/sbml:sbml",
target_namespaces=self.NAMESPACES_L2V4,
new_value=None,
),
]
with self.assertRaises(ValueError):
preprocess_sed_task(task, variables)

def test_exec_sed_task_with_changes_discrete(self):
task = sedml_data_model.Task(
model=sedml_data_model.Model(
source=os.path.join(os.path.dirname(__file__), 'fixtures', 'BIOMD0000000806.xml'),
language=sedml_data_model.ModelLanguage.SBML.value,
changes=[],
),
simulation=sedml_data_model.UniformTimeCourseSimulation(
algorithm=sedml_data_model.Algorithm(
kisao_id='KISAO_0000027',
),
initial_time=0.,
output_start_time=0.,
output_end_time=10.,
number_of_points=10,
),
)
model = task.model
sim = task.simulation

model.changes = [
sedml_data_model.ModelAttributeChange(
id='UnInfected_Tumour_Cells_Xu',
target="/sbml:sbml/sbml:model/sbml:listOfSpecies/sbml:species[@id='{}']/@initialConcentration".format(
'UnInfected_Tumour_Cells_Xu'),
target_namespaces=self.NAMESPACES_L3V1,
new_value=None,
),
sedml_data_model.ModelAttributeChange(
id='r',
target="/sbml:sbml/sbml:model/sbml:listOfParameters/sbml:parameter[@id='{}']/@initialConcentration".format('r'),
target_namespaces=self.NAMESPACES_L3V1,
new_value=None,
),
sedml_data_model.ModelAttributeChange(
id='compartment',
target="/sbml:sbml/sbml:model/sbml:listOfCompartments/sbml:compartment[@id='{}']/@size".format('compartment'),
target_namespaces=self.NAMESPACES_L3V1,
new_value=None,
),
]

variables = [
sedml_data_model.Variable(
id='UnInfected_Tumour_Cells_Xu',
target="/sbml:sbml/sbml:model/sbml:listOfSpecies/sbml:species[@id='{}']".format('UnInfected_Tumour_Cells_Xu'),
target_namespaces=self.NAMESPACES_L3V1,
task=task,
),
]

preprocessed_task = preprocess_sed_task(task, variables)

def test_exec_sed_task_error_handling(self):
task = sedml_data_model.Task(
model=sedml_data_model.Model(

0 comments on commit 0a9a181

Please sign in to comment.