Merge branch 'feature/mult_expiration' into 'master'

MultipleExpriation implementation

See merge request cara/cara!188

cara/apps/calculator/model_generator.py

@@ -537,38 +537,16 @@ class FormData:
         )
 
 
-def build_expiration(expiration_definition) -> models.Expiration:
+def build_expiration(expiration_definition) -> models._ExpirationBase:
     if isinstance(expiration_definition, str):
-        return models.Expiration.types[expiration_definition]
+        return models._ExpirationBase.types[expiration_definition]
     elif isinstance(expiration_definition, dict):
-        return expiration_blend({
-            build_expiration(exp): amount
-            for exp, amount in expiration_definition.items()
-            }
+        return models.MultipleExpiration(
+            tuple([build_expiration(exp) for exp in expiration_definition.keys()]),
+            tuple(expiration_definition.values())
         )
 
 
-def expiration_blend(expiration_weights: typing.Dict[models.Expiration, int]) -> models.Expiration:
-    """
-    Combine together multiple types of Expiration, using a weighted mean to
-    compute their ejection factor and particle sizes.
-
-    """
-    ejection_factor = np.zeros(4)
-    particle_sizes = np.zeros(4)
-
-    total_weight = 0
-    for expiration, weight in expiration_weights.items():
-        total_weight += weight
-        ejection_factor += np.array(expiration.ejection_factor) * weight
-        particle_sizes += np.array(expiration.particle_sizes) * weight
-
-    r_ejection_factor: typing.Tuple[float, float, float, float] = tuple(ejection_factor/total_weight)  # type: ignore
-    r_particle_sizes: typing.Tuple[float, float, float, float] = tuple(particle_sizes/total_weight)  # type: ignore
-
-    return models.Expiration(ejection_factor=r_ejection_factor, particle_sizes=r_particle_sizes)
-
-
 def model_from_form(form: FormData) -> models.ExposureModel:
     # Initializes room with volume either given directly or as product of area and height
     if form.volume_type == 'room_volume_explicit':

cara/models.py

@@ -511,13 +511,31 @@ Mask.types = {
 
 
 @dataclass(frozen=True)
-class Expiration:
+class _ExpirationBase:
+    """
+    Represents the expiration of aerosols by a person.
+    Subclasses of _ExpirationBase represent different models.
+    """
+    #: Pre-populated examples of Masks.
+    types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]]
+
+    def aerosols(self, mask: Mask):
+        # total volume of aerosols expired per volume of air (mL/cm^3).
+        raise NotImplementedError("Subclass must implement")
+
+
+@dataclass(frozen=True)
+class Expiration(_ExpirationBase):
+    """
+    Simple model based on four different sizes of particles emitted,
+    with different ejection factors. See Fig. 4 in L. Morawska et al,
+    Size distribution and sites of origin of droplets expelled from the
+    human respiratory tract during expiratory activities,
+    Aerosol Science 40 (2009) pp. 256 - 269.
+    """
     ejection_factor: typing.Tuple[float, ...]
     particle_sizes: typing.Tuple[float, ...] = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4)  # In cm.
 
-    #: Pre-populated examples of Expiration.
-    types: typing.ClassVar[typing.Dict[str, "Expiration"]]
-
     def aerosols(self, mask: Mask):
         def volume(diameter):
             return (4 * np.pi * (diameter/2)**3) / 3
@@ -529,7 +547,31 @@ class Expiration:
         return total
 
 
-Expiration.types = {
+@dataclass(frozen=True)
+class MultipleExpiration(_ExpirationBase):
+    """
+    Represents an expiration of aerosols.
+    Group together different modes of expiration, that represent
+    each the main expiration mode for a certain fraction of time (given by
+    the weights).
+
+    """
+    expirations: typing.Tuple[_ExpirationBase, ...]
+    weights: typing.Tuple[float, ...]
+
+    def __post_init__(self):
+        if len(self.expirations) != len(self.weights):
+            raise ValueError("expirations and weigths should contain the"
+                             "same number of elements")
+
+    def aerosols(self, mask: Mask):
+        return np.array([
+            weight * expiration.aerosols(mask) / sum(self.weights)
+            for weight,expiration in zip(self.weights,self.expirations)
+        ]).sum(axis=0)
+
+
+_ExpirationBase.types = {
     'Breathing': Expiration((0.084, 0.009, 0.003, 0.002)),
     'Whispering': Expiration((0.11, 0.014, 0.004, 0.002)),
     'Talking': Expiration((0.236, 0.068, 0.007, 0.011)),
@@ -585,7 +627,7 @@ class InfectedPopulation(Population):
     virus: Virus
 
     #: The type of expiration that is being emitted whilst doing the activity.
-    expiration: Expiration
+    expiration: _ExpirationBase
 
     def emission_rate_when_present(self) -> _VectorisedFloat:
         """

cara/tests/apps/calculator/test_model_generator.py

@@ -8,6 +8,7 @@ from cara.apps.calculator.model_generator import minutes_since_midnight
 from cara import models
 from cara import data
 import numpy as np
+import numpy.testing as npt
 
 
 def test_model_from_dict(baseline_form_data):
@@ -24,10 +25,11 @@ def test_model_from_dict_invalid(baseline_form_data):
 def test_blend_expiration():
     blend = {'Breathing': 2, 'Talking': 1}
     r = model_generator.build_expiration(blend)
+    mask = models.Mask.types['Type I']
     expected = models.Expiration(
         (0.13466666666666668, 0.02866666666666667, 0.004333333333333334, 0.005)
     )
-    assert r == expected
+    npt.assert_almost_equal(r.aerosols(mask), expected.aerosols(mask))
 
 
 def test_ventilation_slidingwindow(baseline_form: model_generator.FormData):

cara/tests/test_expiration.py

@@ -0,0 +1,29 @@
+import re
+
+import numpy as np
+import numpy.testing as npt
+import pytest
+
+from cara import models
+
+
+def test_multiple_wrong_weight_size():
+    weights = (1., 2., 3.)
+    e_base = models.Expiration((0.084, 0.009, 0.003, 0.002))
+    with pytest.raises(
+            ValueError,
+            match=re.escape("expirations and weigths should contain the"
+                            "same number of elements")
+    ):
+        e = models.MultipleExpiration([e_base, e_base], weights)
+
+
+def test_multiple():
+    weights = (1., 2.)
+    e1 = models.Expiration((0.03, 0.02, 0.01, 0.005))
+    e2 = models.Expiration((0.05, 0.04, 0.03, 0.01))
+    e = models.MultipleExpiration([e1, e2], weights)
+    assert e.aerosols(models.Mask.types['No mask']) == (
+              e1.aerosols(models.Mask.types['No mask'])/3. +
+            2*e2.aerosols(models.Mask.types['No mask'])/3.
+    )