diff --git a/cara/apps/calculator/__init__.py b/cara/apps/calculator/__init__.py index 1688d494..a60d5f74 100644 --- a/cara/apps/calculator/__init__.py +++ b/cara/apps/calculator/__init__.py @@ -33,7 +33,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 CARA version (found at ``cara.__version__``). -__version__ = "3.0.0" +__version__ = "3.0.1" class BaseRequestHandler(RequestHandler): diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index d1bff9dc..2f55e3eb 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -1,19 +1,15 @@ import concurrent.futures import base64 import dataclasses -from datetime import datetime, timedelta +from datetime import datetime import io +import json import typing import urllib import zlib -import loky import jinja2 -import matplotlib -matplotlib.use('agg') -import matplotlib.pyplot as plt import numpy as np -import qrcode from cara import models from ... import monte_carlo as mc @@ -125,7 +121,7 @@ def calculate_report_data(model: models.ExposureModel): } -def generate_qr_code(base_url, calculator_prefix, form: FormData): +def generate_permalink(base_url, calculator_prefix, form: FormData): form_dict = FormData.to_dict(form, strip_defaults=True) # Generate the calculator URL arguments that would be needed to re-create this @@ -138,20 +134,9 @@ def generate_qr_code(base_url, calculator_prefix, form: FormData): qr_url = f"{base_url}/_c/{compressed_args}" url = f"{base_url}{calculator_prefix}?{args}" - qr = qrcode.QRCode( - version=1, - error_correction=qrcode.constants.ERROR_CORRECT_H, - box_size=10, - border=4, - ) - qr.add_data(qr_url) - qr.make(fit=True) - img = qr.make_image(fill_color="black", back_color="white").convert('RGB') - return { - 'image': img2base64(_img2bytes(img)), 'link': url, - 'qr_url': qr_url, + 'shortened': qr_url, } @@ -162,50 +147,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) - 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 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='Mean concentration') - ax.spines['right'].set_visible(False) - ax.spines['top'].set_visible(False) - - ax.set_xlabel('Time of day') - ax.set_ylabel('Mean concentration ($virions/m^{3}$)') - ax.set_title('Mean concentration of virions') - ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) - - # Plot presence of exposed person - for i, (presence_start, presence_finish) in enumerate(model.exposed.presence.boundaries()): - plt.fill_between( - datetimes, concentrations, 0, - where=(np.array(times) > presence_start) & (np.array(times) < presence_finish), - color="#1f77b4", alpha=0.1, - label="Presence of exposed person(s)" if i == 0 else "" - ) - - # Place a legend outside of the axes itself. - ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') - ax.set_ylim(0) - - return fig - - def minutes_to_time(minutes: int) -> str: minute_string = str(minutes % 60) minute_string = "0" * (2 - len(minute_string)) + minute_string @@ -281,39 +229,7 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp 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) - - dash_styled_scenarios = [ - 'Base scenario with FFP2 masks', - 'Base scenario with HEPA filter', - 'Base scenario with HEPA and FFP2 masks', - ] - - sample_dts = [datetime(1970, 1, 1) + timedelta(hours=time) for time in sample_times] - for name, statistics in scenarios.items(): - concentrations = statistics['concentrations'] - - if name in dash_styled_scenarios: - ax.plot(sample_dts, concentrations, label=name, linestyle='--') - else: - ax.plot(sample_dts, concentrations, label=name, linestyle='-', alpha=0.5) - - # Place a legend outside of the axes itself. - ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') - ax.spines['right'].set_visible(False) - ax.spines['top'].set_visible(False) - 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 { 'probability_of_infection': np.mean(model.infection_probability()), @@ -342,7 +258,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, } @@ -385,7 +300,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['permalink'] = generate_permalink(base_url, self.calculator_prefix, form) context['calculator_prefix'] = self.calculator_prefix context['scale_warning'] = { 'level': 'yellow-2', @@ -405,6 +320,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/css/report.css b/cara/apps/calculator/static/css/report.css index 0e0a7cf7..c9b2b3ab 100644 --- a/cara/apps/calculator/static/css/report.css +++ b/cara/apps/calculator/static/css/report.css @@ -56,7 +56,7 @@ p.notes { margin: 1% } -#pdf-qr-code { +#pdf_qrcode_aref { margin-right: 1%; width: 100pt; } @@ -102,11 +102,6 @@ p.notes { border-radius: 5px; } -.print-button { - margin-left: auto; - margin-right: 1%; -} - /* @media (width: 1200px) { */ @media print { /* #body { @@ -120,7 +115,7 @@ p.notes { #link_reproduce_results { display: none!important; } - #pdf-qr-code { + #pdf_qrcode_aref { visibility: inherit!important; } .collapse { @@ -138,9 +133,6 @@ p.notes { .icon_button { display: none!important; } - .print-button { - display: none!important; - } .card { page-break-inside: avoid; } diff --git a/cara/apps/calculator/static/js/pdf.js b/cara/apps/calculator/static/js/pdf.js deleted file mode 100644 index 6afb4fea..00000000 --- a/cara/apps/calculator/static/js/pdf.js +++ /dev/null @@ -1,29 +0,0 @@ -function generate_pdf_version(qr_link) { - const pdf_version = this.document.getElementById("body"); - - // PDF styling - var opt = { - filename: 'myfile.pdf', - image: { type: 'jpeg', quality: 0.98 }, - html2canvas: { scale: 2, width: 1200, windowWidth: 1200 }, - enableLinks: false, - jsPDF: { - unit: 'pt', - format: 'letter', - orientation: 'portrait', - }, - pagebreak: { mode: '', avoid: '.break-avoid' }, - }; - html2pdf().set(opt).from(pdf_version).toPdf().get('pdf').then(function(pdf) { - var totalPages = pdf.internal.getNumberOfPages(); - pdf.setPage(1); - pdf.link(530, 25, 60, 60, { url: qr_link }); //Hyperlink to reproduce results - - for (i = 1; i <= totalPages; i++) { - pdf.setPage(i); - pdf.setFontSize(10); - pdf.setTextColor(150); - pdf.text('Page ' + i + ' of ' + totalPages, (pdf.internal.pageSize.getWidth() / 2.25), (pdf.internal.pageSize.getHeight() - 10)); - } - }).save(); -}; \ No newline at end of file diff --git a/cara/apps/calculator/static/js/report.js b/cara/apps/calculator/static/js/report.js index d1105bec..b815da2f 100644 --- a/cara/apps/calculator/static/js/report.js +++ b/cara/apps/calculator/static/js/report.js @@ -25,56 +25,16 @@ function draw_concentration_plot(svg_id, times, concentrations, exposed_presence yAxis = d3.axisLeft(yRange); // 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); + plot_scenario_data(vis, data, xTimeRange, yRange, '#1f77b4'); - vis.append('svg:path') - .attr('d', lineFunc(data)) - .attr('stroke', '#1f77b4') - .attr('stroke-width', 2) - .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 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³)'); + // Y axis + plot_y_axis(vis, height, width, margins, yAxis, 'Mean concentration (virions/m³)') // Area representing the presence of exposed person(s). exposed_presence_intervals.forEach(b => { @@ -181,4 +141,178 @@ function draw_concentration_plot(svg_id, times, concentrations, exposed_presence 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; } \ 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 ad04ffae..d06a21bc 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -23,9 +23,8 @@

