diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index ad65795e..1d4f51ec 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -8,7 +8,9 @@ import numpy as np from cara import models from cara import data +import cara.monte_carlo as mc from .. import calculator +from cara.monte_carlo.data import activity_distributions, virus_distributions LOG = logging.getLogger(__name__) @@ -20,7 +22,7 @@ minutes_since_midnight = typing.NewType('minutes_since_midnight', int) # Used to declare when an attribute of a class must have a value provided, and # there should be no default value used. _NO_DEFAULT = object() - +_SAMPLE_SIZE = 50000 @dataclass class FormData: @@ -263,11 +265,15 @@ class FormData: ventilation = models.HVACMechanical( active=always_on, q_air_mech=self.air_supply) + # this is a minimal, always present source of ventilation, due + # to the air infiltration from the outside. + # See CERN-OPEN-2021-004, p. 12. + infiltration_ventilation = models.AirChange(active=always_on, air_exch=0.25) if self.hepa_option: hepa = models.HEPAFilter(active=always_on, q_air_mech=self.hepa_amount) - return models.MultipleVentilation((ventilation, hepa)) + return models.MultipleVentilation((ventilation, hepa, infiltration_ventilation)) else: - return ventilation + return models.MultipleVentilation((ventilation, infiltration_ventilation)) def mask(self) -> models.Mask: # Initializes the mask type if mask wearing is "continuous", otherwise instantiates the mask attribute as @@ -275,9 +281,9 @@ class FormData: mask = models.Mask.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] return mask - def infected_population(self) -> models.InfectedPopulation: + def infected_population(self) -> mc.InfectedPopulation: # Initializes the virus - virus = models.Virus.types[self.virus_type] + virus = virus_distributions[self.virus_type] scenario_activity_and_expiration = { 'office': ( @@ -315,12 +321,12 @@ class FormData: } [activity_defn, expiration_defn] = scenario_activity_and_expiration[self.activity_type] - activity = models.Activity.types[activity_defn] + activity = activity_distributions[activity_defn] expiration = build_expiration(expiration_defn) infected_occupants = self.infected_people - infected = models.InfectedPopulation( + infected = mc.InfectedPopulation( number=infected_occupants, virus=virus, presence=self.infected_present_interval(), @@ -330,7 +336,7 @@ class FormData: ) return infected - def exposed_population(self) -> models.Population: + def exposed_population(self) -> mc.Population: scenario_activity = { 'office': 'Seated', 'controlroom-day': 'Seated', @@ -345,14 +351,14 @@ class FormData: } activity_defn = scenario_activity[self.activity_type] - activity = models.Activity.types[activity_defn] + activity = activity_distributions[activity_defn] infected_occupants = self.infected_people # The number of exposed occupants is the total number of occupants # minus the number of infected occupants. exposed_occupants = self.total_people - infected_occupants - exposed = models.Population( + exposed = mc.Population( number=exposed_occupants, presence=self.exposed_present_interval(), activity=activity, @@ -560,14 +566,14 @@ def model_from_form(form: FormData) -> models.ExposureModel: room = models.Room(volume=volume, humidity=humidity) # Initializes and returns a model with the attributes defined above - return models.ExposureModel( - concentration_model=models.ConcentrationModel( + return mc.ExposureModel( + concentration_model=mc.ConcentrationModel( room=room, ventilation=form.ventilation(), infected=form.infected_population(), ), exposed=form.exposed_population() - ) + ).build_model(size=_SAMPLE_SIZE) def baseline_raw_form_data(): diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index ea22cbb6..4029e54f 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -38,12 +38,13 @@ def calculate_report_data(model: models.ExposureModel): t_start, t_end = model_start_end(model) times = list(np.linspace(t_start, t_end, resolution)) - concentrations = [model.concentration_model.concentration(time) for time in times] + concentrations = [np.mean(model.concentration_model.concentration(time)) + for time in times] highest_const = max(concentrations) - prob = model.infection_probability() - er = model.concentration_model.infected.emission_rate_when_present() + prob = np.mean(model.infection_probability()) + er = np.mean(model.concentration_model.infected.emission_rate_when_present()) exposed_occupants = model.exposed.number - expected_new_cases = model.expected_new_cases() + expected_new_cases = np.mean(model.expected_new_cases()) repeated_events = [] for n in [1, 2, 3, 4, 5]: @@ -52,8 +53,8 @@ def calculate_report_data(model: models.ExposureModel): repeated_events.append( RepeatEvents( repeats=n, - probability_of_infection=repeat_model.infection_probability(), - expected_new_cases=repeat_model.expected_new_cases(), + probability_of_infection=np.mean(repeat_model.infection_probability()), + expected_new_cases=np.mean(repeat_model.expected_new_cases()), ) ) @@ -127,13 +128,13 @@ def plot(times, concentrations, model: models.ExposureModel): fig = plt.figure() ax = fig.add_subplot(1, 1, 1) datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in times] - ax.plot(datetimes, concentrations, lw=2, color='#1f77b4', label='Concentration') + ax.plot(datetimes, concentrations, lw=2, color='#1f77b4', label='Mean concentration') ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.set_xlabel('Time of day') - ax.set_ylabel('Concentration ($q/m^3$)') - ax.set_title('Concentration of infectious quanta') + ax.set_ylabel('Mean concentration ($q/m^3$)') + ax.set_title('Mean concentration of infectious quanta') ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) # Plot presence of exposed person @@ -243,7 +244,8 @@ def comparison_plot(scenarios: typing.Dict[str, models.ExposureModel]): t_start, t_end = model_start_end(model) times = np.linspace(t_start, t_end, resolution) datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in times] - concentrations = [model.concentration_model.concentration(time) for time in times] + concentrations = [np.mean(model.concentration_model.concentration(time)) + for time in times] if name in dash_styled_scenarios: ax.plot(datetimes, concentrations, label=name, linestyle='--') @@ -257,8 +259,8 @@ def comparison_plot(scenarios: typing.Dict[str, models.ExposureModel]): ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) ax.set_xlabel('Time of day') - ax.set_ylabel('Concentration ($q/m^3$)') - ax.set_title('Concentration of infectious quanta') + ax.set_ylabel('Mean concentration ($q/m^3$)') + ax.set_title('Mean concentration of infectious quanta') return fig @@ -267,8 +269,8 @@ def comparison_report(scenarios: typing.Dict[str, models.ExposureModel]): statistics = {} for name, model in scenarios.items(): statistics[name] = { - 'probability_of_infection': model.infection_probability(), - 'expected_new_cases': model.expected_new_cases(), + 'probability_of_infection': np.mean(model.infection_probability()), + 'expected_new_cases': np.mean(model.expected_new_cases()), } return { 'plot': img2base64(_figure2bytes(comparison_plot(scenarios))), diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index 9b618988..cb57d045 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -209,35 +209,12 @@

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: - - - - - - - - - - - {% for repeat_event in repeated_events %} - - - - - - {% endfor %} - -
# of repeated eventsP(i)Expected new cases
{{ repeat_event.repeats }}{{ repeat_event.probability_of_infection | non_zero_percentage }}{{ repeat_event.expected_new_cases | float_format }}
-

