From 9d5243f47f051256fe672d4e0ebd52a73559c973 Mon Sep 17 00:00:00 2001 From: Arthur Goldberg Date: Fri, 8 Jun 2018 09:19:41 -0400 Subject: [PATCH] add Observables to checkpoints and RunResults; dynamically create checkpoints to test RunResults --- examples/simulation_run.py | 2 +- .../multialgorithm/test_dynamic_components.py | 4 +-- .../test_multialgorithm_simulation.py | 17 +++++----- tests/multialgorithm/test_observables.py | 2 +- tests/multialgorithm/test_run_results.py | 31 ++++++++++--------- wc_sim/multialgorithm/dynamic_components.py | 18 ++++++++--- wc_sim/multialgorithm/make_models.py | 4 +-- .../multialgorithm_checkpointing.py | 12 +++++-- .../multialgorithm_simulation.py | 2 +- wc_sim/multialgorithm/observables.py | 11 ++----- wc_sim/multialgorithm/run_results.py | 24 +++++++++++--- 11 files changed, 77 insertions(+), 50 deletions(-) diff --git a/examples/simulation_run.py b/examples/simulation_run.py index 212ccc65..7b85f6e7 100644 --- a/examples/simulation_run.py +++ b/examples/simulation_run.py @@ -26,7 +26,7 @@ # view results print(run_results.get('populations')) -# run_results also contains 'aggregate_states', 'random_states', and 'metadata' +# run_results also contains 'observables', 'aggregate_states', 'random_states', and 'metadata' for component in RunResults.COMPONENTS: print(run_results.get(component)) diff --git a/tests/multialgorithm/test_dynamic_components.py b/tests/multialgorithm/test_dynamic_components.py index 7f7dc250..1cba35f4 100644 --- a/tests/multialgorithm/test_dynamic_components.py +++ b/tests/multialgorithm/test_dynamic_components.py @@ -99,7 +99,7 @@ def read_model(self, model_filename): self.model = Reader().run(model_filename, strict=False) multialgorithm_simulation = MultialgorithmSimulation(self.model, None) dynamic_compartments = multialgorithm_simulation.dynamic_compartments - self.dynamic_model = DynamicModel(self.model, dynamic_compartments) + self.dynamic_model = DynamicModel(self.model, multialgorithm_simulation.local_species_population, dynamic_compartments) # TODO(Arthur): move this proportional test to a utility & use it instead of assertAlmostEqual everywhere def almost_equal_test(self, a, b, frac_diff=1/100): @@ -197,7 +197,7 @@ def test_eval_dynamic_observables(self): model = MakeModels.make_test_model('no reactions') # make a DynamicModel - dyn_mdl = DynamicModel(model, {}) + dyn_mdl = DynamicModel(model, lsp, {}) non_dependent_dynamic_observables = [] for non_dependent_observable in non_dependent_observables: diff --git a/tests/multialgorithm/test_multialgorithm_simulation.py b/tests/multialgorithm/test_multialgorithm_simulation.py index 5519920a..0d4c76e7 100644 --- a/tests/multialgorithm/test_multialgorithm_simulation.py +++ b/tests/multialgorithm/test_multialgorithm_simulation.py @@ -30,6 +30,7 @@ from wc_sim.multialgorithm.config import core as config_core_multialgorithm from wc_sim.multialgorithm.multialgorithm_checkpointing import (MultialgorithmicCheckpointingSimObj, MultialgorithmCheckpoint) +from wc_sim.multialgorithm.run_results import RunResults from wc_sim.core.simulation_checkpoint_object import CheckpointSimulationObject @@ -227,16 +228,16 @@ def test_run_ssa_suite(self): expected_mean_copy_numbers={specie:2000}, delta=50) # species counts, and cell mass and volume steadily decline - previous_ckpt = None + prev_ckpt = None for time in MultialgorithmCheckpoint.list_checkpoints(self.results_dir): ckpt = MultialgorithmCheckpoint.get_checkpoint(self.results_dir, time=time) - if previous_ckpt: - previous_species_pops, previous_aggregate_state = previous_ckpt.state - species_pops, aggregate_state = ckpt.state - self.assertTrue(species_pops[specie] < previous_species_pops[specie]) - self.assertTrue(aggregate_state['cell mass'] < previous_aggregate_state['cell mass']) - self.assertTrue(aggregate_state['cell volume'] < previous_aggregate_state['cell volume']) - previous_ckpt = ckpt + if prev_ckpt: + prev_species_pops, prev_observables, prev_aggregate_state = RunResults.get_state_components(prev_ckpt.state) + species_pops, observables, aggregate_state = RunResults.get_state_components(ckpt.state) + self.assertTrue(species_pops[specie] < prev_species_pops[specie]) + self.assertTrue(aggregate_state['cell mass'] < prev_aggregate_state['cell mass']) + self.assertTrue(aggregate_state['cell volume'] < prev_aggregate_state['cell volume']) + prev_ckpt = ckpt self.perform_ssa_test_run('2 species, 1 reaction', run_time=1000, diff --git a/tests/multialgorithm/test_observables.py b/tests/multialgorithm/test_observables.py index 36917190..8d99b6fa 100644 --- a/tests/multialgorithm/test_observables.py +++ b/tests/multialgorithm/test_observables.py @@ -57,7 +57,7 @@ def setUp(self): self.pseudo_function = PseudoFunction('fun_1', self.tokens, [self.obs_1]) self.model = MakeModels.make_test_model('1 species, 1 reaction') - self.dyn_mdl = DynamicModel(self.model, {}) + self.dyn_mdl = DynamicModel(self.model, self.lsp, {}) def test_dynamic_observable(self): diff --git a/tests/multialgorithm/test_run_results.py b/tests/multialgorithm/test_run_results.py index 266f777f..900ee1ca 100644 --- a/tests/multialgorithm/test_run_results.py +++ b/tests/multialgorithm/test_run_results.py @@ -6,48 +6,49 @@ :License: MIT """ -import os -import getpass import unittest import shutil import tempfile -from argparse import Namespace -import warnings import pandas import numpy +from wc_lang.core import SpeciesType from wc_sim.multialgorithm.multialgorithm_errors import MultialgorithmError from wc_sim.multialgorithm.simulation import Simulation from wc_sim.multialgorithm.run_results import RunResults +from wc_sim.multialgorithm.make_models import MakeModels class TestRunResults(unittest.TestCase): def setUp(self): - # use stored checkpoints and metadata from simulation of 2_species_1_reaction model - self.RESULTS_DIR = os.path.join(os.path.dirname(__file__), 'fixtures', 'results_dir') - self.results_dir = tempfile.mkdtemp() - self.results_copy = os.path.join(self.results_dir, 'results_copy') - shutil.copytree(self.RESULTS_DIR, self.results_copy) + SpeciesType.objects.reset() + # create stored checkpoints and metadata + self.temp_dir = tempfile.mkdtemp() + # create and run simulation + model = MakeModels.make_test_model('2 species, 1 reaction') + simulation = Simulation(model) self.checkpoint_period = 10 self.max_time = 100 + _, self.results_dir = simulation.run(end_time=self.max_time, results_dir=self.temp_dir, + checkpoint_period=self.checkpoint_period) def tearDown(self): - shutil.rmtree(self.results_dir) + shutil.rmtree(self.temp_dir) def test_run_results(self): - run_results_1 = RunResults(self.results_copy) + run_results_1 = RunResults(self.results_dir) # after run_results file created - run_results_2 = RunResults(self.results_copy) + run_results_2 = RunResults(self.results_dir) for component in RunResults.COMPONENTS: component_data = run_results_1.get(component) self.assertTrue(run_results_1.get(component).equals(run_results_2.get(component))) - expected_times = pandas.Float64Index(numpy.linspace(0, self.max_time, 1 + self.max_time/self.checkpoint_period)) - for component in ['populations', 'aggregate_states', 'random_states']: + for component in ['populations', 'observables', 'aggregate_states', 'random_states']: component_data = run_results_1.get(component) + self.assertFalse(component_data.empty) self.assertTrue(component_data.index.equals(expected_times)) # total population is invariant @@ -61,6 +62,6 @@ def test_run_results(self): def test_run_results_errors(self): - run_results = RunResults(self.results_copy) + run_results = RunResults(self.results_dir) with self.assertRaisesRegexp(MultialgorithmError, "component '.*' is not an element of "): run_results.get('not_a_component') diff --git a/wc_sim/multialgorithm/dynamic_components.py b/wc_sim/multialgorithm/dynamic_components.py index a3c8352e..b4a580bc 100644 --- a/wc_sim/multialgorithm/dynamic_components.py +++ b/wc_sim/multialgorithm/dynamic_components.py @@ -123,6 +123,8 @@ class DynamicModel(object): dynamic_compartments (:obj: `dict`): map from compartment ID to `DynamicCompartment`; the simulation's `DynamicCompartment`s, one for each compartment in `model` cellular_dyn_compartments (:obj:`list`): list of the cellular compartments + species_population (:obj:`LocalSpeciesPopulation`): an object that represents + the populations of species in this `DynamicCompartment` dynamic_observables (:obj:`dict` of `DynamicObservable`): the simulation's dynamic observables, indexed by their ids dynamic_functions (:obj:`dict` of `DynamicFunction`): the simulation's dynamic functions, @@ -133,18 +135,18 @@ class DynamicModel(object): a constant water_in_model (:obj:`bool`): if set, the model represents water """ - def __init__(self, model, dynamic_compartments): + def __init__(self, model, species_population, dynamic_compartments): """ Prepare a `DynamicModel` for a discrete-event simulation Args: model (:obj:`Model`): the description of the whole-cell model in `wc_lang` + species_population (:obj:`LocalSpeciesPopulation`): an object that represents + the populations of species in this `DynamicCompartment` dynamic_compartments (:obj: `dict`): the simulation's `DynamicCompartment`s, one for each compartment in `model` """ self.dynamic_compartments = dynamic_compartments - self.dynamic_observables = {} - self.dynamic_functions = {} - self.dynamic_stop_conditions = {} + self.species_population = species_population # Classify compartments into extracellular and cellular; those which are not extracellular are cellular # Assumes at most one extracellular compartment @@ -169,6 +171,14 @@ def __init__(self, model, dynamic_compartments): self.fraction_dry_weight = utils.get_component_by_id(model.get_parameters(), 'fractionDryWeight').value + # create dynamic observables + self.dynamic_observables = {} + for observable in model.observables: + self.dynamic_observables[observable.id] = DynamicObservable(self, self.species_population, observable) + + self.dynamic_functions = {} + self.dynamic_stop_conditions = {} + def cell_mass(self): """ Compute the cell's mass diff --git a/wc_sim/multialgorithm/make_models.py b/wc_sim/multialgorithm/make_models.py index e7eb845a..204e5615 100644 --- a/wc_sim/multialgorithm/make_models.py +++ b/wc_sim/multialgorithm/make_models.py @@ -105,8 +105,8 @@ def add_test_submodel(model, model_type, submodel_num, submodel_type, init_vol, Concentration(species=spec, value=concentration, units=ConcentrationUnit.M.value) else: Concentration(species=spec, value=default_concentration, units=ConcentrationUnit.M.value) - # TODO: fix some problem with SpeciesCoefficient that prevents coefficients from being identical - species_coefficient = spec.species_coefficients.create(coefficient=1. + i) + # TODO: fix some problem with SpeciesCoefficient that prevents coefficient values from being identical + species_coefficient = spec.species_coefficients.create(coefficient=1.5) obs_plain = model.observables.create(id='obs_{}_{}'.format(submodel_num, i)) obs_plain.species.append(species_coefficient) obs_coeff = obs_plain.observable_coefficients.create(coefficient=2.) diff --git a/wc_sim/multialgorithm/multialgorithm_checkpointing.py b/wc_sim/multialgorithm/multialgorithm_checkpointing.py index f9f71981..e3624ae8 100644 --- a/wc_sim/multialgorithm/multialgorithm_checkpointing.py +++ b/wc_sim/multialgorithm/multialgorithm_checkpointing.py @@ -43,10 +43,16 @@ def get_checkpoint_state(self, time): """ Obtain a checkpoint of the biological state Returns: - :obj:`tuple` of (`dict`, `dict`): dictionaries with the species populations and the - cell's aggregate state, respectively + :obj:`dict` of `dict`: dictionaries containing the simulation's state: `population` has + the simulation's species populations, `aggregate_state` contains its aggregrate compartment + states, and `observables` contains all of its observables """ - return (self.local_species_population.read(time), self.dynamic_model.get_aggregate_state()) + state = { + 'population': self.local_species_population.read(time), + 'aggregate_state': self.dynamic_model.get_aggregate_state(), + 'observables': self.dynamic_model.eval_dynamic_observables(time) + } + return state def get_random_state(self): """ Obtain a checkpoint of the random state diff --git a/wc_sim/multialgorithm/multialgorithm_simulation.py b/wc_sim/multialgorithm/multialgorithm_simulation.py index 71e04749..bdb60264 100644 --- a/wc_sim/multialgorithm/multialgorithm_simulation.py +++ b/wc_sim/multialgorithm/multialgorithm_simulation.py @@ -161,7 +161,7 @@ def build_simulation(self): dynamic model """ self.partition_species() - self.dynamic_model = DynamicModel(self.model, self.dynamic_compartments) + self.dynamic_model = DynamicModel(self.model, self.local_species_population, self.dynamic_compartments) if 'results_dir' in self.args and self.args['results_dir']: self.checkpointing_sim_obj = self.create_multialgorithm_checkpointing( self.args['results_dir'], diff --git a/wc_sim/multialgorithm/observables.py b/wc_sim/multialgorithm/observables.py index bdeafd89..e6460b98 100644 --- a/wc_sim/multialgorithm/observables.py +++ b/wc_sim/multialgorithm/observables.py @@ -19,17 +19,12 @@ ''' Also: - Include in checkpoints, and checkpoint processing - expand - calc all observables in dynamic_model - get_checkpoint_state - RunResults(), COMPONENTS, convert_checkpoints, - - Make available to rate law calculations - prohibit observable names that conflict with functions (but could be relaxed with smarter parsing) + Make observalbes available to rate law calculations expression parser for DynamicFunction, generalized for RateLaws too + prohibit observable names that conflict with functions (but could be relaxed with smarter parsing) failing test for wc_lang that shows the str.replace() problem push wc_lang + memoize performance comparison Later: clean up memoize cache file ''' diff --git a/wc_sim/multialgorithm/run_results.py b/wc_sim/multialgorithm/run_results.py index 766990fc..2b2df851 100644 --- a/wc_sim/multialgorithm/run_results.py +++ b/wc_sim/multialgorithm/run_results.py @@ -28,6 +28,7 @@ class RunResults(object): COMPONENTS = { 'populations', # predicted populations of species at all checkpoint times 'aggregate_states', # predicted aggregate states of the cell over the simulation + 'observables', # predicted values of all observables over the simulation 'random_states', # states of the simulation's random number geerators over the simulation 'metadata', # the simulation's global metadata } @@ -51,9 +52,11 @@ def __init__(self, results_dir): else: # create the HDF file containing the run results - population_df, aggregate_states_df, random_states_s = self.convert_checkpoints() + population_df, observables_df, aggregate_states_df, random_states_s = self.convert_checkpoints() # populations population_df.to_hdf(self._hdf_file(), 'populations') + # observables + observables_df.to_hdf(self._hdf_file(), 'observables') # aggregate states aggregate_states_df.to_hdf(self._hdf_file(), 'aggregate_states') # random states @@ -106,21 +109,28 @@ def convert_metadata(self): simulation_metadata = SimulationMetadata.read_metadata(self.results_dir) return pandas.Series(as_dict(simulation_metadata)) + @staticmethod + def get_state_components(state): + return (state['population'], state['observables'], state['aggregate_state']) + def convert_checkpoints(self): """ Convert the data in saved checkpoints into pandas dataframes Returns: :obj:`tuple` of pandas objects: dataframes of the components of a simulation checkpoint history - population_df, aggregate_states_df, random_states_s + population_df, observables_df, aggregate_states_df, random_states_s """ # create pandas objects for species populations, aggregate states and simulation random states checkpoints = Checkpoint.list_checkpoints(self.results_dir) first_checkpoint = Checkpoint.get_checkpoint(self.results_dir, time=0) - species_pop, aggregate_state = first_checkpoint.state + species_pop, observables, aggregate_state = self.get_state_components(first_checkpoint.state) species_ids = species_pop.keys() population_df = pandas.DataFrame(index=checkpoints, columns=species_ids, dtype=numpy.float64) + observable_ids = observables.keys() + observables_df = pandas.DataFrame(index=checkpoints, columns=observable_ids, dtype=numpy.float64) + compartments = list(aggregate_state['compartments'].keys()) properties = list(aggregate_state['compartments'][compartments[0]].keys()) compartment_property_tuples = list(zip(compartments, properties)) @@ -132,10 +142,14 @@ def convert_checkpoints(self): for time in Checkpoint.list_checkpoints(self.results_dir): checkpoint = Checkpoint.get_checkpoint(self.results_dir, time=time) - species_populations, aggregate_state = checkpoint.state + species_populations, observables, aggregate_state = self.get_state_components(checkpoint.state) + for species_id,population in species_populations.items(): population_df.loc[time, species_id] = population + for observable_id, observable in observables.items(): + observables_df.loc[time, observable_id] = observable + compartment_states = aggregate_state['compartments'] for compartment_id,agg_states in compartment_states.items(): for property,value in agg_states.items(): @@ -143,4 +157,4 @@ def convert_checkpoints(self): random_states_s[time] = pickle.dumps(checkpoint.random_state) - return (population_df, aggregate_states_df, random_states_s) + return (population_df, observables_df, aggregate_states_df, random_states_s)