diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index fb34392d..ad65795e 100644 --- a/cara/apps/calculator/model_generator.py +++ b/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': diff --git a/cara/models.py b/cara/models.py index 0943bd87..c2500d68 100644 --- a/cara/models.py +++ b/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: """ diff --git a/cara/tests/apps/calculator/test_model_generator.py b/cara/tests/apps/calculator/test_model_generator.py index 6c5050e4..eb216727 100644 --- a/cara/tests/apps/calculator/test_model_generator.py +++ b/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): diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py new file mode 100644 index 00000000..7638d613 --- /dev/null +++ b/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. + )