+

Alternative scenarios:

@@ -248,7 +225,7 @@ Scenario - P(i) + P(I) Expected new cases @@ -269,7 +246,7 @@

Notes for alternative scenarios:

  1. This graph shows the concentration of infectious quanta in the air. The filtration of Type I and FFP2 masks, if worn, applies not only to the emission rate but also to the individual exposure (i.e. inhalation). - For this reason, scenarios with different types of mask will show the same concentration on the graph but have different Pi values.
    2. + For this reason, scenarios with different types of mask will show the same concentration on the graph but have different absorbed doses and infection probabilities.
  2. If you have selected more sophisticated options, such as HEPA filtration or FFP2 masks, this will be indicated in the plot as the "base scenario", representing the inputs inserted in the form.
    The other alternative scenarios shown for comparison will not include either HEPA filtration or FFP2 masks.
diff --git a/cara/apps/calculator/templates/userguide.html.j2 b/cara/apps/calculator/templates/userguide.html.j2 index 013a1763..f519a9bb 100644 --- a/cara/apps/calculator/templates/userguide.html.j2 +++ b/cara/apps/calculator/templates/userguide.html.j2 @@ -8,7 +8,7 @@

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

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

-

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

+

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

Disclaimer

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

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

Generate Report

-

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

+

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

Report

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

Results

-

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

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

-

Exposure graph

-

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

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

@@ -204,7 +203,7 @@ This allows for:

Conclusion

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

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

