From c27f567f3ab5aa936b528eb3460cfb0de5d48ae7 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 7 Sep 2021 17:00:51 +0200 Subject: [PATCH] merge with master --- cara/apps/calculator/model_generator.py | 7 +- cara/apps/calculator/report_generator.py | 111 ++------ cara/apps/calculator/static/js/report.js | 257 ++++++++++++++---- .../templates/base/calculator.report.html.j2 | 20 +- cara/apps/templates/common_text.md.j2 | 3 + cara/models.py | 1 - cara/monte_carlo/data.py | 12 +- cara/monte_carlo/sampleable.py | 12 + 8 files changed, 264 insertions(+), 159 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index a6a41fab..06c29229 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -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: diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 10bd67ee..23ad9fb1 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -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'. 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: diff --git a/cara/apps/calculator/static/js/report.js b/cara/apps/calculator/static/js/report.js index 00e85892..6266a09b 100644 --- a/cara/apps/calculator/static/js/report.js +++ b/cara/apps/calculator/static/js/report.js @@ -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('Mean concentration of virions'); + 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; } \ No newline at end of file diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index 371788c0..108a9a9c 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -86,12 +86,12 @@ {% endblock report_summary_footnote %}

* The results are based on the parameters and assumptions published in the CERN Open Report CERN-OPEN-2021-004.

- +

@@ -109,8 +109,12 @@
- - + + {% block report_scenarios_summary_table %} @@ -437,4 +441,4 @@ - \ No newline at end of file + diff --git a/cara/apps/templates/common_text.md.j2 b/cara/apps/templates/common_text.md.j2 index 63512dee..ac8829d7 100644 --- a/cara/apps/templates/common_text.md.j2 +++ b/cara/apps/templates/common_text.md.j2 @@ -71,3 +71,6 @@ We wish to thank CERN’s 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.
[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.
[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.
+[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.
+[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.
+[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.
\ No newline at end of file diff --git a/cara/models.py b/cara/models.py index 772c0324..a86457eb 100644 --- a/cara/models.py +++ b/cara/models.py @@ -904,7 +904,6 @@ class ExposureModel: return 0 def inhaled_exposure_between_bounds(self, time: float) -> _VectorisedFloat: - exposure = self.exposure_between_bounds(time) return ( diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 613b830c..4415efbb 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -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)), +} \ No newline at end of file diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 4333c93f..27907e49 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -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