add Observables to checkpoints and RunResults; dynamically create che…

…ckpoints to test RunResults

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))
 

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:

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,

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):
 

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')

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
 

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.)

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

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'],

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
 '''

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,15 +142,19 @@ 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():
                     aggregate_states_df.loc[time, (compartment_id, property)] = value
 
             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)