Merge branch 'feature/monte_carlo_in_calculator' into 'feature/MonteCarlo'

Monte-Carlo in calculator

See merge request cara/cara!195
This commit is contained in:
Andre Henriques 2021-06-07 13:36:53 +00:00
commit ff199020e3
14 changed files with 79 additions and 88 deletions

View file

@ -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():

View file

@ -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))),

View file

@ -209,35 +209,12 @@
<p class="result_title">Results:</p>
<p class="data_text">
{% 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 <b>probability of one exposed occupant getting infected is {{ prob_inf | non_zero_percentage }}</b><a href="#section1">[*]</a> and the <b>expected number of new cases is {{ expected_new_cases | float_format }}</b>.
<p id="section1">[*] The results are based on the parameters and assumptions published in the CERN Open Report <a href="https://cds.cern.ch/record/2756083"> CERN-OPEN-2021-004</a></p>
{% endblock report_summary %}
</p>
<p class="data_title">Exposure graph:</p>
<img id="scenario_concentration_plot" src="{{ scenario_plot_src }}">
<p class="data_title">Repeated events:</p>
<p class="data_text">
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:
<table class="table table-striped w-auto">
<thead class="thead-light">
<tr>
<th># of repeated events</th>
<th>P(i)</th>
<th>Expected new cases</th>
</tr>
</thead>
<tbody>
{% for repeat_event in repeated_events %}
<tr>
<td>{{ repeat_event.repeats }}</td>
<td>{{ repeat_event.probability_of_infection | non_zero_percentage }}</td>
<td style="text-align:right">{{ repeat_event.expected_new_cases | float_format }}</td>
</tr>
{% endfor %}
</tbody>
</table>
</p>
<img id="scenario_concentration_plot" src="{{ scenario_plot_src }}">
<p class="data_title">Alternative scenarios:</p>
<p class="data_text">
@ -248,7 +225,7 @@
<thead class="thead-light">
<tr>
<th>Scenario</th>
<th>P(i)</th>
<th>P(I)</th>
<th>Expected new cases</th>
</tr>
</thead>
@ -269,7 +246,7 @@
<p class="data_text"> <strong> Notes for alternative scenarios: </strong><br>
<ol>
<li>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.</li>
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.</li>
<li>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.<br>
The other alternative scenarios shown for comparison will not include either HEPA filtration or FFP2 masks.</li>
</ol>

View file

@ -8,7 +8,7 @@
<p>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.</p>
<p>For more information on the Airborne Transmission of SARS-CoV-2, feel free to check out the HSE Seminar: <a href="https://cds.cern.ch/record/2743403">https://cds.cern.ch/record/2743403</a></p>
<p>The methodology, mathematical equations and parameters of the model are described here: <a href="https://edms.cern.ch/ui/file/2566402/1/CARA_Deterministic_parameters_2020.pdf">https://edms.cern.ch/ui/file/2566402/1/CARA_Deterministic_parameters_2020.pdf</a></p>
<p>The methodology, mathematical equations and parameters of the model are described here in the CERN Report: <a href="https://cds.cern.ch/record/2756083"> CERN-OPEN-2021-004</a></p>
<h2>Disclaimer</h2>
<p>
@ -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</p>
<p>For the time being only the Type 1 surgical and FFP2 masks can be selected.</p>
<h2>Generate Report</h2>
<p>When you have entered all the necessary information, please click on the Generate Report button to execute the model.</p>
<p>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.</p>
<h1>Report</h1>
<p>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.</p>
<h2>Results</h2>
<p>This part of the report shows the <code>P(i)</code> or probability of one exposed person getting infected.
<p>This part of the report shows the <code>P(I)</code> 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 <code>(1-(1-P(i))^5)</code>.
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 <code>expected number of new cases</code> for the simulation is calculated based on the probability of infection, multiplied by the number of exposed occupants.</p>
<h3>Exposure graph</h3>
<p>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.
<p>The graph shows the variation in the concentration of infectious viruses within the simulated volume.
It is determined by:</p>
<ul>
<li>The presence of the infected person, who emits airborne viruses in the volume.</li>
<li>The emission rate is related to the type of activity of the infected person (sitting, light exercise), their level of vocalisation (breathing, whispering or talking).</li>
<li>The emission rate is related to the type of activity of the infected person (sitting, light exercise), their level of vocalisation (breathing, talking or shouting).</li>
<li>The accumulation of infectious quanta in the volume, which is driven, among other factors, by ventilation (if applicable).<ul>
<li>In a mechanical ventilation scenario, the removal rate is constant, based on fresh airflow supply in and out of the simulated space.</li>
<li>Under natural ventilation conditions, the effectiveness of ventilation relies upon the hourly temperature difference between the inside and outside air temperature.</li>
<li>A HEPA filter removes infectious quanta from the air at a constant rate and is modelled in the same way as mechanical ventilation, however air passed through a HEPA filter is recycled (i.e. it is not fresh air).</li>
<li>A HEPA filter removes infectious virus from the air at a constant rate and is modelled in the same way as mechanical ventilation, however air passed through a HEPA filter is recycled (i.e. it is not fresh air).</li>
</ul>
</li>
</ul>
@ -204,7 +203,7 @@ This allows for:</p>
</ul>
<h1>Conclusion</h1>
<p>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.</p>
If you have any comments on your experience with the app, or feedback for potential improvements, please share them with the development team <a href="mailto:cara-dev@cern.ch">Send email</a>.</p>
{% endblock contents %}

View file

@ -3,7 +3,7 @@
{% block report_preamble %}
<p><strong>Applicable rules: <br>
Please ensure that this scenario conforms to current <a href="https://hse.cern/covid-19-information"> CERN HSE rules</a> (minimum ventilation requirements, mask wearing and the maximum number of people permitted in a space).</strong> <br>
Please ensure that this scenario conforms to current COVID-related <a href="https://hse.cern/covid-19-information"> Health & Safety requirements</a>, under the applicable COVID Scale and measures in force at the time of the CARA assessment.</strong> <br>
The results of this simulation are colour coded according to the risk values authorized at CERN (approved in December 2020):
<ul><li>Events with a <span class="green_bkg">P(i) less than 5%</span> may go ahead without further mitigation measures.</li>
<li>Events with a <span class="yellow_bkg">P(i) between 5% and 15%</span> shall be subject to ALARA principles (see footnote) to minimise the risk before proceeding.</li>

View file

@ -17,7 +17,7 @@ CARA stands for COVID Airborne Risk Assessment and was developed in the spring o
</ul>
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.<br>
<p>The methodology, mathematical equations and parameters of the model are described here: <a href="https://edms.cern.ch/ui/file/2566402/1/CARA_Deterministic_parameters_2020.pdf">https://edms.cern.ch/ui/file/2566402/1/CARA_Deterministic_parameters_2020.pdf</a>.</p>
<p>The methodology, mathematical equations and parameters of the model are described here in the CERN Report: <a href="https://cds.cern.ch/record/2756083"> CERN-OPEN-2021-004</a>.</p>
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.

View file

@ -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()
}

View file

@ -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)

View file

@ -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])

View file

@ -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()

View file

@ -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):

View file

@ -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)

View file

@ -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):

View file

@ -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