diff --git a/caimira/apps/calculator/__init__.py b/caimira/apps/calculator/__init__.py index 8399bd9d..3da5f048 100644 --- a/caimira/apps/calculator/__init__.py +++ b/caimira/apps/calculator/__init__.py @@ -35,7 +35,7 @@ from .user import AuthenticatedUser, AnonymousUser # calculator version. If the calculator needs to make breaking changes (e.g. change # form attributes) then it can also increase its MAJOR version without needing to # increase the overall CAiMIRA version (found at ``caimira.__version__``). -__version__ = "4.7" +__version__ = "4.8" class BaseRequestHandler(RequestHandler): @@ -122,6 +122,10 @@ class ConcentrationModel(BaseRequestHandler): max_workers=self.settings['handler_worker_pool_size'], timeout=300, ) + # Re-generate the report with the conditional probability of infection plot + if self.get_cookie('conditional_plot'): + form.conditional_probability_plot = True if self.get_cookie('conditional_plot') == '1' else False + self.clear_cookie('conditional_plot') # Clears cookie after changing the form value. report_task = executor.submit( report_generator.build_report, base_url, form, executor_factory=functools.partial( diff --git a/caimira/apps/calculator/model_generator.py b/caimira/apps/calculator/model_generator.py index 5931c980..5c0fe232 100644 --- a/caimira/apps/calculator/model_generator.py +++ b/caimira/apps/calculator/model_generator.py @@ -36,6 +36,7 @@ class FormData: specific_breaks: dict precise_activity: dict ceiling_height: float + conditional_probability_plot: bool exposed_coffee_break_option: str exposed_coffee_duration: int exposed_finish: minutes_since_midnight @@ -104,6 +105,7 @@ class FormData: 'precise_activity': '{}', 'calculator_version': _NO_DEFAULT, 'ceiling_height': 0., + 'conditional_probability_plot': False, 'exposed_coffee_break_option': 'coffee_break_0', 'exposed_coffee_duration': 5, 'exposed_finish': '17:30', @@ -900,6 +902,7 @@ def baseline_raw_form_data() -> typing.Dict[str, typing.Union[str, float]]: 'air_changes': '', 'air_supply': '', 'ceiling_height': '', + 'conditional_probability_plot': '0', 'exposed_coffee_break_option': 'coffee_break_4', 'exposed_coffee_duration': '10', 'exposed_finish': '18:00', diff --git a/caimira/apps/calculator/report_generator.py b/caimira/apps/calculator/report_generator.py index 1dcfa00d..39919fac 100644 --- a/caimira/apps/calculator/report_generator.py +++ b/caimira/apps/calculator/report_generator.py @@ -10,6 +10,7 @@ import zlib import jinja2 import numpy as np +import matplotlib.pyplot as plt from caimira import models from caimira.apps.calculator import markdown_tools @@ -130,11 +131,13 @@ def calculate_report_data(form: FormData, model: models.ExposureModel) -> typing for time1, time2 in zip(times[:-1], times[1:]) ]) - prob = np.array(model.infection_probability()).mean() + prob = np.array(model.infection_probability()) + prob_dist_count, prob_dist_bins = np.histogram(prob/100, bins=100, density=True) prob_probabilistic_exposure = np.array(model.total_probability_rule()).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() + uncertainties_plot_src = img2base64(_figure2bytes(uncertainties_plot(model))) if form.conditional_probability_plot else None return { "model_repr": repr(model), @@ -147,11 +150,16 @@ def calculate_report_data(form: FormData, model: models.ExposureModel) -> typing "highest_const": highest_const, "cumulative_doses": list(cumulative_doses), "long_range_cumulative_doses": list(long_range_cumulative_doses), - "prob_inf": prob, + "prob_inf": prob.mean(), + "prob_inf_sd": np.std(prob), + "prob_dist": list(prob), + "prob_hist_count": list(prob_dist_count), + "prob_hist_bins": list(prob_dist_bins), "prob_probabilistic_exposure": prob_probabilistic_exposure, "emission_rate": er, "exposed_occupants": exposed_occupants, "expected_new_cases": expected_new_cases, + "uncertainties_plot_src": uncertainties_plot_src, } @@ -174,6 +182,68 @@ def generate_permalink(base_url, get_root_url, get_root_calculator_url, form: F } +def uncertainties_plot(exposure_model: models.ExposureModel): + fig = plt.figure(figsize=(4, 7), dpi=110) + + viral_loads = np.linspace(2, 10, 600) + pi_means, lower_percentiles, upper_percentiles = [], [], [] + for vl in viral_loads: + model_vl = dataclass_utils.nested_replace( + exposure_model, { + 'concentration_model.infected.virus.viral_load_in_sputum' : 10**vl, + } + ) + pi = model_vl.infection_probability()/100 + + pi_means.append(np.mean(pi)) + lower_percentiles.append(np.quantile(pi, 0.05)) + upper_percentiles.append(np.quantile(pi, 0.95)) + + histogram_data = exposure_model.infection_probability() / 100 + + fig, axs = plt.subplots(2, 3, + gridspec_kw={'width_ratios': [5, 0.5] + [1], + 'height_ratios': [3, 1], 'wspace': 0}, + sharey='row', + sharex='col') + + for y, x in [(0, 1)] + [(1, i + 1) for i in range(2)]: + axs[y, x].axis('off') + + axs[0, 1].set_visible(False) + + axs[0, 0].plot(viral_loads, pi_means, label='Predictive total probability') + axs[0, 0].fill_between(viral_loads, lower_percentiles, upper_percentiles, alpha=0.1, label='5ᵗʰ and 95ᵗʰ percentile') + + axs[0, 2].hist(histogram_data, bins=30, orientation='horizontal') + axs[0, 2].set_xticks([]) + axs[0, 2].set_xticklabels([]) + axs[0, 2].set_facecolor("lightgrey") + + highest_bar = axs[0, 2].get_xlim()[1] + axs[0, 2].set_xlim(0, highest_bar) + + axs[0, 2].text(highest_bar * 0.5, 0.5, + rf"$\bf{np.round(np.mean(histogram_data) * 100, 1)}$%", ha='center', va='center') + axs[1, 0].hist(np.log10(exposure_model.concentration_model.infected.virus.viral_load_in_sputum), + bins=150, range=(2, 10), color='grey') + axs[1, 0].set_facecolor("lightgrey") + axs[1, 0].set_yticks([]) + axs[1, 0].set_yticklabels([]) + axs[1, 0].set_xticks([i for i in range(2, 13, 2)]) + axs[1, 0].set_xticklabels(['$10^{' + str(i) + '}$' for i in range(2, 13, 2)]) + axs[1, 0].set_xlim(2, 10) + axs[1, 0].set_xlabel('Viral load\n(RNA copies)', fontsize=12) + axs[0, 0].set_ylabel('Conditional Probability\nof Infection', fontsize=12) + + axs[0, 0].text(9.5, -0.01, '$(i)$') + axs[1, 0].text(9.5, axs[1, 0].get_ylim()[1] * 0.8, '$(ii)$') + axs[0, 2].set_title('$(iii)$', fontsize=10) + + axs[0, 0].legend() + return fig + + def _img2bytes(figure): # Draw the image img_data = io.BytesIO() @@ -181,6 +251,13 @@ 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, dpi=110) + return img_data + + def img2base64(img_data) -> str: img_data.seek(0) pic_hash = base64.b64encode(img_data.read()).decode('ascii') diff --git a/caimira/apps/calculator/static/js/report.js b/caimira/apps/calculator/static/js/report.js index 26e41cc3..49286914 100644 --- a/caimira/apps/calculator/static/js/report.js +++ b/caimira/apps/calculator/static/js/report.js @@ -1,3 +1,8 @@ +function on_report_load(conditional_probability_plot) { + // Check/uncheck uncertainties image generation + document.getElementById('conditional_probability_plot').checked = conditional_probability_plot +} + /* Generate the concentration plot using d3 library. */ function draw_plot(svg_id) { @@ -537,7 +542,6 @@ function draw_plot(svg_id) { }); } - // 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 @@ -853,6 +857,196 @@ function draw_alternative_scenarios_plot(concentration_plot_svg_id, alternative_ }); } +function draw_histogram(svg_id, prob, prob_sd) { + // Add main SVG element + var plot_div = document.getElementById(svg_id); + var div_width = plot_div.clientWidth; + var div_height = plot_div.clientHeight; + var vis = d3.select(plot_div).append('svg'); + + // set the dimensions and margins of the graph + if (div_width > 900) { + div_width = 900; + var margins = { top: 30, right: 20, bottom: 50, left: 60 }; + var graph_width = div_width * (2/3); + const svg_margins = {'margin-left': '0rem'}; + Object.entries(svg_margins).forEach(([prop,val]) => vis.style(prop,val)); + } + + vis.attr("width", div_width).attr('height', div_height); + + let hist_count = prob_hist_count; + let hist_bins = prob_hist_bins; + + // X axis: scale and draw: + var x = d3.scaleLinear() + .domain([0, d3.max(hist_bins)]) + .range([margins.left, graph_width - margins.right]); + vis.append("svg:g") + .attr("transform", "translate(0," + (graph_height - margins.bottom) + ")") + .call(d3.axisBottom(x)); + + // X axis label. + vis.append('text') + .attr('class', 'x label') + .attr('fill', 'black') + .attr('text-anchor', 'middle') + .text('Probability of Infection') + .attr('x', (graph_width + margins.right) / 2) + .attr('y', graph_height * 0.97); + + // set the parameters for the histogram + var histogram = d3.histogram() + .value(d => d) + .domain(x.domain()) // then the domain of the graphic + .thresholds(x.ticks(100)); // then the numbers of bins + + // And apply this function to data to get the bins + var bins = histogram(prob_dist); + + // Y left axis: scale and draw: + var y_left = d3.scaleLinear() + .range([graph_height - margins.bottom, margins.top]); + y_left.domain([0, d3.max(hist_count)]); // d3.hist has to be called before the Y axis obviously + vis.append("svg:g") + .attr('transform', 'translate(' + margins.left + ',0)') + .call(d3.axisLeft(y_left)); + + // Y left axis label. + vis.append('svg:text') + .attr('class', 'y label') + .attr('fill', 'black') + .attr('text-anchor', 'middle') + .text('Density') + .attr('x', (graph_height * 0.9 + margins.bottom) / 2) + .attr('y', (graph_height + margins.left) * 0.9) + .attr('transform', 'rotate(-90, 0,' + graph_height + ')'); + + // append the bar rectangles to the svg element + vis.selectAll("rect") + .data(bins.slice(0, -1)) + .enter() + .append("rect") + .attr("x", 1) + .attr("transform", function(d, i) { + return "translate(" + x(d.x0) + "," + y_left(hist_count[i]) + ")"; }) + .attr("width", function(d) { return x(d.x1) - x(d.x0) -1 ; }) + .attr("height", function(d, i) { return graph_height - y_left(hist_count[i]) - margins.bottom; }) + .attr('fill', '#1f77b4'); + + // Y right axis: scale and draw: + var y_right = d3.scaleLinear() + .range([graph_height - margins.bottom, margins.top]); + y_right.domain([0, 1]); + vis.append("svg:g") + .attr('transform', 'translate(' + (graph_width - margins.right) + ',0)') + .call(d3.axisRight(y_right)); + + // Y right axis label. + vis.append('svg:text') + .attr('class', 'y label') + .attr('fill', 'black') + .attr('text-anchor', 'middle') + .text('Cumulative Density Function (CDF)') + .attr('transform', 'rotate(-90, 0,' + graph_height + ')') + .attr('x', (graph_height + margins.bottom * 0.55) / 2) + .attr('y', graph_width + 430); + + // CDF Calculation + let count_sum = hist_count.reduce((partialSum, a) => partialSum + a, 0); + let pdf = hist_count.map((el, i) => el/count_sum); + let cdf = pdf.map((sum => value => sum += value)(0)); + // Add the CDF line + vis.append("svg:path") + .datum(cdf) + .attr("fill", "none") + .attr("stroke", "lightblue") + .attr("stroke-width", 1.5) + .attr("d", d3.line() + .x(function(d, i) { return x(hist_bins[i]) }) + .y(function(d) { return y_right(d) }) + ); + + // Add the mean dashed line + vis.append("svg:line") + .attr("fill", "none") + .attr('stroke-width', 2) + .attr('stroke-dasharray', (5, 5)) + .attr("x1", x(prob)) + .attr("y1", y_right(1)) + .attr("x2", x(prob)) + .attr("y2", y_right(0)) + .attr("stroke", "grey"); + + // Plot tile + vis.append("svg:text") + .attr("x", x(0.5)) + .attr("y", 0 + margins.top) + .attr("text-anchor", "middle") + .style("font-size", "16px") + .text(`P(I) -- Mean(SD) = ${prob.toFixed(2)}(${prob_sd.toFixed(2)}) `); + + // Legend for the plot elements + const size = 15; + var legend_x_start = 50; + const space_between_text_icon = 30; + const text_height = 6; + // CDF line icon + vis.append('rect') + .attr('width', 20) + .attr('height', 3) + .style('fill', 'lightblue') + .attr('x', graph_width + legend_x_start) + .attr('y', margins.top + size); + // CDF line text + vis.append('text') + .text('CDF') + .style('font-size', '15px') + .attr('x', graph_width + legend_x_start + space_between_text_icon) + .attr('y', margins.top + size + text_height); + // Hist icon + vis.append('rect') + .attr('width', 20) + .attr('height', 15) + .attr('fill', '#1f77b4') + .attr('x', graph_width + legend_x_start) + .attr('y', margins.top + (2 * size)); + // Hist text + vis.append('text') + .text('Histogram') + .style('font-size', '15px') + .attr('x', graph_width + legend_x_start + space_between_text_icon) + .attr('y', margins.top + 2 * size + text_height*2); + // Mean text + vis.append('line') + .attr('stroke', 'grey') + .attr('stroke-width', 2) + .attr('stroke-dasharray', (5, 5)) + .attr("x1", graph_width + legend_x_start) + .attr("x2", graph_width + legend_x_start + 20) + .attr("y1", margins.top + 3.85 * size) + .attr("y2", margins.top + 3.85 * size); + // Mean line text + vis.append('text') + .text('Mean') + .style('font-size', '15px') + .attr('x', graph_width + legend_x_start + space_between_text_icon) + .attr('y', margins.top + 3 * size + text_height*3); + + // Legend Bbox + vis.append('rect') + .attr('width', 120) + .attr('height', 65) + .attr('stroke', 'lightgrey') + .attr('stroke-width', '2') + .attr('rx', '5px') + .attr('ry', '5px') + .attr('stroke-linejoin', 'round') + .attr('fill', 'none') + .attr('x', graph_width * 1.07) + .attr('y', margins.top * 1.1); +} + function copy_clipboard(shareable_link) { $("#mobile_link").attr('title', 'Copied!') @@ -880,6 +1074,20 @@ function display_rename_column(bool, id) { else document.getElementById(id).style.display = 'none'; } +function conditional_probability_plot(value, is_generated) { + // If the image was previously generated, there is no need to reload the page. + if (value && is_generated == 1) { + document.getElementById('conditional_probability_div').style.display = 'block' + } + else if (value && is_generated == 0) { + document.getElementById('label_conditional_probability_plot').innerHTML = `Loading...`; + document.getElementById('conditional_probability_plot').setAttribute('disabled', true); + document.cookie = `conditional_plot= 1; path=/`; + window.location.reload(); + } + else document.getElementById('conditional_probability_div').style.display = "none"; +} + function export_csv() { // This function generates a CSV file according to the user's input. // It is composed of a list of lists. diff --git a/caimira/apps/templates/base/calculator.form.html.j2 b/caimira/apps/templates/base/calculator.form.html.j2 index 448350b7..8277f829 100644 --- a/caimira/apps/templates/base/calculator.form.html.j2 +++ b/caimira/apps/templates/base/calculator.form.html.j2 @@ -406,6 +406,10 @@ +
+ +
+
diff --git a/caimira/apps/templates/base/calculator.report.html.j2 b/caimira/apps/templates/base/calculator.report.html.j2 index 7c202f6f..cba60719 100644 --- a/caimira/apps/templates/base/calculator.report.html.j2 +++ b/caimira/apps/templates/base/calculator.report.html.j2 @@ -15,7 +15,7 @@ - +