To more faithfully represent initial population distributions, have s…

…ample_copy_num_from_concentration() return floating populations and have DynamicSpeciesState accept them and only round those not modeled by continuous submodels; Addresses issue #57 "Improve handling of DistributionInitConcentration"

tests/submodels/test_dynamic_submodel.py

@@ -67,13 +67,6 @@ def test_get_num_submodels(self):
         for dynamic_submodel in self.dynamic_submodels.values():
             self.assertEqual(dynamic_submodel.get_num_submodels(), 1)
 
-    def expected_molar_conc(self, dynamic_submodel, species_id):
-        species = list(filter(lambda s: s.id == species_id, dynamic_submodel.species))[0]
-        init_volume = species.compartment.init_volume.mean
-        copy_num = ModelUtilities.sample_copy_num_from_concentration(species, init_volume, RandomStateManager.instance())
-        volume = dynamic_submodel.dynamic_compartments[species.compartment.id].volume()
-        return copy_num / (volume * Avogadro)
-
     def test_calc_reaction_rates(self):
         # set standard deviation of initial conc. to 0
         self.setUp(std_init_concentrations=0.)

tests/test_model_utilities.py

@@ -110,12 +110,15 @@ def test_sample_copy_num_from_concentration(self):
         self.assertEqual(copy_number, 10**-6 * conc_value * Avogadro)
         copy_number = conc_to_molecules(species['nanomolar'], species['nanomolar'].compartment.init_volume.mean, random_state)
         self.assertEqual(copy_number, 10**-9 * conc_value * Avogadro)
-        copy_number = conc_to_molecules(species['picomolar'], species['picomolar'].compartment.init_volume.mean, random_state)
+        copy_number = round(conc_to_molecules(species['picomolar'],
+                                              species['picomolar'].compartment.init_volume.mean, random_state))
         self.assertEqual(copy_number, 10**-12 * conc_value * Avogadro)
-        copy_number = conc_to_molecules(species['femtomolar'], species['femtomolar'].compartment.init_volume.mean, random_state)
-        self.assertAlmostEqual(copy_number, 10**-15 * conc_value * Avogadro, delta=1)
-        copy_number = conc_to_molecules(species['attomolar'], species['attomolar'].compartment.init_volume.mean, random_state)
-        self.assertAlmostEqual(copy_number, 10**-18 * conc_value * Avogadro, delta=1)
+        copy_number = round(conc_to_molecules(species['femtomolar'],
+                                              species['femtomolar'].compartment.init_volume.mean, random_state))
+        self.assertEqual(copy_number, 10**-15 * conc_value * Avogadro)
+        copy_number = round(conc_to_molecules(species['attomolar'],
+                                              species['attomolar'].compartment.init_volume.mean, random_state))
+        self.assertEqual(copy_number, 10**-18 * conc_value * Avogadro)
         copy_number = conc_to_molecules(species['no_concentration'],
                                         species['no_concentration'].compartment.init_volume.mean, random_state)
         self.assertEqual(copy_number, 0)

tests/test_species_populations.py

@@ -296,7 +296,8 @@ def setUp(self):
         self.local_species_pop_no_init_pop_slope = \
             LocalSpeciesPopulation('test',
                                    self.init_populations,
-                                   self.molecular_weights)
+                                   self.molecular_weights,
+                                   random_state=RandomStateManager.instance())
 
         self.MODEL_FILENAME = os.path.join(os.path.dirname(__file__), 'fixtures',
                                            'test_model_for_access_species_populations.xlsx')
@@ -309,7 +310,8 @@ def add_initial_continuous_adjustments(self, lsp, cont_submodel_id, init_populat
 
     def test_init(self):
         self.assertEqual(self.local_species_pop_no_init_pop_slope._all_species(), set(self.species_ids))
-        an_LSP = LocalSpeciesPopulation('test', {}, {}, retain_history=False)
+        an_LSP = LocalSpeciesPopulation('test', {}, {}, random_state=RandomStateManager.instance(),
+                                        retain_history=False)
         s1_id = 's1[c]'
         mw = 1.5
         an_LSP.init_cell_state_species(s1_id, 2, mw)
@@ -642,7 +644,8 @@ def test_invalid_weights(self):
         init_populations = dict(zip(species_ids, [1]*num_mws))
         local_species_pop = LocalSpeciesPopulation('test_invalid_weights',
                                                    init_populations,
-                                                   molecular_weights)
+                                                   molecular_weights,
+                                                   random_state=RandomStateManager.instance())
         ids_w_bad_mws = species_ids[:len(bad_molecular_weights)]
         self.assertEqual(local_species_pop.invalid_weights(), set(ids_w_bad_mws))
         ids_w_bad_or_no_mw = ['x'] + ids_w_bad_mws
@@ -785,17 +788,21 @@ def test_dynamic_species_state_init(self):
         with self.assertRaisesRegex(AssertionError, 'cont_submodel_ids must be None or a list'):
             DynamicSpeciesState('s0[c]', self.random_state, pop, cont_submodel_ids=3)
 
-        with self.assertRaisesRegex(AssertionError,
-                                    re.escape("a species population modeled only by discrete submodel(s) "
-                                              "must be a non-negative integer")):
-            DynamicSpeciesState('s0[c]', self.random_state, 1.5)
+    def test_modeled_discretely(self):
+
+        pop = 3.4
+        s0 = DynamicSpeciesState('s0[c]', self.random_state, pop)
+        self.assertIn(s0.last_population, [round(pop), round(pop) + 1])
+        self.assertFalse(s0.modeled_continuously())
 
     def test_modeled_continuously(self):
 
         s0 = DynamicSpeciesState('s0[c]', self.random_state, 3)
         self.assertFalse(s0.modeled_continuously())
 
