diff --git a/README.md b/README.md index 02b95b78..90319c63 100644 --- a/README.md +++ b/README.md @@ -2,15 +2,15 @@ CARA is a risk assessment tool developed to model the concentration of viruses in enclosed spaces, in order to inform space-management decisions. -CARA models the concentration profile of potential virions in enclosed spaces with clear and intuitive graphs. +CARA models the concentration profile of potential virions in enclosed spaces , both as background (room) concentration and during close-proximity interations, with clear and intuitive graphs. The user can set a number of parameters, including room volume, exposure time, activity type, mask-wearing and ventilation. The report generated indicates how to avoid exceeding critical concentrations and chains of airborne transmission in spaces such as individual offices, meeting rooms and labs. -The risk assessment tool simulates the long-range airborne spread SARS-CoV-2 virus in a finite volume, assuming a homogenous mixture, and estimates the risk of COVID-19 infection therein. -The results DO NOT include short-range airborne exposure (where the physical distance is a significant factor) nor the other known modes of SARS-CoV-2 transmission. -Hence, the output from this model is only valid when the other recommended public health & safety instructions are observed, such as adequate physical distancing, good hand hygiene and other barrier measures. +The risk assessment tool simulates the airborne spread SARS-CoV-2 virus in a finite volume, assuming a homogenous mixture and a two-stage exhaled jet model, and estimates the risk of COVID-19 infection therein. +The results DO NOT include the other known modes of SARS-CoV-2 transmission, such as fomite or blood-bound. +Hence, the output from this model is only valid when the other recommended public health & safety instructions are observed, such as good hand hygiene and other barrier measures. -The model used is based on scientific publications relating to airborne transmission of infectious diseases, dose-response exposures and aerosol science, as of February 2021. +The model used is based on scientific publications relating to airborne transmission of infectious diseases, dose-response exposures and aerosol science, as of February 2022. It can be used to compare the effectiveness of different airborne-related risk mitigation measures. Note that this model applies a deterministic approach, i.e., it is assumed at least one person is infected and shedding viruses into the simulated volume. @@ -44,11 +44,13 @@ CARA – COVID Airborne Risk Assessment tool **For use of the model** Henriques A, Mounet N, Aleixo L, Elson P, Devine J, Azzopardi G, Andreini M, Rognlien M, Tarocco N, Tang J. (2022). Modelling airborne transmission of SARS-CoV-2 using CARA: risk assessment for enclosed spaces. _Interface Focus 20210076_. https://doi.org/10.1098/rsfs.2021.0076 +_Note that the short-range component of the model has not yet been published._ + ## Applications ### COVID Calculator -A risk assessment tool which simulates the long range airborne spread of the SARS-CoV-2 virus for space managers. +A risk assessment tool which simulates the long-range airborne spread of the SARS-CoV-2 virus for space managers. ### CARA Expert App diff --git a/cara/apps/calculator/__init__.py b/cara/apps/calculator/__init__.py index e396ea8c..c0c93887 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__ = "4.0.0" +__version__ = "4.1.0" class BaseRequestHandler(RequestHandler): diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index cccce020..ab90fab8 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -3,6 +3,8 @@ import datetime import html import logging import typing +import ast +import json import numpy as np @@ -11,8 +13,8 @@ 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, mask_distributions -from cara.monte_carlo.data import expiration_distribution, expiration_BLO_factors, expiration_distributions +from cara.monte_carlo.data import activity_distributions, virus_distributions, mask_distributions, short_range_distances +from cara.monte_carlo.data import expiration_distribution, expiration_BLO_factors, expiration_distributions, short_range_expiration_distributions LOG = logging.getLogger(__name__) @@ -76,6 +78,8 @@ class FormData: window_width: float windows_number: int window_opening_regime: str + short_range_option: str + short_range_interactions: list #: The default values for undefined fields. Note that the defaults here #: and the defaults in the html form must not be contradictory. @@ -127,6 +131,8 @@ class FormData: 'windows_frequency': 0., 'windows_number': 0, 'window_opening_regime': 'windows_open_permanently', + 'short_range_option': 'short_range_no', + 'short_range_interactions': '[]', } @classmethod @@ -240,14 +246,27 @@ class FormData: humidity = 0.5 room = models.Room(volume=volume, humidity=humidity) + infected_population = self.infected_population() + + short_range = [] + if self.short_range_option == "short_range_yes": + for interaction in self.short_range_interactions: + short_range.append(mc.ShortRangeModel( + expiration=short_range_expiration_distributions[interaction['expiration']], + activity=infected_population.activity, + presence=self.short_range_interval(interaction), + distance=short_range_distances, + )) + # Initializes and returns a model with the attributes defined above return mc.ExposureModel( concentration_model=mc.ConcentrationModel( room=room, ventilation=self.ventilation(), - infected=self.infected_population(), + infected=infected_population, evaporation_factor=0.3, ), + short_range = tuple(short_range), exposed=self.exposed_population(), ) @@ -629,6 +648,11 @@ class FormData: breaks=self.infected_lunch_break_times() + self.infected_coffee_break_times(), ) + def short_range_interval(self, interaction) -> models.SpecificInterval: + start_time = time_string_to_minutes(interaction['start_time']) + duration = float(interaction['duration']) + return models.SpecificInterval(present_times=((start_time/60, (start_time + duration)/60),)) + def exposed_present_interval(self) -> models.Interval: return self.present_interval( self.exposed_start, self.exposed_finish, @@ -636,7 +660,7 @@ class FormData: ) -def build_expiration(expiration_definition) -> models._ExpirationBase: +def build_expiration(expiration_definition) -> mc._ExpirationBase: if isinstance(expiration_definition, str): return expiration_distributions[expiration_definition] elif isinstance(expiration_definition, dict): @@ -645,7 +669,7 @@ def build_expiration(expiration_definition) -> models._ExpirationBase: np.array(expiration_BLO_factors[exp_type]) * weight/total_weight for exp_type, weight in expiration_definition.items() ], axis=0) - return expiration_distribution(tuple(BLO_factors)) + return expiration_distribution(BLO_factors=tuple(BLO_factors)) def baseline_raw_form_data(): @@ -697,7 +721,9 @@ def baseline_raw_form_data(): 'window_type': 'window_sliding', 'window_width': '2', 'windows_number': '1', - 'window_opening_regime': 'windows_open_permanently' + 'window_opening_regime': 'windows_open_permanently', + 'short_range_option': 'short_range_no', + 'short_range_interactions': '[]', } @@ -742,6 +768,14 @@ def time_minutes_to_string(time: int) -> str: return "{0:0=2d}".format(int(time/60)) + ":" + "{0:0=2d}".format(time%60) +def string_to_list(l: str) -> list: + return list(ast.literal_eval(l.replace(""", "\""))) + + +def list_to_string(s: list) -> str: + return json.dumps(s) + + def _safe_int_cast(value) -> int: if isinstance(value, int): return value @@ -773,3 +807,6 @@ for _field in dataclasses.fields(FormData): elif _field.type is bool: _CAST_RULES_FORM_ARG_TO_NATIVE[_field.name] = lambda v: v == '1' _CAST_RULES_NATIVE_TO_FORM_ARG[_field.name] = int + elif _field.type is list: + _CAST_RULES_FORM_ARG_TO_NATIVE[_field.name] = string_to_list + _CAST_RULES_NATIVE_TO_FORM_ARG[_field.name] = list_to_string diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index c06b63ee..b526d5bf 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -97,29 +97,54 @@ def interesting_times(model: models.ExposureModel, approx_n_pts=100) -> typing.L return nice_times -def calculate_report_data(model: models.ExposureModel): - times = interesting_times(model) +def concentrations_with_sr_breathing(form: FormData, model: models.ExposureModel, times: typing.List[float], short_range_intervals: typing.List) -> typing.List[float]: + lower_concentrations = [] + for time in times: + for index, (start, stop) in enumerate(short_range_intervals): + # For visualization issues, add short-range breathing activity to the initial long-range concentrations + if start <= time <= stop and form.short_range_interactions[index]['expiration'] == 'Breathing': + lower_concentrations.append(np.array(model.concentration(float(time))).mean()) + break + lower_concentrations.append(np.array(model.concentration_model.concentration(float(time))).mean()) + return lower_concentrations + +def calculate_report_data(form: FormData, model: models.ExposureModel): + times = interesting_times(model) + short_range_intervals = [interaction.presence.boundaries()[0] for interaction in model.short_range] + short_range_expirations = [interaction['expiration'] for interaction in form.short_range_interactions] if form.short_range_option == "short_range_yes" else [] + concentrations = [ - np.array(model.concentration_model.concentration(float(time))).mean() + np.array(model.concentration(float(time))).mean() for time in times - ] + ] + lower_concentrations = concentrations_with_sr_breathing(form, model, times, short_range_intervals) highest_const = max(concentrations) - prob = np.array(model.infection_probability()).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 = np.cumsum([ np.array(model.deposited_exposure_between_bounds(float(time1), float(time2))).mean() for time1, time2 in zip(times[:-1], times[1:]) ]) + long_range_cumulative_doses = np.cumsum([ + np.array(model.long_range_deposited_exposure_between_bounds(float(time1), float(time2))).mean() + for time1, time2 in zip(times[:-1], times[1:]) + ]) + + prob = np.array(model.infection_probability()).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() return { "times": list(times), "exposed_presence_intervals": [list(interval) for interval in model.exposed.presence.boundaries()], - "cumulative_doses": list(cumulative_doses), + "short_range_intervals": short_range_intervals, + "short_range_expirations": short_range_expirations, "concentrations": concentrations, + "concentrations_zoomed": lower_concentrations, "highest_const": highest_const, + "cumulative_doses": list(cumulative_doses), + "long_range_cumulative_doses": list(long_range_cumulative_doses), "prob_inf": prob, "emission_rate": er, "exposed_occupants": exposed_occupants, @@ -197,51 +222,52 @@ def non_zero_percentage(percentage: int) -> str: def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.ExposureModel]: scenarios = {} + if (form.short_range_option == "short_range_no"): + # Two special option cases - HEPA and/or FFP2 masks. + 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') + 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) + scenarios['Base scenario without HEPA filter, with FFP2 masks'] = noHEPAalternative.build_mc_model() - # Two special option cases - HEPA and/or FFP2 masks. - 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') - 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) - 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) + # and no HEPA filtration. + form = dataclass_utils.replace(form, mask_type='Type I') + if form.hepa_option: + form = dataclass_utils.replace(form, hepa_option=False) - # The remaining scenarios are based on Type I masks (possibly not worn) - # and no HEPA filtration. - form = dataclass_utils.replace(form, mask_type='Type I') - if form.hepa_option: - 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') - with_mask = dataclass_utils.replace(form, mask_wearing_option='mask_on') - 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() + scenarios['Mechanical ventilation without masks'] = without_mask.build_mc_model() - if form.ventilation_type == 'mechanical_ventilation': - #scenarios['Mechanical ventilation with Type I masks'] = with_mask.build_mc_model() - scenarios['Mechanical ventilation without masks'] = without_mask.build_mc_model() + elif form.ventilation_type == 'natural_ventilation': + #scenarios['Windows open with Type I masks'] = with_mask.build_mc_model() + scenarios['Windows open without masks'] = without_mask.build_mc_model() - elif form.ventilation_type == 'natural_ventilation': - #scenarios['Windows open with Type I masks'] = with_mask.build_mc_model() - 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') - 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() + # 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') + 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 scenario_statistics(mc_model: mc.ExposureModel, sample_times: np.ndarray): +def scenario_statistics(mc_model: mc.ExposureModel, sample_times: typing.List[float]): model = mc_model.build_model(size=_DEFAULT_MC_SAMPLE_SIZE) + return { 'probability_of_infection': np.mean(model.infection_probability()), 'expected_new_cases': np.mean(model.expected_new_cases()), 'concentrations': [ - np.mean(model.concentration_model.concentration(time)) + np.mean(model.concentration(time)) for time in sample_times ], } @@ -301,7 +327,7 @@ class ReportGenerator: scenario_sample_times = interesting_times(model) - context.update(calculate_report_data(model)) + context.update(calculate_report_data(form, model)) alternative_scenarios = manufacture_alternative_scenarios(form) context['alternative_scenarios'] = comparison_report( alternative_scenarios, scenario_sample_times, executor_factory=executor_factory, diff --git a/cara/apps/calculator/static/css/report.css b/cara/apps/calculator/static/css/report.css index d2b096e3..3a5a5487 100644 --- a/cara/apps/calculator/static/css/report.css +++ b/cara/apps/calculator/static/css/report.css @@ -144,6 +144,15 @@ p.notes { padding: 15px; page-break-inside: avoid; } + #button_full_exposure, #button_hide_high_concentration { + display: none!important; + } + #long_range_cumulative_checkbox, #lr_cumulative_checkbox_label { + display: none!important; + } + #button_alternative_full_exposure, #button_alternative_hide_high_concentration { + display: none!important; + } } diff --git a/cara/apps/calculator/static/js/form.js b/cara/apps/calculator/static/js/form.js index 5a20f025..fd4f3c73 100644 --- a/cara/apps/calculator/static/js/form.js +++ b/cara/apps/calculator/static/js/form.js @@ -242,6 +242,16 @@ function on_wearing_mask_change() { if (this.checked) { getChildElement($(this)).show(); require_fields(this); + if (this.id == "mask_on") { + $('#short_range_no').click(); + $('input[name="short_range_option"]').attr('disabled', true); + $("#short_range_warning").show(); + } + else { + $('input[name="short_range_option"]').attr('disabled', false); + $("#short_range_warning").hide(); + } + } else { getChildElement($(this)).hide(); @@ -250,6 +260,20 @@ function on_wearing_mask_change() { }) } +function on_short_range_option_change() { + short_range = $('input[type=radio][name=short_range_option]') + short_range.each(function (index){ + if (this.checked) { + getChildElement($(this)).show(); + require_fields(this); + } + else { + getChildElement($(this)).hide(); + require_fields(this); + } + }) +} + /* -------UI------- */ function show_disclaimer() { @@ -379,6 +403,27 @@ function validate_form(form) { } } + // Generate the short-range interactions list + var short_range_interactions = []; + $(".form_field_outer_row").each(function (index, element){ + let obj = {}; + const $element = $(element); + obj.expiration = $element.find("[name='short_range_expiration']").val(); + obj.start_time = $element.find("[name='short_range_start_time']").val(); + obj.duration = $element.find("[name='short_range_duration']").val(); + short_range_interactions.push(JSON.stringify(obj)); + }); + + // Sort list by time + short_range_interactions.sort(function (a, b) { + return JSON.parse(a).start_time.localeCompare(JSON.parse(b).start_time); + }); + $("input[type=text][name=short_range_interactions]").val('[' + short_range_interactions + ']'); + if (short_range_interactions.length == 0) { + $("input[type=radio][id=short_range_no]").prop("checked", true); + on_short_range_option_change(); + } + if (submit) { $("#generate_report").prop("disabled", true); //Add spinner to button @@ -492,6 +537,84 @@ function validateLunchTime(obj) { return true; } +function overlapped_times(obj, start_time, finish_time) { + removeErrorFor($(obj)); + $(obj).removeClass("red_border"); + + let parameter = document.getElementById($(obj).attr('id')); + + if ($(obj).attr('name') == "short_range_duration" && parseFloat($(obj).val()) < 15.0) { + if (!$(obj).hasClass("red_border")) $(parameter).addClass("red_border"); //Adds the red border and error message. + insertErrorFor(parameter, "Must be ≥ 15 min.") + return false; + } + + let simulation_start = parseTimeToMins($("#exposed_start").val()) + let simulation_finish = parseTimeToMins($("#exposed_finish").val()) + var simulation_lunch_start, simulation_lunch_finish; + if ($('input[name=exposed_lunch_option]:checked').val() == 1) { + simulation_lunch_start = parseTimeToMins($("#exposed_lunch_start").val()) + simulation_lunch_finish = parseTimeToMins($("#exposed_lunch_finish").val()) + } else { + simulation_lunch_start = 0 + simulation_lunch_finish = 0 + } + if (start_time < simulation_start || start_time > simulation_finish || + finish_time < simulation_start || finish_time > simulation_finish || + start_time >= simulation_lunch_start && start_time <= simulation_lunch_finish || + finish_time >= simulation_lunch_start && finish_time <= simulation_lunch_finish ) {//If start and finish inputs are out of the simulation period + //Adds the red border and error message. + if (!$(obj).hasClass("red_border")) $(parameter).addClass("red_border"); + insertErrorFor(parameter, "Out of event time."); + return false; + } + let current_interaction = $(obj).closest(".form_field_outer_row"); + var toReturn = true; + $(".form_field_outer_row.row_validated").not(current_interaction).each(function(index, el) { + let current_start_el = $(el).find("input[name='short_range_start_time']"); + let current_duration_el = $(el).find("input[name='short_range_duration']") + start_time_2 = parseTimeToMins(current_start_el.val()) + finish_time_2 = parseTimeToMins(current_start_el.val()) + parseInt(current_duration_el.val()); + if ((start_time > start_time_2 && start_time < finish_time_2) || ( //If hour input is within other time range + finish_time > start_time_2 && finish_time < finish_time_2) || //If finish time input is within other time range + (start_time <= start_time_2 && finish_time >= finish_time_2) || //If start and finish inputs encompass other time range + start_time == start_time_2) { + if (!$(obj).hasClass("red_border")) $(parameter).addClass("red_border"); //Adds the red border and error message. + insertErrorFor(parameter, "Time overlap.") + toReturn = false; + return false; + } + }); + return toReturn; +} + +function validate_sr_time(obj) { + let obj_id = $(obj).attr('id').split('_').slice(-1)[0]; + var start_time, finish_time; + if ($(obj).val() != "") { + if ($('#sr_start_no_' + String(obj_id)).val()) start_time = parseTimeToMins($('#sr_start_no_' + String(obj_id)).val()); + else start = 0. + finish_time = start_time + parseInt($('#sr_duration_no_' + String(obj_id)).val()); + } + return overlapped_times(obj, start_time, finish_time); +}; + +// Check if short-range durations are filled, and if there is no repetitions +function validate_sr_parameter(obj, error_message) { + if ($(obj).val() == "" || $(obj).val() == null) { + if (!$(obj).hasClass("red_border") && !$(obj).prop("disabled")) { + var parameter = document.getElementById($(obj).attr('id')); + insertErrorFor(parameter, error_message) + $(parameter).addClass("red_border"); + } + return false; + } else { + removeErrorFor($(obj)); + $(obj).removeClass("red_border"); + return true; + } +} + function parseValToNumber(val) { return parseInt(val.replace(':',''), 10); } @@ -529,6 +652,22 @@ $(document).ready(function () { elemObj.checked = (value==1); } + // Read short-range from URL + else if (name == 'short_range_interactions') { + let index = 1; + for (const interaction of JSON.parse(value)) { + $("#dialog_sr").append(inject_sr_interaction(index, value = interaction, is_validated="row_validated")) + $('#sr_expiration_no_' + String(index)).val(interaction.expiration).change(); + document.getElementById('sr_expiration_no_' + String(index)).disabled = true; + document.getElementById('sr_start_no_' + String(index)).disabled = true; + document.getElementById('sr_duration_no_' + String(index)).disabled = true; + document.getElementById('edit_row_no_' + String(index)).style.cssText = 'display:inline !important'; + document.getElementById('validate_row_no_' + String(index)).style.cssText = 'display: none !important'; + index++; + } + $("#sr_interactions").text(index - 1); + } + //Ignore 0 (default) values from server side else if (!(elemObj.classList.contains("non_zero") || elemObj.classList.contains("remove_zero")) || (value != "0.0" && value != "0")) { elemObj.value = value; @@ -578,6 +717,13 @@ $(document).ready(function () { // Call the function now to handle forward/back button presses in the browser. on_wearing_mask_change(); + // When the short_range_option changes we want to make its respective + // children show/hide. + $("input[type=radio][name=short_range_option]").change(on_short_range_option_change); + + // Call the function now to handle forward/back button presses in the browser. + on_short_range_option_change(); + // Setup the maximum number of people at page load (to handle back/forward), // and update it when total people is changed. setMaxInfectedPeople(); @@ -690,6 +836,134 @@ $(document).ready(function () { } return selectedSuggestion.text; } + + function inject_sr_interaction(index, value, is_validated) { + return `
+
+ +
+

