Merge branch 'feature/PoI_uncertainties' into 'master'

Prob. of infection uncertainties

See merge request caimira/caimira!436
This commit is contained in:
Andre Henriques 2023-04-03 16:41:15 +02:00
commit d30cfbdac9
6 changed files with 350 additions and 16 deletions

View file

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

View file

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

View file

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

View file

@ -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 = `<span id="loading_spinner" class="spinner-border spinner-border-sm mr-2 mt-0" role="status" aria-hidden="true"></span>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.

View file

@ -406,6 +406,10 @@
</div>
</div>
<div class="form-check d-none">
<input type="checkbox" id="conditional_probability_plot" class="tabbed form-check-input" name="conditional_probability_plot" value="0" disabled>
</div>
<span id="training_limit_error" class="red_text" hidden>Conference/Training activities limited to 1 infected<br></span>
<hr width="80%">

View file

@ -15,7 +15,7 @@
</head>
<body id="body">
<body id="body" onload="on_report_load({{ form.conditional_probability_plot | int }})">
<!-- MODEL REPR - Available in the developer tools once the report is generated. Useful to re-create the model using an interpreter that has CAiMIRA installed:
@ -171,17 +171,55 @@
{% endif %}
<div id="concentration_plot" style="height: 400px"></div>
<script type="application/javascript">
let times = {{ times | JSONify }}
let concentrations_zoomed = {{ concentrations_zoomed | JSONify }}
let concentrations = {{ concentrations | JSONify }}
let cumulative_doses = {{ cumulative_doses | JSONify }}
let long_range_cumulative_doses = {{ long_range_cumulative_doses | JSONify }}
let exposed_presence_intervals = {{ exposed_presence_intervals | JSONify }}
let short_range_intervals = {{ short_range_intervals | JSONify }}
let short_range_expirations = {{ short_range_expirations | JSONify }}
draw_plot("concentration_plot")
let times = {{ times | JSONify }};
let concentrations_zoomed = {{ concentrations_zoomed | JSONify }};
let concentrations = {{ concentrations | JSONify }};
let cumulative_doses = {{ cumulative_doses | JSONify }};
let long_range_cumulative_doses = {{ long_range_cumulative_doses | JSONify }};
let exposed_presence_intervals = {{ exposed_presence_intervals | JSONify }};
let short_range_intervals = {{ short_range_intervals | JSONify }};
let short_range_expirations = {{ short_range_expirations | JSONify }};
draw_plot("concentration_plot");
</script>
</p>
</p>
</div>
</div>
</div>
<div class="card bg-light mb-3" id="results-div">
<div class="card-header"><strong>Result uncertainties </strong>
<button class="icon_button p-0 float-right" data-toggle="collapse" href="#collapseUncertainties" role="button" aria-expanded="true" aria-controls="collapseUncertainties">
<svg xmlns="http://www.w3.org/2000/svg" width="16" height="16" fill="currentColor" class="bi bi-chevron-expand" viewBox="0 0 16 16">
<path fill-rule="evenodd" d="M3.646 9.146a.5.5 0 0 1 .708 0L8 12.793l3.646-3.647a.5.5 0 0 1 .708.708l-4 4a.5.5 0 0 1-.708 0l-4-4a.5.5 0 0 1 0-.708zm0-2.292a.5.5 0 0 0 .708 0L8 3.207l3.646 3.647a.5.5 0 0 0 .708-.708l-4-4a.5.5 0 0 0-.708 0l-4 4a.5.5 0 0 0 0 .708z"/>
</svg>
</button>
</div>
<div class="collapse show" id="collapseUncertainties">
<div class="card-body">
<div class="align-self-center">
<div id="prob_inf_hist" style="height: 400px"></div>
<script type="application/javascript">
let prob_dist = {{ prob_dist | JSONify }};
let prob_hist_count = {{ prob_hist_count | JSONify }};
let prob_hist_bins = {{ prob_hist_bins | JSONify }};
draw_histogram("prob_inf_hist", {{ prob_inf }}, {{ prob_inf_sd }});
</script>
<br>
<div class="form-check">
<input type="checkbox" id="conditional_probability_plot" class="tabbed form-check-input" name="conditional_probability_plot" value="1" onClick="conditional_probability_plot(this.checked, {{ form.conditional_probability_plot | int }});">
<label id="label_conditional_probability_plot" for="conditional_probability_plot" class="form-check-label col-sm-12">Generate full uncertainty data (as function of the viral load)</label>
</div>
{% if form.conditional_probability_plot %}
<div id="conditional_probability_div">
<img src= "{{ uncertainties_plot_src }}" />
<div class="ml-5">
<p>(i) &nbsp;&nbsp;Predictive probability of infection for a given value of the viral load</p>
<p>(ii) &nbsp;Histogram of the viral load data</p>
<p>(iii) Histogram of the conditional probability of infection (result of total predictive probability in the middle)</p>
</div>
</div>
{% endif %}
</div>
</div>
</div>
</div>
@ -206,7 +244,7 @@
{% endif %}
<div id="alternative_scenario_plot" style="height: 400px"></div>
<script type="application/javascript">
let alternative_scenarios = {{ alternative_scenarios.stats | JSONify }}
let alternative_scenarios = {{ alternative_scenarios.stats | JSONify }};
draw_alternative_scenarios_plot("concentration_plot", "alternative_scenario_plot");
</script>
<br>