CARA - CALCULATOR REPORT

Created {{ creation_date }} using CARA calculator version v{{ form.calculator_version }}

- - {# To be replaced by "Generate PDF" #} - + + {% endblock report_header %} @@ -88,9 +87,9 @@

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

@@ -108,8 +107,12 @@
- - + + {% block report_scenarios_summary_table %} @@ -158,10 +161,10 @@
- +

- Click the QR code to regenerate the report and get a shareable link.
Alternatively, scan to regenerate the report.
Mobile-friendly app coming soon! + Click the QR code to regenerate the report and get a shareable link.
Alternatively, scan to regenerate the report.

@@ -429,11 +432,26 @@
{% endblock disclaimer_container %} - + + + diff --git a/cara/apps/calculator/templates/calculator.form.html.j2 b/cara/apps/calculator/templates/calculator.form.html.j2 index b9dcd305..6626dd25 100644 --- a/cara/apps/calculator/templates/calculator.form.html.j2 +++ b/cara/apps/calculator/templates/calculator.form.html.j2 @@ -57,21 +57,16 @@ v{{ calculator_version }} Please sen
? -
- -
-
+

-
- -
- -
+
+ +

@@ -456,13 +451,13 @@ v{{ calculator_version }} Please sen Quick Guide:
This tool simulates the long range airborne spread SARS-CoV-2 virus in a finite volume and estimates the risk of COVID-19 infection. It is based on current scientific data and can be used to compare the effectiveness of different mitigation measures.
Virus data:
- SARS-CoV-2 covers typical strains of the virus and three variants of concern (VOC):
+ SARS-CoV-2 covers the original "wild type" strain of the virus and three variants of concern (VOC):
  • Alpha (also known as B.1.1.7, first identified in UK, Dec 2020),
  • Gamma (also known as P.1, first identified in Brazil/Japan, Jan 2021).
  • Delta (also known as B.1.617.2, first identified in India, Oct 2020).
- Choose variant according to local area prevalence, e.g. for Geneva + Modify the default as necessary, according to local area prevalence e.g. for Geneva or Ain (France).
Ventilation data:
    diff --git a/cara/models.py b/cara/models.py index 53cf52bc..b7885768 100644 --- a/cara/models.py +++ b/cara/models.py @@ -620,7 +620,6 @@ _ExpirationBase.types = { 'Talking': Expiration((1., 1., 1.)), 'Shouting': Expiration((1., 5., 5.)), 'Singing': Expiration((1., 5., 5.)), - 'Superspreading event': Expiration((np.inf, 0., 0.)), } @@ -669,44 +668,21 @@ class Population: @dataclass(frozen=True) -class InfectedPopulation(Population): +class _PopulationWithVirus(Population): #: The virus with which the population is infected. virus: Virus - #: The type of expiration that is being emitted whilst doing the activity. - expiration: _ExpirationBase - @method_cache def emission_rate_when_present(self) -> _VectorisedFloat: """ - The emission rate if the infected population is present. - - Note that the rate is not currently time-dependent. - + The emission rate if the infected population is present + (in virions / h). It should not be a function of time. """ - # Emission Rate (virions / h) - # Note on units: exhalation rate is in m^3/h, aerosols in mL/cm^3 - # and viral load in virus/mL -> 1e6 conversion factor - aerosols = self.expiration.aerosols(self.mask) + raise NotImplementedError("Subclass must implement") - ER = (self.virus.viral_load_in_sputum * - self.activity.exhalation_rate * - 10 ** 6 * - aerosols) - - # For superspreading event, where ejection_factor is infinite we fix the ER - # based on Miller et al. (2020). - if isinstance(aerosols, np.ndarray): - ER[np.isinf(aerosols)] = 970 * self.virus.infectious_dose - elif np.isinf(aerosols): - ER = 970 * self.virus.infectious_dose - - return ER - - def individual_emission_rate(self, time) -> _VectorisedFloat: + def emission_rate(self, time) -> _VectorisedFloat: """ - The emission rate of a single individual in the population. - + The emission rate of the population vs time. """ # Note: The original model avoids time dependence on the emission rate # at the cost of implementing a piecewise (on time) concentration function. @@ -721,19 +697,50 @@ class InfectedPopulation(Population): return self.emission_rate_when_present() - def emission_rate(self, time) -> _VectorisedFloat: - """ - The emission rate of the entire population. +@dataclass(frozen=True) +class EmittingPopulation(_PopulationWithVirus): + #: The emission rate of a single individual, in virions / h. + known_individual_emission_rate: float + + @method_cache + def emission_rate_when_present(self) -> _VectorisedFloat: """ - return self.individual_emission_rate(time) * self.number + The emission rate if the infected population is present. + """ + return self.known_individual_emission_rate * self.number + + +@dataclass(frozen=True) +class InfectedPopulation(_PopulationWithVirus): + #: The type of expiration that is being emitted whilst doing the activity. + expiration: _ExpirationBase + + @method_cache + def emission_rate_when_present(self) -> _VectorisedFloat: + """ + The emission rate if the infected population is present. + Note that the rate is not currently time-dependent. + """ + # Emission Rate (virions / h) + # Note on units: exhalation rate is in m^3/h, aerosols in mL/cm^3 + # and viral load in virus/mL -> 1e6 conversion factor + + aerosols = self.expiration.aerosols(self.mask) + + ER = (self.virus.viral_load_in_sputum * + self.activity.exhalation_rate * + 10 ** 6 * + aerosols) + + return ER * self.number @dataclass(frozen=True) class ConcentrationModel: room: Room ventilation: _VentilationBase - infected: InfectedPopulation + infected: _PopulationWithVirus @property def virus(self): @@ -753,16 +760,22 @@ class ConcentrationModel: ) @method_cache - def _concentration_limit(self, time: float) -> _VectorisedFloat: + def _normed_concentration_limit(self, time: float) -> _VectorisedFloat: """ Provides a constant that represents the theoretical asymptotic value reached by the concentration when time goes to infinity, if all parameters were to stay time-independent. + This is normalized by the emission rate, the latter acting as a + multiplicative constant factor for the concentration model that + can be put back in front of the concentration after the time + dependence has been solved for. """ + if not self.infected.person_present(time): + return 0. V = self.room.volume IVRR = self.infectious_virus_removal_rate(time) - return (self.infected.emission_rate(time)) / (IVRR * V) + return 1. / (IVRR * V) @method_cache def state_change_times(self) -> typing.List[float]: @@ -776,6 +789,14 @@ class ConcentrationModel: state_change_times.update(self.ventilation.transition_times()) return sorted(state_change_times) + @method_cache + def _first_presence_time(self) -> float: + """ + First presence time. Before that, the concentration is zero. + + """ + return self.infected.presence.boundaries()[0][0] + def last_state_change(self, time: float) -> float: """ Find the most recent/previous state change. @@ -807,15 +828,16 @@ class ConcentrationModel: ) @method_cache - def _concentration_cached(self, time: float) -> _VectorisedFloat: - # A cached version of the concentration method. Use this method if you - # expect that there may be multiple concentration calculations for the - # same time (e.g. at state change times). - return self.concentration(time) + def _normed_concentration_cached(self, time: float) -> _VectorisedFloat: + # A cached version of the _normed_concentration method. Use this + # method if you expect that there may be multiple concentration + # calculations for the same time (e.g. at state change times). + return self._normed_concentration(time) - def concentration(self, time: float) -> _VectorisedFloat: + def _normed_concentration(self, time: float) -> _VectorisedFloat: """ - Virus exposure concentration, as a function of time. + Virus exposure concentration, as a function of time, and + normalized by the emission rate. The formulas used here assume that all parameters (ventilation, emission rate) are constant between two state changes - only the value of these parameters at the next state change, are used. @@ -823,27 +845,42 @@ class ConcentrationModel: Note that time is not vectorised. You can only pass a single float to this method. """ - if time == 0: + # The model always starts at t=0, but we avoid running concentration calculations + # before the first presence as an optimisation. + if time <= self._first_presence_time(): return 0.0 next_state_change_time = self._next_state_change(time) IVRR = self.infectious_virus_removal_rate(next_state_change_time) - concentration_limit = self._concentration_limit(next_state_change_time) + conc_limit = self._normed_concentration_limit(next_state_change_time) t_last_state_change = self.last_state_change(time) - concentration_at_last_state_change = self._concentration_cached(t_last_state_change) + conc_at_last_state_change = self._normed_concentration_cached(t_last_state_change) delta_time = time - t_last_state_change fac = np.exp(-IVRR * delta_time) - return concentration_limit * (1 - fac) + concentration_at_last_state_change * fac + return conc_limit * (1 - fac) + conc_at_last_state_change * fac + + def concentration(self, time: float) -> _VectorisedFloat: + """ + Virus exposure concentration, as a function of time. + + Note that time is not vectorised. You can only pass a single float + to this method. + """ + return (self._normed_concentration(time) * + self.infected.emission_rate_when_present()) @method_cache - def integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat: + def normed_integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat: """ - Get the integrated concentration dose between the times start and stop. + Get the integrated concentration of viruses in the air between the times start and stop, + normalized by the emission rate. """ + if stop <= self._first_presence_time(): + return 0.0 state_change_times = self.state_change_times() req_start, req_stop = start, stop - total_concentration = 0. + total_normed_concentration = 0. for interval_start, interval_stop in zip(state_change_times[:-1], state_change_times[1:]): if req_start > interval_stop or req_stop < interval_start: continue @@ -851,17 +888,24 @@ class ConcentrationModel: start = max([interval_start, req_start]) stop = min([interval_stop, req_stop]) - conc_start = self._concentration_cached(start) + conc_start = self._normed_concentration_cached(start) next_conc_state = self._next_state_change(stop) - conc_limit = self._concentration_limit(next_conc_state) + conc_limit = self._normed_concentration_limit(next_conc_state) IVRR = self.infectious_virus_removal_rate(next_conc_state) delta_time = stop - start - total_concentration += ( + total_normed_concentration += ( conc_limit * delta_time + (conc_limit - conc_start) * (np.exp(-IVRR*delta_time)-1) / IVRR ) - return total_concentration + return total_normed_concentration + + def integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat: + """ + Get the integrated concentration of viruses in the air between the times start and stop. + """ + return (self.normed_integrated_concentration(start, stop) * + self.infected.emission_rate_when_present()) @dataclass(frozen=True) @@ -878,14 +922,22 @@ class ExposureModel: #: The fraction of viruses actually deposited in the respiratory tract fraction_deposited: _VectorisedFloat = 0.6 - def exposure(self) -> _VectorisedFloat: - """The number of virus per meter^3.""" - exposure = 0.0 + def _normed_exposure(self) -> _VectorisedFloat: + """ + The number of virus per meter^3, normalized by the emission rate + of the infected population. + """ + normed_exposure = 0.0 for start, stop in self.exposed.presence.boundaries(): - exposure += self.concentration_model.integrated_concentration(start, stop) + normed_exposure += self.concentration_model.normed_integrated_concentration(start, stop) - return exposure * self.repeats + return normed_exposure * self.repeats + + def exposure(self) -> _VectorisedFloat: + """The number of virus per meter^3.""" + return (self._normed_exposure() * + self.concentration_model.infected.emission_rate_when_present()) def infection_probability(self) -> _VectorisedFloat: exposure = self.exposure() diff --git a/cara/tests/apps/calculator/test_webapp.py b/cara/tests/apps/calculator/test_webapp.py index 40b75595..09bc104c 100644 --- a/cara/tests/apps/calculator/test_webapp.py +++ b/cara/tests/apps/calculator/test_webapp.py @@ -4,7 +4,7 @@ import pytest import tornado.testing import cara.apps.calculator -from cara.apps.calculator.report_generator import generate_qr_code +from cara.apps.calculator.report_generator import generate_permalink _TIMEOUT = 20. @@ -97,26 +97,26 @@ class TestOpenApp(tornado.testing.AsyncHTTPTestCase): assert response.code == 404 -async def test_qrcode_urls(http_server_client, baseline_form): +async def test_permalink_urls(http_server_client, baseline_form): base_url = 'proto://hostname/prefix' - qr_data = generate_qr_code(base_url, "/calculator", baseline_form) + permalink_data = generate_permalink(base_url, "/calculator", baseline_form) expected = f'{base_url}/calculator?exposed_coffee_break_option={baseline_form.exposed_coffee_break_option}&' - assert qr_data['link'].startswith(expected) + assert permalink_data['link'].startswith(expected) # We should get a 200 for the link. - response = await http_server_client.fetch(qr_data['link'].replace(base_url, '')) + response = await http_server_client.fetch(permalink_data['link'].replace(base_url, '')) assert response.code == 200 # And a 302 for the QR url itself. The redirected URL should be the same as # in the link. - assert qr_data['qr_url'].startswith(base_url) + assert permalink_data['shortened'].startswith(base_url) response = await http_server_client.fetch( - qr_data['qr_url'].replace(base_url, ''), + permalink_data['shortened'].replace(base_url, ''), max_redirects=0, raise_error=False, ) assert response.code == 302 - assert response.headers['Location'] == qr_data['link'].replace(base_url, '') + assert response.headers['Location'] == permalink_data['link'].replace(base_url, '') async def test_invalid_compressed_url(http_server_client, baseline_form): diff --git a/cara/tests/conftest.py b/cara/tests/conftest.py index 9499dfcf..e828e0cd 100644 --- a/cara/tests/conftest.py +++ b/cara/tests/conftest.py @@ -13,13 +13,15 @@ def baseline_model(): active=models.SpecificInterval(((0., 24.), )), air_exch=30., ), - infected=models.InfectedPopulation( + infected=models.EmittingPopulation( number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((0., 4.), (5., 8.))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light activity'], - expiration=models.Expiration.types['Superspreading event'], + known_individual_emission_rate=970 * 50, + # superspreading event, where ejection factor is fixed based + # on Miller et al. (2020) - 50 represents the infectious dose. ), ) return model diff --git a/cara/tests/models/test_concentration_model.py b/cara/tests/models/test_concentration_model.py index 46e7e88e..9373d4a9 100644 --- a/cara/tests/models/test_concentration_model.py +++ b/cara/tests/models/test_concentration_model.py @@ -121,6 +121,10 @@ def test_next_state_change_time_out_of_range(simple_conc_model: models.Concentra simple_conc_model._next_state_change(3.1) +def test_first_presence_time(simple_conc_model): + assert simple_conc_model._first_presence_time() == 0.5 + + def test_integrated_concentration(simple_conc_model): c1 = simple_conc_model.integrated_concentration(0, 2) c2 = simple_conc_model.integrated_concentration(0, 1) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index cf8cfb6e..ac82dae3 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -11,20 +11,20 @@ from cara.dataclass_utils import replace @dataclass(frozen=True) -class KnownConcentrations(models.ConcentrationModel): +class KnownNormedconcentration(models.ConcentrationModel): """ A ConcentrationModel which is based on pre-known exposure concentrations and which therefore doesn't need other components. Useful for testing. """ - concentration_function: typing.Callable + normed_concentration_function: typing.Callable def infectious_virus_removal_rate(self, time: float) -> models._VectorisedFloat: # very large decay constant -> same as constant concentration return 1.e50 - def _concentration_limit(self, time: float) -> models._VectorisedFloat: - return self.concentration_function(time) + def _normed_concentration_limit(self, time: float) -> models._VectorisedFloat: + return self.normed_concentration_function(time) def state_change_times(self): return [0., 24.] @@ -32,8 +32,8 @@ class KnownConcentrations(models.ConcentrationModel): def _next_state_change(self, time: float): return 24. - def concentration(self, time: float) -> models._VectorisedFloat: # noqa - return self.concentration_function(time) + def _normed_concentration(self, time: float) -> models._VectorisedFloat: # noqa + return self.normed_concentration_function(time) halftime = models.PeriodicInterval(120, 60) @@ -67,7 +67,9 @@ def known_concentrations(func): virus=models.Virus.types['SARS_CoV_2_B117'], expiration=models.Expiration.types['Talking'] ) - return KnownConcentrations(dummy_room, dummy_ventilation, dummy_infected_population, func) + normed_func = lambda x: func(x) / dummy_infected_population.emission_rate_when_present() + return KnownNormedconcentration(dummy_room, dummy_ventilation, + dummy_infected_population, normed_func) @pytest.mark.parametrize( @@ -142,13 +144,15 @@ def conc_model(): return models.ConcentrationModel( models.Room(25), models.AirChange(always, 5), - models.InfectedPopulation( + models.EmittingPopulation( number=1, presence=interesting_times, mask=models.Mask.types['No mask'], activity=models.Activity.types['Seated'], virus=models.Virus.types['SARS_CoV_2'], - expiration=models.Expiration.types['Superspreading event'], + known_individual_emission_rate=970 * 50, + # superspreading event, where ejection factor is fixed based + # on Miller et al. (2020) - 50 represents the infectious dose. ) ) diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 89c0268d..d1db42c2 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -66,13 +66,15 @@ def build_model(interval_duration): active=models.PeriodicInterval(period=120, duration=interval_duration), q_air_mech=500., ), - infected=models.InfectedPopulation( + infected=models.EmittingPopulation( number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((0., 4.), (5., 8.))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light activity'], - expiration=models.Expiration.types['Superspreading event'], + known_individual_emission_rate=970 * 50, + # superspreading event, where ejection factor is fixed based + # on Miller et al. (2020) - 50 represents the infectious dose. ), ) return model @@ -226,13 +228,13 @@ def build_hourly_dependent_model( outside_temp=outside_temp, window_height=1.6, opening_length=0.6, ), - infected=models.InfectedPopulation( + infected=models.EmittingPopulation( number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(intervals_presence_infected), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light activity'], - expiration=models.Expiration.types['Superspreading event'], + known_individual_emission_rate=970 * 50, ), ) return model @@ -247,13 +249,13 @@ def build_constant_temp_model(outside_temp, intervals_open=((7.5, 8.5),)): outside_temp=models.PiecewiseConstant((0., 24.), (outside_temp,)), window_height=1.6, opening_length=0.6, ), - infected=models.InfectedPopulation( + infected=models.EmittingPopulation( number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((0., 4.), (5., 7.5))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light activity'], - expiration=models.Expiration.types['Superspreading event'], + known_individual_emission_rate=970 * 50, ), ) return model @@ -275,13 +277,13 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5 model = models.ConcentrationModel( room=models.Room(volume=75), ventilation=vent, - infected=models.InfectedPopulation( + infected=models.EmittingPopulation( number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((0., 4.), (5., 7.5))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light activity'], - expiration=models.Expiration.types['Superspreading event'], + known_individual_emission_rate=970 * 50, ), ) return model diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 5cf88198..08bd2dad 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -238,7 +238,7 @@ def skagit_chorale_mc(): ["shared_office_mc", 10.7, 0.32, 57.24, 654], ["classroom_mc", 36.1, 6.85, 780.0, 28464], ["ski_cabin_mc", 16.3, 0.49, 35.94, 7404], - ["gym_mc", 2.25, 0.63, 0.7842, 984], + ["gym_mc", 2.25, 0.63, 0.7842, 1968], ["waiting_room_mc", 9.72, 1.36, 34.26, 3534], ["skagit_chorale_mc",29.9, 17.9, 190.0, 141400], ] diff --git a/requirements.txt b/requirements.txt index 457a69b7..724bf5cc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -62,7 +62,6 @@ pyparsing==2.4.7 pyrsistent==0.18.0 python-dateutil==2.8.2 pyzmq==22.1.0 -qrcode==7.2 requests==2.26.0 requests-unixsocket==0.2.0 scikit-learn==0.24.2 diff --git a/setup.py b/setup.py index 9e5eb596..704f923c 100644 --- a/setup.py +++ b/setup.py @@ -20,7 +20,7 @@ REQUIREMENTS: dict = { 'core': [ 'dataclasses; python_version < "3.7"', 'ipykernel', - 'ipympl', + 'ipympl != 0.8.0', 'ipywidgets', 'Jinja2', 'loky', @@ -30,7 +30,6 @@ REQUIREMENTS: dict = { 'numpy', 'psutil', 'python-dateutil', - 'qrcode[pil]', 'scipy', 'sklearn', 'timezonefinder',