1''' Methods for executing COMBINE archives and SED-ML files involving simulations of models encoded in Smoldyn 2 3 Information about encoding Smoldyn simulations into SED-ML files and COMBINE archives is available at 4 `https://github.com/ssandrews/Smoldyn/blob/master/Using-Smoldyn-with-SED-ML-COMBINE-BioSimulators.md <https://github.com/ssandrews/Smoldyn/blob/master/Using-Smoldyn-with-SED-ML-COMBINE-BioSimulators.md>`_. 5''' 6 7__author__ = 'Jonathan Karr' 8__email__ = 'karr@mssm.edu' 9 10__all__ = [ 11 'exec_sedml_docs_in_combine_archive', 12 'exec_sed_task', 13] 14 15from .data_model import (SmoldynCommand, SmoldynOutputFile, SimulationChange, SimulationChangeExecution, AlgorithmParameterType, 16 KISAO_ALGORITHMS_MAP, KISAO_ALGORITHM_PARAMETERS_MAP) 17from biosimulators_utils.combine.exec import exec_sedml_docs_in_archive 18from biosimulators_utils.config import get_config, Config # noqa: F401 19from biosimulators_utils.log.data_model import CombineArchiveLog, TaskLog, StandardOutputErrorCapturerLevel # noqa: F401 20from biosimulators_utils.viz.data_model import VizFormat # noqa: F401 21from biosimulators_utils.report.data_model import ReportFormat, VariableResults, SedDocumentResults # noqa: F401 22from biosimulators_utils.sedml import validation 23from biosimulators_utils.sedml.data_model import (Task, ModelLanguage, ModelAttributeChange, # noqa: F401 24 UniformTimeCourseSimulation, AlgorithmParameterChange, Variable, 25 Symbol) 26from biosimulators_utils.sedml.exec import exec_sed_doc as base_exec_sed_doc 27from biosimulators_utils.utils.core import validate_str_value, parse_value, raise_errors_warnings 28from smoldyn import smoldyn 29import functools 30import os 31import numpy 32import pandas 33import re 34import tempfile 35import types # noqa: F401 36 37__all__ = ['exec_sedml_docs_in_combine_archive', 'exec_sed_task', 'exec_sed_doc', 'preprocess_sed_task'] 38 39 40def exec_sedml_docs_in_combine_archive(archive_filename, out_dir, config=None): 41 ''' Execute the SED tasks defined in a COMBINE/OMEX archive and save the outputs 42 43 Args: 44 archive_filename (:obj:`str`): path to COMBINE/OMEX archive 45 out_dir (:obj:`str`): path to store the outputs of the archive 46 47 * CSV: directory in which to save outputs to files 48 ``{ out_dir }/{ relative-path-to-SED-ML-file-within-archive }/{ report.id }.csv`` 49 * HDF5: directory in which to save a single HDF5 file (``{ out_dir }/reports.h5``), 50 with reports at keys ``{ relative-path-to-SED-ML-file-within-archive }/{ report.id }`` within the HDF5 file 51 52 config (:obj:`Config`, optional): BioSimulators common configuration 53 54 Returns: 55 :obj:`tuple`: 56 57 * :obj:`SedDocumentResults`: results 58 * :obj:`CombineArchiveLog`: log 59 ''' 60 return exec_sedml_docs_in_archive(exec_sed_doc, archive_filename, out_dir, 61 apply_xml_model_changes=False, 62 config=config) 63 64 65def exec_sed_doc(doc, working_dir, base_out_path, rel_out_path=None, 66 apply_xml_model_changes=True, 67 log=None, indent=0, pretty_print_modified_xml_models=False, 68 log_level=StandardOutputErrorCapturerLevel.c, config=None): 69 """ Execute the tasks specified in a SED document and generate the specified outputs 70 71 Args: 72 doc (:obj:`SedDocument` or :obj:`str`): SED document or a path to SED-ML file which defines a SED document 73 working_dir (:obj:`str`): working directory of the SED document (path relative to which models are located) 74 75 base_out_path (:obj:`str`): path to store the outputs 76 77 * CSV: directory in which to save outputs to files 78 ``{base_out_path}/{rel_out_path}/{report.id}.csv`` 79 * HDF5: directory in which to save a single HDF5 file (``{base_out_path}/reports.h5``), 80 with reports at keys ``{rel_out_path}/{report.id}`` within the HDF5 file 81 82 rel_out_path (:obj:`str`, optional): path relative to :obj:`base_out_path` to store the outputs 83 apply_xml_model_changes (:obj:`bool`, optional): if :obj:`True`, apply any model changes specified in the SED-ML file before 84 calling :obj:`task_executer`. 85 log (:obj:`SedDocumentLog`, optional): log of the document 86 indent (:obj:`int`, optional): degree to indent status messages 87 pretty_print_modified_xml_models (:obj:`bool`, optional): if :obj:`True`, pretty print modified XML models 88 log_level (:obj:`StandardOutputErrorCapturerLevel`, optional): level at which to log output 89 config (:obj:`Config`, optional): BioSimulators common configuration 90 simulator_config (:obj:`SimulatorConfig`, optional): tellurium configuration 91 92 Returns: 93 :obj:`tuple`: 94 95 * :obj:`ReportResults`: results of each report 96 * :obj:`SedDocumentLog`: log of the document 97 """ 98 return base_exec_sed_doc(exec_sed_task, doc, working_dir, base_out_path, 99 rel_out_path=rel_out_path, 100 apply_xml_model_changes=apply_xml_model_changes, 101 log=log, 102 indent=indent, 103 pretty_print_modified_xml_models=pretty_print_modified_xml_models, 104 log_level=log_level, 105 config=config) 106 107 108def exec_sed_task(task, variables, preprocessed_task=None, log=None, config=None): 109 ''' Execute a task and save its results 110 111 Args: 112 task (:obj:`Task`): task 113 variables (:obj:`list` of :obj:`Variable`): variables that should be recorded 114 preprocessed_task (:obj:`dict`, optional): preprocessed information about the task, including possible 115 model changes and variables. This can be used to avoid repeatedly executing the same initialization 116 for repeated calls to this method. 117 log (:obj:`TaskLog`, optional): log for the task 118 config (:obj:`Config`, optional): BioSimulators common configuration 119 120 Returns: 121 :obj:`tuple`: 122 123 :obj:`VariableResults`: results of variables 124 :obj:`TaskLog`: log 125 ''' 126 if not config: 127 config = get_config() 128 129 if config.LOG and not log: 130 log = TaskLog() 131 132 sed_model_changes = task.model.changes 133 sed_simulation = task.simulation 134 135 if preprocessed_task is None: 136 preprocessed_task = preprocess_sed_task(task, variables, config=config) 137 sed_model_changes = list(filter(lambda change: change.target in preprocessed_task['sed_smoldyn_simulation_change_map'], 138 sed_model_changes)) 139 140 # read Smoldyn configuration 141 smoldyn_simulation = preprocessed_task['simulation'] 142 143 # apply model changes to the Smoldyn configuration 144 sed_smoldyn_simulation_change_map = preprocessed_task['sed_smoldyn_simulation_change_map'] 145 for change in sed_model_changes: 146 smoldyn_change = sed_smoldyn_simulation_change_map.get(change.target, None) 147 if smoldyn_change is None or smoldyn_change.execution != SimulationChangeExecution.simulation: 148 raise NotImplementedError('Target `{}` can only be changed during simulation preprocessing.'.format(change.target)) 149 apply_change_to_smoldyn_simulation( 150 smoldyn_simulation, change, smoldyn_change) 151 152 # get the Smoldyn representation of the SED uniform time course simulation 153 smoldyn_simulation_run_timecourse_args = get_smoldyn_run_timecourse_args(sed_simulation) 154 155 # execute the simulation 156 smoldyn_run_args = dict( 157 **smoldyn_simulation_run_timecourse_args, 158 **preprocessed_task['simulation_run_alg_param_args'], 159 ) 160 smoldyn_simulation.run(**smoldyn_run_args, overwrite=True, display=False, quit_at_end=False) 161 162 # get the result of each SED variable 163 variable_output_cmd_map = preprocessed_task['variable_output_cmd_map'] 164 smoldyn_output_files = preprocessed_task['output_files'] 165 variable_results = get_variable_results(sed_simulation.number_of_steps, variables, variable_output_cmd_map, smoldyn_output_files) 166 167 # cleanup output files 168 for smoldyn_output_file in smoldyn_output_files.values(): 169 os.remove(smoldyn_output_file.filename) 170 171 # log simulation 172 if config.LOG: 173 log.algorithm = sed_simulation.algorithm.kisao_id 174 log.simulator_details = { 175 'class': 'smoldyn.Simulation', 176 'instanceAttributes': preprocessed_task['simulation_attrs'], 177 'method': 'run', 178 'methodArguments': smoldyn_run_args, 179 } 180 181 # return the values of the variables and log 182 return variable_results, log 183 184 185def preprocess_sed_task(task, variables, config=None): 186 """ Preprocess a SED task, including its possible model changes and variables. This is useful for avoiding 187 repeatedly initializing tasks on repeated calls of :obj:`exec_sed_task`. 188 189 Args: 190 task (:obj:`Task`): task 191 variables (:obj:`list` of :obj:`Variable`): variables that should be recorded 192 config (:obj:`Config`, optional): BioSimulators common configuration 193 194 Returns: 195 :obj:`dict`: preprocessed information about the task 196 """ 197 config = config or get_config() 198 199 sed_model = task.model 200 sed_simulation = task.simulation 201 202 if config.VALIDATE_SEDML: 203 raise_errors_warnings(validation.validate_task(task), 204 error_summary='Task `{}` is invalid.'.format(task.id)) 205 raise_errors_warnings(validation.validate_model_language(sed_model.language, ModelLanguage.Smoldyn), 206 error_summary='Language for model `{}` is not supported.'.format(sed_model.id)) 207 raise_errors_warnings(validation.validate_model_change_types(sed_model.changes, (ModelAttributeChange, )), 208 error_summary='Changes for model `{}` are not supported.'.format(sed_model.id)) 209 raise_errors_warnings(*validation.validate_model_changes(sed_model), 210 error_summary='Changes for model `{}` are invalid.'.format(sed_model.id)) 211 raise_errors_warnings(validation.validate_simulation_type(sed_simulation, (UniformTimeCourseSimulation, )), 212 error_summary='{} `{}` is not supported.'.format(sed_simulation.__class__.__name__, sed_simulation.id)) 213 raise_errors_warnings(*validation.validate_simulation(sed_simulation), 214 error_summary='Simulation `{}` is invalid.'.format(sed_simulation.id)) 215 raise_errors_warnings(*validation.validate_data_generator_variables(variables), 216 error_summary='Data generator variables for task `{}` are invalid.'.format(task.id)) 217 218 if sed_simulation.algorithm.kisao_id not in KISAO_ALGORITHMS_MAP: 219 msg = 'Algorithm `{}` is not supported. The following algorithms are supported:{}'.format( 220 sed_simulation.algorithm.kisao_id, 221 ''.join('\n {}: {}'.format(kisao_id, alg_props['name']) for kisao_id, alg_props in KISAO_ALGORITHMS_MAP.items()) 222 ) 223 raise NotImplementedError(msg) 224 225 # read Smoldyn configuration 226 simulation_configuration = read_smoldyn_simulation_configuration(sed_model.source) 227 normalize_smoldyn_simulation_configuration(simulation_configuration) 228 229 # turn off Smoldyn's graphics 230 disable_smoldyn_graphics_in_simulation_configuration(simulation_configuration) 231 232 # preprocess model changes 233 sed_smoldyn_preprocessing_change_map = {} 234 sed_smoldyn_simulation_change_map = {} 235 for change in sed_model.changes: 236 smoldyn_change = validate_model_change(change) 237 238 if smoldyn_change.execution == SimulationChangeExecution.preprocessing: 239 sed_smoldyn_preprocessing_change_map[change] = smoldyn_change 240 else: 241 sed_smoldyn_simulation_change_map[change.target] = smoldyn_change 242 243 # apply preprocessing-time changes 244 for change, smoldyn_change in sed_smoldyn_preprocessing_change_map.items(): 245 apply_change_to_smoldyn_simulation_configuration( 246 simulation_configuration, change, smoldyn_change) 247 248 # write the modified Smoldyn configuration to a temporary file 249 fid, smoldyn_configuration_filename = tempfile.mkstemp(suffix='.txt') 250 os.close(fid) 251 write_smoldyn_simulation_configuration(simulation_configuration, smoldyn_configuration_filename) 252 253 # initialize a simulation from the Smoldyn file 254 smoldyn_simulation = init_smoldyn_simulation_from_configuration_file(smoldyn_configuration_filename) 255 256 # clean up temporary file 257 os.remove(smoldyn_configuration_filename) 258 259 # apply the SED algorithm parameters to the Smoldyn simulation and to the arguments to its ``run`` method 260 smoldyn_simulation_attrs = {} 261 smoldyn_simulation_run_alg_param_args = {} 262 for sed_algorithm_parameter_change in sed_simulation.algorithm.changes: 263 val = get_smoldyn_instance_attr_or_run_algorithm_parameter_arg(sed_algorithm_parameter_change) 264 if val['type'] == AlgorithmParameterType.run_argument: 265 smoldyn_simulation_run_alg_param_args[val['name']] = val['value'] 266 else: 267 smoldyn_simulation_attrs[val['name']] = val['value'] 268 269 # apply the SED algorithm parameters to the Smoldyn simulation and to the arguments to its ``run`` method 270 for attr_name, value in smoldyn_simulation_attrs.items(): 271 setter = getattr(smoldyn_simulation, attr_name) 272 setter(value) 273 274 # validate SED variables 275 variable_output_cmd_map = validate_variables(variables) 276 277 # Setup Smoldyn output files for the SED variables 278 smoldyn_configuration_dirname = os.path.dirname(smoldyn_configuration_filename) 279 smoldyn_output_files = add_smoldyn_output_files_for_sed_variables( 280 smoldyn_configuration_dirname, variables, variable_output_cmd_map, smoldyn_simulation) 281 282 # return preprocessed information 283 return { 284 'simulation': smoldyn_simulation, 285 'simulation_attrs': smoldyn_simulation_attrs, 286 'simulation_run_alg_param_args': smoldyn_simulation_run_alg_param_args, 287 'sed_smoldyn_simulation_change_map': sed_smoldyn_simulation_change_map, 288 'variable_output_cmd_map': variable_output_cmd_map, 289 'output_files': smoldyn_output_files, 290 } 291 292 293def init_smoldyn_simulation_from_configuration_file(filename): 294 ''' Initialize a simulation for a Smoldyn model from a file 295 296 Args: 297 filename (:obj:`str`): path to model file 298 299 Returns: 300 :obj:`smoldyn.Simulation`: simulation 301 ''' 302 if not os.path.isfile(filename): 303 raise FileNotFoundError('Model source `{}` is not a file.'.format(filename)) 304 305 smoldyn_simulation = smoldyn.Simulation.fromFile(filename) 306 if not smoldyn_simulation.getSimPtr(): 307 error_code, error_msg = smoldyn.getError() 308 msg = 'Model source `{}` is not a valid Smoldyn file.\n\n {}: {}'.format( 309 filename, error_code.name[0].upper() + error_code.name[1:], error_msg.replace('\n', '\n ')) 310 raise ValueError(msg) 311 312 return smoldyn_simulation 313 314 315def read_smoldyn_simulation_configuration(filename): 316 ''' Read a configuration for a Smoldyn simulation 317 318 Args: 319 filename (:obj:`str`): path to model file 320 321 Returns: 322 :obj:`list` of :obj:`str`: simulation configuration 323 ''' 324 with open(filename, 'r') as file: 325 return [line.strip('\n') for line in file] 326 327 328def write_smoldyn_simulation_configuration(configuration, filename): 329 ''' Write a configuration for Smoldyn simulation to a file 330 331 Args: 332 configuration 333 filename (:obj:`str`): path to save configuration 334 ''' 335 with open(filename, 'w') as file: 336 for line in configuration: 337 file.write(line) 338 file.write('\n') 339 340 341def normalize_smoldyn_simulation_configuration(configuration): 342 ''' Normalize a configuration for a Smoldyn simulation 343 344 Args: 345 configuration (:obj:`list` of :obj:`str`): configuration for a Smoldyn simulation 346 ''' 347 # normalize spacing and comments 348 for i_line, line in enumerate(configuration): 349 if '#' in line: 350 cmd, comment = re.split('#+', line, maxsplit=1) 351 cmd = re.sub(' +', ' ', cmd).strip() 352 comment = comment.strip() 353 354 if cmd: 355 if comment: 356 line = cmd + ' # ' + comment 357 else: 358 line = cmd 359 else: 360 if comment: 361 line = '# ' + comment 362 else: 363 line = '' 364 365 else: 366 line = re.sub(' +', ' ', line).strip() 367 368 configuration[i_line] = line 369 370 # remove end_file and following lines 371 for i_line, line in enumerate(configuration): 372 if re.match(r'^end_file( |$)', line): 373 for i_line_remove in range(len(configuration) - i_line): 374 configuration.pop() 375 break 376 377 # remove empty starting lines 378 for line in list(configuration): 379 if not line: 380 configuration.pop(0) 381 else: 382 break 383 384 # remove empty ending lines 385 for line in reversed(configuration): 386 if not line: 387 configuration.pop() 388 else: 389 break 390 391 392def disable_smoldyn_graphics_in_simulation_configuration(configuration): 393 ''' Turn off graphics in the configuration of a Smoldyn simulation 394 395 Args: 396 configuration (:obj:`list` of :obj:`str`): simulation configuration 397 ''' 398 for i_line, line in enumerate(configuration): 399 if line.startswith('graphics '): 400 configuration[i_line] = re.sub(r'^graphics +[a-z_]+', 'graphics none', line) 401 402 403def validate_model_change(sed_model_change): 404 ''' Validate a SED model attribute change to a configuration for a Smoldyn simulation 405 406 ==================================================================== =================== 407 target newValue 408 ==================================================================== =================== 409 ``define {name}`` float 410 ``difc {species}`` float 411 ``difc {species}({state})`` float 412 ``difc_rule {species}({state})`` float 413 ``difm {species}`` float[] 414 ``difm {species}({state})`` float[] 415 ``difm_rule {species}({state})`` float[] 416 ``drift {species}`` float[] 417 ``drift {species}({state})`` float[] 418 ``drift_rule {species}({state})`` float[] 419 ``surface_drift {species}({state}) {surface} {panel-shape}`` float[] 420 ``surface_drift_rule {species}({state}) {surface} {panel-shape}`` float[] 421 ``killmol {species}({state})`` 0 422 ``killmolprob {species}({state}) {prob}`` 0 423 ``killmolincmpt {species}({state}) {compartment}`` 0 424 ``killmolinsphere {species}({state}) {surface}`` 0 425 ``killmoloutsidesystem {species}({state})`` 0 426 ``fixmolcount {species}({state})`` integer 427 ``fixmolcountincmpt {species}({staet}) {compartment}`` integer 428 ``fixmolcountonsurf {species}({state}) {surface}`` integer 429 ==================================================================== =================== 430 431 Args: 432 sed_model_change (:obj:`ModelAttributeChange`): SED model change 433 434 Returns: 435 :obj:`SimulationChange`: Smoldyn representation of the model change 436 437 Raises: 438 :obj:`NotImplementedError`: unsupported model change 439 ''' 440 # TODO: support additional types of model changes 441 442 target_type, _, target = sed_model_change.target.strip().partition(' ') 443 target_type = target_type.strip() 444 target = re.sub(r' +', ' ', target).strip() 445 446 if target_type in [ 447 'killmol', 'killmolprob', 'killmolinsphere', 'killmolincmpt', 'killmoloutsidesystem', 448 ]: 449 # Examples: 450 # killmol red 451 def new_line_func(new_value): 452 return target_type + ' ' + target 453 execution = SimulationChangeExecution.simulation 454 455 elif target_type in [ 456 'fixmolcount', 'fixmolcountonsurf', 'fixmolcountincmpt', 457 ]: 458 # Examples: 459 # fixmolcount red 460 species_name, species_target_sep, target = target.partition(' ') 461 462 def new_line_func(new_value): 463 return target_type + ' ' + species_name + ' ' + new_value + species_target_sep + target 464 465 execution = SimulationChangeExecution.simulation 466 467 elif target_type in [ 468 'define', 'difc', 'difc_rule', 'difm', 'difm_rule', 469 'drift', 'drift_rule', 'surface_drift', 'surface_drift_rule', 470 ]: 471 # Examples: 472 # define K_FWD 0.001 473 # difc S 3 474 # difm red 1 0 0 0 0 0 0 0 2 475 def new_line_func(new_value): 476 return target_type + ' ' + target + ' ' + new_value 477 478 execution = SimulationChangeExecution.preprocessing 479 480 # elif target_type in [ 481 # 'mol', 'compartment_mol', 'surface_mol', 482 # ]: 483 # # Examples: 484 # # mol 5 red u 485 # # compartment_mol 500 S inside 486 # # surface_mol 100 E(front) membrane all all 487 # def new_line_func(new_value): 488 # return target_type + ' ' + new_value + ' ' + target 489 # 490 # execution = SimulationChangeExecution.preprocessing 491 492 # elif target_type in [ 493 # 'dim', 494 # ]: 495 # # Examples: 496 # # dim 1 497 # def new_line_func(new_value): 498 # return target_type + ' ' + new_value 499 # 500 # execution = SimulationChangeExecution.preprocessing 501 502 # elif target_type in [ 503 # 'boundaries', 'low_wall', 'high_wall', 504 # ]: 505 # # Examples: 506 # # low_wall x -10 507 # # high_wall y 10 508 # def new_line_func(new_value): 509 # return target_type + ' ' + target + ' ' + new_value 510 # 511 # execution = SimulationChangeExecution.preprocessing 512 513 else: 514 raise NotImplementedError('Target `{}` is not supported.'.format(sed_model_change.target)) 515 516 return SimulationChange(new_line_func, execution) 517 518 519def apply_change_to_smoldyn_simulation(smoldyn_simulation, sed_change, smoldyn_change): 520 ''' Apply a SED model attribute change to a Smoldyn simulation 521 522 Args: 523 smoldyn_simulation (:obj:`smoldyn.Simulation`): Smoldyn simulation 524 sed_change (:obj:`ModelAttributeChange`): SED model change 525 smoldyn_change (:obj:`SimulationChange`): Smoldyn representation of the model change 526 ''' 527 new_value = str(sed_change.new_value).strip() 528 new_line = smoldyn_change.command(new_value) 529 smoldyn_simulation.addCommand(new_line, 'b') 530 531 532def apply_change_to_smoldyn_simulation_configuration(smoldyn_simulation_configuration, sed_change, smoldyn_change): 533 ''' Apply a SED model attribute change to a configuration for a Smoldyn simulation 534 535 Args: 536 smoldyn_simulation_configuration (:obj:`list` of :obj:`str`): configuration for the Smoldyn simulation 537 sed_change (:obj:`ModelAttributeChange`): SED model change 538 smoldyn_change (:obj:`SimulationChange`): Smoldyn representation of the model change 539 ''' 540 new_value = str(sed_change.new_value).strip() 541 new_line = smoldyn_change.command(new_value) 542 smoldyn_simulation_configuration.insert(0, new_line) 543 544 545def get_smoldyn_run_timecourse_args(sed_simulation): 546 ''' Get the Smoldyn representation of a SED uniform time course simulation 547 548 Args: 549 sed_simulation (:obj:`UniformTimeCourseSimulation`): SED uniform time course simulation 550 551 Returns: 552 :obj:`dict`: dictionary with keys ``start``, ``stop``, and ``dt`` that captures the Smoldyn 553 representation of the time course 554 555 Raises: 556 :obj:`NotImplementedError`: unsupported timecourse 557 ''' 558 number_of_steps = ( 559 ( 560 sed_simulation.output_end_time - sed_simulation.initial_time 561 ) / ( 562 sed_simulation.output_end_time - sed_simulation.output_start_time 563 ) * ( 564 sed_simulation.number_of_steps 565 ) 566 ) 567 if (number_of_steps - int(number_of_steps)) > 1e-8: 568 msg = ( 569 'Simulations must specify an integer number of steps, not {}.' 570 '\n Initial time: {}' 571 '\n Output start time: {}' 572 '\n Output end time: {}' 573 '\n Number of steps (output start to end time): {}' 574 ).format( 575 number_of_steps, 576 sed_simulation.initial_time, 577 sed_simulation.output_start_time, 578 sed_simulation.output_end_time, 579 sed_simulation.number_of_steps, 580 ) 581 raise NotImplementedError(msg) 582 583 dt = (sed_simulation.output_end_time - sed_simulation.output_start_time) / sed_simulation.number_of_steps 584 585 return { 586 'start': sed_simulation.initial_time, 587 'stop': sed_simulation.output_end_time, 588 'dt': dt, 589 } 590 591 592def get_smoldyn_instance_attr_or_run_algorithm_parameter_arg(sed_algorithm_parameter_change): 593 ''' Get the Smoldyn representation of a SED uniform time course simulation 594 595 Args: 596 sed_algorithm_parameter_change (:obj:`AlgorithmParameterChange`): SED change to a parameter of an algorithm 597 598 Returns: 599 :obj:`dict`: dictionary with keys ``type``, ``name``, and ``value`` that captures the Smoldyn representation 600 of the algorithm parameter 601 602 Raises: 603 :obj:`ValueError`: unsupported algorithm parameter value 604 :obj:`NotImplementedError`: unsupported algorithm parameter 605 ''' 606 parameter_props = KISAO_ALGORITHM_PARAMETERS_MAP.get(sed_algorithm_parameter_change.kisao_id, None) 607 if parameter_props: 608 if not validate_str_value(sed_algorithm_parameter_change.new_value, parameter_props['data_type']): 609 msg = '{} ({}) must be a {}, not `{}`.'.format( 610 parameter_props['name'], sed_algorithm_parameter_change.kisao_id, 611 parameter_props['data_type'].name, sed_algorithm_parameter_change.new_value, 612 ) 613 raise ValueError(msg) 614 new_value = parse_value(sed_algorithm_parameter_change.new_value, parameter_props['data_type']) 615 616 return { 617 'type': parameter_props['type'], 618 'name': parameter_props['name'], 619 'value': new_value, 620 } 621 622 else: 623 msg = 'Algorithm parameter `{}` is not supported. The following parameters are supported:{}'.format( 624 sed_algorithm_parameter_change.kisao_id, 625 ''.join('\n {}: {}'.format(kisao_id, parameter_props['name']) 626 for kisao_id, parameter_props in KISAO_ALGORITHM_PARAMETERS_MAP.items()) 627 ) 628 raise NotImplementedError(msg) 629 630 631def add_smoldyn_output_file(configuration_dirname, smoldyn_simulation): 632 ''' Add an output file to a Smoldyn simulation 633 634 Args: 635 configuration_dirname (:obj:`str`): path to the parent directory of the Smoldyn configuration file for the simulation 636 smoldyn_simulation (:obj:`smoldyn.Simulation`): Smoldyn simulation 637 638 Returns: 639 :obj:`SmoldynOutputFile`: output file 640 ''' 641 fid, filename = tempfile.mkstemp(dir=configuration_dirname, suffix='.ssv') 642 os.close(fid) 643 name = os.path.relpath(filename, configuration_dirname) 644 smoldyn_simulation.setOutputFile(name, append=False) 645 smoldyn_simulation.setOutputPath('./') 646 return SmoldynOutputFile(name=name, filename=filename) 647 648 649def add_commands_to_smoldyn_output_file(simulation, output_file, commands): 650 ''' Add commands to a Smoldyn output file 651 652 Args: 653 smoldyn_simulation (:obj:`smoldyn.Simulation`): Smoldyn simulation 654 smoldyn_output_file (:obj:`SmoldynOutputFile`): Smoldyn output file 655 commands (:obj:`list` of :obj`SmoldynCommand`): Smoldyn commands 656 ''' 657 for command in commands: 658 simulation.addCommand(command.command + ' ' + output_file.name, command.type) 659 660 661def validate_variables(variables): 662 ''' Validate SED variables 663 664 ============================================================================================================================================= =========================================================================================================================================== =========================================== 665 Smoldyn output file SED variable target Shape 666 ============================================================================================================================================= =========================================================================================================================================== =========================================== 667 ``molcount`` ``molcount {species}`` (numberOfSteps + 1,) 668 ``molcountspecies {species}({state})`` ``molcountspecies {species}({state})`` (numberOfSteps + 1,) 669 ``molcountspecieslist {species}({state})+`` ``molcountspecies {species}({state})`` (numberOfSteps + 1,) 670 ``molcountinbox {low-x} {hi-x}`` ``molcountinbox {species} {low-x} {hi-x}`` (numberOfSteps + 1,) 671 ``molcountinbox {low-x} {hi-x} {low-y} {hi-y}`` ``molcountinbox {species} {low-x} {hi-x} {low-y} {hi-y}`` (numberOfSteps + 1,) 672 ``molcountinbox {low-x} {hi-x} {low-y} {hi-y} {low-z} {hi-z}`` ``molcountinbox {species} {low-x} {hi-x} {low-y} {hi-y} {low-z} {hi-z}`` (numberOfSteps + 1,) 673 ``molcountincmpt {compartment}`` ``molcountincmpt {species} {compartment}`` (numberOfSteps + 1,) 674 ``molcountincmpts {compartment}+`` ``molcountincmpt {species} {compartment}`` (numberOfSteps + 1,) 675 ``molcountincmpt2 {compartment} {state}`` ``molcountincmpt2 {species} {compartment} {state}`` (numberOfSteps + 1,) 676 ``molcountonsurf {surface}`` ``molcountonsurf {species} {surface}`` (numberOfSteps + 1,) 677 ``molcountspace {species}({state}) {axis} {low} {hi} {bins} 0`` ``molcountspace {species}({state}) {axis} {low} {hi} {bins}`` (numberOfSteps + 1, bins) 678 ``molcountspace {species}({state}) {axis} {low} {hi} {bins} {low} {hi} 0`` ``molcountspace {species}({state}) {axis} {low} {hi} {bins} {low} {hi}`` (numberOfSteps + 1, bins) 679 ``molcountspace {species}({state}) {axis} {low} {hi} {bins} {low} {hi} {low} {hi} 0`` ``molcountspace {species}({state}) {axis} {low} {hi} {bins} {low} {hi} {low} {hi}`` (numberOfSteps + 1, bins) 680 ``molcountspace2d {species}({state}) z {low-x} {hi-x} {bins-x} {low-y} {hi-y} {bins-y} 0`` ``molcountspace2d {species}({state}) z {low-x} {hi-x} {bins-x} {low-y} {hi-y} {bins-y}`` (numberOfSteps + 1, bins-x, bins-y) 681 ``molcountspace2d {species}({state}) {axis} {low-1} {hi-1} {bins-1} {low-2} {hi-2} {bins-2} {low-3} {hi-3} 0`` ``molcountspace2d {species}({state}) {axis} {low-1} {hi-1} {bins-1} {low-2} {hi-2} {bins-3} {low-3} {hi-3}`` (numberOfSteps + 1, bins-1, bins-2) 682 ``molcountspaceradial {species}({state}) {center-x} {radius} {bins} 0`` ``molcountspaceradial {species}({state}) {center-x} {radius} {bins}`` (numberOfSteps + 1, bins) 683 ``molcountspaceradial {species}({state}) {center-x} {center-y} {radius} {bins} 0`` ``molcountspaceradial {species}({state}) {center-x} {center-y} {radius} {bins}`` (numberOfSteps + 1, bins) 684 ``molcountspaceradial {species}({state}) {center-x} {center-y} {center-z} {radius} {bins} 0`` ``molcountspaceradial {species}({state}) {center-x} {center-y} {center-z} {radius} {bins}`` (numberOfSteps + 1, bins) 685 ``molcountspacepolarangle {species}({state}) {center-x} {center-y} {pole-x} {pole-y} {radius-min} {radius-max} {bins} 0`` ``molcountspacepolarangle {species}({state}) {center-x} {center-y} {pole-x} {pole-y} {radius-min} {radius-max} {bins}`` (numberOfSteps + 1, bins) 686 ``molcountspacepolarangle {species}({state}) {center-x} {center-y} {center-z} {pole-x} {pole-y} {pole-z} {radius-min} {radius-max} {bins} 0`` ``molcountspacepolarangle {species}({state}) {center-x} {center-y} {center-z} {pole-x} {pole-y} {pole-z} {radius-min} {radius-max} {bins}`` (numberOfSteps + 1, bins) 687 ``radialdistribution {species-1}({state-1}) {species-2}({state-2}) {radius} {bins} 0`` ``radialdistribution {species-1}({state-1}) {species-2}({state-2}) {radius} {bins}`` (numberOfSteps + 1, bins) 688 ``radialdistribution2 {species-1}({state-1}) {species-2}({state-2}) {low-x} {hi-x} {low-y} {hi-y} {low-z} {hi-z} {radius} {bins} 0`` ``radialdistribution2 {species-1}({state-1}) {species-2}({state-2}) {low-x} {hi-x} {low-y} {hi-y} {low-z} {hi-z} {radius} {bins}`` (numberOfSteps + 1, bins) 689 ============================================================================================================================================= ========================================================================================================================================== =========================================== 690 691 Args: 692 variables (:obj:`list` of :obj:`Variable`): variables that should be recorded 693 694 Returns: 695 :obj:`dict`: dictionary that maps variable targets and symbols to Smoldyn output commands 696 ''' 697 # TODO: support additional kinds of outputs 698 699 variable_output_cmd_map = {} 700 701 invalid_symbols = [] 702 invalid_targets = [] 703 704 for variable in variables: 705 if variable.symbol: 706 if variable.symbol == Symbol.time.value: 707 output_command_args = 'molcount' 708 include_header = True 709 shape = None 710 results_slicer = functools.partial(results_key_slicer, key='time') 711 712 else: 713 invalid_symbols.append('{}: {}'.format(variable.id, variable.symbol)) 714 output_command_args = None 715 include_header = None 716 717 else: 718 output_command, _, output_args = re.sub(r' +', ' ', variable.target).partition(' ') 719 720 if output_command in ['molcount', 'molcountinbox', 'molcountincmpt', 'molcountincmpt2', 'molcountonsurf']: 721 species_name, _, output_args = output_args.partition(' ') 722 output_command_args = output_command + ' ' + output_args 723 include_header = True 724 shape = None 725 results_slicer = functools.partial(results_key_slicer, key=species_name) 726 727 elif output_command in ['molcountspecies']: 728 output_command_args = output_command + ' ' + output_args 729 include_header = False 730 shape = None 731 results_slicer = results_array_slicer 732 733 elif output_command in ['molcountspace', 'molcountspaceradial', 734 'molcountspacepolarangle', 'radialdistribution', 'radialdistribution2']: 735 output_command_args = output_command + ' ' + output_args + ' 0' 736 include_header = False 737 shape = None 738 results_slicer = results_matrix_slicer 739 740 elif output_command in ['molcountspace2d']: 741 output_command_args = output_command + ' ' + output_args + ' 0' 742 include_header = False 743 output_args_list = output_args.split(' ') 744 if len(output_args_list) == 8: 745 shape = (int(output_args_list[-4]), int(output_args_list[-1])) 746 else: 747 shape = (int(output_args_list[-6]), int(output_args_list[-3])) 748 results_slicer = None 749 750 else: 751 invalid_targets.append('{}: {}'.format(variable.id, variable.target)) 752 output_command_args = None 753 754 if output_command_args is not None: 755 output_command_args = output_command_args.strip() 756 variable_output_cmd_map[(variable.target, variable.symbol)] = (output_command_args, include_header, shape, results_slicer) 757 758 if invalid_symbols: 759 msg = '{} symbols cannot be recorded:\n {}\n\nThe following symbols can be recorded:\n {}'.format( 760 len(invalid_symbols), 761 '\n '.join(sorted(invalid_symbols)), 762 '\n '.join(sorted([Symbol.time.value])), 763 ) 764 raise ValueError(msg) 765 766 if invalid_targets: 767 valid_target_output_commands = [ 768 'molcount', 769 770 'molcount', 'molcountinbox', 'molcountincmpt', 'molcountincmpt2', 'molcountonsurf', 771 772 'molcountspace', 'molcountspaceradial', 773 'molcountspacepolarangle', 'radialdistribution', 'radialdistribution2', 774 775 'molcountspace2d', 776 ] 777 778 msg = '{} targets cannot be recorded:\n {}\n\nTargets are supported for the following output commands:\n {}'.format( 779 len(invalid_targets), 780 '\n '.join(sorted(invalid_targets)), 781 '\n '.join(sorted(set(valid_target_output_commands))), 782 ) 783 raise NotImplementedError(msg) 784 785 return variable_output_cmd_map 786 787 788def results_key_slicer(results, key): 789 """ Get the results for a key from a set of results 790 791 Args: 792 results (:obj:`pandas.DataFrame`): set of results 793 794 Returns: 795 :obj:`pandas.DataFrame`: results for a key 796 """ 797 return results.get(key, None) 798 799 800def results_array_slicer(results): 801 """ Extract an array of results from a matrix of time and results 802 803 Args: 804 results (:obj:`pandas.DataFrame`): matrix of time and results 805 806 Returns: 807 :obj:`pandas.DataFrame`: results 808 """ 809 return results.iloc[:, 1] 810 811 812def results_matrix_slicer(results): 813 """ Extract a matrix array of results from a matrix of time and results 814 815 Args: 816 results (:obj:`pandas.DataFrame`): matrix of time and results 817 818 Returns: 819 :obj:`pandas.DataFrame`: results 820 """ 821 return results.iloc[:, 1:] 822 823 824def add_smoldyn_output_files_for_sed_variables(configuration_dirname, variables, variable_output_cmd_map, smoldyn_simulation): 825 ''' Add Smoldyn output files for capturing each SED variable 826 827 Args: 828 configuration_dirname (:obj:`str`): path to the parent directory of the Smoldyn configuration file for the simulation 829 variables (:obj:`list` of :obj:`Variable`): variables that should be recorded 830 variable_output_cmd_map (:obj:`dict`): dictionary that maps variable targets and symbols to Smoldyn output commands 831 smoldyn_simulation (:obj:`smoldyn.Simulation`): Smoldyn simulation 832 833 Returns: 834 :obj:`dict` of :obj:`str` => :obj:`SmoldynOutputFile`: Smoldyn output files 835 ''' 836 smoldyn_output_files = {} 837 838 output_cmds = set() 839 for variable in variables: 840 output_cmds.add(variable_output_cmd_map[(variable.target, variable.symbol)]) 841 842 for command, include_header, _, _ in output_cmds: 843 add_smoldyn_output_file_for_output(configuration_dirname, smoldyn_simulation, 844 command, include_header, 845 smoldyn_output_files) 846 # return output files 847 return smoldyn_output_files 848 849 850def add_smoldyn_output_file_for_output(configuration_dirname, smoldyn_simulation, 851 smoldyn_output_command, include_header, smoldyn_output_files): 852 ''' Add a Smoldyn output file for molecule counts 853 854 Args: 855 configuration_dirname (:obj:`str`): path to the parent directory of the Smoldyn configuration file for the simulation 856 smoldyn_simulation (:obj:`smoldyn.Simulation`): Smoldyn simulation 857 smoldyn_output_command (:obj:`str`): Smoldyn output command (e.g., ``molcount``) 858 include_header (:obj:`bool`): whether to include a header 859 smoldyn_output_files (:obj:`dict` of :obj:`str` => :obj:`SmoldynOutputFile`): Smoldyn output files 860 ''' 861 smoldyn_output_files[smoldyn_output_command] = add_smoldyn_output_file(configuration_dirname, smoldyn_simulation) 862 863 commands = [] 864 if include_header: 865 commands.append(SmoldynCommand('molcountheader', 'B')) 866 commands.append(SmoldynCommand(smoldyn_output_command, 'E')) 867 868 add_commands_to_smoldyn_output_file( 869 smoldyn_simulation, 870 smoldyn_output_files[smoldyn_output_command], 871 commands, 872 ) 873 874 875def get_variable_results(number_of_steps, variables, variable_output_cmd_map, smoldyn_output_files): 876 ''' Get the result of each SED variable 877 878 Args: 879 number_of_steps (:obj:`int`): number of steps 880 variables (:obj:`list` of :obj:`Variable`): variables that should be recorded 881 variable_output_cmd_map (:obj:`dict`): dictionary that maps variable targets and symbols to Smoldyn output commands 882 smoldyn_output_files (:obj:`dict` of :obj:`str` => :obj:`SmoldynOutputFile`): Smoldyn output files 883 884 Returns: 885 :obj:`VariableResults`: result of each SED variable 886 887 Raises: 888 :obj:`ValueError`: unsupported results 889 ''' 890 smoldyn_results = {} 891 892 missing_variables = [] 893 894 variable_results = VariableResults() 895 for variable in variables: 896 output_command_args, _, shape, results_slicer = variable_output_cmd_map[(variable.target, variable.symbol)] 897 variable_result = get_smoldyn_output(output_command_args, True, shape, smoldyn_output_files, smoldyn_results) 898 if results_slicer: 899 variable_result = results_slicer(variable_result) 900 901 if variable_result is None: 902 missing_variables.append('{}: {}: {}'.format(variable.id, 'target', variable.target)) 903 else: 904 if variable_result.ndim == 1: 905 variable_results[variable.id] = variable_result.to_numpy()[-(number_of_steps + 1):, ] 906 elif variable_result.ndim == 2: 907 variable_results[variable.id] = variable_result.to_numpy()[-(number_of_steps + 1):, :] 908 else: 909 variable_results[variable.id] = variable_result[-(number_of_steps + 1):, :, :] 910 911 if missing_variables: 912 msg = '{} variables could not be recorded:\n {}'.format( 913 len(missing_variables), 914 '\n '.join(missing_variables), 915 ) 916 raise ValueError(msg) 917 918 return variable_results 919 920 921def get_smoldyn_output(smoldyn_output_command, has_header, three_d_shape, smoldyn_output_files, smoldyn_results): 922 ''' Get the simulated count of each molecule 923 924 Args: 925 smoldyn_output_command (:obj:`str`): Smoldyn output command (e.g., ``molcount``) 926 has_header (:obj:`bool`): whether to include a header 927 three_d_shape (:obj:`tuple` of :obj:`int`): dimensions of the output 928 smoldyn_output_files (:obj:`dict` of :obj:`str` => :obj:`SmoldynOutputFile`): Smoldyn output files 929 smoldyn_results (:obj:`dict`) 930 931 Returns: 932 :obj:`pandas.DataFrame` or :obj:`numpy.ndarray`: results 933 ''' 934 smoldyn_output_command = smoldyn_output_command.strip() 935 if smoldyn_output_command not in smoldyn_results: 936 smoldyn_output_file = smoldyn_output_files[smoldyn_output_command] 937 if three_d_shape: 938 with open(smoldyn_output_file.filename, 'r') as file: 939 data_list = [] 940 941 i_line = 0 942 for line in file: 943 if i_line % (three_d_shape[1] + 1) == 0: 944 time_point_data = [] 945 else: 946 profile = [int(el) for el in line.split(' ')] 947 time_point_data.append(profile) 948 949 if i_line % (three_d_shape[1] + 1) == three_d_shape[1]: 950 data_list.append(time_point_data) 951 952 i_line += 1 953 954 smoldyn_results[smoldyn_output_command] = numpy.array(data_list).transpose((0, 2, 1)) 955 956 else: 957 smoldyn_results[smoldyn_output_command] = pandas.read_csv(smoldyn_output_file.filename, sep=' ') 958 959 return smoldyn_results[smoldyn_output_command] 960