From 9d4ffdb89053db2715e0ebdde5c9af5e573102bd Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 26 May 2021 18:57:56 +0200 Subject: [PATCH 01/83] Introducing _MaskBase class, of which Mask is a subclass --- cara/models.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/cara/models.py b/cara/models.py index 0943bd87..55989ec5 100644 --- a/cara/models.py +++ b/cara/models.py @@ -472,8 +472,21 @@ Virus.types = { @dataclass(frozen=True) -class Mask: - #: Filtration efficiency. (In %/100) +class _MaskBase: + """ + Represents the filtration of aerosols by a mask, both inward and + outward. + The nature of the various air exchange schemes means that it is expected + for subclasses of _MaskBase to exist. + """ + def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: + # Overall efficiency, including the effect of the leaks. + raise NotImplementedError("Subclass must implement") + + +@dataclass(frozen=True) +class Mask(_MaskBase): + #: Filtration efficiency. η_exhale: _VectorisedFloat #: Leakage through side of masks. From 676fd64844cc7c83069f44ed265ef111c1e6ed4a Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 26 May 2021 23:37:10 +0200 Subject: [PATCH 02/83] Adding tests on Mask (new classes and new methods) --- cara/tests/models/test_mask.py | 58 ++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 cara/tests/models/test_mask.py diff --git a/cara/tests/models/test_mask.py b/cara/tests/models/test_mask.py new file mode 100644 index 00000000..f2069cbd --- /dev/null +++ b/cara/tests/models/test_mask.py @@ -0,0 +1,58 @@ + +import dataclasses + +import numpy as np +import numpy.testing as npt +import pytest + +from cara import models + +@pytest.mark.parametrize( + "η_inhale, expected_inhale_efficiency", + [ + [0.5, 0.5], + [np.array([0.3, 0.5]), np.array([0.3, 0.5])], + ], +) +def test_masks_inhale(η_inhale, expected_inhale_efficiency): + mask = models.Mask(η_inhale=η_inhale,η_exhale=0.95,η_leaks=0.15) + measuredmask = models.MeasuredMask(η_inhale=η_inhale) + npt.assert_equal(mask.inhale_efficiency(), + expected_inhale_efficiency) + npt.assert_equal(measuredmask.inhale_efficiency(), + expected_inhale_efficiency) + + +@pytest.mark.parametrize( + "η_exhale, η_leaks, expected_exhale_efficiency_small, expected_exhale_efficiency_large", + [ + [0.95, 0.15, 0., 0.8075], + [np.array([0.95, 1.]), 0.15, np.zeros(2), np.array([0.8075, 0.85])], + [0.95, np.array([0.15, 0.]), np.zeros(2), np.array([0.8075, 0.95])], + [np.array([0.95, 1.]), np.array([0.15, 0.]), np.zeros(2), np.array([0.8075, 1.])], + ], +) +def test_mask_exhale(η_exhale, η_leaks, expected_exhale_efficiency_small, + expected_exhale_efficiency_large): + mask = models.Mask(η_inhale=0.3,η_exhale=η_exhale,η_leaks=η_leaks) + # we test one small and one large diameter (resp. 1 and 4 microns) + npt.assert_equal(mask.exhale_efficiency(1.e-4), + expected_exhale_efficiency_small) + npt.assert_equal(mask.exhale_efficiency(4.e-4), + expected_exhale_efficiency_large) + + +@pytest.mark.parametrize( + "diameter, expected_exhale_efficiency", + [ + [0.3e-4, 0.], + [0.7e-4, 0.56711], + [1.e-4, 0.7149], + [4.e-4, 0.8167], + ], +) +def test_measuredmask_exhale(diameter, expected_exhale_efficiency): + mask = models.MeasuredMask(η_inhale=0.3) + npt.assert_almost_equal(mask.exhale_efficiency(diameter), + expected_exhale_efficiency) + From 48f0c395819f9738b4ad146b89c2d47a6a3d6aa8 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 26 May 2021 23:39:33 +0200 Subject: [PATCH 03/83] New method inhale_efficiency in _MaskBase; using _MaskBase everywhere needed (also for types); new MeasuredMask with different exhale_efficiency functions; new mask types --- cara/apps/calculator/model_generator.py | 4 +- cara/apps/expert.py | 14 +++--- cara/models.py | 58 +++++++++++++++++++++---- 3 files changed, 58 insertions(+), 18 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index fb34392d..1ec76a77 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -269,10 +269,10 @@ class FormData: else: return ventilation - def mask(self) -> models.Mask: + def mask(self) -> models._MaskBase: # Initializes the mask type if mask wearing is "continuous", otherwise instantiates the mask attribute as # the "No mask"-mask - mask = models.Mask.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] + mask = models._MaskBase.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] return mask def infected_population(self) -> models.InfectedPopulation: diff --git a/cara/apps/expert.py b/cara/apps/expert.py index 31903b79..a0f273e1 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -390,10 +390,10 @@ class ModelWidgets(View): def _build_mask(self, node): mask = node.dcs_instance() - for name, mask_ in models.Mask.types.items(): + for name, mask_ in models._MaskBase.types.items(): if mask == mask_: break - mask_choice = widgets.Select(options=list(models.Mask.types.keys()), value=name) + mask_choice = widgets.Select(options=list(models._MaskBase.types.keys()), value=name) def on_mask_change(change): node.dcs_select(change['new']) @@ -496,7 +496,7 @@ baseline_model = models.ExposureModel( number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((8, 12), (13, 17))), - mask=models.Mask.types['No mask'], + mask=models._MaskBase.types['No mask'], activity=models.Activity.types['Seated'], expiration=models.Expiration.types['Talking'], ), @@ -505,7 +505,7 @@ baseline_model = models.ExposureModel( number=10, presence=models.SpecificInterval(((8, 12), (13, 17))), activity=models.Activity.types['Seated'], - mask=models.Mask.types['No mask'], + mask=models._MaskBase.types['No mask'], ), ) @@ -515,10 +515,10 @@ class CARAStateBuilder(state.StateBuilder): # For example, build_type__VentilationBase is called when dealing with ConcentrationModel # types as it has a ventilation: _VentilationBase field. - def build_type_Mask(self, _: dataclasses.Field): + def build_type__MaskBase(self, _: dataclasses.Field): return state.DataclassStatePredefined( - models.Mask, - choices=models.Mask.types, + models._MaskBase, + choices=models._MaskBase.types, ) def build_type_Virus(self, _: dataclasses.Field): diff --git a/cara/models.py b/cara/models.py index 55989ec5..2987d613 100644 --- a/cara/models.py +++ b/cara/models.py @@ -474,13 +474,20 @@ Virus.types = { @dataclass(frozen=True) class _MaskBase: """ - Represents the filtration of aerosols by a mask, both inward and + Represents the filtration of aerosols by a mask, both inward and outward. The nature of the various air exchange schemes means that it is expected for subclasses of _MaskBase to exist. """ + #: Pre-populated examples of Masks. + types: typing.ClassVar[typing.Dict[str, "_MaskBase"]] + def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: - # Overall efficiency, including the effect of the leaks. + # Overall exhale efficiency, including the effect of the leaks. + raise NotImplementedError("Subclass must implement") + + def inhale_efficiency(self) -> _VectorisedFloat: + # Overall inhale efficiency, including the effect of the leaks. raise NotImplementedError("Subclass must implement") @@ -495,9 +502,6 @@ class Mask(_MaskBase): #: Filtration efficiency of masks when inhaling. η_inhale: _VectorisedFloat - #: Pre-populated examples of Masks. - types: typing.ClassVar[typing.Dict[str, "Mask"]] - def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: # Overall efficiency with the effect of the leaks for aerosol emission # Gammaitoni et al (1997). Diameter is in cm. @@ -507,8 +511,38 @@ class Mask(_MaskBase): eta_out = self.η_exhale * (1 - self.η_leaks) return eta_out + def inhale_efficiency(self) -> _VectorisedFloat: + # Overall inhale efficiency, including the effect of the leaks. + return self.η_inhale -Mask.types = { + +@dataclass(frozen=True) +class MeasuredMask(_MaskBase): + #: Filtration efficiency of masks when inhaling. + η_inhale: _VectorisedFloat + + def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: + # See CERN-OPEN-2021-004 (doi: 10.17181/CERN.1GDQ.5Y75), and Ref. + # therein (Asadi 2020). + # Obtained from measurements of filtration efficiency and of + # the leakage through the sides. + # Diameter is in cm. + if diameter < 0.5e-4: + eta_out = 0. + elif diameter < 0.94614e-4: + eta_out = 0.5893 * diameter * 1e4 + 0.1546 + elif diameter < 3e-4: + eta_out = 0.0509 * diameter * 1e4 + 0.664 + else: + eta_out = 0.8167 + return eta_out + + def inhale_efficiency(self) -> _VectorisedFloat: + # Overall inhale efficiency, including the effect of the leaks. + return self.η_inhale + + +_MaskBase.types = { 'No mask': Mask(0, 0, 0), 'Type I': Mask( η_exhale=0.95, @@ -520,6 +554,12 @@ Mask.types = { η_leaks=0.15, # (same outward effect as type 1 - Asadi 2020) η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) ), + 'Type I measured': MeasuredMask( + η_inhale=0.3, # (Browen 2010) + ), + 'FFP2 measured': MeasuredMask( + η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) + ), } @@ -531,7 +571,7 @@ class Expiration: #: Pre-populated examples of Expiration. types: typing.ClassVar[typing.Dict[str, "Expiration"]] - def aerosols(self, mask: Mask): + def aerosols(self, mask: _MaskBase): def volume(diameter): return (4 * np.pi * (diameter/2)**3) / 3 total = 0 @@ -583,7 +623,7 @@ class Population: presence: Interval #: The kind of mask being worn by the people. - mask: Mask + mask: _MaskBase #: The physical activity being carried out by the people. activity: Activity @@ -806,7 +846,7 @@ class ExposureModel: inf_aero = ( self.exposed.activity.inhalation_rate * - (1 - self.exposed.mask.η_inhale) * + (1 - self.exposed.mask.inhale_efficiency()) * exposure ) From fec6a71c62d6cc0ba53627fa0e415835c707c7ef Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 May 2021 09:20:28 +0200 Subject: [PATCH 04/83] Correction in docstring of _MaskBase --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 2987d613..620e6ca4 100644 --- a/cara/models.py +++ b/cara/models.py @@ -476,7 +476,7 @@ class _MaskBase: """ Represents the filtration of aerosols by a mask, both inward and outward. - The nature of the various air exchange schemes means that it is expected + The nature of the various mask models means that it is expected for subclasses of _MaskBase to exist. """ #: Pre-populated examples of Masks. From beb2fd737175d146967cd0f2c64d87dca894c8a9 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 May 2021 12:33:12 +0200 Subject: [PATCH 05/83] Adding tests for MultipleExpiration --- cara/tests/test_expiration.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 cara/tests/test_expiration.py 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. + ) From a8081d8b0b74a6ad134752d1b92e2620e2aa9909 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 May 2021 12:34:14 +0200 Subject: [PATCH 06/83] Adapting tests for model_generator --- cara/tests/apps/calculator/test_model_generator.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) 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): From af730ccd7c249042a3b5c42924acf05c74dccf99 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 May 2021 13:40:46 +0200 Subject: [PATCH 07/83] Introducing _ExpirationBase and MultipleExpiration classes; adapting tests and model_generator accordingly (removing now obsolete expiration_blend function) --- cara/apps/calculator/model_generator.py | 32 +++------------- cara/models.py | 51 ++++++++++++++++++++++--- 2 files changed, 50 insertions(+), 33 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 1ec76a77..0d7506ae 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 620e6ca4..bb3ece66 100644 --- a/cara/models.py +++ b/cara/models.py @@ -564,13 +564,28 @@ _MaskBase.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: _MaskBase): + # total volume of aerosols expired (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. + """ 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: _MaskBase): def volume(diameter): return (4 * np.pi * (diameter/2)**3) / 3 @@ -582,7 +597,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: _MaskBase): + 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)), @@ -638,7 +677,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: """ From 213dd2fb929d848f0c76d8d88cb4c27df996ff78 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 May 2021 13:44:03 +0200 Subject: [PATCH 08/83] Modifying eta_inhale of Type I measured mask according to CERN report --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 620e6ca4..f95e28ed 100644 --- a/cara/models.py +++ b/cara/models.py @@ -555,7 +555,7 @@ _MaskBase.types = { η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) ), 'Type I measured': MeasuredMask( - η_inhale=0.3, # (Browen 2010) + η_inhale=0.5, # (CERN-OPEN-2021-004) ), 'FFP2 measured': MeasuredMask( η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) From 22539e2f9fbb3e77cbbb7a845de90b4030dabff8 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 May 2021 18:17:59 +0200 Subject: [PATCH 09/83] Adding test on BLO expiration --- cara/tests/test_expiration.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py index 7638d613..5cfe94b9 100644 --- a/cara/tests/test_expiration.py +++ b/cara/tests/test_expiration.py @@ -23,7 +23,15 @@ def test_multiple(): 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. - ) + mask = models.Mask.types['No mask'] + assert e.aerosols(mask) == e1.aerosols(mask)/3. + 2*e2.aerosols(mask)/3. + + +def test_multiple_BLO(): + weights = (1., 1.) + e1 = models.ExpirationBLO((1., 0., 0.)) + e2 = models.ExpirationBLO((4., 5., 5.)) + e_expected = models.ExpirationBLO((2.5, 2.5, 2.5)) + e = models.MultipleExpiration([e1, e2], weights) + mask = models.Mask.types['Type I'] + assert e_expected.aerosols(mask) == e.aerosols(mask) From 23510c8f1a1a1d56f5d16042bdb7a1adcb2525f9 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 May 2021 18:26:37 +0200 Subject: [PATCH 10/83] Adding BLO model for expiration (ExpirationBLO class) --- cara/models.py | 51 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/cara/models.py b/cara/models.py index a7c8d334..d491e0da 100644 --- a/cara/models.py +++ b/cara/models.py @@ -37,6 +37,7 @@ import typing import numpy as np from scipy.interpolate import interp1d +from scipy.integrate import quad if not typing.TYPE_CHECKING: from memoization import cached @@ -621,12 +622,62 @@ class MultipleExpiration(_ExpirationBase): ]).sum(axis=0) +@dataclass(frozen=True) +class ExpirationBLO(_ExpirationBase): + """ + BLO model for the expiration (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). + All diameters are in cm. + """ + #: factors assigned to resp. the B, L and O modes. Depends on the + # kind of expiratory activity (e.g. breathing, speaking, singing, + # or shouting). + BLO_factors: typing.Tuple[float, float, float] + + def aerosols(self, mask: _MaskBase): + + def volume(diameter): + return (np.pi * diameter**3) / 6. + + def _Bmode(diameter: float) -> float: + # B-mode (see ref. above). + d = diameter * 1e4 # microns + return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.26)) * + np.exp(-1 * (np.log(d) - 1.0) ** 2 / (2 * 0.26 ** 2))) + + def _Lmode(diameter: float) -> float: + # L-mode (see ref. above). + d = diameter * 1e4 # microns + return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.5)) * + np.exp(-1 * (np.log(d) - 1.4) ** 2 / (2 * 0.5 ** 2))) + + def _Omode(diameter: float) -> float: + # O-mode (see ref. above). + d = diameter * 1e4 # microns + return ( (1 / d) * (0.001 / (np.sqrt(2 * np.pi) * 0.56)) * + np.exp(-1 * (np.log(d) - 4.98) ** 2 / (2 * 0.56 ** 2))) + + def integrand(diameter: float) -> float: + return (self.BLO_factors[0] * _Bmode(diameter) + + self.BLO_factors[1] * _Lmode(diameter) + + self.BLO_factors[2] * _Omode(diameter) + ) * volume(diameter) * (1 - mask.exhale_efficiency(diameter)) + + return quad(integrand, 0.1, 30)[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)), 'Unmodulated Vocalization': Expiration((0.751, 0.139, 0.0139, 0.059)), 'Superspreading event': Expiration((np.inf, np.inf, np.inf, np.inf)), + 'Breathing BLO': ExpirationBLO((1., 0., 0.)), + 'Talking BLO': ExpirationBLO((1., 1., 1.)), + 'Shouting BLO': ExpirationBLO((5., 5., 5.)), + 'Singing BLO': ExpirationBLO((5., 5., 5.)), } From a9cd36df7529213241d0b55d57bf32cc11c99870 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 May 2021 18:35:54 +0200 Subject: [PATCH 11/83] Minor docstring fix in ExpirationBLO --- cara/models.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cara/models.py b/cara/models.py index d491e0da..2bdade1f 100644 --- a/cara/models.py +++ b/cara/models.py @@ -631,9 +631,9 @@ class ExpirationBLO(_ExpirationBase): https://doi.org/10.1016/j.jaerosci.2011.07.009). All diameters are in cm. """ - #: factors assigned to resp. the B, L and O modes. Depends on the - # kind of expiratory activity (e.g. breathing, speaking, singing, - # or shouting). + #: 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). BLO_factors: typing.Tuple[float, float, float] def aerosols(self, mask: _MaskBase): From 6e49ed4f0213c32a25179593673f863b4ae00bfd Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 28 May 2021 06:56:13 +0200 Subject: [PATCH 12/83] Improving docstrings in expiration classes --- cara/models.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/cara/models.py b/cara/models.py index a7c8d334..a72f1b48 100644 --- a/cara/models.py +++ b/cara/models.py @@ -573,7 +573,7 @@ class _ExpirationBase: types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]] def aerosols(self, mask: _MaskBase): - # total volume of aerosols expired (cm^3). + # total volume of aerosols expired per volume of air (mL/cm^3). raise NotImplementedError("Subclass must implement") @@ -581,7 +581,10 @@ class _ExpirationBase: class Expiration(_ExpirationBase): """ Simple model based on four different sizes of particles emitted, - with different ejection factors. + 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. From cd3f9057f9cc8d35cd1b0712a8e9afd13e4bce63 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 28 May 2021 07:08:23 +0200 Subject: [PATCH 13/83] Dealing better with the units in ExpirationBLO.aerosols --- cara/models.py | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/cara/models.py b/cara/models.py index c1bee285..496b86e2 100644 --- a/cara/models.py +++ b/cara/models.py @@ -632,7 +632,7 @@ class ExpirationBLO(_ExpirationBase): 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). - All diameters are in cm. + Here all diameters (d) are in microns. """ #: factors assigned to resp. the B, L and O modes. They are # charateristics of the kind of expiratory activity (e.g. breathing, @@ -640,33 +640,30 @@ class ExpirationBLO(_ExpirationBase): BLO_factors: typing.Tuple[float, float, float] def aerosols(self, mask: _MaskBase): + # result is in mL.cm^-3 + def volume(d): + return (np.pi * d**3) / 6. - def volume(diameter): - return (np.pi * diameter**3) / 6. - - def _Bmode(diameter: float) -> float: + def _Bmode(d: float) -> float: # B-mode (see ref. above). - d = diameter * 1e4 # microns return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.26)) * np.exp(-1 * (np.log(d) - 1.0) ** 2 / (2 * 0.26 ** 2))) - def _Lmode(diameter: float) -> float: + def _Lmode(d: float) -> float: # L-mode (see ref. above). - d = diameter * 1e4 # microns return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.5)) * np.exp(-1 * (np.log(d) - 1.4) ** 2 / (2 * 0.5 ** 2))) - def _Omode(diameter: float) -> float: + def _Omode(d: float) -> float: # O-mode (see ref. above). - d = diameter * 1e4 # microns return ( (1 / d) * (0.001 / (np.sqrt(2 * np.pi) * 0.56)) * np.exp(-1 * (np.log(d) - 4.98) ** 2 / (2 * 0.56 ** 2))) - def integrand(diameter: float) -> float: - return (self.BLO_factors[0] * _Bmode(diameter) + - self.BLO_factors[1] * _Lmode(diameter) + - self.BLO_factors[2] * _Omode(diameter) - ) * volume(diameter) * (1 - mask.exhale_efficiency(diameter)) + def integrand(d: float) -> float: + return (self.BLO_factors[0] * _Bmode(d) + + self.BLO_factors[1] * _Lmode(d) + + self.BLO_factors[2] * _Omode(d) + ) * volume(d*1e-4) * (1 - mask.exhale_efficiency(d*1e-4)) return quad(integrand, 0.1, 30)[0] From 9322c27af93718df6938672d585590ae4e5d583c Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 26 May 2021 18:57:56 +0200 Subject: [PATCH 14/83] Introducing _MaskBase class, of which Mask is a subclass --- cara/models.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/cara/models.py b/cara/models.py index c2500d68..f5e80d1a 100644 --- a/cara/models.py +++ b/cara/models.py @@ -472,8 +472,21 @@ Virus.types = { @dataclass(frozen=True) -class Mask: - #: Filtration efficiency. (In %/100) +class _MaskBase: + """ + Represents the filtration of aerosols by a mask, both inward and + outward. + The nature of the various air exchange schemes means that it is expected + for subclasses of _MaskBase to exist. + """ + def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: + # Overall efficiency, including the effect of the leaks. + raise NotImplementedError("Subclass must implement") + + +@dataclass(frozen=True) +class Mask(_MaskBase): + #: Filtration efficiency. η_exhale: _VectorisedFloat #: Leakage through side of masks. From 78cdb22798f7d5cad8445dbba2bac01df4693d90 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 26 May 2021 23:37:10 +0200 Subject: [PATCH 15/83] Adding tests on Mask (new classes and new methods) --- cara/tests/models/test_mask.py | 58 ++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 cara/tests/models/test_mask.py diff --git a/cara/tests/models/test_mask.py b/cara/tests/models/test_mask.py new file mode 100644 index 00000000..f2069cbd --- /dev/null +++ b/cara/tests/models/test_mask.py @@ -0,0 +1,58 @@ + +import dataclasses + +import numpy as np +import numpy.testing as npt +import pytest + +from cara import models + +@pytest.mark.parametrize( + "η_inhale, expected_inhale_efficiency", + [ + [0.5, 0.5], + [np.array([0.3, 0.5]), np.array([0.3, 0.5])], + ], +) +def test_masks_inhale(η_inhale, expected_inhale_efficiency): + mask = models.Mask(η_inhale=η_inhale,η_exhale=0.95,η_leaks=0.15) + measuredmask = models.MeasuredMask(η_inhale=η_inhale) + npt.assert_equal(mask.inhale_efficiency(), + expected_inhale_efficiency) + npt.assert_equal(measuredmask.inhale_efficiency(), + expected_inhale_efficiency) + + +@pytest.mark.parametrize( + "η_exhale, η_leaks, expected_exhale_efficiency_small, expected_exhale_efficiency_large", + [ + [0.95, 0.15, 0., 0.8075], + [np.array([0.95, 1.]), 0.15, np.zeros(2), np.array([0.8075, 0.85])], + [0.95, np.array([0.15, 0.]), np.zeros(2), np.array([0.8075, 0.95])], + [np.array([0.95, 1.]), np.array([0.15, 0.]), np.zeros(2), np.array([0.8075, 1.])], + ], +) +def test_mask_exhale(η_exhale, η_leaks, expected_exhale_efficiency_small, + expected_exhale_efficiency_large): + mask = models.Mask(η_inhale=0.3,η_exhale=η_exhale,η_leaks=η_leaks) + # we test one small and one large diameter (resp. 1 and 4 microns) + npt.assert_equal(mask.exhale_efficiency(1.e-4), + expected_exhale_efficiency_small) + npt.assert_equal(mask.exhale_efficiency(4.e-4), + expected_exhale_efficiency_large) + + +@pytest.mark.parametrize( + "diameter, expected_exhale_efficiency", + [ + [0.3e-4, 0.], + [0.7e-4, 0.56711], + [1.e-4, 0.7149], + [4.e-4, 0.8167], + ], +) +def test_measuredmask_exhale(diameter, expected_exhale_efficiency): + mask = models.MeasuredMask(η_inhale=0.3) + npt.assert_almost_equal(mask.exhale_efficiency(diameter), + expected_exhale_efficiency) + From adc11cb527e2a033ad49b8496136d5c0e85b0b9f Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 26 May 2021 23:39:33 +0200 Subject: [PATCH 16/83] New method inhale_efficiency in _MaskBase; using _MaskBase everywhere needed (also for types); new MeasuredMask with different exhale_efficiency functions; new mask types --- cara/apps/calculator/model_generator.py | 4 +- cara/apps/expert.py | 14 +++--- cara/models.py | 58 +++++++++++++++++++++---- 3 files changed, 58 insertions(+), 18 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index ad65795e..0d7506ae 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -269,10 +269,10 @@ class FormData: else: return ventilation - def mask(self) -> models.Mask: + def mask(self) -> models._MaskBase: # Initializes the mask type if mask wearing is "continuous", otherwise instantiates the mask attribute as # the "No mask"-mask - mask = models.Mask.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] + mask = models._MaskBase.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] return mask def infected_population(self) -> models.InfectedPopulation: diff --git a/cara/apps/expert.py b/cara/apps/expert.py index 31903b79..a0f273e1 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -390,10 +390,10 @@ class ModelWidgets(View): def _build_mask(self, node): mask = node.dcs_instance() - for name, mask_ in models.Mask.types.items(): + for name, mask_ in models._MaskBase.types.items(): if mask == mask_: break - mask_choice = widgets.Select(options=list(models.Mask.types.keys()), value=name) + mask_choice = widgets.Select(options=list(models._MaskBase.types.keys()), value=name) def on_mask_change(change): node.dcs_select(change['new']) @@ -496,7 +496,7 @@ baseline_model = models.ExposureModel( number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((8, 12), (13, 17))), - mask=models.Mask.types['No mask'], + mask=models._MaskBase.types['No mask'], activity=models.Activity.types['Seated'], expiration=models.Expiration.types['Talking'], ), @@ -505,7 +505,7 @@ baseline_model = models.ExposureModel( number=10, presence=models.SpecificInterval(((8, 12), (13, 17))), activity=models.Activity.types['Seated'], - mask=models.Mask.types['No mask'], + mask=models._MaskBase.types['No mask'], ), ) @@ -515,10 +515,10 @@ class CARAStateBuilder(state.StateBuilder): # For example, build_type__VentilationBase is called when dealing with ConcentrationModel # types as it has a ventilation: _VentilationBase field. - def build_type_Mask(self, _: dataclasses.Field): + def build_type__MaskBase(self, _: dataclasses.Field): return state.DataclassStatePredefined( - models.Mask, - choices=models.Mask.types, + models._MaskBase, + choices=models._MaskBase.types, ) def build_type_Virus(self, _: dataclasses.Field): diff --git a/cara/models.py b/cara/models.py index f5e80d1a..5dbc134e 100644 --- a/cara/models.py +++ b/cara/models.py @@ -474,13 +474,20 @@ Virus.types = { @dataclass(frozen=True) class _MaskBase: """ - Represents the filtration of aerosols by a mask, both inward and + Represents the filtration of aerosols by a mask, both inward and outward. The nature of the various air exchange schemes means that it is expected for subclasses of _MaskBase to exist. """ + #: Pre-populated examples of Masks. + types: typing.ClassVar[typing.Dict[str, "_MaskBase"]] + def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: - # Overall efficiency, including the effect of the leaks. + # Overall exhale efficiency, including the effect of the leaks. + raise NotImplementedError("Subclass must implement") + + def inhale_efficiency(self) -> _VectorisedFloat: + # Overall inhale efficiency, including the effect of the leaks. raise NotImplementedError("Subclass must implement") @@ -495,9 +502,6 @@ class Mask(_MaskBase): #: Filtration efficiency of masks when inhaling. η_inhale: _VectorisedFloat - #: Pre-populated examples of Masks. - types: typing.ClassVar[typing.Dict[str, "Mask"]] - def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: # Overall efficiency with the effect of the leaks for aerosol emission # Gammaitoni et al (1997). Diameter is in cm. @@ -507,8 +511,38 @@ class Mask(_MaskBase): eta_out = self.η_exhale * (1 - self.η_leaks) return eta_out + def inhale_efficiency(self) -> _VectorisedFloat: + # Overall inhale efficiency, including the effect of the leaks. + return self.η_inhale -Mask.types = { + +@dataclass(frozen=True) +class MeasuredMask(_MaskBase): + #: Filtration efficiency of masks when inhaling. + η_inhale: _VectorisedFloat + + def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: + # See CERN-OPEN-2021-004 (doi: 10.17181/CERN.1GDQ.5Y75), and Ref. + # therein (Asadi 2020). + # Obtained from measurements of filtration efficiency and of + # the leakage through the sides. + # Diameter is in cm. + if diameter < 0.5e-4: + eta_out = 0. + elif diameter < 0.94614e-4: + eta_out = 0.5893 * diameter * 1e4 + 0.1546 + elif diameter < 3e-4: + eta_out = 0.0509 * diameter * 1e4 + 0.664 + else: + eta_out = 0.8167 + return eta_out + + def inhale_efficiency(self) -> _VectorisedFloat: + # Overall inhale efficiency, including the effect of the leaks. + return self.η_inhale + + +_MaskBase.types = { 'No mask': Mask(0, 0, 0), 'Type I': Mask( η_exhale=0.95, @@ -520,6 +554,12 @@ Mask.types = { η_leaks=0.15, # (same outward effect as type 1 - Asadi 2020) η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) ), + 'Type I measured': MeasuredMask( + η_inhale=0.3, # (Browen 2010) + ), + 'FFP2 measured': MeasuredMask( + η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) + ), } @@ -549,7 +589,7 @@ class Expiration(_ExpirationBase): ejection_factor: typing.Tuple[float, ...] particle_sizes: typing.Tuple[float, ...] = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4) # In cm. - def aerosols(self, mask: Mask): + def aerosols(self, mask: _MaskBase): def volume(diameter): return (4 * np.pi * (diameter/2)**3) / 3 total = 0 @@ -625,7 +665,7 @@ class Population: presence: Interval #: The kind of mask being worn by the people. - mask: Mask + mask: _MaskBase #: The physical activity being carried out by the people. activity: Activity @@ -848,7 +888,7 @@ class ExposureModel: inf_aero = ( self.exposed.activity.inhalation_rate * - (1 - self.exposed.mask.η_inhale) * + (1 - self.exposed.mask.inhale_efficiency()) * exposure ) From 4aa819bea475f2547782e81ad52bf0ada7270c1b Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 May 2021 09:20:28 +0200 Subject: [PATCH 17/83] Correction in docstring of _MaskBase --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 5dbc134e..087704e8 100644 --- a/cara/models.py +++ b/cara/models.py @@ -476,7 +476,7 @@ class _MaskBase: """ Represents the filtration of aerosols by a mask, both inward and outward. - The nature of the various air exchange schemes means that it is expected + The nature of the various mask models means that it is expected for subclasses of _MaskBase to exist. """ #: Pre-populated examples of Masks. From 95640cb5be37aeeb8fb0aebd1e6dcb638c3ac37f Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 May 2021 13:44:03 +0200 Subject: [PATCH 18/83] Modifying eta_inhale of Type I measured mask according to CERN report --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 087704e8..6efa29db 100644 --- a/cara/models.py +++ b/cara/models.py @@ -555,7 +555,7 @@ _MaskBase.types = { η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) ), 'Type I measured': MeasuredMask( - η_inhale=0.3, # (Browen 2010) + η_inhale=0.5, # (CERN-OPEN-2021-004) ), 'FFP2 measured': MeasuredMask( η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) From de5e96fd0f86898f8cb2934be7db7e4315cc06fa Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 May 2021 13:40:46 +0200 Subject: [PATCH 19/83] Introducing _ExpirationBase and MultipleExpiration classes; adapting tests and model_generator accordingly (removing now obsolete expiration_blend function) --- cara/models.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cara/models.py b/cara/models.py index 6efa29db..50eaaec2 100644 --- a/cara/models.py +++ b/cara/models.py @@ -572,8 +572,8 @@ class _ExpirationBase: #: 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). + def aerosols(self, mask: _MaskBase): + # total volume of aerosols expired (cm^3). raise NotImplementedError("Subclass must implement") @@ -617,7 +617,7 @@ class MultipleExpiration(_ExpirationBase): raise ValueError("expirations and weigths should contain the" "same number of elements") - def aerosols(self, mask: Mask): + def aerosols(self, mask: _MaskBase): return np.array([ weight * expiration.aerosols(mask) / sum(self.weights) for weight,expiration in zip(self.weights,self.expirations) From b04dfbad01806cf18dc742510bf68e6d965afb2b Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 07:56:06 +0200 Subject: [PATCH 20/83] Fixing previous merge --- cara/models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cara/models.py b/cara/models.py index 6efa29db..a72f1b48 100644 --- a/cara/models.py +++ b/cara/models.py @@ -572,7 +572,7 @@ class _ExpirationBase: #: Pre-populated examples of Masks. types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]] - def aerosols(self, mask: Mask): + def aerosols(self, mask: _MaskBase): # total volume of aerosols expired per volume of air (mL/cm^3). raise NotImplementedError("Subclass must implement") @@ -617,7 +617,7 @@ class MultipleExpiration(_ExpirationBase): raise ValueError("expirations and weigths should contain the" "same number of elements") - def aerosols(self, mask: Mask): + def aerosols(self, mask: _MaskBase): return np.array([ weight * expiration.aerosols(mask) / sum(self.weights) for weight,expiration in zip(self.weights,self.expirations) From 947bd013e0c07131eb8ddc33c4c876f2bbead2e8 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 19:25:39 +0200 Subject: [PATCH 21/83] Modifying tests on mask to comply with a single Mask class --- cara/tests/models/test_mask.py | 43 +++++++--------------------------- 1 file changed, 9 insertions(+), 34 deletions(-) diff --git a/cara/tests/models/test_mask.py b/cara/tests/models/test_mask.py index f2069cbd..7f7cb242 100644 --- a/cara/tests/models/test_mask.py +++ b/cara/tests/models/test_mask.py @@ -1,6 +1,3 @@ - -import dataclasses - import numpy as np import numpy.testing as npt import pytest @@ -15,44 +12,22 @@ from cara import models ], ) def test_masks_inhale(η_inhale, expected_inhale_efficiency): - mask = models.Mask(η_inhale=η_inhale,η_exhale=0.95,η_leaks=0.15) - measuredmask = models.MeasuredMask(η_inhale=η_inhale) + mask = models.Mask(η_inhale=η_inhale) npt.assert_equal(mask.inhale_efficiency(), expected_inhale_efficiency) - npt.assert_equal(measuredmask.inhale_efficiency(), - expected_inhale_efficiency) @pytest.mark.parametrize( - "η_exhale, η_leaks, expected_exhale_efficiency_small, expected_exhale_efficiency_large", + "diameter, factor_exhale, expected_exhale_efficiency", [ - [0.95, 0.15, 0., 0.8075], - [np.array([0.95, 1.]), 0.15, np.zeros(2), np.array([0.8075, 0.85])], - [0.95, np.array([0.15, 0.]), np.zeros(2), np.array([0.8075, 0.95])], - [np.array([0.95, 1.]), np.array([0.15, 0.]), np.zeros(2), np.array([0.8075, 1.])], + [0.3e-4, 1., 0.], + [0.7e-4, 0.3, 0.56711*0.3], + [1.e-4, 1., 0.7149], + [4.e-4, 0.5, 0.8167*0.5], + [5.e-4, 0., 0.], ], ) -def test_mask_exhale(η_exhale, η_leaks, expected_exhale_efficiency_small, - expected_exhale_efficiency_large): - mask = models.Mask(η_inhale=0.3,η_exhale=η_exhale,η_leaks=η_leaks) - # we test one small and one large diameter (resp. 1 and 4 microns) - npt.assert_equal(mask.exhale_efficiency(1.e-4), - expected_exhale_efficiency_small) - npt.assert_equal(mask.exhale_efficiency(4.e-4), - expected_exhale_efficiency_large) - - -@pytest.mark.parametrize( - "diameter, expected_exhale_efficiency", - [ - [0.3e-4, 0.], - [0.7e-4, 0.56711], - [1.e-4, 0.7149], - [4.e-4, 0.8167], - ], -) -def test_measuredmask_exhale(diameter, expected_exhale_efficiency): - mask = models.MeasuredMask(η_inhale=0.3) +def test_mask_exhale(diameter, factor_exhale, expected_exhale_efficiency): + mask = models.Mask(η_inhale=0.3, factor_exhale=factor_exhale) npt.assert_almost_equal(mask.exhale_efficiency(diameter), expected_exhale_efficiency) - From fa0550233c8369f81158ad6abf1cb4bc8bf12817 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 19:28:51 +0200 Subject: [PATCH 22/83] Removing _MaskBase and old Mask class, which is replaced by MeasuredMask --- cara/models.py | 100 ++++++++++++++----------------------------------- 1 file changed, 28 insertions(+), 72 deletions(-) diff --git a/cara/models.py b/cara/models.py index a72f1b48..736e0992 100644 --- a/cara/models.py +++ b/cara/models.py @@ -472,61 +472,25 @@ Virus.types = { @dataclass(frozen=True) -class _MaskBase: - """ - Represents the filtration of aerosols by a mask, both inward and - outward. - The nature of the various mask models means that it is expected - for subclasses of _MaskBase to exist. - """ +class Mask: + #: Filtration efficiency of masks when inhaling. + η_inhale: _VectorisedFloat + + #: Global factor applied to filtration efficiency of masks when exhaling. + factor_exhale: _VectorisedFloat = 1. + #: Pre-populated examples of Masks. - types: typing.ClassVar[typing.Dict[str, "_MaskBase"]] + types: typing.ClassVar[typing.Dict[str, "Mask"]] def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: - # Overall exhale efficiency, including the effect of the leaks. - raise NotImplementedError("Subclass must implement") - - def inhale_efficiency(self) -> _VectorisedFloat: - # Overall inhale efficiency, including the effect of the leaks. - raise NotImplementedError("Subclass must implement") - - -@dataclass(frozen=True) -class Mask(_MaskBase): - #: Filtration efficiency. - η_exhale: _VectorisedFloat - - #: Leakage through side of masks. - η_leaks: _VectorisedFloat - - #: Filtration efficiency of masks when inhaling. - η_inhale: _VectorisedFloat - - def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: - # Overall efficiency with the effect of the leaks for aerosol emission - # Gammaitoni et al (1997). Diameter is in cm. - if diameter < 3e-4: - eta_out = 0. - else: - eta_out = self.η_exhale * (1 - self.η_leaks) - return eta_out - - def inhale_efficiency(self) -> _VectorisedFloat: - # Overall inhale efficiency, including the effect of the leaks. - return self.η_inhale - - -@dataclass(frozen=True) -class MeasuredMask(_MaskBase): - #: Filtration efficiency of masks when inhaling. - η_inhale: _VectorisedFloat - - def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: - # See CERN-OPEN-2021-004 (doi: 10.17181/CERN.1GDQ.5Y75), and Ref. - # therein (Asadi 2020). - # Obtained from measurements of filtration efficiency and of - # the leakage through the sides. - # Diameter is in cm. + """ + Overall exhale efficiency, including the effect of the leaks. + See CERN-OPEN-2021-004 (doi: 10.17181/CERN.1GDQ.5Y75), and Ref. + therein (Asadi 2020). + Obtained from measurements of filtration efficiency and of + the leakage through the sides. + Diameter is in cm. + """ if diameter < 0.5e-4: eta_out = 0. elif diameter < 0.94614e-4: @@ -535,29 +499,21 @@ class MeasuredMask(_MaskBase): eta_out = 0.0509 * diameter * 1e4 + 0.664 else: eta_out = 0.8167 - return eta_out + return eta_out*self.factor_exhale def inhale_efficiency(self) -> _VectorisedFloat: - # Overall inhale efficiency, including the effect of the leaks. + """ + Overall inhale efficiency, including the effect of the leaks. + """ return self.η_inhale -_MaskBase.types = { - 'No mask': Mask(0, 0, 0), +Mask.types = { + 'No mask': Mask(0, 0), 'Type I': Mask( - η_exhale=0.95, - η_leaks=0.15, # (Huang 2007) - η_inhale=0.3, # (Browen 2010) - ), - 'FFP2': Mask( - η_exhale=0.95, # (same outward effect as type 1 - Asadi 2020) - η_leaks=0.15, # (same outward effect as type 1 - Asadi 2020) - η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) - ), - 'Type I measured': MeasuredMask( η_inhale=0.5, # (CERN-OPEN-2021-004) ), - 'FFP2 measured': MeasuredMask( + 'FFP2': Mask( η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) ), } @@ -569,10 +525,10 @@ class _ExpirationBase: Represents the expiration of aerosols by a person. Subclasses of _ExpirationBase represent different models. """ - #: Pre-populated examples of Masks. + #: Pre-populated examples of Expirations. types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]] - def aerosols(self, mask: _MaskBase): + def aerosols(self, mask: Mask): # total volume of aerosols expired per volume of air (mL/cm^3). raise NotImplementedError("Subclass must implement") @@ -589,7 +545,7 @@ class Expiration(_ExpirationBase): ejection_factor: typing.Tuple[float, ...] particle_sizes: typing.Tuple[float, ...] = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4) # In cm. - def aerosols(self, mask: _MaskBase): + def aerosols(self, mask: Mask): def volume(diameter): return (4 * np.pi * (diameter/2)**3) / 3 total = 0 @@ -617,7 +573,7 @@ class MultipleExpiration(_ExpirationBase): raise ValueError("expirations and weigths should contain the" "same number of elements") - def aerosols(self, mask: _MaskBase): + def aerosols(self, mask: Mask): return np.array([ weight * expiration.aerosols(mask) / sum(self.weights) for weight,expiration in zip(self.weights,self.expirations) @@ -665,7 +621,7 @@ class Population: presence: Interval #: The kind of mask being worn by the people. - mask: _MaskBase + mask: Mask #: The physical activity being carried out by the people. activity: Activity From 26075cd53cc22a0aba35185d28c8b49e496e4108 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 19:30:25 +0200 Subject: [PATCH 23/83] Modifying infected populations tests and concentration tests, according to the new Mask class --- cara/tests/models/test_concentration_model.py | 12 +++++------- cara/tests/test_infected_population.py | 9 +++------ 2 files changed, 8 insertions(+), 13 deletions(-) diff --git a/cara/tests/models/test_concentration_model.py b/cara/tests/models/test_concentration_model.py index 48ef0b66..f7230430 100644 --- a/cara/tests/models/test_concentration_model.py +++ b/cara/tests/models/test_concentration_model.py @@ -1,6 +1,7 @@ import re import numpy as np +import numpy.testing as npt import pytest from cara import models @@ -13,8 +14,7 @@ from cara import models {'air_change': np.array([100, 120])}, {'viral_load_in_sputum': np.array([5e8, 1e9])}, {'quantum_infectious_dose': np.array([50, 20])}, - {'η_exhale': np.array([0.92, 0.95])}, - {'η_leaks': np.array([0.15, 0.20])}, + {'factor_exhale': np.array([0.92, 0.95])}, ] ) def test_concentration_model_vectorisation(override_params): @@ -24,8 +24,7 @@ def test_concentration_model_vectorisation(override_params): 'air_change': 100, 'viral_load_in_sputum': 1e9, 'quantum_infectious_dose': 50, - 'η_exhale': 0.95, - 'η_leaks': 0.15, + 'factor_exhale': 0.95, } defaults.update(override_params) @@ -37,8 +36,7 @@ def test_concentration_model_vectorisation(override_params): number=1, presence=always, mask=models.Mask( - η_exhale=defaults['η_exhale'], - η_leaks=defaults['η_leaks'], + factor_exhale=defaults['factor_exhale'], η_inhale=0.3, ), activity=models.Activity( @@ -128,4 +126,4 @@ def test_integrated_concentration(simple_conc_model): c2 = simple_conc_model.integrated_concentration(0, 1) c3 = simple_conc_model.integrated_concentration(1, 2) assert c1 != 0 - assert c1 == c2 + c3 + npt.assert_almost_equal(c1, c2 + c3, decimal=15) diff --git a/cara/tests/test_infected_population.py b/cara/tests/test_infected_population.py index e87ddb20..8ef94729 100644 --- a/cara/tests/test_infected_population.py +++ b/cara/tests/test_infected_population.py @@ -8,8 +8,7 @@ import cara.models "override_params", [ {'viral_load_in_sputum': np.array([5e8, 1e9])}, {'quantum_infectious_dose': np.array([50, 20])}, - {'η_exhale': np.array([0.92, 0.95])}, - {'η_leaks': np.array([0.15, 0.20])}, + {'factor_exhale': np.array([0.92, 0.95])}, {'exhalation_rate': np.array([0.75, 0.81])}, ] ) @@ -17,8 +16,7 @@ def test_infected_population_vectorisation(override_params): defaults = { 'viral_load_in_sputum': 1e9, 'quantum_infectious_dose': 50, - 'η_exhale': 0.95, - 'η_leaks': 0.15, + 'factor_exhale': 0.95, 'exhalation_rate': 0.75, } defaults.update(override_params) @@ -28,8 +26,7 @@ def test_infected_population_vectorisation(override_params): number=1, presence=office_hours, mask=cara.models.Mask( - η_exhale=defaults['η_exhale'], - η_leaks=defaults['η_leaks'], + factor_exhale=defaults['factor_exhale'], η_inhale=0.3, ), activity=cara.models.Activity( From 0087bff41ed69b9dcb561189978b6e1d623a6e1e Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 19:31:13 +0200 Subject: [PATCH 24/83] Propagating the change in Mask class to the calculator and expoert apps --- cara/apps/calculator/model_generator.py | 4 ++-- cara/apps/expert.py | 14 +++++++------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 0d7506ae..ad65795e 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -269,10 +269,10 @@ class FormData: else: return ventilation - def mask(self) -> models._MaskBase: + def mask(self) -> models.Mask: # Initializes the mask type if mask wearing is "continuous", otherwise instantiates the mask attribute as # the "No mask"-mask - mask = models._MaskBase.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] + mask = models.Mask.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] return mask def infected_population(self) -> models.InfectedPopulation: diff --git a/cara/apps/expert.py b/cara/apps/expert.py index a0f273e1..31903b79 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -390,10 +390,10 @@ class ModelWidgets(View): def _build_mask(self, node): mask = node.dcs_instance() - for name, mask_ in models._MaskBase.types.items(): + for name, mask_ in models.Mask.types.items(): if mask == mask_: break - mask_choice = widgets.Select(options=list(models._MaskBase.types.keys()), value=name) + mask_choice = widgets.Select(options=list(models.Mask.types.keys()), value=name) def on_mask_change(change): node.dcs_select(change['new']) @@ -496,7 +496,7 @@ baseline_model = models.ExposureModel( number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((8, 12), (13, 17))), - mask=models._MaskBase.types['No mask'], + mask=models.Mask.types['No mask'], activity=models.Activity.types['Seated'], expiration=models.Expiration.types['Talking'], ), @@ -505,7 +505,7 @@ baseline_model = models.ExposureModel( number=10, presence=models.SpecificInterval(((8, 12), (13, 17))), activity=models.Activity.types['Seated'], - mask=models._MaskBase.types['No mask'], + mask=models.Mask.types['No mask'], ), ) @@ -515,10 +515,10 @@ class CARAStateBuilder(state.StateBuilder): # For example, build_type__VentilationBase is called when dealing with ConcentrationModel # types as it has a ventilation: _VentilationBase field. - def build_type__MaskBase(self, _: dataclasses.Field): + def build_type_Mask(self, _: dataclasses.Field): return state.DataclassStatePredefined( - models._MaskBase, - choices=models._MaskBase.types, + models.Mask, + choices=models.Mask.types, ) def build_type_Virus(self, _: dataclasses.Field): From d7bd53c288416a41c229ed8ddefd87523733a7d0 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 19:31:56 +0200 Subject: [PATCH 25/83] Starting to modify exposure tests according to new Mask class --- cara/tests/models/test_exposure_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 464f25a9..c1e9090d 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -43,7 +43,7 @@ populations = [ ), # A population with some array component for η_inhale. models.Population( - 10, halftime, models.Mask(0.95, 0.15, np.array([0.3, 0.35])), + 10, halftime, models.Mask(np.array([0.3, 0.35])), models.Activity.types['Standing'], ), # A population with some array component for inhalation_rate. From ffb85c6b72b39b7fe62e1cdfa420282452960650 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 06:32:17 +0200 Subject: [PATCH 26/83] Fixing the remaining tests, with the new masks --- cara/tests/models/test_exposure_model.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index c1e9090d..009cf092 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -60,10 +60,10 @@ populations = [ np.array([14.4, 14.4]), np.array([99.6803184113, 99.5181053773])], [populations[2], KnownConcentrations(lambda t: 1.2), - np.array([14.4, 14.4]), np.array([99.4146994564, 99.6803184113])], + np.array([14.4, 14.4]), np.array([97.4574432074, 98.3493482895])], [populations[0], KnownConcentrations(lambda t: np.array([1.2, 2.4])), - np.array([14.4, 28.8]), np.array([99.6803184113, 99.9989780368])], + np.array([14.4, 28.8]), np.array([98.3493482895, 99.9727534893])], [populations[1], KnownConcentrations(lambda t: np.array([1.2, 2.4])), np.array([14.4, 28.8]), np.array([99.6803184113, 99.9976777757])], @@ -123,22 +123,22 @@ def conc_model(): models.InfectedPopulation( number=1, presence=interesting_times, - mask=models.Mask.types['Type I'], + mask=models.Mask.types['No mask'], activity=models.Activity.types['Seated'], virus=models.Virus.types['SARS_CoV_2'], - expiration=models.Expiration.types['Breathing'], + expiration=models.Expiration.types['Superspreading event'], ) ) # expected quanta were computed with a trapezoidal integration, using # a mesh of 10'000 pts per exposed presence interval. @pytest.mark.parametrize("exposed_time_interval, expected_quanta", [ - [(0, 1), 0.0055680845], - [(1, 1.01), 6.4960491e-05], - [(1.01, 1.02), 6.3187723e-05], - [(12, 12.01), 1.9307359e-06], - [(12, 24), 0.079347465], - [(0, 24), 0.086122050], + [(0, 1), 5.4869151], + [(1, 1.01), 0.064013521], + [(1.01, 1.02), 0.062266596], + [(12, 12.01), 0.0019025904], + [(12, 24), 78.190763], + [(0, 24), 84.866592], ] ) def test_exposure_model_integral_accuracy(exposed_time_interval, From bfdb3e322e1b6f1de291fb64adeb4f1591974a10 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 06:55:07 +0200 Subject: [PATCH 27/83] Mask exhale_efficiency use now a diameter in microns --- cara/models.py | 14 +++++++------- cara/tests/models/test_mask.py | 12 ++++++------ 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/cara/models.py b/cara/models.py index 736e0992..56f663f4 100644 --- a/cara/models.py +++ b/cara/models.py @@ -489,14 +489,14 @@ class Mask: therein (Asadi 2020). Obtained from measurements of filtration efficiency and of the leakage through the sides. - Diameter is in cm. + Diameter is in microns. """ - if diameter < 0.5e-4: + if diameter < 0.5: eta_out = 0. - elif diameter < 0.94614e-4: - eta_out = 0.5893 * diameter * 1e4 + 0.1546 - elif diameter < 3e-4: - eta_out = 0.0509 * diameter * 1e4 + 0.664 + elif diameter < 0.94614: + eta_out = 0.5893 * diameter + 0.1546 + elif diameter < 3.: + eta_out = 0.0509 * diameter + 0.664 else: eta_out = 0.8167 return eta_out*self.factor_exhale @@ -551,7 +551,7 @@ class Expiration(_ExpirationBase): total = 0 for diameter, factor in zip(self.particle_sizes, self.ejection_factor): contribution = (volume(diameter) * factor * - (1 - mask.exhale_efficiency(diameter))) + (1 - mask.exhale_efficiency(diameter*1e4))) total += contribution return total diff --git a/cara/tests/models/test_mask.py b/cara/tests/models/test_mask.py index 7f7cb242..4a71d13c 100644 --- a/cara/tests/models/test_mask.py +++ b/cara/tests/models/test_mask.py @@ -11,7 +11,7 @@ from cara import models [np.array([0.3, 0.5]), np.array([0.3, 0.5])], ], ) -def test_masks_inhale(η_inhale, expected_inhale_efficiency): +def test_mask_inhale(η_inhale, expected_inhale_efficiency): mask = models.Mask(η_inhale=η_inhale) npt.assert_equal(mask.inhale_efficiency(), expected_inhale_efficiency) @@ -20,11 +20,11 @@ def test_masks_inhale(η_inhale, expected_inhale_efficiency): @pytest.mark.parametrize( "diameter, factor_exhale, expected_exhale_efficiency", [ - [0.3e-4, 1., 0.], - [0.7e-4, 0.3, 0.56711*0.3], - [1.e-4, 1., 0.7149], - [4.e-4, 0.5, 0.8167*0.5], - [5.e-4, 0., 0.], + [0.3, 1., 0.], + [0.7, 0.3, 0.56711*0.3], + [1., 1., 0.7149], + [4., 0.5, 0.8167*0.5], + [5., 0., 0.], ], ) def test_mask_exhale(diameter, factor_exhale, expected_exhale_efficiency): From 9325eea09e5d339c103df48dfda2ba773cecda14 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 08:55:55 +0200 Subject: [PATCH 28/83] Fixing docstring --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 56f663f4..88e5493a 100644 --- a/cara/models.py +++ b/cara/models.py @@ -529,7 +529,7 @@ class _ExpirationBase: types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]] def aerosols(self, mask: Mask): - # total volume of aerosols expired per volume of air (mL/cm^3). + # total volume of aerosols expired (cm^3). raise NotImplementedError("Subclass must implement") From 6fd8d6531b768631c68b118bd3fbdf6ad37100a1 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 07:56:06 +0200 Subject: [PATCH 29/83] Fixing previous merge --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 50eaaec2..a72f1b48 100644 --- a/cara/models.py +++ b/cara/models.py @@ -573,7 +573,7 @@ class _ExpirationBase: types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]] def aerosols(self, mask: _MaskBase): - # total volume of aerosols expired (cm^3). + # total volume of aerosols expired per volume of air (mL/cm^3). raise NotImplementedError("Subclass must implement") From 3226bb04e8d71847757fee8bd0f35f89f2eb4bfd Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 19:25:39 +0200 Subject: [PATCH 30/83] Modifying tests on mask to comply with a single Mask class --- cara/tests/models/test_mask.py | 43 +++++++--------------------------- 1 file changed, 9 insertions(+), 34 deletions(-) diff --git a/cara/tests/models/test_mask.py b/cara/tests/models/test_mask.py index f2069cbd..7f7cb242 100644 --- a/cara/tests/models/test_mask.py +++ b/cara/tests/models/test_mask.py @@ -1,6 +1,3 @@ - -import dataclasses - import numpy as np import numpy.testing as npt import pytest @@ -15,44 +12,22 @@ from cara import models ], ) def test_masks_inhale(η_inhale, expected_inhale_efficiency): - mask = models.Mask(η_inhale=η_inhale,η_exhale=0.95,η_leaks=0.15) - measuredmask = models.MeasuredMask(η_inhale=η_inhale) + mask = models.Mask(η_inhale=η_inhale) npt.assert_equal(mask.inhale_efficiency(), expected_inhale_efficiency) - npt.assert_equal(measuredmask.inhale_efficiency(), - expected_inhale_efficiency) @pytest.mark.parametrize( - "η_exhale, η_leaks, expected_exhale_efficiency_small, expected_exhale_efficiency_large", + "diameter, factor_exhale, expected_exhale_efficiency", [ - [0.95, 0.15, 0., 0.8075], - [np.array([0.95, 1.]), 0.15, np.zeros(2), np.array([0.8075, 0.85])], - [0.95, np.array([0.15, 0.]), np.zeros(2), np.array([0.8075, 0.95])], - [np.array([0.95, 1.]), np.array([0.15, 0.]), np.zeros(2), np.array([0.8075, 1.])], + [0.3e-4, 1., 0.], + [0.7e-4, 0.3, 0.56711*0.3], + [1.e-4, 1., 0.7149], + [4.e-4, 0.5, 0.8167*0.5], + [5.e-4, 0., 0.], ], ) -def test_mask_exhale(η_exhale, η_leaks, expected_exhale_efficiency_small, - expected_exhale_efficiency_large): - mask = models.Mask(η_inhale=0.3,η_exhale=η_exhale,η_leaks=η_leaks) - # we test one small and one large diameter (resp. 1 and 4 microns) - npt.assert_equal(mask.exhale_efficiency(1.e-4), - expected_exhale_efficiency_small) - npt.assert_equal(mask.exhale_efficiency(4.e-4), - expected_exhale_efficiency_large) - - -@pytest.mark.parametrize( - "diameter, expected_exhale_efficiency", - [ - [0.3e-4, 0.], - [0.7e-4, 0.56711], - [1.e-4, 0.7149], - [4.e-4, 0.8167], - ], -) -def test_measuredmask_exhale(diameter, expected_exhale_efficiency): - mask = models.MeasuredMask(η_inhale=0.3) +def test_mask_exhale(diameter, factor_exhale, expected_exhale_efficiency): + mask = models.Mask(η_inhale=0.3, factor_exhale=factor_exhale) npt.assert_almost_equal(mask.exhale_efficiency(diameter), expected_exhale_efficiency) - From d5fb86d694a9b7cb33fd6e305ce89b3e765feafc Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 19:28:51 +0200 Subject: [PATCH 31/83] Removing _MaskBase and old Mask class, which is replaced by MeasuredMask --- cara/models.py | 100 ++++++++++++++----------------------------------- 1 file changed, 28 insertions(+), 72 deletions(-) diff --git a/cara/models.py b/cara/models.py index a72f1b48..736e0992 100644 --- a/cara/models.py +++ b/cara/models.py @@ -472,61 +472,25 @@ Virus.types = { @dataclass(frozen=True) -class _MaskBase: - """ - Represents the filtration of aerosols by a mask, both inward and - outward. - The nature of the various mask models means that it is expected - for subclasses of _MaskBase to exist. - """ +class Mask: + #: Filtration efficiency of masks when inhaling. + η_inhale: _VectorisedFloat + + #: Global factor applied to filtration efficiency of masks when exhaling. + factor_exhale: _VectorisedFloat = 1. + #: Pre-populated examples of Masks. - types: typing.ClassVar[typing.Dict[str, "_MaskBase"]] + types: typing.ClassVar[typing.Dict[str, "Mask"]] def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: - # Overall exhale efficiency, including the effect of the leaks. - raise NotImplementedError("Subclass must implement") - - def inhale_efficiency(self) -> _VectorisedFloat: - # Overall inhale efficiency, including the effect of the leaks. - raise NotImplementedError("Subclass must implement") - - -@dataclass(frozen=True) -class Mask(_MaskBase): - #: Filtration efficiency. - η_exhale: _VectorisedFloat - - #: Leakage through side of masks. - η_leaks: _VectorisedFloat - - #: Filtration efficiency of masks when inhaling. - η_inhale: _VectorisedFloat - - def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: - # Overall efficiency with the effect of the leaks for aerosol emission - # Gammaitoni et al (1997). Diameter is in cm. - if diameter < 3e-4: - eta_out = 0. - else: - eta_out = self.η_exhale * (1 - self.η_leaks) - return eta_out - - def inhale_efficiency(self) -> _VectorisedFloat: - # Overall inhale efficiency, including the effect of the leaks. - return self.η_inhale - - -@dataclass(frozen=True) -class MeasuredMask(_MaskBase): - #: Filtration efficiency of masks when inhaling. - η_inhale: _VectorisedFloat - - def exhale_efficiency(self, diameter: float) -> _VectorisedFloat: - # See CERN-OPEN-2021-004 (doi: 10.17181/CERN.1GDQ.5Y75), and Ref. - # therein (Asadi 2020). - # Obtained from measurements of filtration efficiency and of - # the leakage through the sides. - # Diameter is in cm. + """ + Overall exhale efficiency, including the effect of the leaks. + See CERN-OPEN-2021-004 (doi: 10.17181/CERN.1GDQ.5Y75), and Ref. + therein (Asadi 2020). + Obtained from measurements of filtration efficiency and of + the leakage through the sides. + Diameter is in cm. + """ if diameter < 0.5e-4: eta_out = 0. elif diameter < 0.94614e-4: @@ -535,29 +499,21 @@ class MeasuredMask(_MaskBase): eta_out = 0.0509 * diameter * 1e4 + 0.664 else: eta_out = 0.8167 - return eta_out + return eta_out*self.factor_exhale def inhale_efficiency(self) -> _VectorisedFloat: - # Overall inhale efficiency, including the effect of the leaks. + """ + Overall inhale efficiency, including the effect of the leaks. + """ return self.η_inhale -_MaskBase.types = { - 'No mask': Mask(0, 0, 0), +Mask.types = { + 'No mask': Mask(0, 0), 'Type I': Mask( - η_exhale=0.95, - η_leaks=0.15, # (Huang 2007) - η_inhale=0.3, # (Browen 2010) - ), - 'FFP2': Mask( - η_exhale=0.95, # (same outward effect as type 1 - Asadi 2020) - η_leaks=0.15, # (same outward effect as type 1 - Asadi 2020) - η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) - ), - 'Type I measured': MeasuredMask( η_inhale=0.5, # (CERN-OPEN-2021-004) ), - 'FFP2 measured': MeasuredMask( + 'FFP2': Mask( η_inhale=0.865, # (94% penetration efficiency + 8% max inward leakage -> EN 149) ), } @@ -569,10 +525,10 @@ class _ExpirationBase: Represents the expiration of aerosols by a person. Subclasses of _ExpirationBase represent different models. """ - #: Pre-populated examples of Masks. + #: Pre-populated examples of Expirations. types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]] - def aerosols(self, mask: _MaskBase): + def aerosols(self, mask: Mask): # total volume of aerosols expired per volume of air (mL/cm^3). raise NotImplementedError("Subclass must implement") @@ -589,7 +545,7 @@ class Expiration(_ExpirationBase): ejection_factor: typing.Tuple[float, ...] particle_sizes: typing.Tuple[float, ...] = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4) # In cm. - def aerosols(self, mask: _MaskBase): + def aerosols(self, mask: Mask): def volume(diameter): return (4 * np.pi * (diameter/2)**3) / 3 total = 0 @@ -617,7 +573,7 @@ class MultipleExpiration(_ExpirationBase): raise ValueError("expirations and weigths should contain the" "same number of elements") - def aerosols(self, mask: _MaskBase): + def aerosols(self, mask: Mask): return np.array([ weight * expiration.aerosols(mask) / sum(self.weights) for weight,expiration in zip(self.weights,self.expirations) @@ -665,7 +621,7 @@ class Population: presence: Interval #: The kind of mask being worn by the people. - mask: _MaskBase + mask: Mask #: The physical activity being carried out by the people. activity: Activity From 8305144b1d31dbcd1f2b369f5a6277cdc0544ef3 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 19:30:25 +0200 Subject: [PATCH 32/83] Modifying infected populations tests and concentration tests, according to the new Mask class --- cara/tests/models/test_concentration_model.py | 12 +++++------- cara/tests/test_infected_population.py | 9 +++------ 2 files changed, 8 insertions(+), 13 deletions(-) diff --git a/cara/tests/models/test_concentration_model.py b/cara/tests/models/test_concentration_model.py index 48ef0b66..f7230430 100644 --- a/cara/tests/models/test_concentration_model.py +++ b/cara/tests/models/test_concentration_model.py @@ -1,6 +1,7 @@ import re import numpy as np +import numpy.testing as npt import pytest from cara import models @@ -13,8 +14,7 @@ from cara import models {'air_change': np.array([100, 120])}, {'viral_load_in_sputum': np.array([5e8, 1e9])}, {'quantum_infectious_dose': np.array([50, 20])}, - {'η_exhale': np.array([0.92, 0.95])}, - {'η_leaks': np.array([0.15, 0.20])}, + {'factor_exhale': np.array([0.92, 0.95])}, ] ) def test_concentration_model_vectorisation(override_params): @@ -24,8 +24,7 @@ def test_concentration_model_vectorisation(override_params): 'air_change': 100, 'viral_load_in_sputum': 1e9, 'quantum_infectious_dose': 50, - 'η_exhale': 0.95, - 'η_leaks': 0.15, + 'factor_exhale': 0.95, } defaults.update(override_params) @@ -37,8 +36,7 @@ def test_concentration_model_vectorisation(override_params): number=1, presence=always, mask=models.Mask( - η_exhale=defaults['η_exhale'], - η_leaks=defaults['η_leaks'], + factor_exhale=defaults['factor_exhale'], η_inhale=0.3, ), activity=models.Activity( @@ -128,4 +126,4 @@ def test_integrated_concentration(simple_conc_model): c2 = simple_conc_model.integrated_concentration(0, 1) c3 = simple_conc_model.integrated_concentration(1, 2) assert c1 != 0 - assert c1 == c2 + c3 + npt.assert_almost_equal(c1, c2 + c3, decimal=15) diff --git a/cara/tests/test_infected_population.py b/cara/tests/test_infected_population.py index e87ddb20..8ef94729 100644 --- a/cara/tests/test_infected_population.py +++ b/cara/tests/test_infected_population.py @@ -8,8 +8,7 @@ import cara.models "override_params", [ {'viral_load_in_sputum': np.array([5e8, 1e9])}, {'quantum_infectious_dose': np.array([50, 20])}, - {'η_exhale': np.array([0.92, 0.95])}, - {'η_leaks': np.array([0.15, 0.20])}, + {'factor_exhale': np.array([0.92, 0.95])}, {'exhalation_rate': np.array([0.75, 0.81])}, ] ) @@ -17,8 +16,7 @@ def test_infected_population_vectorisation(override_params): defaults = { 'viral_load_in_sputum': 1e9, 'quantum_infectious_dose': 50, - 'η_exhale': 0.95, - 'η_leaks': 0.15, + 'factor_exhale': 0.95, 'exhalation_rate': 0.75, } defaults.update(override_params) @@ -28,8 +26,7 @@ def test_infected_population_vectorisation(override_params): number=1, presence=office_hours, mask=cara.models.Mask( - η_exhale=defaults['η_exhale'], - η_leaks=defaults['η_leaks'], + factor_exhale=defaults['factor_exhale'], η_inhale=0.3, ), activity=cara.models.Activity( From 58e40eff273a637f844b2d291217d5310c65c4bc Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 19:31:13 +0200 Subject: [PATCH 33/83] Propagating the change in Mask class to the calculator and expoert apps --- cara/apps/calculator/model_generator.py | 4 ++-- cara/apps/expert.py | 14 +++++++------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 0d7506ae..ad65795e 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -269,10 +269,10 @@ class FormData: else: return ventilation - def mask(self) -> models._MaskBase: + def mask(self) -> models.Mask: # Initializes the mask type if mask wearing is "continuous", otherwise instantiates the mask attribute as # the "No mask"-mask - mask = models._MaskBase.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] + mask = models.Mask.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] return mask def infected_population(self) -> models.InfectedPopulation: diff --git a/cara/apps/expert.py b/cara/apps/expert.py index a0f273e1..31903b79 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -390,10 +390,10 @@ class ModelWidgets(View): def _build_mask(self, node): mask = node.dcs_instance() - for name, mask_ in models._MaskBase.types.items(): + for name, mask_ in models.Mask.types.items(): if mask == mask_: break - mask_choice = widgets.Select(options=list(models._MaskBase.types.keys()), value=name) + mask_choice = widgets.Select(options=list(models.Mask.types.keys()), value=name) def on_mask_change(change): node.dcs_select(change['new']) @@ -496,7 +496,7 @@ baseline_model = models.ExposureModel( number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((8, 12), (13, 17))), - mask=models._MaskBase.types['No mask'], + mask=models.Mask.types['No mask'], activity=models.Activity.types['Seated'], expiration=models.Expiration.types['Talking'], ), @@ -505,7 +505,7 @@ baseline_model = models.ExposureModel( number=10, presence=models.SpecificInterval(((8, 12), (13, 17))), activity=models.Activity.types['Seated'], - mask=models._MaskBase.types['No mask'], + mask=models.Mask.types['No mask'], ), ) @@ -515,10 +515,10 @@ class CARAStateBuilder(state.StateBuilder): # For example, build_type__VentilationBase is called when dealing with ConcentrationModel # types as it has a ventilation: _VentilationBase field. - def build_type__MaskBase(self, _: dataclasses.Field): + def build_type_Mask(self, _: dataclasses.Field): return state.DataclassStatePredefined( - models._MaskBase, - choices=models._MaskBase.types, + models.Mask, + choices=models.Mask.types, ) def build_type_Virus(self, _: dataclasses.Field): From cb1b31c820d2ba1a676394fd6f01f8d6766ab07b Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sun, 30 May 2021 19:31:56 +0200 Subject: [PATCH 34/83] Starting to modify exposure tests according to new Mask class --- cara/tests/models/test_exposure_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 464f25a9..c1e9090d 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -43,7 +43,7 @@ populations = [ ), # A population with some array component for η_inhale. models.Population( - 10, halftime, models.Mask(0.95, 0.15, np.array([0.3, 0.35])), + 10, halftime, models.Mask(np.array([0.3, 0.35])), models.Activity.types['Standing'], ), # A population with some array component for inhalation_rate. From 421495241e47e704614a0180f4f97ec61c0a69bd Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 06:32:17 +0200 Subject: [PATCH 35/83] Fixing the remaining tests, with the new masks --- cara/tests/models/test_exposure_model.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index c1e9090d..009cf092 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -60,10 +60,10 @@ populations = [ np.array([14.4, 14.4]), np.array([99.6803184113, 99.5181053773])], [populations[2], KnownConcentrations(lambda t: 1.2), - np.array([14.4, 14.4]), np.array([99.4146994564, 99.6803184113])], + np.array([14.4, 14.4]), np.array([97.4574432074, 98.3493482895])], [populations[0], KnownConcentrations(lambda t: np.array([1.2, 2.4])), - np.array([14.4, 28.8]), np.array([99.6803184113, 99.9989780368])], + np.array([14.4, 28.8]), np.array([98.3493482895, 99.9727534893])], [populations[1], KnownConcentrations(lambda t: np.array([1.2, 2.4])), np.array([14.4, 28.8]), np.array([99.6803184113, 99.9976777757])], @@ -123,22 +123,22 @@ def conc_model(): models.InfectedPopulation( number=1, presence=interesting_times, - mask=models.Mask.types['Type I'], + mask=models.Mask.types['No mask'], activity=models.Activity.types['Seated'], virus=models.Virus.types['SARS_CoV_2'], - expiration=models.Expiration.types['Breathing'], + expiration=models.Expiration.types['Superspreading event'], ) ) # expected quanta were computed with a trapezoidal integration, using # a mesh of 10'000 pts per exposed presence interval. @pytest.mark.parametrize("exposed_time_interval, expected_quanta", [ - [(0, 1), 0.0055680845], - [(1, 1.01), 6.4960491e-05], - [(1.01, 1.02), 6.3187723e-05], - [(12, 12.01), 1.9307359e-06], - [(12, 24), 0.079347465], - [(0, 24), 0.086122050], + [(0, 1), 5.4869151], + [(1, 1.01), 0.064013521], + [(1.01, 1.02), 0.062266596], + [(12, 12.01), 0.0019025904], + [(12, 24), 78.190763], + [(0, 24), 84.866592], ] ) def test_exposure_model_integral_accuracy(exposed_time_interval, From be5aea72f0d288d663847e5a48002687fee39ac7 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 06:55:07 +0200 Subject: [PATCH 36/83] Mask exhale_efficiency use now a diameter in microns --- cara/models.py | 14 +++++++------- cara/tests/models/test_mask.py | 12 ++++++------ 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/cara/models.py b/cara/models.py index 736e0992..56f663f4 100644 --- a/cara/models.py +++ b/cara/models.py @@ -489,14 +489,14 @@ class Mask: therein (Asadi 2020). Obtained from measurements of filtration efficiency and of the leakage through the sides. - Diameter is in cm. + Diameter is in microns. """ - if diameter < 0.5e-4: + if diameter < 0.5: eta_out = 0. - elif diameter < 0.94614e-4: - eta_out = 0.5893 * diameter * 1e4 + 0.1546 - elif diameter < 3e-4: - eta_out = 0.0509 * diameter * 1e4 + 0.664 + elif diameter < 0.94614: + eta_out = 0.5893 * diameter + 0.1546 + elif diameter < 3.: + eta_out = 0.0509 * diameter + 0.664 else: eta_out = 0.8167 return eta_out*self.factor_exhale @@ -551,7 +551,7 @@ class Expiration(_ExpirationBase): total = 0 for diameter, factor in zip(self.particle_sizes, self.ejection_factor): contribution = (volume(diameter) * factor * - (1 - mask.exhale_efficiency(diameter))) + (1 - mask.exhale_efficiency(diameter*1e4))) total += contribution return total diff --git a/cara/tests/models/test_mask.py b/cara/tests/models/test_mask.py index 7f7cb242..4a71d13c 100644 --- a/cara/tests/models/test_mask.py +++ b/cara/tests/models/test_mask.py @@ -11,7 +11,7 @@ from cara import models [np.array([0.3, 0.5]), np.array([0.3, 0.5])], ], ) -def test_masks_inhale(η_inhale, expected_inhale_efficiency): +def test_mask_inhale(η_inhale, expected_inhale_efficiency): mask = models.Mask(η_inhale=η_inhale) npt.assert_equal(mask.inhale_efficiency(), expected_inhale_efficiency) @@ -20,11 +20,11 @@ def test_masks_inhale(η_inhale, expected_inhale_efficiency): @pytest.mark.parametrize( "diameter, factor_exhale, expected_exhale_efficiency", [ - [0.3e-4, 1., 0.], - [0.7e-4, 0.3, 0.56711*0.3], - [1.e-4, 1., 0.7149], - [4.e-4, 0.5, 0.8167*0.5], - [5.e-4, 0., 0.], + [0.3, 1., 0.], + [0.7, 0.3, 0.56711*0.3], + [1., 1., 0.7149], + [4., 0.5, 0.8167*0.5], + [5., 0., 0.], ], ) def test_mask_exhale(diameter, factor_exhale, expected_exhale_efficiency): From a61e900178faa46b331490dd1546c3b05ff01e77 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 08:55:55 +0200 Subject: [PATCH 37/83] Fixing docstring --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 56f663f4..88e5493a 100644 --- a/cara/models.py +++ b/cara/models.py @@ -529,7 +529,7 @@ class _ExpirationBase: types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]] def aerosols(self, mask: Mask): - # total volume of aerosols expired per volume of air (mL/cm^3). + # total volume of aerosols expired (cm^3). raise NotImplementedError("Subclass must implement") From 8e68bd7a2317b0ccabcbf91298dd639c7c5deda9 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 11:07:53 +0200 Subject: [PATCH 38/83] Removing old Expiration and ref. to ExpirationBLO in test_expiration --- cara/tests/test_expiration.py | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py index 5cfe94b9..ff08798b 100644 --- a/cara/tests/test_expiration.py +++ b/cara/tests/test_expiration.py @@ -9,7 +9,7 @@ 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)) + e_base = models.Expiration((1., 0., 0.)) with pytest.raises( ValueError, match=re.escape("expirations and weigths should contain the" @@ -19,19 +19,10 @@ def test_multiple_wrong_weight_size(): 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) - mask = models.Mask.types['No mask'] - assert e.aerosols(mask) == e1.aerosols(mask)/3. + 2*e2.aerosols(mask)/3. - - -def test_multiple_BLO(): weights = (1., 1.) - e1 = models.ExpirationBLO((1., 0., 0.)) - e2 = models.ExpirationBLO((4., 5., 5.)) - e_expected = models.ExpirationBLO((2.5, 2.5, 2.5)) + e1 = models.Expiration((1., 0., 0.)) + e2 = models.Expiration((4., 5., 5.)) + e_expected = models.Expiration((2.5, 2.5, 2.5)) e = models.MultipleExpiration([e1, e2], weights) mask = models.Mask.types['Type I'] - assert e_expected.aerosols(mask) == e.aerosols(mask) + npt.assert_almost_equal(e_expected.aerosols(mask), e.aerosols(mask)) From ecf9ce8655a48c6186a8d80e9fc176fd290db756 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 11:08:42 +0200 Subject: [PATCH 39/83] Removing old Expiration class; replaced by previously named ExpirationBLO class --- cara/models.py | 87 +++++++++++++++++--------------------------------- 1 file changed, 30 insertions(+), 57 deletions(-) diff --git a/cara/models.py b/cara/models.py index 284acdd0..f88107c9 100644 --- a/cara/models.py +++ b/cara/models.py @@ -538,53 +538,6 @@ class _ExpirationBase: @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. - - def aerosols(self, mask: Mask): - def volume(diameter): - return (4 * np.pi * (diameter/2)**3) / 3 - total = 0 - for diameter, factor in zip(self.particle_sizes, self.ejection_factor): - contribution = (volume(diameter) * factor * - (1 - mask.exhale_efficiency(diameter*1e4))) - total += contribution - return total - - -@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) - - -@dataclass(frozen=True) -class ExpirationBLO(_ExpirationBase): """ BLO model for the expiration (G. Johnson et al., Modality of human expired aerosol size distributions, Journal of Aerosol Science, @@ -597,7 +550,7 @@ class ExpirationBLO(_ExpirationBase): # speaking, singing, or shouting). BLO_factors: typing.Tuple[float, float, float] - def aerosols(self, mask: _MaskBase): + def aerosols(self, mask: Mask): """ Result is in mL.cm^-3 """ def volume(d): return (np.pi * d**3) / 6. @@ -626,16 +579,36 @@ class ExpirationBLO(_ExpirationBase): return quad(integrand, 0.1, 30)[0] +@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)), - 'Unmodulated Vocalization': Expiration((0.751, 0.139, 0.0139, 0.059)), - 'Superspreading event': Expiration((np.inf, np.inf, np.inf, np.inf)), - 'Breathing BLO': ExpirationBLO((1., 0., 0.)), - 'Talking BLO': ExpirationBLO((1., 1., 1.)), - 'Shouting BLO': ExpirationBLO((5., 5., 5.)), - 'Singing BLO': ExpirationBLO((5., 5., 5.)), + 'Breathing': Expiration((1., 0., 0.)), + 'Talking': Expiration((1., 1., 1.)), + 'Shouting': Expiration((5., 5., 5.)), + 'Singing': Expiration((5., 5., 5.)), + 'Superspreading event': Expiration((np.inf, 0., 0.)), } From d12e809d8011cb4b4759594df78ce2239ac71784 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 11:29:49 +0200 Subject: [PATCH 40/83] Better import of scipy.integrate --- cara/models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cara/models.py b/cara/models.py index f88107c9..24e99e4e 100644 --- a/cara/models.py +++ b/cara/models.py @@ -37,7 +37,7 @@ import typing import numpy as np from scipy.interpolate import interp1d -from scipy.integrate import quad +import scipy.integrate if not typing.TYPE_CHECKING: from memoization import cached @@ -576,7 +576,7 @@ class Expiration(_ExpirationBase): self.BLO_factors[2] * _Omode(d) ) * volume(d*1e-4) * (1 - mask.exhale_efficiency(d*1e-4)) - return quad(integrand, 0.1, 30)[0] + return scipy.integrate.quad(integrand, 0.1, 30)[0] @dataclass(frozen=True) From 612a4c8e08102546bcbe7fd60f22a22c256f67ea Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 16:38:00 +0200 Subject: [PATCH 41/83] New expiration test to make sure expiration rates are correct --- cara/tests/test_expiration.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py index ff08798b..1a6fdbb7 100644 --- a/cara/tests/test_expiration.py +++ b/cara/tests/test_expiration.py @@ -26,3 +26,18 @@ def test_multiple(): e = models.MultipleExpiration([e1, e2], weights) mask = models.Mask.types['Type I'] npt.assert_almost_equal(e_expected.aerosols(mask), e.aerosols(mask)) + + +# expected values obtained from separate model +@pytest.mark.parametrize( + "BLO_weights, expected_aerosols", + [ + [(1.,0.,0.), 1.38924e-12], + [(1.,1.,1.), 1.07129e-10], + [(1.,5.,5.), 5.30088e-10], + ], +) +def test_expiration_aerosols(BLO_weights, expected_aerosols): + mask = models.Mask.types['No mask'] + e = models.Expiration(BLO_weights) + npt.assert_allclose(e.aerosols(mask), expected_aerosols, rtol=1e-4) From 3a390435229756cd81bf84ee139d5fbddcb79974 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 16:41:41 +0200 Subject: [PATCH 42/83] More accurate B, L and O modes; correcting units --- cara/models.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/cara/models.py b/cara/models.py index 24e99e4e..a5c63b2e 100644 --- a/cara/models.py +++ b/cara/models.py @@ -557,26 +557,27 @@ class Expiration(_ExpirationBase): def _Bmode(d: float) -> float: # B-mode (see ref. above). - return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.26)) * - np.exp(-1 * (np.log(d) - 1.0) ** 2 / (2 * 0.26 ** 2))) + return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.262364)) * + np.exp(-1 * (np.log(d) - 0.989541) ** 2 / (2 * 0.262364 ** 2))) def _Lmode(d: float) -> float: # L-mode (see ref. above). - return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.5)) * - np.exp(-1 * (np.log(d) - 1.4) ** 2 / (2 * 0.5 ** 2))) + return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.506818)) * + np.exp(-1 * (np.log(d) - 1.38629) ** 2 / (2 * 0.506818 ** 2))) def _Omode(d: float) -> float: # O-mode (see ref. above). - return ( (1 / d) * (0.001 / (np.sqrt(2 * np.pi) * 0.56)) * - np.exp(-1 * (np.log(d) - 4.98) ** 2 / (2 * 0.56 ** 2))) + return ( (1 / d) * (0.0010008 / (np.sqrt(2 * np.pi) * 0.585005)) * + np.exp(-1 * (np.log(d) - 4.97673) ** 2 / (2 * 0.585005 ** 2))) def integrand(d: float) -> float: return (self.BLO_factors[0] * _Bmode(d) + self.BLO_factors[1] * _Lmode(d) + self.BLO_factors[2] * _Omode(d) - ) * volume(d*1e-4) * (1 - mask.exhale_efficiency(d*1e-4)) + ) * volume(d) * (1 - mask.exhale_efficiency(d)) - return scipy.integrate.quad(integrand, 0.1, 30)[0] + # final result converted from microns^3/cm3 to mL/cm^3 + return scipy.integrate.quad(integrand, 0.1, 30.)[0]*1e-12 @dataclass(frozen=True) @@ -606,8 +607,8 @@ class MultipleExpiration(_ExpirationBase): _ExpirationBase.types = { 'Breathing': Expiration((1., 0., 0.)), 'Talking': Expiration((1., 1., 1.)), - 'Shouting': Expiration((5., 5., 5.)), - 'Singing': Expiration((5., 5., 5.)), + 'Shouting': Expiration((1., 5., 5.)), + 'Singing': Expiration((1., 5., 5.)), 'Superspreading event': Expiration((np.inf, 0., 0.)), } From 62a4880a5b92227b1b4887337f0c63e523d7ca4a Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 17:12:18 +0200 Subject: [PATCH 43/83] Avoiding vectorisation of factor_exhale, and changing test_infected_population accordingly (incompatibility with scipy quad integration) --- cara/models.py | 2 +- cara/tests/test_infected_population.py | 8 ++------ 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/cara/models.py b/cara/models.py index a5c63b2e..481074cd 100644 --- a/cara/models.py +++ b/cara/models.py @@ -478,7 +478,7 @@ class Mask: η_inhale: _VectorisedFloat #: Global factor applied to filtration efficiency of masks when exhaling. - factor_exhale: _VectorisedFloat = 1. + factor_exhale: float = 1. #: Pre-populated examples of Masks. types: typing.ClassVar[typing.Dict[str, "Mask"]] diff --git a/cara/tests/test_infected_population.py b/cara/tests/test_infected_population.py index 8ef94729..7acef7f7 100644 --- a/cara/tests/test_infected_population.py +++ b/cara/tests/test_infected_population.py @@ -8,7 +8,6 @@ import cara.models "override_params", [ {'viral_load_in_sputum': np.array([5e8, 1e9])}, {'quantum_infectious_dose': np.array([50, 20])}, - {'factor_exhale': np.array([0.92, 0.95])}, {'exhalation_rate': np.array([0.75, 0.81])}, ] ) @@ -16,7 +15,6 @@ def test_infected_population_vectorisation(override_params): defaults = { 'viral_load_in_sputum': 1e9, 'quantum_infectious_dose': 50, - 'factor_exhale': 0.95, 'exhalation_rate': 0.75, } defaults.update(override_params) @@ -26,7 +24,7 @@ def test_infected_population_vectorisation(override_params): number=1, presence=office_hours, mask=cara.models.Mask( - factor_exhale=defaults['factor_exhale'], + factor_exhale=0.95, η_inhale=0.3, ), activity=cara.models.Activity( @@ -37,9 +35,7 @@ def test_infected_population_vectorisation(override_params): viral_load_in_sputum=defaults['viral_load_in_sputum'], quantum_infectious_dose=defaults['quantum_infectious_dose'], ), - expiration=cara.models.Expiration( - ejection_factor=(0.084, 0.009, 0.003, 0.002), - ), + expiration=cara.models.Expiration((1., 0., 0.)), ) emission_rate = infected.emission_rate(10) assert isinstance(emission_rate, np.ndarray) From 4a211aa7fc453aa6daa8a03c6b8cb106e9bb1ee3 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 23:09:22 +0200 Subject: [PATCH 44/83] Fixing tests with the change of Expiration, in particular: removing useless tests on aerosols in test_known_quantities; removing reference to Unmodulated vocalization expiration type; avoiding the use of factor_exhale as array; changing some integration test values on r0; replacing some integration tests on r0 by integration tests on quanta_exposure with a simpler model; changing values in some concentration tests --- cara/tests/conftest.py | 12 ++- cara/tests/models/test_concentration_model.py | 8 +- cara/tests/test_expiration.py | 2 +- cara/tests/test_known_quantities.py | 77 ++++++++----------- cara/tests/test_monte_carlo.py | 2 +- 5 files changed, 39 insertions(+), 62 deletions(-) diff --git a/cara/tests/conftest.py b/cara/tests/conftest.py index 24983cf1..44b25423 100644 --- a/cara/tests/conftest.py +++ b/cara/tests/conftest.py @@ -7,11 +7,9 @@ import pytest def baseline_model(): model = models.ConcentrationModel( room=models.Room(volume=75), - ventilation=models.SlidingWindow( - active=models.PeriodicInterval(period=120, duration=120), - inside_temp=models.PiecewiseConstant((0,24),(293,)), - outside_temp=models.PiecewiseConstant((0,24),(283,)), - window_height=1.6, opening_length=0.6, + ventilation=models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=30., ), infected=models.InfectedPopulation( number=1, @@ -19,7 +17,7 @@ def baseline_model(): presence=models.SpecificInterval(((0, 4), (5, 8))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light activity'], - expiration=models.Expiration.types['Unmodulated Vocalization'], + expiration=models.Expiration.types['Superspreading event'], ), ) return model @@ -30,7 +28,7 @@ def baseline_exposure_model(baseline_model): return models.ExposureModel( baseline_model, exposed=models.Population( - number=10, + number=1000, presence=baseline_model.infected.presence, activity=baseline_model.infected.activity, mask=baseline_model.infected.mask, diff --git a/cara/tests/models/test_concentration_model.py b/cara/tests/models/test_concentration_model.py index f7230430..26707236 100644 --- a/cara/tests/models/test_concentration_model.py +++ b/cara/tests/models/test_concentration_model.py @@ -14,7 +14,6 @@ from cara import models {'air_change': np.array([100, 120])}, {'viral_load_in_sputum': np.array([5e8, 1e9])}, {'quantum_infectious_dose': np.array([50, 20])}, - {'factor_exhale': np.array([0.92, 0.95])}, ] ) def test_concentration_model_vectorisation(override_params): @@ -24,7 +23,6 @@ def test_concentration_model_vectorisation(override_params): 'air_change': 100, 'viral_load_in_sputum': 1e9, 'quantum_infectious_dose': 50, - 'factor_exhale': 0.95, } defaults.update(override_params) @@ -36,7 +34,7 @@ def test_concentration_model_vectorisation(override_params): number=1, presence=always, mask=models.Mask( - factor_exhale=defaults['factor_exhale'], + factor_exhale=0.95, η_inhale=0.3, ), activity=models.Activity( @@ -47,9 +45,7 @@ def test_concentration_model_vectorisation(override_params): viral_load_in_sputum=defaults['viral_load_in_sputum'], quantum_infectious_dose=defaults['quantum_infectious_dose'], ), - expiration=models.Expiration( - ejection_factor=(0.084, 0.009, 0.003, 0.002), - ), + expiration=models.Expiration((1., 0., 0.)), ) ) concentrations = c_model.concentration(10) diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py index 1a6fdbb7..434caf8c 100644 --- a/cara/tests/test_expiration.py +++ b/cara/tests/test_expiration.py @@ -28,7 +28,7 @@ def test_multiple(): npt.assert_almost_equal(e_expected.aerosols(mask), e.aerosols(mask)) -# expected values obtained from separate model +# expected values obtained from another code @pytest.mark.parametrize( "BLO_weights, expected_aerosols", [ diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index ea8a1ad0..a3b1e5db 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -6,21 +6,12 @@ import cara.models as models import cara.data as data -def test_no_mask_aerosols(baseline_model): - exp = models.Expiration.types['Unmodulated Vocalization'] - npt.assert_allclose( - exp.aerosols(models.Mask.types['No mask']), - 6.077541e-12, - rtol=1e-5, - ) - - -def test_no_mask_emission_rate(baseline_model): - rate = 151.938514 +def test_no_mask_superspeading_emission_rate(baseline_model): + expected_rate = 970. npt.assert_allclose( [baseline_model.infected.emission_rate(t) for t in [0, 1, 4, 4.5, 5, 8, 9]], - [0, rate, rate, 0, 0, rate, 0], - rtol=1e-5 + [0, expected_rate, expected_rate, 0, 0, expected_rate, 0], + rtol=1e-12 ) @@ -48,23 +39,24 @@ def baseline_periodic_hepa(): def test_concentrations(baseline_model): + # expected concentrations were computed analytically ts = [0, 4, 5, 7, 10] concentrations = [baseline_model.concentration(t) for t in ts] npt.assert_allclose( concentrations, - [0.000000e+00, 2.619538e-01, 1.146999e-04, 2.619537e-01, 5.022289e-08], - rtol=1e-5 + [0.000000e+00, 0.4189594, 1.6422648e-14, 0.4189594, 6.4374587e-28], + rtol=1e-6 ) def test_smooth_concentrations(baseline_model): # We don't care about the actual concentrations in this test, but rather # that the curve itself is smooth. - dx = 0.1 - dy_limit = dx * 2 # Anything more than this is a bit steep. + dx = 0.002 + dy_limit = 0.2 # Anything more than this (in relative) is a bit steep. ts = np.arange(0, 10, dx) concentrations = [baseline_model.concentration(t) for t in ts] - assert np.abs(np.diff(concentrations)).max() < dy_limit + assert np.abs(np.diff(concentrations)).max()/np.mean(concentrations) < dy_limit def build_model(interval_duration): @@ -80,7 +72,7 @@ def build_model(interval_duration): presence=models.SpecificInterval(((0, 4), (5, 8))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light activity'], - expiration=models.Expiration.types['Unmodulated Vocalization'], + expiration=models.Expiration.types['Superspreading event'], ), ) return model @@ -98,8 +90,8 @@ def test_concentrations_startup(baseline_model): def test_r0(baseline_exposure_model): # expected r0 was computed with a trapezoidal integration, using # a mesh of 100'000 pts per exposed presence interval. - p = baseline_exposure_model.infection_probability() - npt.assert_allclose(p, 89.001748) + r0 = baseline_exposure_model.reproduction_number() + npt.assert_allclose(r0, 973.535888) def test_periodic_window(baseline_periodic_window, baseline_room): @@ -183,15 +175,6 @@ def test_multiple_ventilation_HEPA_HVAC_AirChange(volume, expected_value): expected_value,rtol=1e-5) -def test_expiration_aerosols(): - mask = models.Mask.types['Type I'] - exp1 = models.Expiration((0.751, 0.139, 0.0139, 0.059), - particle_sizes = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4)) - exp2 = models.Expiration((0.059, 0.0139, 0.751, 0.139), - particle_sizes = (5.5e-4, 3.5e-4, 0.8e-4, 1.8e-4)) - npt.assert_allclose(exp1.aerosols(mask), exp2.aerosols(mask), rtol=1e-5) - - @pytest.mark.parametrize( "time, expected_value", [ @@ -239,7 +222,7 @@ def build_hourly_dependent_model(month, intervals_open=((7.5, 8.5),), presence=models.SpecificInterval(intervals_presence_infected), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light activity'], - expiration=models.Expiration.types['Unmodulated Vocalization'], + expiration=models.Expiration.types['Superspreading event'], ), ) return model @@ -260,7 +243,7 @@ def build_constant_temp_model(outside_temp, intervals_open=((7.5, 8.5),)): presence=models.SpecificInterval(((0, 4), (5, 7.5))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light activity'], - expiration=models.Expiration.types['Unmodulated Vocalization'], + expiration=models.Expiration.types['Superspreading event'], ), ) return model @@ -287,7 +270,7 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5 presence=models.SpecificInterval(((0, 4), (5, 7.5))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light activity'], - expiration=models.Expiration.types['Unmodulated Vocalization'], + expiration=models.Expiration.types['Superspreading event'], ), ) return model @@ -370,16 +353,16 @@ def build_exposure_model(concentration_model): ) -# expected r0 were computed with a trapezoidal integration, using +# expected quanta were computed with a trapezoidal integration, using # a mesh of 100'000 pts per exposed presence interval. @pytest.mark.parametrize( - "month, expected_r0", + "month, expected_quanta", [ - ['Jan', 86.258533], - ['Jun', 99.972103], + ['Jan', 10.136783], + ['Jun', 41.800377], ], ) -def test_r0_hourly_dep(month,expected_r0): +def test_quanta_hourly_dep(month,expected_quanta): m = build_exposure_model( build_hourly_dependent_model( month, @@ -387,19 +370,19 @@ def test_r0_hourly_dep(month,expected_r0): intervals_presence_infected=((8, 12), (13, 17)) ) ) - p = m.infection_probability() - npt.assert_allclose(p, expected_r0) + quanta = m.quanta_exposure() + npt.assert_allclose(quanta, expected_quanta) -# expected r0 were computed with a trapezoidal integration, using +# expected quanta were computed with a trapezoidal integration, using # a mesh of 100'000 pts per exposed presence interval. @pytest.mark.parametrize( - "month, expected_r0", + "month, expected_quanta", [ - ['Jan', 86.422736], - ['Jun', 99.982323], + ['Jan', 10.19818], + ['Jun', 44.130683], ], ) -def test_r0_hourly_dep_refined(month,expected_r0): +def test_quanta_hourly_dep_refined(month,expected_quanta): m = build_exposure_model( build_hourly_dependent_model( month, @@ -408,5 +391,5 @@ def test_r0_hourly_dep_refined(month,expected_r0): temperatures=data.GenevaTemperatures, ) ) - p = m.infection_probability() - npt.assert_allclose(p, expected_r0) + quanta = m.quanta_exposure() + npt.assert_allclose(quanta, expected_quanta) diff --git a/cara/tests/test_monte_carlo.py b/cara/tests/test_monte_carlo.py index f3d9fb08..c3bfc2fe 100644 --- a/cara/tests/test_monte_carlo.py +++ b/cara/tests/test_monte_carlo.py @@ -53,7 +53,7 @@ def baseline_mc_model() -> cara.monte_carlo.ConcentrationModel: presence=cara.models.SpecificInterval(((0, 4), (5, 8))), mask=cara.models.Mask.types['No mask'], activity=cara.models.Activity.types['Light activity'], - expiration=cara.models.Expiration.types['Unmodulated Vocalization'], + expiration=cara.models.Expiration.types['Breathing'], ), ) return mc_model From 390d2299445296465075b7981c23dd34ed5955f4 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 31 May 2021 23:23:42 +0200 Subject: [PATCH 45/83] Caching aerosols computation --- cara/models.py | 1 + 1 file changed, 1 insertion(+) diff --git a/cara/models.py b/cara/models.py index 481074cd..cf2ffb8f 100644 --- a/cara/models.py +++ b/cara/models.py @@ -550,6 +550,7 @@ class Expiration(_ExpirationBase): # speaking, singing, or shouting). BLO_factors: typing.Tuple[float, float, float] + @cached() def aerosols(self, mask: Mask): """ Result is in mL.cm^-3 """ def volume(d): From 8193a51f99d6f265704720b9673002bb791bc8d7 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Tue, 1 Jun 2021 09:21:09 +0200 Subject: [PATCH 46/83] Adding tests on fraction_deposited (ExposureModel) - including vectorisation; modifying existing exposure model tests --- cara/tests/conftest.py | 3 ++- cara/tests/models/test_exposure_model.py | 19 +++++++++++-------- cara/tests/test_known_quantities.py | 3 ++- 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/cara/tests/conftest.py b/cara/tests/conftest.py index 44b25423..023e4e5a 100644 --- a/cara/tests/conftest.py +++ b/cara/tests/conftest.py @@ -32,5 +32,6 @@ def baseline_exposure_model(baseline_model): presence=baseline_model.infected.presence, activity=baseline_model.infected.activity, mask=baseline_model.infected.mask, - ) + ), + fraction_deposited = 1., ) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 009cf092..35c5ac0c 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -55,22 +55,25 @@ populations = [ @pytest.mark.parametrize( - "population, cm, expected_exposure, expected_probability",[ - [populations[1], KnownConcentrations(lambda t: 1.2), + "population, cm, f_dep, expected_exposure, expected_probability",[ + [populations[1], KnownConcentrations(lambda t: 1.2), 1., np.array([14.4, 14.4]), np.array([99.6803184113, 99.5181053773])], - [populations[2], KnownConcentrations(lambda t: 1.2), + [populations[2], KnownConcentrations(lambda t: 1.2), 1., np.array([14.4, 14.4]), np.array([97.4574432074, 98.3493482895])], - [populations[0], KnownConcentrations(lambda t: np.array([1.2, 2.4])), + [populations[0], KnownConcentrations(lambda t: np.array([1.2, 2.4])), 1., np.array([14.4, 28.8]), np.array([98.3493482895, 99.9727534893])], - [populations[1], KnownConcentrations(lambda t: np.array([1.2, 2.4])), + [populations[1], KnownConcentrations(lambda t: np.array([1.2, 2.4])), 1., np.array([14.4, 28.8]), np.array([99.6803184113, 99.9976777757])], + + [populations[0], KnownConcentrations(lambda t: 2.4), np.array([0.5, 1.]), + 28.8, np.array([98.3493482895, 99.9727534893])], ]) -def test_exposure_model_ndarray(population, cm, +def test_exposure_model_ndarray(population, cm, f_dep, expected_exposure, expected_probability): - model = ExposureModel(cm, population) + model = ExposureModel(cm, population, fraction_deposited = f_dep) np.testing.assert_almost_equal( model.quanta_exposure(), expected_exposure ) @@ -148,5 +151,5 @@ def test_exposure_model_integral_accuracy(exposed_time_interval, 10, presence_interval, models.Mask.types['Type I'], models.Activity.types['Standing'], ) - model = ExposureModel(conc_model, population) + model = ExposureModel(conc_model, population, fraction_deposited=1.) np.testing.assert_allclose(model.quanta_exposure(), expected_quanta) diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index a3b1e5db..4af1c06f 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -349,7 +349,8 @@ def build_exposure_model(concentration_model): presence=infected.presence, activity=infected.activity, mask=infected.mask, - ) + ), + fraction_deposited = 1., ) From 02cbc2be303b24001c527b84c5c687483e985848 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Tue, 1 Jun 2021 09:22:15 +0200 Subject: [PATCH 47/83] Adding fraction_deposited to ExposureModel (vectorised); improving docstrings --- cara/models.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index cf2ffb8f..794ed21b 100644 --- a/cara/models.py +++ b/cara/models.py @@ -616,7 +616,10 @@ _ExpirationBase.types = { @dataclass(frozen=True) class Activity: + #: Inhalation rate in m^3/h inhalation_rate: _VectorisedFloat + + #: Exhalation rate in m^3/h exhalation_rate: _VectorisedFloat #: Pre-populated examples of activities. @@ -671,6 +674,8 @@ class InfectedPopulation(Population): """ # Emission Rate (infectious quantum / h) + # Note on units: exhalation rate is in m^3/h, aerosols in mL/cm^3 + # and viral load in virus/mL -> 1e6 conversion factor aerosols = self.expiration.aerosols(self.mask) ER = (self.virus.viral_load_in_sputum * @@ -855,6 +860,9 @@ class ExposureModel: #: The number of times the exposure event is repeated (default 1). repeats: int = 1 + #: The fraction of viruses actually deposited in the respiratory tract + fraction_deposited: _VectorisedFloat = 0.6 + def quanta_exposure(self) -> _VectorisedFloat: """The number of virus quanta per meter^3.""" exposure = 0.0 @@ -870,7 +878,7 @@ class ExposureModel: inf_aero = ( self.exposed.activity.inhalation_rate * (1 - self.exposed.mask.inhale_efficiency()) * - exposure + exposure * self.fraction_deposited ) # Probability of infection. From bef2d4949ebc3a097bceccc503bc6cf711555623 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Tue, 1 Jun 2021 21:39:01 +0200 Subject: [PATCH 48/83] Adding tests on different kind of sampleable distributions --- cara/tests/test_sampleable_distribution.py | 81 ++++++++++++++++++++++ 1 file changed, 81 insertions(+) create mode 100644 cara/tests/test_sampleable_distribution.py diff --git a/cara/tests/test_sampleable_distribution.py b/cara/tests/test_sampleable_distribution.py new file mode 100644 index 00000000..b53f9ad6 --- /dev/null +++ b/cara/tests/test_sampleable_distribution.py @@ -0,0 +1,81 @@ +import numpy as np +import numpy.testing as npt +import pytest + +from cara.monte_carlo import sampleable + +@pytest.mark.parametrize( + "mean, std",[ + [1., 0.5], + ] +) +def test_normal(mean, std): + # test that the sample has approximately the right mean, + # std deviation and distribution function. + sample_size = 2000000 + samples = sampleable.Normal(mean, std).generate_samples(sample_size) + histogram, bins = np.histogram(samples,bins=100, density=True) + x = (bins[1:]+bins[:-1])/2 + exact_dist = 1/(np.sqrt(2*np.pi)*std) * np.exp(-((x-mean)/std)**2/2) + + assert len(samples) == sample_size + npt.assert_allclose([samples.mean(), samples.std()], [mean, std], atol=mean/100) + npt.assert_allclose(histogram, exact_dist, atol=exact_dist.mean()/20) + + +@pytest.mark.parametrize( + "mean_gaussian, std_gaussian",[ + [-0.6872121723362303, 0.10498338229297108], + ] +) +def test_lognormal(mean_gaussian, std_gaussian): + # test that the sample has approximately the right mean, + # std deviation and distribution function. + sample_size = 2000000 + samples = sampleable.LogNormal(mean_gaussian, std_gaussian + ).generate_samples(sample_size) + histogram, bins = np.histogram(samples,bins=50, density=True) + x = (bins[1:]+bins[:-1])/2 + exact_dist = ( 1/(x*np.sqrt(2*np.pi)*std_gaussian) * + np.exp(-((np.log(x)-mean_gaussian)/std_gaussian)**2/2) ) + exact_mean = np.exp(mean_gaussian + std_gaussian**2/2) + exact_std = np.sqrt( (np.exp(std_gaussian**2)-1) * + np.exp(2*mean_gaussian + std_gaussian**2) ) + + assert len(samples) == sample_size + npt.assert_allclose([samples.mean(), samples.std()], + [exact_mean, exact_std], atol=exact_mean/100) + npt.assert_allclose(histogram, exact_dist, atol=exact_dist.mean()/20) + + +@pytest.mark.parametrize( + "use_kernel", + [False, True], +) +def test_custom(use_kernel): + # test that the sample has approximately the right distribution + # function, with both Custom and CustomKernel method. The latter + # is less accurate for smooth functions. + # the distribution function is an inverted parabola, with maximum 0.15, + # which is 0 at the bounds (0,10) (normalized) + norm = 500/3. + function = lambda x: (-(5 - x)**2 + 25)/norm + max_function = 0.15 + sample_size = 2000000 + + if use_kernel: + variable = np.linspace(0.1,9.9,100) + frequencies = function(variable) + samples = sampleable.CustomKernel(variable, frequencies, + kernel_bandwidth=0.1 + ).generate_samples(sample_size) + tolerance = max_function/10 + else: + samples = sampleable.Custom((0, 10), function, max_function + ).generate_samples(sample_size) + tolerance = max_function/50 + + histogram, bins = np.histogram(samples, bins=100, density=True) + correct_dist = function((bins[1:]+bins[:-1])/2) + assert len(samples) == sample_size + npt.assert_allclose(histogram, correct_dist, atol=tolerance) From 0cc1e85ececa08ce1a69c2ac28a7f63a81a2381b Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Tue, 1 Jun 2021 21:40:07 +0200 Subject: [PATCH 49/83] Adding new sampleable distributions: lognormal and custom (with two different methods) - thanks to Markus Rognlien --- cara/monte_carlo/sampleable.py | 70 ++++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 4ed49d82..df7a7646 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -1,6 +1,7 @@ import typing import numpy as np +from sklearn.neighbors import KernelDensity # type: ignore import cara.models @@ -16,6 +17,9 @@ class SampleableDistribution: class Normal(SampleableDistribution): + """ + Defines a normal (i.e. Gaussian) distribution + """ def __init__(self, mean: float, standard_deviation: float): self.mean = mean self.standard_deviation = standard_deviation @@ -24,6 +28,72 @@ class Normal(SampleableDistribution): return np.random.normal(self.mean, self.standard_deviation, size=size) +class LogNormal(SampleableDistribution): + """ + Defines a lognormal distribution (i.e. Gaussian distribution vs. the + natural logarithm of the random variable) + """ + + def __init__(self, mean_gaussian: float, standard_deviation_gaussian: float): + # these are resp. the mean and std. deviation of the underlying + # Gaussian distribution + self.mean_gaussian = mean_gaussian + self.standard_deviation_gaussian = standard_deviation_gaussian + + def generate_samples(self, size: int) -> float_array_size_n: + return np.random.lognormal(self.mean_gaussian, + self.standard_deviation_gaussian, + size=size) + + +class Custom(SampleableDistribution): + """ + Defines a distribution which follows a custom curve vs. the random + variable. Uses a simple algorithm. This is appropriate for a smooth + distribution function (one should know its maximum). + """ + def __init__(self, bounds: typing.Tuple[float, float], + function: typing.Callable, max_function: float): + self.bounds = bounds + self.function = function + self.max_function = max_function + + def generate_samples(self, size: int) -> float_array_size_n: + fvalue = np.random.uniform(0,self.max_function,size) + x = np.random.uniform(*self.bounds,size) + invalid = np.where(fvalue>self.function(x))[0] + while len(invalid)>0: + fvalue[invalid] = np.random.uniform(0,self.max_function,len(invalid)) + x[invalid] = np.random.uniform(*self.bounds,len(invalid)) + invalid = np.where(fvalue>self.function(x))[0] + + return x + + +class CustomKernel(SampleableDistribution): + """ + Defines a distribution which follows a custom curve vs. the + random variable. Uses a Gaussian kernel density fit. This is more + appropriate for a noisy distribution function. + """ + def __init__(self, variable: float_array_size_n, + frequencies: float_array_size_n, + kernel_bandwidth: float): + # these are resp. the random variable, the distribution + # frequencies at these values, and the bandwidth of the Gaussian + # kernel + self.variable = variable + self.frequencies = frequencies + self.kernel_bandwidth = kernel_bandwidth + + def generate_samples(self, size: int) -> float_array_size_n: + kde_model = KernelDensity(kernel='gaussian', + bandwidth=self.kernel_bandwidth) + kde_model.fit(self.variable.reshape(-1, 1), + sample_weight=self.frequencies) + return kde_model.sample(n_samples=size)[:, 0] + + _VectorisedFloatOrSampleable = typing.Union[ SampleableDistribution, cara.models._VectorisedFloat, ] From 867536295019f0cf4ba9aa9dee40b2f549966faa Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Tue, 1 Jun 2021 21:48:04 +0200 Subject: [PATCH 50/83] Minor fix in docstring --- cara/monte_carlo/sampleable.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index df7a7646..95374370 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -79,7 +79,7 @@ class CustomKernel(SampleableDistribution): def __init__(self, variable: float_array_size_n, frequencies: float_array_size_n, kernel_bandwidth: float): - # these are resp. the random variable, the distribution + # these are resp. the random variable, the distribution # frequencies at these values, and the bandwidth of the Gaussian # kernel self.variable = variable From 185d57639744e8f4bcb68635b8f8b6d2a370568e Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 2 Jun 2021 07:16:00 +0200 Subject: [PATCH 51/83] Adding sklearn to the dependencies --- setup.py | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.py b/setup.py index e99d3185..3d5fc17b 100644 --- a/setup.py +++ b/setup.py @@ -29,6 +29,7 @@ REQUIREMENTS: dict = { 'numpy', 'qrcode[pil]', 'scipy', + 'sklearn', 'tornado', 'voila >=0.2.4', ], From b7f4ad8c3e1af27ce7d30d2b55869e5e96301ddb Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 2 Jun 2021 09:06:38 +0200 Subject: [PATCH 52/83] Adding a TODO in sampleable tests --- cara/tests/test_sampleable_distribution.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cara/tests/test_sampleable_distribution.py b/cara/tests/test_sampleable_distribution.py index b53f9ad6..b3887763 100644 --- a/cara/tests/test_sampleable_distribution.py +++ b/cara/tests/test_sampleable_distribution.py @@ -3,6 +3,8 @@ import numpy.testing as npt import pytest from cara.monte_carlo import sampleable +# TODO: seed deterministically the random number generators, here (to +# avoid random issues with tests) @pytest.mark.parametrize( "mean, std",[ From 61b06fc1f0bd80f96472fba85a31b34d20f4d055 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 2 Jun 2021 12:43:10 +0200 Subject: [PATCH 53/83] Adding LogCustomKernel sampleable distribution (curve is vs. a log variable) --- cara/monte_carlo/sampleable.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 95374370..4333c93f 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -94,6 +94,30 @@ class CustomKernel(SampleableDistribution): return kde_model.sample(n_samples=size)[:, 0] +class LogCustomKernel(SampleableDistribution): + """ + Defines a distribution which follows a custom curve vs. the log + (in base 10) of the random variable. Uses a Gaussian kernel density + fit. This is more appropriate for a noisy distribution function. + """ + def __init__(self, log_variable: float_array_size_n, + frequencies: float_array_size_n, + kernel_bandwidth: float): + # these are resp. the log of the random variable, the distribution + # frequencies at these values, and the bandwidth of the Gaussian + # kernel + self.log_variable = log_variable + self.frequencies = frequencies + self.kernel_bandwidth = kernel_bandwidth + + def generate_samples(self, size: int) -> float_array_size_n: + kde_model = KernelDensity(kernel='gaussian', + bandwidth=self.kernel_bandwidth) + kde_model.fit(self.log_variable.reshape(-1, 1), + sample_weight=self.frequencies) + return 10 ** kde_model.sample(n_samples=size)[:, 0] + + _VectorisedFloatOrSampleable = typing.Union[ SampleableDistribution, cara.models._VectorisedFloat, ] From 387fda67d844ce875f529d326b8665b562975629 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 2 Jun 2021 12:45:49 +0200 Subject: [PATCH 54/83] Improving tests on sampleable distributions (more robust; using seed); adding test for LogCustomKernel --- cara/tests/test_sampleable_distribution.py | 57 ++++++++++++++++------ 1 file changed, 42 insertions(+), 15 deletions(-) diff --git a/cara/tests/test_sampleable_distribution.py b/cara/tests/test_sampleable_distribution.py index b3887763..e3801c7e 100644 --- a/cara/tests/test_sampleable_distribution.py +++ b/cara/tests/test_sampleable_distribution.py @@ -3,8 +3,8 @@ import numpy.testing as npt import pytest from cara.monte_carlo import sampleable -# TODO: seed deterministically the random number generators, here (to -# avoid random issues with tests) +# TODO: seed better the random number generators +np.random.seed(2000) @pytest.mark.parametrize( "mean, std",[ @@ -17,12 +17,14 @@ def test_normal(mean, std): sample_size = 2000000 samples = sampleable.Normal(mean, std).generate_samples(sample_size) histogram, bins = np.histogram(samples,bins=100, density=True) - x = (bins[1:]+bins[:-1])/2 - exact_dist = 1/(np.sqrt(2*np.pi)*std) * np.exp(-((x-mean)/std)**2/2) + selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip( + (bins[1:]+bins[:-1])/2,histogram) if b>=0.25 and b<=1.75]) + exact_dist = 1/(np.sqrt(2*np.pi)*std) * np.exp( + -((np.array(selected_bins)-mean)/std)**2/2) assert len(samples) == sample_size - npt.assert_allclose([samples.mean(), samples.std()], [mean, std], atol=mean/100) - npt.assert_allclose(histogram, exact_dist, atol=exact_dist.mean()/20) + npt.assert_allclose([samples.mean(), samples.std()], [mean, std], rtol=0.01) + npt.assert_allclose(selected_histogram, exact_dist, rtol=0.02) @pytest.mark.parametrize( @@ -37,17 +39,19 @@ def test_lognormal(mean_gaussian, std_gaussian): samples = sampleable.LogNormal(mean_gaussian, std_gaussian ).generate_samples(sample_size) histogram, bins = np.histogram(samples,bins=50, density=True) - x = (bins[1:]+bins[:-1])/2 - exact_dist = ( 1/(x*np.sqrt(2*np.pi)*std_gaussian) * - np.exp(-((np.log(x)-mean_gaussian)/std_gaussian)**2/2) ) + selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip( + (bins[1:]+bins[:-1])/2,histogram) if b>=0.4 and b<=0.6]) + selected_bins = np.array(selected_bins) + exact_dist = ( 1/(selected_bins*np.sqrt(2*np.pi)*std_gaussian) * + np.exp(-((np.log(selected_bins)-mean_gaussian)/std_gaussian)**2/2) ) exact_mean = np.exp(mean_gaussian + std_gaussian**2/2) exact_std = np.sqrt( (np.exp(std_gaussian**2)-1) * np.exp(2*mean_gaussian + std_gaussian**2) ) assert len(samples) == sample_size npt.assert_allclose([samples.mean(), samples.std()], - [exact_mean, exact_std], atol=exact_mean/100) - npt.assert_allclose(histogram, exact_dist, atol=exact_dist.mean()/20) + [exact_mean, exact_std], rtol=0.01) + npt.assert_allclose(selected_histogram, exact_dist, rtol=0.02) @pytest.mark.parametrize( @@ -71,13 +75,36 @@ def test_custom(use_kernel): samples = sampleable.CustomKernel(variable, frequencies, kernel_bandwidth=0.1 ).generate_samples(sample_size) - tolerance = max_function/10 else: samples = sampleable.Custom((0, 10), function, max_function ).generate_samples(sample_size) - tolerance = max_function/50 histogram, bins = np.histogram(samples, bins=100, density=True) - correct_dist = function((bins[1:]+bins[:-1])/2) + selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip( + (bins[1:]+bins[:-1])/2,histogram) if b>=1 and b<=9]) + correct_dist = function(np.array(selected_bins)) assert len(samples) == sample_size - npt.assert_allclose(histogram, correct_dist, atol=tolerance) + npt.assert_allclose(selected_histogram, correct_dist, rtol=0.05) + + +def test_logcustomkernel(): + # test that the sample has approximately the right distribution + # function, for the LogCustomKernel. + # the distribution function is an inverted parabola vs. the log of + # the variable (normalized) + norm = 500/3. + function = lambda x: (-(5 - x)**2 + 25)/norm + sample_size = 2000000 + + log_variable = np.linspace(0.1,9.9,100) + frequencies = function(log_variable) + samples = sampleable.LogCustomKernel(log_variable, frequencies, + kernel_bandwidth=0.1 + ).generate_samples(sample_size) + + histogram, bins = np.histogram(np.log10(samples), bins=100, density=True) + selected_bins,selected_histogram = zip(*[(b,h) for b,h in zip( + (bins[1:]+bins[:-1])/2,histogram) if b>=1 and b<=9]) + correct_dist = function(np.array(selected_bins)) + assert len(samples) == sample_size + npt.assert_allclose(selected_histogram, correct_dist, rtol=0.05) From cc628cd1e6f280cab157e1db11229a5aff16a6df Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 3 Jun 2021 10:23:19 +0200 Subject: [PATCH 55/83] Adding tests on predefined distributions --- cara/tests/test_predefined_distributions.py | 47 +++++++++++++++++++++ cara/tests/test_sampleable_distribution.py | 2 + 2 files changed, 49 insertions(+) create mode 100644 cara/tests/test_predefined_distributions.py diff --git a/cara/tests/test_predefined_distributions.py b/cara/tests/test_predefined_distributions.py new file mode 100644 index 00000000..1e066432 --- /dev/null +++ b/cara/tests/test_predefined_distributions.py @@ -0,0 +1,47 @@ +import numpy as np +import numpy.testing as npt +import pytest + +from cara.monte_carlo.data import activity_distributions, virus_distributions + +# TODO: seed better the random number generators +np.random.seed(0) + + +# mean & std deviations from CERN-OPEN-2021-04 (Table 4) +# NOTE: a mistake was corrected for the std deviation of the +# "Moderate exercise" case (0.37 in the report, but should be 0.34) +@pytest.mark.parametrize( + "distribution, mean, std",[ + ['Seated', 0.51, 0.053], + + ['Standing', 0.57, 0.053], + + ['Light activity', 1.24, 0.12], + + ['Moderate activity', 1.77, 0.34], + + ['Heavy exercise', 3.28, 0.72,], + ] +) +def test_activity_distributions(distribution, mean, std): + activity = activity_distributions[distribution].build_model(size=1000000) + npt.assert_allclose(activity.inhalation_rate.mean(), mean, atol=0.01) + npt.assert_allclose(activity.inhalation_rate.std(), std, atol=0.01) + + +# mean & std deviations from CERN-OPEN-2021-04 (Table 4) - with a refined +# precision on the values +@pytest.mark.parametrize( + "distribution, mean, std",[ + ['SARS_CoV_2', 6.59, 1.74], + + ['SARS_CoV_2_B117', 6.59, 1.74], + + ['SARS_CoV_2_P1', 6.59, 1.74], + ] +) +def test_viral_load_logdistribution(distribution, mean, std): + virus = virus_distributions[distribution].build_model(size=1000000) + npt.assert_allclose(np.log10(virus.viral_load_in_sputum).mean(), mean, atol=0.01) + npt.assert_allclose(np.log10(virus.viral_load_in_sputum).std(), std, atol=0.01) diff --git a/cara/tests/test_sampleable_distribution.py b/cara/tests/test_sampleable_distribution.py index e3801c7e..585b4dde 100644 --- a/cara/tests/test_sampleable_distribution.py +++ b/cara/tests/test_sampleable_distribution.py @@ -3,9 +3,11 @@ import numpy.testing as npt import pytest from cara.monte_carlo import sampleable + # TODO: seed better the random number generators np.random.seed(2000) + @pytest.mark.parametrize( "mean, std",[ [1., 0.5], From d246530704ae98959b6e799d195a762e89b39b38 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 3 Jun 2021 10:35:17 +0200 Subject: [PATCH 56/83] Adding data.py with pre-defined distributions for breathing rates and viral loads --- cara/monte_carlo/data.py | 56 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 cara/monte_carlo/data.py diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py new file mode 100644 index 00000000..ad1f4088 --- /dev/null +++ b/cara/monte_carlo/data.py @@ -0,0 +1,56 @@ +import numpy as np + +import cara.monte_carlo as mc +from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel + + +# From CERN-OPEN-2021-04 and refererences therein +activity_distributions = { + 'Seated': mc.Activity(LogNormal(-0.6872121723362303, 0.10498338229297108), + LogNormal(-0.6872121723362303, 0.10498338229297108)), + + 'Standing': mc.Activity(LogNormal(-0.5742377578494785, 0.09373162411398223), + LogNormal(-0.5742377578494785, 0.09373162411398223)), + + 'Light activity': mc.Activity(LogNormal(0.21380242785625422,0.09435378091059601), + LogNormal(0.21380242785625422,0.09435378091059601)), + + 'Moderate activity': mc.Activity(LogNormal(0.551771330362601, 0.1894616357138137), + LogNormal(0.551771330362601, 0.1894616357138137)), + + 'Heavy exercise': mc.Activity(LogNormal(1.1644665696723049, 0.21744554768657565), + LogNormal(1.1644665696723049, 0.21744554768657565)), +} + + +# From CERN-OPEN-2021-04 and refererences therein +log_symptomatic_vl_frequencies = LogCustomKernel( + np.array((2.46032, 2.67431, 2.85434, 3.06155, 3.25856, 3.47256, 3.66957, 3.85979, 4.09927, 4.27081, + 4.47631, 4.66653, 4.87204, 5.10302, 5.27456, 5.46478, 5.6533, 5.88428, 6.07281, 6.30549, + 6.48552, 6.64856, 6.85407, 7.10373, 7.30075, 7.47229, 7.66081, 7.85782, 8.05653, 8.27053, + 8.48453, 8.65607, 8.90573, 9.06878, 9.27429, 9.473, 9.66152, 9.87552)), + np.array((0.001206885, 0.007851618, 0.008078144, 0.01502491, 0.013258014, 0.018528495, 0.020053765, + 0.021896167, 0.022047184, 0.018604005, 0.01547796, 0.018075445, 0.021503523, 0.022349217, + 0.025097721, 0.032875078, 0.030594727, 0.032573045, 0.034717482, 0.034792991, + 0.033267721, 0.042887485, 0.036846816, 0.03876473, 0.045016819, 0.040063473, 0.04883754, + 0.043944602, 0.048142864, 0.041588741, 0.048762031, 0.027921732, 0.033871788, + 0.022122693, 0.016927718, 0.008833228, 0.00478598, 0.002807662)), + kernel_bandwidth=0.1 +) + + +# From CERN-OPEN-2021-04 and refererences therein +virus_distributions = { + 'SARS_CoV_2': mc.SARSCoV2( + viral_load_in_sputum=log_symptomatic_vl_frequencies, + quantum_infectious_dose=100, + ), + 'SARS_CoV_2_B117': mc.SARSCoV2( + viral_load_in_sputum=log_symptomatic_vl_frequencies, + quantum_infectious_dose=60, + ), + 'SARS_CoV_2_P1': mc.SARSCoV2( + viral_load_in_sputum=log_symptomatic_vl_frequencies, + quantum_infectious_dose=100/2.25, + ), +} From 944db9bf8acfdb869c78d47c8c80c21b55aec616 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 3 Jun 2021 10:47:24 +0200 Subject: [PATCH 57/83] Minor (transparent) fix in some variable name --- cara/monte_carlo/data.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index ad1f4088..79e25517 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -24,7 +24,7 @@ activity_distributions = { # From CERN-OPEN-2021-04 and refererences therein -log_symptomatic_vl_frequencies = LogCustomKernel( +symptomatic_vl_frequencies = LogCustomKernel( np.array((2.46032, 2.67431, 2.85434, 3.06155, 3.25856, 3.47256, 3.66957, 3.85979, 4.09927, 4.27081, 4.47631, 4.66653, 4.87204, 5.10302, 5.27456, 5.46478, 5.6533, 5.88428, 6.07281, 6.30549, 6.48552, 6.64856, 6.85407, 7.10373, 7.30075, 7.47229, 7.66081, 7.85782, 8.05653, 8.27053, @@ -42,15 +42,15 @@ log_symptomatic_vl_frequencies = LogCustomKernel( # From CERN-OPEN-2021-04 and refererences therein virus_distributions = { 'SARS_CoV_2': mc.SARSCoV2( - viral_load_in_sputum=log_symptomatic_vl_frequencies, + viral_load_in_sputum=symptomatic_vl_frequencies, quantum_infectious_dose=100, ), 'SARS_CoV_2_B117': mc.SARSCoV2( - viral_load_in_sputum=log_symptomatic_vl_frequencies, + viral_load_in_sputum=symptomatic_vl_frequencies, quantum_infectious_dose=60, ), 'SARS_CoV_2_P1': mc.SARSCoV2( - viral_load_in_sputum=log_symptomatic_vl_frequencies, + viral_load_in_sputum=symptomatic_vl_frequencies, quantum_infectious_dose=100/2.25, ), } From 4fe0fbefbeb36ade9b0356ad42dfeb916022c0a6 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 05:45:25 +0200 Subject: [PATCH 58/83] Modifying vg for gravitational settlement, to comply with the CERN note --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 794ed21b..d92d76e7 100644 --- a/cara/models.py +++ b/cara/models.py @@ -732,7 +732,7 @@ class ConcentrationModel: def infectious_virus_removal_rate(self, time: float) -> _VectorisedFloat: # Particle deposition on the floor - vg = 1 * 10 ** -4 + vg = 1.88e-4 # Height of the emission source to the floor - i.e. mouth/nose (m) h = 1.5 # Deposition rate (h^-1) From cc98bf003b7d2e7f918526595b07a20f1e10827a Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 09:23:48 +0200 Subject: [PATCH 59/83] Fix the tests which depend on the gravitational settlement (tests on concentration, its integration, r0, etc.); set a common seed for all tests --- cara/tests/models/test_exposure_model.py | 12 ++++++------ cara/tests/test_known_quantities.py | 12 ++++++------ cara/tests/test_predefined_distributions.py | 2 +- cara/tests/test_sampleable_distribution.py | 2 +- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 35c5ac0c..1c6c6d3d 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -136,12 +136,12 @@ def conc_model(): # expected quanta were computed with a trapezoidal integration, using # a mesh of 10'000 pts per exposed presence interval. @pytest.mark.parametrize("exposed_time_interval, expected_quanta", [ - [(0, 1), 5.4869151], - [(1, 1.01), 0.064013521], - [(1.01, 1.02), 0.062266596], - [(12, 12.01), 0.0019025904], - [(12, 24), 78.190763], - [(0, 24), 84.866592], + [(0, 1), 5.3334352], + [(1, 1.01), 0.061759078], + [(1.01, 1.02), 0.060016487], + [(12, 12.01), 0.0019012647], + [(12, 24), 75.513005], + [(0, 24), 81.956988], ] ) def test_exposure_model_integral_accuracy(exposed_time_interval, diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 4af1c06f..ab24bb08 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -44,7 +44,7 @@ def test_concentrations(baseline_model): concentrations = [baseline_model.concentration(t) for t in ts] npt.assert_allclose( concentrations, - [0.000000e+00, 0.4189594, 1.6422648e-14, 0.4189594, 6.4374587e-28], + [0.000000e+00, 0.41611256, 1.3205628e-14, 0.41611256, 4.1909001e-28], rtol=1e-6 ) @@ -91,7 +91,7 @@ def test_r0(baseline_exposure_model): # expected r0 was computed with a trapezoidal integration, using # a mesh of 100'000 pts per exposed presence interval. r0 = baseline_exposure_model.reproduction_number() - npt.assert_allclose(r0, 973.535888) + npt.assert_allclose(r0, 972.880852) def test_periodic_window(baseline_periodic_window, baseline_room): @@ -359,8 +359,8 @@ def build_exposure_model(concentration_model): @pytest.mark.parametrize( "month, expected_quanta", [ - ['Jan', 10.136783], - ['Jun', 41.800377], + ['Jan', 9.930854], + ['Jun', 37.962708], ], ) def test_quanta_hourly_dep(month,expected_quanta): @@ -379,8 +379,8 @@ def test_quanta_hourly_dep(month,expected_quanta): @pytest.mark.parametrize( "month, expected_quanta", [ - ['Jan', 10.19818], - ['Jun', 44.130683], + ['Jan', 9.989881], + ['Jun', 39.99636], ], ) def test_quanta_hourly_dep_refined(month,expected_quanta): diff --git a/cara/tests/test_predefined_distributions.py b/cara/tests/test_predefined_distributions.py index 1e066432..4197b2fb 100644 --- a/cara/tests/test_predefined_distributions.py +++ b/cara/tests/test_predefined_distributions.py @@ -5,7 +5,7 @@ import pytest from cara.monte_carlo.data import activity_distributions, virus_distributions # TODO: seed better the random number generators -np.random.seed(0) +np.random.seed(2000) # mean & std deviations from CERN-OPEN-2021-04 (Table 4) diff --git a/cara/tests/test_sampleable_distribution.py b/cara/tests/test_sampleable_distribution.py index 585b4dde..c30fcbf7 100644 --- a/cara/tests/test_sampleable_distribution.py +++ b/cara/tests/test_sampleable_distribution.py @@ -53,7 +53,7 @@ def test_lognormal(mean_gaussian, std_gaussian): assert len(samples) == sample_size npt.assert_allclose([samples.mean(), samples.std()], [exact_mean, exact_std], rtol=0.01) - npt.assert_allclose(selected_histogram, exact_dist, rtol=0.02) + npt.assert_allclose(selected_histogram, exact_dist, rtol=0.03) @pytest.mark.parametrize( From 3d06228366e977f9c4bbb48f246f5826dad8769d Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 09:24:38 +0200 Subject: [PATCH 60/83] Adding end-to-end integration tests of the full Monte-Carlo approach (comparison with CERN report cases and a few more) --- cara/tests/test_monte_carlo_full_models.py | 296 +++++++++++++++++++++ 1 file changed, 296 insertions(+) create mode 100644 cara/tests/test_monte_carlo_full_models.py diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py new file mode 100644 index 00000000..8a43afd6 --- /dev/null +++ b/cara/tests/test_monte_carlo_full_models.py @@ -0,0 +1,296 @@ +import numpy as np +import numpy.testing as npt +import pytest + +import cara.monte_carlo as mc +from cara import models,data +from cara.monte_carlo.data import activity_distributions, virus_distributions + +# TODO: seed better the random number generators +np.random.seed(2000) +SAMPLE_SIZE = 20000 +TOLERANCE = 0.05 + +# references values for infection_probability and expected new cases +# in the following tests, were obtained from the feature/mc branch + +def test_shared_office(): + # corresponds to the 1st line of Table 5 in CERN-OPEN-2021-04, but + # speaking 30% of the time (instead of 1/3) + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=50, humidity=0.3), + ventilation=models.MultipleVentilation( + ( + models.SlidingWindow( + active=models.PeriodicInterval(period=120, duration=10), + inside_temp=models.PiecewiseConstant((0, 24), (293,)), + outside_temp=models.PiecewiseConstant((0, 24), (283,)), + window_height=1.6, opening_length=0.6, + ), + models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=0.25, + ), + ), + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2_B117'], + 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)), + ), + ) + exposure_mc = mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=3, + presence=concentration_mc.infected.presence, + activity=models.Activity.types['Seated'], + mask=concentration_mc.infected.mask, + ), + ) + exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) + npt.assert_allclose(exposure_model.infection_probability().mean(), + 10.7, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.expected_new_cases().mean(), + 0.32, rtol=TOLERANCE) + + +def test_classroom(): + # corresponds to the 2nd line of Table 5 in CERN-OPEN-2021-04 + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=160, humidity=0.3), + ventilation=models.MultipleVentilation( + ( + models.SlidingWindow( + active=models.PeriodicInterval(period=120, duration=10), + inside_temp=models.PiecewiseConstant((0, 24), (293,)), + outside_temp=models.PiecewiseConstant((0, 24), (283,)), + window_height=1.6, opening_length=0.6, + ), + models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=0.25, + ), + ), + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2_B117'], + 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'], + ), + ) + exposure_mc = mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=19, + presence=concentration_mc.infected.presence, + activity=models.Activity.types['Seated'], + mask=concentration_mc.infected.mask, + ), + ) + exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) + npt.assert_allclose(exposure_model.infection_probability().mean(), + 36.1, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.expected_new_cases().mean(), + 6.87, rtol=TOLERANCE) + + +def test_skicabin(): + # corresponds to the 3rd line of Table 5 in CERN-OPEN-2021-04 + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=10, humidity=0.5), + ventilation=models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=0, + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2_B117'], + presence=mc.SpecificInterval(((0, 1/3),)), + mask=models.Mask(η_inhale=0.3), + activity=activity_distributions['Moderate activity'], + expiration=models.Expiration.types['Talking'], + ), + ) + exposure_mc = mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=3, + presence=concentration_mc.infected.presence, + activity=models.Activity.types['Moderate activity'], + mask=concentration_mc.infected.mask, + ), + ) + exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) + npt.assert_allclose(exposure_model.infection_probability().mean(), + 16.2, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.expected_new_cases().mean(), + 0.49, rtol=TOLERANCE) + + +def test_gym(): + # corresponds to the 4th line of Table 5 in CERN-OPEN-2021-04, + # but there the expected number of cases is wrongly reported as 0.56 + # while it should be 0.63. + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=300, humidity=0.5), + ventilation=models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=6, + ), + infected=mc.InfectedPopulation( + number=2, + virus=virus_distributions['SARS_CoV_2_B117'], + presence=mc.SpecificInterval(((0, 1),)), + mask=models.Mask.types["No mask"], + activity=activity_distributions['Heavy exercise'], + expiration=models.Expiration.types['Breathing'], + ), + ) + exposure_mc = mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=28, + presence=concentration_mc.infected.presence, + activity=models.Activity.types['Heavy exercise'], + mask=concentration_mc.infected.mask, + ), + ) + exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) + npt.assert_allclose(exposure_model.infection_probability().mean(), + 2.2, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.expected_new_cases().mean(), + 0.63, rtol=TOLERANCE) + + +def test_waiting_room(): + # corresponds to the 5th line of Table 5 in CERN-OPEN-2021-04, but + # speaking 30% of the time (instead of 20%) + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=100, humidity=0.5), + ventilation=models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=0.25, + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2_B117'], + 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)), + ), + ) + exposure_mc = mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=14, + presence=concentration_mc.infected.presence, + activity=models.Activity.types['Seated'], + mask=concentration_mc.infected.mask, + ), + ) + exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) + npt.assert_allclose(exposure_model.infection_probability().mean(), + 9.7, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.expected_new_cases().mean(), + 1.36, rtol=TOLERANCE) + + +def test_skagit_chorale(): + # corresponds to the 6th line of Table 5 in CERN-OPEN-2021-04, but + # infection probability should be 29.8% instead of 32%, and number of + # new cases 17.9 instead of 21. + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=810, humidity=0.5), + ventilation=models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=0.7, + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2'], + presence=mc.SpecificInterval(((0, 2.5),)), + mask=models.Mask.types["No mask"], + activity=activity_distributions['Light activity'], + expiration=models.Expiration((5., 5., 5.)), + ), + ) + exposure_mc = mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=60, + presence=concentration_mc.infected.presence, + activity=models.Activity.types['Moderate activity'], + mask=concentration_mc.infected.mask, + ), + ) + exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) + npt.assert_allclose(exposure_model.infection_probability().mean(), + 29.8, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.expected_new_cases().mean(), + 17.9, rtol=TOLERANCE) + + +@pytest.mark.parametrize( + "mask_type, month, expected_probability", + [ + ["No mask", "Jul", 30], + ["Type I", "Jul", 10.2], + ["FFP2", "Jul", 4], + ["Type I", "Feb", 4.3], + ], +) +def test_small_shared_office_Geneva(mask_type,month,expected_probability): + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=33, humidity=0.5), + ventilation=models.MultipleVentilation( + ( + models.SlidingWindow( + active=models.SpecificInterval(((0,24),)), + inside_temp=models.PiecewiseConstant((0, 24), (293,)), + outside_temp=data.GenevaTemperatures[month], + window_height=1.5, opening_length=0.2, + ), + models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=0.25, + ), + ), + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2_B117'], + 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)), + ), + ) + exposure_mc = mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=1, + presence=concentration_mc.infected.presence, + activity=activity_distributions['Seated'], + mask=concentration_mc.infected.mask, + ), + ) + exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) + npt.assert_allclose(exposure_model.infection_probability().mean(), + expected_probability, rtol=TOLERANCE) From 82667698d56a845f9385ebebed2eb351cf5c6dc5 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 10:24:02 +0200 Subject: [PATCH 61/83] Adding comment to refer to CERN report for the gravitational settlement --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index d92d76e7..1561b09d 100644 --- a/cara/models.py +++ b/cara/models.py @@ -731,7 +731,7 @@ class ConcentrationModel: return self.infected.virus def infectious_virus_removal_rate(self, time: float) -> _VectorisedFloat: - # Particle deposition on the floor + # Particle deposition on the floor (value from CERN-OPEN-2021-04) vg = 1.88e-4 # Height of the emission source to the floor - i.e. mouth/nose (m) h = 1.5 From 2f61cf1a967570c550b04f4f388d0126ff7d4ab1 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 11:58:23 +0200 Subject: [PATCH 62/83] Adding more quantities to test in integration tests (qR, dose); improving test structure --- cara/tests/test_monte_carlo_full_models.py | 144 ++++++++++++--------- 1 file changed, 82 insertions(+), 62 deletions(-) diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 8a43afd6..7c08757f 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -8,15 +8,18 @@ from cara.monte_carlo.data import activity_distributions, virus_distributions # TODO: seed better the random number generators np.random.seed(2000) -SAMPLE_SIZE = 20000 +SAMPLE_SIZE = 50000 TOLERANCE = 0.05 # references values for infection_probability and expected new cases # in the following tests, were obtained from the feature/mc branch -def test_shared_office(): - # corresponds to the 1st line of Table 5 in CERN-OPEN-2021-04, but - # speaking 30% of the time (instead of 1/3) +@pytest.fixture +def shared_office_mc(): + """ + Corresponds to the 1st line of Table 5 in CERN-OPEN-2021-04, but + speaking 30% of the time (instead of 1/3) + """ concentration_mc = mc.ConcentrationModel( room=models.Room(volume=50, humidity=0.3), ventilation=models.MultipleVentilation( @@ -45,7 +48,7 @@ def test_shared_office(): weights=(0.3, 0.7)), ), ) - exposure_mc = mc.ExposureModel( + return mc.ExposureModel( concentration_model=concentration_mc, exposed=mc.Population( number=3, @@ -54,15 +57,13 @@ def test_shared_office(): mask=concentration_mc.infected.mask, ), ) - exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) - npt.assert_allclose(exposure_model.infection_probability().mean(), - 10.7, rtol=TOLERANCE) - npt.assert_allclose(exposure_model.expected_new_cases().mean(), - 0.32, rtol=TOLERANCE) -def test_classroom(): - # corresponds to the 2nd line of Table 5 in CERN-OPEN-2021-04 +@pytest.fixture +def classroom_mc(): + """ + Corresponds to the 2nd line of Table 5 in CERN-OPEN-2021-04 + """ concentration_mc = mc.ConcentrationModel( room=models.Room(volume=160, humidity=0.3), ventilation=models.MultipleVentilation( @@ -88,7 +89,7 @@ def test_classroom(): expiration=models.Expiration.types['Talking'], ), ) - exposure_mc = mc.ExposureModel( + return mc.ExposureModel( concentration_model=concentration_mc, exposed=mc.Population( number=19, @@ -97,15 +98,13 @@ def test_classroom(): mask=concentration_mc.infected.mask, ), ) - exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) - npt.assert_allclose(exposure_model.infection_probability().mean(), - 36.1, rtol=TOLERANCE) - npt.assert_allclose(exposure_model.expected_new_cases().mean(), - 6.87, rtol=TOLERANCE) -def test_skicabin(): - # corresponds to the 3rd line of Table 5 in CERN-OPEN-2021-04 +@pytest.fixture +def ski_cabin_mc(): + """ + Corresponds to the 3rd line of Table 5 in CERN-OPEN-2021-04 + """ concentration_mc = mc.ConcentrationModel( room=models.Room(volume=10, humidity=0.5), ventilation=models.AirChange( @@ -121,7 +120,7 @@ def test_skicabin(): expiration=models.Expiration.types['Talking'], ), ) - exposure_mc = mc.ExposureModel( + return mc.ExposureModel( concentration_model=concentration_mc, exposed=mc.Population( number=3, @@ -130,17 +129,15 @@ def test_skicabin(): mask=concentration_mc.infected.mask, ), ) - exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) - npt.assert_allclose(exposure_model.infection_probability().mean(), - 16.2, rtol=TOLERANCE) - npt.assert_allclose(exposure_model.expected_new_cases().mean(), - 0.49, rtol=TOLERANCE) -def test_gym(): - # corresponds to the 4th line of Table 5 in CERN-OPEN-2021-04, - # but there the expected number of cases is wrongly reported as 0.56 - # while it should be 0.63. +@pytest.fixture +def gym_mc(): + """ + Corresponds to the 4th line of Table 5 in CERN-OPEN-2021-04, + but there the expected number of cases is wrongly reported as 0.56 + while it should be 0.63. + """ concentration_mc = mc.ConcentrationModel( room=models.Room(volume=300, humidity=0.5), ventilation=models.AirChange( @@ -156,7 +153,7 @@ def test_gym(): expiration=models.Expiration.types['Breathing'], ), ) - exposure_mc = mc.ExposureModel( + return mc.ExposureModel( concentration_model=concentration_mc, exposed=mc.Population( number=28, @@ -165,16 +162,14 @@ def test_gym(): mask=concentration_mc.infected.mask, ), ) - exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) - npt.assert_allclose(exposure_model.infection_probability().mean(), - 2.2, rtol=TOLERANCE) - npt.assert_allclose(exposure_model.expected_new_cases().mean(), - 0.63, rtol=TOLERANCE) -def test_waiting_room(): - # corresponds to the 5th line of Table 5 in CERN-OPEN-2021-04, but - # speaking 30% of the time (instead of 20%) +@pytest.fixture +def waiting_room_mc(): + """ + Corresponds to the 5th line of Table 5 in CERN-OPEN-2021-04, but + speaking 30% of the time (instead of 20%) + """ concentration_mc = mc.ConcentrationModel( room=models.Room(volume=100, humidity=0.5), ventilation=models.AirChange( @@ -193,7 +188,7 @@ def test_waiting_room(): weights=(0.3, 0.7)), ), ) - exposure_mc = mc.ExposureModel( + return mc.ExposureModel( concentration_model=concentration_mc, exposed=mc.Population( number=14, @@ -202,17 +197,15 @@ def test_waiting_room(): mask=concentration_mc.infected.mask, ), ) - exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) - npt.assert_allclose(exposure_model.infection_probability().mean(), - 9.7, rtol=TOLERANCE) - npt.assert_allclose(exposure_model.expected_new_cases().mean(), - 1.36, rtol=TOLERANCE) -def test_skagit_chorale(): - # corresponds to the 6th line of Table 5 in CERN-OPEN-2021-04, but - # infection probability should be 29.8% instead of 32%, and number of - # new cases 17.9 instead of 21. +@pytest.fixture +def skagit_chorale_mc(): + """ + Corresponds to the 6th line of Table 5 in CERN-OPEN-2021-04, but + infection probability should be 29.8% instead of 32%, and number of + new cases 17.9 instead of 21. + """ concentration_mc = mc.ConcentrationModel( room=models.Room(volume=810, humidity=0.5), ventilation=models.AirChange( @@ -228,7 +221,7 @@ def test_skagit_chorale(): expiration=models.Expiration((5., 5., 5.)), ), ) - exposure_mc = mc.ExposureModel( + return mc.ExposureModel( concentration_model=concentration_mc, exposed=mc.Population( number=60, @@ -237,23 +230,45 @@ def test_skagit_chorale(): mask=concentration_mc.infected.mask, ), ) - exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) - npt.assert_allclose(exposure_model.infection_probability().mean(), - 29.8, rtol=TOLERANCE) - npt.assert_allclose(exposure_model.expected_new_cases().mean(), - 17.9, rtol=TOLERANCE) @pytest.mark.parametrize( - "mask_type, month, expected_probability", + "mc_model, expected_pi, expected_new_cases, expected_dose, expected_qR", [ - ["No mask", "Jul", 30], - ["Type I", "Jul", 10.2], - ["FFP2", "Jul", 4], - ["Type I", "Feb", 4.3], + ["shared_office_mc", 10.7, 0.32, 0.954, 10.9], + ["classroom_mc", 36.1, 6.85, 13.0, 474.4], + ["ski_cabin_mc", 16.3, 0.49, 0.599, 123.4], + ["gym_mc", 2.25, 0.63, 0.01307, 16.4], + ["waiting_room_mc", 9.72, 1.36, 0.571, 58.9], + ["skagit_chorale_mc",29.9, 17.9, 1.90, 1414], + ] +) +def test_report_models(mc_model, expected_pi, expected_new_cases, + expected_dose, expected_qR, request): + mc_model = request.getfixturevalue(mc_model) + exposure_model = mc_model.build_model(size=SAMPLE_SIZE) + npt.assert_allclose(exposure_model.infection_probability().mean(), + expected_pi, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.expected_new_cases().mean(), + expected_new_cases, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.quanta_exposure().mean(), + expected_dose, rtol=TOLERANCE) + npt.assert_allclose( + exposure_model.concentration_model.infected.emission_rate_when_present().mean(), + expected_qR, rtol=TOLERANCE) + + +@pytest.mark.parametrize( + "mask_type, month, expected_pi, expected_dose, expected_qR", + [ + ["No mask", "Jul", 30.0, 6.764, 64.9], + ["Type I", "Jul", 10.2, 1.223, 11.7], + ["FFP2", "Jul", 4.0, 1.223, 11.7], + ["Type I", "Feb", 4.25, 0.357, 11.7], ], ) -def test_small_shared_office_Geneva(mask_type,month,expected_probability): +def test_small_shared_office_Geneva(mask_type, month, expected_pi, + expected_dose, expected_qR): concentration_mc = mc.ConcentrationModel( room=models.Room(volume=33, humidity=0.5), ventilation=models.MultipleVentilation( @@ -293,4 +308,9 @@ def test_small_shared_office_Geneva(mask_type,month,expected_probability): ) exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) npt.assert_allclose(exposure_model.infection_probability().mean(), - expected_probability, rtol=TOLERANCE) + expected_pi, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.quanta_exposure().mean(), + expected_dose, rtol=TOLERANCE) + npt.assert_allclose( + exposure_model.concentration_model.infected.emission_rate_when_present().mean(), + expected_qR, rtol=TOLERANCE) From 84733984031a1b38006caeaa90956cfe001d7f62 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 13:43:35 +0200 Subject: [PATCH 63/83] Using Monte-Carlo approach in calculator; hard-coding the sample size for now --- cara/apps/calculator/model_generator.py | 24 +++++++++++++----------- cara/apps/calculator/report_generator.py | 20 +++++++++++--------- 2 files changed, 24 insertions(+), 20 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index ad65795e..68b060ab 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -8,7 +8,9 @@ import numpy as np from cara import models from cara import data +import cara.monte_carlo as mc from .. import calculator +from cara.monte_carlo.data import activity_distributions, virus_distributions LOG = logging.getLogger(__name__) @@ -20,7 +22,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() - +_SAMPLE_SIZE = 50000 @dataclass class FormData: @@ -275,9 +277,9 @@ class FormData: mask = models.Mask.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] return mask - def infected_population(self) -> models.InfectedPopulation: + def infected_population(self) -> mc.InfectedPopulation: # Initializes the virus - virus = models.Virus.types[self.virus_type] + virus = virus_distributions[self.virus_type] scenario_activity_and_expiration = { 'office': ( @@ -315,12 +317,12 @@ class FormData: } [activity_defn, expiration_defn] = scenario_activity_and_expiration[self.activity_type] - activity = models.Activity.types[activity_defn] + activity = activity_distributions[activity_defn] expiration = build_expiration(expiration_defn) infected_occupants = self.infected_people - infected = models.InfectedPopulation( + infected = mc.InfectedPopulation( number=infected_occupants, virus=virus, presence=self.infected_present_interval(), @@ -330,7 +332,7 @@ class FormData: ) return infected - def exposed_population(self) -> models.Population: + def exposed_population(self) -> mc.Population: scenario_activity = { 'office': 'Seated', 'controlroom-day': 'Seated', @@ -345,14 +347,14 @@ class FormData: } activity_defn = scenario_activity[self.activity_type] - activity = models.Activity.types[activity_defn] + activity = activity_distributions[activity_defn] infected_occupants = self.infected_people # The number of exposed occupants is the total number of occupants # minus the number of infected occupants. exposed_occupants = self.total_people - infected_occupants - exposed = models.Population( + exposed = mc.Population( number=exposed_occupants, presence=self.exposed_present_interval(), activity=activity, @@ -560,14 +562,14 @@ def model_from_form(form: FormData) -> models.ExposureModel: room = models.Room(volume=volume, humidity=humidity) # Initializes and returns a model with the attributes defined above - return models.ExposureModel( - concentration_model=models.ConcentrationModel( + return mc.ExposureModel( + concentration_model=mc.ConcentrationModel( room=room, ventilation=form.ventilation(), infected=form.infected_population(), ), exposed=form.exposed_population() - ) + ).build_model(size=_SAMPLE_SIZE) def baseline_raw_form_data(): diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index ea22cbb6..6cead28c 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -38,12 +38,13 @@ def calculate_report_data(model: models.ExposureModel): t_start, t_end = model_start_end(model) times = list(np.linspace(t_start, t_end, resolution)) - concentrations = [model.concentration_model.concentration(time) for time in times] + concentrations = [np.mean(model.concentration_model.concentration(time)) + for time in times] highest_const = max(concentrations) - prob = model.infection_probability() - er = model.concentration_model.infected.emission_rate_when_present() + prob = np.mean(model.infection_probability()) + er = np.mean(model.concentration_model.infected.emission_rate_when_present()) exposed_occupants = model.exposed.number - expected_new_cases = model.expected_new_cases() + expected_new_cases = np.mean(model.expected_new_cases()) repeated_events = [] for n in [1, 2, 3, 4, 5]: @@ -52,8 +53,8 @@ def calculate_report_data(model: models.ExposureModel): repeated_events.append( RepeatEvents( repeats=n, - probability_of_infection=repeat_model.infection_probability(), - expected_new_cases=repeat_model.expected_new_cases(), + probability_of_infection=np.mean(repeat_model.infection_probability()), + expected_new_cases=np.mean(repeat_model.expected_new_cases()), ) ) @@ -243,7 +244,8 @@ def comparison_plot(scenarios: typing.Dict[str, models.ExposureModel]): t_start, t_end = model_start_end(model) times = np.linspace(t_start, t_end, resolution) datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in times] - concentrations = [model.concentration_model.concentration(time) for time in times] + concentrations = [np.mean(model.concentration_model.concentration(time)) + for time in times] if name in dash_styled_scenarios: ax.plot(datetimes, concentrations, label=name, linestyle='--') @@ -267,8 +269,8 @@ def comparison_report(scenarios: typing.Dict[str, models.ExposureModel]): statistics = {} for name, model in scenarios.items(): statistics[name] = { - 'probability_of_infection': model.infection_probability(), - 'expected_new_cases': model.expected_new_cases(), + 'probability_of_infection': np.mean(model.infection_probability()), + 'expected_new_cases': np.mean(model.expected_new_cases()), } return { 'plot': img2base64(_figure2bytes(comparison_plot(scenarios))), From 3a8c20a38a840b3a72975cddba0eafacdb9e1756 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 13:45:18 +0200 Subject: [PATCH 64/83] Adding sklearn to requirements.txt --- requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/requirements.txt b/requirements.txt index 1bbecfe1..f2133fa2 100644 --- a/requirements.txt +++ b/requirements.txt @@ -62,6 +62,7 @@ python-dateutil==2.8.1 pyzmq==22.0.3 qrcode==6.1 scipy==1.5.4 +sklearn==0.24.1 Send2Trash==1.5.0 six==1.15.0 sniffio==1.2.0 From 71b59357fad2cd7d1b6eb972f2dc7129c4c9fef2 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 14:16:35 +0200 Subject: [PATCH 65/83] Adding 0.25 ACH to all ventilation schemes in calculator; changing modele_generator tests accordingly --- cara/apps/calculator/model_generator.py | 8 ++++++-- cara/tests/apps/calculator/test_model_generator.py | 10 +++++----- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 68b060ab..1d4f51ec 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -265,11 +265,15 @@ class FormData: ventilation = models.HVACMechanical( active=always_on, q_air_mech=self.air_supply) + # this is a minimal, always present source of ventilation, due + # to the air infiltration from the outside. + # See CERN-OPEN-2021-004, p. 12. + infiltration_ventilation = models.AirChange(active=always_on, air_exch=0.25) if self.hepa_option: hepa = models.HEPAFilter(active=always_on, q_air_mech=self.hepa_amount) - return models.MultipleVentilation((ventilation, hepa)) + return models.MultipleVentilation((ventilation, hepa, infiltration_ventilation)) else: - return ventilation + return models.MultipleVentilation((ventilation, infiltration_ventilation)) def mask(self) -> models.Mask: # Initializes the mask type if mask wearing is "continuous", otherwise instantiates the mask attribute as diff --git a/cara/tests/apps/calculator/test_model_generator.py b/cara/tests/apps/calculator/test_model_generator.py index eb216727..53045564 100644 --- a/cara/tests/apps/calculator/test_model_generator.py +++ b/cara/tests/apps/calculator/test_model_generator.py @@ -50,7 +50,7 @@ def test_ventilation_slidingwindow(baseline_form: model_generator.FormData): baseline_form.opening_distance = 0.6 ts = np.linspace(8, 16, 100) - np.testing.assert_allclose([window.air_exchange(room, t) for t in ts], + np.testing.assert_allclose([window.air_exchange(room, t)+0.25 for t in ts], [baseline_form.ventilation().air_exchange(room, t) for t in ts]) @@ -73,7 +73,7 @@ def test_ventilation_hingedwindow(baseline_form: model_generator.FormData): baseline_form.opening_distance = 0.6 ts = np.linspace(8, 16, 100) - np.testing.assert_allclose([window.air_exchange(room, t) for t in ts], + np.testing.assert_allclose([window.air_exchange(room, t)+0.25 for t in ts], [baseline_form.ventilation().air_exchange(room, t) for t in ts]) @@ -88,7 +88,7 @@ def test_ventilation_mechanical(baseline_form: model_generator.FormData): baseline_form.air_supply = 500. ts = np.linspace(8, 16, 100) - np.testing.assert_allclose([mech.air_exchange(room, t) for t in ts], + np.testing.assert_allclose([mech.air_exchange(room, t)+0.25 for t in ts], [baseline_form.ventilation().air_exchange(room, t) for t in ts]) @@ -103,7 +103,7 @@ def test_ventilation_airchanges(baseline_form: model_generator.FormData): baseline_form.air_changes = 3. ts = np.linspace(8, 16, 100) - np.testing.assert_allclose([airchange.air_exchange(room, t) for t in ts], + np.testing.assert_allclose([airchange.air_exchange(room, t)+0.25 for t in ts], [baseline_form.ventilation().air_exchange(room, t) for t in ts]) @@ -131,7 +131,7 @@ def test_ventilation_window_hepa(baseline_form: model_generator.FormData): baseline_form.hepa_option = True ts = np.linspace(9, 17, 100) - np.testing.assert_allclose([ventilation.air_exchange(room, t) for t in ts], + np.testing.assert_allclose([ventilation.air_exchange(room, t)+0.25 for t in ts], [baseline_form.ventilation().air_exchange(room, t) for t in ts]) From f672ee4208225b37c9bd8994db1eabf283ee538f Mon Sep 17 00:00:00 2001 From: Andre Henriques Date: Fri, 4 Jun 2021 12:25:29 +0000 Subject: [PATCH 66/83] change CERN theme text for COVID measures --- .../calculator/themes/cern/templates/calculator.report.html.j2 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/apps/calculator/themes/cern/templates/calculator.report.html.j2 b/cara/apps/calculator/themes/cern/templates/calculator.report.html.j2 index 4d6769f8..85d8860b 100644 --- a/cara/apps/calculator/themes/cern/templates/calculator.report.html.j2 +++ b/cara/apps/calculator/themes/cern/templates/calculator.report.html.j2 @@ -3,7 +3,7 @@ {% block report_preamble %}

Applicable rules:
- Please ensure that this scenario conforms to current CERN HSE rules (minimum ventilation requirements, mask wearing and the maximum number of people permitted in a space).
+ Please ensure that this scenario conforms to current COVID-related Health & Safety requirements, under the applicable COVID Scale and measures in force at the time of the CARA assessment.
The results of this simulation are colour coded according to the risk values authorized at CERN (approved in December 2020):

  • Events with a P(i) less than 5% may go ahead without further mitigation measures.
  • Events with a P(i) between 5% and 15% shall be subject to ALARA principles (see footnote) to minimise the risk before proceeding.
    • From c7c2969bb891d2329df83553554f3bf15b55971f Mon Sep 17 00:00:00 2001 From: Andre Henriques Date: Fri, 4 Jun 2021 12:38:36 +0000 Subject: [PATCH 67/83] update of report text --- cara/apps/calculator/templates/base/calculator.report.html.j2 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index 9b618988..2ad82bfa 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -209,7 +209,8 @@

    Results:

    {% block report_summary %} - In this scenario, the estimated probability of one exposed occupant getting infected P(i) is {{ prob_inf | non_zero_percentage }} and the expected number of new cases is {{ expected_new_cases | float_format }}. + Taking into account the uncertainties tied to the model variables, in this scenario, the probability of one exposed occupant getting infected P(i) is {{ prob_inf | non_zero_percentage }}[*] and the expected number of new cases is {{ expected_new_cases | float_format }}. +

    [*] The results are based on the parameters and assumptions published in the CERN Open Report CERN-OPEN-2021-004

    {% endblock report_summary %}

    Exposure graph:

    From 1d34db3a81e7fb99ce08dfaa1ca1efcd1766c0bf Mon Sep 17 00:00:00 2001 From: Andre Henriques Date: Fri, 4 Jun 2021 12:40:37 +0000 Subject: [PATCH 68/83] delete repeated events section --- .../templates/base/calculator.report.html.j2 | 24 ------------------- 1 file changed, 24 deletions(-) diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index 2ad82bfa..4896ebf0 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -216,30 +216,6 @@

    Exposure graph:

    -

    Repeated events:

    -

    - The P(i) and expected number of new cases if repeating this scenario event - provided the infected person emits the same amount of viruses each day and the exposed person is subject to the same daily exposure time: - - - - - - - - - - - {% for repeat_event in repeated_events %} - - - - - - {% endfor %} - -
    # of repeated eventsP(i)Expected new cases
    {{ repeat_event.repeats }}{{ repeat_event.probability_of_infection | non_zero_percentage }}{{ repeat_event.expected_new_cases | float_format }}
    -

    -

    Alternative scenarios:

    From 8d55eaa632045ddb7c20299197562e551166a9f3 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 14:55:22 +0200 Subject: [PATCH 69/83] Improving the reference values used for the test on the temperature time discretization (using a very fine mesh to compute the reference); decreasing the level of accuracy required for this test --- cara/tests/test_known_quantities.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index ab24bb08..d8597cf8 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -375,12 +375,13 @@ def test_quanta_hourly_dep(month,expected_quanta): npt.assert_allclose(quanta, expected_quanta) # expected quanta were computed with a trapezoidal integration, using -# a mesh of 100'000 pts per exposed presence interval. +# a mesh of 100'000 pts per exposed presence interval and 25 pts per hour +# for the temperature discretization. @pytest.mark.parametrize( "month, expected_quanta", [ - ['Jan', 9.989881], - ['Jun', 39.99636], + ['Jan', 9.993842], + ['Jun', 40.151985], ], ) def test_quanta_hourly_dep_refined(month,expected_quanta): @@ -393,4 +394,4 @@ def test_quanta_hourly_dep_refined(month,expected_quanta): ) ) quanta = m.quanta_exposure() - npt.assert_allclose(quanta, expected_quanta) + npt.assert_allclose(quanta, expected_quanta, rtol=0.02) From 50a36140ab87d76d70376427d80190dde6e2fb18 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 15:00:02 +0200 Subject: [PATCH 70/83] Decreasing (still to an acceptable level) the accuracy of the mode - less points per hour for the temperature, and less sample points in model_generator --- cara/apps/calculator/model_generator.py | 2 +- cara/data.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 1d4f51ec..4baddaf7 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -22,7 +22,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() -_SAMPLE_SIZE = 50000 +_SAMPLE_SIZE = 10000 @dataclass class FormData: diff --git a/cara/data.py b/cara/data.py index f84e627e..61f103b9 100644 --- a/cara/data.py +++ b/cara/data.py @@ -38,7 +38,7 @@ GenevaTemperatures_hourly = { } # same temperatures on a finer temperature mesh GenevaTemperatures = { - month: GenevaTemperatures_hourly[month].refine(refine_factor=10) + month: GenevaTemperatures_hourly[month].refine(refine_factor=4) for month,temperatures in Geneva_hourly_temperatures_celsius_per_hour.items() } From 6e223050fbcd4ddd57ebe9137c66ab72970650b1 Mon Sep 17 00:00:00 2001 From: Andre Henriques Date: Fri, 4 Jun 2021 13:08:03 +0000 Subject: [PATCH 71/83] update C(t) plot labels --- cara/apps/calculator/report_generator.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 6cead28c..8ae90534 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -128,13 +128,13 @@ def plot(times, concentrations, model: models.ExposureModel): fig = plt.figure() ax = fig.add_subplot(1, 1, 1) datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in times] - ax.plot(datetimes, concentrations, lw=2, color='#1f77b4', label='Concentration') + ax.plot(datetimes, concentrations, lw=2, color='#1f77b4', label='Mean concentration') ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.set_xlabel('Time of day') - ax.set_ylabel('Concentration ($q/m^3$)') - ax.set_title('Concentration of infectious quanta') + ax.set_ylabel('Mean concentration ($q/m^3$)') + ax.set_title('Mean concentration of infectious quanta') ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) # Plot presence of exposed person From 8dc8e8f42ab64d6bf3e62e8abe35682b4e3c1653 Mon Sep 17 00:00:00 2001 From: Andre Henriques Date: Fri, 4 Jun 2021 13:17:45 +0000 Subject: [PATCH 72/83] delete 'Exposure Graph' title --- .../calculator/templates/base/calculator.report.html.j2 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index 4896ebf0..174a85eb 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -210,11 +210,11 @@

    {% block report_summary %} Taking into account the uncertainties tied to the model variables, in this scenario, the probability of one exposed occupant getting infected P(i) is {{ prob_inf | non_zero_percentage }}[*] and the expected number of new cases is {{ expected_new_cases | float_format }}. -

    [*] The results are based on the parameters and assumptions published in the CERN Open Report CERN-OPEN-2021-004

    +

    [*] The results are based on the parameters and assumptions published in the CERN Open Report CERN-OPEN-2021-004

    {% endblock report_summary %}

    -

    Exposure graph:

    - + +

    Alternative scenarios:

    From 30923e490fc4ecf0d0942a1dd8d0f3e1431b1137 Mon Sep 17 00:00:00 2001 From: Andre Henriques Date: Fri, 4 Jun 2021 13:25:52 +0000 Subject: [PATCH 73/83] text edit in alternative scenario notes --- cara/apps/calculator/templates/base/calculator.report.html.j2 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index 174a85eb..58a4096e 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -246,7 +246,7 @@

    Notes for alternative scenarios:

    1. This graph shows the concentration of infectious quanta in the air. The filtration of Type I and FFP2 masks, if worn, applies not only to the emission rate but also to the individual exposure (i.e. inhalation). - For this reason, scenarios with different types of mask will show the same concentration on the graph but have different Pi values.
      2. + For this reason, scenarios with different types of mask will show the same concentration on the graph but have different absorbed doses and infection probabilities.
    2. If you have selected more sophisticated options, such as HEPA filtration or FFP2 masks, this will be indicated in the plot as the "base scenario", representing the inputs inserted in the form.
      The other alternative scenarios shown for comparison will not include either HEPA filtration or FFP2 masks.
    From 4241b1c061f4efb7847c43133f7592a9e2c76402 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 17:53:50 +0200 Subject: [PATCH 74/83] Increasing the timeout in tests to 20 seconds; reverting to more accurate sample size & temperature refinement factor (resp. 50000 and 10); changing a few strings in the tests on the report text, that was modified by previous commits --- cara/apps/calculator/model_generator.py | 2 +- cara/data.py | 2 +- .../tests/apps/calculator/test_report_generator.py | 2 +- cara/tests/apps/calculator/test_webapp.py | 14 ++++++++------ cara/tests/test_known_quantities.py | 2 +- 5 files changed, 12 insertions(+), 10 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 4baddaf7..1d4f51ec 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -22,7 +22,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() -_SAMPLE_SIZE = 10000 +_SAMPLE_SIZE = 50000 @dataclass class FormData: diff --git a/cara/data.py b/cara/data.py index 61f103b9..f84e627e 100644 --- a/cara/data.py +++ b/cara/data.py @@ -38,7 +38,7 @@ GenevaTemperatures_hourly = { } # same temperatures on a finer temperature mesh GenevaTemperatures = { - month: GenevaTemperatures_hourly[month].refine(refine_factor=4) + month: GenevaTemperatures_hourly[month].refine(refine_factor=10) for month,temperatures in Geneva_hourly_temperatures_celsius_per_hour.items() } diff --git a/cara/tests/apps/calculator/test_report_generator.py b/cara/tests/apps/calculator/test_report_generator.py index da648635..75733d5f 100644 --- a/cara/tests/apps/calculator/test_report_generator.py +++ b/cara/tests/apps/calculator/test_report_generator.py @@ -11,7 +11,7 @@ def test_generate_report(baseline_form): # generate a report for it. Because this is what happens in the cara # calculator, we confirm that the generation happens within a reasonable # time threshold. - time_limit: float = 5.0 # seconds + time_limit: float = 20.0 # seconds start = time.perf_counter() diff --git a/cara/tests/apps/calculator/test_webapp.py b/cara/tests/apps/calculator/test_webapp.py index d8a8d6d7..9ee1a703 100644 --- a/cara/tests/apps/calculator/test_webapp.py +++ b/cara/tests/apps/calculator/test_webapp.py @@ -45,11 +45,12 @@ class TestBasicApp(tornado.testing.AsyncHTTPTestCase): def get_app(self): return cara.apps.calculator.make_app() + @tornado.testing.gen_test(timeout=20) def test_report(self): - response = self.fetch('/calculator/baseline-model/result') + response = yield self.http_client.fetch(self.get_url('/calculator/baseline-model/result')) self.assertEqual(response.code, 200) - assert 'CERN HSE rules' not in response.body.decode() - assert 'the expected number of new cases is' in response.body.decode() + assert 'CERN HSE' not in response.body.decode() + assert 'expected number of new cases is' in response.body.decode() class TestCernApp(tornado.testing.AsyncHTTPTestCase): @@ -57,11 +58,12 @@ class TestCernApp(tornado.testing.AsyncHTTPTestCase): cern_theme = Path(cara.apps.calculator.__file__).parent / 'themes' / 'cern' return cara.apps.calculator.make_app(theme_dir=cern_theme) + @tornado.testing.gen_test(timeout=20) def test_report(self): - response = self.fetch('/calculator/baseline-model/result') + response = yield self.http_client.fetch(self.get_url('/calculator/baseline-model/result')) self.assertEqual(response.code, 200) - assert 'CERN HSE rules' in response.body.decode() - assert 'the expected number of new cases is' in response.body.decode() + assert 'CERN HSE' in response.body.decode() + assert 'expected number of new cases is' in response.body.decode() async def test_qrcode_urls(http_server_client, baseline_form): diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index d8597cf8..2577c837 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -394,4 +394,4 @@ def test_quanta_hourly_dep_refined(month,expected_quanta): ) ) quanta = m.quanta_exposure() - npt.assert_allclose(quanta, expected_quanta, rtol=0.02) + npt.assert_allclose(quanta, expected_quanta, rtol=0.01) From 5b74538fb6c1aed393597bab2f2b12c0d238451b Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 4 Jun 2021 18:09:37 +0200 Subject: [PATCH 75/83] Setting timeout to 30 seconds --- cara/tests/apps/calculator/test_report_generator.py | 2 +- cara/tests/apps/calculator/test_webapp.py | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/cara/tests/apps/calculator/test_report_generator.py b/cara/tests/apps/calculator/test_report_generator.py index 75733d5f..44131d3a 100644 --- a/cara/tests/apps/calculator/test_report_generator.py +++ b/cara/tests/apps/calculator/test_report_generator.py @@ -11,7 +11,7 @@ def test_generate_report(baseline_form): # generate a report for it. Because this is what happens in the cara # calculator, we confirm that the generation happens within a reasonable # time threshold. - time_limit: float = 20.0 # seconds + time_limit: float = 30.0 # seconds start = time.perf_counter() diff --git a/cara/tests/apps/calculator/test_webapp.py b/cara/tests/apps/calculator/test_webapp.py index 9ee1a703..fcb8c9b4 100644 --- a/cara/tests/apps/calculator/test_webapp.py +++ b/cara/tests/apps/calculator/test_webapp.py @@ -6,6 +6,7 @@ import tornado.testing import cara.apps.calculator from cara.apps.calculator.report_generator import generate_qr_code +_TIMEOUT = 30. @pytest.fixture def app(): @@ -45,7 +46,7 @@ class TestBasicApp(tornado.testing.AsyncHTTPTestCase): def get_app(self): return cara.apps.calculator.make_app() - @tornado.testing.gen_test(timeout=20) + @tornado.testing.gen_test(timeout=_TIMEOUT) def test_report(self): response = yield self.http_client.fetch(self.get_url('/calculator/baseline-model/result')) self.assertEqual(response.code, 200) @@ -58,7 +59,7 @@ class TestCernApp(tornado.testing.AsyncHTTPTestCase): cern_theme = Path(cara.apps.calculator.__file__).parent / 'themes' / 'cern' return cara.apps.calculator.make_app(theme_dir=cern_theme) - @tornado.testing.gen_test(timeout=20) + @tornado.testing.gen_test(timeout=_TIMEOUT) def test_report(self): response = yield self.http_client.fetch(self.get_url('/calculator/baseline-model/result')) self.assertEqual(response.code, 200) From 56ecf150a2cddc7791eff5ac839c75d603b9b2ce Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sat, 5 Jun 2021 18:36:06 +0200 Subject: [PATCH 76/83] Fixing on Monte-Carlo test (type issue) --- cara/tests/test_monte_carlo.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/cara/tests/test_monte_carlo.py b/cara/tests/test_monte_carlo.py index c3bfc2fe..23e139fd 100644 --- a/cara/tests/test_monte_carlo.py +++ b/cara/tests/test_monte_carlo.py @@ -76,7 +76,9 @@ def test_build_concentration_model(baseline_mc_model: cara.monte_carlo.Concentra model = baseline_mc_model.build_model(7) assert isinstance(model, cara.models.ConcentrationModel) assert isinstance(model.concentration(time=0), float) - assert model.concentration(time=1).shape == (7, ) + conc = model.concentration(time=1) + assert isinstance(conc, np.ndarray) + assert conc.shape == (7, ) def test_build_exposure_model(baseline_mc_exposure_model: cara.monte_carlo.ExposureModel): From 17956d766daba582ea093d9a9a1f9365b3b0ace4 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sat, 5 Jun 2021 18:38:57 +0200 Subject: [PATCH 77/83] Reverting to a refinement factor of 4 in the time-dependence of Geneva temperature; optimizing the caching; decreasing the timeout value in tests --- cara/data.py | 2 +- cara/models.py | 8 ++++---- cara/tests/apps/calculator/test_report_generator.py | 2 +- cara/tests/apps/calculator/test_webapp.py | 2 +- cara/tests/test_known_quantities.py | 2 +- 5 files changed, 8 insertions(+), 8 deletions(-) diff --git a/cara/data.py b/cara/data.py index f84e627e..61f103b9 100644 --- a/cara/data.py +++ b/cara/data.py @@ -38,7 +38,7 @@ GenevaTemperatures_hourly = { } # same temperatures on a finer temperature mesh GenevaTemperatures = { - month: GenevaTemperatures_hourly[month].refine(refine_factor=10) + month: GenevaTemperatures_hourly[month].refine(refine_factor=4) for month,temperatures in Geneva_hourly_temperatures_celsius_per_hour.items() } diff --git a/cara/models.py b/cara/models.py index 1561b09d..ad4493ba 100644 --- a/cara/models.py +++ b/cara/models.py @@ -711,7 +711,6 @@ class InfectedPopulation(Population): return self.emission_rate_when_present() - @cached() def emission_rate(self, time) -> _VectorisedFloat: """ The emission rate of the entire population. @@ -741,7 +740,6 @@ class ConcentrationModel: return k + self.virus.decay_constant(self.room.humidity ) + self.ventilation.air_exchange(self.room, time) - @cached() def _concentration_limit(self, time: float) -> _VectorisedFloat: """ Provides a constant that represents the theoretical asymptotic @@ -753,7 +751,6 @@ class ConcentrationModel: return (self.infected.emission_rate(time)) / (IVRR * V) - @cached() def state_change_times(self): """ All time dependent entities on this model must provide information about @@ -798,6 +795,9 @@ class ConcentrationModel: return (self.last_state_change(stop) <= start) @cached() + def _concentration_at_state_change(self, time: float) -> _VectorisedFloat: + return self.concentration(time) + def concentration(self, time: float) -> _VectorisedFloat: """ Virus quanta concentration, as a function of time. @@ -816,7 +816,7 @@ class ConcentrationModel: concentration_limit = self._concentration_limit(next_state_change_time) t_last_state_change = self.last_state_change(time) - concentration_at_last_state_change = self.concentration(t_last_state_change) + concentration_at_last_state_change = self._concentration_at_state_change(t_last_state_change) delta_time = time - t_last_state_change fac = np.exp(-IVRR * delta_time) diff --git a/cara/tests/apps/calculator/test_report_generator.py b/cara/tests/apps/calculator/test_report_generator.py index 44131d3a..75733d5f 100644 --- a/cara/tests/apps/calculator/test_report_generator.py +++ b/cara/tests/apps/calculator/test_report_generator.py @@ -11,7 +11,7 @@ def test_generate_report(baseline_form): # generate a report for it. Because this is what happens in the cara # calculator, we confirm that the generation happens within a reasonable # time threshold. - time_limit: float = 30.0 # seconds + time_limit: float = 20.0 # seconds start = time.perf_counter() diff --git a/cara/tests/apps/calculator/test_webapp.py b/cara/tests/apps/calculator/test_webapp.py index fcb8c9b4..0a22fe59 100644 --- a/cara/tests/apps/calculator/test_webapp.py +++ b/cara/tests/apps/calculator/test_webapp.py @@ -6,7 +6,7 @@ import tornado.testing import cara.apps.calculator from cara.apps.calculator.report_generator import generate_qr_code -_TIMEOUT = 30. +_TIMEOUT = 20. @pytest.fixture def app(): diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 2577c837..d8597cf8 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -394,4 +394,4 @@ def test_quanta_hourly_dep_refined(month,expected_quanta): ) ) quanta = m.quanta_exposure() - npt.assert_allclose(quanta, expected_quanta, rtol=0.01) + npt.assert_allclose(quanta, expected_quanta, rtol=0.02) From cf817d0eb72e5aab05f6cbb1bb2a4d4e9181d07c Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sat, 5 Jun 2021 19:08:46 +0200 Subject: [PATCH 78/83] Modifying sklearn version in requirements --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index f2133fa2..d7f8affc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -62,7 +62,7 @@ python-dateutil==2.8.1 pyzmq==22.0.3 qrcode==6.1 scipy==1.5.4 -sklearn==0.24.1 +sklearn==0.23.1 Send2Trash==1.5.0 six==1.15.0 sniffio==1.2.0 From 20a288a8876a8a7bf172a80d69b8317b0fa5d215 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sat, 5 Jun 2021 17:18:14 +0000 Subject: [PATCH 79/83] Update requirements.txt --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index d7f8affc..eac04f3e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -62,7 +62,7 @@ python-dateutil==2.8.1 pyzmq==22.0.3 qrcode==6.1 scipy==1.5.4 -sklearn==0.23.1 +scikit_learn==0.23.1 Send2Trash==1.5.0 six==1.15.0 sniffio==1.2.0 From 263330796aa899eeeffde8eea38d1bb39456e946 Mon Sep 17 00:00:00 2001 From: Andre Henriques Date: Sat, 5 Jun 2021 19:34:30 +0000 Subject: [PATCH 80/83] update text format in report --- .../calculator/templates/base/calculator.report.html.j2 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index 58a4096e..cb57d045 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -209,8 +209,8 @@

    Results:

    {% block report_summary %} - Taking into account the uncertainties tied to the model variables, in this scenario, the probability of one exposed occupant getting infected P(i) is {{ prob_inf | non_zero_percentage }}[*] and the expected number of new cases is {{ expected_new_cases | float_format }}. -

    [*] The results are based on the parameters and assumptions published in the CERN Open Report CERN-OPEN-2021-004

    + Taking into account the uncertainties tied to the model variables, in this scenario, the probability of one exposed occupant getting infected is {{ prob_inf | non_zero_percentage }}[*] and the expected number of new cases is {{ expected_new_cases | float_format }}. +

    [*] The results are based on the parameters and assumptions published in the CERN Open Report CERN-OPEN-2021-004

    {% endblock report_summary %}

    @@ -225,7 +225,7 @@ Scenario - P(i) + P(I) Expected new cases From 527a6263970b1c7b8e9e48af164d746b46e932a5 Mon Sep 17 00:00:00 2001 From: Andre Henriques Date: Sat, 5 Jun 2021 19:36:25 +0000 Subject: [PATCH 81/83] update labels in 'alternative scenarios' plot --- cara/apps/calculator/report_generator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 8ae90534..4029e54f 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -259,8 +259,8 @@ def comparison_plot(scenarios: typing.Dict[str, models.ExposureModel]): ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) ax.set_xlabel('Time of day') - ax.set_ylabel('Concentration ($q/m^3$)') - ax.set_title('Concentration of infectious quanta') + ax.set_ylabel('Mean concentration ($q/m^3$)') + ax.set_title('Mean concentration of infectious quanta') return fig From a6bb522de9b601c06b5cb7b125ba72f7f3639af1 Mon Sep 17 00:00:00 2001 From: Andre Henriques Date: Mon, 7 Jun 2021 12:19:44 +0000 Subject: [PATCH 82/83] Update userguide for mc --- .../apps/calculator/templates/userguide.html.j2 | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/cara/apps/calculator/templates/userguide.html.j2 b/cara/apps/calculator/templates/userguide.html.j2 index 013a1763..f519a9bb 100644 --- a/cara/apps/calculator/templates/userguide.html.j2 +++ b/cara/apps/calculator/templates/userguide.html.j2 @@ -8,7 +8,7 @@

    This is a guide to help you use the calculator app. If you are using the expert version of the tool, you should look at the expert notes.

    For more information on the Airborne Transmission of SARS-CoV-2, feel free to check out the HSE Seminar: https://cds.cern.ch/record/2743403

    -

    The methodology, mathematical equations and parameters of the model are described here: https://edms.cern.ch/ui/file/2566402/1/CARA_Deterministic_parameters_2020.pdf

    +

    The methodology, mathematical equations and parameters of the model are described here in the CERN Report: CERN-OPEN-2021-004

    Disclaimer

    @@ -172,26 +172,25 @@ Please check what are the applicable rules, before deciding which assumptions ar Please confirm what are the applicable rules, before deciding which assumptions are used for the simulation

    For the time being only the Type 1 surgical and FFP2 masks can be selected.

    Generate Report

    -

    When you have entered all the necessary information, please click on the Generate Report button to execute the model.

    +

    When you have entered all the necessary information, please click on the Generate Report button to execute the model. With the implementation of Monte Carlo simulations, the browser might take a few secounds to react.

    Report

    The report will open in your web browser. It contains a summary of all the input data, which will allow the simulation to be repeated if required in the future as we improve the model.

    Results

    -

    This part of the report shows the P(i) or probability of one exposed person getting infected. +

    This part of the report shows the P(I) or probability of one exposed person getting infected. It is estimated based on the emission rate of virus into the simulated volume, and the amount which is inhaled by exposed individuals. -This probability is valid for the simulation duration - i.e. if you have simulated one day and plan to work 5 days in these conditions and the infected person emits the same amount of virus each day, the cumulative probability of infection is (1-(1-P(i))^5). +This probability is valid for the simulation duration - i.e. the start and end time. If you are using the natural ventilation option, the simulation is only valid for the selected month, because the following or preceding month will have a different average temperature profile. The expected number of new cases for the simulation is calculated based on the probability of infection, multiplied by the number of exposed occupants.

    -

    Exposure graph

    -

    The graph shows the variation in the concentration of infectious quanta (one quanta is the amount of inhaled virus that can cause infection in 63% of the exposed occupants) within the simulated volume. +

    The graph shows the variation in the concentration of infectious viruses within the simulated volume. It is determined by:

    • The presence of the infected person, who emits airborne viruses in the volume.
      • -
    • The emission rate is related to the type of activity of the infected person (sitting, light exercise), their level of vocalisation (breathing, whispering or talking).
      • +
    • The emission rate is related to the type of activity of the infected person (sitting, light exercise), their level of vocalisation (breathing, talking or shouting).
    • The accumulation of infectious quanta in the volume, which is driven, among other factors, by ventilation (if applicable).
      • In a mechanical ventilation scenario, the removal rate is constant, based on fresh airflow supply in and out of the simulated space.
      • Under natural ventilation conditions, the effectiveness of ventilation relies upon the hourly temperature difference between the inside and outside air temperature.
        • -
      • A HEPA filter removes infectious quanta from the air at a constant rate and is modelled in the same way as mechanical ventilation, however air passed through a HEPA filter is recycled (i.e. it is not fresh air).
        • +
      • A HEPA filter removes infectious virus from the air at a constant rate and is modelled in the same way as mechanical ventilation, however air passed through a HEPA filter is recycled (i.e. it is not fresh air).
    @@ -204,7 +203,7 @@ This allows for:

