diff --git a/cara/apps/calculator/__init__.py b/cara/apps/calculator/__init__.py index 6bfcb9d7..b9572896 100644 --- a/cara/apps/calculator/__init__.py +++ b/cara/apps/calculator/__init__.py @@ -33,7 +33,7 @@ from .user import AuthenticatedUser, AnonymousUser # calculator version. If the calculator needs to make breaking changes (e.g. change # form attributes) then it can also increase its MAJOR version without needing to # increase the overall CARA version (found at ``cara.__version__``). -__version__ = "3.3.0" +__version__ = "4.0.0" class BaseRequestHandler(RequestHandler): diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 917cacf5..c6205519 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -138,7 +138,7 @@ def calculate_report_data(model: models.ExposureModel): exposed_occupants = model.exposed.number expected_new_cases = np.array(model.expected_new_cases()).mean() cumulative_doses = np.cumsum([ - np.array(model.exposure_between_bounds(float(time1), float(time2))).mean() + np.array(model.deposited_exposure_between_bounds(float(time1), float(time2))).mean() for time1, time2 in zip(times[:-1], times[1:]) ]) @@ -337,11 +337,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', - 'incidence_rate': 'higher or equal to 100 new cases per 100 000 inhabitants', - 'onsite_access': 'lower than 4000', - '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': '2%' # '2%' - '' - '' - '' } return context diff --git a/cara/apps/templates/cern/calculator.report.html.j2 b/cara/apps/templates/cern/calculator.report.html.j2 index 78e73b8f..cea08ada 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,41 +23,41 @@ {% 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" %} {% 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.

