diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 06c29229..dfd5c581 100644 --- a/cara/apps/calculator/model_generator.py +++ b/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(): diff --git a/cara/models.py b/cara/models.py index 16cf0201..b91e2407 100644 --- a/cara/models.py +++ b/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([ diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index e033a558..71498ee0 100644 --- a/cara/monte_carlo/data.py +++ b/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() } diff --git a/cara/tests/apps/calculator/test_model_generator.py b/cara/tests/apps/calculator/test_model_generator.py index ba4f3a6a..2c61c63e 100644 --- a/cara/tests/apps/calculator/test_model_generator.py +++ b/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): diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py index e3497efe..8a600086 100644 --- a/cara/tests/test_expiration.py +++ b/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']