Preventing MultipleExpiration to be used with distributions

cara/apps/calculator/model_generator.py

@@ -12,6 +12,7 @@ import cara.data.weather
 import cara.monte_carlo as mc
 from .. import calculator
 from cara.monte_carlo.data import activity_distributions, virus_distributions, mask_distributions
+from cara.monte_carlo.data import expiration_distribution, expiration_BLO_factors, expiration_distributions
 
 
 LOG = logging.getLogger(__name__)
@@ -628,12 +629,14 @@ class FormData:
 
 def build_expiration(expiration_definition) -> models._ExpirationBase:
     if isinstance(expiration_definition, str):
-        return models._ExpirationBase.types[expiration_definition]
+        return expiration_distributions[expiration_definition]
     elif isinstance(expiration_definition, dict):
-        return models.MultipleExpiration(
-            tuple([build_expiration(exp) for exp in expiration_definition.keys()]),
-            tuple(expiration_definition.values())
-        )
+        total_weight = sum(expiration_definition.values())
+        BLO_factors = np.sum([
+            np.array(expiration_BLO_factors[exp_type]) * weight/total_weight
+            for exp_type, weight in expiration_definition.items()
+            ], axis=0)
+        return expiration_distribution(tuple(BLO_factors))
 
 
 def baseline_raw_form_data():

cara/models.py

@@ -580,7 +580,8 @@ class MultipleExpiration(_ExpirationBase):
     Group together different modes of expiration, that represent
     each the main expiration mode for a certain fraction of time (given by
     the weights).
-
+    Obsolet class that can only be used with single diameters (it cannot
+    be used with diameter distributions).
     """
     expirations: typing.Tuple[_ExpirationBase, ...]
     weights: typing.Tuple[float, ...]
@@ -589,6 +590,8 @@ class MultipleExpiration(_ExpirationBase):
         if len(self.expirations) != len(self.weights):
             raise ValueError("expirations and weigths should contain the"
                              "same number of elements")
+        if not all(np.isscalar(e.diameter) for e in self.expirations):
+            raise ValueError("diameters should all be scalars")
 
     def aerosols(self, mask: Mask):
         return np.array([

cara/monte_carlo/data.py

@@ -131,10 +131,10 @@ mask_distributions = {
 }
 
 
-def expiration_distribution(BLO_factors: typing.Tuple[float, float, float]):
+def expiration_distribution(BLO_factors):
     """
     Returns an Expiration with an aerosol diameter distribution, defined
-    by the BLO factors.
+    by the BLO factors (a length-3 tuple).
     The total concentration of aerosols is computed by integrating
     the distribution between 0.1 and 30 microns - these boundaries are
     an historical choice based on previous implementations of the model
@@ -146,9 +146,15 @@ def expiration_distribution(BLO_factors: typing.Tuple[float, float, float]):
                 BLOmodel(BLO_factors).integrate(0.1, 30.))
 
 
-expiration_distributions = {
-    'Breathing': expiration_distribution((1., 0., 0.)),
-    'Talking':   expiration_distribution((1., 1., 1.)),
-    'Singing':   expiration_distribution((1., 5., 5.)),
-    'Shouting':  expiration_distribution((1., 5., 5.)),
+expiration_BLO_factors = {
+    'Breathing': (1., 0., 0.),
+    'Talking':   (1., 1., 1.),
+    'Singing':   (1., 5., 5.),
+    'Shouting':  (1., 5., 5.),
+}
+
+
+expiration_distributions = {
+    exp_type: expiration_distribution(BLO_factors)
+    for exp_type,BLO_factors in expiration_BLO_factors.items()
 }

cara/tests/apps/calculator/test_model_generator.py

@@ -9,7 +9,13 @@ from cara.apps.calculator import model_generator
 from cara.apps.calculator.model_generator import _hours2timestring
 from cara.apps.calculator.model_generator import minutes_since_midnight
 from cara import models
+from cara.monte_carlo.data import expiration_distributions
 
+# TODO: seed better the random number generators
+# this is only for test_blend_expiration
+np.random.seed(2000)
+SAMPLE_SIZE = 500000
+TOLERANCE = 0.01
 
 def test_model_from_dict(baseline_form_data):
     form = model_generator.FormData.from_dict(baseline_form_data)
@@ -31,11 +37,11 @@ def test_model_from_dict_invalid(baseline_form_data):
 )
 def test_blend_expiration(mask_type):
     blend = {'Breathing': 2, 'Talking': 1}
-    r = model_generator.build_expiration(blend)
+    r = model_generator.build_expiration(blend).build_model(SAMPLE_SIZE)
     mask = models.Mask.types[mask_type]
-    expected = (models.Expiration.types['Breathing'].aerosols(mask)*2/3. +
-                models.Expiration.types['Talking'].aerosols(mask)/3.)
-    npt.assert_allclose(r.aerosols(mask), expected)
+    expected = (expiration_distributions['Breathing'].build_model(SAMPLE_SIZE).aerosols(mask).mean()*2/3. +
+                expiration_distributions['Talking'].build_model(SAMPLE_SIZE).aerosols(mask).mean()/3.)
+    npt.assert_allclose(r.aerosols(mask).mean(), expected, rtol=TOLERANCE)
 
 
 def test_ventilation_slidingwindow(baseline_form: model_generator.FormData):

cara/tests/test_expiration.py

@@ -18,6 +18,18 @@ def test_multiple_wrong_weight_size():
         e = models.MultipleExpiration([e_base, e_base], weights)
 
 
+def test_multiple_wrong_diameters():
+    weights = (1., 2., 3.)
+    e1 = models.Expiration(np.array([1., 1.]))
+    e2 = models.Expiration(1.)
+    e3 = models.Expiration(2.)
+    with pytest.raises(
+            ValueError,
+            match=re.escape("diameters should all be scalars")
+    ):
+        e = models.MultipleExpiration([e1, e2, e3], weights)
+
+
 def test_multiple():
     weights = (1., 1.)
     mask = models.Mask.types['Type I']