merge with master

This commit is contained in:
Luis Aleixo 2021-09-07 17:00:51 +02:00
parent 0b83c6491e
commit c27f567f3a
8 changed files with 264 additions and 159 deletions

View file

@ -11,7 +11,7 @@ from cara import data
import cara.data.weather
import cara.monte_carlo as mc
from .. import calculator
from cara.monte_carlo.data import activity_distributions, virus_distributions
from cara.monte_carlo.data import activity_distributions, virus_distributions, mask_distributions
LOG = logging.getLogger(__name__)
@ -359,7 +359,10 @@ class FormData:
def mask(self) -> models.Mask:
# Initializes the mask type if mask wearing is "continuous", otherwise instantiates the mask attribute as
# the "No mask"-mask
mask = models.Mask.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask']
if self.mask_wearing_option == 'mask_on':
mask = mask_distributions[self.mask_type]
else:
mask = models.Mask.types['No mask']
return mask
def infected_population(self) -> mc.InfectedPopulation:

View file

@ -1,10 +1,3 @@
from ... import dataclass_utils
from .model_generator import FormData, _DEFAULT_MC_SAMPLE_SIZE
from ... import monte_carlo as mc
from cara import models
import qrcode
import numpy as np
import matplotlib.pyplot as plt
import concurrent.futures
import base64
import dataclasses
@ -16,9 +9,14 @@ import zlib
import loky
import jinja2
import matplotlib
from numpy.lib.function_base import append, quantile
matplotlib.use('agg')
import numpy as np
import qrcode
import json
from cara import models
from ... import monte_carlo as mc
from .model_generator import FormData, _DEFAULT_MC_SAMPLE_SIZE
from ... import dataclass_utils
def model_start_end(model: models.ExposureModel):
@ -96,8 +94,7 @@ def interesting_times(model: models.ExposureModel, approx_n_pts=100) -> typing.L
# Expand the times list to ensure that we have a maximum gap size between
# the key times.
nice_times = fill_big_gaps(times, gap_size=(
max(times) - min(times)) / approx_n_pts)
nice_times = fill_big_gaps(times, gap_size=(max(times) - min(times)) / approx_n_pts)
return nice_times
@ -110,8 +107,7 @@ def calculate_report_data(model: models.ExposureModel):
]
highest_const = max(concentrations)
prob = np.array(model.infection_probability()).mean()
er = np.array(
model.concentration_model.infected.emission_rate_when_present()).mean()
er = np.array(model.concentration_model.infected.emission_rate_when_present()).mean()
exposed_occupants = model.exposed.number
expected_new_cases = np.array(model.expected_new_cases()).mean()
cumulative_doses = [
@ -169,21 +165,12 @@ def _img2bytes(figure):
return img_data
def _figure2bytes(figure):
# Draw the image
img_data = io.BytesIO()
figure.savefig(img_data, format='png',
bbox_inches="tight", transparent=True)
return img_data
def img2base64(img_data) -> str:
plt.close()
img_data.seek(0)
pic_hash = base64.b64encode(img_data.read()).decode('ascii')
# A src suitable for a tag such as f'<img id="scenario_concentration_plot" src="{result}">.
return f'data:image/png;base64,{pic_hash}'
def minutes_to_time(minutes: int) -> str:
minute_string = str(minutes % 60)
@ -224,18 +211,14 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp
scenarios = {}
# Two special option cases - HEPA and/or FFP2 masks.
FFP2_being_worn = bool(form.mask_wearing_option ==
'mask_on' and form.mask_type == 'FFP2')
FFP2_being_worn = bool(form.mask_wearing_option == 'mask_on' and form.mask_type == 'FFP2')
if FFP2_being_worn and form.hepa_option:
FFP2andHEPAalternative = dataclass_utils.replace(
form, mask_type='Type I')
FFP2andHEPAalternative = dataclass_utils.replace(form, mask_type='Type I')
scenarios['Base scenario with HEPA filter and Type I masks'] = FFP2andHEPAalternative.build_mc_model()
if not FFP2_being_worn and form.hepa_option:
noHEPAalternative = dataclass_utils.replace(form, mask_type='FFP2')
noHEPAalternative = dataclass_utils.replace(
noHEPAalternative, mask_wearing_option='mask_on')
noHEPAalternative = dataclass_utils.replace(
noHEPAalternative, hepa_option=False)
noHEPAalternative = dataclass_utils.replace(form, mask_type = 'FFP2')
noHEPAalternative = dataclass_utils.replace(noHEPAalternative, mask_wearing_option = 'mask_on')
noHEPAalternative = dataclass_utils.replace(noHEPAalternative, hepa_option=False)
scenarios['Base scenario without HEPA filter, with FFP2 masks'] = noHEPAalternative.build_mc_model()
# The remaining scenarios are based on Type I masks (possibly not worn)
@ -245,8 +228,7 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp
form = dataclass_utils.replace(form, hepa_option=False)
with_mask = dataclass_utils.replace(form, mask_wearing_option='mask_on')
without_mask = dataclass_utils.replace(
form, mask_wearing_option='mask_off')
without_mask = dataclass_utils.replace(form, mask_wearing_option='mask_off')
if form.ventilation_type == 'mechanical_ventilation':
#scenarios['Mechanical ventilation with Type I masks'] = with_mask.build_mc_model()
@ -257,62 +239,17 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp
scenarios['Windows open without masks'] = without_mask.build_mc_model()
# No matter the ventilation scheme, we include scenarios which don't have any ventilation.
with_mask_no_vent = dataclass_utils.replace(
with_mask, ventilation_type='no_ventilation')
without_mask_or_vent = dataclass_utils.replace(
without_mask, ventilation_type='no_ventilation')
with_mask_no_vent = dataclass_utils.replace(with_mask, ventilation_type='no_ventilation')
without_mask_or_vent = dataclass_utils.replace(without_mask, ventilation_type='no_ventilation')
scenarios['No ventilation with Type I masks'] = with_mask_no_vent.build_mc_model()
scenarios['Neither ventilation nor masks'] = without_mask_or_vent.build_mc_model()
return scenarios
def comparison_plot(scenarios: typing.Dict[str, dict], sample_times: typing.List[float]):
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax1 = ax.twinx()
dash_styled_scenarios = [
'Base scenario with FFP2 masks',
'Base scenario with HEPA filter',
'Base scenario with HEPA and FFP2 masks',
]
datetimes = [datetime(1970, 1, 1) + timedelta(hours=time)
for time in sample_times]
for name, statistics in scenarios.items():
model = statistics['model']
concentrations = statistics['concentrations']
# Plot concentrations
if name in dash_styled_scenarios:
ax.plot(datetimes, concentrations, label=name, linestyle='--')
else:
ax.plot(datetimes, concentrations,
label=name, linestyle='-', alpha=0.5)
# Place a legend outside of the axes itself.
fig.legend(bbox_to_anchor=(1.05, 0.95), loc='upper left')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel('Time of day', fontsize=14)
ax.set_ylabel('Mean viral concentration\n(virion m$^{-3}$)', fontsize=14)
ax.set_title('Concentration profile')
ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M"))
ax.set_xlabel('Time of day')
ax.set_ylabel('Mean concentration ($virions/m^{3}$)')
ax.set_title('Mean concentration of virions')
return fig
def scenario_statistics(mc_model: mc.ExposureModel, sample_times: typing.List[float]):
def scenario_statistics(mc_model: mc.ExposureModel, sample_times: np.ndarray):
model = mc_model.build_model(size=_DEFAULT_MC_SAMPLE_SIZE)
return {
'model': model,
'probability_of_infection': np.mean(model.infection_probability()),
'expected_new_cases': np.mean(model.expected_new_cases()),
'concentrations': [
@ -339,7 +276,6 @@ def comparison_report(
for (name, model), model_stats in zip(scenarios.items(), results):
statistics[name] = model_stats
return {
'plot': img2base64(_figure2bytes(comparison_plot(statistics, sample_times))),
'stats': statistics,
}
@ -356,8 +292,7 @@ class ReportGenerator:
executor_factory: typing.Callable[[], concurrent.futures.Executor],
) -> str:
model = form.build_model()
context = self.prepare_context(
base_url, model, form, executor_factory=executor_factory)
context = self.prepare_context(base_url, model, form, executor_factory=executor_factory)
return self.render(context)
def prepare_context(
@ -383,8 +318,7 @@ class ReportGenerator:
context['alternative_scenarios'] = comparison_report(
alternative_scenarios, scenario_sample_times, executor_factory=executor_factory,
)
context['qr_code'] = generate_qr_code(
base_url, self.calculator_prefix, form)
context['qr_code'] = generate_qr_code(base_url, self.calculator_prefix, form)
context['calculator_prefix'] = self.calculator_prefix
context['scale_warning'] = {
'level': 'yellow-2',
@ -404,6 +338,7 @@ class ReportGenerator:
env.filters['minutes_to_time'] = minutes_to_time
env.filters['float_format'] = "{0:.2f}".format
env.filters['int_format'] = "{:0.0f}".format
env.filters['JSONify'] = json.dumps
return env
def render(self, context: dict) -> str:

View file

@ -27,70 +27,18 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses
yCumulatedAxis = d3.axisRight(yCumulatedRange);
// Plot tittle.
vis.append('svg:foreignObject')
.attr("background-color", "transparent")
.attr('width', width)
.attr('height', margins.top)
.style('text-align', 'center')
.html('<b>Mean concentration of virions</b>');
plot_title(vis, width, margins.top, 'Mean concentration of virions');
// Line representing the mean concentration.
var lineFunc = d3.line()
.defined(d => !isNaN(d.concentration))
.x(d => xTimeRange(d.time))
.y(d => yRange(d.concentration))
.curve(d3.curveBasis);
vis.append('svg:path')
.attr('d', lineFunc(data))
.attr('stroke', '#1f77b4')
.attr('stroke-width', 2)
.attr('fill', 'none');
plot_scenario_data(vis, data, xTimeRange, yRange, '#1f77b4');
// Line representing the cumulative concentration.
var lineCumulativeFunc = d3.line()
.defined(d => !isNaN(d.cumulative_doses))
.x(d => xTimeRange(d.time))
.y(d => yCumulatedRange(d.cumulative_doses))
.curve(d3.curveBasis);
plot_cumulative_data(vis, data, xTimeRange, yCumulatedRange, '#1f77b4');
vis.append('svg:path')
.attr('d', lineCumulativeFunc(data))
.attr('stroke', '#1f77b4')
.attr('stroke-width', 2)
.style("stroke-dasharray", "5 5")
.attr('fill', 'none');
// X axis.
plot_x_axis(vis, height, width, margins, xAxis, 'Time of day');
// X axis declaration.
vis.append('svg:g')
.attr('class', 'x axis')
.attr('transform', 'translate(0,' + (height - margins.bottom) + ')')
.call(xAxis);
// X axis label.
vis.append('text')
.attr('class', 'x label')
.attr('fill', 'black')
.attr('text-anchor', 'middle')
.attr('x', (width + margins.right) / 2)
.attr('y', height * 0.97)
.text('Time of day')
// Y concentration axis declaration.
vis.append('svg:g')
.attr('class', 'y axis')
.attr('transform', 'translate(' + margins.left + ',0)')
.call(yAxis);
// Y concentration axis label.
vis.append('svg:text')
.attr('class', 'y label')
.attr('fill', 'black')
.attr('transform', 'rotate(-90, 0,' + height + ')')
.attr('text-anchor', 'middle')
.attr('x', (height + margins.bottom) / 2)
.attr('y', (height + margins.left) * 0.92)
.text('Mean concentration (virions/m³)');
// Y axis
plot_y_axis(vis, height, width, margins, yAxis, 'Mean concentration (virions/m³)')
// Y cumulative concentration axis declaration.
vis.append('svg:g')
@ -230,4 +178,195 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses
focus.select('#tooltip-time').text('x = ' + time_format(d.hour));
focus.select('#tooltip-concentration').text('y = ' + d.concentration.toFixed(2));
}
}
// Generate the alternative scenarios plot using d3 library.
// 'alternative_scenarios' is a dictionary with all the alternative scenarios
// 'times' is a list of times for all the scenarios
function draw_alternative_scenarios_plot(svg_id, width, height, alternative_scenarios, times) {
// H:M format
var time_format = d3.timeFormat('%H:%M');
// D3 array of ten categorical colors represented as RGB hexadecimal strings.
var colors = d3.schemeAccent;
// Variable for the highest concentration for all the scenarios
var highest_concentration = 0.
var data_for_scenarios = {}
for (scenario in alternative_scenarios) {
scenario_concentrations = alternative_scenarios[scenario].concentrations
highest_concentration = Math.max(highest_concentration, Math.max(...scenario_concentrations))
var data = []
times.map((time, index) => data.push({ 'time': time, 'hour': new Date().setHours(Math.trunc(time), (time - Math.trunc(time)) * 60), 'concentration': scenario_concentrations[index] }))
// Add data into lines dictionary
data_for_scenarios[scenario] = data
}
// We need one scenario to get the time range
var first_scenario = Object.values(data_for_scenarios)[0]
var vis = d3.select(svg_id),
width = width,
height = height,
margins = { top: 30, right: 20, bottom: 50, left: 50 },
// H:M time format for x axis.
xRange = d3.scaleTime().range([margins.left, width - margins.right]).domain([first_scenario[0].hour, first_scenario[first_scenario.length - 1].hour]),
xTimeRange = d3.scaleLinear().range([margins.left, width - margins.right]).domain([times[0], times[times.length - 1]]),
yRange = d3.scaleLinear().range([height - margins.bottom, margins.top]).domain([0., highest_concentration]),
xAxis = d3.axisBottom(xRange).tickFormat(d => time_format(d)),
yAxis = d3.axisLeft(yRange);
// Plot title.
plot_title(vis, width, margins.top, 'Mean concentration of virions');
// Line representing the mean concentration for each scenario.
for (const [scenario_name, data] of Object.entries(data_for_scenarios)) {
var scenario_index = Object.keys(data_for_scenarios).indexOf(scenario_name)
// Line representing the mean concentration.
plot_scenario_data(vis, data, xTimeRange, yRange, colors[scenario_index])
// Legend for the plot elements - lines.
var size = 20 * (scenario_index + 1)
vis.append('rect')
.attr('x', width + 20)
.attr('y', margins.top + size)
.attr('width', 20)
.attr('height', 3)
.style('fill', colors[scenario_index]);
vis.append('text')
.attr('x', width + 3 * 20)
.attr('y', margins.top + size)
.text(scenario_name)
.style('font-size', '15px')
.attr('alignment-baseline', 'central');
}
// X axis.
plot_x_axis(vis, height, width, margins, xAxis, "Time of day");
// Y axis declaration.
vis.append('svg:g')
.attr('class', 'y axis')
.attr('transform', 'translate(' + margins.left + ',0)')
.call(yAxis);
// Y axis label.
vis.append('svg:text')
.attr('class', 'y label')
.attr('fill', 'black')
.attr('transform', 'rotate(-90, 0,' + height + ')')
.attr('text-anchor', 'middle')
.attr('x', (height + margins.bottom) / 2)
.attr('y', (height + margins.left) * 0.92)
.text('Mean concentration (virions/m³)');
// Legend bounding box.
vis.append('rect')
.attr('width', 275)
.attr('height', 25 * (Object.keys(data_for_scenarios).length))
.attr('x', width * 1.005)
.attr('y', margins.top + 5)
.attr('stroke', 'lightgrey')
.attr('stroke-width', '2')
.attr('rx', '5px')
.attr('ry', '5px')
.attr('stroke-linejoin', 'round')
.attr('fill', 'none');
}
// Functions used to build the plots' components
function plot_title(vis, width, margin_top, title) {
vis.append('svg:foreignObject')
.attr('width', width)
.attr('height', margin_top)
.attr('fill', 'none')
.append('xhtml:div')
.style('text-align', 'center')
.html(title);
return vis;
}
function plot_x_axis(vis, height, width, margins, xAxis, label) {
// X axis declaration
vis.append('svg:g')
.attr('class', 'x axis')
.attr('transform', 'translate(0,' + (height - margins.bottom) + ')')
.call(xAxis);
// X axis label.
vis.append('text')
.attr('class', 'x label')
.attr('fill', 'black')
.attr('text-anchor', 'middle')
.attr('x', (width + margins.right) / 2)
.attr('y', height * 0.97)
.text(label);
return vis;
}
function plot_y_axis(vis, height, width, margins, yAxis, label) {
// Y axis declaration.
vis.append('svg:g')
.attr('class', 'y axis')
.attr('transform', 'translate(' + margins.left + ',0)')
.call(yAxis);
// Y axis label.
vis.append('svg:text')
.attr('class', 'y label')
.attr('fill', 'black')
.attr('transform', 'rotate(-90, 0,' + height + ')')
.attr('text-anchor', 'middle')
.attr('x', (height + margins.bottom) / 2)
.attr('y', (height + margins.left) * 0.92)
.text(label);
return vis;
}
function plot_scenario_data(vis, data, xTimeRange, yRange, line_color) {
var lineFunc = d3.line()
.defined(d => !isNaN(d.concentration))
.x(d => xTimeRange(d.time))
.y(d => yRange(d.concentration))
.curve(d3.curveBasis);
vis.append('svg:path')
.attr('d', lineFunc(data))
.attr("stroke", line_color)
.attr('stroke-width', 2)
.attr('fill', 'none');
return vis;
}
function plot_cumulative_data(vis, data, xTimeRange, yCumulativeRange, line_color) {
var lineCumulativeFunc = d3.line()
.defined(d => !isNaN(d.cumulative_doses))
.x(d => xTimeRange(d.time))
.y(d => yCumulativeRange(d.cumulative_doses))
.curve(d3.curveBasis);
vis.append('svg:path')
.attr('d', lineCumulativeFunc(data))
.attr('stroke', line_color)
.attr('stroke-width', 2)
.style("stroke-dasharray", "5 5")
.attr('fill', 'none');
return vis;
}

View file

@ -86,12 +86,12 @@
{% endblock report_summary_footnote %}
</div>
<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>
<svg id="result_plot" width="1000" height="400"></svg>
<svg id="result_plot" width="900" height="400"></svg>
<script type="application/javascript">
var times = {{times}}
var concentrations = {{concentrations}}
var cumulative_doses = {{cumulative_doses}}
var exposed_presence_intervals = {{exposed_presence_intervals}}
var times = {{ times | JSONify }}
var concentrations = {{ concentrations | JSONify }}
var cumulative_doses = {{ cumulative_doses | JSONify }}
var exposed_presence_intervals = {{ exposed_presence_intervals | JSONify }}
draw_concentration_plot("#result_plot", times, concentrations, cumulative_doses, exposed_presence_intervals);
</script>
</p>
@ -109,8 +109,12 @@
<div class="collapse" id="collapseAlternativeScenarios">
<div class="card-body">
<div>
<img id="scenario_concentration_plot" src="{{ alternative_scenarios.plot }}" />
<svg id="alternative_scenario_plot" width="900" height="400"></svg>
<script type="application/javascript">
var alternative_scenarios = {{ alternative_scenarios.stats | JSONify }}
var times = {{ times | JSONify }}
draw_alternative_scenarios_plot("#alternative_scenario_plot", width=600, height=400, alternative_scenarios, times);
</script>
{% block report_scenarios_summary_table %}
<table class="table w-auto">
<thead class="thead-light">
@ -437,4 +441,4 @@
<script src="https://cdnjs.cloudflare.com/ajax/libs/html2pdf.js/0.9.2/html2pdf.bundle.js"></script>
</body>
</html>
</html>

View file

@ -71,3 +71,6 @@ We wish to thank CERNs HSE Unit, Beams Department, Experimental Physics Depar
[54] Leung, N.H.L et al. Respiratory virus shedding in exhaled breath and efficacy of face masks. Nat Med (2020). 10.1038/s41591-020-0843-2.<br>
[55] Asadi, S., Cappa, C.D., Barreda, S. et al. Efficacy of masks and face coverings in controlling outward aerosol particle emission from expiratory activities. Sci Rep 10, 15665 (2020). https://doi.org/10.1038/s41598-020-72798-7.<br>
[56] Endo A, Abbott S et al. Estimating the overdispersion in COVID-19 transmission using outbreak sizes outside China [version 3; peer review: 2 approved]. Wellcome Open Res 2020, 5:67. doi:10.12688/wellcomeopenres.15842.3.<br>
[57] Jin Pan, Charbel Harb, Weinan Leng & Linsey C. Marr (2021) Inward and outward effectiveness of cloth masks, a surgical mask, and a face shield, Aerosol Science and Technology, 55:6, 718-733, doi: 10.1080/02786826.2021.1890687.<br>
[58] C. Makison Booth, M. Clayton, B. Crook, J.M. Gawn, Effectiveness of surgical masks against influenza bioaerosols, Journal of Hospital Infection, Volume 84, Issue 1, 2013, Pages 22-26, https://doi.org/10.1016/j.jhin.2013.02.007.<br>
[59] Riediker, M., Monn, C. (2021). Simulation of SARS-CoV-2 Aerosol Emissions in the Infected Population and Resulting Airborne Exposures in Different Indoor Scenarios. Aerosol Air Qual. Res. 21, 200531. https://doi.org/10.4209/aaqr.2020.08.0531.<br>

View file

@ -904,7 +904,6 @@ class ExposureModel:
return 0
def inhaled_exposure_between_bounds(self, time: float) -> _VectorisedFloat:
exposure = self.exposure_between_bounds(time)
return (

View file

@ -1,7 +1,7 @@
import numpy as np
import cara.monte_carlo as mc
from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel
from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel, Uniform
# From CERN-OPEN-2021-04 and refererences therein
@ -58,3 +58,13 @@ virus_distributions = {
infectious_dose=60/1.6,
),
}
# From:
# https://doi.org/10.1080/02786826.2021.1890687
# https://doi.org/10.1016/j.jhin.2013.02.007
# https://doi.org/10.4209/aaqr.2020.08.0531
mask_distributions = {
'Type I': mc.Mask(Uniform(0.25, 0.80)),
'FFP2': mc.Mask(Uniform(0.83, 0.91)),
}

View file

@ -28,6 +28,18 @@ class Normal(SampleableDistribution):
return np.random.normal(self.mean, self.standard_deviation, size=size)
class Uniform(SampleableDistribution):
"""
Defines a continuous uniform distribution
"""
def __init__(self, low: float, high: float):
self.low = low
self.high = high
def generate_samples(self, size: int) -> float_array_size_n:
return np.random.uniform(self.low, self.high, size=size)
class LogNormal(SampleableDistribution):
"""
Defines a lognormal distribution (i.e. Gaussian distribution vs. the