{% endblock contents %} 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 mathematical and physical model simulate the long-range airborne spread of SARS-CoV-2 virus in a finite volume, assuming a homogenous mixture, and estimates the risk of COVID-19 infection therein. The results DO NOT include short-range airborne exposure (where the physical distance plays a factor) nor the other known modes of SARS-CoV-2 transmission. Hence, the output from this model is only valid when the other recommended public health & safety instructions are observed, such as adequate physical distancing, good hand hygiene and other barrier measures.
-

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

+

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

The model used is based on scientific publications relating to airborne transmission of infectious diseases, virology, epidemiology and aerosol science. It can be used to compare the effectiveness of different airborne-related risk mitigation measures. diff --git a/cara/data.py b/cara/data.py index f84e627e..61f103b9 100644 --- a/cara/data.py +++ b/cara/data.py @@ -38,7 +38,7 @@ GenevaTemperatures_hourly = { } # same temperatures on a finer temperature mesh GenevaTemperatures = { - month: GenevaTemperatures_hourly[month].refine(refine_factor=10) + month: GenevaTemperatures_hourly[month].refine(refine_factor=4) for month,temperatures in Geneva_hourly_temperatures_celsius_per_hour.items() } diff --git a/cara/models.py b/cara/models.py index 1561b09d..ad4493ba 100644 --- a/cara/models.py +++ b/cara/models.py @@ -711,7 +711,6 @@ class InfectedPopulation(Population): return self.emission_rate_when_present() - @cached() def emission_rate(self, time) -> _VectorisedFloat: """ The emission rate of the entire population. @@ -741,7 +740,6 @@ class ConcentrationModel: return k + self.virus.decay_constant(self.room.humidity ) + self.ventilation.air_exchange(self.room, time) - @cached() def _concentration_limit(self, time: float) -> _VectorisedFloat: """ Provides a constant that represents the theoretical asymptotic @@ -753,7 +751,6 @@ class ConcentrationModel: return (self.infected.emission_rate(time)) / (IVRR * V) - @cached() def state_change_times(self): """ All time dependent entities on this model must provide information about @@ -798,6 +795,9 @@ class ConcentrationModel: return (self.last_state_change(stop) <= start) @cached() + def _concentration_at_state_change(self, time: float) -> _VectorisedFloat: + return self.concentration(time) + def concentration(self, time: float) -> _VectorisedFloat: """ Virus quanta concentration, as a function of time. @@ -816,7 +816,7 @@ class ConcentrationModel: concentration_limit = self._concentration_limit(next_state_change_time) t_last_state_change = self.last_state_change(time) - concentration_at_last_state_change = self.concentration(t_last_state_change) + concentration_at_last_state_change = self._concentration_at_state_change(t_last_state_change) delta_time = time - t_last_state_change fac = np.exp(-IVRR * delta_time) diff --git a/cara/tests/apps/calculator/test_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/test_known_quantities.py b/cara/tests/test_known_quantities.py index ab24bb08..d8597cf8 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -375,12 +375,13 @@ def test_quanta_hourly_dep(month,expected_quanta): npt.assert_allclose(quanta, expected_quanta) # expected quanta were computed with a trapezoidal integration, using -# a mesh of 100'000 pts per exposed presence interval. +# a mesh of 100'000 pts per exposed presence interval and 25 pts per hour +# for the temperature discretization. @pytest.mark.parametrize( "month, expected_quanta", [ - ['Jan', 9.989881], - ['Jun', 39.99636], + ['Jan', 9.993842], + ['Jun', 40.151985], ], ) def test_quanta_hourly_dep_refined(month,expected_quanta): @@ -393,4 +394,4 @@ def test_quanta_hourly_dep_refined(month,expected_quanta): ) ) quanta = m.quanta_exposure() - npt.assert_allclose(quanta, expected_quanta) + npt.assert_allclose(quanta, expected_quanta, rtol=0.02) diff --git a/cara/tests/test_monte_carlo.py b/cara/tests/test_monte_carlo.py index c3bfc2fe..23e139fd 100644 --- a/cara/tests/test_monte_carlo.py +++ b/cara/tests/test_monte_carlo.py @@ -76,7 +76,9 @@ def test_build_concentration_model(baseline_mc_model: cara.monte_carlo.Concentra model = baseline_mc_model.build_model(7) assert isinstance(model, cara.models.ConcentrationModel) assert isinstance(model.concentration(time=0), float) - assert model.concentration(time=1).shape == (7, ) + conc = model.concentration(time=1) + assert isinstance(conc, np.ndarray) + assert conc.shape == (7, ) def test_build_exposure_model(baseline_mc_exposure_model: cara.monte_carlo.ExposureModel): 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