From cbb3846eb7df1971669a45532765efd4449c49ba Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 1 Sep 2021 14:59:13 +0200 Subject: [PATCH 1/5] Added uniform distribution for masks --- cara/apps/calculator/model_generator.py | 7 +++++-- cara/monte_carlo/data.py | 9 ++++++++- cara/monte_carlo/sampleable.py | 12 ++++++++++++ 3 files changed, 25 insertions(+), 3 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index a4b139e7..cafc8077 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -10,7 +10,7 @@ from cara import models from cara import data import cara.monte_carlo as mc from .. import calculator -from cara.monte_carlo.data import activity_distributions, virus_distributions +from cara.monte_carlo.data import activity_distributions, virus_distributions, mask_distributions LOG = logging.getLogger(__name__) @@ -301,7 +301,10 @@ class FormData: def mask(self) -> models.Mask: # Initializes the mask type if mask wearing is "continuous", otherwise instantiates the mask attribute as # the "No mask"-mask - mask = models.Mask.types[self.mask_type if self.mask_wearing_option == "mask_on" else 'No mask'] + if self.mask_wearing_option == 'mask_on': + mask = mask_distributions[self.mask_type] + else: + mask = models.Mask.types['No mask'] return mask def infected_population(self) -> mc.InfectedPopulation: diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 613b830c..250bd6e6 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -1,7 +1,7 @@ import numpy as np import cara.monte_carlo as mc -from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel +from cara.monte_carlo.sampleable import Normal,LogNormal,LogCustomKernel, Uniform # From CERN-OPEN-2021-04 and refererences therein @@ -58,3 +58,10 @@ virus_distributions = { infectious_dose=60/1.6, ), } + + +# From CERN-OPEN-2021-04 and refererences therein +mask_distributions = { + 'Type I': mc.Mask(Uniform(0.25, 0.80)), + 'FFP2': mc.Mask(Uniform(0.83, 0.91)), +} \ No newline at end of file diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 4333c93f..27907e49 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -28,6 +28,18 @@ class Normal(SampleableDistribution): return np.random.normal(self.mean, self.standard_deviation, size=size) +class Uniform(SampleableDistribution): + """ + Defines a continuous uniform distribution + """ + def __init__(self, low: float, high: float): + self.low = low + self.high = high + + def generate_samples(self, size: int) -> float_array_size_n: + return np.random.uniform(self.low, self.high, size=size) + + class LogNormal(SampleableDistribution): """ Defines a lognormal distribution (i.e. Gaussian distribution vs. the From db02132853009295c74f3751ad8c8df606823b07 Mon Sep 17 00:00:00 2001 From: James Devine Date: Thu, 2 Sep 2021 11:58:32 +0200 Subject: [PATCH 2/5] Allow location selection as part of the calculator, with temp profile from local weather station --- README.md | 26 ++- cara/.gitattributes | 2 + cara/apps/calculator/__init__.py | 2 +- cara/apps/calculator/model_generator.py | 79 +++++++++- cara/apps/calculator/static/js/form.js | 135 +++++++++++++++- .../templates/base/calculator.report.html.j2 | 6 +- .../templates/calculator.form.html.j2 | 23 ++- cara/{data.py => data/__init__.py} | 2 + cara/data/global_weather_set.json | 3 + .../hadisd_station_fullinfo_v311_202001p.txt | 3 + cara/data/weather.py | 149 ++++++++++++++++++ .../apps/calculator/test_model_generator.py | 128 ++++++++++----- cara/tests/data/__init__.py | 0 cara/tests/data/test_weather.py | 74 +++++++++ requirements.txt | 1 + setup.cfg | 4 + setup.py | 5 + 17 files changed, 573 insertions(+), 69 deletions(-) create mode 100644 cara/.gitattributes rename cara/{data.py => data/__init__.py} (96%) create mode 100644 cara/data/global_weather_set.json create mode 100644 cara/data/hadisd_station_fullinfo_v311_202001p.txt create mode 100644 cara/data/weather.py create mode 100644 cara/tests/data/__init__.py create mode 100644 cara/tests/data/test_weather.py diff --git a/README.md b/README.md index 1be26e55..7853027d 100644 --- a/README.md +++ b/README.md @@ -58,19 +58,6 @@ CARA has not undergone review, approval or certification by competent authoritie The software is provided "as is", without warranty of any kind, express or implied, including but not limited to the warranties of merchantability, fitness for a particular purpose and non-infringement. In no event shall the authors or copyright holders be liable for any claim, damages or other liability, whether in an action of contract, tort or otherwise, arising from, out of or in connection with the software or the use or other dealings in the software. -## Adapting CARA to your location - -The default weather data (average hourly outdoor temperature in Celcius for each month of the year) used in CARA is for Geneva, Switzerland. -In order for the natural ventilation option to work correctly for other geographic locations, the outdoor temperatures must be updated. -There are some scripts to help download and process the temperature data from your nearest weather station in the https://gitlab.cern.ch/cara/climatology-data repository. -Once you have used the scripts, the hourly temperature data for your location should be added to the file `data.py` in place of the default values for Geneva. The temperature values for your locations should be pasted into the `Geneva_hourly_temperatures_celsius_per_hour` variable, **without changing the variable name** in the following format: - - `'Jan': [0.2, -0.3, -0.5, -0.9, -1.1, -1.4, -1.5, -1.5, -1.1, 0.1, 1.5, - 2.8, 3.8, 4.4, 4.5, 4.4, 4.4, 3.9, 3.1, 2.7, 2.2, 1.7, 1.5, 1.1], - 'Feb': [0.9, 0.3, 0.0, -0.5, -0.7, -1.1, -1.2, -1.1, -0.7, 0.8, 2.5, - 4.2, 5.4, 6.2, 6.3, 6.2, 6.1, 5.5, 4.5, 4.1, 3.5, 2.8, 2.5, 2.0],...` - -CARA currently supports **only one geographic location for weather data per instance**. ## Running CARA locally @@ -85,10 +72,20 @@ This will start a local version of CARA, which can be visited at http://localhos ## Development guide +The CARA repository makes use of Git's Large File Storage (LFS) feature. +You will need a working installation of git-lfs in order to run CARA in development mode. +See https://git-lfs.github.com/ for installation instructions. + +### Installing CARA in editable mode + +``` +git lfs pull # Fetch the data from LFS +pip install -e . # At the root of the repository +``` + ### Running the COVID calculator app in development mode ``` -pip install -e . # At the root of the repository python -m cara.apps.calculator ``` @@ -107,7 +104,6 @@ python -m cara.apps.calculator --prefix=/mycalc ### Running the CARA Expert-App app in development mode ``` -pip install -e . # At the root of the repository voila cara/apps/expert/cara.ipynb --port=8080 ``` diff --git a/cara/.gitattributes b/cara/.gitattributes new file mode 100644 index 00000000..c1609315 --- /dev/null +++ b/cara/.gitattributes @@ -0,0 +1,2 @@ +global_weather_set.json filter=lfs diff=lfs merge=lfs -text +hadisd_station_fullinfo_v311_202001p.txt filter=lfs diff=lfs merge=lfs -text diff --git a/cara/apps/calculator/__init__.py b/cara/apps/calculator/__init__.py index dbfa15ef..1688d494 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__ = "2.1.0" +__version__ = "3.0.0" class BaseRequestHandler(RequestHandler): diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index a4b139e7..a6a41fab 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -1,5 +1,5 @@ import dataclasses -from dataclasses import dataclass +import datetime import html import logging import typing @@ -8,6 +8,7 @@ import numpy as np from cara import models 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 @@ -25,7 +26,7 @@ _NO_DEFAULT = object() _DEFAULT_MC_SAMPLE_SIZE = 50000 -@dataclass +@dataclasses.dataclass class FormData: activity_type: str air_changes: float @@ -50,6 +51,9 @@ class FormData: infected_lunch_start: minutes_since_midnight #Used if infected_dont_have_breaks_with_exposed infected_people: int infected_start: minutes_since_midnight + location_name: str + location_latitude: float + location_longitude: float mask_type: str mask_wearing_option: str mechanical_ventilation_type: str @@ -100,6 +104,9 @@ class FormData: 'infected_lunch_start': '12:30', 'infected_people': _NO_DEFAULT, 'infected_start': '08:30', + 'location_latitude': _NO_DEFAULT, + 'location_longitude': _NO_DEFAULT, + 'location_name': _NO_DEFAULT, 'mask_type': 'Type I', 'mask_wearing_option': 'mask_off', 'mechanical_ventilation_type': 'not-applicable', @@ -197,7 +204,8 @@ class FormData: ('virus_type', VIRUS_TYPES), ('volume_type', VOLUME_TYPES), ('window_opening_regime', WINDOWS_OPENING_REGIMES), - ('window_type', WINDOWS_TYPES)] + ('window_type', WINDOWS_TYPES), + ('event_month', MONTH_NAMES)] for attr_name, valid_set in validation_tuples: if getattr(self, attr_name) not in valid_set: raise ValueError(f"{getattr(self, attr_name)} is not a valid value for {attr_name}") @@ -244,6 +252,52 @@ class FormData: def build_model(self, sample_size=_DEFAULT_MC_SAMPLE_SIZE) -> models.ExposureModel: return self.build_mc_model().build_model(size=sample_size) + def tz_name_and_utc_offset(self) -> typing.Tuple[str, float]: + """ + Return the timezone name (e.g. CET), and offset, in hours, that need to + be *added* to UTC to convert to the form location's timezone. + + """ + month = MONTH_NAMES.index(self.event_month) + 1 + timezone = cara.data.weather.timezone_at( + latitude=self.location_latitude, longitude=self.location_longitude, + ) + # We choose the first of the month for the current year. + date = datetime.datetime(datetime.datetime.now().year, month, 1) + name = timezone.tzname(date) + assert isinstance(name, str) + utc_offset_td = timezone.utcoffset(date) + assert isinstance(utc_offset_td, datetime.timedelta) + utc_offset_hours = utc_offset_td.total_seconds() / 60 / 60 + return name, utc_offset_hours + + def outside_temp(self) -> models.PiecewiseConstant: + """ + Return the outside temperature as a PiecewiseConstant in the destination + timezone. + + """ + month = MONTH_NAMES.index(self.event_month) + 1 + + wx_station = self.nearest_weather_station() + temp_profile = cara.data.weather.mean_hourly_temperatures(wx_station[0], month) + + _, utc_offset = self.tz_name_and_utc_offset() + + # Offset the source times according to the difference from UTC (as a + # result the first data value may no longer be a midnight, and the hours + # no longer ordered modulo 24). + source_times = np.arange(24) + utc_offset + times, temp_profile = cara.data.weather.refine_hourly_data( + source_times, + temp_profile, + npts=24*10, # 10 steps per hour => 6 min steps + ) + outside_temp = models.PiecewiseConstant( + tuple(float(t) for t in times), tuple(float(t) for t in temp_profile), + ) + return outside_temp + def ventilation(self) -> models._VentilationBase: always_on = models.PeriodicInterval(period=120, duration=120) # Initializes a ventilation instance as a window if 'natural_ventilation' is selected, or as a HEPA-filter otherwise @@ -253,10 +307,8 @@ class FormData: else: window_interval = always_on - month = self.event_month[:3] - + outside_temp = self.outside_temp() inside_temp = models.PiecewiseConstant((0, 24), (293,)) - outside_temp = data.GenevaTemperatures[month] ventilation: models.Ventilation if self.window_type == 'window_sliding': @@ -298,6 +350,12 @@ class FormData: else: return models.MultipleVentilation((ventilation, infiltration_ventilation)) + def nearest_weather_station(self) -> cara.data.weather.WxStationRecordType: + """Return the nearest weather station (which has valid data) for this form""" + return cara.data.weather.nearest_wx_station( + longitude=self.location_longitude, latitude=self.location_latitude + ) + def mask(self) -> models.Mask: # Initializes the mask type if mask wearing is "continuous", otherwise instantiates the mask attribute as # the "No mask"-mask @@ -601,6 +659,9 @@ def baseline_raw_form_data(): 'infected_lunch_start': '12:30', 'infected_people': '1', 'infected_start': '09:00', + 'location_latitude': 46.20833, + 'location_longitude': 6.14275, + 'location_name': 'Geneva', 'mask_type': 'Type I', 'mask_wearing_option': 'mask_off', 'mechanical_ventilation_type': '', @@ -637,6 +698,11 @@ WINDOWS_TYPES = {'window_sliding', 'window_hinged', 'not-applicable'} COFFEE_OPTIONS_INT = {'coffee_break_0': 0, 'coffee_break_1': 1, 'coffee_break_2': 2, 'coffee_break_4': 4} +MONTH_NAMES = [ + 'January', 'February', 'March', 'April', 'May', 'June', 'July', + 'August', 'September', 'October', 'November', 'December', +] + def _hours2timestring(hours: float): # Convert times like 14.5 to strings, like "14:30" @@ -692,4 +758,3 @@ 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 - diff --git a/cara/apps/calculator/static/js/form.js b/cara/apps/calculator/static/js/form.js index b15f2fc6..abb888c5 100644 --- a/cara/apps/calculator/static/js/form.js +++ b/cara/apps/calculator/static/js/form.js @@ -311,7 +311,7 @@ function validate_form(form) { var lunch_finish = document.getElementById(activity+"_lunch_finish"); lunch_mins = parseTimeToMins(lunch_finish.value) - parseTimeToMins(lunch_start.value); } - + var coffee_breaks = parseInt(document.querySelector('input[name="'+activity+'_coffee_break_option"]:checked').value); var coffee_duration = parseInt(document.getElementById(activity+"_coffee_duration").value); var coffee_mins = coffee_breaks * coffee_duration; @@ -328,6 +328,24 @@ function validate_form(form) { }); } + // Validate location input. + if (submit) { + // We make the non-visible location inputs mandatory, without marking them as "required" inputs. + // See https://stackoverflow.com/q/22148080/741316 for motivation. + var locationSelectObj= document.getElementById("location_select"); + removeErrorFor(locationSelectObj); + $("input[name*='location']").each(function() { + el = $(this); + if ($.trim(el.val()) == ''){ + submit = false; + } + }); + + if (!submit) { + insertErrorFor(locationSelectObj, "Please select a location"); + } + } + //Validate all non zero values $("input[required].non_zero").each(function() { if (!validateValue(this)) { @@ -335,7 +353,6 @@ function validate_form(form) { } }); - //Validate window venting duration < venting frequency if (!$("#windows_duration").hasClass("disabled")) { var windowsDurationObj = document.getElementById("windows_duration"); @@ -464,9 +481,9 @@ function parseTimeToMins(cTime) { /* -------On Load------- */ $(document).ready(function () { - + var url = new URL(decodeURIComponent(window.location.href)); //Pre-fill form with known values - (new URL(decodeURIComponent(window.location.href))).searchParams.forEach((value, name) => { + url.searchParams.forEach((value, name) => { //If element exists if(document.getElementsByName(name).length > 0) { @@ -484,6 +501,7 @@ $(document).ready(function () { else if (elemObj.type === 'checkbox') { elemObj.checked = (value==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; @@ -492,6 +510,21 @@ $(document).ready(function () { } }); + // Handle default URL values if they are not explicitly defined. + if (Array.from(url.searchParams).length > 0) { + if (!url.searchParams.has('location_name')) { + $('[name="location_name"]').val('Geneva') + $('[name="location_select"]').val('Geneva') + } + if (!url.searchParams.has('location_latitude')) { + $('[name="location_latitude"]').val('46.20833') + } + if (!url.searchParams.has('location_longitude')) { + $('[name="location_longitude"]').val('6.14275') + } + } + + // When the document is ready, deal with the fact that we may be here // as a result of a forward/back browser action. If that is the case, update // the visibility of some of our inputs. @@ -531,8 +564,102 @@ $(document).ready(function () { $(".start_time[data-lunch-for]").each(function() {validateLunchBreak($(this).data('time-group'))}); $("[data-lunch-for]").change(function() {validateLunchBreak($(this).data('time-group'))}); $("[data-lunch-break]").change(function() {validateLunchBreak($(this).data('lunch-break'))}); + + $("#location_select").select2({ + ajax: { + // Docs for the geocoding service at: + // https://developers.arcgis.com/rest/geocode/api-reference/geocoding-service-output.htm + url: "https://geocode.arcgis.com/arcgis/rest/services/World/GeocodeServer/suggest", + dataType: 'json', + delay: 250, + data: function(params) { + return { + text: params.term, // search term + f: 'json', + page: params.page, + maxSuggestions: 20, + }; + }, + processResults: function(data, params) { + // Enable infinite scrolling + params.page = params.page || 1; + return { + results: data.suggestions.map(function(suggestion) { + return { + id: suggestion.magicKey, // The unique reference to this result. + text: suggestion.text, + magicKey: suggestion.magicKey + } + }), + pagination: { + more: (params.page * 10) < data.suggestions.length + } + }; + }, + cache: true + }, + placeholder: 'Search for a location', + minimumInputLength: 1, + templateResult: formatlocation, + templateSelection: formatLocationSelection + }); + + function formatlocation(suggestedLocation) { + // Function is called for each location from the geocoding API. + + if (suggestedLocation.loading) { + // Update the first message in the search results to show the + // "Searching..." message. + return suggestedLocation.text; + } + + // Create a container for this location (to be added to the DOM by the select2 + // library when returned). + // This will become one of many search results in the dropdown. + var $container = $( + "
" + + "
" + + "
" + suggestedLocation.text + "
" + + "
" + + "
" + ); + return $container; + } + + function formatLocationSelection(selectedSuggestion) { + // Function is called when a selection is made in the search result dropdown. + + // ID may be empty, for example when the page is refreshed or back button pressed. + if (selectedSuggestion.id != "") { + + // Turn the suggestion into a proper location (so that we can get its latitude & longitude). + $.ajax({ + dataType: "json", + url: 'https://geocode.arcgis.com/arcgis/rest/services/World/GeocodeServer/findAddressCandidates', + data: { + magicKey: selectedSuggestion.magicKey, + outFields: 'country, location', + f: "json" + }, + success: function (locations) { + // If there isn't precisely one result something is very wrong. + geocoded_loc = locations.candidates[0]; + $('input[name="location_name"]').val(selectedSuggestion.text); + $('input[name="location_latitude"]').val(geocoded_loc.location.y.toPrecision(7)); + $('input[name="location_longitude"]').val(geocoded_loc.location.x.toPrecision(7)); + } + }); + + } else if ($('input[name="location_name"]').val() != "") { + // If we have no selection AND the location_name is available, use that in the search bar. + // This means that we preserve the location through refresh/back button. + return $('input[name="location_name"]').val(); + } + return selectedSuggestion.text; + } }); + /* -------Debugging------- */ function debug_submit(form) { diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index 8a36f52b..ad04ffae 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -199,6 +199,10 @@

  • Room Volume: {{ model.concentration_model.room.volume }} m³

  • Room Central Heating: {{ "On" if form.room_heating_option else "Off" }}

    • +

  • Geographic Location: {{ form.location_name }}

    • + {% if form.ventilation_type == "natural_ventilation" %} +

  • Nearest weather station: {{ form.nearest_weather_station()[1].strip().title() }}

    • + {% endif %} @@ -245,7 +249,7 @@ {% endif %}

    -

    When using the natural ventilation option, air flows are calculated using averaged hourly temperatures for the Geneva region, based on historical data for the month selected.

    +

    When using the natural ventilation option, air flows are calculated using averaged hourly temperatures for the region {{ form.location_name }}, based on historical data for the month selected.

    {% else %} No

    {% endif %} diff --git a/cara/apps/calculator/templates/calculator.form.html.j2 b/cara/apps/calculator/templates/calculator.form.html.j2 index 712432fb..c119d4e5 100644 --- a/cara/apps/calculator/templates/calculator.form.html.j2 +++ b/cara/apps/calculator/templates/calculator.form.html.j2 @@ -6,11 +6,13 @@ {% block extra_headers %} + {% endblock extra_headers %} {% block body_scripts %} + {% endblock body_scripts %} @@ -102,6 +104,23 @@ v{{ calculator_version }} Please sen       + +
    + +
    + ? +
    + +
    + + + + +
    +
    @@ -206,10 +225,10 @@ v{{ calculator_version }} Please sen
    - Exposed person(s) presence:
    + Exposed person(s) presence:
    Start:    Finish:
    - Infected person(s) presence:
    + Infected person(s) presence:
    Start:    Finish:
    diff --git a/cara/data.py b/cara/data/__init__.py similarity index 96% rename from cara/data.py rename to cara/data/__init__.py index 8e844dae..353992e6 100644 --- a/cara/data.py +++ b/cara/data/__init__.py @@ -1,6 +1,8 @@ import numpy as np from cara import models +# TODO: The values in this module to be removed and instead use the cara.data.weather functionality. + # average temperature of each month, hour per hour (from midnight to 11 pm) Geneva_hourly_temperatures_celsius_per_hour = { 'Jan': [0.2, -0.3, -0.5, -0.9, -1.1, -1.4, -1.5, -1.5, -1.1, 0.1, 1.5, diff --git a/cara/data/global_weather_set.json b/cara/data/global_weather_set.json new file mode 100644 index 00000000..ba1765ee --- /dev/null +++ b/cara/data/global_weather_set.json @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c6b41b08c350c543bced4d15b670851042cf1ca9135551b2ed2afb2d99ec63e8 +size 13803155 diff --git a/cara/data/hadisd_station_fullinfo_v311_202001p.txt b/cara/data/hadisd_station_fullinfo_v311_202001p.txt new file mode 100644 index 00000000..1747636a --- /dev/null +++ b/cara/data/hadisd_station_fullinfo_v311_202001p.txt @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4843d34b6e4c26d4382860e011451d5f32157b9a3660830f8d2894a11d022298 +size 772370 diff --git a/cara/data/weather.py b/cara/data/weather.py new file mode 100644 index 00000000..af5f27ec --- /dev/null +++ b/cara/data/weather.py @@ -0,0 +1,149 @@ +import datetime +import functools +import json +from pathlib import Path +import typing + +import dateutil.tz +import numpy as np +from scipy.spatial import cKDTree +from timezonefinder import TimezoneFinder + + +WX_DATA_LOCATION = Path(__file__).absolute().parent +WxStationIdType = str +MonthType = str +# HourlyTempType - 24 temperatures, one for each hour of the day (the average for the given month). +HourlyTempType = typing.List[float] +WxStationRecordType = typing.Tuple[WxStationIdType, str, float, float] + + +@functools.lru_cache() +def wx_data() -> typing.Dict[WxStationIdType, typing.Dict[MonthType, HourlyTempType]]: + """ + Load the weather data (temperature in kelvin). + + The data is structured by station location, and for each station location, by month. + + """ + with (WX_DATA_LOCATION / 'global_weather_set.json').open("r") as json_file: + data = json.load(json_file) + + for station in list(data.keys()): + for month in list(data[station].keys()): + if not np.any(np.isnan(data[station][month])): + data[station][month] = tuple( + 273.15 + np.array(data[station][month])) + return data + + +@functools.lru_cache() +def wx_station_data() -> typing.Dict[WxStationIdType, WxStationRecordType]: + """ + Return a dictionary of ``station-id: station records``, where station records + are of the form ``(station-id, station-name, station-latitude, station-longitude)``. + + The stations returned are guaranteed to have valid weather data. + + """ + weather_data = wx_data() + station_data = {} + fixed_delimits = [0, 12, 13, 44, 51, 60, 69, 90, 91] + station_file = WX_DATA_LOCATION / 'hadisd_station_fullinfo_v311_202001p.txt' + + for line in station_file.open('rt'): + start_end_positions = zip(fixed_delimits[:-1], fixed_delimits[1:]) + split_vals = [line[start:end] for start, end in start_end_positions] + station_location = ( + split_vals[0], split_vals[2], float(split_vals[3]), float(split_vals[4]), + ) + # We only consider stations with weather data, don't include the rest. + if split_vals[0] in weather_data: + station_data[split_vals[0]] = station_location + return station_data + + +@functools.lru_cache() +def _wx_station_kdtree() -> cKDTree: + """Build a kd-tree of wx station longitude & latitudes (note the coordinate order)""" + station_data = wx_station_data().values() + coords = np.array([(stn_record[3], stn_record[2]) + for stn_record in station_data]) + return cKDTree(coords) + + +def mean_hourly_temperatures(wx_station: str, month: int) -> HourlyTempType: + """ + Return the mean monthly temperature for the given weather station and month. + + Returns + ------- + + temperatures: List[24 floats] + A list containing 24 temperature values, one for each hour, in kelvin. + Index 0 of the result corresponds to hour 00:00 (UTC), and index 23 (the last) to 23:00 (UTC). + + """ + # Note that the current dataset encodes month number as a string. + return wx_data()[wx_station][str(month)] + + +def timezone_at(*, latitude: float, longitude: float) -> datetime.tzinfo: + """Find a timezone for the given location, or raise.""" + tf = TimezoneFinder() + tz_name = tf.timezone_at(lat=latitude, lng=longitude) + tz = dateutil.tz.gettz(tz_name) + if tz_name is None or tz is None: + raise ValueError( + f"Unable to determine the timezone of given location " + f"(lat={latitude}, lng={longitude})" + ) + return tz + + +def refine_hourly_data(source_times, hourly_data, npts): + """ + Given times (in hours), where each data point is on the hour, + interpolate the data to mid-point of the returned boundaries. + + For example: + + >>> time_bounds, data = refine_hourly_data(list(range(24)), list(range(24)), 24) + >>> len(time_bounds), len(data) + (25, 24) + >>> time_bounds + array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., + 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24.]) + >>> data + array([ 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, + 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, + 22.5, 11.5]) + + The source times need not be monotonic, which allows for data to be + time-offset shifted. For example: + + >>> time_bounds, data = refine_hourly_data( + ... list(range(6, 28)) + [4, 5], list(range(24)), 24) + >>> data + array([18.5, 19.5, 20.5, 21.5, 22.5, 11.5, 0.5, 1.5, 2.5, 3.5, 4.5, + 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, + 16.5, 17.5]) + + """ + target_time_boundaries, step = np.linspace( + 0, 24, npts + 1, retstep=True, endpoint=True, + ) + target_times = target_time_boundaries[:-1] + step / 2 + data = np.interp(target_times, np.array(source_times), hourly_data, period=24) + return target_time_boundaries, data + + +def nearest_wx_station(*, longitude: float, latitude: float) -> WxStationRecordType: + """ + Given a latitude & longitude, return the nearest station with valid weather data. + + """ + ktree = _wx_station_kdtree() + station_data = list(wx_station_data().values()) + dd, ii = ktree.query((longitude, latitude), k=[1]) + return station_data[ii[0]] diff --git a/cara/tests/apps/calculator/test_model_generator.py b/cara/tests/apps/calculator/test_model_generator.py index 53045564..6ee821c0 100644 --- a/cara/tests/apps/calculator/test_model_generator.py +++ b/cara/tests/apps/calculator/test_model_generator.py @@ -1,14 +1,14 @@ import dataclasses +import typing +import numpy as np +import numpy.testing as npt import pytest from cara.apps.calculator import model_generator from cara.apps.calculator.model_generator import _hours2timestring from cara.apps.calculator.model_generator import minutes_since_midnight from cara import models -from cara import data -import numpy as np -import numpy.testing as npt def test_model_from_dict(baseline_form_data): @@ -33,13 +33,6 @@ def test_blend_expiration(): def test_ventilation_slidingwindow(baseline_form: model_generator.FormData): - room = models.Room(75) - window = models.SlidingWindow( - active=models.PeriodicInterval(period=120, duration=10), - inside_temp=models.PiecewiseConstant((0, 24), (293,)), - outside_temp=data.GenevaTemperatures['Dec'], - window_height=1.6, opening_length=0.6, - ) baseline_form.ventilation_type = 'natural_ventilation' baseline_form.windows_duration = 10 baseline_form.windows_frequency = 120 @@ -49,19 +42,28 @@ def test_ventilation_slidingwindow(baseline_form: model_generator.FormData): baseline_form.window_height = 1.6 baseline_form.opening_distance = 0.6 - ts = np.linspace(8, 16, 100) - np.testing.assert_allclose([window.air_exchange(room, t)+0.25 for t in ts], - [baseline_form.ventilation().air_exchange(room, t) for t in ts]) + baseline_vent = baseline_form.ventilation() + assert isinstance(baseline_vent, models.MultipleVentilation) + baseline_window = baseline_vent.ventilations[0] + assert isinstance(baseline_window, models.SlidingWindow) + + window = models.SlidingWindow( + active=models.PeriodicInterval(period=120, duration=10), + inside_temp=models.PiecewiseConstant((0, 24), (293,)), + outside_temp=baseline_window.outside_temp, + window_height=1.6, opening_length=0.6, + ) + + ach = models.AirChange( + active=models.PeriodicInterval(period=120, duration=120), + air_exch=0.25, + ) + ventilation = models.MultipleVentilation((window, ach)) + + assert ventilation == baseline_vent def test_ventilation_hingedwindow(baseline_form: model_generator.FormData): - room = models.Room(75) - window = models.HingedWindow( - active=models.PeriodicInterval(period=120, duration=10), - inside_temp=models.PiecewiseConstant((0, 24), (293,)), - outside_temp=data.GenevaTemperatures['Dec'], - window_height=1.6, window_width=1., opening_length=0.6, - ) baseline_form.ventilation_type = 'natural_ventilation' baseline_form.windows_duration = 10 baseline_form.windows_frequency = 120 @@ -72,9 +74,24 @@ def test_ventilation_hingedwindow(baseline_form: model_generator.FormData): baseline_form.window_width = 1. baseline_form.opening_distance = 0.6 - ts = np.linspace(8, 16, 100) - np.testing.assert_allclose([window.air_exchange(room, t)+0.25 for t in ts], - [baseline_form.ventilation().air_exchange(room, t) for t in ts]) + baseline_vent = baseline_form.ventilation() + assert isinstance(baseline_vent, models.MultipleVentilation) + baseline_window = baseline_vent.ventilations[0] + assert isinstance(baseline_window, models.HingedWindow) + + window = models.HingedWindow( + active=models.PeriodicInterval(period=120, duration=10), + inside_temp=models.PiecewiseConstant((0, 24), (293,)), + outside_temp=baseline_window.outside_temp, + window_height=1.6, window_width=1., opening_length=0.6, + ) + ach = models.AirChange( + active=models.PeriodicInterval(period=120, duration=120), + air_exch=0.25, + ) + ventilation = models.MultipleVentilation((window, ach)) + + assert ventilation == baseline_vent def test_ventilation_mechanical(baseline_form: model_generator.FormData): @@ -108,19 +125,6 @@ def test_ventilation_airchanges(baseline_form: model_generator.FormData): def test_ventilation_window_hepa(baseline_form: model_generator.FormData): - room = models.Room(75) - window = models.SlidingWindow( - active=models.PeriodicInterval(period=120, duration=10), - inside_temp=models.PiecewiseConstant((0, 24), (293,)), - outside_temp=data.GenevaTemperatures['Dec'], - window_height=1.6, opening_length=0.6, - ) - hepa = models.HEPAFilter( - active=models.PeriodicInterval(period=120, duration=120), - q_air_mech=250., - ) - ventilation = models.MultipleVentilation((window,hepa)) - baseline_form.ventilation_type = 'natural_ventilation' baseline_form.windows_duration = 10 baseline_form.windows_frequency = 120 @@ -130,9 +134,29 @@ def test_ventilation_window_hepa(baseline_form: model_generator.FormData): baseline_form.opening_distance = 0.6 baseline_form.hepa_option = True - ts = np.linspace(9, 17, 100) - np.testing.assert_allclose([ventilation.air_exchange(room, t)+0.25 for t in ts], - [baseline_form.ventilation().air_exchange(room, t) for t in ts]) + baseline_vent = baseline_form.ventilation() + assert isinstance(baseline_vent, models.MultipleVentilation) + baseline_window = baseline_vent.ventilations[0] + assert isinstance(baseline_window, models.SlidingWindow) + + # Now build the equivalent ventilation instance directly, and compare. + window = models.SlidingWindow( + active=models.PeriodicInterval(period=120, duration=10), + inside_temp=models.PiecewiseConstant((0, 24), (293,)), + outside_temp=baseline_window.outside_temp, + window_height=1.6, opening_length=0.6, + ) + hepa = models.HEPAFilter( + active=models.PeriodicInterval(period=120, duration=120), + q_air_mech=250., + ) + ach = models.AirChange( + active=models.PeriodicInterval(period=120, duration=120), + air_exch=0.25, + ) + ventilation = models.MultipleVentilation((window, hepa, ach)) + + assert ventilation == baseline_vent def present_times(interval: models.Interval) -> models.BoundarySequence_t: @@ -410,6 +434,12 @@ def test_key_validation_mech_ventilation_type_na(baseline_form_data): model_generator.FormData.from_dict(baseline_form_data) +def test_key_validation_event_month(baseline_form_data): + baseline_form_data['event_month'] = 'invalid month' + with pytest.raises(ValueError, match='invalid month is not a valid value for event_month'): + model_generator.FormData.from_dict(baseline_form_data) + + def test_default_types(): # Validate that FormData._DEFAULTS are complete and of the correct type. # Validate that we have the right types and matching attributes to the DEFAULTS. @@ -444,3 +474,23 @@ def test_form_to_dict(baseline_form): # If we set the value to the default one, it should no longer turn up in the dictionary. baseline_form.exposed_coffee_break_option = model_generator.FormData._DEFAULTS['exposed_coffee_break_option'] assert 'exposed_coffee_break_option' not in baseline_form.to_dict(baseline_form, strip_defaults=True) + + +@pytest.mark.parametrize( + ["longitude", "latitude", "month", "expected_tz_name", "expected_offset"], + [ + [6.14275, 46.20833, "January", 'CET', 1], # Geneva, winter + [6.14275, 46.20833, "May", 'CEST', 2], # Geneva, summer + [144.96751, -37.81739, "January", 'AEDT', 11], # Melbourne, summer time + [144.96751, -37.81739, "June", 'AEST', 10], # Melbourne, winter time + [-176.433333, -44.033333, 'August', '+1245', 12.75], # Chatham Islands + ] +) +def test_form_timezone(baseline_form_data, longitude, latitude, month, expected_tz_name, expected_offset): + baseline_form_data['location_latitude'] = latitude + baseline_form_data['location_longitude'] = longitude + baseline_form_data['event_month'] = month + form = model_generator.FormData.from_dict(baseline_form_data) + name, offset = form.tz_name_and_utc_offset() + assert name == expected_tz_name + assert offset == expected_offset diff --git a/cara/tests/data/__init__.py b/cara/tests/data/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/cara/tests/data/test_weather.py b/cara/tests/data/test_weather.py new file mode 100644 index 00000000..2a667df8 --- /dev/null +++ b/cara/tests/data/test_weather.py @@ -0,0 +1,74 @@ +import datetime + +import dateutil.tz +import numpy as np +import numpy.testing +import pytest + +import cara.data.weather as wx + + +def test_nearest_wx_station(): + melbourne_lat, melbourne_lon = -37.81739, 144.96751 + station_rec = wx.nearest_wx_station(longitude=melbourne_lon, latitude=melbourne_lat) + + station_name = station_rec[1].strip() + # Note: For Melbourne, the nearest station is 'MELBOURNE REGIONAL OFFICE', + # but the nearest location with suitable wx data is 'MELBOURNE ESSENDON' + assert station_name == 'MELBOURNE ESSENDON' + + +def test_refine(): + source_times = [0, 3, 6, 9, 12, 15, 18, 21] + data = [0, 30, 60, 90, 120, 90, 60, 30] + + time_bounds, data = wx.refine_hourly_data(source_times, data, 4) + + # Notice that the expected data falls in the mid-point of the + # expected time bounds. + np.testing.assert_array_equal(time_bounds, [0., 6., 12., 18., 24.]) + np.testing.assert_array_equal(data, [30., 90., 90., 30.]) + + +def test_refine_offset(): + source_times = [14, 20, 26, 32] + data = [200., 182, 168, 192] + + time_bounds, data = wx.refine_hourly_data(source_times, data, 6) + + # Notice that the expected data falls in the mid-point of the + # expected time bounds. + np.testing.assert_array_equal(time_bounds, [0., 4., 8., 12., 16., 20., 24.]) + np.testing.assert_array_almost_equal(data, [168., 184., 194.666667, 200., 188., 177.333333]) + + +def test_refine_non_monotonic(): + source_times = [14, 20, 2, 8] + data = [200., 182, 168, 192] + + time_bounds, data = wx.refine_hourly_data(source_times, data, 6) + + # Notice that the expected data falls in the mid-point of the + # expected time bounds. + np.testing.assert_array_equal(time_bounds, [0., 4., 8., 12., 16., 20., 24.]) + np.testing.assert_array_almost_equal(data, [168., 184., 194.666667, 200., 188., 177.333333]) + + + +def test_timezone_at__out_of_range(): + with pytest.raises(ValueError, match='out of bounds'): + wx.timezone_at(latitude=88, longitude=181) + + +@pytest.mark.parametrize( + ["latitude", "longitude", "expected_tz_name"], + [ + [6.14275, 46.20833, 'Europe/Zurich'], # Geneva + [144.96751, -37.81739, "Australia/Melbourne"], # Melbourne + [-176.433333, -44.033333, 'Pacific/Chatham'], # Chatham Islands + ] +) +def test_timezone_at__expected(latitude, longitude, expected_tz_name): + assert wx.timezone_at(latitude=longitude, longitude=latitude) == dateutil.tz.gettz(expected_tz_name) + assert wx.timezone_at(latitude=0, longitude=-175) == dateutil.tz.gettz('Etc/GMT+12') + assert wx.timezone_at(latitude=89.8, longitude=-170) == dateutil.tz.gettz('Etc/GMT+11') diff --git a/requirements.txt b/requirements.txt index 302de847..457a69b7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -74,6 +74,7 @@ sniffio==1.2.0 terminado==0.10.1 testpath==0.5.0 threadpoolctl==2.2.0 +timezonefinder==5.2.0 tornado==6.1 traitlets==5.0.5 urllib3==1.26.6 diff --git a/setup.cfg b/setup.cfg index 21bbf8d4..c2f968b6 100644 --- a/setup.cfg +++ b/setup.cfg @@ -25,3 +25,7 @@ ignore_missing_imports = True [mypy-scipy.*] ignore_missing_imports = True + +[mypy-timezonefinder.*] +ignore_missing_imports = True + diff --git a/setup.py b/setup.py index a70d5237..9e5eb596 100644 --- a/setup.py +++ b/setup.py @@ -29,9 +29,11 @@ REQUIREMENTS: dict = { 'mistune', 'numpy', 'psutil', + 'python-dateutil', 'qrcode[pil]', 'scipy', 'sklearn', + 'timezonefinder', 'tornado', 'voila >=0.2.4', ], @@ -42,6 +44,7 @@ REQUIREMENTS: dict = { 'pytest-tornasync', 'numpy-stubs @ git+https://github.com/numpy/numpy-stubs.git', 'types-dataclasses', + 'types-python-dateutil', ], 'dev': [ 'jupyterlab', @@ -84,5 +87,7 @@ setup( 'apps/*/*/*', 'apps/*/*/*/*', 'apps/*/*/*/*/*', + 'data/*.json', + 'data/*.txt', ]}, ) From bdc51d81b8e3315010e7bd9c136ff57957bd2edd Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Thu, 2 Sep 2021 12:11:22 +0200 Subject: [PATCH 3/5] Changed relevant references for mask distributions --- cara/monte_carlo/data.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 250bd6e6..4415efbb 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -60,7 +60,10 @@ virus_distributions = { } -# From CERN-OPEN-2021-04 and refererences therein +# From: +# https://doi.org/10.1080/02786826.2021.1890687 +# https://doi.org/10.1016/j.jhin.2013.02.007 +# https://doi.org/10.4209/aaqr.2020.08.0531 mask_distributions = { 'Type I': mc.Mask(Uniform(0.25, 0.80)), 'FFP2': mc.Mask(Uniform(0.83, 0.91)), From ab153ca7812fdceda7befb5a5de0bd56663c8f9d Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Thu, 2 Sep 2021 15:23:49 +0200 Subject: [PATCH 4/5] added inward mask efficiency references --- cara/apps/templates/common_text.md.j2 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/cara/apps/templates/common_text.md.j2 b/cara/apps/templates/common_text.md.j2 index 63512dee..ac8829d7 100644 --- a/cara/apps/templates/common_text.md.j2 +++ b/cara/apps/templates/common_text.md.j2 @@ -71,3 +71,6 @@ We wish to thank CERN’s HSE Unit, Beams Department, Experimental Physics Depar [54] Leung, N.H.L et al. Respiratory virus shedding in exhaled breath and efficacy of face masks. Nat Med (2020). 10.1038/s41591-020-0843-2.
    [55] Asadi, S., Cappa, C.D., Barreda, S. et al. Efficacy of masks and face coverings in controlling outward aerosol particle emission from expiratory activities. Sci Rep 10, 15665 (2020). https://doi.org/10.1038/s41598-020-72798-7.
    [56] Endo A, Abbott S et al. Estimating the overdispersion in COVID-19 transmission using outbreak sizes outside China [version 3; peer review: 2 approved]. Wellcome Open Res 2020, 5:67. doi:10.12688/wellcomeopenres.15842.3.
    +[57] Jin Pan, Charbel Harb, Weinan Leng & Linsey C. Marr (2021) Inward and outward effectiveness of cloth masks, a surgical mask, and a face shield, Aerosol Science and Technology, 55:6, 718-733, doi: 10.1080/02786826.2021.1890687.
    +[58] C. Makison Booth, M. Clayton, B. Crook, J.M. Gawn, Effectiveness of surgical masks against influenza bioaerosols, Journal of Hospital Infection, Volume 84, Issue 1, 2013, Pages 22-26, https://doi.org/10.1016/j.jhin.2013.02.007.
    +[59] Riediker, M., Monn, C. (2021). Simulation of SARS-CoV-2 Aerosol Emissions in the Infected Population and Resulting Airborne Exposures in Different Indoor Scenarios. Aerosol Air Qual. Res. 21, 200531. https://doi.org/10.4209/aaqr.2020.08.0531.
    \ No newline at end of file From 250b7c0582ecedf58b3e90eb4d7b3f752dcbdd04 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 7 Sep 2021 06:43:22 +0200 Subject: [PATCH 5/5] Plot alternative scenarios in report using D3.js, speeding up the response time and improving user interactivity --- cara/apps/calculator/report_generator.py | 77 +----- cara/apps/calculator/static/js/report.js | 226 ++++++++++++++---- .../templates/base/calculator.report.html.j2 | 14 +- 3 files changed, 192 insertions(+), 125 deletions(-) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index d1bff9dc..773481f3 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -9,11 +9,9 @@ import zlib import loky import jinja2 -import matplotlib -matplotlib.use('agg') -import matplotlib.pyplot as plt import numpy as np import qrcode +import json from cara import models from ... import monte_carlo as mc @@ -162,50 +160,13 @@ def _img2bytes(figure): return img_data -def _figure2bytes(figure): - # Draw the image - img_data = io.BytesIO() - figure.savefig(img_data, format='png', bbox_inches="tight", transparent=True) - return img_data - - def img2base64(img_data) -> str: - plt.close() img_data.seek(0) pic_hash = base64.b64encode(img_data.read()).decode('ascii') # A src suitable for a tag such as f'. return f'data:image/png;base64,{pic_hash}' -def plot(times, concentrations, model: models.ExposureModel): - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1) - datetimes = [datetime(1970, 1, 1) + timedelta(hours=time) for time in times] - ax.plot(datetimes, concentrations, lw=2, color='#1f77b4', label='Mean concentration') - ax.spines['right'].set_visible(False) - ax.spines['top'].set_visible(False) - - ax.set_xlabel('Time of day') - ax.set_ylabel('Mean concentration ($virions/m^{3}$)') - ax.set_title('Mean concentration of virions') - ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) - - # Plot presence of exposed person - for i, (presence_start, presence_finish) in enumerate(model.exposed.presence.boundaries()): - plt.fill_between( - datetimes, concentrations, 0, - where=(np.array(times) > presence_start) & (np.array(times) < presence_finish), - color="#1f77b4", alpha=0.1, - label="Presence of exposed person(s)" if i == 0 else "" - ) - - # Place a legend outside of the axes itself. - ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') - ax.set_ylim(0) - - return fig - - def minutes_to_time(minutes: int) -> str: minute_string = str(minutes % 60) minute_string = "0" * (2 - len(minute_string)) + minute_string @@ -281,39 +242,7 @@ def manufacture_alternative_scenarios(form: FormData) -> typing.Dict[str, mc.Exp return scenarios -def comparison_plot(scenarios: typing.Dict[str, dict], sample_times: typing.List[float]): - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1) - - dash_styled_scenarios = [ - 'Base scenario with FFP2 masks', - 'Base scenario with HEPA filter', - 'Base scenario with HEPA and FFP2 masks', - ] - - sample_dts = [datetime(1970, 1, 1) + timedelta(hours=time) for time in sample_times] - for name, statistics in scenarios.items(): - concentrations = statistics['concentrations'] - - if name in dash_styled_scenarios: - ax.plot(sample_dts, concentrations, label=name, linestyle='--') - else: - ax.plot(sample_dts, concentrations, label=name, linestyle='-', alpha=0.5) - - # Place a legend outside of the axes itself. - ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') - ax.spines['right'].set_visible(False) - ax.spines['top'].set_visible(False) - ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter("%H:%M")) - - ax.set_xlabel('Time of day') - ax.set_ylabel('Mean concentration ($virions/m^{3}$)') - ax.set_title('Mean concentration of virions') - - return fig - - -def scenario_statistics(mc_model: mc.ExposureModel, sample_times: typing.List[float]): +def scenario_statistics(mc_model: mc.ExposureModel, sample_times: np.ndarray): model = mc_model.build_model(size=_DEFAULT_MC_SAMPLE_SIZE) return { 'probability_of_infection': np.mean(model.infection_probability()), @@ -342,7 +271,6 @@ def comparison_report( for (name, model), model_stats in zip(scenarios.items(), results): statistics[name] = model_stats return { - 'plot': img2base64(_figure2bytes(comparison_plot(statistics, sample_times))), 'stats': statistics, } @@ -405,6 +333,7 @@ class ReportGenerator: env.filters['minutes_to_time'] = minutes_to_time env.filters['float_format'] = "{0:.2f}".format env.filters['int_format'] = "{:0.0f}".format + env.filters['JSONify'] = json.dumps return env def render(self, context: dict) -> str: diff --git a/cara/apps/calculator/static/js/report.js b/cara/apps/calculator/static/js/report.js index d1105bec..b815da2f 100644 --- a/cara/apps/calculator/static/js/report.js +++ b/cara/apps/calculator/static/js/report.js @@ -25,56 +25,16 @@ function draw_concentration_plot(svg_id, times, concentrations, exposed_presence yAxis = d3.axisLeft(yRange); // Plot tittle. - vis.append('svg:foreignObject') - .attr("background-color", "transparent") - .attr('width', width) - .attr('height', margins.top) - .style('text-align', 'center') - .html('Mean concentration of virions'); + plot_title(vis, width, margins.top, 'Mean concentration of virions'); // Line representing the mean concentration. - var lineFunc = d3.line() - .defined(d => !isNaN(d.concentration)) - .x(d => xTimeRange(d.time)) - .y(d => yRange(d.concentration)) - .curve(d3.curveBasis); + plot_scenario_data(vis, data, xTimeRange, yRange, '#1f77b4'); - vis.append('svg:path') - .attr('d', lineFunc(data)) - .attr('stroke', '#1f77b4') - .attr('stroke-width', 2) - .attr('fill', 'none'); + // X axis. + plot_x_axis(vis, height, width, margins, xAxis, 'Time of day'); - // X axis declaration. - vis.append('svg:g') - .attr('class', 'x axis') - .attr('transform', 'translate(0,' + (height - margins.bottom) + ')') - .call(xAxis); - - // X axis label. - vis.append('text') - .attr('class', 'x label') - .attr('fill', 'black') - .attr('text-anchor', 'middle') - .attr('x', (width + margins.right) / 2) - .attr('y', height * 0.97) - .text('Time of day') - - // Y axis declaration. - vis.append('svg:g') - .attr('class', 'y axis') - .attr('transform', 'translate(' + margins.left + ',0)') - .call(yAxis); - - // Y axis label. - vis.append('svg:text') - .attr('class', 'y label') - .attr('fill', 'black') - .attr('transform', 'rotate(-90, 0,' + height + ')') - .attr('text-anchor', 'middle') - .attr('x', (height + margins.bottom) / 2) - .attr('y', (height + margins.left) * 0.92) - .text('Mean concentration (virions/m³)'); + // Y axis + plot_y_axis(vis, height, width, margins, yAxis, 'Mean concentration (virions/m³)') // Area representing the presence of exposed person(s). exposed_presence_intervals.forEach(b => { @@ -181,4 +141,178 @@ function draw_concentration_plot(svg_id, times, concentrations, exposed_presence focus.select('#tooltip-time').text('x = ' + time_format(d.hour)); focus.select('#tooltip-concentration').text('y = ' + d.concentration.toFixed(2)); } +} + +// Generate the alternative scenarios plot using d3 library. +// 'alternative_scenarios' is a dictionary with all the alternative scenarios +// 'times' is a list of times for all the scenarios +function draw_alternative_scenarios_plot(svg_id, width, height, alternative_scenarios, times) { + // H:M format + var time_format = d3.timeFormat('%H:%M'); + // D3 array of ten categorical colors represented as RGB hexadecimal strings. + var colors = d3.schemeAccent; + + // Variable for the highest concentration for all the scenarios + var highest_concentration = 0. + + var data_for_scenarios = {} + for (scenario in alternative_scenarios) { + scenario_concentrations = alternative_scenarios[scenario].concentrations + + highest_concentration = Math.max(highest_concentration, Math.max(...scenario_concentrations)) + + var data = [] + times.map((time, index) => data.push({ 'time': time, 'hour': new Date().setHours(Math.trunc(time), (time - Math.trunc(time)) * 60), 'concentration': scenario_concentrations[index] })) + + // Add data into lines dictionary + data_for_scenarios[scenario] = data + } + + // We need one scenario to get the time range + var first_scenario = Object.values(data_for_scenarios)[0] + + var vis = d3.select(svg_id), + width = width, + height = height, + margins = { top: 30, right: 20, bottom: 50, left: 50 }, + + // H:M time format for x axis. + xRange = d3.scaleTime().range([margins.left, width - margins.right]).domain([first_scenario[0].hour, first_scenario[first_scenario.length - 1].hour]), + xTimeRange = d3.scaleLinear().range([margins.left, width - margins.right]).domain([times[0], times[times.length - 1]]), + + yRange = d3.scaleLinear().range([height - margins.bottom, margins.top]).domain([0., highest_concentration]), + + xAxis = d3.axisBottom(xRange).tickFormat(d => time_format(d)), + yAxis = d3.axisLeft(yRange); + + // Plot title. + plot_title(vis, width, margins.top, 'Mean concentration of virions'); + + // Line representing the mean concentration for each scenario. + for (const [scenario_name, data] of Object.entries(data_for_scenarios)) { + var scenario_index = Object.keys(data_for_scenarios).indexOf(scenario_name) + + // Line representing the mean concentration. + plot_scenario_data(vis, data, xTimeRange, yRange, colors[scenario_index]) + + // Legend for the plot elements - lines. + var size = 20 * (scenario_index + 1) + vis.append('rect') + .attr('x', width + 20) + .attr('y', margins.top + size) + .attr('width', 20) + .attr('height', 3) + .style('fill', colors[scenario_index]); + + vis.append('text') + .attr('x', width + 3 * 20) + .attr('y', margins.top + size) + .text(scenario_name) + .style('font-size', '15px') + .attr('alignment-baseline', 'central'); + + } + + // X axis. + plot_x_axis(vis, height, width, margins, xAxis, "Time of day"); + + // Y axis declaration. + vis.append('svg:g') + .attr('class', 'y axis') + .attr('transform', 'translate(' + margins.left + ',0)') + .call(yAxis); + + // Y axis label. + vis.append('svg:text') + .attr('class', 'y label') + .attr('fill', 'black') + .attr('transform', 'rotate(-90, 0,' + height + ')') + .attr('text-anchor', 'middle') + .attr('x', (height + margins.bottom) / 2) + .attr('y', (height + margins.left) * 0.92) + .text('Mean concentration (virions/m³)'); + + // Legend bounding box. + vis.append('rect') + .attr('width', 275) + .attr('height', 25 * (Object.keys(data_for_scenarios).length)) + .attr('x', width * 1.005) + .attr('y', margins.top + 5) + .attr('stroke', 'lightgrey') + .attr('stroke-width', '2') + .attr('rx', '5px') + .attr('ry', '5px') + .attr('stroke-linejoin', 'round') + .attr('fill', 'none'); +} + + +// Functions used to build the plots' components + +function plot_title(vis, width, margin_top, title) { + vis.append('svg:foreignObject') + .attr('width', width) + .attr('height', margin_top) + .attr('fill', 'none') + .append('xhtml:div') + .style('text-align', 'center') + .html(title); + + return vis; +} + +function plot_x_axis(vis, height, width, margins, xAxis, label) { + // X axis declaration + vis.append('svg:g') + .attr('class', 'x axis') + .attr('transform', 'translate(0,' + (height - margins.bottom) + ')') + .call(xAxis); + + // X axis label. + vis.append('text') + .attr('class', 'x label') + .attr('fill', 'black') + .attr('text-anchor', 'middle') + .attr('x', (width + margins.right) / 2) + .attr('y', height * 0.97) + .text(label); + + return vis; +} + +function plot_y_axis(vis, height, width, margins, yAxis, label) { + // Y axis declaration. + vis.append('svg:g') + .attr('class', 'y axis') + .attr('transform', 'translate(' + margins.left + ',0)') + .call(yAxis); + + // Y axis label. + vis.append('svg:text') + .attr('class', 'y label') + .attr('fill', 'black') + .attr('transform', 'rotate(-90, 0,' + height + ')') + .attr('text-anchor', 'middle') + .attr('x', (height + margins.bottom) / 2) + .attr('y', (height + margins.left) * 0.92) + .text(label); + + return vis; + +} + +function plot_scenario_data(vis, data, xTimeRange, yRange, line_color) { + var lineFunc = d3.line() + .defined(d => !isNaN(d.concentration)) + .x(d => xTimeRange(d.time)) + .y(d => yRange(d.concentration)) + .curve(d3.curveBasis); + + vis.append('svg:path') + .attr('d', lineFunc(data)) + .attr("stroke", line_color) + .attr('stroke-width', 2) + .attr('fill', 'none'); + + return vis; } \ No newline at end of file diff --git a/cara/apps/calculator/templates/base/calculator.report.html.j2 b/cara/apps/calculator/templates/base/calculator.report.html.j2 index ad04ffae..146478dc 100644 --- a/cara/apps/calculator/templates/base/calculator.report.html.j2 +++ b/cara/apps/calculator/templates/base/calculator.report.html.j2 @@ -88,9 +88,9 @@

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

    @@ -108,8 +108,12 @@
    - - + + {% block report_scenarios_summary_table %}