Conclusion

This tool provides informative comparisons for COVID-19 (long-range) airborne risk only - see Disclaimer -If you have any comments on your experience with the app, or feedback for potential improvements, please share them with the development team at cara-dev@cern.ch.

+If you have any comments on your experience with the app, or feedback for potential improvements, please share them with the development team Send email.

{% endblock contents %} From 6a3751c9e6756bd2bab0724d8168030834be0b01 Mon Sep 17 00:00:00 2001 From: Andre Henriques Date: Mon, 7 Jun 2021 12:22:57 +0000 Subject: [PATCH 83/83] Update about page for mc --- cara/apps/templates/about.html.j2 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/apps/templates/about.html.j2 b/cara/apps/templates/about.html.j2 index 3ff6222d..56d0a307 100644 --- a/cara/apps/templates/about.html.j2 +++ b/cara/apps/templates/about.html.j2 @@ -17,7 +17,7 @@ CARA stands for COVID Airborne Risk Assessment and was developed in the spring o The mathematical and physical model simulate the long-range airborne spread of SARS-CoV-2 virus in a finite volume, assuming a homogenous mixture, and estimates the risk of COVID-19 infection therein. The results DO NOT include short-range airborne exposure (where the physical distance plays a factor) nor the other known modes of SARS-CoV-2 transmission. Hence, the output from this model is only valid when the other recommended public health & safety instructions are observed, such as adequate physical distancing, good hand hygiene and other barrier measures.
-

The methodology, mathematical equations and parameters of the model are described here: https://edms.cern.ch/ui/file/2566402/1/CARA_Deterministic_parameters_2020.pdf.

+

The methodology, mathematical equations and parameters of the model are described here in the CERN Report: CERN-OPEN-2021-004.

The model used is based on scientific publications relating to airborne transmission of infectious diseases, virology, epidemiology and aerosol science. It can be used to compare the effectiveness of different airborne-related risk mitigation measures.