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 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:
-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: - -
| # of repeated events | -P(i) | -Expected new cases | -
|---|---|---|
| {{ repeat_event.repeats }} | -{{ repeat_event.probability_of_infection | non_zero_percentage }} | -{{ repeat_event.expected_new_cases | float_format }} | -
Alternative scenarios:
@@ -248,7 +225,7 @@
@@ -269,7 +246,7 @@
Scenario
- P(i)
+ P(I)
Expected new cases
Notes for alternative scenarios:
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
@@ -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.
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.
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.
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.
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:
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 %} 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):
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. 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 c2500d68..ad4493ba 100644 --- a/cara/models.py +++ b/cara/models.py @@ -37,6 +37,7 @@ import typing import numpy as np from scipy.interpolate import interp1d +import scipy.integrate if not typing.TYPE_CHECKING: from memoization import cached @@ -473,38 +474,47 @@ Virus.types = { @dataclass(frozen=True) class Mask: - #: Filtration efficiency. (In %/100) - η_exhale: _VectorisedFloat - - #: Leakage through side of masks. - η_leaks: _VectorisedFloat - #: Filtration efficiency of masks when inhaling. η_inhale: _VectorisedFloat + #: Global factor applied to filtration efficiency of masks when exhaling. + factor_exhale: float = 1. + #: 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. - if diameter < 3e-4: + """ + 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 microns. + """ + if diameter < 0.5: eta_out = 0. + elif diameter < 0.94614: + eta_out = 0.5893 * diameter + 0.1546 + elif diameter < 3.: + eta_out = 0.0509 * diameter + 0.664 else: - eta_out = self.η_exhale * (1 - self.η_leaks) - return eta_out + eta_out = 0.8167 + return eta_out*self.factor_exhale + + def inhale_efficiency(self) -> _VectorisedFloat: + """ + Overall inhale efficiency, including the effect of the leaks. + """ + return self.η_inhale Mask.types = { - 'No mask': Mask(0, 0, 0), + 'No mask': Mask(0, 0), 'Type I': Mask( - η_exhale=0.95, - η_leaks=0.15, # (Huang 2007) - η_inhale=0.3, # (Browen 2010) + η_inhale=0.5, # (CERN-OPEN-2021-004) ), '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) ), } @@ -516,35 +526,59 @@ 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: Mask): - # total volume of aerosols expired per volume of air (mL/cm^3). + """ + total volume of aerosols expired per volume of air (mL/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. 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. + 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). + Here all diameters (d) are in microns. """ - ejection_factor: typing.Tuple[float, ...] - particle_sizes: typing.Tuple[float, ...] = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4) # In cm. + #: 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] + @cached() 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))) - total += contribution - return total + """ Result is in mL.cm^-3 """ + def volume(d): + return (np.pi * d**3) / 6. + + def _Bmode(d: float) -> float: + # B-mode (see ref. above). + 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.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.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) * (1 - mask.exhale_efficiency(d)) + + # 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) @@ -572,17 +606,20 @@ class MultipleExpiration(_ExpirationBase): _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': Expiration((1., 0., 0.)), + 'Talking': Expiration((1., 1., 1.)), + 'Shouting': Expiration((1., 5., 5.)), + 'Singing': Expiration((1., 5., 5.)), + 'Superspreading event': Expiration((np.inf, 0., 0.)), } @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. @@ -637,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 * @@ -672,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. @@ -692,8 +730,8 @@ class ConcentrationModel: return self.infected.virus def infectious_virus_removal_rate(self, time: float) -> _VectorisedFloat: - # Particle deposition on the floor - vg = 1 * 10 ** -4 + # 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 # Deposition rate (h^-1) @@ -702,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 @@ -714,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 @@ -759,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. @@ -777,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) @@ -821,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 @@ -835,8 +877,8 @@ class ExposureModel: inf_aero = ( self.exposed.activity.inhalation_rate * - (1 - self.exposed.mask.η_inhale) * - exposure + (1 - self.exposed.mask.inhale_efficiency()) * + exposure * self.fraction_deposited ) # Probability of infection. diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py new file mode 100644 index 00000000..79e25517 --- /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 +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=symptomatic_vl_frequencies, + quantum_infectious_dose=100, + ), + 'SARS_CoV_2_B117': mc.SARSCoV2( + viral_load_in_sputum=symptomatic_vl_frequencies, + quantum_infectious_dose=60, + ), + 'SARS_CoV_2_P1': mc.SARSCoV2( + viral_load_in_sputum=symptomatic_vl_frequencies, + quantum_infectious_dose=100/2.25, + ), +} diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 4ed49d82..4333c93f 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,96 @@ 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] + + +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, ] 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]) 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..0a22fe59 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 = 20. @pytest.fixture def app(): @@ -45,11 +46,12 @@ class TestBasicApp(tornado.testing.AsyncHTTPTestCase): def get_app(self): return cara.apps.calculator.make_app() + @tornado.testing.gen_test(timeout=_TIMEOUT) 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 +59,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=_TIMEOUT) 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/conftest.py b/cara/tests/conftest.py index 24983cf1..023e4e5a 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,9 +28,10 @@ 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, - ) + ), + fraction_deposited = 1., ) diff --git a/cara/tests/models/test_concentration_model.py b/cara/tests/models/test_concentration_model.py index 48ef0b66..26707236 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,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])}, - {'η_exhale': np.array([0.92, 0.95])}, - {'η_leaks': np.array([0.15, 0.20])}, ] ) def test_concentration_model_vectorisation(override_params): @@ -24,8 +23,6 @@ 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, } defaults.update(override_params) @@ -37,8 +34,7 @@ def test_concentration_model_vectorisation(override_params): number=1, presence=always, mask=models.Mask( - η_exhale=defaults['η_exhale'], - η_leaks=defaults['η_leaks'], + factor_exhale=0.95, η_inhale=0.3, ), activity=models.Activity( @@ -49,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) @@ -128,4 +122,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/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 464f25a9..1c6c6d3d 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. @@ -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), - np.array([14.4, 14.4]), np.array([99.4146994564, 99.6803184113])], + [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])), - np.array([14.4, 28.8]), np.array([99.6803184113, 99.9989780368])], + [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 ) @@ -123,22 +126,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.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, @@ -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/models/test_mask.py b/cara/tests/models/test_mask.py new file mode 100644 index 00000000..4a71d13c --- /dev/null +++ b/cara/tests/models/test_mask.py @@ -0,0 +1,33 @@ +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_mask_inhale(η_inhale, expected_inhale_efficiency): + mask = models.Mask(η_inhale=η_inhale) + npt.assert_equal(mask.inhale_efficiency(), + expected_inhale_efficiency) + + +@pytest.mark.parametrize( + "diameter, factor_exhale, expected_exhale_efficiency", + [ + [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): + mask = models.Mask(η_inhale=0.3, factor_exhale=factor_exhale) + npt.assert_almost_equal(mask.exhale_efficiency(diameter), + expected_exhale_efficiency) diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py index 7638d613..434caf8c 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,11 +19,25 @@ 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)) + weights = (1., 1.) + 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) - 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['Type I'] + npt.assert_almost_equal(e_expected.aerosols(mask), e.aerosols(mask)) + + +# expected values obtained from another code +@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) diff --git a/cara/tests/test_infected_population.py b/cara/tests/test_infected_population.py index e87ddb20..7acef7f7 100644 --- a/cara/tests/test_infected_population.py +++ b/cara/tests/test_infected_population.py @@ -8,8 +8,6 @@ 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])}, {'exhalation_rate': np.array([0.75, 0.81])}, ] ) @@ -17,8 +15,6 @@ def test_infected_population_vectorisation(override_params): defaults = { 'viral_load_in_sputum': 1e9, 'quantum_infectious_dose': 50, - 'η_exhale': 0.95, - 'η_leaks': 0.15, 'exhalation_rate': 0.75, } defaults.update(override_params) @@ -28,8 +24,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=0.95, η_inhale=0.3, ), activity=cara.models.Activity( @@ -40,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) diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index ea8a1ad0..d8597cf8 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.41611256, 1.3205628e-14, 0.41611256, 4.1909001e-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, 972.880852) 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 @@ -366,20 +349,21 @@ def build_exposure_model(concentration_model): presence=infected.presence, activity=infected.activity, mask=infected.mask, - ) + ), + fraction_deposited = 1., ) -# 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', 9.930854], + ['Jun', 37.962708], ], ) -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 +371,20 @@ 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 -# a mesh of 100'000 pts per exposed presence interval. +# expected quanta 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_r0", + "month, expected_quanta", [ - ['Jan', 86.422736], - ['Jun', 99.982323], + ['Jan', 9.993842], + ['Jun', 40.151985], ], ) -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 +393,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, rtol=0.02) diff --git a/cara/tests/test_monte_carlo.py b/cara/tests/test_monte_carlo.py index f3d9fb08..23e139fd 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 @@ -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): 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..7c08757f --- /dev/null +++ b/cara/tests/test_monte_carlo_full_models.py @@ -0,0 +1,316 @@ +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 = 50000 +TOLERANCE = 0.05 + +# references values for infection_probability and expected new cases +# in the following tests, were obtained from the feature/mc branch + +@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( + ( + 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)), + ), + ) + return 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, + ), + ) + + +@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( + ( + 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'], + ), + ) + return 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, + ), + ) + + +@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( + 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'], + ), + ) + return 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, + ), + ) + + +@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( + 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'], + ), + ) + return 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, + ), + ) + + +@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( + 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)), + ), + ) + return 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, + ), + ) + + +@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( + 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.)), + ), + ) + return 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, + ), + ) + + +@pytest.mark.parametrize( + "mc_model, expected_pi, expected_new_cases, expected_dose, expected_qR", + [ + ["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_pi, + expected_dose, expected_qR): + 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_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) diff --git a/cara/tests/test_predefined_distributions.py b/cara/tests/test_predefined_distributions.py new file mode 100644 index 00000000..4197b2fb --- /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(2000) + + +# 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 new file mode 100644 index 00000000..c30fcbf7 --- /dev/null +++ b/cara/tests/test_sampleable_distribution.py @@ -0,0 +1,112 @@ +import numpy as np +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], + ] +) +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) + 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], rtol=0.01) + npt.assert_allclose(selected_histogram, exact_dist, rtol=0.02) + + +@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) + 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], rtol=0.01) + npt.assert_allclose(selected_histogram, exact_dist, rtol=0.03) + + +@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) + else: + samples = sampleable.Custom((0, 10), function, max_function + ).generate_samples(sample_size) + + histogram, bins = np.histogram(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) + + +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) diff --git a/requirements.txt b/requirements.txt index 1bbecfe1..eac04f3e 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 +scikit_learn==0.23.1 Send2Trash==1.5.0 six==1.15.0 sniffio==1.2.0 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', ],