-        s1 = DynamicSpeciesState('s1[c]', self.random_state, 1.5, cont_submodel_ids=['ode'])
+        pop = 3.4
+        s1 = DynamicSpeciesState('s1[c]', self.random_state, pop, cont_submodel_ids=['ode'])
+        self.assertEqual(pop, s1.last_population)
         self.assertTrue(s1.modeled_continuously())
 
     def test__all_slopes_set(self):
@@ -1176,7 +1183,7 @@ def test_get_population_with_interpolation_on_or_off(self):
         config_multialgorithm['interpolate'] = existing_interpolate
 
     def test_temp_populations(self):
-        ds = DynamicSpeciesState('s[c]', None, 0)
+        ds = DynamicSpeciesState('s[c]', self.random_state, 0)
         self.assertEqual(ds.get_temp_population_value(), None)
         population = 3
         ds.set_temp_population_value(population)

wc_sim/model_utilities.py

@@ -43,7 +43,7 @@ def find_private_species(model, return_ids=False):
             return_ids (:obj:`boolean`, optional): if set, return object ids rather than references
 
         Returns:
-            dict: a dict that maps each submodel to a set containing the species
+            :obj:`dict`: a dict that maps each submodel to a set containing the species
                 modeled by only the submodel.
         """
         species_to_submodels = collections.defaultdict(list)
@@ -78,7 +78,7 @@ def find_shared_species(model, return_ids=False):
             return_ids (:obj:`boolean`, optional): if set, return object ids rather than references
 
         Returns:
-            set: a set containing the shared species.
+            :obj:`set`: a set containing the shared species.
         """
         all_species = model.get_species()
 
@@ -95,8 +95,8 @@ def find_shared_species(model, return_ids=False):
     def sample_copy_num_from_concentration(species, volume, random_state):
         """ Provide the initial copy number of `species` from its specified value
 
-        The initial copy number is sampled from a specified distribution in molecules or molarity.
-        Copy number is rounded to the closest integer to avoid truncating small populations.
+        The initial copy number is sampled from a specified distribution whose mean is given
+        in molecules or molarity.
 
         Args:
             species (:obj:`Species`): a `Species` instance; the `species.concentration.units` must
@@ -106,7 +106,7 @@ def sample_copy_num_from_concentration(species, volume, random_state):
                 concentrations
 
         Returns:
-            `int`: the `species'` copy number
+            :obj:`float`: the `species'` copy number
 
         Raises:
             :obj:`ValueError`: if the concentration uses illegal or unsupported units
@@ -129,14 +129,13 @@ def sample_copy_num_from_concentration(species, volume, random_state):
 
             try:
                 scale = units.to(unit_registry.parse_units('molecule'))
-                return round(scale.magnitude * conc)
+                return scale.magnitude * conc
             except pint.DimensionalityError:
                 pass
 
             try:
                 scale = units.to(unit_registry.parse_units('M'))
-                # population must be rounded to the closest integer to avoid truncating small populations
-                return int(round(scale.magnitude * conc * volume * Avogadro))
+                return scale.magnitude * conc * volume * Avogadro
             except pint.DimensionalityError as error:
                 pass
 
@@ -153,7 +152,7 @@ def get_species_types(species_ids):
             species_ids (:obj:`iterator`): an iterator that provides specie ids
 
         Returns:
-            `list`: an iterator over the specie type ids in `species_ids`
+            :obj:`list`: an iterator over the specie type ids in `species_ids`
         """
         species_types = set()
         species_types_list = []

wc_sim/species_populations.py

@@ -602,7 +602,7 @@ def __init__(self, name, initial_population, molecular_weights, cont_submodel_id
                                              sim_time=self.time)
 
     def init_cell_state_species(self, species_id, population, molecular_weight, cont_submodel_ids=None):
-        """ Initialize a species with the given population and ids of continuous submodels that model it
+        """ Initialize a species with the given population and, optionally, the ids of continuous submodels that model it
 
         Add a species to the cell state.
 
@@ -1499,6 +1499,8 @@ def __init__(self, species_name, random_state, initial_population, cont_submodel
                  record_history=False, default_rounding=None):
         """ Initialize a species object, defaulting to a simulation time start time of 0
 
+        If a species is not modeled continuously then its initial population is stochastically rounded to an integer.
+
         Args:
             species_name (:obj:`str`): the species' name; not logically needed, but helpful for error
                 reporting, logging, debugging, etc.
@@ -1517,12 +1519,13 @@ def __init__(self, species_name, random_state, initial_population, cont_submodel
         assert cont_submodel_ids is None or isinstance(cont_submodel_ids, list), \
             (f"DynamicSpeciesState '{species_name}': cont_submodel_ids must be None or a list "
              f"but is a(n) {type(cont_submodel_ids).__name__}")
-        # if a species is not modeled continuously then its initial population must be a non-negative integer
-        assert cont_submodel_ids is not None or float(initial_population).is_integer(), \
-            (f"DynamicSpeciesState '{species_name}': a species population modeled only by discrete submodel(s) "
-             f"must be a non-negative integer, but {initial_population} isn't")
+
         self.random_state = random_state
+
         self.last_population = initial_population
+        if cont_submodel_ids is None:
+            self.last_population = self.random_state.round(initial_population)
+
         if cont_submodel_ids is not None:
             # initialize population_slopes to None for each submodel
             self.population_slopes = dict.fromkeys(cont_submodel_ids, None)