+

+
+
+ +
+
+

+
+ +
+
+

+
+ +
+ + + +
+
+
` + } + + // Add one empty row if none. + $("#set_interactions_button").on("click", e => { + if ($(".form_field_outer").find(".form_field_outer_row").length == 0) $(".add_node_btn_frm_field").click(); + }); + + // When short_range_yes option is selected, we want to inject rows for each expiractory activity, start_time and duration. + $("body").on("click", ".add_node_btn_frm_field", function(e) { + let last_row = $(".form_field_outer").find(".form_field_outer_row"); + if (last_row.length == 0) $("#dialog_sr").append(inject_sr_interaction(1, value = { activity: "", start_time: "", duration: "15" })); + else { + last_index = last_row.last().find(".short_range_option").prop("id").split("_").slice(-1)[0]; + index = parseInt(last_index) + 1; + $("#dialog_sr").append(inject_sr_interaction(index, value = { activity: "", start_time: "", duration: "15" })); + } + }); + + // Validate row button (Save button) + $("body").on("click", ".validate_node_btn_frm_field", function() { + var index = $(this).attr('id').split('_').slice(-1)[0]; + let activity = validate_sr_parameter('#sr_expiration_no_' + String(index)[0], "Required input."); + let start = validate_sr_parameter('#sr_start_no_' + String(index)[0], "Required input."); + let duration = validate_sr_parameter('#sr_duration_no_' + String(index)[0], "Required input."); + if (activity && start && duration) { + if (validate_sr_time('#sr_start_no_' + String(index)) && validate_sr_time('#sr_duration_no_' + String(index))) { + document.getElementById('sr_expiration_no_' + String(index)).disabled = true; + document.getElementById('sr_start_no_' + String(index)).disabled = true; + document.getElementById('sr_duration_no_' + String(index)).disabled = true; + document.getElementById('edit_row_no_' + String(index)).style.cssText = 'display:inline !important'; + $(this).closest(".form_field_outer_row").addClass("row_validated"); + $(this).hide(); + index = index + 1; + } + } + // On save, check open/unvalidated rows. + $(".validate_node_btn_frm_field").not(".row_validated").not(this).each(function( index ) { + index = $(this).attr('id').split('_').slice(-1)[0]; + if ($('#sr_start_no_' + String(index)[0]).val() != "") { + validate_sr_parameter('#sr_start_no_' + String(index)[0], "Required input.") + validate_sr_time('#sr_start_no_' + String(index)); + }; + if ($('#sr_duration_no_' + String(index)[0]).val() != "") { + validate_sr_parameter('#sr_duration_no_' + String(index)[0], "Required input."); + validate_sr_time('#sr_duration_no_' + String(index)); + } + }); + }); + + //Edit short-range activity type + $("body").on("click", ".edit_node_btn_frm_field", function() { + $(this).closest(".form_field_outer_row").removeClass("row_validated"); + $(this).hide(); + let id = $(this).attr('id').split('_').slice(-1)[0]; + document.getElementById('sr_expiration_no_' + String(id)).disabled = false; + document.getElementById('sr_start_no_' + String(id)).disabled = false; + document.getElementById('sr_duration_no_' + String(id)).disabled = false; + document.getElementById('validate_row_no_' + String(id)).style.cssText = 'display:inline !important'; + }) + + //Remove short-range interaction (modal field row). + $("body").on("click", ".remove_node_btn_frm_field", function() { + $(this).closest(".form_field_outer_row").remove(); + // On delete, check open/unvalidated rows. + $(".validate_node_btn_frm_field").not(".row_validated").not(this).each(function( index ) { + index = $(this).attr('id').split('_').slice(-1)[0]; + if ($('#sr_start_no_' + String(index)[0]).val() != "") { + validate_sr_parameter('#sr_start_no_' + String(index)[0], "Required input.") + validate_sr_time('#sr_start_no_' + String(index)); + }; + if ($('#sr_duration_no_' + String(index)[0]).val() != "") { + validate_sr_parameter('#sr_duration_no_' + String(index)[0], "Required input."); + validate_sr_time('#sr_duration_no_' + String(index)); + } + }); + }); + + //Short-range modal - close and save button + $("body").on("click", ".close_btn_frm_field", function() { + $(".validate_node_btn_frm_field").click(); + if ($(".form_field_outer").find(".form_field_outer_row.row_validated").length == $(".form_field_outer").find(".form_field_outer_row").length) { + $("#sr_interactions").text($(".form_field_outer").find(".form_field_outer_row.row_validated").length); + $(".form_field_outer_row").not(".row_validated").remove(); + $('#short_range_dialog').modal('hide'); + } + }); + + //Short-range modal - reset button + $("body").on("click", ".dismiss_btn_frm_field", function() { + $(".form_field_outer_row").remove(); + $("#sr_interactions").text(0); + $('input[type=radio][id=short_range_no]').prop("checked", true); + on_short_range_option_change(); + }); + }); diff --git a/cara/apps/calculator/static/js/report.js b/cara/apps/calculator/static/js/report.js index be601146..56ba86bc 100644 --- a/cara/apps/calculator/static/js/report.js +++ b/cara/apps/calculator/static/js/report.js @@ -1,56 +1,38 @@ /* Generate the concentration plot using d3 library. */ -function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses, exposed_presence_intervals) { +function draw_plot(svg_id) { - console.log(cumulative_doses) - - var time_format = d3.timeFormat('%H:%M'); + // Used for controlling the short-range interactions + let button_full_exposure = document.getElementById("button_full_exposure"); + let button_hide_high_concentration = document.getElementById("button_hide_high_concentration"); + let long_range_checkbox = document.getElementById('long_range_cumulative_checkbox') + let show_sr_legend = short_range_expirations.length > 0; var data_for_graphs = { 'concentrations': [], 'cumulative_doses': [], + 'long_range_cumulative_doses': [], } times.map((time, index) => data_for_graphs.concentrations.push({ 'time': time, 'hour': new Date().setHours(Math.trunc(time), (time - Math.trunc(time)) * 60), 'concentration': concentrations[index]})); times.map((time, index) => data_for_graphs.cumulative_doses.push({ 'time': time, 'hour': new Date().setHours(Math.trunc(time), (time - Math.trunc(time)) * 60), 'concentration': cumulative_doses[index]})); + times.map((time, index) => data_for_graphs.long_range_cumulative_doses.push({ 'time': time, 'hour': new Date().setHours(Math.trunc(time), (time - Math.trunc(time)) * 60), 'concentration': long_range_cumulative_doses[index]})); + + const tooltip_data_for_graphs = Object.fromEntries(Object.entries(data_for_graphs).filter(([key]) => !key.includes('long_range_cumulative_doses'))); // Add main SVG element var plot_div = document.getElementById(svg_id); var vis = d3.select(plot_div).append('svg'); + var time_format = d3.timeFormat('%H:%M'); // H:M time format for x axis. xRange = d3.scaleTime().domain([data_for_graphs.concentrations[0].hour, data_for_graphs.concentrations[data_for_graphs.concentrations.length - 1].hour]), xTimeRange = d3.scaleLinear().domain([data_for_graphs.concentrations[0].time, data_for_graphs.concentrations[data_for_graphs.concentrations.length - 1].time]), bisecHour = d3.bisector((d) => { return d.hour; }).left, - yRange = d3.scaleLinear().domain([0., Math.max(...concentrations)]), - yCumulativeRange = d3.scaleLinear().domain([0., Math.max(...cumulative_doses)*1.1]), - - xAxis = d3.axisBottom(xRange).tickFormat(d => time_format(d)), - yAxis = d3.axisLeft(yRange).ticks(4), - yCumulativeAxis = d3.axisRight(yCumulativeRange).ticks(4); - - // Line representing the mean concentration. - var lineFunc = d3.line(); - var draw_line = vis.append('svg:path') - .attr('stroke', '#1f77b4') - .attr('stroke-width', 2) - .attr('fill', 'none'); - - var lineCumulative = d3.line(); - var draw_cumulative_line = vis.append('svg:path') - .attr('stroke', '#1f77b4') - .attr('stroke-width', 2) - .style("stroke-dasharray", "5 5") - .attr('fill', 'none'); - - // Area representing the presence of exposed person(s). - var exposedArea = {}; - var drawArea = {}; - exposed_presence_intervals.forEach((b, index) => { - exposedArea[index] = d3.area(); - drawArea[index] = vis.append('svg:path') - .attr('fill', '#1f77b4') - .attr('fill-opacity', '0.1'); - }); + yRange = d3.scaleLinear(), + yCumulativeRange = d3.scaleLinear(), + + yAxis = d3.axisLeft(), + yCumulativeAxis = d3.axisRight(); // X axis declaration. var xAxisEl = vis.append('svg:g') @@ -77,7 +59,6 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses // Y cumulative concentration axis declaration. var yAxisCumEl = vis.append('svg:g') .attr('class', 'y axis') - .style('font-size', 14) .style("stroke-dasharray", "5 5"); // Y cumulated concentration axis label. @@ -88,41 +69,76 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses .text('Mean cumulative dose (infectious virus)'); // Legend for the plot elements - line and area. + + // Concentration line icon var legendLineIcon = vis.append('rect') .attr('width', 20) .attr('height', 3) .style('fill', '#1f77b4'); + // Concentration line text + var legendLineText = vis.append('text') + .text('Mean concentration') + .style('font-size', '15px'); + // Cumulative dose line icon var legendCumulativeIcon = vis.append('line') .style("stroke-dasharray", "5 5") //dashed array for line .attr('stroke-width', '2') .style("stroke", '#1f77b4'); - - var legendAreaIcon = vis.append('rect') - .attr('width', 20) - .attr('height', 20) - .attr('fill', '#1f77b4') - .attr('fill-opacity', '0.1'); - - var legendLineText = vis.append('text') - .text('Mean concentration') - .style('font-size', '15px') - .attr('alignment-baseline', 'central'); - + // Cumulative dose line text var legendCumutiveText = vis.append('text') .text('Cumulative dose') - .style('font-size', '15px') - .attr('alignment-baseline', 'central'); + .style('font-size', '15px'); + // Area line icon + var legendAreaIcon = vis.append('rect') + .attr('width', 20) + .attr('height', 15) + .attr('fill', '#1f77b4') + .attr('fill-opacity', '0.1'); + // Area line text var legendAreaText = vis.append('text') .text('Presence of exposed person(s)') - .style('font-size', '15px') - .attr('alignment-baseline', 'central'); + .style('font-size', '15px'); + + sr_unique_activities = [...new Set(short_range_expirations)] + if (show_sr_legend) { + // Long range cumulative dose line legend - line and area + var legendLongCumulativeIcon = vis.append('line') + .style("stroke-dasharray", "5 5") //dashed array for line + .attr('stroke-width', '2') + .style("stroke", 'purple') + .attr('opacity', 0); + var legendLongCumutiveText = vis.append('text') + .text('Long-range cumulative dose') + .style('font-size', '15px') + .attr('opacity', 0); + // Short range area icon + var legendShortRangeAreaIcon = {}; + sr_unique_activities.forEach((b, index) => { + legendShortRangeAreaIcon[index] = vis.append('rect') + .attr('width', 20) + .attr('height', 15); + // Short range area icon colors + if (sr_unique_activities[index] == 'Breathing') legendShortRangeAreaIcon[index].attr('fill', 'red').attr('fill-opacity', '0.2'); + else if (sr_unique_activities[index] == 'Speaking') legendShortRangeAreaIcon[index].attr('fill', 'green').attr('fill-opacity', '0.1'); + else legendShortRangeAreaIcon[index].attr('fill', 'blue').attr('fill-opacity', '0.1'); + }); + // Short range area text + var legendShortRangeText = {}; + sr_unique_activities.forEach((b, index) => { + legendShortRangeText[index] = vis.append('text') + .text('Short-range - ' + sr_unique_activities[index]) + .style('font-size', '15px'); + }); + } // Legend bounding + if (show_sr_legend) legendBBox_height = 68 + 20 * sr_unique_activities.length; + else legendBBox_height = 68; var legendBBox = vis.append('rect') .attr('width', 255) - .attr('height', 70) + .attr('height', legendBBox_height) .attr('stroke', 'lightgrey') .attr('stroke-width', '2') .attr('rx', '5px') @@ -130,9 +146,64 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses .attr('stroke-linejoin', 'round') .attr('fill', 'none'); + var clip = vis.append("defs").append("svg:clipPath") + .attr("id", "clip") + .append("svg:rect"); + + var draw_area = vis.append('svg:g') + .attr('clip-path', 'url(#clip)'); + // Line representing the mean concentration. + var lineFunc = d3.line(); + draw_area.append('svg:path') + .attr('class', 'line') + .attr('stroke', '#1f77b4') + .attr('stroke-width', 2) + .attr('fill', 'none'); + + // Line representing the cumulative concentration. + var lineCumulative = d3.line(); + var draw_cumulative_line = draw_area.append('svg:path') + .attr('stroke', '#1f77b4') + .attr('stroke-width', 2) + .style("stroke-dasharray", "5 5") + .attr('fill', 'none'); + + // Line representing the long-range cumulative concentration. + if (show_sr_legend) { + var longRangeCumulative = d3.line(); + var draw_long_range_cumulative_line = draw_area.append('svg:path') + .attr('stroke', 'purple') + .attr('stroke-width', 2) + .style("stroke-dasharray", "5 5") + .attr('fill', 'none') + .attr('opacity', 0); + } + + // Area representing the presence of exposed person(s). + var exposedArea = {}; + var drawArea = {}; + exposed_presence_intervals.forEach((b, index) => { + exposedArea[index] = d3.area(); + drawArea[index] = draw_area.append('svg:path') + .attr('fill', '#1f77b4') + .attr('fill-opacity', '0.1'); + }); + + // Area representing the short-range interaction(s). + var shortRangeArea = {}; + var drawShortRangeArea = {}; + short_range_intervals.forEach((b, index) => { + shortRangeArea[index] = d3.area(); + drawShortRangeArea[index] = draw_area.append('svg:path'); + + if (short_range_expirations[index] == 'Breathing') drawShortRangeArea[index].attr('fill', 'red').attr('fill-opacity', '0.2'); + else if (short_range_expirations[index] == 'Speaking') drawShortRangeArea[index].attr('fill', 'green').attr('fill-opacity', '0.1'); + else drawShortRangeArea[index].attr('fill', 'blue').attr('fill-opacity', '0.1'); + }); + // Tooltip. var focus = {}, tooltip_rect = {}, tooltip_time = {}, tooltip_concentration = {}, toolBox = {}; - for (const [concentration, data] of Object.entries(data_for_graphs)) { + for (const [concentration, data] of Object.entries(tooltip_data_for_graphs)) { focus[concentration] = vis.append('svg:g') .style('display', 'none'); @@ -145,7 +216,6 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses .attr('stroke', '#000') .attr('width', 85) .attr('height', 50) - .attr('x', 10) .attr('y', -22) .attr('rx', 4) .attr('ry', 4); @@ -165,11 +235,66 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses .attr('pointer-events', 'all') .on('mouseover', () => { for (const [concentration, data] of Object.entries(focus)) focus[concentration].style('display', null); }) .on('mouseout', () => { for (const [concentration, data] of Object.entries(focus)) focus[concentration].style('display', 'none'); }) - .on('mousemove', mousemove); + .on('mousemove', mousemove);; } - var graph_width; - var graph_height; + function update_concentration_plot(concentration_data, cumulative_data) { + yRange.domain([0., Math.max(...concentration_data)*1.1]); + yAxisEl.transition().duration(1000).call(yAxis); + + yCumulativeRange.domain([0., Math.max(...cumulative_data)*1.1]); + yAxisCumEl.transition().duration(1000).call(yCumulativeAxis) + + // Concentration line + lineFunc.defined(d => !isNaN(d.concentration)) + .x(d => xTimeRange(d.time)) + .y(d => yRange(d.concentration)); + draw_area.select('.line') + .transition() + .duration(1000) + .attr("d", lineFunc(data_for_graphs.concentrations)); + + // Cumulative line. + lineCumulative.defined(d => !isNaN(d.concentration)) + .x(d => xTimeRange(d.time)) + .y(d => yCumulativeRange(d.concentration)); + draw_cumulative_line.transition() + .duration(1000) + .attr("d", lineCumulative(data_for_graphs.cumulative_doses)); + + // Long-range cumulative line. + if (show_sr_legend) { + longRangeCumulative.defined(d => !isNaN(d.concentration)) + .x(d => xTimeRange(d.time)) + .y(d => yCumulativeRange(d.concentration)); + draw_long_range_cumulative_line.transition() + .duration(1000) + .attr("d", lineCumulative(data_for_graphs.long_range_cumulative_doses)); + } + + // Area. + exposed_presence_intervals.forEach((b, index) => { + exposedArea[index].x(d => xTimeRange(d.time)) + .y0(graph_height - 50) + .y1(d => yRange(d.concentration) + ); + drawArea[index].transition().duration(1000).attr('d', exposedArea[index](data_for_graphs.concentrations.filter(d => { + return d.time >= b[0] && d.time <= b[1] + }))); + }); + + // Short-Range Area. + short_range_intervals.forEach((b, index) => { + shortRangeArea[index].x(d => xTimeRange(d.time)) + .y0(graph_height - 50) + .y1(d => yRange(d.concentration)); + + drawShortRangeArea[index].transition().duration(1000).attr('d', shortRangeArea[index](data_for_graphs.concentrations.filter(d => { + return d.time >= b[0] && d.time <= b[1] + }))); + }); + + } function redraw() { @@ -182,7 +307,7 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses var margins = { top: 30, right: 20, bottom: 50, left: 60 }; div_width = 900; graph_width = div_width * (2/3); - const svg_margins = {'margin-left': '0rem', 'margin-top': '0rem'}; + const svg_margins = {'margin-left': '0rem'}; Object.entries(svg_margins).forEach(([prop,val]) => vis.style(prop,val)); } else { @@ -190,7 +315,7 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses div_width = div_width * 1.1 graph_width = div_width * .9; graph_height = div_height * 0.65; // On mobile screen sizes we want the legend to be on the bottom of the graph. - const svg_margins = {'margin-left': '-1rem', 'margin-top': '3rem'}; + const svg_margins = {'margin-left': '-1rem'}; Object.entries(svg_margins).forEach(([prop,val]) => vis.style(prop,val)); }; @@ -200,45 +325,29 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses // SVG components according to the width and height. + // clipPath: everything out of this area won't be drawn. + clip.attr("x", margins.left) + .attr("y", margins.top) + .attr("width", graph_width - margins.right - margins.left) + .attr("height", graph_height - margins.top - margins.bottom); + // Axis ranges. xRange.range([margins.left, graph_width - margins.right]); xTimeRange.range([margins.left, graph_width - margins.right]); yRange.range([graph_height - margins.bottom, margins.top]); yCumulativeRange.range([graph_height - margins.bottom, margins.top]); - // Line. - lineFunc.defined(d => !isNaN(d.concentration)) - .x(d => xTimeRange(d.time)) - .y(d => yRange(d.concentration)); - draw_line.attr("d", lineFunc(data_for_graphs.concentrations)); - - // Cumulative line - lineCumulative.defined(d => !isNaN(d.concentration)) - .x(d => xTimeRange(d.time)) - .y(d => yCumulativeRange(d.concentration)); - draw_cumulative_line.attr("d", lineCumulative(data_for_graphs.cumulative_doses)); - - // Area. - exposed_presence_intervals.forEach((b, index) => { - exposedArea[index].x(d => xTimeRange(d.time)) - .y0(graph_height - margins.bottom) - .y1(d => yRange(d.concentration)); - - drawArea[index].attr('d', exposedArea[index](data_for_graphs.concentrations.filter(d => { - return d.time >= b[0] && d.time <= b[1] - }))); - }); - // Axis. var xAxis = d3.axisBottom(xRange).tickFormat(d => time_format(d)); - var yAxis = d3.axisLeft(yRange); + yAxis.scale(yRange); + yCumulativeAxis.scale(yCumulativeRange); xAxisEl.attr('transform', 'translate(0,' + (graph_height - margins.bottom) + ')') .call(xAxis); xAxisLabelEl.attr('x', (graph_width + margins.right) / 2) .attr('y', graph_height * 0.97); - - yAxisEl.attr('transform', 'translate(' + margins.left + ',0)').call(yAxis); + + yAxisEl.attr('transform', 'translate(' + margins.left + ',0)'); yAxisLabelEl.attr('x', (graph_height * 0.9 + margins.bottom) / 2) .attr('y', (graph_height + margins.left) * 0.9) .attr('transform', 'rotate(-90, 0,' + graph_height + ')'); @@ -258,64 +367,129 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses .attr('y', graph_width + 290); } - // Legend on right side. const size = 20; + var legend_x_start = 50; + const space_between_text_icon = 30; + const text_height = 6; + // Legend on right side. if (plot_div.clientWidth >= 900) { - legendLineIcon.attr('x', graph_width + size * 2.5) - .attr('y', margins.top + size); - legendLineText.attr('x', graph_width + 4 * size) + legendLineIcon.attr('x', graph_width + legend_x_start) .attr('y', margins.top + size); + legendLineText.attr('x', graph_width + legend_x_start + space_between_text_icon) + .attr('y', margins.top + size + text_height); + + legendCumulativeIcon.attr("x1", graph_width + legend_x_start) + .attr("x2", graph_width + legend_x_start + 20) + .attr("y1", margins.top + 2 * size) + .attr("y2", margins.top + 2 * size); + legendCumutiveText.attr('x', graph_width + legend_x_start + space_between_text_icon) + .attr('y', margins.top + 2 * size + text_height); - legendCumulativeIcon.attr("x1", graph_width + size + 30) - .attr("x2", graph_width + 2 * size + 32) - .attr("y1", 3.5 * size) - .attr("y2", 3.5 * size); - legendCumutiveText.attr('x', graph_width + 2.5 * size + 30) - .attr('y', margins.top + 2 * size); + legendAreaIcon.attr('x', graph_width + legend_x_start) + .attr('y', margins.top + (3 * size) - 15/2); + legendAreaText.attr('x', graph_width + legend_x_start + space_between_text_icon) + .attr('y', margins.top + 3 * size + text_height); - legendAreaIcon.attr('x', graph_width + size * 2.5) - .attr('y', margins.top + 2.5 * size); - legendAreaText.attr('x', graph_width + 4 * size) - .attr('y', margins.top + 3 * size); + if (show_sr_legend) { + sr_unique_activities.forEach((b, index) => { + legendShortRangeAreaIcon[index].attr('x', graph_width + legend_x_start) + .attr('y', margins.top + (4 + index) * size - 15/2); + legendShortRangeText[index].attr('x', graph_width + legend_x_start + space_between_text_icon) + .attr('y', margins.top + (4 + index) * size + text_height); + }); + legendLongCumulativeIcon.attr("x1", graph_width + legend_x_start) + .attr("x2", graph_width + legend_x_start + 20) + .attr("y1", margins.top + (4 + sr_unique_activities.length) * size) + .attr("y2", margins.top + (4 + sr_unique_activities.length) * size); + legendLongCumutiveText.attr('x', graph_width + legend_x_start + space_between_text_icon) + .attr('y', margins.top + (4 + sr_unique_activities.length) * size + + text_height); + } legendBBox.attr('x', graph_width * 1.07) .attr('y', margins.top * 1.2); } // Legend on the bottom. else { - legendLineIcon.attr('x', size * 0.5) - .attr('y', graph_height * 1.05); - legendLineText.attr('x', 2 * size) - .attr('y', graph_height * 1.05); + legend_x_start = 10; - legendCumulativeIcon.attr("x1", size * 0.5) - .attr("x2", size * 1.55) - .attr("y1", graph_height * 1.05 + size) - .attr("y2", graph_height * 1.05 + size); - legendCumutiveText.attr('x', 2 * size) - .attr('y', graph_height + 1.65 * size); + legendLineIcon.attr('x', legend_x_start) + .attr('y', graph_height + size); + legendLineText.attr('x', legend_x_start + space_between_text_icon) + .attr('y', graph_height + size + text_height); - legendAreaIcon.attr('x', size * 0.50) - .attr('y', graph_height * 1.09 + size); - legendAreaText.attr('x', 2 * size) - .attr('y', graph_height + 2.7 * size); + legendCumulativeIcon.attr("x1", legend_x_start) + .attr("x2", legend_x_start + 20) + .attr("y1", graph_height + 2 * size) + .attr("y2", graph_height + 2 * size); + legendCumutiveText.attr('x', legend_x_start + space_between_text_icon) + .attr('y', graph_height + 2 * size + text_height); + + legendAreaIcon.attr('x', legend_x_start) + .attr('y', graph_height + 3 * size - 15/2); + legendAreaText.attr('x', legend_x_start + space_between_text_icon) + .attr('y', graph_height + 3 * size + text_height); + + if (show_sr_legend) { + sr_unique_activities.forEach((b, index) => { + legendShortRangeAreaIcon[index].attr('x', legend_x_start) + .attr('y', graph_height + 4 * size - 15/2); + legendShortRangeText[index].attr('x', legend_x_start + space_between_text_icon) + .attr('y', graph_height + 4 * size + text_height); + }); + legendLongCumulativeIcon.attr("x1", legend_x_start) + .attr("x2", legend_x_start + 20) + .attr("y1", graph_height + (4 + sr_unique_activities.length) * size) + .attr("y2", graph_height + (4 + sr_unique_activities.length) * size) + legendLongCumutiveText.attr('x', legend_x_start + space_between_text_icon) + .attr('y', graph_height + (4 + sr_unique_activities.length) * size + text_height); + } legendBBox.attr('x', 1) - .attr('y', graph_height); + .attr('y', graph_height + 6); } - + // ToolBox. - for (const [concentration, data] of Object.entries(data_for_graphs)) { + for (const [concentration, data] of Object.entries(tooltip_data_for_graphs)) { toolBox[concentration].attr('width', graph_width - margins.right) .attr('height', graph_height); } + } - // Draw for the first time to initialize. - redraw(); + if (show_sr_legend) { + long_range_checkbox.addEventListener("click", () => { + if (long_range_checkbox.checked) { + draw_long_range_cumulative_line.transition().duration(1000).attr("opacity", 1); + legendBBox.transition().duration(1000).attr("height", legendBBox_height + 20); + legendLongCumulativeIcon.transition().duration(1000).attr("opacity", 1); + legendLongCumutiveText.transition().duration(1000).attr("opacity", 1); + } + else { + draw_long_range_cumulative_line.transition().duration(1000).attr("opacity", 0); + legendBBox.transition().duration(1000).attr("height", legendBBox_height); + legendLongCumulativeIcon.transition().duration(1000).attr("opacity", 0); + legendLongCumutiveText.transition().duration(1000).attr("opacity", 0); + } + }); + }; + + if (button_full_exposure) { + button_full_exposure.addEventListener("click", () => { + update_concentration_plot(concentrations, cumulative_doses); + button_full_exposure.disabled = true; + button_hide_high_concentration.disabled = false; + }); + } + if (button_hide_high_concentration) { + button_hide_high_concentration.addEventListener("click", () => { + update_concentration_plot(concentrations_zoomed, long_range_cumulative_doses); + button_full_exposure.disabled = false; + button_hide_high_concentration.disabled = true; + }); + } function mousemove() { - for (const [scenario, data] of Object.entries(data_for_graphs)) { + for (const [scenario, data] of Object.entries(tooltip_data_for_graphs)) { if (d3.pointer(event)[0] < graph_width / 2) { tooltip_rect[scenario].attr('x', 10) tooltip_time[scenario].attr('x', 18) @@ -351,27 +525,39 @@ function draw_concentration_plot(svg_id, times, concentrations, cumulative_doses } } - // Redraw based on the new size whenever the browser window is resized. - window.addEventListener("resize", redraw); + // Draw for the first time to initialize. + redraw(); + update_concentration_plot(concentrations, cumulative_doses); - + // Redraw based on the new size whenever the browser window is resized. + window.addEventListener("resize", e => { + redraw(); + if (button_full_exposure && button_full_exposure.disabled) update_concentration_plot(concentrations, cumulative_doses); + else update_concentration_plot(concentrations_zoomed, long_range_cumulative_doses) + }); } + // 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(concentration_plot_svg_id, alternative_plot_svg_id, times, alternative_scenarios) { - // H:M format +// The method is prepared to consider short range interactions if needed. +function draw_alternative_scenarios_plot(concentration_plot_svg_id, alternative_plot_svg_id) { + // 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; + // Used for controlling the short-range interactions + let button_full_exposure = document.getElementById("button_alternative_full_exposure"); + let button_hide_high_concentration = document.getElementById("button_alternative_hide_high_concentration"); + // 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 + scenario_concentrations = alternative_scenarios[scenario].concentrations; highest_concentration = Math.max(highest_concentration, Math.max(...scenario_concentrations)) @@ -393,34 +579,9 @@ function draw_alternative_scenarios_plot(concentration_plot_svg_id, alternative_ var xTimeRange = d3.scaleLinear().domain([times[0], times[times.length - 1]]); var bisecHour = d3.bisector((d) => { return d.hour; }).left; - var yRange = d3.scaleLinear().domain([0., highest_concentration]); + var yRange = d3.scaleLinear(); + var yAxis = d3.axisLeft(); - // Line representing the mean concentration for each scenario. - var lineFuncs = {}, draw_lines = {}, label_icons = {}, label_text = {}; - 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. - lineFuncs[scenario_name] = d3.line(); - - draw_lines[scenario_name] = vis.append('svg:path') - .attr("stroke", colors[scenario_index]) - .attr('stroke-width', 2) - .attr('fill', 'none'); - - // Legend for the plot elements - lines. - label_icons[scenario_name] = vis.append('rect') - .attr('width', 20) - .attr('height', 3) - .style('fill', colors[scenario_index]); - - label_text[scenario_name] = vis.append('text') - .text(scenario_name) - .style('font-size', '15px') - .attr('alignment-baseline', 'central'); - - } - // X axis. var xAxisEl = vis.append('svg:g') .attr('class', 'x axis'); @@ -454,6 +615,38 @@ function draw_alternative_scenarios_plot(concentration_plot_svg_id, alternative_ .attr('stroke-linejoin', 'round') .attr('fill', 'none'); + var clip = vis.append("defs").append("svg:clipPath") + .attr("id", "clip") + .append("svg:rect"); + + var draw_area = vis.append('svg:g') + .attr('clip-path', 'url(#clip)'); + + // Line representing the mean concentration for each scenario. + var lineFuncs = {}, draw_lines = {}, label_icons = {}, label_text = {}; + 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. + lineFuncs[scenario_name] = d3.line(); + + draw_lines[scenario_name] = draw_area.append('svg:path') + .attr("stroke", colors[scenario_index]) + .attr('stroke-width', 2) + .attr('fill', 'none'); + + // Legend for the plot elements - lines. + label_icons[scenario_name] = vis.append('rect') + .attr('width', 20) + .attr('height', 3) + .style('fill', colors[scenario_index]); + + label_text[scenario_name] = vis.append('text') + .text(scenario_name) + .style('font-size', '15px'); + + } + // Tooltip. var focus = {}, tooltip_rect = {}, tooltip_time = {}, tooltip_concentration = {}, toolBox = {}; for (const [scenario_name, data] of Object.entries(data_for_scenarios)) { @@ -489,6 +682,26 @@ function draw_alternative_scenarios_plot(concentration_plot_svg_id, alternative_ .on('mousemove', mousemove); } + function update_alternative_concentration_plot(concentration_data) { + var highest_concentration = 0. + + for (scenario in alternative_scenarios) { + scenario_concentrations = alternative_scenarios[scenario][concentration_data]; + highest_concentration = Math.max(highest_concentration, Math.max(...scenario_concentrations)); + } + + yRange.domain([0., highest_concentration]); + yAxisEl.transition().duration(1000).call(yAxis); + + for (const [scenario_name, data] of Object.entries(data_for_scenarios)) { + // Lines. + lineFuncs[scenario_name].defined(d => !isNaN(d.concentration)) + .x(d => xTimeRange(d.time)) + .y(d => yRange(d.concentration)); + draw_lines[scenario_name].transition().duration(1000).attr("d", lineFuncs[scenario_name](data)); + } + } + var graph_width; var graph_height; @@ -519,41 +732,45 @@ function draw_alternative_scenarios_plot(concentration_plot_svg_id, alternative_ .attr('height', div_height); // SVG components according to the width and height. + // clipPath: everything out of this area won't be drawn. + clip.attr("x", margins.left) + .attr("y", margins.top) + .attr("width", graph_width - margins.right - margins.left) + .attr("height", graph_height - margins.top - margins.bottom); // Axis ranges. xRange.range([margins.left, graph_width - margins.right]); xTimeRange.range([margins.left, graph_width - margins.right]); yRange.range([graph_height - margins.bottom, margins.top]); + var legend_x_start = 25; + const space_between_text_icon = 30; + const text_height = 6; for (const [scenario_name, data] of Object.entries(data_for_scenarios)) { var scenario_index = Object.keys(data_for_scenarios).indexOf(scenario_name) - // Lines. - lineFuncs[scenario_name].defined(d => !isNaN(d.concentration)) - .x(d => xTimeRange(d.time)) - .y(d => yRange(d.concentration)); - draw_lines[scenario_name].attr("d", lineFuncs[scenario_name](data)); - // Legend on right side. var size = 20 * (scenario_index + 1); if (document.getElementById(concentration_plot_svg_id).clientWidth >= 900) { - label_icons[scenario_name].attr('x', graph_width + 20) - .attr('y', margins.top + size); - label_text[scenario_name].attr('x', graph_width + 3 * 20) + label_icons[scenario_name].attr('x', graph_width + legend_x_start) .attr('y', margins.top + size); + label_text[scenario_name].attr('x', graph_width + legend_x_start + space_between_text_icon) + .attr('y', margins.top + size + text_height); } // Legend on the bottom. else { - label_icons[scenario_name].attr('x', margins.left * 0.3) - .attr('y', graph_height + size); - label_text[scenario_name].attr('x', margins.left * 1.4) + legend_x_start = 10; + + label_icons[scenario_name].attr('x', legend_x_start) .attr('y', graph_height + size); + label_text[scenario_name].attr('x', legend_x_start + space_between_text_icon) + .attr('y', graph_height + size + text_height); } } // Axis. var xAxis = d3.axisBottom(xRange).tickFormat(d => time_format(d)); - var yAxis = d3.axisLeft(yRange); + yAxis.scale(yRange); xAxisEl.attr('transform', 'translate(0,' + (graph_height - margins.bottom) + ')') .call(xAxis); @@ -585,8 +802,20 @@ function draw_alternative_scenarios_plot(concentration_plot_svg_id, alternative_ } } - // Draw for the first time to initialize. - redraw(); + if (button_full_exposure) { + button_full_exposure.addEventListener("click", () => { + update_alternative_concentration_plot('concentrations'); + button_full_exposure.disabled = true; + button_hide_high_concentration.disabled = false; + }); + } + if (button_hide_high_concentration) { + button_hide_high_concentration.addEventListener("click", () => { + update_alternative_concentration_plot('concentrations_zoomed'); + button_full_exposure.disabled = false; + button_hide_high_concentration.disabled = true; + }); + } function mousemove() { for (const [scenario_name, data] of Object.entries(data_for_scenarios)) { @@ -611,8 +840,16 @@ function draw_alternative_scenarios_plot(concentration_plot_svg_id, alternative_ } } + // Draw for the first time to initialize. + redraw(); + update_alternative_concentration_plot('concentrations'); + // Redraw based on the new size whenever the browser window is resized. - window.addEventListener("resize", redraw); + window.addEventListener("resize", e => { + redraw(); + if (button_full_exposure && button_full_exposure.disabled) update_alternative_concentration_plot('concentrations'); + else update_alternative_concentration_plot('concentrations') + }); } function copy_clipboard(shareable_link) { diff --git a/cara/apps/expert.py b/cara/apps/expert.py index bba12b6b..84beed27 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -503,6 +503,7 @@ baseline_model = models.ExposureModel( ), evaporation_factor=0.3, ), + short_range=(), exposed=models.Population( number=10, presence=models.SpecificInterval(((8., 12.), (13., 17.))), diff --git a/cara/apps/templates/about.html.j2 b/cara/apps/templates/about.html.j2 index 50d2a8bf..47e29072 100644 --- a/cara/apps/templates/about.html.j2 +++ b/cara/apps/templates/about.html.j2 @@ -17,8 +17,9 @@ CARA stands for COVID Airborne Risk Assessment and was developed in the spring o
  • CARA expert app
    • -The mathematical and physical model simulate the long-range airborne spread of SARS-CoV-2 virus in a finite volume, assuming a homogenous mixture, and estimates the risk of COVID-19 airborne transmission therein. The results DO NOT include (for now) short-range airborne exposure (where the physical distance plays a factor) nor the other known modes of SARS-CoV-2 transmission. Hence, the output from this model is only valid when the other recommended public health & safety instructions are observed, such as adequate physical distancing, good hand hygiene and other barrier measures.
    +The mathematical and physical model simulate the airborne spread of SARS-CoV-2 virus in a finite volume, assuming a homogenous mixture and a two-stage exhaled jet model, and estimates the risk of COVID-19 airborne transmission therein. The results DO NOT include other known modes of SARS-CoV-2 transmission. Hence, the output from this model is only valid when the other recommended public health & safety instructions are observed, such as good hand hygiene and other barrier measures.

    The methodology, mathematical equations and parameters of the model are published here in the CARA paper: Modelling airborne transmission of SARS-CoV-2 using CARA: risk assessment for enclosed spaces.

    +

    Note that the short-range component of the model has not yet been published.

    The model used is based on scientific publications relating to airborne transmission of infectious diseases, virology, epidemiology and aerosol science. It can be used to compare the effectiveness of different airborne-related risk mitigation measures. diff --git a/cara/apps/templates/base/calculator.form.html.j2 b/cara/apps/templates/base/calculator.form.html.j2 index 84b3e987..499fc761 100644 --- a/cara/apps/templates/base/calculator.form.html.j2 +++ b/cara/apps/templates/base/calculator.form.html.j2 @@ -366,6 +366,47 @@
    +
    +
    Short-range interactions (without masks):
    +
    + + + + +
    +
    +

    The use of masks mitigates exposure at short-range. The analytical model with short-range interactions does not take mask wearing into account.

    + +
    +
    + +

    0 short-range interactions.

    +
    +
    + + + +
    +
    @@ -515,7 +556,7 @@
    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.
    + 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 the original "wild type" strain of the virus and three variants of concern (VOC):
      diff --git a/cara/apps/templates/base/calculator.report.html.j2 b/cara/apps/templates/base/calculator.report.html.j2 index 52d8ad6c..45a8ff0d 100644 --- a/cara/apps/templates/base/calculator.report.html.j2 +++ b/cara/apps/templates/base/calculator.report.html.j2 @@ -24,7 +24,7 @@

      REPORT

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

    - +
    @@ -78,7 +78,7 @@
    {% block report_summary %} {% endblock report_summary %}
    @@ -89,71 +89,89 @@ {% endblock report_summary_footnote %}

    * The results are based on the parameters and assumptions published in the CARA publication: doi.org/10.1098/rsfs.2021.0076.


    - + {% if form.short_range_option == "short_range_yes" %} + {% if 'Speaking' in form.short_range_interactions|string or 'Shouting' in form.short_range_interactions|string %} + + + {% endif %} + + {% endif %}

    -
    -
    Alternative scenarios - -
    -
    -
    -
    -
    - -
    - {% block report_scenarios_summary_table %} - - - - - - - - - - {% for scenario_name, scenario_stats in alternative_scenarios.stats.items() %} - - - - - - {% endfor %} - -
    ScenarioP(I)Expected new cases
    {{ scenario_name }} {{ scenario_stats.probability_of_infection | non_zero_percentage }}{{ scenario_stats.expected_new_cases | float_format }}
    - {% endblock report_scenarios_summary_table %} + + {% if form.short_range_option == "short_range_no" %} +
    +
    Alternative scenarios + +
    +
    +
    +
    + {% if form.short_range_option == "short_range_yes" %} + {% if 'Speaking' in form.short_range_interactions|string or 'Shouting' in form.short_range_interactions|string %} + + + {% endif %} + {% endif %} +
    + +
    + {% block report_scenarios_summary_table %} + + + + + + + + + + {% for scenario_name, scenario_stats in alternative_scenarios.stats.items() %} + + + + + + {% endfor %} + +
    ScenarioP(I)Expected new cases
    {{ scenario_name }} {{ scenario_stats.probability_of_infection | non_zero_percentage }}{{ scenario_stats.expected_new_cases | float_format }}
    + {% endblock report_scenarios_summary_table %} +
    +
    +

    Notes for alternative scenarios:
    +

      +
    1. This graph shows the concentration of infectious quanta in the air. The filtration of Type I and FFP2 masks, if worn, applies not only to the emission rate but also to the individual exposure (i.e. inhalation). + For this reason, scenarios with different types of mask will show the same concentration on the graph but have different absorbed doses and infection probabilities.
      2. +
    2. If you have selected more sophisticated options, such as HEPA filtration or FFP2 masks, this will be indicated in the plot as the "base scenario", representing the inputs inserted in the form.
      + The other alternative scenarios shown for comparison will not include either HEPA filtration or FFP2 masks.
      3. +
    +
    +

    -
    -

    Notes for alternative scenarios:
    -

      -
    1. This graph shows the concentration of infectious quanta in the air. The filtration of Type I and FFP2 masks, if worn, applies not only to the emission rate but also to the individual exposure (i.e. inhalation). - For this reason, scenarios with different types of mask will show the same concentration on the graph but have different absorbed doses and infection probabilities.
      2. -
    2. If you have selected more sophisticated options, such as HEPA filtration or FFP2 masks, this will be indicated in the plot as the "base scenario", representing the inputs inserted in the form.
      - The other alternative scenarios shown for comparison will not include either HEPA filtration or FFP2 masks.
      3. -
    -
    -

    -
    + {% endif %} {% endblock report_results %} {% block report_footer %} @@ -311,6 +329,18 @@ Gym = For comparison only, all persons doing heavy physical exercise, breathing and not speaking. {% endif %}

    + {% if form.short_range_option == "short_range_yes" %} +

  • + Short-range interactions: {{ form.short_range_interactions|length }} +

    • +
      + {% for interaction in form.short_range_interactions %} +
    • Expiratory activity {{ loop.index if form.short_range_interactions|length > 1 }}: {{ interaction.expiration }}
      • +
    • Start time {{ loop.index if form.short_range_interactions|length > 1 }}: {{ interaction.start_time }}
      • +
    • Duration {{ loop.index if form.short_range_interactions|length > 1 }}: {{ interaction.duration }} {{ "minutes" if interaction.duration|float > 1 else "minute" }}
      • + {% endfor %} +
    + {% endif %}

  • Exposed occupant(s) activity time:

    • Start time: {{ form.exposed_start | minutes_to_time }}

      • diff --git a/cara/apps/templates/cern/calculator.report.html.j2 b/cara/apps/templates/cern/calculator.report.html.j2 index 4d8c5d90..11c6f5c2 100644 --- a/cara/apps/templates/cern/calculator.report.html.j2 +++ b/cara/apps/templates/cern/calculator.report.html.j2 @@ -28,22 +28,22 @@ {% endblock warning_animation %} {% block report_summary %} + {% set report_message = "Taking into account the uncertainties tied to the model variables, in this scenario and assuming all occupants are exposed equally, the probability of one exposed occupant getting infected is " + prob_inf | non_zero_percentage + " and the expected number of new cases is " + expected_new_cases | float_format + "*." %}
      {% if scale_warning == 'red' %} {% elif scale_warning == 'orange' %} {% elif scale_warning == 'green' %} {% endif %} diff --git a/cara/models.py b/cara/models.py index f667ba1f..eb7324d5 100644 --- a/cara/models.py +++ b/cara/models.py @@ -37,7 +37,6 @@ import typing import numpy as np from scipy.interpolate import interp1d -import scipy.integrate if not typing.TYPE_CHECKING: from memoization import cached @@ -634,6 +633,12 @@ class _ExpirationBase: """ raise NotImplementedError("Subclass must implement") + def jet_origin_concentration(self): + """ + concentration of viruses at the jet origin (mL/m3). + """ + raise NotImplementedError("Subclass must implement") + @dataclass(frozen=True) class Expiration(_ExpirationBase): @@ -667,6 +672,14 @@ class Expiration(_ExpirationBase): return self.cn * (volume(self.diameter) * (1 - mask.exhale_efficiency(self.diameter))) * 1e-12 + @cached() + def jet_origin_concentration(self): + def volume(d): + return (np.pi * d**3) / 6. + + # final result converted from microns^3/cm3 to mL/m3 + return self.cn * volume(self.diameter) * 1e-6 + @dataclass(frozen=True) class MultipleExpiration(_ExpirationBase): @@ -886,7 +899,7 @@ class InfectedPopulation(_PopulationWithVirus): class ConcentrationModel: room: Room ventilation: _VentilationBase - infected: _PopulationWithVirus + infected: InfectedPopulation #: evaporation factor: the particles' diameter is multiplied by this # factor as soon as they are in the air (but AFTER going out of the, @@ -987,7 +1000,7 @@ class ConcentrationModel: def _normed_concentration(self, time: float) -> _VectorisedFloat: """ - Virus exposure concentration, as a function of time, and + Virus long-range 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 @@ -1013,18 +1026,18 @@ class ConcentrationModel: def concentration(self, time: float) -> _VectorisedFloat: """ - Virus exposure concentration, as a function of time. + Virus long-range 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) * + return (self._normed_concentration_cached(time) * self.infected.emission_rate_when_present()) @method_cache def normed_integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat: """ - Get the integrated concentration of viruses in the air between the times start and stop, + Get the integrated long-range concentration of viruses in the air between the times start and stop, normalized by the emission rate. """ if stop <= self._first_presence_time(): @@ -1059,6 +1072,117 @@ class ConcentrationModel: self.infected.emission_rate_when_present()) +@dataclass(frozen=True) +class ShortRangeModel: + #: Expiration type + expiration: _ExpirationBase + + #: Activity type + activity: Activity + + #: Short-range expiration and respective presence + presence: SpecificInterval + + #: Interpersonal distances + distance: _VectorisedFloat + + def dilution_factor(self) -> _VectorisedFloat: + ''' + The dilution factor for the respective expiratory activity type. + ''' + # Average mouth diameter + D = 0.02 + # Convert Breathing rate from m3/h to m3/s + BR = np.array(self.activity.exhalation_rate/3600.) + # Area of the mouth assuming a perfect circle + Am = np.pi*(D**2)/4 + # Initial velocity from the division of the Breathing rate with the area + u0 = np.array(BR/Am) + + tstar = 2.0 + Cr1 = 0.18 + Cr2 = 0.2 + Cx1 = 2.4 + + # The expired flow rate during the expiration period, m^3/s + Q0 = u0 * np.pi/4*D**2 + # Parameters in the jet-like stage + x01 = D/2/Cr1 + # Time of virtual origin + t01 = (x01/Cx1)**2 * (Q0*u0)**(-0.5) + # The transition point, m + xstar = np.array(Cx1*(Q0*u0)**0.25*(tstar + t01)**0.5 - x01) + # Dilution factor at the transition point xstar + Sxstar = np.array(2*Cr1*(xstar+x01)/D) + + distances = np.array(self.distance) + + factors = np.empty(distances.shape, dtype=np.float64) + factors[distances < xstar] = 2*Cr1*(distances[distances < xstar] + + x01)/D + factors[distances >= xstar] = Sxstar[distances >= xstar]*(1 + + Cr2*(distances[distances >= xstar] - + xstar[distances >= xstar])/Cr1/(xstar[distances >= xstar] + + x01))**3 + return factors + + def _normed_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: + """ + Virus short-range exposure concentration, as a function of time. + + If the given time falls within a short-range interval it returns the + short-range concentration normalized by the virus viral load. Otherwise + it returns 0. + """ + start, stop = self.presence.boundaries()[0] + # Verifies if the given time falls within a short-range interaction + if start <= time <= stop: + dilution = self.dilution_factor() + jet_origin_concentration = self.expiration.jet_origin_concentration() + # Long-range concentration normalized by the virus viral load + long_range_normed_concentration = (concentration_model.concentration(time) / + concentration_model.virus.viral_load_in_sputum) + + # The long-range concentration values are then approximated using interpolation: + # The set of points where we want the interpolated values are the short-range particle diameters (given the current expiration); + # The set of points with a known value are the long-range particle diameters (given the initial expiration); + # The set of known values are the long-range concentration values normalized by the viral load. + long_range_normed_concentration_interpolated=np.interp(self.expiration.particle.diameter, + concentration_model.infected.particle.diameter, long_range_normed_concentration) + + # Short-range concentration formula. The long-range concentration is added in the concentration method (ExposureModel). + return ((1/dilution)*(jet_origin_concentration - long_range_normed_concentration_interpolated)) + return 0. + + def short_range_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: + """ + Virus short-range exposure concentration, as a function of time. + """ + return (self._normed_concentration(concentration_model, time) * + concentration_model.virus.viral_load_in_sputum) + + @method_cache + def _normed_short_range_concentration_cached(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat: + # A cached version of the _normed_concentration method. Use this + # method if you expect that there may be multiple short-range concentration + # calculations for the same time (e.g. at state change times). + return self._normed_concentration(concentration_model, time) + + def normed_exposure_between_bounds(self, concentration_model: ConcentrationModel, time1: float, time2: float): + """ + Get the integrated short-range concentration of viruses in the air between the times start and stop, + normalized by the virus viral load. + """ + start_bound, stop_bound = self.presence.boundaries()[0] + + jet_origin_integrated = self.expiration.jet_origin_concentration() + dilution = self.dilution_factor() + + total_normed_concentration = -(concentration_model.integrated_concentration(start_bound, stop_bound)/concentration_model.virus.viral_load_in_sputum/dilution) + total_normed_concentration_interpolated = np.interp(self.expiration.particle.diameter, concentration_model.infected.particle.diameter, total_normed_concentration) + return (jet_origin_integrated/dilution * (stop_bound - start_bound)) + total_normed_concentration_interpolated + + @dataclass(frozen=True) class ExposureModel: """ @@ -1071,13 +1195,16 @@ class ExposureModel: #: The virus concentration model which this exposure model should consider. concentration_model: ConcentrationModel + #: The list of short-range models which this exposure model should consider. + short_range: typing.Tuple[ShortRangeModel, ...] + #: The population of non-infected people to be used in the model. exposed: Population #: The number of times the exposure event is repeated (default 1). repeats: int = 1 - def fraction_deposited(self) -> _VectorisedFloat: + def long_range_fraction_deposited(self) -> _VectorisedFloat: """ The fraction of particles actually deposited in the respiratory tract (over the total number of particles). It depends on the @@ -1086,7 +1213,7 @@ class ExposureModel: return self.concentration_model.infected.particle.fraction_deposited( self.concentration_model.evaporation_factor) - def _normed_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: + def _long_range_normed_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: """The number of virions per meter^3 between any two times, normalized by the emission rate of the infected population""" exposure = 0. @@ -1104,28 +1231,26 @@ class ExposureModel: elif time1 <= start and stop < time2: exposure += self.concentration_model.normed_integrated_concentration(start, stop) return exposure - - def _normed_exposure(self) -> _VectorisedFloat: - """ - The number of virions per meter^3, normalized by the emission rate - of the infected population. - """ - normed_exposure = 0.0 - for start, stop in self.exposed.presence.boundaries(): - normed_exposure += self.concentration_model.normed_integrated_concentration(start, stop) - - return normed_exposure * self.repeats - - def deposited_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: + def concentration(self, time: float) -> _VectorisedFloat: """ - The number of virus per m^3 deposited on the respiratory tract - between any two times. + Virus exposure concentration, as a function of time. + + It considers the long-range concentration with the + contribution of the short-range concentration. """ + concentration = self.concentration_model.concentration(time) + for interaction in self.short_range: + concentration += interaction.short_range_concentration(self.concentration_model, time) + return concentration + + def long_range_deposited_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: + deposited_exposure = 0. + emission_rate_per_aerosol = self.concentration_model.infected.emission_rate_per_aerosol_when_present() aerosols = self.concentration_model.infected.aerosols() - fdep = self.fraction_deposited() f_inf = self.concentration_model.infected.fraction_of_infectious_virus() + fdep = self.long_range_fraction_deposited() diameter = self.concentration_model.infected.particle.diameter @@ -1134,20 +1259,74 @@ class ExposureModel: # to perform properly the Monte-Carlo integration over # particle diameters (doing things in another order would # lead to wrong results). - dep_exposure_integrated = np.array(self._normed_exposure_between_bounds(time1, time2) * - aerosols * - fdep).mean() + dep_exposure_integrated = np.array(self._long_range_normed_exposure_between_bounds(time1, time2) * + aerosols * + fdep).mean() else: # in the case of a single diameter or no diameter defined, # one should not take any mean at this stage. - dep_exposure_integrated = self._normed_exposure_between_bounds(time1, time2)*aerosols*fdep + dep_exposure_integrated = self._long_range_normed_exposure_between_bounds(time1, time2)*aerosols*fdep # then we multiply by the diameter-independent quantity emission_rate_per_aerosol, - # and parameters of the vD equation (i.e. f_inf, BR_k and n_in). - return (dep_exposure_integrated * emission_rate_per_aerosol * - f_inf * self.exposed.activity.inhalation_rate * + # and parameters of the vD equation (i.e. BR_k and n_in). + deposited_exposure += (dep_exposure_integrated * emission_rate_per_aerosol * + self.exposed.activity.inhalation_rate * (1 - self.exposed.mask.inhale_efficiency())) + # In the end we multiply the final results by the fraction of infectious virus of the vD equation. + return deposited_exposure * f_inf + + def deposited_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat: + """ + The number of virus per m^3 deposited on the respiratory tract + between any two times. + + Considers a contribution between the short-range and long-range exposures: + It calculates the deposited exposure given a short-range interaction (if any). + Then, the deposited exposure given the long-range interactions is added to the + initial deposited exposure. + """ + deposited_exposure = 0. + for interaction in self.short_range: + start, stop = interaction.presence.boundaries()[0] + if stop < time1: + continue + elif start > time2: + break + elif start <= time1 and time2<= stop: + start_bound, stop_bound = time1, time2 + elif start <= time1 and stop < time2: + start_bound, stop_bound = time1, stop + elif time1 < start and time2 <= stop: + start_bound, stop_bound = start, time2 + elif time1 <= start and stop < time2: + start_bound, stop_bound = start, stop + short_range_exposure = interaction.normed_exposure_between_bounds(self.concentration_model, start_bound, stop_bound) + + fdep = interaction.expiration.particle.fraction_deposited(evaporation_factor=1.0) + diameter = interaction.expiration.particle.diameter + + # Aerosols not considered given the formula for the initial concentration at mouth/nose. + if diameter is not None and not np.isscalar(diameter): + # we compute first the mean of all diameter-dependent quantities + # to perform properly the Monte-Carlo integration over + # particle diameters (doing things in another order would + # lead to wrong results). + deposited_exposure += np.array(short_range_exposure * + fdep).mean() + else: + # in the case of a single diameter or no diameter defined, + # one should not take any mean at this stage. + deposited_exposure += short_range_exposure*fdep + + # then we multiply by the diameter-independent quantity virus viral load + deposited_exposure *= self.concentration_model.virus.viral_load_in_sputum + # long-range concentration + f_inf = self.concentration_model.infected.fraction_of_infectious_virus() + deposited_exposure += self.long_range_deposited_exposure_between_bounds(time1, time2)/f_inf + + return deposited_exposure * f_inf + def deposited_exposure(self) -> _VectorisedFloat: """ The number of virus per m^3 deposited on the respiratory tract. @@ -1162,7 +1341,7 @@ class ExposureModel: def infection_probability(self) -> _VectorisedFloat: # viral dose (vD) vD = self.deposited_exposure() - + # oneoverln2 multiplied by ID_50 corresponds to ID_63. infectious_dose = oneoverln2 * self.concentration_model.virus.infectious_dose diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 290f5920..2c6268e6 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -5,11 +5,13 @@ import numpy as np from scipy import special as sp import cara.monte_carlo as mc -from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel,CustomKernel,Uniform +from cara.monte_carlo.sampleable import LogNormal,LogCustomKernel,CustomKernel,Uniform + sqrt2pi = np.sqrt(2.*np.pi) sqrt2 = np.sqrt(2.) + @dataclass(frozen=True) class BLOmodel: """ @@ -65,7 +67,7 @@ class BLOmodel: return result -# From https://doi.org/10.1101/2021.10.14.21264988 and refererences therein +# From https://doi.org/10.1101/2021.10.14.21264988 and references therein activity_distributions = { 'Seated': mc.Activity(LogNormal(-0.6872121723362303, 0.10498338229297108), LogNormal(-0.6872121723362303, 0.10498338229297108)), @@ -84,7 +86,7 @@ activity_distributions = { } -# From https://doi.org/10.1101/2021.10.14.21264988 and refererences therein +# From https://doi.org/10.1101/2021.10.14.21264988 and references therein symptomatic_vl_frequencies = LogCustomKernel( np.array((2.46032, 2.67431, 2.85434, 3.06155, 3.25856, 3.47256, 3.66957, 3.85979, 4.09927, 4.27081, 4.47631, 4.66653, 4.87204, 5.10302, 5.27456, 5.46478, 5.6533, 5.88428, 6.07281, 6.30549, @@ -157,7 +159,10 @@ mask_distributions = { } -def expiration_distribution(BLO_factors): +def expiration_distribution( + BLO_factors, + d_max=30., +) -> mc.Expiration: """ Returns an Expiration with an aerosol diameter distribution, defined by the BLO factors (a length-3 tuple). @@ -166,10 +171,15 @@ def expiration_distribution(BLO_factors): an historical choice based on previous implementations of the model (it limits the influence of the O-mode). """ - dscan = np.linspace(0.1, 30. ,3000) - return mc.Expiration(CustomKernel(dscan, - BLOmodel(BLO_factors).distribution(dscan),kernel_bandwidth=0.1), - cn=BLOmodel(BLO_factors).integrate(0.1, 30.)) + dscan = np.linspace(0.1, d_max, 3000) + return mc.Expiration( + CustomKernel( + dscan, + BLOmodel(BLO_factors).distribution(dscan), + kernel_bandwidth=0.1, + ), + cn=BLOmodel(BLO_factors).integrate(0.1, d_max), + ) expiration_BLO_factors = { @@ -182,5 +192,15 @@ expiration_BLO_factors = { expiration_distributions = { exp_type: expiration_distribution(BLO_factors) - for exp_type,BLO_factors in expiration_BLO_factors.items() + for exp_type, BLO_factors in expiration_BLO_factors.items() } + + +short_range_expiration_distributions = { + exp_type: expiration_distribution(BLO_factors, d_max=100) + for exp_type, BLO_factors in expiration_BLO_factors.items() +} + + +# Fit from Fig 8 a) "stand-stand" in https://www.mdpi.com/1660-4601/17/4/1445/htm +short_range_distances = LogNormal(-0.269359136417347, 0.4728300188814934) diff --git a/cara/monte_carlo/models.py b/cara/monte_carlo/models.py index 6f09023e..f47d3416 100644 --- a/cara/monte_carlo/models.py +++ b/cara/monte_carlo/models.py @@ -37,7 +37,7 @@ class MCModelBase(typing.Generic[_ModelType]): def build_model(self, size: int) -> _ModelType: """ - Turn this MCModelBase subclass into a cara.models Model instance + Turn this MCModelBase subclass into a cara.model Model instance from which you can then run the model. """ @@ -72,6 +72,9 @@ def _build_mc_model(model: _ModelType) -> typing.Type[MCModelBase[_ModelType]]: elif new_field.type == typing.Tuple[cara.models._ExpirationBase, ...]: EB = getattr(sys.modules[__name__], "_ExpirationBase") field_type = typing.Tuple[typing.Union[cara.models._ExpirationBase, EB], ...] + elif new_field.type == typing.Tuple[cara.models.SpecificInterval, ...]: + SI = getattr(sys.modules[__name__], "SpecificInterval") + field_type = typing.Tuple[typing.Union[cara.models.SpecificInterval, SI], ...] else: # Check that we don't need to do anything with this type. for item in new_field.type.__args__: diff --git a/cara/tests/conftest.py b/cara/tests/conftest.py index 27ce9f1d..3cce3eb3 100644 --- a/cara/tests/conftest.py +++ b/cara/tests/conftest.py @@ -6,7 +6,7 @@ import pytest @pytest.fixture -def baseline_model(): +def baseline_concentration_model(): model = models.ConcentrationModel( room=models.Room(volume=75), ventilation=models.AirChange( @@ -30,14 +30,20 @@ def baseline_model(): @pytest.fixture -def baseline_exposure_model(baseline_model): +def baseline_sr_model(): + return () + + +@pytest.fixture +def baseline_exposure_model(baseline_concentration_model, baseline_sr_model): return models.ExposureModel( - baseline_model, + baseline_concentration_model, + baseline_sr_model, exposed=models.Population( number=1000, - presence=baseline_model.infected.presence, - activity=baseline_model.infected.activity, - mask=baseline_model.infected.mask, + presence=baseline_concentration_model.infected.presence, + activity=baseline_concentration_model.infected.activity, + mask=baseline_concentration_model.infected.mask, host_immunity=0., ), ) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 238e2316..44979024 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -90,8 +90,8 @@ def known_concentrations(func): np.array([40.91708675, 91.46172332]), np.array([51.6749232285, 80.3196524031])], ]) def test_exposure_model_ndarray(population, cm, - expected_exposure, expected_probability): - model = ExposureModel(cm, population) + expected_exposure, expected_probability, sr_model): + model = ExposureModel(cm, sr_model, population) np.testing.assert_almost_equal( model.deposited_exposure(), expected_exposure ) @@ -110,10 +110,10 @@ def test_exposure_model_ndarray(population, cm, [populations[1], np.array([2.13410688, 1.98167067])], [populations[2], np.array([1.36390289, 1.52436206])], ]) -def test_exposure_model_ndarray_and_float_mix(population, expected_deposited_exposure): +def test_exposure_model_ndarray_and_float_mix(population, expected_deposited_exposure, sr_model): cm = known_concentrations( lambda t: 0. if np.floor(t) % 2 else np.array([1.2, 1.2])) - model = ExposureModel(cm, population) + model = ExposureModel(cm, sr_model, population) np.testing.assert_almost_equal( model.deposited_exposure(), expected_deposited_exposure @@ -128,17 +128,17 @@ def test_exposure_model_ndarray_and_float_mix(population, expected_deposited_exp [populations[1], np.array([2.13410688, 1.98167067])], [populations[2], np.array([1.36390289, 1.52436206])], ]) -def test_exposure_model_vector(population, expected_deposited_exposure): +def test_exposure_model_vector(population, expected_deposited_exposure, sr_model): cm_array = known_concentrations(lambda t: np.array([1.2, 1.2])) - model_array = ExposureModel(cm_array, population) + model_array = ExposureModel(cm_array, sr_model, population) np.testing.assert_almost_equal( model_array.deposited_exposure(), np.array(expected_deposited_exposure) ) -def test_exposure_model_scalar(): +def test_exposure_model_scalar(sr_model): cm_scalar = known_concentrations(lambda t: 1.2) - model_scalar = ExposureModel(cm_scalar, populations[0]) + model_scalar = ExposureModel(cm_scalar, sr_model, populations[0]) expected_deposited_exposure = 1.52436206 np.testing.assert_almost_equal( model_scalar.deposited_exposure(), expected_deposited_exposure @@ -169,6 +169,11 @@ def conc_model(): ) +@pytest.fixture +def sr_model(): + return () + + # Expected deposited exposure were computed with a trapezoidal integration, using # a mesh of 10'000 pts per exposed presence interval. @pytest.mark.parametrize( @@ -183,17 +188,17 @@ def conc_model(): ] ) def test_exposure_model_integral_accuracy(exposed_time_interval, - expected_deposited_exposure, conc_model): + expected_deposited_exposure, conc_model, sr_model): presence_interval = models.SpecificInterval((exposed_time_interval,)) population = models.Population( 10, presence_interval, models.Mask.types['Type I'], models.Activity.types['Standing'], 0., ) - model = ExposureModel(conc_model, population) + model = ExposureModel(conc_model, sr_model, population) np.testing.assert_allclose(model.deposited_exposure(), expected_deposited_exposure) -def test_infectious_dose_vectorisation(): +def test_infectious_dose_vectorisation(sr_model): infected_population = models.InfectedPopulation( number=1, presence=halftime, @@ -216,7 +221,7 @@ def test_infectious_dose_vectorisation(): 10, presence_interval, models.Mask.types['Type I'], models.Activity.types['Standing'], 0., ) - model = ExposureModel(cm, population) + model = ExposureModel(cm, sr_model, population) inf_probability = model.infection_probability() assert isinstance(inf_probability, np.ndarray) assert inf_probability.shape == (3, ) diff --git a/cara/tests/models/test_short_range_model.py b/cara/tests/models/test_short_range_model.py new file mode 100644 index 00000000..a07a656b --- /dev/null +++ b/cara/tests/models/test_short_range_model.py @@ -0,0 +1,105 @@ +import typing + +import numpy as np +import pytest + +from cara import models +import cara.monte_carlo as mc_models +from cara.apps.calculator.model_generator import build_expiration +from cara.monte_carlo.data import short_range_expiration_distributions, short_range_distances, activity_distributions + +# TODO: seed better the random number generators +np.random.seed(2000) + +@pytest.fixture +def concentration_model() -> mc_models.ConcentrationModel: + return mc_models.ConcentrationModel( + room=models.Room(volume=75), + ventilation=models.AirChange( + active=models.SpecificInterval(present_times=((8.5, 12.5), (13.5, 17.5))), + air_exch=10_000_000., + ), + infected=mc_models.InfectedPopulation( + number=1, + virus=models.Virus.types['SARS_CoV_2'], + presence=models.SpecificInterval(present_times=((8.5, 12.5), (13.5, 17.5))), + mask=models.Mask.types['No mask'], + activity=models.Activity.types['Light activity'], + expiration=build_expiration({'Speaking': 0.33, 'Breathing': 0.67}), + host_immunity=0., + ), + evaporation_factor=0.3, + ) + + +@pytest.fixture +def short_range_model(): + return mc_models.ShortRangeModel(expiration=short_range_expiration_distributions['Breathing'], + activity=activity_distributions['Seated'], + presence=models.SpecificInterval(present_times=((10.5, 11.0),)), + distance=short_range_distances) + + +def test_short_range_model_ndarray(concentration_model, short_range_model): + concentration_model = concentration_model.build_model(250_000) + model = short_range_model.build_model(250_000) + assert isinstance(model._normed_concentration(concentration_model, 10.75), np.ndarray) + assert isinstance(model.short_range_concentration(concentration_model, 10.75), np.ndarray) + assert isinstance(model.normed_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray) + assert isinstance(model.short_range_concentration(concentration_model, 14.0), float) + + +@pytest.mark.parametrize( + "activity, expected_dilution", [ + ["Seated", 176.04075727780327], + ["Standing", 157.12965288170005], + ["Light activity", 69.06672998536413], + ["Moderate activity", 47.165817446310115], + ["Heavy exercise", 23.759992220217875], + ] +) +def test_dilution_factor(activity, expected_dilution): + model = models.ShortRangeModel(expiration="Breathing", + activity=models.Activity.types[activity], + presence=models.SpecificInterval(present_times=((10.5, 11.0),)), + distance=0.854) + assert isinstance(model.dilution_factor(), np.ndarray) + np.testing.assert_almost_equal( + model.dilution_factor(), expected_dilution, decimal=10 + ) + + +@pytest.mark.parametrize( + "time, expected_short_range_concentration", [ + [8.5, 0.], + [10.5, 15.24806213], + [10.6, 15.24806213], + [11.0, 15.24806213], + [12.0, 0.], + ] +) +def test_short_range_concentration(time, expected_short_range_concentration, concentration_model, short_range_model): + concentration_model = concentration_model.build_model(250_000) + model = short_range_model.build_model(250_000) + np.testing.assert_allclose( + np.array(model.short_range_concentration(concentration_model, time)).mean(), + expected_short_range_concentration, rtol=0.01 + ) + +@pytest.mark.parametrize( + "start, stop, expected_exposure", [ + [8.5, 12.5, 7.875963317294013e-09], + [10.5, 11.0, 7.875963317294013e-09], + [10.4, 11.1, 7.875963317294013e-09], + [10.5, 11.1, 7.875963317294013e-09], + [10.6, 11.1, 7.66539809488759e-09], + [10.4, 10.9, 7.66539809488759e-09], + + ] +) +def test_normed_exposure_between_bounds(start, stop, expected_exposure, concentration_model, short_range_model): + concentration_model = concentration_model.build_model(250_000) + model = short_range_model.build_model(250_000) + np.testing.assert_almost_equal( + model.normed_exposure_between_bounds(concentration_model, start, stop).mean(), expected_exposure + ) diff --git a/cara/tests/test_full_algorithm.py b/cara/tests/test_full_algorithm.py new file mode 100644 index 00000000..56264d24 --- /dev/null +++ b/cara/tests/test_full_algorithm.py @@ -0,0 +1,383 @@ +from dataclasses import dataclass +import typing + +import numpy as np +from scipy.integrate import quad +from scipy.special import erf +import numpy.testing as npt +import pytest + +import cara.monte_carlo as mc +from cara import models,data +from cara.utils import method_cache +from cara.models import _VectorisedFloat,Interval,SpecificInterval +from cara.monte_carlo.sampleable import LogNormal +from cara.monte_carlo.data import (expiration_distributions, + expiration_BLO_factors,short_range_expiration_distributions, + short_range_distances) + +# TODO: seed better the random number generators +np.random.seed(2000) +SAMPLE_SIZE = 500000 +TOLERANCE = 0.05 + +sqrt2pi = np.sqrt(2.*np.pi) +sqrt2 = np.sqrt(2.) +ln2 = np.log(2) + +@dataclass(frozen=True) +class SimpleConcentrationModel: + """ + Simple model for the background (long-range) concentration, without + all the flexibility of cara.models.ConcentrationModel. + For independent, end-to-end testing purposes. + This assumes no mask wearing, and the same ventilation rate at all + times. + """ + + #: infected people presence interval + infected_presence: Interval + + #: viral load (RNA copies / mL) + viral_load: _VectorisedFloat + + #: breathing rate (m^3/h) + breathing_rate: _VectorisedFloat + + #: room volume (m^3) + room_volume: _VectorisedFloat + + #: ventilation rate (air changes per hour) - including HEPA + lambda_ventilation: _VectorisedFloat + + #: BLO factors + BLO_factors: typing.Tuple[float, float, float] + + #: number of infected people + num_infected: int = 1 + + #: relative humidity RH + humidity: float = 0.3 + + #: minimum particle diameter considered (microns) + diameter_min: float = 0.1 + + #: maximum particle diameter considered (microns) + diameter_max: float = 30. + + #: evaporation factor + evaporation: float = 0.3 + + #: cn (cm^-3) for resp. the B, L and O modes. Corresponds to the + # total concentration of aerosols for each mode. + cn: typing.Tuple[float, float, float] = (0.06, 0.2, 0.0010008) + + # mean of the underlying normal distributions (represents the log of a + # diameter in microns), for resp. the B, L and O modes. + mu: typing.Tuple[float, float, float] = (0.989541, 1.38629, 4.97673) + + # std deviation of the underlying normal distribution, for resp. + # the B, L and O modes. + sigma: typing.Tuple[float, float, float] = (0.262364, 0.506818, 0.585005) + + def removal_rate(self) -> _VectorisedFloat: + """ + removal rate lambda in h^-1, excluding the deposition rate. + """ + return (self.lambda_ventilation + + ln2/(6.43 if self.humidity<=0.4 else 1.1) ) + + @method_cache + def deposition_removal_coefficient(self) -> float: + """ + coefficient in front of gravitational deposition rate, in h^-1.microns^-2 + Note: 0.4512 = 1.88e-4 * 3600 / 1.5 + """ + return 0.4512*(self.evaporation/2.5)**2 + + @method_cache + def aerosol_volume(self,diameter: float) -> float: + """ + particle volume in microns^3 + """ + return 4*np.pi/3. * (diameter/2.)**3 + + @method_cache + def Np(self,diameter: float) -> float: + """ + number of emitted particles per unit volume (BLO model) + in cm^-3.ln(micron)^-1 + """ + result = 0. + for cn,mu,sigma,famp in zip(self.cn,self.mu,self.sigma, + self.BLO_factors): + result += ( (cn * famp)/sigma * + np.exp(-(np.log(diameter)-mu)**2/(2*sigma**2))) + return result/(diameter*sqrt2pi) + + def vR(self,diameter: float) -> float: + """ + emission rate per unit diameter, in RNA copies / h / micron + """ + return (self.viral_load * self.breathing_rate * self.Np(diameter) + * self.aerosol_volume(diameter) * 1e-6) + + @method_cache + def f(self, removal_rate: _VectorisedFloat, deltat: float) -> _VectorisedFloat: + """ + A general function to compute the main integral over diameters + """ + def integrand(diameter): + # function to return the integrand + a = self.deposition_removal_coefficient() + a_dsquare = a*diameter**2 + return (self.vR(diameter)/(a_dsquare + removal_rate) + * np.exp(-a_dsquare*deltat)) + + return quad(integrand,self.diameter_min,self.diameter_max)[0] + + def concentration(self,t: float) -> _VectorisedFloat: + """ + concentration at a given time t + """ + trans_times = sorted(self.infected_presence.transition_times()) + if t==trans_times[0]: + return 0. + + lambda_rate = self.removal_rate() + # transition_times[i] < t <= transition_times[i]+1 + i: int = np.searchsorted(trans_times,t) - 1 # type: ignore + ti = trans_times[i] + Pim1 = self.infected_presence.triggered((trans_times[i-1]+ti)/2.) + Pi = self.infected_presence.triggered((ti+trans_times[i+1])/2.) + + result = (0 if not Pim1 else self.f(lambda_rate,t-ti)) + result -= (0 if not Pi else self.f(lambda_rate,t-ti)) + + for k,tk in enumerate(trans_times[:i]): + Pkm1 = self.infected_presence.triggered((trans_times[k-1]+tk)/2.) + Pk = self.infected_presence.triggered((tk+trans_times[k+1])/2.) + s = np.sum([lambda_rate*(trans_times[l]-trans_times[l-1]) + for l in range(k+1,i+1)]) + result += ( (0 if not Pkm1 else self.f(lambda_rate,t-tk)) + -(0 if not Pk else self.f(lambda_rate,t-tk)) + ) * np.exp(-s) + + return ( ( (0 if not self.infected_presence.triggered(t) + else self.f(lambda_rate,0)) + + result * np.exp(-lambda_rate*(t-ti)) ) + * self.num_infected/self.room_volume) + + +@dataclass(frozen=True) +class SimpleShortRangeModel: + """ + Simple model for the short-range concentration, without + all the flexibility of cara.models.ShortRangeModel. + For independent, end-to-end testing purposes. + This assumes no mask wearing. + """ + + #: time intervals in which a short-range interaction occurs + interaction_interval: SpecificInterval + + #: tuple with interpersonal distanced from infected person (m) + distance : _VectorisedFloat = 0.854 + + #: breathing rate (m^3/h) + breathing_rate: _VectorisedFloat = 0.51 + + #: tuple with BLO factors + BLO_factors: typing.Tuple[float, float, float] = (1,0,0) + + #: minimum diameter for integration (short-range only) (microns) + diameter_min: float = 0.1 + + #: maximum diameter for integration (short-range only) (microns) + diameter_max: float = 100. + + #: mouth opening diameter (m) + D: float = 0.02 + + #: duration of the expiration (s) + tstar: float = 2. + + #: Streamwise and radial penetration coefficients + Cr1: float = 0.18 + Cx1: float = 2.4 + Cr2: float = 0.2 + Cx2: float = 2.2 + + @method_cache + def dilution_factor(self) -> _VectorisedFloat: + """ + computes dilution factor at a certain distance x + based on Wei JIA matlab script. + """ + x = np.array(self.distance) + dilution = np.empty(x.shape, dtype=np.float64) + # expired flow rate during the expiration period, m^3/s + Q0 = np.array(self.breathing_rate/3600) + # the expired flow velocity at the noozle (mouth opening), m/s + u0 = np.array(Q0/(np.pi/4. * self.D**2)) + # parameters in the jet-like stage + # position of virtual origin + x01 = self.D/2/self.Cr1 + # time of virtual origin + t01 = (x01/self.Cx1)**2 * (Q0*u0)**(-0.5) + # transition point (in m) + xstar = np.array(self.Cx1*(Q0*u0)**0.25*(self.tstar + t01)**0.5 + - x01) + # dilution factor at the transition point xstar + Sxstar = np.array(2.*self.Cr1*(xstar+x01)/self.D) + + # calculate dilution factor at the short-range distance x + dilution[x <= xstar] = 2.*self.Cr1*(x[x <= xstar] + x01)/self.D + dilution[x > xstar] = Sxstar[x > xstar]*(1. + self.Cr2*(x[x > xstar] + - xstar[x > xstar]) + /self.Cr1/(xstar[x > xstar] + x01))**3 + + return dilution + + @method_cache + def jet_concentration(self,conc_model: SimpleConcentrationModel) -> _VectorisedFloat: + """ + virion concentration at the origin of the jet (close to + the mouth of the infected person), in m^-3 + we perform the integral of Np(d)*V(d) over diameter analytically + """ + vl = conc_model.viral_load + dmin = self.diameter_min + dmax = self.diameter_max + result = 0. + for cn,mu,sigma,famp in zip(conc_model.cn,conc_model.mu,conc_model.sigma, + self.BLO_factors): + d0 = np.exp(mu) + ymin = (np.log(dmin)-mu)/(sqrt2*sigma)-3.*sigma/sqrt2 + ymax = (np.log(dmax)-mu)/(sqrt2*sigma)-3.*sigma/sqrt2 + result += ( (cn * famp * d0**3)/2. * np.exp(9*sigma**2/2.) * + (erf(ymax) - erf(ymin)) ) + return vl * 1e-6 * result * np.pi/6. + + def concentration(self, conc_model: SimpleConcentrationModel, time: float) -> _VectorisedFloat: + """ + compute the short-range part of the concentration, and add it + to the background concentration + """ + if self.interaction_interval.triggered(time): + background_concentration = conc_model.concentration(time) + S = self.dilution_factor() + return (self.jet_concentration(conc_model) + - background_concentration) / S + else: + return 0. + + +presence = models.SpecificInterval(present_times=((8.5, 12), (13, 17.5))) +interaction_intervals = (models.SpecificInterval(present_times=((10.5, 11.0),)), + models.SpecificInterval(present_times=((14.5, 15.0),)) + ) + + +@pytest.fixture +def conc_model() -> mc.ConcentrationModel: + return mc.ConcentrationModel( + room=models.Room(volume=50, humidity=0.3), + ventilation=models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=1.), + infected=mc.InfectedPopulation( + number=1, + presence=presence, + virus=models.Virus.types['SARS_CoV_2_DELTA'], + mask=models.Mask.types['No mask'], + activity=models.Activity.types['Seated'], + expiration=expiration_distributions['Breathing'], + host_immunity=0., + ), + evaporation_factor=0.3, + ).build_model(SAMPLE_SIZE) + + +@pytest.fixture +def simple_conc_model() -> SimpleConcentrationModel: + return SimpleConcentrationModel( + infected_presence = presence, + viral_load = models.Virus.types['SARS_CoV_2_DELTA'].viral_load_in_sputum, + breathing_rate = models.Activity.types['Seated'].exhalation_rate, + room_volume = 50., + lambda_ventilation= 1., + BLO_factors = expiration_BLO_factors['Breathing'], + ) + + +@pytest.fixture +def sr_models() -> typing.Tuple[mc.ShortRangeModel, ...]: + return ( + mc.ShortRangeModel( + expiration = short_range_expiration_distributions['Breathing'], + activity = models.Activity.types['Seated'], + presence = interaction_intervals[0], + distance = 0.854, + ).build_model(SAMPLE_SIZE), + mc.ShortRangeModel( + expiration = short_range_expiration_distributions['Speaking'], + activity = models.Activity.types['Seated'], + presence = interaction_intervals[1], + distance = 0.854, + ).build_model(SAMPLE_SIZE), + ) + + +@pytest.fixture +def simple_sr_models() -> typing.Tuple[SimpleShortRangeModel, ...]: + return ( + SimpleShortRangeModel( + interaction_interval = interaction_intervals[0], + distance = 0.854, + breathing_rate = models.Activity.types['Seated'].exhalation_rate, + BLO_factors = expiration_BLO_factors['Breathing'], + ), + SimpleShortRangeModel( + interaction_interval = interaction_intervals[1], + distance = 0.854, + breathing_rate = models.Activity.types['Seated'].exhalation_rate, + BLO_factors = expiration_BLO_factors['Speaking'], + ) + ) + + +@pytest.mark.parametrize( + "time", np.linspace(8.5,17.5,12), +) +def test_background_concentration(time,conc_model,simple_conc_model): + npt.assert_allclose( + conc_model.concentration(time).mean(), + simple_conc_model.concentration(time), rtol=TOLERANCE + ) + +@pytest.mark.parametrize( + "time", [10, 10.7, 11., 12.5, 14.75, 14.9, 17] +) +def test_shortrange_concentration(time,conc_model,simple_conc_model, + sr_models,simple_sr_models): + result_sr_model = np.sum([np.array( + sr_mod.short_range_concentration(conc_model,time)).mean() + for sr_mod in sr_models]) + result_simple_sr_model = np.sum([np.array( + sr_mod.concentration(simple_conc_model,time)).mean() + for sr_mod in simple_sr_models]) + npt.assert_allclose( + result_sr_model,result_simple_sr_model,rtol=TOLERANCE + ) + + +#times = np.linspace(8.5,17.5,120) +#plt.plot(times,[conc_model.concentration(t).mean() for t in times]) +#plt.plot(times,[simple_conc_model.concentration(t) for t in times]) + +#plt.plot(times,[np.sum([np.array(sr_mod.short_range_concentration(conc_model,t)).mean() +# for sr_mod in sr_models]) +# + conc_model.concentration(t).mean() for t in times]) +#plt.plot(times,[np.sum([np.array(sr_mod.concentration(simple_conc_model,t)).mean() +# for sr_mod in simple_sr_models]) +# + simple_conc_model.concentration(t) +# for t in times]) diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 2ab80fc8..78c4541f 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -6,10 +6,10 @@ import cara.models as models import cara.data as data -def test_no_mask_superspeading_emission_rate(baseline_model): +def test_no_mask_superspeading_emission_rate(baseline_concentration_model): expected_rate = 48500. npt.assert_allclose( - [baseline_model.infected.emission_rate(float(t)) for t in [0, 1, 4, 4.5, 5, 8, 9]], + [baseline_concentration_model.infected.emission_rate(float(t)) for t in [0, 1, 4, 4.5, 5, 8, 9]], [0, expected_rate, expected_rate, 0, 0, expected_rate, 0], rtol=1e-12 ) @@ -38,10 +38,10 @@ def baseline_periodic_hepa(): ) -def test_concentrations(baseline_model): +def test_concentrations(baseline_concentration_model): # expected concentrations were computed analytically ts = [0, 4, 5, 7, 10] - concentrations = [baseline_model.concentration(float(t)) for t in ts] + concentrations = [baseline_concentration_model.concentration(float(t)) for t in ts] npt.assert_allclose( concentrations, [0.000000e+00, 20.805628, 6.602814e-13, 20.805628, 2.09545e-26], @@ -49,13 +49,13 @@ def test_concentrations(baseline_model): ) -def test_smooth_concentrations(baseline_model): +def test_smooth_concentrations(baseline_concentration_model): # We don't care about the actual concentrations in this test, but rather # that the curve itself is smooth. dx = 0.002 dy_limit = 0.2 # Anything more than this (in relative) is a bit steep. ts = np.arange(0, 10, dx) - concentrations = [baseline_model.concentration(float(t)) for t in ts] + concentrations = [baseline_concentration_model.concentration(float(t)) for t in ts] assert np.abs(np.diff(concentrations)).max()/np.mean(concentrations) < dy_limit @@ -367,10 +367,11 @@ def test_concentrations_refine_times(time): npt.assert_allclose(m1.concentration(time), m2.concentration(time), rtol=1e-8) -def build_exposure_model(concentration_model): +def build_exposure_model(concentration_model, short_range_model): infected = concentration_model.infected return models.ExposureModel( concentration_model=concentration_model, + short_range=short_range_model, exposed=models.Population( number=10, presence=infected.presence, @@ -390,13 +391,13 @@ def build_exposure_model(concentration_model): ['Jun', 1721.03336729], ], ) -def test_exposure_hourly_dep(month,expected_deposited_exposure): +def test_exposure_hourly_dep(month,expected_deposited_exposure, baseline_sr_model): m = build_exposure_model( build_hourly_dependent_model( month, intervals_open=((0., 24.), ), intervals_presence_infected=((8., 12.), (13., 17.)) - ) + ), baseline_sr_model ) deposited_exposure = m.deposited_exposure() npt.assert_allclose(deposited_exposure, expected_deposited_exposure) @@ -411,14 +412,14 @@ def test_exposure_hourly_dep(month,expected_deposited_exposure): ['Jun', 1799.17597184], ], ) -def test_exposure_hourly_dep_refined(month,expected_deposited_exposure): +def test_exposure_hourly_dep_refined(month,expected_deposited_exposure, baseline_sr_model): m = build_exposure_model( build_hourly_dependent_model( month, intervals_open=((0., 24.),), intervals_presence_infected=((8., 12.), (13., 17.)), temperatures=data.GenevaTemperatures, - ) + ), baseline_sr_model ) deposited_exposure = m.deposited_exposure() npt.assert_allclose(deposited_exposure, expected_deposited_exposure, rtol=0.02) diff --git a/cara/tests/test_monte_carlo.py b/cara/tests/test_monte_carlo.py index 8a8b0271..a398e716 100644 --- a/cara/tests/test_monte_carlo.py +++ b/cara/tests/test_monte_carlo.py @@ -38,7 +38,7 @@ def test_type_annotations(): @pytest.fixture -def baseline_mc_model() -> cara.monte_carlo.ConcentrationModel: +def baseline_mc_concentration_model() -> cara.monte_carlo.ConcentrationModel: mc_model = cara.monte_carlo.ConcentrationModel( room=cara.monte_carlo.Room(volume=cara.monte_carlo.sampleable.Normal(75, 20)), ventilation=cara.monte_carlo.SlidingWindow( @@ -62,21 +62,27 @@ def baseline_mc_model() -> cara.monte_carlo.ConcentrationModel: @pytest.fixture -def baseline_mc_exposure_model(baseline_mc_model) -> cara.monte_carlo.ExposureModel: +def baseline_mc_sr_model() -> cara.monte_carlo.ShortRangeModel: + return () + + +@pytest.fixture +def baseline_mc_exposure_model(baseline_mc_concentration_model, baseline_mc_sr_model) -> cara.monte_carlo.ExposureModel: return cara.monte_carlo.ExposureModel( - baseline_mc_model, + baseline_mc_concentration_model, + baseline_mc_sr_model, exposed=cara.models.Population( number=10, - presence=baseline_mc_model.infected.presence, - activity=baseline_mc_model.infected.activity, - mask=baseline_mc_model.infected.mask, + presence=baseline_mc_concentration_model.infected.presence, + activity=baseline_mc_concentration_model.infected.activity, + mask=baseline_mc_concentration_model.infected.mask, host_immunity=0., ) ) -def test_build_concentration_model(baseline_mc_model: cara.monte_carlo.ConcentrationModel): - model = baseline_mc_model.build_model(7) +def test_build_concentration_model(baseline_mc_concentration_model: cara.monte_carlo.ConcentrationModel): + model = baseline_mc_concentration_model.build_model(7) assert isinstance(model, cara.models.ConcentrationModel) assert isinstance(model.concentration(time=0.), float) conc = model.concentration(time=1.) diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index b4aa247d..ecdb3374 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -9,8 +9,8 @@ from cara.apps.calculator.model_generator import build_expiration # TODO: seed better the random number generators np.random.seed(2000) -SAMPLE_SIZE = 250000 -TOLERANCE = 0.05 +SAMPLE_SIZE = 250_000 +TOLERANCE = 0.06 # Load the weather data (temperature in kelvin) for Toronto. toronto_coordinates = (43.667, 79.400) @@ -39,7 +39,6 @@ TorontoTemperatures = { # references values for infection_probability and expected new cases # in the following tests, were obtained from the feature/mc branch - @pytest.fixture def shared_office_mc(): """ @@ -72,6 +71,7 @@ def shared_office_mc(): ) return mc.ExposureModel( concentration_model=concentration_mc, + short_range=(), exposed=mc.Population( number=3, presence=mc.SpecificInterval(present_times=((0, 3.5), (4.5, 9))), @@ -114,6 +114,7 @@ def classroom_mc(): ) return mc.ExposureModel( concentration_model=concentration_mc, + short_range=(), exposed=mc.Population( number=19, presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), @@ -147,6 +148,7 @@ def ski_cabin_mc(): ) return mc.ExposureModel( concentration_model=concentration_mc, + short_range=(), exposed=mc.Population( number=3, presence=models.SpecificInterval(((0, 20/60),)), @@ -186,6 +188,7 @@ def skagit_chorale_mc(): ) return mc.ExposureModel( concentration_model=concentration_mc, + short_range=(), exposed=mc.Population( number=60, presence=models.SpecificInterval(((0, 2.5), )), @@ -225,6 +228,7 @@ def bus_ride_mc(): ) return mc.ExposureModel( concentration_model=concentration_mc, + short_range=(), exposed=mc.Population( number=67, presence=models.SpecificInterval(((0, 1.67), )), @@ -259,6 +263,7 @@ def gym_mc(): ) return mc.ExposureModel( concentration_model=concentration_mc, + short_range=(), exposed=mc.Population( number=28, presence=concentration_mc.infected.presence, @@ -293,6 +298,7 @@ def waiting_room_mc(): ) return mc.ExposureModel( concentration_model=concentration_mc, + short_range=(), exposed=mc.Population( number=14, presence=concentration_mc.infected.presence, @@ -370,6 +376,7 @@ def test_small_shared_office_Geneva(mask_type, month, expected_pi, ) exposure_mc = mc.ExposureModel( concentration_model=concentration_mc, + short_range=(), exposed=mc.Population( number=1, presence=concentration_mc.infected.presence,