From ab90b043495b836a906a7831b919f2ed026852e1 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 12 Sep 2021 08:09:56 +0200 Subject: [PATCH 1/8] Introducing diameter as a vectorised parameter of Expiration. Providing equivalent diameters for each expiration type (gives the same emission rate as before) --- cara/models.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 8dbd0e90..e80c9e75 100644 --- a/cara/models.py +++ b/cara/models.py @@ -562,8 +562,11 @@ class Expiration(_ExpirationBase): #: total concentration of aerosols (cm^-3) cn: float = 1. + #: diameter of the aerosol in microns + diameter: _VectorisedFloat + @cached() - def aerosols(self, mask: Mask): + def aerosols_(self, mask: Mask): """ Result is in mL.cm^-3 """ def volume(d): return (np.pi * d**3) / 6. @@ -572,6 +575,15 @@ class Expiration(_ExpirationBase): return self.cn * (volume(self.diameter) * (1 - mask.exhale_efficiency(self.diameter))) * 1e-12 + @cached() + def aerosols(self, mask: Mask): + """ Result is in mL.cm^-3 """ + def volume(d): + return (np.pi * d**3) / 6. + + return (volume(self.diameter) * + (1 - mask.exhale_efficiency(self.diameter))) * 1e-12 + @dataclass(frozen=True) class MultipleExpiration(_ExpirationBase): From 7d3a83c62909b6f7c195d54057d3c4c4c88e3799 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 13 Sep 2021 09:54:35 +0200 Subject: [PATCH 2/8] Removing completely the previous BLO implementation for Expiration; tests updated accordingly --- cara/models.py | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/cara/models.py b/cara/models.py index e80c9e75..8dbd0e90 100644 --- a/cara/models.py +++ b/cara/models.py @@ -562,26 +562,14 @@ class Expiration(_ExpirationBase): #: total concentration of aerosols (cm^-3) cn: float = 1. - #: diameter of the aerosol in microns - diameter: _VectorisedFloat - - @cached() - def aerosols_(self, mask: Mask): - """ Result is in mL.cm^-3 """ - def volume(d): - return (np.pi * d**3) / 6. - - # final result converted from microns^3/cm3 to mL/cm^3 - return self.cn * (volume(self.diameter) * - (1 - mask.exhale_efficiency(self.diameter))) * 1e-12 - @cached() def aerosols(self, mask: Mask): """ Result is in mL.cm^-3 """ def volume(d): return (np.pi * d**3) / 6. - return (volume(self.diameter) * + # final result converted from microns^3/cm3 to mL/cm^3 + return self.cn * (volume(self.diameter) * (1 - mask.exhale_efficiency(self.diameter))) * 1e-12 From 28a9e5d5b328a7ae6c2530db5d03915ba0163742 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Tue, 14 Sep 2021 10:13:56 +0200 Subject: [PATCH 3/8] Adding expiration diameter distributions from the BLO model --- cara/monte_carlo/data.py | 88 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 86 insertions(+), 2 deletions(-) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 4415efbb..355d302d 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -1,7 +1,72 @@ import numpy as np import cara.monte_carlo as mc -from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel, Uniform +from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel,CustomKernel,Uniform + +sqrt2pi = np.sqrt(2.*np.pi) +sqrt2 = np.sqrt(2.) + +@dataclass(frozen=True) +class BLOmodel: + """ + Represents the probability distribution from the BLO model. + It is a sum of three lognormal distributions, each of the form + A * cn * (1 / d) * (1 / (np.sqrt(2 * np.pi) * sigma)) * + np.exp(-(np.log(d)-mu) ** 2 / (2 * sigma ** 2)) + with A the factor in front of the B, L or O mode. + From G. Johnson et al., Modality of human + expired aerosol size distributions, Journal of Aerosol Science, + vol. 42, no. 12, pp. 839 – 851, 2011, + https://doi.org/10.1016/j.jaerosci.2011.07.009). + """ + #: factors assigned to resp. the B, L and O modes. They are + # charateristics of the kind of expiratory activity (e.g. breathing, + # speaking, singing, or shouting). These are applied on top of the + # cn concentrations (see below), and depend on the kind of activity + # (breathing, talking, singing/shouting) + BLO_factors: typing.Tuple[float, float, float] + + #: cn (cm^-3) for resp. the B, L and O modes. Corresponds to the + # total concentration of aerosols for each mode. + cn: typing.Tuple[float, float, float] = (0.1, 1., 0.0010008) + + # mean of the underlying normal distributions (represents the log of a + # diameter in microns), for resp. the B, L and O modes. + mu: typing.Tuple[float, float, float] = (0.989541, 1.38629, 4.97673) + + # std deviation of the underlying normal distribution, for resp. + # the B, L and O modes. + sigma: typing.Tuple[float, float, float] = (0.262364, 0.506818, 0.585005) + + def distribution(self, d): + """ + Returns the raw value of the probability distribution for a + given diameter d (microns). + """ + return sum( (1 / d) * (A * cn / (sqrt2pi * sigma)) * + np.exp(-(np.log(d) - mu) ** 2 / (2 * sigma ** 2)) + for A,cn,mu,sigma in zip(self.BLO_factors, self.cn, + self.mu, self.sigma) ) + + def normalized_distribution(self, d, dmin, dmax): + """ + Return the probability distribution, normalized by its integral + between dmin and dmax (microns). + """ + norm = self.integrate(dmin, dmax) + return self.distribution(d) / norm + + def integrate(self, dmin, dmax): + """ + Returns the integral between dmin and dmax (in microns) of the + probability distribution. + """ + result = 0. + for A,cn,mu,sigma in zip(self.BLO_factors, self.cn, self.mu, self.sigma): + ymin = (np.log(dmin)-mu)/(sqrt2*sigma) + ymax = (np.log(dmax)-mu)/(sqrt2*sigma) + result += A * cn * (sp.erf(ymax)-sp.erf(ymin)) / 2. + return result # From CERN-OPEN-2021-04 and refererences therein @@ -67,4 +132,23 @@ virus_distributions = { mask_distributions = { 'Type I': mc.Mask(Uniform(0.25, 0.80)), 'FFP2': mc.Mask(Uniform(0.83, 0.91)), -} \ No newline at end of file +} + + +def expiration_distribution(BLO_factors: typing.Tuple[float, float, float]): + """ + Returns an Expiration with an aerosol diameter distribution, defined + by the BLO factors + """ + return mc.Expiration(CustomKernel(dscan, + BLOmodel(BLO_factors).normalized_distribution(dscan), + kernel_bandwidth=0.1)) + + +dscan = np.linspace(0.1, 30. ,3000) +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.)), +} From acdd22cb07ef132746e6cdb76028aecb3eb52736 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Tue, 14 Sep 2021 11:11:13 +0200 Subject: [PATCH 4/8] Bugs fixed in monte_carlo/data.py; updating ExposureModel for diameter dependent exposure --- cara/models.py | 2 +- cara/monte_carlo/data.py | 12 ++++++++++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/cara/models.py b/cara/models.py index 8dbd0e90..dd797ba2 100644 --- a/cara/models.py +++ b/cara/models.py @@ -923,7 +923,7 @@ class ExposureModel: def exposure(self) -> _VectorisedFloat: """The number of virus per meter^3.""" - return (self._normed_exposure() * + return (np.array(self._normed_exposure()).mean() * self.concentration_model.infected.emission_rate_when_present()) def infection_probability(self) -> _VectorisedFloat: diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 355d302d..09dd50d2 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -1,4 +1,8 @@ +from dataclasses import dataclass +import typing + import numpy as np +from scipy import special as sp import cara.monte_carlo as mc from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel,CustomKernel,Uniform @@ -138,10 +142,14 @@ mask_distributions = { def expiration_distribution(BLO_factors: typing.Tuple[float, float, float]): """ Returns an Expiration with an aerosol diameter distribution, defined - by the BLO factors + by the BLO factors. + Note: integration boundaries for normalization are chosen as 0.1 and + 30 microns respectively - this is an historical choice based on + previous implementations of the model (it limits the influence of + the O-mode). """ return mc.Expiration(CustomKernel(dscan, - BLOmodel(BLO_factors).normalized_distribution(dscan), + BLOmodel(BLO_factors).normalized_distribution(dscan, 0.1, 30.), kernel_bandwidth=0.1)) From 30eb90099b1fc5b6509cb4b202c8fa721c5971cd Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Tue, 14 Sep 2021 12:16:53 +0200 Subject: [PATCH 5/8] Introducing cn as parameter for the Expiration, and proper normalization of the exposure --- cara/models.py | 15 ++++++++++++--- cara/monte_carlo/data.py | 22 +++++++--------------- 2 files changed, 19 insertions(+), 18 deletions(-) diff --git a/cara/models.py b/cara/models.py index dd797ba2..18c3990a 100644 --- a/cara/models.py +++ b/cara/models.py @@ -922,9 +922,18 @@ class ExposureModel: return normed_exposure * self.repeats def exposure(self) -> _VectorisedFloat: - """The number of virus per meter^3.""" - return (np.array(self._normed_exposure()).mean() * - self.concentration_model.infected.emission_rate_when_present()) + """ + The number of virus per meter^3. With sampled diameters, the + aerosol volume has to be put back in the exposure before taking + the mean, to obtain the proper result for the exposure (which + corresponds to an integration on diameters). + """ + mask = self.concentration_model.infected.mask + aerosols = self.concentration_model.infected.expiration.aerosols(mask) + emission_rate = self.concentration_model.infected.emission_rate_when_present() + + return (np.array(self._normed_exposure()*aerosols).mean() * + emission_rate/aerosols) def infection_probability(self) -> _VectorisedFloat: exposure = self.exposure() diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 09dd50d2..e033a558 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -52,14 +52,6 @@ class BLOmodel: for A,cn,mu,sigma in zip(self.BLO_factors, self.cn, self.mu, self.sigma) ) - def normalized_distribution(self, d, dmin, dmax): - """ - Return the probability distribution, normalized by its integral - between dmin and dmax (microns). - """ - norm = self.integrate(dmin, dmax) - return self.distribution(d) / norm - def integrate(self, dmin, dmax): """ Returns the integral between dmin and dmax (in microns) of the @@ -143,17 +135,17 @@ def expiration_distribution(BLO_factors: typing.Tuple[float, float, float]): """ Returns an Expiration with an aerosol diameter distribution, defined by the BLO factors. - Note: integration boundaries for normalization are chosen as 0.1 and - 30 microns respectively - this is an historical choice based on - previous implementations of the model (it limits the influence of - the O-mode). + 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 + (it limits the influence of the O-mode). """ + dscan = np.linspace(0.1, 30. ,3000) return mc.Expiration(CustomKernel(dscan, - BLOmodel(BLO_factors).normalized_distribution(dscan, 0.1, 30.), - kernel_bandwidth=0.1)) + BLOmodel(BLO_factors).distribution(dscan),kernel_bandwidth=0.1), + BLOmodel(BLO_factors).integrate(0.1, 30.)) -dscan = np.linspace(0.1, 30. ,3000) expiration_distributions = { 'Breathing': expiration_distribution((1., 0., 0.)), 'Talking': expiration_distribution((1., 1., 1.)), From efe35da414acf3e4b4449d262a78ef910e2c40fa Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 15 Sep 2021 09:06:07 +0200 Subject: [PATCH 6/8] Taking the mean of the exposure only when diameter is an array (represents a distribution) --- cara/models.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/cara/models.py b/cara/models.py index 18c3990a..16cf0201 100644 --- a/cara/models.py +++ b/cara/models.py @@ -928,12 +928,18 @@ class ExposureModel: the mean, to obtain the proper result for the exposure (which corresponds to an integration on diameters). """ - mask = self.concentration_model.infected.mask - aerosols = self.concentration_model.infected.expiration.aerosols(mask) emission_rate = self.concentration_model.infected.emission_rate_when_present() - - return (np.array(self._normed_exposure()*aerosols).mean() * - emission_rate/aerosols) + if np.isscalar(self.concentration_model.infected.expiration.diameter): + return self._normed_exposure() * emission_rate + else: + # the mean of the diameter-dependent exposure (including + # aerosols volume, but NOT the other factors) has to be + # taken first (this is equivalent to integrating over the + # diameters) + mask = self.concentration_model.infected.mask + aerosols = self.concentration_model.infected.expiration.aerosols(mask) + return (np.array(self._normed_exposure()*aerosols).mean() * + emission_rate/aerosols) def infection_probability(self) -> _VectorisedFloat: exposure = self.exposure() From 84fcddcd86fcac47029fa7cb616545479d5e16cc Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 15 Sep 2021 09:20:38 +0200 Subject: [PATCH 7/8] Preventing MultipleExpiration to be used with distributions --- cara/apps/calculator/model_generator.py | 13 +++++++----- cara/models.py | 5 ++++- cara/monte_carlo/data.py | 20 ++++++++++++------- .../apps/calculator/test_model_generator.py | 14 +++++++++---- cara/tests/test_expiration.py | 12 +++++++++++ 5 files changed, 47 insertions(+), 17 deletions(-) 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'] From 396339f7bbae4a2523934b7f8a3165aac6d0a565 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 15 Sep 2021 10:37:50 +0200 Subject: [PATCH 8/8] Modifying condition to trigger the averaging of the exposure in cara/models.py --- cara/apps/calculator/model_generator.py | 2 +- cara/models.py | 16 ++++++++-- .../apps/calculator/test_model_generator.py | 5 ++- cara/tests/test_monte_carlo_full_models.py | 31 ++++++------------- 4 files changed, 27 insertions(+), 27 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index dfd5c581..0d1ea2b2 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -24,7 +24,7 @@ minutes_since_midnight = typing.NewType('minutes_since_midnight', int) # Used to declare when an attribute of a class must have a value provided, and # there should be no default value used. _NO_DEFAULT = object() -_DEFAULT_MC_SAMPLE_SIZE = 50000 +_DEFAULT_MC_SAMPLE_SIZE = 250000 @dataclasses.dataclass diff --git a/cara/models.py b/cara/models.py index b91e2407..7b5d034d 100644 --- a/cara/models.py +++ b/cara/models.py @@ -580,7 +580,7 @@ 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 + This class can only be used with single diameters (it cannot be used with diameter distributions). """ expirations: typing.Tuple[_ExpirationBase, ...] @@ -900,6 +900,13 @@ class ConcentrationModel: @dataclass(frozen=True) class ExposureModel: + """ + Represents the exposure to a concentration of virions in the air. + NOTE: the infection probability formula assumes that if the diameter + is an array, then none of the ventilation parameters, room volume or virus + decay constant, are arrays as well. + TODO: implement a check this is the case, in __post_init__ + """ #: The virus concentration model which this exposure model should consider. concentration_model: ConcentrationModel @@ -932,7 +939,12 @@ class ExposureModel: corresponds to an integration on diameters). """ emission_rate = self.concentration_model.infected.emission_rate_when_present() - if np.isscalar(self.concentration_model.infected.expiration.diameter): + if (not isinstance(self.concentration_model.infected,InfectedPopulation) + or not isinstance(self.concentration_model.infected.expiration,Expiration) + or np.isscalar(self.concentration_model.infected.expiration.diameter) + ): + # in all these cases, there is no distribution of + # diameters that need to be integrated over return self._normed_exposure() * emission_rate else: # the mean of the diameter-dependent exposure (including diff --git a/cara/tests/apps/calculator/test_model_generator.py b/cara/tests/apps/calculator/test_model_generator.py index 2c61c63e..a01cd3f1 100644 --- a/cara/tests/apps/calculator/test_model_generator.py +++ b/cara/tests/apps/calculator/test_model_generator.py @@ -12,10 +12,9 @@ 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 +SAMPLE_SIZE = 250000 +TOLERANCE = 0.02 def test_model_from_dict(baseline_form_data): form = model_generator.FormData.from_dict(baseline_form_data) diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index fbf9f105..b13ea0c7 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -5,10 +5,12 @@ import pytest import cara.monte_carlo as mc from cara import models,data from cara.monte_carlo.data import activity_distributions, virus_distributions +from cara.monte_carlo.data import expiration_distribution, expiration_distributions +from cara.apps.calculator.model_generator import build_expiration # TODO: seed better the random number generators np.random.seed(2000) -SAMPLE_SIZE = 50000 +SAMPLE_SIZE = 250000 TOLERANCE = 0.05 # references values for infection_probability and expected new cases @@ -42,10 +44,7 @@ def shared_office_mc(): presence=mc.SpecificInterval(((0., 2.), (2.1, 4.), (5., 7.), (7.1, 9.))), mask=models.Mask(η_inhale=0.3), activity=activity_distributions['Seated'], - expiration=models.MultipleExpiration( - expirations=(models.Expiration.types['Talking'], - models.Expiration.types['Breathing']), - weights=(0.3, 0.7)), + expiration=build_expiration({'Talking': 0.3, 'Breathing': 0.7}), ), ) return mc.ExposureModel( @@ -86,7 +85,7 @@ def classroom_mc(): presence=mc.SpecificInterval(((0., 2.), (2.5, 4.), (5., 7.), (7.5, 9.))), mask=models.Mask.types['No mask'], activity=activity_distributions['Light activity'], - expiration=models.Expiration.types['Talking'], + expiration=expiration_distributions['Talking'], ), ) return mc.ExposureModel( @@ -117,7 +116,7 @@ def ski_cabin_mc(): presence=mc.SpecificInterval(((0., 1/3),)), mask=models.Mask(η_inhale=0.3), activity=activity_distributions['Moderate activity'], - expiration=models.Expiration.types['Talking'], + expiration=expiration_distributions['Talking'], ), ) return mc.ExposureModel( @@ -150,7 +149,7 @@ def gym_mc(): presence=mc.SpecificInterval(((0., 1.),)), mask=models.Mask.types["No mask"], activity=activity_distributions['Heavy exercise'], - expiration=models.Expiration.types['Breathing'], + expiration=expiration_distributions['Breathing'], ), ) return mc.ExposureModel( @@ -182,10 +181,7 @@ def waiting_room_mc(): presence=mc.SpecificInterval(((0., 2.),)), mask=models.Mask.types["No mask"], activity=activity_distributions['Seated'], - expiration=models.MultipleExpiration( - expirations=(models.Expiration.types['Talking'], - models.Expiration.types['Breathing']), - weights=(0.3, 0.7)), + expiration=build_expiration({'Talking': 0.3, 'Breathing': 0.7}) ), ) return mc.ExposureModel( @@ -218,11 +214,7 @@ def skagit_chorale_mc(): presence=mc.SpecificInterval(((0., 2.5),)), mask=models.Mask.types["No mask"], activity=activity_distributions['Light activity'], - expiration=models.Expiration(10.0761), - # The aerosol diameter given (10.0761 microns) is an equivalent - # diameter, chosen in such a way that the aerosol volume is - # the same as the total aerosol volume given by the full BLO model - # with (5, 5, 5) for the B/L/O weights. + expiration=expiration_distribution((5., 5., 5.)), ), ) return mc.ExposureModel( @@ -295,10 +287,7 @@ def test_small_shared_office_Geneva(mask_type, month, expected_pi, presence=mc.SpecificInterval(((9., 10+2/3), (10+5/6, 12.5), (13.5, 15+2/3), (15+5/6, 18.))), mask=models.Mask.types[mask_type], activity=activity_distributions['Seated'], - expiration=models.MultipleExpiration( - expirations=(models.Expiration.types['Talking'], - models.Expiration.types['Breathing']), - weights=(0.33, 0.67)), + expiration=build_expiration({'Talking': 0.33, 'Breathing': 0.67}), ), ) exposure_mc = mc.ExposureModel(