From 8afa8bac9d9812512aabff9e660c27a87872e2c8 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 12 Jan 2022 17:50:27 +0100 Subject: [PATCH 01/29] Stylistic changes --- cara/tests/test_monte_carlo_full_models.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 71a630b0..e5c96b82 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -22,7 +22,7 @@ def shared_office_mc(): """ concentration_mc = mc.ConcentrationModel( room=models.Room(volume=50, humidity=0.5), - ventilation = models.MultipleVentilation( + ventilation=models.MultipleVentilation( ventilations=( models.SlidingWindow( active=models.PeriodicInterval(period=120, duration=120), @@ -36,7 +36,7 @@ def shared_office_mc(): ), infected=mc.InfectedPopulation( number=1, - presence=mc.SpecificInterval(present_times = ((0, 3.5), (4.5, 9))), + presence=mc.SpecificInterval(present_times=((0, 3.5), (4.5, 9))), virus=virus_distributions['SARS_CoV_2_DELTA'], mask=models.Mask.types['No mask'], activity=activity_distributions['Seated'], @@ -49,7 +49,7 @@ def shared_office_mc(): concentration_model=concentration_mc, exposed=mc.Population( number=3, - presence=mc.SpecificInterval(present_times = ((0, 3.5), (4.5, 9))), + presence=mc.SpecificInterval(present_times=((0, 3.5), (4.5, 9))), activity=activity_distributions['Seated'], mask=models.Mask.types['No mask'], host_immunity=0., @@ -64,7 +64,7 @@ def classroom_mc(): """ concentration_mc = mc.ConcentrationModel( room=models.Room(volume=160, humidity=0.3), - ventilation = models.MultipleVentilation( + ventilation=models.MultipleVentilation( ventilations=( models.SlidingWindow( active=models.PeriodicInterval(period=120, duration=120), From e26249327d86581f5716ee2e275231cd75c61172 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Thu, 13 Jan 2022 11:18:02 +0100 Subject: [PATCH 02/29] updated test_r0 assert value --- cara/tests/test_known_quantities.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 9fc10597..6aca5e70 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -95,7 +95,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, 776.9419902161412) + npt.assert_allclose(r0, 776.941990) def test_periodic_window(baseline_periodic_window, baseline_room): From 8d2d846492545135e471df4ff0e69576655a0150 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Thu, 13 Jan 2022 11:23:54 +0100 Subject: [PATCH 03/29] oneoverln2 declaration changed location on models.py --- cara/models.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/cara/models.py b/cara/models.py index 3cedbf22..dcf4eb78 100644 --- a/cara/models.py +++ b/cara/models.py @@ -50,8 +50,6 @@ from .utils import method_cache from .dataclass_utils import nested_replace -oneoverln2 = 1 / np.log(2) - # Define types for items supporting vectorisation. In the future this may be replaced # by ``np.ndarray[]`` once/if that syntax is supported. Note that vectorization # implies 1d arrays: multi-dimensional arrays are not supported. @@ -1109,6 +1107,7 @@ class ExposureModel: ) # oneoverln2 multiplied by ID_50 corresponds to ID_63. + oneoverln2 = 1 / np.log(2) infectious_dose = oneoverln2 * self.concentration_model.virus.infectious_dose # Probability of infection. From ad74d68db41d8e3a40e8bfc55694d4c2d2fc4ace Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Thu, 13 Jan 2022 11:25:49 +0100 Subject: [PATCH 04/29] removed fraction_deposited comment on text_exposure_model --- 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 8b572daf..858b0ada 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -206,7 +206,7 @@ def test_infectious_dose_vectorisation(): 10, presence_interval, models.Mask.types['Type I'], models.Activity.types['Standing'], 0., ) - model = ExposureModel(cm, population) #, fraction_deposited=1.0 + model = ExposureModel(cm, population) inf_probability = model.infection_probability() assert isinstance(inf_probability, np.ndarray) assert inf_probability.shape == (3, ) From f5f0365725772f073281911f0de2490bf5dce618 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 17 Jan 2022 14:01:40 +0100 Subject: [PATCH 05/29] Added cn reference to the code --- cara/monte_carlo/data.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 3b9172d2..290f5920 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -161,7 +161,7 @@ def expiration_distribution(BLO_factors): """ Returns an Expiration with an aerosol diameter distribution, defined by the BLO factors (a length-3 tuple). - The total concentration of aerosols is computed by integrating + The total concentration of aerosols, cn, is computed by integrating the distribution between 0.1 and 30 microns - these boundaries are an historical choice based on previous implementations of the model (it limits the influence of the O-mode). @@ -169,7 +169,7 @@ def expiration_distribution(BLO_factors): dscan = np.linspace(0.1, 30. ,3000) return mc.Expiration(CustomKernel(dscan, BLOmodel(BLO_factors).distribution(dscan),kernel_bandwidth=0.1), - BLOmodel(BLO_factors).integrate(0.1, 30.)) + cn=BLOmodel(BLO_factors).integrate(0.1, 30.)) expiration_BLO_factors = { From 0a6472bdbfc2cbf3f9b599ba382d17b15d5c605f Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 17 Jan 2022 14:09:24 +0100 Subject: [PATCH 06/29] oneoverln2 back on the header of models.py --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index dcf4eb78..24110fa6 100644 --- a/cara/models.py +++ b/cara/models.py @@ -50,6 +50,7 @@ from .utils import method_cache from .dataclass_utils import nested_replace +oneoverln2 = 1 / np.log(2) # Define types for items supporting vectorisation. In the future this may be replaced # by ``np.ndarray[]`` once/if that syntax is supported. Note that vectorization # implies 1d arrays: multi-dimensional arrays are not supported. @@ -1107,7 +1108,6 @@ class ExposureModel: ) # oneoverln2 multiplied by ID_50 corresponds to ID_63. - oneoverln2 = 1 / np.log(2) infectious_dose = oneoverln2 * self.concentration_model.virus.infectious_dose # Probability of infection. From 9da5b5b929447e8d272bae21e241528c882748a5 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Tue, 25 Jan 2022 22:53:41 +0100 Subject: [PATCH 07/29] Improving docstring for cn in Expiration class --- cara/models.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 24110fa6..167d2a67 100644 --- a/cara/models.py +++ b/cara/models.py @@ -589,7 +589,9 @@ class Expiration(_ExpirationBase): #: diameter of the aerosol in microns diameter: _VectorisedFloat - #: total concentration of aerosols (cm^-3) + #: total concentration of aerosols per unit volume of expired air + # (in cm^-3), integrated over all aerosol diameters (corresponding + # to c_n,i in Eq. (4) of https://doi.org/10.1101/2021.10.14.21264988) cn: float = 1. @cached() From c9329da4d1a349cd109003db31c166101c10ebb1 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 26 Jan 2022 13:03:37 +0100 Subject: [PATCH 08/29] models.py: some re-structuring to removing some code-smells (introduction of a Particle class where settling velocity and deposited fraction are computed) --- cara/models.py | 108 +++++++++++++++++++++++++++++++++++-------------- 1 file changed, 78 insertions(+), 30 deletions(-) diff --git a/cara/models.py b/cara/models.py index 167d2a67..88d66555 100644 --- a/cara/models.py +++ b/cara/models.py @@ -563,6 +563,56 @@ Mask.types = { } +@dataclass(frozen=True) +class Particle: + """ + Represents an aerosol particle. + """ + + #: diameter of the aerosol in microns + diameter: typing.Union[None,_VectorisedFloat] = None + + def settling_velocity(self, evaporation_factor: float=0.3) -> _VectorisedFloat: + """ + Settling velocity (i.e. speed of deposition on the floor due + to gravity), for aerosols, in m/s. Diameter-dependent expression + from https://doi.org/10.1101/2021.10.14.21264988 + When an aerosol-diameter is not given, returns + the default value of 1.88e-4 m/s (corresponds to diameter of + 2.5 microns, i.e. geometric average of the breathing + expiration distribution, taking evaporation into account, see + https://doi.org/10.1101/2021.10.14.21264988) + evaporation_factor represents the factor applied to the diameter, + due to instantaneous evaporation of the particle in the air. + """ + if self.diameter is None: + vg = 1.88e-4 + else: + vg = 1.88e-4 * (self.diameter*evaporation_factor / 2.5)**2 + return vg + + def fraction_deposited(self, evaporation_factor: float=0.3): + """ + The fraction of particles actually deposited in the respiratory tract. + From W. C. Hinds, New York, Wiley, 1999 (pp. 233 – 259). + evaporation_factor represents the factor applied to the diameter, + due to instantaneous evaporation of the particle in the air. + """ + if self.diameter is None: + # model is not evaluated for specific values of aerosol + # diameters - we choose a single "average" deposition factor, + # as in https://doi.org/10.1101/2021.10.14.21264988. + fdep = 0.6 + else: + # deposition fraction depends on aerosol particle diameter. + d = (self.diameter * evaporation_factor) + IFrac = 1 - 0.5 * (1 - (1 / (1 + (0.00076*(d**2.8))))) + fdep = IFrac * (0.0587 + + (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d)))) + + (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d))))) + return fdep + + @dataclass(frozen=True) class _ExpirationBase: """ @@ -572,6 +622,12 @@ class _ExpirationBase: #: Pre-populated examples of Expirations. types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]] + def particle(self) -> Particle: + """ + the Particle object representing the aerosol - here the default one + """ + return Particle() + def aerosols(self, mask: Mask): """ total volume of aerosols expired per volume of air (mL/cm^3). @@ -594,6 +650,12 @@ class Expiration(_ExpirationBase): # to c_n,i in Eq. (4) of https://doi.org/10.1101/2021.10.14.21264988) cn: float = 1. + def particle(self) -> Particle: + """ + the Particle object representing the aerosol + """ + return Particle(diameter=self.diameter) + @cached() def aerosols(self, mask: Mask): """ Result is in mL.cm^-3 """ @@ -732,6 +794,13 @@ class _PopulationWithVirus(Population): return self.emission_rate_when_present() + def particle(self) -> Particle: + """ + the Particle object representing the aerosol expired by the + population - here we take the default Particle object + """ + return Particle() + @dataclass(frozen=True) class EmittingPopulation(_PopulationWithVirus): @@ -778,6 +847,12 @@ class InfectedPopulation(_PopulationWithVirus): return ER * self.number + def particle(self) -> Particle: + """ + the Particle object representing the aerosol - here the default one + """ + return self.expiration.particle() + @dataclass(frozen=True) class ConcentrationModel: @@ -796,20 +871,7 @@ class ConcentrationModel: def infectious_virus_removal_rate(self, time: float) -> _VectorisedFloat: # Equilibrium velocity of particle motion toward the floor - if (isinstance(self.infected, InfectedPopulation) - and isinstance(self.infected.expiration, Expiration)): - d = self.infected.expiration.diameter * self.evaporation_factor - vg = 1.88e-4 * (d / 2.5)**2 - # see https://doi.org/10.1101/2021.10.14.21264988 - # (velocity of 1.88e-4 corresponds to diameter of 2.5 microns) - else: - # model is not evaluated for specific values of aerosol - # diameters - we choose a single velocity value - # corresponding to that obtained with a diameter of 2.5 microns - # (geometric average of the breathing expiration distribution, - # taking evaporation into account, see - # https://doi.org/10.1101/2021.10.14.21264988) - vg = 1.88e-4 + vg = self.infected.particle().settling_velocity(self.evaporation_factor) # Height of the emission source to the floor - i.e. mouth/nose (m) h = 1.5 # Deposition rate (h^-1) @@ -990,23 +1052,9 @@ class ExposureModel: def fraction_deposited(self): """ The fraction of viruses actually deposited in the respiratory tract. - From W. C. Hinds, New York, Wiley, 1999 (pp. 233 – 259). """ - if not (isinstance(self.concentration_model.infected,InfectedPopulation) - and isinstance(self.concentration_model.infected.expiration,Expiration)): - # model is not evaluated for specific values of aerosol - # diameters - we choose a single "average" deposition factor, - # as in https://doi.org/10.1101/2021.10.14.21264988. - fdep = 0.6 - else: - # deposition factor depends on aerosol particle diameter. - d = (self.concentration_model.infected.expiration.diameter - * self.concentration_model.evaporation_factor) - IFrac = 1 - 0.5 * (1 - (1 / (1 + (0.00076*(d**2.8))))) - fdep = IFrac * (0.0587 - + (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d)))) - + (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d))))) - return fdep + return self.concentration_model.infected.particle().fraction_deposited( + self.concentration_model.evaporation_factor) def _normed_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: """The number of virions per meter^3 between any two times, normalized From 9b1a324369e91cd781fce92385adefb85a830cf8 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 26 Jan 2022 14:23:31 +0100 Subject: [PATCH 09/29] Minor change in models: particle methods changed into properties --- cara/models.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/cara/models.py b/cara/models.py index 88d66555..d98baf02 100644 --- a/cara/models.py +++ b/cara/models.py @@ -622,6 +622,7 @@ class _ExpirationBase: #: Pre-populated examples of Expirations. types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]] + @property def particle(self) -> Particle: """ the Particle object representing the aerosol - here the default one @@ -650,6 +651,7 @@ class Expiration(_ExpirationBase): # to c_n,i in Eq. (4) of https://doi.org/10.1101/2021.10.14.21264988) cn: float = 1. + @property def particle(self) -> Particle: """ the Particle object representing the aerosol @@ -794,6 +796,7 @@ class _PopulationWithVirus(Population): return self.emission_rate_when_present() + @property def particle(self) -> Particle: """ the Particle object representing the aerosol expired by the @@ -847,11 +850,12 @@ class InfectedPopulation(_PopulationWithVirus): return ER * self.number + @property def particle(self) -> Particle: """ the Particle object representing the aerosol - here the default one """ - return self.expiration.particle() + return self.expiration.particle @dataclass(frozen=True) @@ -871,7 +875,7 @@ class ConcentrationModel: def infectious_virus_removal_rate(self, time: float) -> _VectorisedFloat: # Equilibrium velocity of particle motion toward the floor - vg = self.infected.particle().settling_velocity(self.evaporation_factor) + vg = self.infected.particle.settling_velocity(self.evaporation_factor) # Height of the emission source to the floor - i.e. mouth/nose (m) h = 1.5 # Deposition rate (h^-1) @@ -1053,7 +1057,7 @@ class ExposureModel: """ The fraction of viruses actually deposited in the respiratory tract. """ - return self.concentration_model.infected.particle().fraction_deposited( + return self.concentration_model.infected.particle.fraction_deposited( self.concentration_model.evaporation_factor) def _normed_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: From 10565bd448579efe3e7a3ae7545af3f8d6e8c9a0 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 26 Jan 2022 17:14:52 +0100 Subject: [PATCH 10/29] Updated covid level text --- cara/apps/calculator/report_generator.py | 4 ++-- cara/apps/templates/cern/calculator.report.html.j2 | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index e703e225..ee33fbda 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -309,8 +309,8 @@ class ReportGenerator: context['calculator_prefix'] = self.calculator_prefix context['scale_warning'] = { 'level': 'red-4', - 'incidence_rate': 'higher or equal to 100 new cases per 100 000 inhabitants', - 'onsite_access': 'lower than 4000', + 'risk': 'strong', + 'onsite_access': '4’000', 'threshold': '5%' } return context diff --git a/cara/apps/templates/cern/calculator.report.html.j2 b/cara/apps/templates/cern/calculator.report.html.j2 index 78e73b8f..b01aad11 100644 --- a/cara/apps/templates/cern/calculator.report.html.j2 +++ b/cara/apps/templates/cern/calculator.report.html.j2 @@ -45,19 +45,19 @@
{% if scale_warning.level == "green-1" %} {% elif scale_warning.level == "yellow-2" %} {% elif scale_warning.level == "orange-3" %} {% elif scale_warning.level == "red-4" %} {% else %}

Note: The CERN COVID Level is not specified.

From 6b8f71048acfa8f748d5ed92e7ed97d141d877a2 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 26 Jan 2022 17:46:54 +0100 Subject: [PATCH 11/29] Added covid level variable comments on the code for future reference. --- cara/apps/calculator/report_generator.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index ee33fbda..1a688cc5 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -307,11 +307,13 @@ class ReportGenerator: ) context['permalink'] = generate_permalink(base_url, self.calculator_prefix, form) context['calculator_prefix'] = self.calculator_prefix + + # For further information about these values visit https://gitlab.cern.ch/cara/cara/-/merge_requests/321. context['scale_warning'] = { - 'level': 'red-4', - 'risk': 'strong', - 'onsite_access': '4’000', - 'threshold': '5%' + 'level': 'red-4', # 'red-4' - 'orange-3' - 'yellow-2' - 'green-1' + 'risk': 'strong', # 'strong' - 'medium' - 'reduced' - '' + 'onsite_access': '4’000', # '4’000' - '5’000' - '6’500' - '8’000' + 'threshold': '5%' # '5%' - '' - '' - '' } return context From 5302001ed94f2b7be13ba961154460f8bf1effda Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Wed, 26 Jan 2022 20:47:28 +0100 Subject: [PATCH 12/29] Restructuring exposure class to avoid smelly code (cara/models.py) --- cara/models.py | 146 ++++++++++++++++++++++++++++--------------------- 1 file changed, 84 insertions(+), 62 deletions(-) diff --git a/cara/models.py b/cara/models.py index d09f9f69..35c85dd3 100644 --- a/cara/models.py +++ b/cara/models.py @@ -673,8 +673,8 @@ class MultipleExpiration(_ExpirationBase): Group together different modes of expiration, that represent each the main expiration mode for a certain fraction of time (given by the weights). - This class can only be used with single diameters (it cannot - be used with diameter distributions). + This class can only be used with single diameters defined in each + expiration (it cannot be used with diameter distributions). """ expirations: typing.Tuple[_ExpirationBase, ...] weights: typing.Tuple[float, ...] @@ -768,13 +768,28 @@ class _PopulationWithVirus(Population): """ return 1. + def aerosols(self): + """ + total volume of aerosols expired per volume of air (mL/cm^3). + """ + raise NotImplementedError("Subclass must implement") + + def emission_rate_per_aerosol_when_present(self) -> _VectorisedFloat: + """ + The emission rate of virions per fraction of aerosol volume in + the expired air, if the infected population is present, in cm^3/(mL.h). + It should not be a function of time. + """ + raise NotImplementedError("Subclass must implement") + @method_cache def emission_rate_when_present(self) -> _VectorisedFloat: """ The emission rate if the infected population is present - (in virions / h). It should not be a function of time. + (in virions / h). """ - raise NotImplementedError("Subclass must implement") + return (self.emission_rate_per_aerosol_when_present() * + self.aerosols()) def emission_rate(self, time) -> _VectorisedFloat: """ @@ -807,10 +822,19 @@ class EmittingPopulation(_PopulationWithVirus): #: The emission rate of a single individual, in virions / h. known_individual_emission_rate: float - @method_cache - def emission_rate_when_present(self) -> _VectorisedFloat: + def aerosols(self): """ - The emission rate if the infected population is present. + total volume of aerosols expired per volume of air (mL/cm^3). + Here arbitrarily set to 1 as the full emission rate is known. + """ + return 1. + + @method_cache + def emission_rate_per_aerosol_when_present(self) -> _VectorisedFloat: + """ + The emission rate of virions per fraction of aerosol volume in + the expired air, if the infected population is present, in cm^3/(mL.h). + It should not be a function of time. """ return self.known_individual_emission_rate * self.number @@ -824,26 +848,27 @@ class InfectedPopulation(_PopulationWithVirus): def fraction_of_infectious_virus(self) -> _VectorisedFloat: """ The fraction of infectious virus. - """ return self.virus.viable_to_RNA_ratio * (1 - self.host_immunity) - @method_cache - def emission_rate_when_present(self) -> _VectorisedFloat: + def aerosols(self): """ - The emission rate if the infected population is present. - Note that the rate is not currently time-dependent. + total volume of aerosols expired per volume of air (mL/cm^3). """ - # Emission Rate (virions / 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 + return self.expiration.aerosols(self.mask) - aerosols = self.expiration.aerosols(self.mask) + @method_cache + def emission_rate_per_aerosol_when_present(self) -> _VectorisedFloat: + """ + The emission rate of virions per fraction of aerosol volume in + the expired air, if the infected population is present, in cm^3/(mL.h). + It should not be a function of time. + """ + # Note on units: exhalation rate is in m^3/h -> 1e6 conversion factor ER = (self.virus.viral_load_in_sputum * self.activity.exhalation_rate * - 10 ** 6 * - aerosols) + 10 ** 6) return ER * self.number @@ -1095,57 +1120,54 @@ class ExposureModel: def exposure(self) -> _VectorisedFloat: """ - The number of virus per meter^3. With sampled diameters, the - aerosol volume has to be put back in the exposure before taking - the mean, to obtain the proper result for the exposure (which - corresponds to an integration on diameters). + The number of virus per meter^3. """ - emission_rate = self.concentration_model.infected.emission_rate_when_present() - if (not isinstance(self.concentration_model.infected,InfectedPopulation) - or not isinstance(self.concentration_model.infected.expiration,Expiration) - or np.isscalar(self.concentration_model.infected.expiration.diameter) - ): - # in all these cases, there is no distribution of - # diameters that need to be integrated over - return self._normed_exposure() * emission_rate + emission_rate_per_aerosol = self.concentration_model.infected.emission_rate_per_aerosol_when_present() + aerosols = self.concentration_model.infected.aerosols() + + diameter = self.concentration_model.infected.particle.diameter + + if not np.isscalar(diameter) and not diameter is None: + # we compute first the mean of all diameter-dependent quantities + # to perform properly the Monte-Carlo integration over + # particle diameters (doing things in another order would + # lead to wrong results). + exposure_integrated = np.array(self._normed_exposure()*aerosols).mean() else: - # the mean of the diameter-dependent exposure (including - # aerosols volume, but NOT the other factors) has to be - # taken first (this is equivalent to integrating over the - # diameters) - mask = self.concentration_model.infected.mask - aerosols = self.concentration_model.infected.expiration.aerosols(mask) - return (np.array(self._normed_exposure()*aerosols).mean() * - emission_rate/aerosols) + # in the case of a single diameter or no diameter defined, + # one should not take any mean at this stage. + exposure_integrated = self._normed_exposure()*aerosols + + # then we multiply by the diameter-independent quantity + # emission_rate_per_aerosol + return exposure_integrated * emission_rate_per_aerosol def _deposited_exposure(self) -> _VectorisedFloat: """ - The number of virus per meter^3 deposited on the respiratory - tract. As in the exposure method, with sampled diameters the - aerosol volume and the deposited fraction, have to be put - in the deposited exposure before taking the mean, to obtain the - proper result (which corresponds to an integration on diameters). + The number of virus per m^3 deposited on the respiratory tract. """ - emission_rate = self.concentration_model.infected.emission_rate_when_present() - if (not isinstance(self.concentration_model.infected,InfectedPopulation) - or not isinstance(self.concentration_model.infected.expiration,Expiration) - or np.isscalar(self.concentration_model.infected.expiration.diameter) - ): - # in all these cases, there is no distribution of - # diameters that need to be integrated over - return (self._normed_exposure() * - self.fraction_deposited() * - emission_rate) + emission_rate_per_aerosol = self.concentration_model.infected.emission_rate_per_aerosol_when_present() + aerosols = self.concentration_model.infected.aerosols() + fdep = self.concentration_model.infected.particle.fraction_deposited() + + diameter = self.concentration_model.infected.particle.diameter + + if not np.isscalar(diameter) and not diameter is None: + # we compute first the mean of all diameter-dependent quantities + # to perform properly the Monte-Carlo integration over + # particle diameters (doing things in another order would + # lead to wrong results). + dep_exposure_integrated = np.array(self._normed_exposure() * + aerosols * + fdep).mean() else: - # the mean of the diameter-dependent exposure (including - # aerosols volume, but NOT the other factors) has to be - # taken first (this is equivalent to integrating over the - # diameters) - mask = self.concentration_model.infected.mask - aerosols = self.concentration_model.infected.expiration.aerosols(mask) - return (np.array(self._normed_exposure() * aerosols * - self.fraction_deposited()).mean() * - emission_rate/aerosols) + # in the case of a single diameter or no diameter defined, + # one should not take any mean at this stage. + dep_exposure_integrated = self._normed_exposure()*aerosols*fdep + + # then we multiply by the diameter-independent quantity + # emission_rate_per_aerosol + return dep_exposure_integrated * emission_rate_per_aerosol def infection_probability(self) -> _VectorisedFloat: deposited_exposure = self._deposited_exposure() From 46146b2a451176daf05c12480c4dfa887f645293 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 Jan 2022 10:14:43 +0100 Subject: [PATCH 13/29] Improving fraction_deposited docstrings (in cara/models.py) --- cara/models.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/cara/models.py b/cara/models.py index 35c85dd3..424f8797 100644 --- a/cara/models.py +++ b/cara/models.py @@ -590,7 +590,9 @@ class Particle: def fraction_deposited(self, evaporation_factor: float=0.3): """ - The fraction of particles actually deposited in the respiratory tract. + The fraction of particles actually deposited in the respiratory + tract (over the total number of particles). It depends on the + particle diameter. From W. C. Hinds, New York, Wiley, 1999 (pp. 233 – 259). evaporation_factor represents the factor applied to the diameter, due to instantaneous evaporation of the particle in the air. @@ -1077,7 +1079,9 @@ class ExposureModel: def fraction_deposited(self): """ - The fraction of viruses actually deposited in the respiratory tract. + The fraction of particles actually deposited in the respiratory + tract (over the total number of particles). It depends on the + particle diameter. """ return self.concentration_model.infected.particle.fraction_deposited( self.concentration_model.evaporation_factor) From 8b485e80a7acf240d49975a91c54ecdb033c730c Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 Jan 2022 10:39:15 +0100 Subject: [PATCH 14/29] Bug fix in Exposure class (cara/models.py) --- cara/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/models.py b/cara/models.py index 424f8797..3dc7b226 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1152,7 +1152,7 @@ class ExposureModel: """ emission_rate_per_aerosol = self.concentration_model.infected.emission_rate_per_aerosol_when_present() aerosols = self.concentration_model.infected.aerosols() - fdep = self.concentration_model.infected.particle.fraction_deposited() + fdep = self.fraction_deposited() diameter = self.concentration_model.infected.particle.diameter From 7af1abcde62f470aed1ea64b10c5a14be0d820c2 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Thu, 27 Jan 2022 15:38:25 +0100 Subject: [PATCH 15/29] Minor fixes in cara/models.py --- cara/models.py | 18 +++++++++--------- cara/tests/test_expiration.py | 4 ++-- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/cara/models.py b/cara/models.py index 3dc7b226..1d60ccb7 100644 --- a/cara/models.py +++ b/cara/models.py @@ -151,9 +151,9 @@ class PiecewiseConstant: def __post_init__(self): if len(self.transition_times) != len(self.values)+1: - raise ValueError("transition_times should contain one more element than values") + raise ValueError("transition_times must contain one more element than values") if tuple(sorted(set(self.transition_times))) != self.transition_times: - raise ValueError("transition_times should not contain duplicated elements and should be sorted") + raise ValueError("transition_times must not contain duplicated elements and must be sorted") shapes = [np.array(v).shape for v in self.values] if not all(shapes[0] == shape for shape in shapes): raise ValueError("All values must have the same shape") @@ -588,7 +588,7 @@ class Particle: vg = 1.88e-4 * (self.diameter*evaporation_factor / 2.5)**2 return vg - def fraction_deposited(self, evaporation_factor: float=0.3): + def fraction_deposited(self, evaporation_factor: float=0.3) -> _VectorisedFloat: """ The fraction of particles actually deposited in the respiratory tract (over the total number of particles). It depends on the @@ -683,10 +683,10 @@ class MultipleExpiration(_ExpirationBase): def __post_init__(self): if len(self.expirations) != len(self.weights): - raise ValueError("expirations and weigths should contain the" + raise ValueError("expirations and weigths must contain the" "same number of elements") if not all(np.isscalar(e.diameter) for e in self.expirations): - raise ValueError("diameters should all be scalars") + raise ValueError("diameters must all be scalars") def aerosols(self, mask: Mask): return np.array([ @@ -844,7 +844,7 @@ class EmittingPopulation(_PopulationWithVirus): @dataclass(frozen=True) class InfectedPopulation(_PopulationWithVirus): #: The type of expiration that is being emitted whilst doing the activity. - expiration: _ExpirationBase + expiration: _ExpirationBase @method_cache def fraction_of_infectious_virus(self) -> _VectorisedFloat: @@ -1077,7 +1077,7 @@ class ExposureModel: #: The number of times the exposure event is repeated (default 1). repeats: int = 1 - def fraction_deposited(self): + def fraction_deposited(self) -> _VectorisedFloat: """ The fraction of particles actually deposited in the respiratory tract (over the total number of particles). It depends on the @@ -1131,7 +1131,7 @@ class ExposureModel: diameter = self.concentration_model.infected.particle.diameter - if not np.isscalar(diameter) and not diameter is None: + if not np.isscalar(diameter) and diameter is not None: # we compute first the mean of all diameter-dependent quantities # to perform properly the Monte-Carlo integration over # particle diameters (doing things in another order would @@ -1156,7 +1156,7 @@ class ExposureModel: diameter = self.concentration_model.infected.particle.diameter - if not np.isscalar(diameter) and not diameter is None: + if not np.isscalar(diameter) and diameter is not None: # we compute first the mean of all diameter-dependent quantities # to perform properly the Monte-Carlo integration over # particle diameters (doing things in another order would diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py index 93c5d4f4..cc82e694 100644 --- a/cara/tests/test_expiration.py +++ b/cara/tests/test_expiration.py @@ -13,7 +13,7 @@ def test_multiple_wrong_weight_size(): e_base = models.Expiration(2.5) with pytest.raises( ValueError, - match=re.escape("expirations and weigths should contain the" + match=re.escape("expirations and weigths must contain the" "same number of elements") ): e = models.MultipleExpiration([e_base, e_base], weights) @@ -26,7 +26,7 @@ def test_multiple_wrong_diameters(): e3 = models.Expiration(2.) with pytest.raises( ValueError, - match=re.escape("diameters should all be scalars") + match=re.escape("diameters must all be scalars") ): e = models.MultipleExpiration([e1, e2, e3], weights) From 68d495ee8e8aeac26a7e523a984002f792719bea Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 28 Jan 2022 11:18:47 +0100 Subject: [PATCH 16/29] Move toronto temperatures into test_monte_carlo_full_models where it is used --- cara/data/__init__.py | 20 ---------------- cara/tests/test_monte_carlo_full_models.py | 27 +++++++++++++++++++++- 2 files changed, 26 insertions(+), 21 deletions(-) diff --git a/cara/data/__init__.py b/cara/data/__init__.py index 6f1ddb0b..31195ba7 100644 --- a/cara/data/__init__.py +++ b/cara/data/__init__.py @@ -22,10 +22,6 @@ geneva_coordinates = (46.204391, 6.143158) local_hourly_temperatures_celsius_per_hour = get_hourly_temperatures_celsius_per_hour( geneva_coordinates) -# Load the weather data (temperature in kelvin) for Toronto. -toronto_coordinates = (43.667, 79.400) -toronto_hourly_temperatures_celsius_per_hour = get_hourly_temperatures_celsius_per_hour( - toronto_coordinates) # Geneva hourly temperatures as piecewise constant function (in Kelvin). GenevaTemperatures_hourly = { @@ -38,25 +34,9 @@ GenevaTemperatures_hourly = { for month, temperatures in local_hourly_temperatures_celsius_per_hour.items() } -# Toronto hourly temperatures as piecewise constant function (in Kelvin). -TorontoTemperatures_hourly = { - month: models.PiecewiseConstant( - # NOTE: It is important that the time type is float, not np.float, in - # order to allow hashability (for caching). - tuple(float(time) for time in range(25)), - tuple(273.15 + np.array(temperatures)), - ) - for month, temperatures in toronto_hourly_temperatures_celsius_per_hour.items() -} # Same Geneva temperatures on a finer temperature mesh (every 6 minutes). GenevaTemperatures = { month: GenevaTemperatures_hourly[month].refine(refine_factor=10) for month, temperatures in local_hourly_temperatures_celsius_per_hour.items() } - -# Same Toronto temperatures on a finer temperature mesh (every 6 minutes). -TorontoTemperatures = { - month: TorontoTemperatures_hourly[month].refine(refine_factor=10) - for month, temperatures in toronto_hourly_temperatures_celsius_per_hour.items() -} diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index e5c96b82..407ea37d 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -12,6 +12,31 @@ np.random.seed(2000) SAMPLE_SIZE = 250000 TOLERANCE = 0.05 +# Load the weather data (temperature in kelvin) for Toronto. +toronto_coordinates = (43.667, 79.400) +toronto_hourly_temperatures_celsius_per_hour = data.get_hourly_temperatures_celsius_per_hour( + toronto_coordinates) + + +# Toronto hourly temperatures as piecewise constant function (in Kelvin). +TorontoTemperatures_hourly = { + month: models.PiecewiseConstant( + # NOTE: It is important that the time type is float, not np.float, in + # order to allow hashability (for caching). + tuple(float(time) for time in range(25)), + tuple(273.15 + np.array(temperatures)), + ) + for month, temperatures in toronto_hourly_temperatures_celsius_per_hour.items() +} + + +# Same Toronto temperatures on a finer temperature mesh (every 6 minutes). +TorontoTemperatures = { + month: TorontoTemperatures_hourly[month].refine(refine_factor=10) + for month, temperatures in toronto_hourly_temperatures_celsius_per_hour.items() +} + + # references values for infection_probability and expected new cases # in the following tests, were obtained from the feature/mc branch @@ -69,7 +94,7 @@ def classroom_mc(): models.SlidingWindow( active=models.PeriodicInterval(period=120, duration=120), inside_temp=models.PiecewiseConstant((0., 24.), (293,)), - outside_temp=data.TorontoTemperatures['Dec'], + outside_temp=TorontoTemperatures['Dec'], window_height=1.6, opening_length=0.2, ), From 3e4aa09eee0d9b37b8713ee3b9eecf1583f1103f Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 1 Feb 2022 09:30:41 +0100 Subject: [PATCH 17/29] Moved SAMPLE_SIZE and TOLERANCE to the respective test where they are used --- cara/tests/apps/calculator/test_model_generator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cara/tests/apps/calculator/test_model_generator.py b/cara/tests/apps/calculator/test_model_generator.py index b9c5a852..5cc9e409 100644 --- a/cara/tests/apps/calculator/test_model_generator.py +++ b/cara/tests/apps/calculator/test_model_generator.py @@ -13,8 +13,6 @@ from cara.monte_carlo.data import expiration_distributions # TODO: seed better the random number generators np.random.seed(2000) -SAMPLE_SIZE = 250000 -TOLERANCE = 0.02 def test_model_from_dict(baseline_form_data): form = model_generator.FormData.from_dict(baseline_form_data) @@ -35,6 +33,8 @@ def test_model_from_dict_invalid(baseline_form_data): ] ) def test_blend_expiration(mask_type): + SAMPLE_SIZE = 250000 + TOLERANCE = 0.02 blend = {'Breathing': 2, 'Speaking': 1} r = model_generator.build_expiration(blend).build_model(SAMPLE_SIZE) mask = models.Mask.types[mask_type] From 02468f42fabb7bb906cbff0b148dc26fbe6289d3 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 1 Feb 2022 12:14:49 +0100 Subject: [PATCH 18/29] updated threshold values --- .../templates/cern/calculator.report.html.j2 | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/cara/apps/templates/cern/calculator.report.html.j2 b/cara/apps/templates/cern/calculator.report.html.j2 index 78e73b8f..07c03361 100644 --- a/cara/apps/templates/cern/calculator.report.html.j2 +++ b/cara/apps/templates/cern/calculator.report.html.j2 @@ -7,9 +7,9 @@ {% endblock report_preamble_navtab %} {% block warning_animation %} - {% if ((prob_inf > 15) or (expected_new_cases >= 1)) %} {% set warning_color= 'bg-danger' %} - {% elif (5 <= prob_inf <= 15) %} {% set warning_color = 'bg-warning' %} - {% elif (prob_inf < 5) %} {% set warning_color = 'bg-success' %} + {% if ((prob_inf > 10) or (expected_new_cases >= 1)) %} {% set warning_color= 'bg-danger' %} + {% elif (2 <= prob_inf <= 10) %} {% set warning_color = 'bg-warning' %} + {% elif (prob_inf < 2) %} {% set warning_color = 'bg-success' %} {% endif %}
@@ -23,25 +23,25 @@ {% block report_summary %}
- {% if ((prob_inf > 15) or (expected_new_cases >= 1)) %} + {% if ((prob_inf > 10) or (expected_new_cases >= 1)) %} - {% elif 5 <= prob_inf <= 15 %} + {% elif 2 <= prob_inf <= 10 %} - {% elif prob_inf < 5 %} + {% elif prob_inf < 2 %} {% endif %} - {% if (prob_inf > 5) %} + {% if (prob_inf > 2) %}
{% if scale_warning.level == "green-1" %}