@@ -67,13 +67,13 @@ {% endblock report_summary %} {% block report_summary_footnote %} - {% if ((prob_inf > 15) or (expected_new_cases >= 1)) %} + {% if ((prob_inf > 10) or (expected_new_cases >= 1)) %} This exceeds the authorised risk threshold or number of expected new cases. The risk level must be reduced before this activity can be undertaken. - {% elif (5 <= prob_inf <= 15) %} + {% elif (2 <= prob_inf <= 10) %} This activity has an elevated level of risk, ALARA principles must be applied to minimise the level of risk before undertaking the activity. See the footnotes for more details on the ALARA principles. - {% elif (prob_inf < 5) %} + {% elif (prob_inf < 2) %} This level of risk is within acceptable parameters, no further actions are required. {% endif %} {% endblock report_summary_footnote %} @@ -124,9 +124,9 @@
- - - + + +
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/models.py b/cara/models.py index 30aeaafd..6b5e6540 100644 --- a/cara/models.py +++ b/cara/models.py @@ -51,7 +51,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. @@ -152,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") @@ -460,8 +459,8 @@ class SARSCoV2(Virus): piecewise constant model (for more details see A. Henriques et al, CERN-OPEN-2021-004, DOI: 10.17181/CERN.1GDQ.5Y75) """ - # Taken from Morris et al (https://doi.org/10.7554/eLife.65902) data at T = 22°C and RH = 40 % and - # from Doremalen et al (https://www.nejm.org/doi/10.1056/NEJMc2004973). + # Taken from Morris et al (https://doi.org/10.7554/eLife.65902) data at T = 22°C and RH = 40 %, + # and from Doremalen et al (https://www.nejm.org/doi/10.1056/NEJMc2004973). return np.piecewise(humidity, [humidity <= 0.4, humidity > 0.4], [6.43, 1.1]) @@ -561,6 +560,58 @@ 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) -> _VectorisedFloat: + """ + 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. + """ + 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: """ @@ -570,6 +621,13 @@ 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 + """ + return Particle() + def aerosols(self, mask: Mask): """ total volume of aerosols expired per volume of air (mL/cm^3). @@ -587,9 +645,18 @@ 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. + @property + 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 """ @@ -608,18 +675,18 @@ 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, ...] 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([ @@ -709,13 +776,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: """ @@ -734,16 +816,33 @@ class _PopulationWithVirus(Population): return self.emission_rate_when_present() + @property + 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): #: 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 @@ -757,29 +856,37 @@ 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 + @property + def particle(self) -> Particle: + """ + the Particle object representing the aerosol - here the default one + """ + return self.expiration.particle + @dataclass(frozen=True) class ConcentrationModel: @@ -798,20 +905,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) @@ -989,26 +1083,14 @@ 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 viruses actually deposited in the respiratory tract. - From W. C. Hinds, New York, Wiley, 1999 (pp. 233 – 259). + The fraction of particles actually deposited in the respiratory + tract (over the total number of particles). It depends on the + particle diameter. """ - 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 @@ -1028,11 +1110,6 @@ class ExposureModel: elif time1 <= start and stop < time2: exposure += self.concentration_model.normed_integrated_concentration(start, stop) return exposure - - def exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: - """The number of virions per meter^3 between any two times.""" - return (self._normed_exposure_between_bounds(time1, time2) * - self.concentration_model.infected.emission_rate_when_present()) def _normed_exposure(self) -> _VectorisedFloat: """ @@ -1046,76 +1123,57 @@ class ExposureModel: return normed_exposure * self.repeats - def exposure(self) -> _VectorisedFloat: + def deposited_exposure_between_bounds(self, time1: float, time2: float) -> _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 m^3 deposited on the respiratory tract + between any two times. """ - 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 - else: - # the mean of the diameter-dependent exposure (including - # aerosols volume, but NOT the other factors) has to be - # taken first (this is equivalent to integrating over the - # diameters) - mask = self.concentration_model.infected.mask - aerosols = self.concentration_model.infected.expiration.aerosols(mask) - return (np.array(self._normed_exposure()*aerosols).mean() * - emission_rate/aerosols) - - def _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). - """ - 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) - 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) - - def infection_probability(self) -> _VectorisedFloat: - deposited_exposure = self._deposited_exposure() - + emission_rate_per_aerosol = self.concentration_model.infected.emission_rate_per_aerosol_when_present() + aerosols = self.concentration_model.infected.aerosols() + fdep = self.fraction_deposited() f_inf = self.concentration_model.infected.fraction_of_infectious_virus() - inf_aero = ( - self.exposed.activity.inhalation_rate * - (1 - self.exposed.mask.inhale_efficiency()) * - deposited_exposure * f_inf - ) + diameter = self.concentration_model.infected.particle.diameter + + 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 + # lead to wrong results). + dep_exposure_integrated = np.array(self._normed_exposure_between_bounds(time1, time2) * + aerosols * + fdep).mean() + else: + # 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_between_bounds(time1, time2)*aerosols*fdep + + # then we multiply by the diameter-independent quantity emission_rate_per_aerosol, + # and parameters of the vD equation (i.e. f_inf, BR_k and n_in). + return (dep_exposure_integrated * emission_rate_per_aerosol * + f_inf * self.exposed.activity.inhalation_rate * + (1 - self.exposed.mask.inhale_efficiency())) + + def deposited_exposure(self) -> _VectorisedFloat: + """ + The number of virus per m^3 deposited on the respiratory tract. + """ + deposited_exposure = 0.0 + + for start, stop in self.exposed.presence.boundaries(): + deposited_exposure += self.deposited_exposure_between_bounds(start, stop) + + return deposited_exposure * self.repeats + + def infection_probability(self) -> _VectorisedFloat: + # viral dose (vD) + vD = self.deposited_exposure() # oneoverln2 multiplied by ID_50 corresponds to ID_63. infectious_dose = oneoverln2 * self.concentration_model.virus.infectious_dose # Probability of infection. - return (1 - np.exp(-((inf_aero * (1 - self.exposed.host_immunity))/(infectious_dose * + return (1 - np.exp(-((vD * (1 - self.exposed.host_immunity))/(infectious_dose * self.concentration_model.virus.transmissibility_factor)))) * 100 def expected_new_cases(self) -> _VectorisedFloat: diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index c4e1dfb2..9dee5d49 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -165,15 +165,15 @@ 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 distribution between 0.1 and D microns - these boundaries are + 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). """ 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.)) def dilution_factor(distance, D=0.02): 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] diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 5b674207..15d20986 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -77,25 +77,25 @@ def known_concentrations(func): @pytest.mark.parametrize( "population, cm, expected_exposure, expected_probability", [ [populations[1], known_concentrations(lambda t: 36.), - np.array([432, 432]), np.array([67.9503762594, 65.2366759251])], + np.array([64.02320633, 59.45012016]), np.array([67.9503762594, 65.2366759251])], [populations[2], known_concentrations(lambda t: 36.), - np.array([432, 432]), np.array([51.6749232285, 55.6374622042])], + np.array([40.91708675, 45.73086166]), np.array([51.6749232285, 55.6374622042])], [populations[0], known_concentrations(lambda t: np.array([36., 72.])), - np.array([432, 864]), np.array([55.6374622042, 80.3196524031])], + np.array([45.73086166, 91.46172332]), np.array([55.6374622042, 80.3196524031])], [populations[1], known_concentrations(lambda t: np.array([36., 72.])), - np.array([432, 864]), np.array([67.9503762594, 87.9151129926])], + np.array([64.02320633, 118.90024032]), np.array([67.9503762594, 87.9151129926])], [populations[2], known_concentrations(lambda t: np.array([36., 72.])), - np.array([432, 864]), np.array([51.6749232285, 80.3196524031])], + np.array([40.91708675, 91.46172332]), np.array([51.6749232285, 80.3196524031])], ]) def test_exposure_model_ndarray(population, cm, expected_exposure, expected_probability): model = ExposureModel(cm, population) np.testing.assert_almost_equal( - model.exposure(), expected_exposure + model.deposited_exposure(), expected_exposure ) np.testing.assert_almost_equal( model.infection_probability(), expected_probability, decimal=10 @@ -107,33 +107,43 @@ def test_exposure_model_ndarray(population, cm, assert model.expected_new_cases().shape == (2,) -@pytest.mark.parametrize("population", populations) -def test_exposure_model_ndarray_and_float_mix(population): +@pytest.mark.parametrize("population, expected_deposited_exposure", [ + [populations[0], np.array([1.52436206, 1.52436206])], + [populations[1], np.array([2.13410688, 1.98167067])], + [populations[2], np.array([1.36390289, 1.52436206])], + ]) +def test_exposure_model_ndarray_and_float_mix(population, expected_deposited_exposure): cm = known_concentrations( lambda t: 0. if np.floor(t) % 2 else np.array([1.2, 1.2])) model = ExposureModel(cm, population) - expected_exposure = np.array([14.4, 14.4]) np.testing.assert_almost_equal( - model.exposure(), expected_exposure + model.deposited_exposure(), expected_deposited_exposure ) assert isinstance(model.infection_probability(), np.ndarray) assert isinstance(model.expected_new_cases(), np.ndarray) -@pytest.mark.parametrize("population", populations) -def test_exposure_model_compare_scalar_vector(population): - cm_scalar = known_concentrations(lambda t: 1.2) +@pytest.mark.parametrize("population, expected_deposited_exposure", [ + [populations[0], np.array([1.52436206, 1.52436206])], + [populations[1], np.array([2.13410688, 1.98167067])], + [populations[2], np.array([1.36390289, 1.52436206])], + ]) +def test_exposure_model_vector(population, expected_deposited_exposure): cm_array = known_concentrations(lambda t: np.array([1.2, 1.2])) - model_scalar = ExposureModel(cm_scalar, population) model_array = ExposureModel(cm_array, population) - expected_exposure = 14.4 np.testing.assert_almost_equal( - model_scalar.exposure(), expected_exposure + model_array.deposited_exposure(), np.array(expected_deposited_exposure) ) + + +def test_exposure_model_scalar(): + cm_scalar = known_concentrations(lambda t: 1.2) + model_scalar = ExposureModel(cm_scalar, populations[0]) + expected_deposited_exposure = 1.52436206 np.testing.assert_almost_equal( - model_array.exposure(), np.array([expected_exposure]*2) + model_scalar.deposited_exposure(), expected_deposited_exposure ) @@ -163,28 +173,28 @@ def conc_model(): ) -# Expected exposure were computed with a trapezoidal integration, using +# Expected deposited exposure were computed with a trapezoidal integration, using # a mesh of 10'000 pts per exposed presence interval. @pytest.mark.parametrize( - ["exposed_time_interval", "expected_exposure"], + ["exposed_time_interval", "expected_deposited_exposure"], [ - [(0., 1.), 266.67176], - [(1., 1.01), 3.0879539], - [(1.01, 1.02), 3.00082435], - [(12., 12.01), 0.095063235], - [(12., 24.), 3775.65025], - [(0., 24.), 4097.8494], + [(0., 1.), 45.6008710], + [(1., 1.01), 0.5280401], + [(1.01, 1.02), 0.51314096385], + [(12., 12.01), 0.016255813185], + [(12., 24.), 645.63619275], + [(0., 24.), 700.7322474], ] ) def test_exposure_model_integral_accuracy(exposed_time_interval, - expected_exposure, conc_model): + expected_deposited_exposure, conc_model): presence_interval = models.SpecificInterval((exposed_time_interval,)) population = models.Population( 10, presence_interval, models.Mask.types['Type I'], models.Activity.types['Standing'], 0., ) model = ExposureModel(conc_model, population) - np.testing.assert_allclose(model.exposure(), expected_exposure) + np.testing.assert_allclose(model.deposited_exposure(), expected_deposited_exposure) def test_infectious_dose_vectorisation(): @@ -212,7 +222,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, ) 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) diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 9910a848..b94e6acd 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -97,7 +97,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): @@ -389,16 +389,16 @@ def build_exposure_model(concentration_model): ) -# expected exposure were computed with a trapezoidal integration, using +# expected deposited exposure were computed with a trapezoidal integration, using # a mesh of 100'000 pts per exposed presence interval. @pytest.mark.parametrize( - "month, expected_exposure", + "month, expected_deposited_exposure", [ - ['Jan', 503.254087759], - ['Jun', 2294.71115639], + ['Jan', 377.440565819], + ['Jun', 1721.03336729], ], ) -def test_exposure_hourly_dep(month,expected_exposure): +def test_exposure_hourly_dep(month,expected_deposited_exposure): m = build_exposure_model( build_hourly_dependent_model( month, @@ -406,20 +406,20 @@ def test_exposure_hourly_dep(month,expected_exposure): intervals_presence_infected=((8., 12.), (13., 17.)) ) ) - exposure = m.exposure() - npt.assert_allclose(exposure, expected_exposure) + deposited_exposure = m.deposited_exposure() + npt.assert_allclose(deposited_exposure, expected_deposited_exposure) -# expected exposure were computed with a trapezoidal integration, using +# expected deposited exposure were computed with a trapezoidal integration, using # a mesh of 100'000 pts per exposed presence interval and 25 pts per hour # for the temperature discretization. @pytest.mark.parametrize( - "month, expected_exposure", + "month, expected_deposited_exposure", [ - ['Jan', 511.118941481], - ['Jun', 2398.90129579], + ['Jan', 383.339206111], + ['Jun', 1799.17597184], ], ) -def test_exposure_hourly_dep_refined(month,expected_exposure): +def test_exposure_hourly_dep_refined(month,expected_deposited_exposure): m = build_exposure_model( build_hourly_dependent_model( month, @@ -428,5 +428,5 @@ def test_exposure_hourly_dep_refined(month,expected_exposure): temperatures=data.GenevaTemperatures, ) ) - exposure = m.exposure() - npt.assert_allclose(exposure, expected_exposure, rtol=0.02) + deposited_exposure = m.deposited_exposure() + npt.assert_allclose(deposited_exposure, expected_deposited_exposure, rtol=0.02) diff --git a/cara/tests/test_monte_carlo.py b/cara/tests/test_monte_carlo.py index c51bcbeb..339612df 100644 --- a/cara/tests/test_monte_carlo.py +++ b/cara/tests/test_monte_carlo.py @@ -89,6 +89,6 @@ def test_build_concentration_model(baseline_mc_model: cara.monte_carlo.Concentra def test_build_exposure_model(baseline_mc_exposure_model: cara.monte_carlo.ExposureModel): model = baseline_mc_exposure_model.build_model(7) assert isinstance(model, cara.models.ExposureModel) - prob = model.exposure() + prob = model.deposited_exposure() assert isinstance(prob, np.ndarray) assert prob.shape == (7, ) diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 90466d96..4a68b9e1 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 @@ -22,7 +47,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 +61,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'], @@ -51,7 +76,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., @@ -66,12 +91,12 @@ 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), inside_temp=models.PiecewiseConstant((0., 24.), (293,)), - outside_temp=data.TorontoTemperatures['Dec'], + outside_temp=TorontoTemperatures['Dec'], window_height=1.6, opening_length=0.2, ), @@ -295,13 +320,13 @@ def waiting_room_mc(): @pytest.mark.parametrize( "mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER", [ - ["shared_office_mc", 6.03, 0.18, 27.22, 809], - ["classroom_mc", 9.5, 1.85, 79.98, 5624], - ["ski_cabin_mc", 16.0, 0.5, 40.25, 7966], - ["skagit_chorale_mc",65.7, 40.0, 241.28, 190422], - ["bus_ride_mc", 12.0, 8.0, 63.79, 5419], - ["gym_mc", 0.45, 0.13, 0.4852, 1145], - ["waiting_room_mc", 1.59, 0.22, 7.23, 737], + ["shared_office_mc", 6.03, 0.18, 3.198, 809], + ["classroom_mc", 9.5, 1.85, 9.478, 5624], + ["ski_cabin_mc", 16.0, 0.5, 17.315, 7966], + ["skagit_chorale_mc",65.7, 40.0, 102.213, 190422], + ["bus_ride_mc", 12.0, 8.0, 7.65, 5419], + ["gym_mc", 0.45, 0.13, 0.208, 1145], + ["waiting_room_mc", 1.59, 0.22, 0.821, 737], ] ) def test_report_models(mc_model, expected_pi, expected_new_cases, @@ -312,7 +337,7 @@ def test_report_models(mc_model, expected_pi, expected_new_cases, expected_pi, rtol=TOLERANCE) npt.assert_allclose(exposure_model.expected_new_cases().mean(), expected_new_cases, rtol=TOLERANCE) - npt.assert_allclose(exposure_model.exposure().mean(), + npt.assert_allclose(exposure_model.deposited_exposure().mean(), expected_dose, rtol=TOLERANCE) npt.assert_allclose( exposure_model.concentration_model.infected.emission_rate_when_present().mean(), @@ -322,10 +347,10 @@ def test_report_models(mc_model, expected_pi, expected_new_cases, @pytest.mark.parametrize( "mask_type, month, expected_pi, expected_dose, expected_ER", [ - ["No mask", "Jul", 9.52, 84.54, 809], - ["Type I", "Jul", 1.7, 15.64, 149], - ["FFP2", "Jul", 0.51, 15.64, 149], - ["Type I", "Feb", 0.57, 4.59, 162], + ["No mask", "Jul", 9.52, 9.920, 809], + ["Type I", "Jul", 1.7, 0.913, 149], + ["FFP2", "Jul", 0.51, 0.239, 149], + ["Type I", "Feb", 0.57, 0.272, 162], ], ) def test_small_shared_office_Geneva(mask_type, month, expected_pi, @@ -372,7 +397,7 @@ def test_small_shared_office_Geneva(mask_type, month, expected_pi, exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) npt.assert_allclose(exposure_model.infection_probability().mean(), expected_pi, rtol=TOLERANCE) - npt.assert_allclose(exposure_model.exposure().mean(), + npt.assert_allclose(exposure_model.deposited_exposure().mean(), expected_dose, rtol=TOLERANCE) npt.assert_allclose( exposure_model.concentration_model.infected.emission_rate_when_present().mean(),